I am going through the various code generation options for SymPy. Here are some notes on them.

Note that I plan to improve the code generation, so this document may become out-of-date (I am writing this on October, 2015). Everything here is based on the development version of SymPy (unreleased 1.0).

Here is a link to the relevant mailing list discussion https://groups.google.com/d/msg/sympy/omYuHN68-_U/-hrGmsUgBwAJ.

-- Aaron Meurer

There are several submodules in SymPy relating to code generation, or otherwise compiling SymPy expressions down to something that can be used to numerically evaluate them fast.

• Created by Øyvind Jensen for his 2010 GSoC project.

• Main functions are autowrap, binary_function, and ufuncify.

• All automatically generates code, writes it to disk, compiles it, and imports it into the current session.

• It also automatically converts expressions containing Indexed objects into summations (tensors).

• autowrap creates a wrapper using f2py or Cython and creates a numerical function.

• binary_function is similar to autowrap except the it creates a symbolic function, which is evaluated numerically using evalf.

• ufuncify creates a numpy ufunc. The f2py and Cython backends create "fake" ufuncs (no type conversion; only works on equal length 1-dim arrays). The numpy backend (default) creates a real NumPy ufunc.

This module deals with creating compilable code from SymPy expressions. This is lower level than autowrap, as it doesn't actually attempt to compile the code, but higher level than the printers, as it generates compilable files (including header files), rather than just code snippets.

• The main entry points here are codegen and make_routine.

• codegen takes a list of (variable, expression) pairs and a language (C, F95, and Octave/Matlab are supported). It returns, as strings, a code file and a header file (for relevant languages). The variables are created as functions that return the value of the expression as output.

  Example:

  >>> [(c_name, c_code), (h_name, c_header)] = codegen([('y', sin(x))], 'C')
  >>> print(c_name)
  y.c
  >>> print(c_code)
  /******************************************************************************
   *                    Code generated with sympy 0.7.7.dev                     *
   *                                                                            *
   *              See http://www.sympy.org/ for more information.               *
   *                                                                            *
   *                       This file is part of 'project'                       *
   ******************************************************************************/
  #include "y.h"
  #include <math.h>
  
  double y(double x) {
  
     double y_result;
     y_result = sin(x);
     return y_result;
  
  }
  
  >>> print(h_name)
  y.h
  >>> print(c_header)
  /******************************************************************************
   *                    Code generated with sympy 0.7.7.dev                     *
   *                                                                            *
   *              See http://www.sympy.org/ for more information.               *
   *                                                                            *
   *                       This file is part of 'project'                       *
   ******************************************************************************/
  
  
  #ifndef PROJECT__Y__H
  #define PROJECT__Y__H
  
  double y(double x);
  
  #endif

• Various flags to codegen let you modify things like the project name and global variables.

• You can also subclass the CodeGen classes to customize this further.

• make_routine creates a Routine object, which represents an evaluation routine for a set of expressions. As far as I can tell, this is only good for internal use by the CodeGen objects, as an intermediate representation from SymPy expression to generated code.

Routine seems like it would be better named Function. It requires a name and a list of input arguments and return values (perhaps the naming was inspired by Fortran). You generally don't want to try to make a Routine object yourself. You should instead use make_routine, which will call the routine method of the relevant CodeGen object for your language. This creates the internal objects representing assignments and so on, and creates the Routine class with them.

The various codegen objects such as Routine and Variable aren't SymPy objects (they don't subclass from Basic).

• This is where the meat of code generation is, the translation of SymPy expressions to specific code.

