GSoC 2019 Application Ankit Pandey: Extending Codegen - sympy/sympy GitHub Wiki

I’m a third year Computer Science and Mathematics major from Kathmandu, Nepal currently studying at Grinnell College in Grinnell, Iowa. My programming environment is Emacs on Debian GNU/Linux, as I prefer Emacs for its amazing plugins (such as magit and helm) and customizability through Emacs Lisp. I also consider myself to be generally familiar with git.

I’ve been using Python as a general purpose scripting language for about 5 years now. I’ve occasionally used SymPy over the past 2 years, though most of my experience with computer algebra systems has been with proprietary software like Maple. The most advanced Python feature I have used is async and await, as they require a very different way of reasoning about how programs are executed.

Outside of code generation, my favorite feature of of sympy is nsimplify. For example, if I define a sequence such that \(x_1 = 1\) and

\[ x_{m + 1} = \frac{1}{m + 3/2} ∑_{k = 1}^m x_k x_{m + 1 - k} \]

but I have no idea how to find the closed form of

\[ lim_{m → ∞}\frac{x_m}{x_{m + 1}} \]

I can get an idea of what it might be by manually computing the ratio for a large value of \(m\) and plugging it into nsimplify:

>>> from sympy.abc import x, k, m
>>> from sympy import nsimplify
>>> eq = 1 / (m + Rational(3, 2)) * Sum(x[k] * x[m + 1 - k], (k, 1, m))
>>> f = lambdify([x, m], eq)
>>> def nth_elm(m):
...    x = [0] * (m + 1) # lambdify assumes array indices starts at 1
...    x[1] = 1
...    for j in range(2, m + 1):
...        x[j] = f(x, j - 1)
...    return x[-1]
>>> nsimplify(nth_elm(50) / nth_elm(51), [pi])
  2
pi
---
 2

I chose SymPy’s code generation because it is both an open source project that I was familiar with and aligned well with my interest in programming languages and computer algebra systems.

In the summer of 2018, I worked on an undergraduate research project to generate working Haskell functions from type signatures, which in very general terms also involved generating code from mathematical expressions (though the actual methods of accomplishing this were very different from those that SymPy uses). I also thoroughly enjoyed the courses I’ve taken in programming language implementation, and I think that these experiences give me the familiarity with how programming languages are constructed I need for this project.

Some of the math courses that I have taken required the extensive use of computer algebra systems. As part of one of my projects for these courses, I had to manually port C code from expressions generated by Maple, a large part of which in retrospect could have probably been automated with SymPy’s codegen. While the current real-world usage of codegen seems to be more advanced than what you would find in an undergraduate level project, I think that this gives me basic familiarity with how a CAS and code from other programming languages may be used in tandem.

So far, I’ve opened the following pull requests on SymPy’s Github. I’ve tried to keep them relevant to SymPy’s codegen and printing:

• (Merged) #16331 Make var names in generated Cython code unique
• (Open) #16503: Adjust alignment of the Sum pretty printer
• (Open) #16369: Implicit output for matrix power, mult, and addition in codegen

The purpose of this project is to extend SymPy’s ability to generate code involving matrix expressions by allowing the transformation of SymPy’s AST before generation. Through extensions to the codegen AST, this interface will allow SymPy to be easily extended to generate highly optimized library calls (such as to BLAS/LAPACK).

A second portion of this project is to integrate SymPy with the LFortran project by generating intermediate representations understandable by the tool. Code generation for Fortran can then be offloaded to the LFortran backend.

By the end of this project, SymPy will be able to consistently generate code involving matrix expressions, either as native expressions in the target language or, with proper implementation, external library calls. SymPy will also be able to generate Fortran code through LFortran, enabling further opportunities for interactive running/interpretation of Fortran code.

SymPy is capable of outputting code in a target programming language that is able to numerically evaluate a given expression. This functionality is primarily exposed through the code generation classes given in sympy.printing.codegen. SymPy additionally provides convenience functions that do not require instantiation of these classes.

While SymPy’s code generation facilities are sufficient for handling basic matrix expressions in some languages, the printers are unable to handle all expressions across all languages, either producing invalid expressions in the target language or failing entirely.

As a basic example, suppose that we want to generate C code using ccode to to generate a function that multiplies two matrices, storing it into a new matrix variable. This works easily enough for matrices of a known size, based on the assumption that the given variables have already been initialized with an appropriate size:

