

In [1]:

    
using Symata

Interface with Julia

This notebook gives examples of Symata interacting with its host language Julia.

Writing and using Julia functions

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.



In [2]:

    
x1 = Range(10.0^3)
y1 = Range(10.0^3);

Define a procedural function in Symata to compute the sum.



In [3]:

    
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.



In [4]:

    
g(x1,y1)

resultS1 = Timing(g(x1,y1))









    Out[4]:





$$  \left[ 0.069080808,1.6716708333325002e17 \right]  $$

In general, it is faster to use mapping and functional methods in Symata. We can compute the sum like this.



In [5]:

    
Apply(Plus, x1^2 / y1^3)

resultS2 = Timing(Apply(Plus, x1^2 / y1^3))









    Out[5]:





$$  \left[ 0.041056992,7.485470860550343 \right]  $$

The second method is indeed a bit faster.



In [6]:

    
resultS1[1]/resultS2[1]









    Out[6]:





$$ 1.68255891712671 $$

Writing a Julia function while in Symata mode

Define a Julia function to do the sum. We choose a functional method.



In [7]:

    
jfunc = J((x,y) -> sum(u -> u[1]^2 / u[2]^(3), zip(x,y)));



In [8]:

    
jfunc(x1,y1)
Timing(jfunc(x1,y1))

resultJ = Timing(jfunc(x1,y1))









    Out[8]:





$$  \left[ 0.000671538,7.485470860550343 \right]  $$

The Julia function is much faster.



In [9]:

    
[resultS1[1], resultS2[1] ] / resultJ[1]









    Out[9]:





$$  \left[ 102.86954424023658,61.138747174396684 \right]  $$

(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.



In [10]:

    
jfunc([a+b, c+d],[u+v,y+z])









    Out[10]:





$$ \frac{ \left( a + b \right) ^{2}}{ \left( u + v \right) ^{3}} + \frac{ \left( c + d \right) ^{2}}{ \left( y + z \right) ^{3}} $$



In [11]:

    
ClearAll(x,a)

Evaluating Symata expressions from Julia

The macro @sym evaluates Symata code while in Julia.

First, we switch to Julia mode.



In [12]:

    
Julia();

Create a Symata expression and bind it to the Julia variable expr.



In [13]:

    
expr = @sym a + b









    Out[13]:





:a + :b

There are julia functions corresponding to many Symata expression heads.



In [14]:

    
Expand(expr^2)









    Out[14]:





:a^2 + :b^2 + 2:a*:b

Return to Symata mode



In [15]:

    
isymata()

Calling an existing Julia function

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.



In [16]:

    
Julia();

Everything we type will be interpreted as Julia language expressions. We write the Julia function. We will explain later how the function works.



In [17]:

    
fj(x,y) = sum(u -> u[1]^2 / u[2]^(3), zip(x,y));

Return to Symata mode.



In [18]:

    
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)



In [19]:

    
fj = J( Main.fj );

Timing(fj(x1,y1))









    Out[19]:





$$  \left[ 0.034632529,7.485470860550343 \right]  $$

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.



In [20]:

    
J(time)()









    Out[20]:





$$ 1.532087489793153e9 $$

Julia functions for Symata

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.



In [21]:

    
y2 = J(range(1, stop=1000.0, length=1000))









    Out[21]:





$$ 1.0:1.0:1000.0 $$

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.



In [22]:

    
jfunc(x1,y2)
resultJ2 = Timing(jfunc(x1,y2))









    Out[22]:





$$  \left[ 0.000486682,7.485470860550343 \right]  $$

Notice that summing over the two types of arrays is a bit slower in this case than using two Symata arrays.



In [23]:

    
resultJ[1]/resultJ2[1]









    Out[23]:





$$ 1.3798291286712885 $$

Now we call jfunc on two Julia abstract arrays.



In [24]:

    
jfunc(y2,y2)
resultJ3 = Timing(jfunc(y2,y2))









    Out[24]:





$$  \left[ 0.00012182,7.485470860550343 \right]  $$

Operating on these Julia arrays is about 10 times faster in this case than including a Symata array.



In [25]:

    
resultJ2[1]/resultJ3[1]









    Out[25]:





$$ 3.995091118043014 $$

Compile and SymataCall

Use Compile to compile a Symata expression to Julia



In [26]:

    
g = Compile([x], x^2);
g(3)









    Out[26]:





$$ 9 $$

