Load Symata like this
In [1]:
using Symata;
Symata is a computer language for symbolic mathematics written in Julia. After typing using Symata
, you are in Symata mode and input is interpreted as Symata language expressions.
Here is a Symata expression
In [2]:
(x+y)^3
Out[2]:
In a Jupyter notebook, you can exit Symata mode and enter Julia mode by entering Julia()
. Expressions are then Julia language expressions.
In [3]:
Julia()
In [4]:
length(zeros(10)) == 10 # A Julia expression
Out[4]:
In Jupyter, type isymata()
to leave Julia mode and enter Symata
mode. At the command line REPL, type =
at the beginning of a line to enter Symata
mode.
In [5]:
isymata() # we enter Symata mode again
Note: to leave Symata
mode and return to Julia
mode, type Julia()
in IJulia
, or backspace at the command line REPL.
Now we can enter Symata
expressions.
In [6]:
expr = Cos(π * x)
Out[6]:
A variable is set like this.
In [7]:
x = 1/3
Out[7]:
Now the cosine is evaluated
In [8]:
expr
Out[8]:
In [9]:
x = 1/6
Out[9]:
In [10]:
expr
Out[10]:
Clear the value of x
and the cosine can no longer be reduced.
In [11]:
Clear(x)
expr
Out[11]:
It is clear what happened above. If x
is not set to a value, then Cos(π * x)
can't be written in a simpler form. If we set x
to some particular values, then Cos(π * x)
can be reduced to a simpler form.
(you can skip the following the first time through)
But, the reason Symata
understands this is a consquence of the procedure it follows in evaluating (almost) all expressions. Symata
evaluates expressions to a fixed point. More precisely, when an expression is given as input, Symata
descends depth-first evaluating each subexpression to a fixed point and finally the top-level expression to a fixed point. When Cos(π * x)
is first evaluated, each of π
and x
evaluates to itself so that π * x
is already at a fixed point. Since there is no rule for evaluating Cos(π * x)
for fixed π * x
, Cos(π * x)
is also at a fixed point.
The expression x=1/3
means that, whenever x
is encountered, it should evaluate to 1/3
. The expression Out(4)
evaluates to the fourth output cell, which is Cos(π * x)
. Then π
evaluates to iteself, x
evaluates to 1
, so that π * x
evaluates to π/3
. There is a rule saying that Cos(π/3)
evaluates to 1/2
.
Clear(x)
says that x
should once again evaluate to itself. Then evaluating Out(4)
follows the same evaluation sequence, leading to Cos(π * x)
There are several kinds of assignment in Symata. The two most common are =
(or Set
) and :=
(or SetDelayed
).
Set
immediatley evaluates the right hand side and binds the left hand side to the result. SetDelayed
does not evaluate the right hand side when the assignment is made. It evaluates the right hand side every time the left hand side is subsequently evalutated and then binds the result.
The following demonstrates the difference.
In [12]:
x = 1
a := x
b = x
c = a
d := a
In [13]:
[x,a,b,c,d]
Out[13]:
In [14]:
ClearAll(x)
In [15]:
[x,a,b,c,d]
Out[15]:
In [16]:
(a = z, [x,a,b,c,d])
Out[16]:
In [17]:
ClearAll(x,a,b,c,d)
Assign two variables at once
In [18]:
[a,b] = [x,y]
Out[18]:
In [19]:
a
Out[19]:
In [20]:
b
Out[20]:
In [21]:
[a,b]
Out[21]:
Swap two values
In [22]:
[a,b] = [b,a]
Out[22]:
In [23]:
[a,b]
Out[23]:
Every expression has a Head
. For function-like expressions, the Head
is the function name. For atomic expressions, the Head
usually is a data type.
In [24]:
Map(Head, [x, x + y, [x,y], Cos(x), f(x), 3, 3.0, BI(3), BF(3)])
Out[24]:
In [25]:
expr = Expand((x+y)^3)
Out[25]:
In [26]:
FullForm(expr) # This shows the internal form. The tree is explicit
Out[26]:
In [27]:
expr[2,2,1] # Return a part of the expression by index into the tree
Out[27]:
In [28]:
expr[2,2,1] = z; # Replace a part of the expression
In [29]:
expr
Out[29]:
In [30]:
Part(expr,2,2,1) # You can do the same thing with Part
Out[30]:
In [31]:
Expand((x+y)^3)[4,1] # You can get parts of expressions directly
Out[31]:
In [32]:
expr = Expand((x+y)^20); # The semi-colon suppresses printing the return value.
In [33]:
expr[14:18:2] # Parts 14 through 18 with step 2
Out[33]:
In [34]:
ClearAll(expr)
Define a function that collects an expression's head and arguments in a list.
In [35]:
headargs(f_(args__)) := [f,args]
In [36]:
headargs(a + b^2 + 3)
Out[36]:
In [37]:
Integrate(f(x),x)
Out[37]:
In [38]:
headargs(Integrate(f(x),x))
Out[38]:
In [39]:
ClearAll(a,b,c,d)
Rotate the head and arguments to make a new expression.
In [40]:
rotheadargs(f_(args__)) := (Last([args])(f,Splat(Most([args]))))
In [41]:
rotheadargs( a + b + c + d)
Out[41]:
In [42]:
rotheadargs( a + b + c + d + g(x))
Out[42]:
In [43]:
ClearAll(x,a,b,c,d) # delete definitions from the previous example
In [44]:
a = 1
Out[44]:
In [45]:
Definition(a)
In [46]:
a := x
Definition(a) # This overwrites the previous definition
In [47]:
f(x_) := x^2
f(x_, y_) := x + y
Definition(f)
In [48]:
Definition(f)
In [49]:
ClearAll(f,a)
In [50]:
Timing((Range(10^6), Null )) # time a single expression
Out[50]:
Time
toggles the timing of every expression entered. Memory allocated and the number of attempts to apply a user defined rule are also printed.
In [51]:
Time(True) # Enable timing all evaluation. Returns the previous value
Out[51]:
In [52]:
Range(10^6);
In [53]:
Time(False); # disable timing all evaluation.
In [54]:
Trace(True); # Trace evaluation
In [55]:
(a+b)*(a+b)
Out[55]:
In [56]:
Trace(False);
In [57]:
? LeafCount
In [58]:
LeafCount(Expand((a+b)^3))
Out[58]:
In [59]:
ByteCount(Expand((a+b)^3))
Out[59]:
In [60]:
? Depth
In [61]:
Depth(Expand((a+b)^3))
Out[61]:
In [62]:
FullForm(Expand((a+b)^3)) # Examine the tree
Out[62]:
In [63]:
Expand((a+b)^3)[2,2,1] # One of the deepest parts
Out[63]:
In [64]:
Integrate( (1+x^2)^(-1), x)
Out[64]:
In [65]:
expr = 1/(1+x^2)
Integrate(expr, x)
Out[65]:
In [66]:
f(x_) := 1/(1+x^2)
Integrate(expr, x)
Out[66]:
In [67]:
g(x_) = expr # Note we do not use ":="
Out[67]:
In [68]:
ClearAll(expr) # We did not use SetDelay, so we can delete expr
In [69]:
Integrate(g(y),y)
Out[69]:
Note: Trying to use a compiled (Julia) function h = J( x -> 1/(1+x^2))
will not work.
In [70]:
ClearAll(f,g,expr)
# Integrate(f(y), y) # The integral can no longer be reduced
In [71]:
MatchQ(z,_) # Blank matchs z
Out[71]:
In [72]:
Map(MatchQ(_), [1,"string", a+b, 1/3]) # MatchQ does Currying with the first argument
Out[72]:
_head
is a Blank
that only matches expressions with Head
equal to head
.
In [73]:
FullForm(_Integer) # underscore is shorthand for Blank
Out[73]:
In [74]:
MatchQ(1, _Integer)
Out[74]:
Use Currying to define a predicate function
In [75]:
myintq = MatchQ(_Integer);
Not all rational numbers are integers
In [76]:
Map(myintq, Range(1/2,5,1/2))
Out[76]:
In [77]:
MatchQ(b^2, _^2) # Match power with exponent equal to 2
Out[77]:
In [78]:
MatchQ(b^3, _^_) # Match any power
Out[78]:
In [79]:
MatchQ((b+c)^3, _^_)
Out[79]:
In [80]:
MatchQ(b^1, _^_)
Out[80]:
This failed because b^1
evaluates to b
, which does not have the structure of a power
The pattern can be complex with blanks deep in an expression.
In [81]:
Map(MatchQ(f(x_^2)), [f(b^2), f(b^3), g(b^2)])
Out[81]:
Specify a "function" Head
that must match
In [82]:
Map(MatchQ(_gg), [gg(x+y), gg(x), g(x)])
Out[82]:
Define a predicate for positive integers
In [83]:
m = MatchQ(_Integer`Positive`)
Map(m, [1,100, 0, -1, 1.0, x])
Out[83]:
We can also put a Condition
on a Pattern
. This matches pairs with the first element smaller than the second.
In [84]:
m = MatchQ(Condition([x_, y_], x < y))
Map(m, [[2, 1], [1, 2], [1,2,3], 1])
Out[84]:
Pattern
s can include Alternative
s.
In [85]:
m = MatchQ(_Integer | _String)
Map(m, [1, "zebra", 1.0])
Out[85]:
Repeated(expr)
matches one or more occurences of expr
.
In [86]:
MatchQ([a,a,a,b], [Repeated(a),b])
Out[86]:
RepeatedNull
matches zero or more occurences.
In [87]:
MatchQ([b], [RepeatedNull(a),b])
Out[87]:
In [88]:
ClearAll(m)
Rule
s are used for many things in Symata, including replacement. Replacement is a key ingredient in the implementation of functions.
When applied, this Rule
matches and does a replacement on any expression with Head
f
and a List
of two elements as the sole argument.
In [89]:
f([x_, y_]) => p(x + y)
Out[89]:
In [90]:
expr = f([x + y, y]) + f(c) + g([a, b]);
In [91]:
ReplaceAll(expr, f([x_, y_]) => p(x + y))
Out[91]:
There are several things to note here.
The pattern x_
puts no restrictions on the match; any expression will match. The name of the pattern x
only serves to identify it later during a replacement.
Here x_
has matched x+y
, but these two uses of x
are not confused in the result. That is, in x_
, the symbol x
is a dummy variable.
The expression f(c)
has a matching Head
, but not matching arguments, so f(c)
fails to match. Likewise, the expression g([a,b])
has matching arguments, but not matching head.
The expression f([x+y,y])
matches, and the replacement is made in (a copy of) expr
. But, Symata alays evaluates expressions to a fixed point. So y+y
is replaced by 2y
, and the terms in expr
are rearranged into the canonical order.
Again, be aware that matching is structural.
In [92]:
ReplaceAll([a/b, 1/b^2, 2/b^2] , b^n_ => d(n))
Out[92]:
In [93]:
ClearAll(expr)
Named patterns that appear in more than one place must match the same expression.
In [94]:
ReplaceAll([b, a, [a, b]], [x_, y_, [x_, y_]] => 1 ) # This does not match
Out[94]:
In [95]:
ReplaceAll([a, b, [a, b]] , [x_, y_, [x_, y_]] => 1 ) # This does match
Out[95]:
This example uses Alternative
s.
In [96]:
ReplaceAll( [a, b, c, d, a, b, b, b], a | b => x)
Out[96]:
The arguments of Sequence
are spliced into expressions during evaluation.
In [97]:
[1,2,Sequence(a,b)]
Out[97]:
An unmatched alternative is replaced by Sequence()
. Upon evaluation to fixed point, this empty sequence is removed.
In [98]:
f(x_, x_ | y_String) := [x,y]
In [99]:
f(2,2) # `y` does not match, so it is removed.
Out[99]:
In [100]:
f(2,"cat") # Here the second Alternative matches
Out[100]:
In [101]:
f(2,3) # Here the Pattern fails to match.
Out[101]:
Alternative
s, and Pattern
s in general, can be explicit expressions, with no Blank
s.
In [102]:
h(a | b) := p
[h(a), h(b), h(c), h(d)]
Out[102]:
ReplaceAll
replaces all matching subexpressions. We can also specify the levels. This matches at level 2 and deeper.
In [103]:
Replace(1 + a + f(a) + g(f(a)), a => b, 2)
Out[103]:
This replaces only at level 2.
In [104]:
Replace(1 + a + f(a) + g(f(a)), a => b, [2]) == 1 + a + f(b) + g(f(a))
Out[104]:
Rule
evaluates the right hand side once, when it is first evaluated.
In [105]:
ReplaceAll([x, x, x, x, x], x => RandomReal() )
Out[105]:
RuleDelayed
evaluates the right hand side every time it is applie
In [106]:
ReplaceAll([x, x, x, x, x], RuleDelayed(x ,RandomReal()))
Out[106]:
Except
matches everything except expressions that match its argument. The following applies the replacement at level 1.
In [107]:
# FIXME: This example is broken.
Replace([1, 7, "Hi", 3, Indeterminate], Except(_`Numeric`) => 0, 1)
Each Rule
in a List
of Rule
s is tried in turn. Matching stops after the first match. ReplaceRepeated
continues applying Rules
until the expression reaches a fixed point.
In [108]:
ReplaceRepeated(x^2 + y^6 , List(x => 2 + a, a => 3))
Out[108]:
Up to this point, we have used named patterns only with a single blank, for example b_
. But, we may associate a name with any pattern expression, including a complex (compound) expression.
In [109]:
ReplaceAll( b^c, a::(_^_) => g(a))
Out[109]:
Patterns
are used to implement optional arguments.
In [110]:
f(x_, y_:3) := x + y
[f(a+b,z), f(a+b)]
Out[110]:
Condition
may be used in definitions like this:
In [111]:
ClearAll(f)
f(x_) := Condition(x^2, x > 3)
In [112]:
[f(4),f(3)]
Out[112]:
We can match and replace with a Pattern
with Head
equal to Plus
In [113]:
ReplaceAll(z * y + b , x_ + y_ => x * y)
Out[113]:
In [114]:
ReplaceAll(z * y + b + c, x_ + y_ => x * y)
Out[114]:
This failed because Plus
with two terms does not match Plus
with three terms. But, we actually do want this to match. Implementing associative-commuatative matching is a major outstanding problem in Symata. Anyone want to give it a try ?
In [115]:
ClearAll(f,h,a,b,x,y)
Define a compiled function to a built-in or user-defined Julia function like this
In [116]:
mylog = J(log )
mylog(2, 2)
Out[116]:
You can also easily write compiled code like this
In [117]:
f = J((x,y) -> x^2 + y^2)
f(3.0,4.0)
Out[117]:
Note that we did not specify the data types. Is this really high-performance compiled code ? Yes it is. The function was compiled after we called it with two floating point numbers. If we call the function with two integers, a version (called a method) that is optimized for integers is compiled. A version optimized for an integer and a rational number or any combination of arguments can also be compiled.
In [118]:
[f(3,4), f(3, 1/2)]
Out[118]:
We define two versions of the same function to see the difference in performance between compiled functions and functions defined via Rule
s.
In [119]:
(a = Range(0.0,100.,.01), ccossq = J(x -> cos(x)^2), cossq(x_) := cos(x)^2);
In [120]:
Time(True);
We run each test twice because compilation time is included in the first run.
In [121]:
Map(cossq, a);
In [122]:
Map(cossq, a);
In [123]:
Map(ccossq, a);
In [124]:
Map(ccossq, a);
In [125]:
(Time(False), ClearAll(f,a,ccossq,cossq,mylog))
The compiled function is about $500$ times faster in this example. Symata Pattern
s can be compiled (automatically) as well, but this has been removed during a rewriting of the Pattern
code that is still underway. With compiled Pattern
s, the factor might be closer to $50$.
In this example we calculate an expression and compile it. The compiled code is as efficient hand-coded Julia.
In [126]:
expr = Integrate( x^2 * Exp(x)* Cos(x), x)
Out[126]:
In [127]:
expr = Collect(expr, Exp(x))
Out[127]:
Now we use Compile
. By default, Compile
does not evauate its arguments. So we ask explicitly for evaluation.
In [128]:
cexpr = Compile(Evaluate(expr))
a = Range(0.0,10.0,.01);
In [129]:
Timing((Map(cexpr, a), Null)) # big performance regressioin
Out[129]:
In [130]:
Timing((Map(cexpr, a), Null))
Out[130]:
In [131]:
cexpr(2.0)
Out[131]:
In [132]:
ClearAll(a,expr,cexpr)
In [133]:
VersionInfo()
In [134]:
InputForm(Now())
Out[134]: