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]
(The Julia function is a factor of 5 times slower here than in Symata 0.3.0-dev.9)
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 );
Timing(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(range(1, stop=1000.0, length=1000))
Notice that 1000 numbers were not printed. range
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], 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
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(); # return to Symata mode
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(range(1.0, stop=1000, length=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 pack 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.
Symata-mode help does work.
Switch from Symata to Julia mode
Julia();
We will use the following functions and macros
? @sym
? symparsestring
search: symparsestring
? symeval
search: symeval symparseeval symtranseval
? symtranseval
search: symtranseval
? symparseeval
search: symparseeval
? getsymata
search: getsymata
? setsymata
search: setsymata setsymval
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(range(1, stop=4, length=3))
[1.0,2.5,4.0]
Parse a string of Symata code
scode = Meta.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 = 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),2*2^(1/2),3]
isymata() # return to Symata mode
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.
ExportJ(s2);
J(x = 3, y = 2);
J(eval(Meta.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)
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);
copyto!(a, 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.
We verify this by inspecting the lowered code.
import InteractiveUtils
InteractiveUtils.code_native(mmul,(Int,Int))
.text ; ┌ @ arithmetic.jl:60 within `mmul' movq %rdi, %rax ; │┌ @ int.jl:88 within `*' imulq %rsi, %rax ; │└ retq nopl (%rax,%rax) ; └
InteractiveUtils.code_native(*, (Int,Int))
.text ; ┌ @ int.jl:88 within `*' movq %rdi, %rax imulq %rsi, %rax retq nopl (%rax,%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.4.6 Julia version 1.6.0-DEV.58 Python version 3.8.3 SymPy version 1.5.1
InputForm(Now())
2020-05-29T22:35:59.933