Welcome to the PyCAS Crash Course! Use SHIFT+ENTER to jump through the Notebook.
We wanted to make working with PyPropagate as easy as possible - even for inexperienced programmers - without loosing performance to inefficiently written python scripts. The solution was to use a symbolic library to define all simulation parameters. These are then optimized and compiled to low-level code on the fly. There is only one well-known symbolic library for python: sympy.
In [1]:
import sympy as sp
However, working with sympy turned out to be quite complicated and soon some serious issues arose. For example, lets try to define a material with the complex refractive index $1+10^{-14}\;i$ in sympy:
In [2]:
index_of_refraction = sp.S(1 + 10**-14 * 1j)
index_of_refraction.simplify()
Out[2]:
Yikes! It completely forgot about the imaginary part (at least in the current stable release: 0.7.6.1). As it turns out, this and other bugs make using sympy for serious tasks quite risky. Unfortunately, there are no known alternatives for python, so we had to write our own: the symbolic library Expresso and minimal computer algebra system PyCAS were born!
In [3]:
import expresso.pycas as pc
In [4]:
index_of_refraction = pc.S(1 + 10**-14 * 1j)
index_of_refraction.evaluate()
Out[4]:
In [5]:
wavelength = pc.Symbol('lambda') # wavelength because lambda is a reserved keyword in python
x,y,z = pc.symbols('x,y,z')
We can create expressions by using some built-in python operators such as +-*/ or by using predefined functions in PyCAS.
In [6]:
pc.sqrt( pc.sin(x) + z / (wavelength + y) )
Out[6]:
We can also use numbers in our expressions, which will automatically be converted to PyCAS symbols. Be aware that operations on numbers will still first be evaluated in python before they are used in a PyCAS function. To explicitly convert a number to PyCAS we can use the pc.S function.
In [7]:
x - z * 2 / 3
Out[7]:
In [8]:
x + 2 / 3 # Here the expression 1/3 is first evaluated by python (and integer division yields 0)
Out[8]:
In [9]:
x + pc.S(2) / 3 # The whole expression is converted to PyCAS
Out[9]:
In [ ]:
n = ### <- insert integer here
pc.S(1./n)
In [10]:
expr = ((-2*x**(z)*(y+2*y/6)/x**(z-1)-y*(x+x/3)+42*3*x**2/(3*x*y**((x-2*x)/x)))+4*x*y)/(x*y)
expr
Out[10]:
In [11]:
expr.evaluate()
Out[11]:
In [12]:
expr = pc.derivative(pc.exp(x*y)*x**2,x)
expr.evaluate()
Out[12]:
To access the arguments of a function we can use the expression.args member of an expression. For commuatative functions the arguments are not necessarily in the same order as in the definition or when its printed. To get the callable function itself we can use the expression.function member.
In [13]:
expr.evaluate().args
Out[13]:
In [14]:
expr.evaluate().args[1].function(x,y)
Out[14]:
It is possible to define your own functions using the pc.Function class.
In [15]:
f = pc.Function('f')
f(x**2+y)
Out[15]:
Values inside expressions can be substituted using the expression.subs method. It is possible to perform multiple substitutions at once.
In [16]:
(x*y**2+y*wavelength).subs(y,wavelength).evaluate()
Out[16]:
In [17]:
(x*(x+y)).subs([(x,y),(y,x)])
Out[17]:
expression.args method multiple times to extract the $x>0$ condition from piecewise_defined_function? Note that in the special case of piecewise there is a hidden outer function.
In [ ]:
a,b = ### <- define symbols here
piecewise_defined_function = pc.piecewise((a,x>0),(b,True))
piecewise_defined_function
In [ ]:
piecewise_defined_function.subs(x,0).evaluate()
In [ ]:
piecewise_defined_function.args[0]
In [ ]:
In [18]:
(pc.sqrt(2)).N(100)
Out[18]:
In [19]:
(4*pc.atan(1)).N(1000) # accurate to the first 1000 digits of pi
Out[19]:
Expressions that contain symbols can be compiled to numpy (or c) functions that can operate on huge arrays very efficiently. The compiled function will be applied element-wise on the arrays. If the return value of the function cannot be deduced, the function result is casted to complex (this can be overwritten by defining the restype argument of pc.numpyfy). The compiled functions take the name of the symbols contained in the expression as keyword arguments.
In [20]:
# import numpy and a pllotting library
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [21]:
func = pc.piecewise((pc.sin(2*x)**5,(pc.sin(x)>0) ),(pc.sin(4*x)**2,True))
func
Out[21]:
In [22]:
numpy_function = pc.numpyfy(func)
plt.plot(numpy_function(x = np.linspace(0,3*np.pi,1000)).real)
Out[22]:
In [23]:
func = pc.sin(x)*pc.sin(y)
func
Out[23]:
In [24]:
numpy_function = pc.numpyfy(func,restype=float)
ny,nx = np.meshgrid(np.linspace(0,3*np.pi,1000),np.linspace(0,3*np.pi,1000))
plt.imshow(numpy_function(x = nx,y = ny))
plt.colorbar()
Out[24]:
In [ ]:
In [ ]:
We can embed numpy arrays in our expressions using the pc.array(name,array) function definition. This will become useful in case an expression is too complicated or computationally intensive to be evaluated explicitly in the expression. Arra access operations are created by calling the array function just lika an ordinary function. If the argument is less than 0 or larger or equal ot the size of the array, the array access returns 0.
In [25]:
rand = pc.array('rand',np.random.rand(100))
rand(x*100)
Out[25]:
In [26]:
plt.plot(pc.numpyfy(rand(x*10),restype=float)(x=np.linspace(-2,12,1000)))
Out[26]:
In [ ]:
def zone_radius(i,wavelength = 1,F = 1):
return (i*wavelength*F + i**2*wavelength**2/4)**0.5
i = pc.Symbol('i')
zone_radius(i).evaluate()
In [ ]: