FAST AND CORRECT AND GENERATING CODE WITH SYMPY Jim Crist ME@UMN // github.com/jcrist // jcrist.github.io

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

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

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)

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)]

Submodules (just a few…) • Calculus • Linear Algebra • Set Theory • Combinatorics • Differential Geometry • Tensors • Physics • etc…

Better Presentations on learning SymPy Available: SciPy 2014 Tutorial

Not what this talk is about…

WHY USE SYMPY???

WHY USE SYMPY??? I

Expressions Can Get Big… • Inverted Pendulum, 1 Link

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

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

DERIVE EQUATIONS IN SYMPY Solution:

github.com/pydy/pydy

Added Bonus: Code Generation MAGIC Symbolic Expression Numeric Code

SymPy Expression Trees >>> Lambda((a, b, c), sin(a) + cos(b**2)*c) Created with Graphviz

SymPy Expression Trees Two Rules: 1. .expr.args holds node children 2. “args invariant”* >>> expr.func(*expr.args) == expr *Some exceptions…

Operations are just Tree Manipulations

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)

Python Abstract Syntax Trees >>> lambda a, b, c: sin(a) + cos(b**2)*c Created with Graphviz

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

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)

BUILDING A MATHEMATICS COMPILER

Need new tree elements • Routine • Argument • DataType • Assignment

Need transformations • Type Inferral (assumption system) • Indexed Operations → For Loops • Assignment Propagation

Need scoped code printer • Easy with Visitor Pattern

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

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)

Printers Symbolic Expression Numeric Code Generators Wrappers • C • Fortran • Javascript • Mathematica • Theano • C • Fortran • Cython • F2Py • Numpy ufuncs

Demo (try it at github.com/jcrist/codegen_talk)

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

Correct • No silly math/typing errors • Program declaratively • Validate expressions symbolically before compilation

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)

Questions???