>>> from sympy.codegen import *
>>> A = MatrixSymbol('A', 2, 2)
>>> B = MatrixSymbol('B', 2, 2)
>>> X = MatrixSymbol('X', 2, 2)
>>> print(ccode(Assignment(X, A * B)))
X[0] = A[0]*B[0] + A[1]*B[2];
X[1] = A[0]*B[1] + A[1]*B[3];
X[2] = A[2]*B[0] + A[3]*B[2];
X[3] = A[2]*B[1] + A[3]*B[3];

However, the C code generator incorrectly renders more complex matrix expressions, such as transpositions and inverses. The generated expression will not compile, as the Python syntax is directly printed:

>>> A = MatrixSymbol('A', 2, 2)
>>> B = MatrixSymbol('B', 2, 2)
>>> X = MatrixSymbol('X', 2, 2)
>>> print(ccode(Assignment(X, (A.T).I * B)))
X[0] = (A.T**(-1))[0]*B[0] + (A.T**(-1))[1]*B[2];
X[1] = (A.T**(-1))[0]*B[1] + (A.T**(-1))[1]*B[3];
X[2] = (A.T**(-1))[2]*B[0] + (A.T**(-1))[3]*B[2];
X[3] = (A.T**(-1))[2]*B[1] + (A.T**(-1))[3]*B[3];

The same issue also occurs in target languages that support built-in matrix operations, such as Octave/MATLAB:

>>> print(octave_code(Assignment(X, (A.T).I * B)))
X = A.T**(-1)*B;

