3.2k

# 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
[email protected] // 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')
>>> 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…

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

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

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???