• Supported languages are C (ccode), Fortran <95 and >95 (fcode), Javascript (jscode), Mathematica (mathematica_code), Octave/Matlab (octave_code), Python (print_python, which is actually more like a lightweight version of codegen for Python, and lambdarepr, which supports many libraries (like NumPy), and theano, if that counts (theano_code).

• The code printers are special cases of the other prints in SymPy (str printer, pretty printer, etc.).

• An important distinction is that the code printer has to deal with assignments (using the sympy.printing.codeprinter.Assignment object).

• The various printers also tend to support Indexed objects well.

There are three levels

autowrap
    |
codegen
    |
code printers

Autowrap uses codegen, and codegen uses the code printers. Autowrap does everything: it lets you go from SymPy expression to numerical function in the same Python process in one step. codegen is actual code generation, i.e., to compile and use later, or to include in some larger project.

The code printers translate the SymPy objects into actual code, like abs(x) -> fabs(x) (for C).

In terms of potential optimizations, there is a lot that can be done at all levels.

The code printers don't print optimal code in many cases. An example of this is powers in C. x**2 prints as pow(x, 2) instead of x*x. Other optimizations (like mathematical simplifications) should happen before the code printers.

cse is not applied anywhere in this chain. It ideally happen at the codegen level, or somewhere above it.

At the autowrap level could go more knowledge of compiler flags, (autowrap currently defers this knowledge to f2py). Autowrap isn't always the best for every use-case, and the automatic summation can be annoying. Perhaps another module that lives at the same level as autowrap, but interfaces more directly with compiles (and maybe isn't intended directly for generating a Python module) is in order.

autowrap is also annoying in that f2py itself tends to have issues.

https://github.com/paralab/SymPyGR

This is a weird little submodule. It seems to date back to 2008. It isn't included in the online documentation, and I'm not sure if anyone uses it. It uses ctypes and tcc to generate fast numerical functions. It's similar to autowrap, except it seems to be mostly self-contained (it uses lambdify and a 16 line heuristic function to generate the C code).

The module only works in Linux (I believe tcc itself only works on Linux), so I was not able to test it.

It hasn't been updated at all since it was written (except for some code cleanups). I have created a pull request to remove it https://github.com/sympy/sympy/pull/10018.

There are a lot of external tools that use SymPy to do code generation. Some common themes that I've seen are:

• Most tools are domain specific. This is to be expected since the end goal is to solve a real problem, but, with a couple of exceptions, most tools build out some layers of code generation that could be generic without any separation from the domain specific parts.

• Given the three layers outlined above (autowrap, codegen, code printers), most tools work at the second layer. However, nothing I've found so far actually uses the codegen module, indicating that this module is unknown, doesn't solve people's problems, or both. Instead, what I see is tools writing their own ad-hoc version of this layer. This generally involves some kind of templating, or otherwise patching together code to create a final C/Python/etc. module.

• The exception is Björn Dahlgren's code (https://github.com/bjodah and https://github.com/chemreac), which uses his pycodeexport module. This module is basically an alternative to SymPy's codegen. It also uses templates, but is written in a generic way (separated from domain specific code), and uses a real templating system (Mako). To quote the source code:

  The approach taken here differs from sympy.utilities.codegen.codegen through the use of Mako templates and classes which controlls the how they are rendered. Which method is best depends on the problem at hand (personal opinion).

To list out some of the stuff I've seen:

https://github.com/pydy/pydy

A unique example, in that PyDy makes as much use of stuff that is already implemented in SymPy (although it does not use codegen). This is because the PyDy developers implement things in SymPy itself when relevant. With that being said, PyDy does also have its own version of autowrap (pydy.codegen.cython_code).

The PyDy code is mostly focused on generating code for matrices (systems of ODEs), which is the primary reason (I believe) why code exists in PyDy that does things that SymPy itself would otherwise do. The PyDy code generators use predefined templates specifically written for this. It does make full use of the code printers (even subclassing the printer for some special case code).

https://github.com/csu-hmc/opty

Also has a custom code generator based on SymPy and PyDy for evaluation of big matrices in a tight loop.

https://github.com/cdsousa/SymPyBotics

The latest cse improvements were motivated by this use case. Similar to PyDy.

https://github.com/chemreac

(Björn says: Actually, chemreac does not use SymPy for code generation, it uses SymPy in its unittests as a sanity check, and it uses pycodeexport as an alternative pre-processor for the C++ code, the latter is probably suboptimal design but it got the job done, a better example for templated codegeneration is then: fastinverse, you might want to look here for a minimal example)