Extending matrix printing support for languages that support native matrix operations is a relatively simple task, as there is a clear correspondence between the SymPy AST and the target expression (in the above case, the output should be inv(A') * B). This is just a matter of extending the code generation classes to for the appropriate matrix expressions for each language.

Solving this same issue for a lower-level language would require greater involvement. Suppose that we want to generate a C expression for an operation involving matrices that have a variable shape:

>>> from sympy.abc import n, m
>>> A = MatrixSymbol('A', n, m)
>>> B = MatrixSymbol('B', m, n)
>>> X = MatrixSymbol('X', n, n)
>>> print(ccode(Assignment(X, A * B)))
TypeError: 'Symbol' object cannot be interpreted as an integer

By comparison, the MATLAB/Octave code printer is able to represent this assignment with a single expression, as MATLAB/Octave internally handles the size of the matrix.

>>> print(octave_code(Assignment(X, A * B)))
X = A*B;

However, the matrix expression may be transformed into an expression that valid code may be generated from by utilizing the code generation extensions introduced by Björn Dahlgreen. If we rewrite the product of the matrices in terms of IndexedBases, we can produce working code (however, the generated C code is not a single expression, i.e. it cannot be embedded in an if statement):

>>> n = symbols('n', integer=True)
>>> A, B, X = IndexedBase('A'), IndexedBase('B'), IndexedBase('X')
>>> i, j, k = Idx('i', n), Idx('j', n), Idx('k', n)
>>> print(ccode(Assignment(X[i, j], A[i, k] * B[k, j])))
for (int i=0; i<n; i++){
   for (int j=0; j<n; j++){
      X[n*i + j] = 0;
   }
}
for (int i=0; i<n; i++){
   for (int j=0; j<n; j++){
      for (int k=0; k<n; k++){
         X[n*i + j] = A[n*i + k]*B[n*k + j] + X[n*i + j];
      }
   }
}

Here’s a naive implementation of a function that transforms a two-dimensional matrix multiplication into an expression using IndexedBase:

def matmul_translate(e):
    # Only works for an assignment of the form X := A * B, i.e. the RHS must be
    # composed by multiplying two matrices. The shapes must also be symbols of
    # type "integer".
    X_tens = IndexedBase(e.lhs.name)
    A, B = e.rhs.args
    A_tens, B_tens = IndexedBase(A.name), IndexedBase(B.name)
    m = A.shape[0]              # Rows in A
    n = A.shape[1]              # Columns in A, rows in B
    p = B.shape[1]              # Columns in B
    i, j, k = Idx('i', m), Idx('j', n), Idx('p')
    return Assignment(X_tens[i, j], A_tens[i, k] * B_tens[k, j])

The above function may be generalized to multiple matrix multiplications by pre-conversion into multiple assignments inside a single CodeBlock. The translator would however also need to address issues of memory management (how to allocate storage for intermediate values), whether elements are row- or column-major, and the order of evaluation that minimizes memory usage. For example, in an expression such as \(A(BC)D\), if \(B\) has a shape of \(1000 × 1\) and \(C\) has a shape of \(1 × 1000\), we would end up with a \(1000 × 1000\) matrix as an intermediate value. Ideally, our code generator would be able to calculate the order based on known values to find out which order is the optimal.

The fact that we can use a pre-defined set of rules to transform an expression containing a matrix expression into one that the code printer supports suggests that we can extend the code printers to support transformations on the abstract syntax tree. Transformations have the additional advantage that they may be generalized across target languages – both the C and Fortran code printers may use the same transformers for matrix expressions, allowing for easier extensibility and maintenance.

The codegen module may be extended to include a base Transformation class. This class may follow the same visitor pattern as the CodePrinter class, with the only difference being that each method outputs a modified AST node instead of a string. For example, a MatrixTransformation class, which would extend the Transformation class, might have a method (something like _transform_MatMult) that transforms a matrix multiplication expression into an appropriate expression that a lower-level language might be able to handle.

The base code generation class may be extended to include a set of transformation objects that are to be applied to each input expression. The language-specific code printers may then specify the appropriate set of transformers to be applied. As an example, both the C and Fortran code generation classes may have an instance of the MatrixTransformation class as one of the transformers to be applied, since these languages require low-level matrix manipulation. To account for language-specific details, methods in these two instances may be overridden if needed.

A matrix expression may be transformed into an appropriate intermediate representation of an operation upon an array instead of the direct translation discussed in the previous section. The required representations for this already exists within codegen.array_utils and is actively being developed by Francesco Bonazzi (see, for example, #15427) though these are currently not utilized by any of the code printers. Converting these intermediate representations to a final form usable by the code generators might be delegated to language-specific transformers.

Implementing AST transformations in this manner could allow for the extension of the code printers for additional unsupported features, since an additional module could be defined which could have descriptions of how the code within it should be generated. A possible solution to #12172 utilizing this feature would be a transformation class that changes a falling factorial to an equivalent gamma expression.

We would also have the ability to use and extend the algorithms provided in codegen.algorithms, in “default” implementations of functions for more advanced expressions. For example, an integral could be translated into either a function definition utilizing numerical integration algorithms defined in codegen.algorithms and an associated function call or an external library call depending on the user’s preference.

The Transformation class may provide a set of generic transformation and traversal operations to enable easy implementation of transformation rules. An incredibly simple example would be insertion to the top-most level of the AST, which might be implemented as a method called with signature _insert_match(e1, e2). Future implementations of autowrap based on the new sympy.codegen module (a feature which, based on discussions, is eventually planned) might use this function to insert implicit arguments, an issue which I addressed in #1639.

Further possible generic transformations might include:

• Traversal
• Deletion/Insertion
• Pattern matching

While simpler matrix operations would be sufficient for most uses such as generating matrix operations that are known to be for a small size, we may also generate code as calls to external libraries that perform operations on matrices.

The AST must be extended to support a more advanced representation of language antagonistic function and library calls. These extensions might include how variables are passed in (whether through values or reference in C, or the intent in Fortran) and the associated types (using Björn’s work in 2017). Implementation of these extensions would allow relatively simple future addition of a module that translates matrix expressions into appropriate library calls.

The Linear Algebra Package (LAPACK) and Basic Linear Algebra Subprograms (BLAS) are a pair of libraries dedicated to high performance numerical linear algebra. These libraries expose a large number of matrix operations on matrix elements, some of which are specific to assumptions on the input matrix. With the AST extensions for library calls, a basic BLAS/LAPACK generator could be implemented with relative ease.

Matthew Rocklin’s computations presented at SciPy 2013 extends SymPy to generate BLAS and LAPACK calls based on symbolic knowledge of matrix expressions. His dissertation describes it as a method of rewriting of terms, which could be expressed in terms of the transformation API. The addition of LAPACK support would for the most part consist of encoding these translation rules in terms of codegen’s transformation support, enabling a future integration of BLAS/LAPACK and other tools having similar functionality in terms of optimization. BLAS/LAPACK functions also differ based on the types of the contents of matrices, and so type support would most likely be helpful.

LFortran is a modern Fortran compiler that leverages LLVM to compile to modern architectures and enable interactive execution. LFortran internally uses two levels of representation: an Abstract Syntax Tree that represents a collection of Fortran expressions and an Abstract Semantic Representation (ASR) that provides a further level of abstraction by representing semantic aspects of code, such as available variables within a function scope.

Once the existing transformation APIs are in place, we can use them to translate a SymPy expression into an LFortran ASR. Direct translation of a SymPy AST to an LFortran ASR will enable generation of Fortran code through LFortran instead of SymPy, allowing optimizations present in the LFortran project to be added. Additionally, using LFortran’s ASR as an intermediate step would greatly simplify the translation aspect as well as testing and maintenance, since complex expressions could be constructed by LFortran’s ASR builder. This would also eliminate the possibility of bugs in generated Fortran code, since LFortran’s ASR requires the program to be semantically valid. Addtionally, the actual syntax-heavy portions of generation could be handled by LFortran, simplifying testing.

As a simple example, suppose that we want to generate a working Fortran function from the expression \(2x^2 + 11\). We can manually build the ASR representation of a function that computes this using the Python REPL (adapted from one of the builder tests of LFortran):

>>> from lfortran.asr.builder import *
>>> from lfortran.asr.asr import *
>>> tr = make_translation_unit()
>>> m = translation_unit_make_module(tr, "mod1")
>>> integer = make_type_integer()
>>> x = Variable(name="a", intent="in", type=integer) # Our input var
>>> ret = Variable(name="ret", type=integer) # Our return var
>>> expr = BinOp(BinOp(BinOp(x, Mul(), x, integer), Mul(), Num(n=2,
...                                                            type=integer),
...                    integer), Add(),
...              Num(n=17,
...                  type=integer),
...              integer) # The actual expression
>>> assign = asr.Assignment(ret, expr) # ret := 2x^2 + 11
>>> f = scope_add_function(m.symtab, "f", args=[x], return_var=ret)
>>> f.body.extend([assign])
>>> from lfortran.asr.asr_check import verify_asr
>>> verify_asr(tr) # We can check the semantic validity of an expression
>>> print(ast_to_src(asr_to_ast.asr_to_ast(tr)))
module mod1
implicit none


contains

integer function f(a) result(ret)
integer, intent(in) :: a
ret = ((a) * (a)) * (2) + 17
end function

end module

The added layer of abstraction that the ASR provides us means that we can insert variables in the symbol table without worrying about their explicit declarations, allowing us to focus on the expression itself and the associated types. There is also a clear correspondence between the LFortran ASR representation of the expression and its associated SymPy expressions in the SymPy AST, meaning transformation would be relatively straightforward for most expressions that can be represented as a single expression in Fortran.

I intend to work on this project full-time starting May 23, and will be able to invest a full 40 hours a week as I will have no other commitments during the period of the program. I might have to travel during which I will be unavailable for 1-3 days. If I do need to travel, I will give notice at least a month in advance and will be able to make up for the lost time appropriately.

I’d really like to get input on the design of the code generation from other contributors who have worked on this previously (if they are still available), who might include Jim Crist, who extended autowrap to matrix expressions in 2014 and Björn Dahlgreen, who worked on the most recent project extending code generation, and other others who might regularly be using these features. There is also ongoing work on a intermediate expressions for array expressions which I would like to find more about. Since the semester ends on May 22 for me after a full week of exams, I might not be available entirely during this period.

• Add matrix operation printers for languages that natively support matrix expressions such as numpy, MATLAB/Octave and Julia
• Expand the documentation of the existing codegen-specific classes, such as undocumented AST nodes (Assignment, Variable, etc.)
• File first pull request

• Add base classes for allowing the transformation of the abstract syntax tree
• Implement basic transformation API functions

• Extend the AST to include functions to perform an external library call
  • These will probably be designed with the intention of future integration of BLAS/LAPACK
• Add associated tests for transformations
• File second pull request

• Implement transformers of matrix expressions to intermediate expressions.
  • I suspect this will take a great amount of time because of complexities associated with optimizing matrix expressions, in addition to implementing all expressions, such as transposition and inverses.
  • Additional optimizations (though not on the level of BLAS) might be possible based on assumptions
• Refine the Transformation API
  • Since this is the first time the transformers will be used, issues might become apparent that would require this interface to be refined.

• Add code printers for the final output of intermediate representations, focusing on languages with non-native support for matrices
• Consideration of any additional features, such as memory management
• File third pull request

See the related issue on Gitlab.

• Basics of a SymPy-LFortran ASR translator
  • The basis for this might be the existing Fortran code transformer, and would cover translations into LFortran objects that are direct and straightforward.

• More advanced portions of SymPy-LFortran ASR translator
  • This might include changes that take advantage of LFortran-specific representations
• Move towards completeness of Fortran code generation
  • The ASR \(→\) AST translator currently does not implement some operations (such as division) that the Fortran code generator is able to handle. These should be implemented for completeness.
• File fourth pull request, in addition an upstream request in LFortran

• Final write-up, work on additional remaining bug fixes and documentation

• Björn Dahlgreen’s final report (2017)
• 2017 Codegen Gitter discussions
• Aaron Muerer’s notes on Codegen (2015)
• Jim Crist’s work on code printing (2014)
• Oyvind Jensen’s Report (2010)
• Developer Tutorial - LFortran Documentation
• SymPy integration Issue on LFortran