Use SymataCall to create a Julia "callback" to Symata.



In [27]:

    
g = SymataCall(x, x^2 );
g(3)









    Out[27]:





$$ 9 $$

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.

@symExpr macro for writing Julia code with Symata

Here is a Symata expression that evaluates to an expression that could be Julia code.



In [28]:

    
ex1 = Together(PolyLog(-1,z),(1-z))









    Out[28]:





$$ \frac{z}{ \left( -1 + z \right) ^{2}} $$



In [29]:

    
ex1 ./ ( z => 3)









    Out[29]:





$$ \frac{3}{4} $$

Switch to Julia mode



In [30]:

    
Julia()

The macro @symExpr evaluates the Symata expression, translates the result to a Julia expression and inserts it into the surrounding Julia expression.



In [31]:

    
? @symExpr



In [32]:

    
f1(z) = @symExpr Together(PolyLog(-1,z),(1-z))









    Out[32]:





f1 (generic function with 1 method)



In [33]:

    
f1(3)









    Out[33]:




3//4

The function method definition above is equivalent to

f1(z) = z/(z-1)^2



In [34]:

    
isymata();  # return to Symata mode

`Unpack`: Importing Julia arrays into Symata

Recall the abstract Julia array that we created above



In [35]:

    
y2









    Out[35]:





$$ 1.0:1.0:1000.0 $$

We import y2 into Symata using Unpack. (Unpack currently works only with one dimensional arrays.)



In [36]:

    
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.



In [37]:

    
[Head(y3), Length(y3), y3[1], y3[-1]]









    Out[37]:





$$  \left[ \text{List},1000,1.0,1000.0 \right]  $$

y3 is equal to y1, which was created with Range.



In [38]:

    
y3 == y1 == Range(10.0^3)









    Out[38]:





$$ \text{True} $$

y3 is not a packed array, but an ordinary Symata array, a list.



In [39]:

    
y3[1] = "cat"









    Out[39]:





$$ \text{"cat"} $$



In [40]:

    
y3[1:10]









    Out[40]:





$$  \left[ \text{"cat"},2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0 \right]  $$

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.



In [41]:

    
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



In [42]:

    
[Head(y4), y4[1]]









    Out[42]:





$$  \left[ \text{Array{Float64,1}},1.0 \right]  $$

You can also construct a Symata expression directly in Julia, like this.



In [43]:

    
y5 = J(mxpr(:List, Any[collect(range(1.0, stop=1000, length=1000))...]));



In [44]:

    
[Head(y5), y5[-1]]









    Out[44]:





$$  \left[ \text{List},1000.0 \right]  $$

`Pack`: Converting Symata arrays to Julia arrays.

Recall that y1 and y3 are both Symata lists.



In [45]:

    
[Head(y1), Head(y3)]









    Out[45]:





$$  \left[ \text{List},\text{List} \right]  $$

We pack these lists to Julia arrays like this.



In [46]:

    
y6 = Pack(y1)
y7 = Pack(y3);

What type of object was created ?



In [47]:

    
[Head(y6), Head(y7)]









    Out[47]:





$$  \left[ \text{Array{Float64,1}},\text{Array{Any,1}} \right]  $$

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.

Using Symata in Julia code

NOTE Julia-mode help is currently broken (in both IJulia, and the Symata-IJulia-Julia mode)

Symata-mode help does work.

Switch from Symata to Julia mode



In [48]:

    
Julia();

We will use the following functions and macros



In [49]:

    
? @sym



In [50]:

    
? symparsestring









    



search: symparsestring



In [51]:

    
? symeval









    



search: symeval symparseeval symtranseval



In [52]:

    
? symtranseval









    



search: symtranseval



In [53]:

    
? symparseeval









    



search: symparseeval



In [54]:

    
? getsymata









    



search: getsymata



In [55]:

    
? setsymata









    



search: setsymata setsymval

The macro @sym interprets and evaluates its argument in Symata.



In [56]:

    
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



In [57]:

    
@sym z = "cat"
@sym z









    Out[57]:





"cat"

But @sym will not work inside a function. Use getsymata and setsymata instead



In [58]:

    
getsymata(:z)









    Out[58]:





"cat"



In [59]:

    
setsymata(:z, "dog")
getsymata(:z)









    Out[59]:





"dog"

unpacktoList converts a Julia array to a Symat list.



