Jim Crist
September 11, 2014
3.4k

# Generating Fast and Correct Code with SymPy

All about the why and how of code generation in SymPy. Discusses the current state of converting symbolic math into fast numeric code, and gives a brief overview of current and future improvements.

Demo ipython notebook available on github here: https://github.com/jcrist/codegen_talk

## Jim Crist

September 11, 2014

## Transcript

1. ### FAST AND CORRECT AND GENERATING CODE WITH SYMPY Jim Crist

ME@UMN // github.com/jcrist // jcrist.github.io
2. ### What is SymPy??? • A Computer Algebra System (CAS) written

in pure Python • Allows you to create expressions, and do things with them • SymPy knows about all the math needed
3. ### Create Expressions >>> from sympy import * # Create Symbols

>>> a, b, c = symbols('a, b, c') # Combine Symbols into expressions >>> expr = sin(a) + cos(b**2)*c >>> expr sin(a) + cos(b**2)*c
4. ### Manipulate Expressions # subs performs substitutions on expressions >>> expr

= expr.subs(b, a**(1/2)) >>> expr sin(a) + sin(a)*c # simplify reduces the size of expressions >>> expr = expr.simplify() >>> expr sin(a)*(1 + c)
5. ### Solve Expressions # Create another symbol x >>> x =

symbols('x') # Quadratic equation >>> lhs = a*x**2 + b*x + c >>> expr = Eq(lhs, 0) >>> expr a*x**2 + b*x + c == 0 # Solve returns solutions to expressions >>> solve(expr, x) [(-b + sqrt(-4*a*c + b**2))/(2*a), -(b + sqrt(-4*a*c + b**2))/(2*a)]
6. ### Submodules (just a few…) • Calculus • Linear Algebra •

Set Theory • Combinatorics • Differential Geometry • Tensors • Physics • etc…

12. ### Expressions Can Get Big… • Inverted Pendulum, 3 Links: Rendered

in size 2.5 font
13. ### Expressions Can Get Big… • Inverted Pendulum, 3 Links: Rendered

in size 2.5 font

17. ### SymPy Expression Trees >>> Lambda((a, b, c), sin(a) + cos(b**2)*c)

Created with Graphviz
18. ### SymPy Expression Trees Two Rules: 1. .expr.args holds node children

2. “args invariant”* >>> expr.func(*expr.args) == expr *Some exceptions…

20. ### Example: subs def crawl(expr, func, *args, **kwargs): """Crawl the expression

tree, and apply func to every node""" val = func(expr, *args, **kwargs) if val is not None: return val new_args = (crawl(arg, func, *args, **kwargs) for arg in expr.args) return expr.func(*new_args) def sub_func(expr, sub_dict): """Perform direct matching substitution""" if expr in sub_dict: return sub_dict[expr] elif not expr.args: return expr # Compose the subs function subs = lambda x, sub_dict: crawl(x, sub_func, sub_dict)
21. ### Python Abstract Syntax Trees >>> lambda a, b, c: sin(a)

+ cos(b**2)*c Created with Graphviz
22. ### Optimization = Tree Manipulation • Loop Optimizations • Fusion •

Fission • Peephole Optimizations • Pattern matching and replacement • Function Inlining • Inline code inside larger function
23. ### Optimization = Tree Manipulation • Loop Optimizations • Fusion (Tensor

Contraction) • Fission • Peephole Optimizations • Pattern matching and replacement (subs, replace, xreplace) • Function Inlining • Inline code inside larger function (subs)

• Assignment
26. ### Need transformations • Type Inferral (assumption system) • Indexed Operations

→ For Loops • Assignment Propagation

28. ### Visitor Pattern 1. .Printer has methods for each node type

2. .Printer._print dispatches to appropriate print method 3. Recursively call Printer._print for all args in tree
29. ### Example: _print_For class CCodePrinter(CodePrinter): ... def _print_For(self, expr): target =

self._print(expr.target) start, stop, step = expr.iterable.args body = '\n'.join(self._print(i) for i in expr.body) return ('for ({target} = {start}; {target} < {stop}; ' '{target} += {step}) {{\n{body}\n}}').format( target=target, start=start, stop=stop, step=step, body=body)
30. ### Printers Symbolic Expression Numeric Code Generators Wrappers • C •

Fortran • Javascript • Mathematica • Theano • C • Fortran • Cython • F2Py • Numpy ufuncs

32. ### Fast • Compile to native code • Can be smarter

than a normal compiler • Domain knowledge >>> peephole optimizations • Change algorithmic complexity • ufuncs → Broadcasting and type conversion for free
33. ### Correct • No silly math/typing errors • Program declaratively •

Validate expressions symbolically before compilation
34. ### The Future • More Optimizations • Horner’s Method • Common

sub-expression elimination • BLAS/LAPACK method selection • Boolean Support • Expressive Iteration (map, reduce, filter, …) • More Modular • Templating • github.com/jcrist/symcc (shameless plug)