Contents:
help() more helpful (subslides)To Do:
In [1]:
import sympy
In [2]:
print(sympy.__doc__)
In [3]:
help(sympy)
Here are some examples of basic symbol operations:
In [4]:
x = sympy.Symbol('x')
y = x
x, x*2+1, y, type(y), x == y
Out[4]:
In [5]:
try:
x*y+z
except NameError as e:
print(e)
In [6]:
sympy.symbols('x5:10'), sympy.symbols('x:z')
Out[6]:
In [7]:
X = sympy.numbered_symbols('variable')
[ next(X) for i in range(5) ]
Out[7]:
SymPy also handles expressions:
In [8]:
e = sympy.sympify('x*(x-1)+(x-1)')
e, sympy.factor(e), sympy.expand(e)
Out[8]:
But most work of interest is more than single expressions, so here is a helper function to handle systems of equations.
In [9]:
from process import parseExpr
help(parseExpr)
In [10]:
import inspect
print(inspect.getsource(parseExpr))
In [11]:
inputs='x y theta'
outputs="x' y'"
expr='''
x' = x*cos(theta) - y*sin(theta)
y' = x*sin(theta) + y*cos(theta)
'''
ins = sympy.symbols(inputs)
outs = sympy.symbols(outputs)
eqn = dict(parseExpr(expr))
# No quote marks, the dictionary keys are SymPy symbols
ins, outs, eqn
Out[11]:
SymPy maintains a tree for all expressions. Everything in SymPy has .args and .func arguments that allow the expression (at that point in the tree) to be reconstructed. The .func argument is essentially the same as calling type and specifies whether the node is a add, multiply, cosine, some other function or a symbol. As you might expect, the .args attribute for leaves is empty.
In [12]:
expr_inputs = set()
expr_functs = set()
for arg in sympy.preorder_traversal(eqn[outs[0]]):
print(arg.func, '\t', arg, '\t', arg.args)
if arg.is_Symbol:
expr_inputs.add(arg)
elif arg.is_Function:
expr_functs.add(arg.func)
expr_inputs, expr_functs
Out[12]:
In [13]:
from process import Block
print(inspect.getdoc(Block))
b = Block(expr, '2-D Rotate', inputs, outputs)
In [14]:
# spoiler alert!
print(inspect.getsource(Block.__init__))
In [15]:
print('\n'.join( str(k)+' = '+sympy.pretty(v) for k,v in eqn.items() ))
print()
print('\n'.join( str(k)+' = '+str(v) for k,v in eqn.items() ))
So our Block instance can do the same thing, it needs to have pretty and __str__ functions defined. A __repr__ function could also be used to return a separate representation of the object.
In [16]:
b.pretty()
print(b)
In [17]:
print(inspect.getsource(b.pretty))
print(inspect.getsource(b.__str__))
In [18]:
# Generate a LaTeX string for the Jupyter notebook to render
print(' \\\\\n'.join([ str(k)+' = '+sympy.latex(v) for k, v in eqn.items() ]))
In [19]:
%%latex
$
x' = x \cos{\left (\theta \right )} - y \sin{\left (\theta \right )} \\
y' = x \sin{\left (\theta \right )} + y \cos{\left (\theta \right )}
$
For our Block instance, the latex function doesn't need any arguments, so it can be handled as an attribute
In [20]:
print(b.latex)
In [21]:
%%latex
$
\underline{\verb;Block: 2-D Rotate;} \\
x' = x \cos{\left (\theta \right )} - y \sin{\left (\theta \right )} \\
y' = x \sin{\left (\theta \right )} + y \cos{\left (\theta \right )}
$
In [22]:
f = inspect.getsource(Block).split('def ')
for i,s in enumerate(f):
if s.startswith('latex'):
# grab the last 2 lines from the previous code block, suppress the final newline
print('\n'.join(f[i-1].rsplit('\n',2)[-2:]), end="")
print('def '+s.strip())
As a property, the latex function above is implicitly called instead of returning the function itself. Attempting to use inspect.getsource results in a TypeError since the LaTeX output isn't source code.
In [23]:
try:
inspect.getsource(getattr(Block,'latex'))
except TypeError as e:
print(type(e),' : ', e)
We've seen a couple SymPy output formats. The init_printing function provides a lot of additional control over how symbols and expressions are shown, including LaTeX.
In [24]:
b.eqn
Out[24]:
In [25]:
sympy.init_printing(use_latex=True)
#sympy.init_printing(use_latex=False)
In [26]:
b.eqn
Out[26]:
In [27]:
sympy.Matrix(b.eqn)
Out[27]:
In [28]:
help(sympy.init_printing)
In [29]:
inputs='x y theta'
outputs="x' y'"
expr='''
x' = x*c - y*s
y' = x*s + y*c
c = cos(theta)
s = sin(theta)
'''
ins2 = sympy.symbols(inputs)
outs2 = sympy.symbols(outputs)
hidden2 = sympy.symbols('c s')
eqn2 = tuple(sympy.Eq(k,v) for k, v in parseExpr(expr))
ins2, outs2, eqn2
Out[29]:
In [30]:
sympy.solve(eqn2, outs2+hidden2)
Out[30]:
In [31]:
help(sympy.solve)
In [32]:
x = sympy.Symbol("x'")
eqn[x].subs(zip(ins, (1, 1, 45)))
Out[32]:
In [33]:
eqn[x].evalf(4, subs=dict(zip(ins, (1, 1, 45))))
Out[33]:
In [34]:
# sanity check with NumPy
import numpy as np
-1*np.sin(45) + np.cos(45)
Out[34]:
In [35]:
e = sympy.sympify('sqrt(x)')
print(e.evalf(subs={'x':2}))
print(e.evalf(60, subs={'x':2}))
In [36]:
print(sympy.pi.evalf())
print(sympy.pi.evalf(100))
In [37]:
help(eqn[x].subs)
In [38]:
help(eqn[x].evalf)
For efficiency, sympy.lambdify is preferred for numerical analysis. It supports mathematical functions from math, sympy.Function, or mpmath. Since these library functions are compiled Python, C, or even Fortran, they are significantly faster than sympy.evalf.
In [39]:
rad=np.linspace(0, np.pi, 8+1)
f = sympy.lambdify(ins, eqn[x], 'numpy')
%timeit f(1,0,rad)
In [40]:
%%timeit
for i in rad: # evalf doesn't support arrays!
eqn[x].evalf(subs={'x':1.0,'y':0.0,'theta':i})
SymPy also supports uFuncify for generating binary functions (using f2py and Cython) and Theano for GPU support. These options are discussed in the SymPy documentation.
In [41]:
help(sympy.lambdify)
Note that sympy.lambdify also supports custom functions (e.g. conditional operations or reshaping arrays). These custom functions can also be optimized for computing as easily as including a @jit decorator from the numba library.
In [42]:
def my_sample(x):
r = len(x) >> 1
return np.reshape(x[:2*r], (r,2))
x = np.linspace(0, np.pi, 9+1)
print(x*180/np.pi) # degrees
In [43]:
e = sympy.sympify('1+sample(x)')
print(e)
f = sympy.lambdify(sympy.Symbol('x'), e, {'sample':my_sample})
f(x)
Out[43]:
In [44]:
f = sympy.lambdify(sympy.Symbol('x'), # inputs
sympy.sympify('sample(cos(x))'), # expressions
('numpy', {'sample':my_sample})) # functions
# Note: Can also use empty function dictionary for consistency
f(x)
Out[44]:
However, the documentation could be better. (Earlier versions didn't even include the expresion.)
In [45]:
help(f)
For more control over the output function, including the __doc__, __name__, and __signature__ attributes used by the help function, let's revisit the Block class to wrap the functionality seen so far.
In [46]:
b2 = Block('''c = cos(theta)
s = sin(theta)
x' = x*c - y*s
y' = x*s + y*c
''', '2-D Rotate', 'x y theta', "x' y'")
print('\n',vars(b2))
When it comes to function arguments, Python 3.5 has had a lot of development beyond Python 2.7. Since the Block code was originally developed in Python 2.7, let's replace those _Dummy arguments from sympy.lambdify.
In [47]:
f = b2.lambdify()
In [48]:
help(f["x'"])
In [49]:
# All lambdify calls are made with the full set of inputs
p = []
# Update the signatures for each compiled expression
for k in f.keys():
sig=inspect.signature(f[k])
if len(p) == 0:
for i,s in zip(sig.parameters.values(), map(str, b2.ins)):
p.append(i.replace(name=s))
f[k].__signature__ = sig.replace(parameters=p)
help(f[k])
In [50]:
import types
# backup functions before overwriting them
Block.__lambdify = Block.lambdify
b2.__lambdify = b2.lambdify
In [51]:
# Rename the old function
b2._lambdify = b2.lambdify
Block._lambdify = Block.lambdify
#test._lambdify = types.MethodType(test.lambdify, test)
#Block._lambdify = types.MethodType(Block.lambdify, Block)
In [52]:
def lambdify(self, *args, **kwargs):
# Call the base function with any arguments
self._lambdify(*args, **kwargs)
# Update the signatures for each compiled expression
for k in self.lambdas.keys():
sig=inspect.signature(self.lambdas[k])
p = []
for i,s in zip(sig.parameters.values(), map(str, self.ins)):
p.append(i.replace(name=s))
self.lambdas[k].__signature__ = sig.replace(parameters=p)
return self.lambdas
# update the old __doc__ string
lambdify.__doc__ = Block._lambdify.__doc__ + ' and set the arguments to the variable names'
# use the new function in the class and existing instance
b2.lambdify = types.MethodType(lambdify, b2)
Block.lambdify = types.MethodType(lambdify, Block)
In [53]:
f = b2.lambdify()
In [54]:
help(f["x'"])
In [55]:
help(inspect.Signature)
And since it's coming up, let's go ahead and improve the Block.eval function used to call the various lambdify functions. Extra credit to anyone who figures out why I couldn't do the same sort of __signature__ update without getting a NameError exception.
In [56]:
help(b2._eval)
help(b2.eval)
In [57]:
# Backup the old functions
b2._eval_save = b2.eval
Block._wrap_eval_save = Block._wrap_eval
In [58]:
# instead of wrapping _eval, just update its help info
def wrap_eval(self):
'''create _eval() wrapper function with useful docstring'''
def tuple_repr(i):
if type(i) == tuple:
o = i
s = repr(i)
else:
ins = (i,)
s = '('+repr(i)+')'
return o, s
ins, s_ins = tuple_repr(self.ins)
outs, s_outs = tuple_repr(self.outs)
f = self._eval
f.__func__.__name__ = str(self.__class__.__name__)+'('+self.name+').eval'
s = 'Given inputs '+s_ins+', calculate outputs '+s_outs+' as follows:'
for e in self.eqn:
s += '\n '+str(e.lhs)+' = '+str(e.rhs)
s += '\n\nBlock.outputMode = '+repr(Block.outputMode)+'\n'
f.__func__.__doc__ = s
return f
b2.eval = wrap_eval(b2)
# use the new function in the class
Block._wrap_eval = types.MethodType(wrap_eval, Block)
In [59]:
help(b2.eval)
b2.eval(1.0,2.0,3.0)
Out[59]:
The outputs were specified when the Block instance was created, so lambdify all of those.
In [60]:
f = b2.lambdify()
f, b2.lambdas
Out[60]:
In [61]:
help(f["x'"])
The functions were made using sympy.lambdify, so they get all the benefits like support for NumPy arrays. The Block class also has an eval function to use that calls all the lambdify results.
In [62]:
out = f["x'"](1,0,np.linspace(0, np.pi, 8+1))
print(out, type(out))
In [63]:
help(b2.eval)
In [64]:
b2.eval(1,0,np.linspace(0, np.pi, 8+1))
Out[64]:
Potential Updates:
In [65]:
from matplotlib import pyplot as plt
from matplotlib import cm
from matplotlib import colors
# for IPython plotting
%pylab inline
In [66]:
# random (x,y) points
xy1 = np.random.normal(0,0.5,(8,2))
x,y = np.hsplit(xy1,2)
# rotate points on a log-scale
a = [0]+np.logspace(-2,1,20,base=5)
b2.verbose = False
xy2 = b2.eval(x,y,a)
xy2["2-D Rotate.x'"].shape
Out[66]:
In [78]:
def example(xy1, xy2):
'''Move plot code into a function for cleaner presentation'''
# split original points into separate x/y arrays
x,y = np.hsplit(xy1,2)
# Fetch a matplotlib color map to use for line colors
c = cm.get_cmap('gist_heat')
cNorm = matplotlib.colors.Normalize(vmin=0, vmax=20)
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=c)
# calculate "unit" circles for each point
a = np.reshape(np.linspace(0,2*np.pi,60), (60,1))
unit = np.cos(a) + 1j*np.sin(a)
circ = np.abs(x + 1j*y) * unit.T
# set up the figure
plt.figure(figsize=(12,9))
x = xy2["2-D Rotate.x'"].T
y = xy2["2-D Rotate.y'"].T
# "unit" circles for background
for i in range(circ.shape[0]):
s = '|iq{:d}|={:0.3f}'.format(i, np.abs(circ[i][0]))
plt.plot(circ[i].real, circ[i].imag, ls=':', color='k', label=s)
#plt.text(circ[i][0].real, circ[i][0].imag, str(i)+' ')
# plot each rotation point
for i in range(x.shape[0]):
plt.plot(x[i], y[i], color=scalarMap.to_rgba(i), marker='o', ls='none')
# label the first rotation points
for i in range(x.shape[1]):
plt.text(x[0][i]+0.04, y[0][i]+0.04, str(i), color='blue')
plt.axis('equal')
plt.grid()
plt.legend()
In [79]:
example(xy1, xy2)
In [67]:
# Fetch a matplotlib color map to use for line colors
c = cm.get_cmap('gist_heat')
cNorm = matplotlib.colors.Normalize(vmin=0, vmax=20)
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=c)
# calculate "unit" circles for each point
a = np.reshape(np.linspace(0,2*np.pi,60), (60,1))
unit = np.cos(a) + 1j*np.sin(a)
circ = np.abs(x + 1j*y) * unit.T
# set up x/y
x = xy2["2-D Rotate.x'"].T
y = xy2["2-D Rotate.y'"].T
circ.shape
Out[67]:
In [86]:
plt.figure(figsize=(12,6))
for i in range(circ.shape[0]): # "unit" circles for background
s = '|iq{:d}|={:0.3f}'.format(i, np.abs(circ[i][0]))
plt.plot(circ[i].real, circ[i].imag, ls=':', color='k', label=s)
for i in range(x.shape[0]): # plot each rotation point
plt.plot(x[i], y[i], color=scalarMap.to_rgba(i), marker='o', ls='none')
for i in range(x.shape[1]): # label the first rotation points
plt.text(x[0][i]+0.04, y[0][i]+0.04, str(i), color='blue')
plt.axis('equal')
plt.grid()
i = plt.legend()
In [70]:
print(type(b.eqn[0]),'\n')
help(b.eqn[0].subs)
In [71]:
help(sympy.lambdify)
I originally made this notebook using Python 2.7. Fortunately, there's a Python library for that too lib2to3, although the canned 2to3 application worked just fine for me. Although this doesn't usually change functionality, there may be some slight changes as seen when running help on the sympy.lambdify output functions.
In [72]:
%%file process.py
import sympy
import numpy as np
def parseExpr(expr=''):
'''Helper function to iterate through a list of equations'''
err = 'Malformed expression! Does not match "y = f(x)"\n {0:s}'
for s in expr.strip().split('\n'):
# Parse anything that looks like an equation and isn't commented out
if ('=' in s) and (s[0] != '#'):
# convert expression to sympy
y, f = map(str.strip, s.split('=',1))
y = sympy.Symbol(y)
f = sympy.sympify(f)
assert type(y) == sympy.symbol.Symbol, err.format(s)
yield (y, f)
class Block(object):
'''Block(expr, inputs='', outputs='', functions={})
Generic processing block that performs specified calculations on given inputs
in such a way to ease documentation
'''
verbose = True # Enable verbose output
outputMode = 'dict' # eval() return type ('dict' or 'tuple')
def __init__(self, expr, name='', inputs='', outputs='', functions={}):
'''Create a Block instance'''
if Block.verbose:
s = 'Creating Block(expr={:s}, name={:s}, inputs={:s}, outputs={:s}, functions='
s = s.format(*map(repr,[expr, name, inputs, outputs]))+'{'
if functions:
s_functions = []
for k, v in functions.iteritems():
s2 = k+':'+v.__name__
args, varargs, keywords, defaults = getargspec(v)
if varargs != None:
args += ['*args']
if keywords != None:
args += ['*kwargs']
s2 += '('+','.join(args)+')'
s_functions.append(s2)
s += ','.join(s_functions)
print s+'})'
self.name = name
self.user = functions
# save in list form for ordered outputs
eqn = tuple(parseExpr(expr))
self.eqn = tuple([ sympy.Eq(k, v) for k, v in eqn ])
# placeholder for compiled expressions
self.lambdas = {}
# Extract inputs and functions used
expr_inputs = set()
expr_functs = set()
for k, v in eqn:
for arg in sympy.preorder_traversal(v):
if arg.is_Symbol:
expr_inputs.add(arg)
elif arg.is_Function:
expr_functs.add(arg.func)
# save SymPy style .args and .func attributes
self.args = tuple(expr_inputs)
self.func = tuple(expr_functs)
if inputs:
self.ins = sympy.symbols(inputs)
else: # extract inputs from expr
self.ins = tuple(self.args)
outs = tuple([ i[0] for i in eqn ])
if outputs:
self.outs = sympy.symbols(outputs)
print 'outs='+repr(outs)
self.hidden = tuple([ i for i in outs if i not in self.outs ])
else: # extract inputs from expr
self.outs = outs
self.hidden = tuple()
if Block.verbose:
print ' Extracting outputs:', self.outs
# create _eval() wrapper function with useful docstring
self.eval = self._wrap_eval()
def _wrap_eval(self):
'''create _eval() wrapper function with useful docstring'''
if type(self.ins) == tuple:
ins = self.ins
s = repr(self.ins)
else:
ins = (self.ins,)
s = '('+repr(self.ins)+')'
# not sure how to do a decorator with non-generic dynamic args
func = sympy.lambdify(self.ins, '_eval'+s, {'_eval':self._eval})
func.__name__ = 'Block.eval'
s = repr(self.outs)
if type(self.outs) != tuple:
s = '('+s+')'
func.__doc__ = 'Calculate outputs '+s+' as follows:\n '
func.__doc__ += '\n '.join(map(str,self.eqn))
func.__doc__ += '\n\nBlock.outputMode = '+repr(Block.outputMode)+'\n'
return func
def __str__(self):
'''Support str calls on Block instances'''
s = 'Block: '+str(self.name)+'\n'
s += '\n'.join( str(e).replace('==','=') for e in self.eqn )
return s
def __repr__(self):
'''return representation that sympy will automatically show as LaTeX'''
s = 'Block'
if self.name:
s += '('+self.name +')'
return repr({s:[ sympy.Eq(sympy.Symbol(k), v) for k,v in self.eqn.items() ]})
# method doesn't work currently
__repr__ = __str__
def pretty(self):
'''Show a SymPy pretty print version of Block equations'''
print '\nBlock:', self.name
print '\n'.join( sympy.pretty(e) for e in self.eqn )
@property
def latex(self):
'''generate a latex version of Block equations'''
# \verb;*; leaves contents unformatted
s = r'\underline{\verb;Block: '+str(self.name)+';} \\\\ \n'
s += ' \\\\ \n'.join( sympy.latex(e) for e in self.eqn )
return s
def lambdify(self, unknowns=None):
'''generate an executable function using NumPy'''
lookup='numpy'
# by default, unknowns are everything but outputs
if unknowns == None:
unknowns = self.outs+self.hidden
if self.verbose:
print 'Compiling Block '+repr(self.name)+' expressions:'
# check for missing functions
defined = set(dir(np) + self.user.keys())
missing = [ f for f in self.func if str(f) not in defined ]
if missing:
s = 'Unable to find functions '+str(missing)+' in '+str(lookup)
s += ' or user-defined functions ' + self.user.keys()
raise LookupError, s
# solve equations in terms of unknowns
eqn = sympy.solve(self.eqn, unknowns, dict=True)
if self.verbose:
print 'Solving equations for unknowns:', unknowns
print self.eqn
print ' =', eqn, '\n'
s = " {0:s} = sympy.lambdify({1:s}, {2:s}, '{3:s}')" # verbose string
for k, v in eqn[0].iteritems():
if k in self.hidden: # skip intermediate values
continue
k = str(k) # convert output variable from sympy.Symbol
if self.verbose:
print s.format(*map(str, (k, self.ins, v, (lookup, self.user.keys()) )))
f = sympy.lambdify(self.ins, v, (lookup, self.user))
# Update the function name and doc string from "<lambda> lambda x, y, theta"
f.__name__ = "Block"
if self.name:
f.__name__ += '('+self.name+')'
f.__name__ += '.lambdas['+k+']'
f.__doc__ = k+' = '+str(v)
self.lambdas[k] = f
return self.lambdas
def _eval(self, *args, **kwargs):
'''evaluate outputs for given inputs'''
# make sure in/ouputs are iterable
if hasattr(self.ins, '__iter__'):
ins = self.ins
else:
ins = [self.ins]
# make sure all inputs are given
assert len(args) == len(ins)
# default kwargs values
if not kwargs.has_key('output'):
output=Block.outputMode
# compile expressions if needed
if not self.lambdas:
self.lambdify()
if self.verbose:
s = '\nBlock.eval('
s += ', '.join( map(lambda k, v: str(k)+'='+str(v), ins, args) )
s += ', output='+repr(output)
print s+')\n'
if output == 'dict':
if self.name:
prefix=self.name+'.'
else:
prefix=''
return { prefix+k:v(*args) for k, v in self.lambdas.iteritems() }
else: # output tuple in order of .outs
return tuple( self.lambdas[k](*args) for k in self.outs )
Let's test the conversion of the rewritten Python 2.7 code
In [73]:
!2to3 -w process.py
Success!
In [74]:
from process import *
In [75]:
for s in (parseExpr, Block):
help(s)
Python and Jupyter Notebooks
Notebook reveal-based slideshow tutorial
A brief tour of the IPython notebook: Same presentation, just later on
2to3 - Automated Python 2 to 3 code translation
But do I really want to use Python 3? Yes, it's been around since 2008
SymPy and LaTeX
Scientific Python Lectures: Symbolic algebra
Wolfram Alpha-style search engine powered by SymPy SymPy Tutorial
SymPy solvers module: Algebraic/Differential/Recurrence/Diophantine equations, Utilities, Systems of Polynomial Equations Equations, Inequalities
SymPy Simplification Functions: Function/Trigonometric/Powers/Exponentials/logarithms Simplification, Special Functions
SymPy Basic Operations: Substitution, string inputs, evalf, lambdify
SymPy Expression Trees and Manipulations
Embed LaTeX math equations into Microsoft Word
SciPy and matplotlib
SciPy Signal Processing module: Convolution, Filters, Linear Systems, Function Generators, Windowing, Wavelets, Peak finding, Spectral Analysis
SciPy Special Function module: Elliptic, Bessel, Gamma, Error, Legendre, Hypergeometric, Mathieu, Wave, and other functions
SciPy Cookbook: BPSK: A simple BPSK example with AWGN (no coding)
Matplotlib Gallery: Example plots
This presentation courtesy of Python
jupyter nbconvert --to slides Process.ipynb --post serve
[NbConvertApp] Converting notebook Process.ipynb to slides
[NbConvertApp] Writing 341651 bytes to Process.slides.html
[NbConvertApp] Redirecting reveal.js requests to https://cdn.jsdelivr.net/reveal.js/2.6.2
Serving your slides at http://127.0.0.1:8000/Process.slides.html
Use Control-C to stop this server
WARNING:tornado.access:404 GET /custom.css (127.0.0.1) 4.00ms
WARNING:tornado.access:404 GET /custom.css (127.0.0.1) 3.00ms
In [77]:
%%file slides.bat
jupyter nbconvert --to slides Process.ipynb --post serve