In [60]:

    
unpacktoList(range(1, stop=4, length=3))









    Out[60]:





[1.0,2.5,4.0]

Parse a string of Symata code



In [61]:

    
scode = Meta.parse("Sqrt(a)")









    Out[61]:





:(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.



In [62]:

    
res = symtranseval(scode)









    Out[62]:





3^(1//2)

Print this as Symata would



In [63]:

    
symprintln(res)









    



3^(1/2)

Put this together in a function



In [64]:

    
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]



In [65]:

    
isymata()  # return to Symata mode

You can translate Symata to Julia like this



In [66]:

    
ToJuliaString(3 * x^2 * y^3 + Cos(1))









    Out[66]:





$$ \text{"3 * x ^ 2 * y ^ 3 + cos(1)"} $$

The following gives code that will only work when the Symata module is loaded.



In [67]:

    
s2 = ToJuliaString( 3*x^2*y^3 + Cos(1), NoSymata => False)









    Out[67]:





$$ \text{"mplus(mmul(3, mpow(x, 2), mpow(y, 3)), Cos(1))"} $$

We can parse and evaluate this string in Julia.



In [68]:

    
ExportJ(s2);
J(x = 3, y = 2);
J(eval(Meta.parse(Main.s2)))









    Out[68]:





$$ 216 + \text{Cos} \!  \left( 1 \right)  $$

Symata expressions are of type `Mxpr`

We assigned a value to a in symata



In [69]:

    
Julia()



In [70]:

    
a = getsymata(:a)









    Out[70]:





[1,2,3,4,5,6,7,8,9]



In [71]:

    
typeof(a)









    Out[71]:





Mxpr{:List}

mhead returns the head of a symata expression.



In [72]:

    
mhead(a)









    Out[72]:





:List

margs returns the arguments of a Symata expression.



In [73]:

    
margs(a)









    Out[73]:





9-element Array{Any,1}:
 1
 2
 3
 4
 5
 6
 7
 8
 9

Use `mxpr` to construct a Symata expression



In [74]:

    
ex = mxpr(:Cos, mxpr(:Times, :Pi, 2))









    Out[74]:





Cos(:Pi*2)

mxpr creates an object but does not evaluate it. Evaluate ex with symeval.



In [75]:

    
symeval(ex)









    Out[75]:




1

It may be more efficient to create a Symata expression by first filling an array of arguments



In [76]:

    
a = newargs(5);



In [77]:

    
copyto!(a, 1:5)









    Out[77]:





5-element Array{Any,1}:
 1
 2
 3
 4
 5



In [78]:

    
mxpra(:List,a)









    Out[78]:





[1,2,3,4,5]



In [79]:

    
mxpra(:Plus,a)









    Out[79]:





1 + 2 + 3 + 4 + 5



In [80]:

    
symeval(mxpra(:Plus,a))









    Out[80]:




15

Instead of mxpr, we used mxpra, which does not copy the array of arguments.

We can create Symata expressions like this



In [81]:

    
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 inspectin the lowered code.



In [82]:

    
import InteractiveUtils  
InteractiveUtils.code_native(mmul,(Int,Int))









    



	.text
; Function mmul {
; Location: arithmetic.jl:60
; Function *; {
; Location: julia_level.jl:155
; Function *; {
; Location: arithmetic.jl:60
	imulq	%rsi, %rdi
;}}
	movq	%rdi, %rax
	retq
	nopl	(%rax,%rax)
;}



In [83]:

    
InteractiveUtils.code_native(*, (Int,Int))









    



	.text
; Function * {
; Location: int.jl:54
	imulq	%rsi, %rdi
	movq	%rdi, %rax
	retq
	nopl	(%rax,%rax)
;}

Many Symata functions have equivalents in Julia. These functions both construct Symata expressions and evaluate them.



In [84]:

    
Cos(mmul(2,Pi))









    Out[84]:




1

symmatamath() defines methods allowing you to use * for mmul, etc.



In [85]:

    
symatamath()



In [86]:

    
Cos(2Pi)









    Out[86]:




1

Note that Cos is replaced by efficient Julia methods when possible.



In [87]:

    
isymata()

Version and date



In [88]:

    
VersionInfo()









    



Symata version     0.4.1-dev.3
Julia version      0.7.0-beta2.1
Python version     2.7.14+
SymPy version      1.0



In [89]:

    
InputForm(Now())









    Out[89]:





2018-07-20T13:51:52.167