Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Generating Fast and Correct Code with SymPy

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

Jim Crist

September 11, 2014
Tweet

More Decks by Jim Crist

Other Decks in Programming

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. SymPy Expression Trees Two Rules: 1. .expr.args holds node children

    2. “args invariant”* >>> expr.func(*expr.args) == expr *Some exceptions…
  8. 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)
  9. Python Abstract Syntax Trees >>> lambda a, b, c: sin(a)

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

    Fission • Peephole Optimizations • Pattern matching and replacement • Function Inlining • Inline code inside larger function
  11. 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)
  12. 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
  13. 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)
  14. Printers Symbolic Expression Numeric Code Generators Wrappers • C •

    Fortran • Javascript • Mathematica • Theano • C • Fortran • Cython • F2Py • Numpy ufuncs
  15. 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
  16. Correct • No silly math/typing errors • Program declaratively •

    Validate expressions symbolically before compilation
  17. 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)