Pro Yearly is on sale from $80 to $50! »

Generating Fast and Correct Code with SymPy

85bba1ca66eb909a289448a90e88f53a?s=47 Jim Crist
September 11, 2014

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

85bba1ca66eb909a289448a90e88f53a?s=128

Jim Crist

September 11, 2014
Tweet

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