Björn Dahlgren's library (this along with some other code of his). This makes full use of his pycodeexport library, as an alternative to codegen. This means large C files with templates (example: https://github.com/chemreac/chemreac/blob/master/src/chemreac_template.cpp). This is definitely easier to read than spread out code with lots of lines.append(stuff) and code.format(**vars), but I'm not sure if it's easier to maintain. For instance, a codegen-type approach (or at least what I'd envision; I haven't yet reviewed codegen itself in depth) would generate an AST (using SymPy expressions or otherwise) of the generated code. The advantage here is that the AST itself could be simplified or modified before dumping it to code.

• The code generation itself is done in pycodeexport and uses the code printers.

• Question: I see that pycodeexport supports using cse, but I can't tell how it works, and I don't see any examples that use it. (Björn says: -It's not too well documented I'm afraid, here's one example)

• Final note: Björn has another library, pycompilation, which might be considered an alternative to autowrap (except, AFAICT, with more of a focus on reusability). It appears to be a smarter distutils with better support for things that use code generation (it does not use code generation directly itself; pycodeexport uses pycompilation, not the other way around). (Björn says: yes, I was unable to cache object files with distutils so I wrote pycompilation which I also gave some convenience functions for abstracting away compilation and import of the resulting python module).

Just a little academic code from a paper I found (http://cpc.cs.qub.ac.uk/summaries/AESD_v1_0.html). You can find the paper online but unfortunately you have to have a journal subscription to download the code.

The code itself isn't that great. The code generation is to Python/NumPy/mpmath (rather than C). The code doesn't make use of existing stuff in SymPy (like lambdarepr; it achieves the equivalent of lambdarepr with a Python header with import as lines). The architecture is not terrible, but not great either (basically writing out a template file line-by-line). The code quality is quite low; it is littered with anti-patterns (trust me, you aren't missing much by not being able to access this code). There is virtually no separation of concerns between the generation and the domain specific knowledge in this code.

The main thing I got from this, aside from being reminded by how bad academic code can be, is a reminder that horner is a useful optimization that we can apply. I suppose it also reinforces the importance that SymPy's str representation be valid Python to recreate the expression, since people don't know that things like lambdarepr exist.

Anthony Scopatz's (and others') nuclear engineering code. The two examples are https://github.com/pyne/pyne/blob/develop/src/decaygen.py and https://github.com/pyne/pyne/blob/develop/pyne/apigen/enrich_multi_sym.py.

The first, decaygen, doesn't actually use SymPy, but shows some of the stuff SymPy could do in terms of optimization. There is a paper on this ("Binary Formulation of Decay Equations", Anthony M. Scopatz, Cameron R. Bates, Paul P. H. Wilson). The idea is that they sped up some existing code for a problem with a closed-form solution by precomputing constants and using exp2 instead of exp. The module does do code generation, but does not use SymPy since the problem is static and does not require much mathematical derivation.

The second, enrich_multi_sym, also has a paper ("First & Second Order Approximations to Stage Numbers in Multicomponent Enrichment Cascades", Anthony Scopatz). Here, a problem is solved that doesn't have a closed-form solution, but a second order approximation is used, which does have a closed-form solution.

SymPy is used, but a lot of things are done manually (e.g., the second order approximation is computed directly (second derivative divided by 2), and solved by hand (the quadradic equation). Also, the fact that the approximation is second order is hard-coded (first order was not accurate enough, but second order was). The fact that the quadradic equation is hard-coded is partly due to the fact that a particular form of it is required for numerical accuracy.

The code also makes use of SymPy's cse() to reduce the code (this is apparently needed because the compiler cse is too slow, requiring the use of -O0).

The code is all templated into hard-coded C header, footer, etc.

One note here: this code does use codegen (in pyne.apigen.utils.cse_to_c), but it really should be using ccode. It generates code from codegen and uses regular expressions to extract it. The code itself does not even work on the latest version of SymPy because of changes in the codegen output.

Theano is a Python library that allows you to define, optimize, and evaluate mathematical expressions involving multi-dimensional arrays efficiently. It can use GPUs and perform efficient symbolic differentiation. It is heavily used in machine learning.

https://github.com/IgnitionProject/ignition.

A library from Andy Terrel designed around writing DSLs to generate code. It uses SymPy.

There is a paper "From equations to code: Automated Scientific Computing" (Andy Terrel, 2011) that describes int_gen, one of the submodules of Ignition.

Looking at the code for ignition, I doubt it works any more (with the current SymPy).

It seems to focus more on the DSL aspect of code generation. It does use ccode, but implements its own version of the codegen layer. This appears to use an AST + visitor pattern approach (the first example I've seen that does). See ignition.code_tools. However (if I'm reading the code correctly), not everything uses it (the sfl module does, but int_gen doesn't seem to).

This has some ideas for optimizations (like precomputing constants), but not many ideas for an optimization pipeline. Each DSL hard-codes its own optimizations (e.g., int_gen knows that it is dealing with integrals, and precomputes constant integrals and evaluates constants, and has a lot of logic around quadrature).

https://github.com/cosmo-ethz/hope

This library lets you add a @jit decorator to functions, and it parses the ast and compiles it to C++. It converts to and from SymPy as an intermediary representation so that it can use simplify and cse for optimization.

It doesn't look as fleshed out as numba, although the fact that they do use SymPy to do some optimizations is interesting.

The pipeline appears to be

Python Code -> Python AST -> SymPy expression -> sympy.simplify()/sympy.cse() -> Python AST -> C++ -> Compile an extension and import

In particular, they go to and from SymPy just long enough to use simplify and cse, and then go right back to Python AST (they do not use any of the code printing modules in SymPy).

Obviously, I'd prefer to see it use more SymPy, but it's an interesting idea. Like virtually every other example, though, the optimizations are hard-coded and cannot be modified (only optionally turned completely off).

Matt Turk, et al. A meta-solver for chemical reaction networks and cooling processes. Makes heavy use of SymPy code generation (ccode, again, no codegen). It fills out predefined C template files using jinja2 templates.

I don't know if there is any sort of symbolic optimization happening here. I didn't see anything, but I also know nothing about the domain.

https://github.com/casadi/casadi

Casadi uses symbolic graphs to generate fast numerical optimization routines.

A code generation packages for Maple to C with a focus on ODEs. The code is at http://www.cecm.sfu.ca/~wittkopf/ToExternal.html and there is a paper at http://www.jnaiam.org/new/uploads/files/16985fffb53018456cf3506db1c5e42b.pdf. It predates Maple's own Compiler module.

An interesting thing here is that Maple has automatic uniquification, making common subexpression elimination (at least in the trivial sense) simpler, because common subexpressions already share the same memory address.

There are also a few optimizations of note here, which I have added to the Code-Generation-Optimizations page.

https://github.com/tensorflow/tensorflow

A Google create lib that allows you to define mathematical ops with graphs that get translated into computation.

https://github.com/coneoproject/COFFEE

COFFEE is a COmpiler For Fast Expression Evaluation. Used as a backend for solving PDEs.

https://github.com/coneoproject/SLOPE

SLOPE is a run-time system for fusing and tiling sequences of loops with indirect memory accesses (e.g., A[B[i]]). It also provides a python interface to ease the integration of the library with scientific codes based on python and run-time code generation.

http://openrave.org/docs/0.8.2/openravepy/ikfast/

ikfast uses sympy to analytically compute the inverse kinematics of robotic manipulators. It then generates C++ code to compute the inverse kinematics at top speed. It uses SymPy's cse method.

https://github.com/PyLinX-Project/PyLinX

Simulink clone that generates code.

https://github.com/jcrist/symcc

Experimental new code printing architecture for SymPy

http://jitcode.readthedocs.io/en/latest/

Uses sympy to do automatic diff of ode's.

https://github.com/opesci/devito. https://arxiv.org/abs/1608.08658

Automated fast finite difference computation. Uses SymPy.

Uses SymPy to parse model descriptions and then manipulate the expressions, utlimately code generation is done with SciPy's weave. This solves neural network problems.

http://briansimulator.org/