Tweet
Tweet

Slide 1

Slide 1 text

SymPy Code Generation Aaron Meurer Anthony Scopatz ERGS, University of South Carolina SciPy 2016, Austin, TX

Slide 2

Slide 2 text

SymPy • Powerful computer algebra system (CAS) written in pure Python • BSD licensed • Usable as a library • Just released version 1.0

Slide 3

Slide 3 text

No content

Slide 4

Slide 4 text

Code Generation • Translate SymPy expressions to another language • Example: translate |sin(π⋅x)| to C >>> ccode(abs(sin(pi*x))) 'fabs(sin(M_PI*x))'

Slide 5

Slide 5 text

Languages supported • C • Fortran • MATLAB/Octave • Python (NumPy/SciPy) • Julia • Mathematica • Javascript • LLVM • Rust • Theano • Easy to extend to others

Slide 6

Slide 6 text

Code generation workflow Derive mathematical formulas symbolically Convert symbolic expressions to numeric function Use numeric function to solve some problem SymPy Code generation Domain specific

Slide 7

Slide 7 text

Code generation workflow Derive mathematical formulas symbolically Convert symbolic expressions to numeric function Use numeric function to solve some problem SymPy Code generation Domain specific

Slide 8

Slide 8 text

Code generation levels of abstraction Code printers codegen ufuncify lambdify 'fabs(sin(M_PI*x))' #include "f.h" #include double f(double x) { double f_result; f_result = fabs(sin(M_PI*x)); return f_result; } f = ufuncify(x, abs(sin(pi*x))) Expression Larger block of code Python callable

Slide 9

Slide 9 text

Why do code generation?

Slide 10

Slide 10 text

Why do code generation?

Slide 11

Slide 11 text

Why do code generation? 1. SymPy can deal with mathematical expressions in a high-level way. For example, it can take symbolic derivatives.

Slide 12

Slide 12 text

Why do code generation? 1. SymPy can deal with mathematical expressions in a high-level way. For example, it can take symbolic derivatives. 2. Using code generation avoids mistakes that come from translating mathematics into low level code.

Slide 13

Slide 13 text

Why do code generation? 1. SymPy can deal with mathematical expressions in a high-level way. For example, it can take symbolic derivatives. 2. Using code generation avoids mistakes that come from translating mathematics into low level code. 3. It’s possible to deal with expressions that are otherwise too large to write by hand.

Slide 14

Slide 14 text

Why do code generation? 1. SymPy can deal with mathematical expressions in a high-level way. For example, it can take symbolic derivatives. 2. Using code generation avoids mistakes that come from translating mathematics into low level code. 3. It’s possible to deal with expressions that are otherwise too large to write by hand. 4. Some “mathematical” optimizations are possible, which a normal compiler would not be able to do.

Slide 15

Slide 15 text

Example: Iodine-131

Slide 16

Slide 16 text

Iodine-131 • Iodine-131 has a half-life of 8.0197 days • Cells in the thyroid absorb Iodine • Radioactive Iodine-131 destroys thyroid cells by short-range beta radiation

Slide 17

Slide 17 text

• I-131 decays with a half-life of 8.02 days 131I ⟶ β- + 131*Xe 131I ⟶ β- + 131Xe 131*Xe ⟶ 131Xe @xi @t = ixi( t ) + X i6=j j j!ixj( t )

Slide 18

Slide 18 text

No content

Slide 19

Slide 19 text

No content

Slide 20

Slide 20 text

No content

Slide 21

Slide 21 text

No content

Slide 22

Slide 22 text

No content

Slide 23

Slide 23 text

• This is a simple example, because the decay of 131I only results in three species, 131I, 131*Xe, and 131Xe • A typical decay chain may result in hundreds of species • With SymPy, we can avoid mistakes by representing the decay equations in a high level way, and deriving the low level representations

Slide 24

Slide 24 text

Example: n-link pendulum on cart (PyDy) Thanks to Jason Moore

Slide 25

Slide 25 text

