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
   

   View Slide

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
   

   View Slide

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
   

   View Slide

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)
   

   View Slide

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

   View Slide

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

   View Slide

7. Better Presentations on learning
   SymPy Available:
   SciPy 2014 Tutorial
   

   View Slide

8. Not what this talk is about…
   

   View Slide

9. WHY USE
   SYMPY???
   

   View Slide

10. WHY USE
    SYMPY???
    I
    

    View Slide

11. Expressions Can Get Big…
    • Inverted Pendulum, 1 Link
    

    View Slide

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

    View Slide

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

    View Slide

14. DERIVE
    EQUATIONS
    IN SYMPY
    Solution:
    

    View Slide

15. github.com/pydy/pydy
    

    View Slide

16. Added Bonus: Code Generation
    MAGIC
    Symbolic
    Expression
    Numeric
    Code
    

    View Slide

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

    View Slide

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

    View Slide

19. Operations are just Tree Manipulations
    

    View Slide

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)
    

    View Slide

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

    View Slide

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

    View Slide

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)
    

    View Slide

24. BUILDING A
    MATHEMATICS
    COMPILER
    

    View Slide

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

    View Slide

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

    View Slide

27. Need scoped code printer
    • Easy with Visitor Pattern
    

    View Slide

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
    

    View Slide

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)
    

    View Slide

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

    View Slide

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

    View Slide

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
    

    View Slide

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

    View Slide

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)
    

    View Slide

35. Questions???
    

    View Slide