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

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…
7. Better Presentations on learning SymPy Available: SciPy 2014 Tutorial
8. Not what this talk is about…
9. WHY USE SYMPY???
10. WHY USE SYMPY??? I
11. Expressions Can Get Big… • Inverted Pendulum, 1 Link
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
14. DERIVE EQUATIONS IN SYMPY Solution:
15. github.com/pydy/pydy
16. Added Bonus: Code Generation MAGIC Symbolic Expression Numeric Code
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…
19. Operations are just Tree Manipulations
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)
24. BUILDING A MATHEMATICS COMPILER
25. Need new tree elements • Routine • Argument • DataType • Assignment
26. Need transformations • Type Inferral (assumption system) • Indexed Operations → For Loops • Assignment Propagation
27. Need scoped code printer • Easy with Visitor Pattern
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
31. Demo (try it at github.com/jcrist/codegen_talk)
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)
35. Questions???