PyDy • Multibody dynamics • SymPy is used to derive the equations of motion for a mechanical system • The resulting equations are a large system of nonlinear ODEs which need to be integrated (F=Ma) • Code generation allows creating fast callbacks which can be integrated over many time steps

Slide 26

Slide 26 text

• Needs to be fast because: • Realtime simulation • Optimal control • Stiff systems require more time steps

Slide 27

Slide 27 text

n-link pendulum on a cart

Slide 28

Slide 28 text

No content

Slide 29

Slide 29 text

No content

Slide 30

Slide 30 text

No content

Slide 31

Slide 31 text

• trigsimp() can simplify the equations of motion significantly • In this example, 
 sin(x)⋅cos(y) - sin(y)⋅cos(x) ⟶ sin(x - y) • The equations must be evaluated at each time step, so this can make a significant difference in performance

Slide 32

Slide 32 text

• PyDy automatically generates a fast ODE solve callback using SymPy code generation and Cython

Slide 33

Slide 33 text

No content

Slide 34

Slide 34 text

Animation (code not shown)

Slide 35

Slide 35 text

Animation (code not shown)

Slide 36

Slide 36 text

We can control the pendulum by forcing the cart

Slide 37

Slide 37 text

Apply the force f to keep the pendulum upright

Slide 38

Slide 38 text

No content

Slide 39

Slide 39 text

No content

Slide 40

Slide 40 text

No content

Slide 41

Slide 41 text

Animation

Slide 42

Slide 42 text

Animation

Slide 43

Slide 43 text

Yes this is possible https://www.youtube.com/watch?v=cyN-CRNrb3E

Slide 44

Slide 44 text

We can increase the number of links (n=6)

Slide 45

Slide 45 text

No content

Slide 46

Slide 46 text

No content

Slide 47

Slide 47 text

No content

Slide 48

Slide 48 text

• More details on this problem are at https:// github.com/pydy/pydy/blob/master/examples/ npendulum/n-pendulum-control.ipynb

Slide 49

Slide 49 text

How does it work? • SymPy expressions are stored as trees • Example: x2 + xy Add Mul Pow Symbol(‘x’) 2 Symbol(‘y’) Symbol(‘x’)

Slide 50

Slide 50 text

How does it work? • Every expression stores its children expression in .args >>> (x**2 + x*y).args (x**2, x*y) >>> (x**2 + x*y).args[0].args (x, 2)

Slide 51

Slide 51 text

How does it work? class CCodePrinter(CodePrinter): def _print_Rational(self, expr): p, q = int(expr.p), int(expr.q) return '%d.0L/%d.0L' % (p, q) def _print_Exp1(self, expr): return "M_E" def _print_Pi(self, expr): return 'M_PI'

Slide 52

Slide 52 text

How does it work? • Printer subclasses walk the expression tree and call methods corresponding to children (visitor pattern) • Subclass CodePrinter, define methods for the expression types to code generate • Easy to write your own code printers, or to extend existing code printers to do the things you need

Slide 53

Slide 53 text

Some other libraries that use SymPy code generation • Chemreac • python library for solving chemical kinetics problems with possible diffusion and drift contributions • SymPyBotics • Symbolic Framework for Modeling and Identification of Robot Dynamics

Slide 54

Slide 54 text

Takeaways 1. SymPy can deal with mathematical expressions in a high-level way. For example, it can take symbolic derivatives. 2. Using code generation avoids mistakes that come from translating mathematics into low level code. 3. It’s possible to deal with expressions that are otherwise too large to write by hand. 4. Some “mathematical” optimizations are possible, which a normal compiler would not be able to do.

Slide 55

Slide 55 text

• Mailing list: http://groups.google.com/group/ sympy • https://github.com/sympy/sympy • @asmeurer, @SymPy on Twitter • These slides are at https://github.com/asmeurer/ SciPy-2016-Talk • I’ll be at the sprints (and other SymPy developers)

Slide 56

Slide 56 text

Questions