using Symata
This notebook gives examples of Symata interacting with its host language Julia.
It is easy to write elegant, fast, powerful, flexible, Julia code from within Symata.
As a first example, we multiply powers of elements in two arrays and sum the results. We first try two different methods using pure Symata. Then we write a Julia function to do the same thing.
Here are the arrays.
x1 = Range(10.0^3)
y1 = Range(10.0^3);
Define a procedural function in Symata to compute the sum.
g(x_, y_) := Module([s=0],
begin
For(i=1, i<=Length(x), i += 1, s += x[i]^2 * y[i]^(-3))
s
end)
Apply this function and time the result.
g(x1,y1)
resultS1 = Timing(g(x1,y1))
In general, it is faster to use mapping and functional methods in Symata. We can compute the sum like this.
Apply(Plus, x1^2 / y1^3)
resultS2 = Timing(Apply(Plus, x1^2 / y1^3))
The second method is indeed a bit faster.
resultS1[1]/resultS2[1]
Define a Julia function to do the sum. We choose a functional method.
jfunc = J( (x,y) -> sum(u -> u[1]^2 / u[2]^(3), zip(x,y)) );
jfunc(x1,y1)
Timing(jfunc(x1,y1))
resultJ = Timing(jfunc(x1,y1))
The Julia function is much faster.
[resultS1[1], resultS2[1] ] / resultJ[1]
We will explain later why it is possible to write such simple and fast code that operates on Symata expressions.
The function f
also works on symbolic expressions.
jfunc([a+b, c+d],[u+v,y+z])
ClearAll(x,a)
The macro @sym
evaluates Symata code while in Julia.
First, we switch to Julia mode.
Julia();
Create a Symata expression and bind it to the Julia variable expr
.
expr = @sym a + b
:a + :b
There are julia functions corresponding to many Symata expression heads.
Expand(expr^2)
:a^2 + :b^2 + 2:a*:b
Return to Symata mode
isymata()
To define the Julia function, we used the Symata function J()
. The arguments of J
are interpreted as pure Julia code, with no translation. It is as if we temporarily enter Julia mode. In fact, we could have defined the function in Julia. Let's try that.
First, we enter Julia mode.
Julia();
Everything we type will be interpreted as Julia language expressions. We write the Julia function. We will explain later how the function works.
fj(x,y) = sum(u -> u[1]^2 / u[2]^(3), zip(x,y));
Return to Symata mode.
isymata();
We set the Symata variable fj
to the Julia function fj
. The Julia function was written in the Main
module. (NB we may change this so Julia functions are evaluated in the Symata module)
fj = J( Main.fj );
fj(x1,y1)
We can use J()
in this way to call any existing Julia function...
time()
Get the system time in seconds since the epoch, with fairly high (typically,
microsecond) resolution.
J(time)()
How does the Julia function
jfunc = (x,y) -> sum(u -> u[1]^2 / u[2]^(3), zip(x,y))
work ?
As in Symata, x -> body
defines a pure, or anonymous, function. zip
returns a list of pairs of elements from two lists. In fact, it returns a virtual list, called an iterator, which is more efficient. These pairs are supplied sequentially to the to the function, and the results are summed. No intermediate arrays are formed.
Symata expressions are iterable objects in Julia. Most Julia code that operates on iterable objects will work with Symata expressions. zip
takes two iterable objects and returns an iterable object.
In Julia, Symata expressions are of type Mxpr
. Notice that we did not write Mxpr
anywhere in the code. The first time jfunc
is called with Symata expressions, Julia compiles a method to handle just this case. The compiler is typically very good at writing code optimized for the input type.
All of this means that the author of Symata wrote no code to implement zip
or sum
for Symata expressions.
In fact jfunc
can be called with many types objects. To demonstrate this, we perform the sum operation on a Symata list and a Julia Array.
We set the Symata variable y2
to a Julia Array
of 1000 numbers.
y2 = J(linspace(1,1000.0,1000))
Notice that 1000 numbers were not printed. linspace
returns a virtual array, that is an iterator. We call jfunc
twice. The first time, Julia compiles a method for the input types which takes some (not much) time.
jfunc(x1,y2)
resultJ2 = Timing(jfunc(x1,y2))
Notice that summing over the two types of arrays is a bit slower in this case than using two Symata arrays.
resultJ[1]/resultJ2[1]
Now we call jfunc
on two Julia abstract arrays.
jfunc(y2,y2)
resultJ3 = Timing(jfunc(y2,y2))
Operating on these Julia arrays is about 10 times faster in this case than including a Symata array.
resultJ2[1]/resultJ3[1]
Use Compile
to compile a Symata expression to Julia
g = Compile( x^2 );
g(3)
Use SymataCall
to create a Julia "callback" to Symata.
g = SymataCall(x, x^2 );
g(3)
But, Compile
and SymataCall
create different functions. Compile
translates Symata to Julia and compiles the result. SymataCall
wraps the Symata expression in a Julia function that sets the variable and then Symata-evaluates the result. The function returned by Compile
executes faster in general. But SymataCall
works with expressions that cannot be translated and compiled.
NIntegrate
effectively uses SymataCall
.
Here is a Symata expression that evaluates to an expression that could be Julia code.
ex1 = Together(PolyLog(-1,z),(1-z))
ex1 ./ ( z => 3)
Switch to Julia mode
Julia()
The macro @symExpr
evaluates the Symata expression, translates the result to a Julia expression and inserts it into the surrounding Julia expression.
? @symExpr
@symExpr(expr)
Symata-evaluates expr
, translates the result to a Julia expression, and inserts it into the surrounding code. This works just like any Julia macro for inserting code, except that the code is generated by evaluating a Symata expression.
The following Symata code returns a rational expression, which is used as the body of a Julia function.
f1(z) = @symExpr Together(PolyLog(-1,z),(1-z))
f1(z) = @symExpr Together(PolyLog(-1,z),(1-z))
f1 (generic function with 1 method)
f1(3)
3//4
The function method definition above is equivalent to
f1(z) = z/(z-1)^2
isymata()
Unpack
: Importing Julia arrays into Symata¶Recall the abstract Julia array that we created above
y2
We import y2
into Symata using Unpack
. (Unpack
currently works only with one dimensional arrays.)
y3 = Unpack(y2);
The result is a Symata list. We check that it is indeed a list of the expected length and the first and last elements.
[Head(y3), Length(y3), y3[1], y3[-1]]
y3
is equal to y1
, which was created with Range.
y3 == y1 == Range(10.0^3)
y3
is not a packed
array, but an ordinary Symata array, a list.
y3[1] = "cat"
y3[1:10]
Of course, we can unpack not just abstract Julia arrays, but physical arrays as well. The Julia function for converting an abstract array to a physical array is collect
.
y4 = J(collect)(y2);
Notice what we did here. J(collect)
gets a Symata reference to the Julia function collect
. We then call the imported function on the Symata variable y2
, which refers to an abstract Julia array.
Now we have a physical Julia array
[Head(y4), y4[1]]
You can also construct a Symata expression directly in Julia, like this.
y5 = J(mxpr(:List, Any[collect(linspace(1.0,1000,1000))...]));
[Head(y5), y5[-1]]
Pack
: Converting Symata arrays to Julia arrays.¶Recall that y1
and y3
are both Symata lists.
[Head(y1), Head(y3)]
We unpack these lists to Julia arrays like this.
y6 = Pack(y1)
y7 = Pack(y3);
What type of object was created ?
[Head(y6), Head(y7)]
The first array is of type Float64
and the second of type Any
. They are different because, while all elements of y1
are floating point numbers, we set the first element of y7
to "cat"
. When copying arrays, Julia creates a container of the most specific type that will contain all elements.
# isymata()
Switch from Symata to Julia mode
Julia();
We will use the following functions and macros
? @sym
@sym expr
Embed symata expression expr
in Julia code.
Read and evaluate expr
embedded in Julia code.
julia> a = 1 # Julia symbol `a`
julia> @sym a = 2 # Symata symbol
julia> a
1
julia> @sym a
2
? symparsestring
search: symparsestring
symparsestring(s::String)
parses s
into one or more Julia expressions, translates each one to a Symata expression, and returns an array of the Symata expressions.
Note that the phrases Symata expression and Julia expression here include numbers, symbols, etc.
? symeval
search: symeval symparseeval symtranseval
symeval(expr::Any)
send expr
through the Symata evaluation sequence. expr
is an Mxpr
or number, symbol, etc.
In particular, an Expr
(i.e. not translated to Symata) Symata-evaluated (whic in this case means returned unchanged) by the Symata evaluation sequence.
? symtranseval
search: symtranseval
symtranseval(expr::Any)
translate expr
from Expr
to Mxpr
and send to the top level of the Symata evaluation sequence. expr
may also be a number, symbol, string, etc.
? symparseeval
search: symparseeval
symparseeval(s::String)
parses s
into Julia expressions, translates them to Symata expressions, and Symata-evaluates each one, returning the value returned by the final evaluation.
? getsymata
search: getsymata
getsymata(s::Symbol)
gets the value that s
is bound to in Symata.
? setsymata
search: setsymata
setsymata(s::Symbol, val)
binds val
to s
in Symata.
The macro @sym
interprets and evaluates its argument in Symata.
a = 1 # set a to 1 in Julia
@sym a = 3 # set a to 3 in Symata
println("a in Julia is $a") # print a in Julia
@sym Println("a in Symata is $a") # print a in Symata
a in Julia is 1 a in Symata is 3
You can assign a value in Symata by using @sym
@sym z = "cat"
@sym z
"cat"
But @sym
will not work inside a function. Use getsymata
and setsymata
instead
getsymata(:z)
"cat"
setsymata(:z, "dog")
getsymata(:z)
"dog"
unpacktoList
converts a Julia array to a Symat list.
unpacktoList(linspace(1,4,3))
[1.0,2.5,4.0]
Parse a string of Symata code
scode = parse("Sqrt(a)")
:(Sqrt(a))
This is a valid Julia expression, although evaluating it in Julia may cause an error. We translate the expression to Symata and send it through the Symata evaluation sequence.
res = symtranseval(scode)
3^(1//2)
Print this as Symata would
symprintln(res)
3^(1/2)
Put this together in a function
function squareroots()
a = [i for i in 1:9]
setsymata(:a, unpacktoList(a))
symprintln(symparseeval("Sqrt(a)"))
nothing
end
squareroots()
[1,2^(1/2),3^(1/2),2,5^(1/2),2^(1/2)*3^(1/2),7^(1/2),22^(1/2),3]
isymata()
You can translate Symata to Julia like this
ToJuliaString( 3*x^2*y^3 + Cos(1))
The following gives code that will only work when the Symata module is loaded.
s2 = ToJuliaString( 3*x^2*y^3 + Cos(1), NoSymata => False)
We can parse and evaluate this string in Julia.
SetJ(s2,s2);
J( x = 3, y = 2);
J( eval(parse(Main.s2)))
Mxpr
¶We assigned a value to a
in symata
Julia()
a = getsymata(:a)
[1,2,3,4,5,6,7,8,9]
typeof(a)
Symata.Mxpr{:List}
mhead
returns the head of a symata expression.
mhead(a)
:List
margs
returns the arguments of a Symata expression.
margs(a)
9-element Array{Any,1}: 1 2 3 4 5 6 7 8 9
mxpr
to construct a Symata expression¶ex = mxpr(:Cos, mxpr(:Times, :Pi, 2))
Cos(:Pi*2)
mxpr
creates an object but does not evaluate it. Evaluate ex
with symeval
.
symeval(ex)
1
It may be more efficient to create a Symata expression by first filling an array of arguments
a = newargs(5);
copy!(a, [i for i in 1:5])
5-element Array{Any,1}: 1 2 3 4 5
mxpra(:List,a)
[1,2,3,4,5]
mxpra(:Plus,a)
1 + 2 + 3 + 4 + 5
symeval(mxpra(:Plus,a))
15
Instead of mxpr
, we used mxpra
, which does not copy the array of arguments.
We can create Symata expressions like this
println(mmul(3,2), ", " , mmul(:b, :a))
6, :b*:a
mmul
, mpow
, mplus
, mminus
are arithemtic methods that can create Symata expressions. When called inside a function with numerical arguments, the compiler will replace them with an efficent Julia method
code_native(mmul,(Int,Int))
.text Filename: arithmetic.jl pushq %rbp movq %rsp, %rbp Source line: 61 imulq %rsi, %rdi movq %rdi, %rax popq %rbp retq nopl (%rax)
code_native(*, (Int,Int))
.text Filename: int.jl pushq %rbp movq %rsp, %rbp Source line: 33 imulq %rsi, %rdi movq %rdi, %rax popq %rbp retq nopl (%rax)
Many Symata functions have equivalents in Julia. These functions both construct Symata expressions and evaluate them.
Cos(mmul(2,Pi))
1
symmatamath()
defines methods allowing you to use *
for mmul
, etc.
symatamath()
Cos(2Pi)
1
Note that Cos
is replaced by efficient Julia methods when possible.
isymata()
VersionInfo()
symata version 0.3.0-dev.9 julia version 0.6.0-dev.435 python version 2.7.12 sympy version 1.0
Now()