The Harder Way: C Code generation, Custom Printers, and CSE [1 hour]

One of the most common low level programming languages in use is C. Compiled C code can be optimized for execution speed for many different computers. Python is written in C as well as many of the vectorized operations in NumPy and numerical algorithms in SciPy. It is often necessary to translate a complex mathematical expression into C for optimal execution speeds and memory management. In this notebook you will learn how to automatically translate a complex SymPy expression into C, compile the code, and run the program.

We will continue examining the complex chemical kinetic reaction ordinary differential equation introduced in the previous lesson.

Learning Objectives

After this lesson you will be able to:

  • use a code printer class to convert a SymPy expression to compilable C code
  • use an array compatible assignment to print valid C array code
  • subclass the printer class and modify it to provide custom behavior
  • utilize common sub expression elimination to simplify and speed up the code execution

Import SymPy and enable mathematical printing in the Jupyter notebook.

In [1]:
import sympy as sym
In [2]:
sym.init_printing()

Ordinary Differential Equations

The previously generated ordinary differential equations that describe chemical kinetic reactions are loaded below. These expressions describe the right hand side of this mathematical equation:

$$\frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y}(t))$$

where the state vector $\mathbf{y}(t)$ is made up of 14 states, i.e. $\mathbf{y}(t) \in \mathbb{R}^{14}$.

Below the variable rhs_of_odes represents $\mathbf{f}(\mathbf{y}(t))$ and states represents $\mathbf{y}(t)$.

From now own we will simply use $\mathbf{y}$ instead of $\mathbf{y}(t)$ and assume an implicit function of $t$.

In [3]:
from scipy2017codegen.chem import load_large_ode
In [4]:
rhs_of_odes, states = load_large_ode()

Exercise [2 min]

Display the expressions (rhs_of_odes and states), inspect them, and find out their types and dimensions. What are some of the characteristics of the equations (type of mathematical expressions, linear or non-linear, etc)?

Double Click For Solution

In [5]:
# write your solution here

Compute the Jacobian

As has been shown in the previous lesson the Jacobian of the right hand side of the differential equations is often very useful for computations, such as integration and optimization. With:

$$\frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y})$$

the Jacobian is defined as:

$$\mathbf{J}(\mathbf{y}) = \frac{\partial\mathbf{f}(\mathbf{y})}{\partial\mathbf{y}}$$

SymPy can compute the Jacobian of matrix objects with the Matrix.jacobian() method.

Exercise [3 min]

Look up the Jacobian in the SymPy documentation then compute the Jacobian and store the result in the variable jac_of_odes. Inspect the resulting Jacobian for dimensionality, type, and the symbolic form.

Double Click For Solution

In [6]:
# write your answer here

C Code Printing

The two expressions are large and will likely have to be excuted many thousands of times to compute the desired numerical values, so we want them to execute as fast as possible. We can use SymPy to print these expressions as C code.

We will design a double precision C function that evaluates both $\mathbf{f}(\mathbf{y})$ and $\mathbf{J}(\mathbf{y})$ simultaneously given the values of the states $\mathbf{y}$. Below is a basic template for a C program that includes such a function, evaluate_odes(). Our job is to populate the function with the C version of the SymPy expressions.

#include <math.h>
#include <stdio.h>

void evaluate_odes(const double state_vals[14], double rhs_result[14], double jac_result[196])
{
      // We need to fill in the code here using SymPy.
}

int main() {

    // initialize the state vector with some values
    double state_vals[14] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14};
    // create "empty" 1D arrays to hold the results of the computation
    double rhs_result[14];
    double jac_result[196];

    // call the function
    evaluate_odes(state_vals, rhs_result, jac_result);

    // print the computed values to the terminal
    int i;

    printf("The right hand side of the equations evaluates to:\n");
    for (i=0; i < 14; i++) {
        printf("%lf\n", rhs_result[i]);
    }

    printf("\nThe Jacobian evaluates to:\n");
    for (i=0; i < 196; i++) {
        printf("%lf\n", jac_result[i]);
    }

    return 0;
}

Instead of using the ccode convenience function you learned earlier let's use the underlying code printer class to do the printing. This will allow us to modify the class to for custom printing further down.

In [7]:
from sympy.printing.ccode import C99CodePrinter

All printing classes have to be instantiated and then the .doprint() method can be used to print SymPy expressions. Let's try to print the right hand side of the differential equations.

In [8]:
printer = C99CodePrinter()
In [9]:
print(printer.doprint(rhs_of_odes))
// Not supported in C:
// ImmutableDenseMatrix
Matrix([
[                                                                                                                                                                                      -14520000.0*y0**2 - 0.0158*y0*y1 - 20900000.0*y0*y11 - 27600000.0*y0*y4 - 35500000.0*y0*y5 - 13600000.0*y0*y6 - 22900000.0*y0*y7 - 13000000.0*y0*y8 - 13000000.0*y0*y9 + 24400.0*y2*y4 + 5.83*y4 + 0.0854],
[-14520000.0*y0**2 - 0.0158*y0*y1 - 27600000.0*y0*y4 - 26000000.0*y0*y8 - 1270.0*y1*y10 - 1270.0*y1*y12 - 4.58e-8*y1*y4 - 0.000155*y1*y8 - 2.12e-5*y1 + 4070000.0*y10*y5 + 118000000.0*y11*y2 + 24400.0*y2*y4 + 13300000.0*y2*y5 + 13300000.0*y2*y6 + 13300000.0*y2*y9 + 39500.0*y3*y5 + 10900000.0*y4*y5 + 36500.0*y4*y6 + 29200.0*y5*y6 + 8840000.0*y5*y9 + 2.6e-7*y6 - 101000.0*y8*y9 - 0.131],
[                            14520000.0*y0**2 + 0.0158*y0*y1 + 27600000.0*y0*y4 + 35500000.0*y0*y5 + 13600000.0*y0*y6 + 26000000.0*y0*y8 + 1270.0*y1*y10 + 1270.0*y1*y12 + 0.000155*y1*y8 + 2.12e-5*y1 + 786000.0*y10*y12 - 118000000.0*y11*y2 + 128000.0*y12*y3 + 4070000.0*y12*y6 - 24400.0*y2*y4 - 13300000.0*y2*y5 - 13300000.0*y2*y6 - 13300000.0*y2*y9 + 11000000.0*y5*y8 + 101000.0*y8*y9],
[                                                                                                                                                                                                                                                                                7260000.0*y0**2 + 27600000.0*y0*y4 + 4.58e-8*y1*y4 - 128000.0*y12*y3 - 39500.0*y3*y5 + 5140000.0*y4**2 + 0.0136],
[                                                                                                                                        0.0158*y0*y1 + 20900000.0*y0*y11 - 27600000.0*y0*y4 - 4.58e-8*y1*y4 + 128000.0*y12*y3 - 24400.0*y2*y4 + 39500.0*y3*y5 - 10280000.0*y4**2 - 10900000.0*y4*y5 - 36500.0*y4*y6 - 13100000.0*y4*y7 - 11400000.0*y4*y8 - 11400000.0*y4*y9 - 5.83*y4 + 0.0188],
[                                                                                              -35500000.0*y0*y5 + 13600000.0*y0*y6 + 1270.0*y1*y12 + 4.58e-8*y1*y4 - 4070000.0*y10*y5 + 50200000.0*y11*y12 - 13300000.0*y2*y5 - 39500.0*y3*y5 - 10900000.0*y4*y5 + 36500.0*y4*y6 + 22800000.0*y4*y9 - 9620000.0*y5**2 - 29200.0*y5*y6 - 11000000.0*y5*y8 - 8840000.0*y5*y9 - 0.0943*y5 + 0.0871],
[                                                                                                                                                                      -13600000.0*y0*y6 + 13000000.0*y0*y8 + 1270.0*y1*y10 + 50200000.0*y10*y11 - 4070000.0*y12*y6 - 13300000.0*y2*y6 - 36500.0*y4*y6 + 4810000.0*y5**2 - 29200.0*y5*y6 - 0.09430026*y6 + 101000.0*y8*y9 + 840.0*y9**2 + 0.0221],
[                                                                                                                                                                                                                                 -22900000.0*y0*y7 - 3750000.0*y12*y7 + 2620.0*y13 - 13100000.0*y4*y7 + 11000000.0*y5*y8 + 8840000.0*y5*y9 + 1.3e-7*y6 + 101000.0*y8*y9 + 840.0*y9**2 + 8.19e-7],
[                                                                                                                                                                                         22900000.0*y0*y7 - 13000000.0*y0*y8 - 0.000155*y1*y8 + 786000.0*y10*y12 + 4070000.0*y10*y5 - 50200000.0*y11*y8 + 13300000.0*y2*y9 - 11400000.0*y4*y8 - 11000000.0*y5*y8 - 101000.0*y8*y9 + 773000.0*y9],
[                                                                                                                                                                             -13000000.0*y0*y9 + 0.000155*y1*y8 + 50200000.0*y11*y8 + 4070000.0*y12*y6 - 13300000.0*y2*y9 + 13100000.0*y4*y7 - 11400000.0*y4*y9 + 29200.0*y5*y6 - 8840000.0*y5*y9 - 101000.0*y8*y9 - 1680.0*y9**2 - 773000.0*y9],
[                                                                                                                                                                                                                                                  13000000.0*y0*y9 - 1270.0*y1*y10 - 50200000.0*y10*y11 - 786000.0*y10*y12 - 4070000.0*y10*y5 + 13300000.0*y2*y6 + 11400000.0*y4*y8 + 0.0943*y6],
[                                                                                                                                                                                                                    -20900000.0*y0*y11 + 2.12e-5*y1 - 50200000.0*y10*y11 - 50200000.0*y11*y12 - 118000000.0*y11*y2 - 50200000.0*y11*y8 + 5.83*y4 + 0.0943*y5 + 0.0943*y6 + 773000.0*y9 + 0.0854],
[                                                                                                                                                                                                                                     -1270.0*y1*y12 - 786000.0*y10*y12 - 50200000.0*y11*y12 - 128000.0*y12*y3 - 4070000.0*y12*y6 - 3750000.0*y12*y7 + 2620.0*y13 + 13300000.0*y2*y5 + 0.0943*y5],
[                                                                                                                                                                                                                                                                                                                                                                  3750000.0*y12*y7 - 2620.0*y13]])

In this case, the C code printer does not do what we desire. It does not support printing a SymPy Matrix (see the first line of the output). In C, on possible representation of a matrix is an array type. The array type in C stores contigous values, e.g. doubles, in a chunk of memory. You can declare an array of doubles in C like:

double my_array[10];

The word double is the data type of the individual values in the array which must all be the same. The word my_array is the variable name we choose to name the array and the [10] is the syntax to declare that this array will have 10 values.

The array is "empty" when first declared and can be filled with values like so:

my_array[0] = 5;
my_array[1] = 6.78;
my array[2] = my_array[0] * 12;

or like:

my_array = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};

It is possible to declare multidimensional arrays in C that could map more directly to the indices of our two dimensional matrix, but in this case we will map our two dimensional matrix to a one dimenasional array using C contingous row ordering.

The code printers are capable of dealing with this need through the assign_to keyword argument in the .doprint() method but we must define a SymPy object that is appropriate to be assigned to. In our case, since we want to assign a Matrix we need to use an appropriately sized Matrix symbol.

In [10]:
rhs_result = sym.MatrixSymbol('rhs_result', 14, 1)
In [11]:
print(rhs_result)
rhs_result
In [12]:
print(rhs_result[0])
rhs_result[0, 0]
In [13]:
print(printer.doprint(rhs_of_odes, assign_to=rhs_result))
rhs_result[0] = -14520000.0*pow(y0, 2) - 0.0158*y0*y1 - 20900000.0*y0*y11 - 27600000.0*y0*y4 - 35500000.0*y0*y5 - 13600000.0*y0*y6 - 22900000.0*y0*y7 - 13000000.0*y0*y8 - 13000000.0*y0*y9 + 24400.0*y2*y4 + 5.83*y4 + 0.0854;
rhs_result[1] = -14520000.0*pow(y0, 2) - 0.0158*y0*y1 - 27600000.0*y0*y4 - 26000000.0*y0*y8 - 1270.0*y1*y10 - 1270.0*y1*y12 - 4.58e-8*y1*y4 - 0.000155*y1*y8 - 2.12e-5*y1 + 4070000.0*y10*y5 + 118000000.0*y11*y2 + 24400.0*y2*y4 + 13300000.0*y2*y5 + 13300000.0*y2*y6 + 13300000.0*y2*y9 + 39500.0*y3*y5 + 10900000.0*y4*y5 + 36500.0*y4*y6 + 29200.0*y5*y6 + 8840000.0*y5*y9 + 2.6e-7*y6 - 101000.0*y8*y9 - 0.131;
rhs_result[2] = 14520000.0*pow(y0, 2) + 0.0158*y0*y1 + 27600000.0*y0*y4 + 35500000.0*y0*y5 + 13600000.0*y0*y6 + 26000000.0*y0*y8 + 1270.0*y1*y10 + 1270.0*y1*y12 + 0.000155*y1*y8 + 2.12e-5*y1 + 786000.0*y10*y12 - 118000000.0*y11*y2 + 128000.0*y12*y3 + 4070000.0*y12*y6 - 24400.0*y2*y4 - 13300000.0*y2*y5 - 13300000.0*y2*y6 - 13300000.0*y2*y9 + 11000000.0*y5*y8 + 101000.0*y8*y9;
rhs_result[3] = 7260000.0*pow(y0, 2) + 27600000.0*y0*y4 + 4.58e-8*y1*y4 - 128000.0*y12*y3 - 39500.0*y3*y5 + 5140000.0*pow(y4, 2) + 0.0136;
rhs_result[4] = 0.0158*y0*y1 + 20900000.0*y0*y11 - 27600000.0*y0*y4 - 4.58e-8*y1*y4 + 128000.0*y12*y3 - 24400.0*y2*y4 + 39500.0*y3*y5 - 10280000.0*pow(y4, 2) - 10900000.0*y4*y5 - 36500.0*y4*y6 - 13100000.0*y4*y7 - 11400000.0*y4*y8 - 11400000.0*y4*y9 - 5.83*y4 + 0.0188;
rhs_result[5] = -35500000.0*y0*y5 + 13600000.0*y0*y6 + 1270.0*y1*y12 + 4.58e-8*y1*y4 - 4070000.0*y10*y5 + 50200000.0*y11*y12 - 13300000.0*y2*y5 - 39500.0*y3*y5 - 10900000.0*y4*y5 + 36500.0*y4*y6 + 22800000.0*y4*y9 - 9620000.0*pow(y5, 2) - 29200.0*y5*y6 - 11000000.0*y5*y8 - 8840000.0*y5*y9 - 0.0943*y5 + 0.0871;
rhs_result[6] = -13600000.0*y0*y6 + 13000000.0*y0*y8 + 1270.0*y1*y10 + 50200000.0*y10*y11 - 4070000.0*y12*y6 - 13300000.0*y2*y6 - 36500.0*y4*y6 + 4810000.0*pow(y5, 2) - 29200.0*y5*y6 - 0.09430026*y6 + 101000.0*y8*y9 + 840.0*pow(y9, 2) + 0.0221;
rhs_result[7] = -22900000.0*y0*y7 - 3750000.0*y12*y7 + 2620.0*y13 - 13100000.0*y4*y7 + 11000000.0*y5*y8 + 8840000.0*y5*y9 + 1.3e-7*y6 + 101000.0*y8*y9 + 840.0*pow(y9, 2) + 8.19e-7;
rhs_result[8] = 22900000.0*y0*y7 - 13000000.0*y0*y8 - 0.000155*y1*y8 + 786000.0*y10*y12 + 4070000.0*y10*y5 - 50200000.0*y11*y8 + 13300000.0*y2*y9 - 11400000.0*y4*y8 - 11000000.0*y5*y8 - 101000.0*y8*y9 + 773000.0*y9;
rhs_result[9] = -13000000.0*y0*y9 + 0.000155*y1*y8 + 50200000.0*y11*y8 + 4070000.0*y12*y6 - 13300000.0*y2*y9 + 13100000.0*y4*y7 - 11400000.0*y4*y9 + 29200.0*y5*y6 - 8840000.0*y5*y9 - 101000.0*y8*y9 - 1680.0*pow(y9, 2) - 773000.0*y9;
rhs_result[10] = 13000000.0*y0*y9 - 1270.0*y1*y10 - 50200000.0*y10*y11 - 786000.0*y10*y12 - 4070000.0*y10*y5 + 13300000.0*y2*y6 + 11400000.0*y4*y8 + 0.0943*y6;
rhs_result[11] = -20900000.0*y0*y11 + 2.12e-5*y1 - 50200000.0*y10*y11 - 50200000.0*y11*y12 - 118000000.0*y11*y2 - 50200000.0*y11*y8 + 5.83*y4 + 0.0943*y5 + 0.0943*y6 + 773000.0*y9 + 0.0854;
rhs_result[12] = -1270.0*y1*y12 - 786000.0*y10*y12 - 50200000.0*y11*y12 - 128000.0*y12*y3 - 4070000.0*y12*y6 - 3750000.0*y12*y7 + 2620.0*y13 + 13300000.0*y2*y5 + 0.0943*y5;
rhs_result[13] = 3750000.0*y12*y7 - 2620.0*y13;

Notice that we have proper array value assignment and valid lines of C code that can be used in our function.

Excercise [5 min]

Print out valid C code for the Jacobian matrix.

Double Click For Solution

In [14]:
# write your answer here

Changing the Behavior of the Printer

The SymPy code printers are relatively easy to extend. They are designed such that if you want to change how a particularly SymPy object prints, for example a Symbol, then you only need to modify the _print_Symbol method of the printer. In general, the code printers have a method for every SymPy object and also many builtin types. Use tab completion with C99CodePrinter._print_ to see all of the options.

Once you find the method you want to modify, it is often useful to look at the existing impelementation of the print method to see how the code is written.

In [15]:
C99CodePrinter._print_Symbol??

Below is a simple example of overiding the Symbol printer method. Note that you should use the self._print() method instead of simply returning the string so that the proper printer, self._print_str(), is dispatched. This is most important if you are printing non-singletons, i.e. expressions that are made up of multiple singletons.

In [16]:
C99CodePrinter._print_str??
In [17]:
class MyCodePrinter(C99CodePrinter):
    def _print_Symbol(self, expr):
        return self._print("No matter what symbol you pass in I will always print:\n\nNi!")
In [18]:
my_printer = MyCodePrinter()
In [19]:
theta = sym.symbols('theta')
theta
Out[19]:
$$\theta$$
In [20]:
print(my_printer.doprint(theta))
No matter what symbol you pass in I will always print:

Ni!

Exercise [10 min]

One issue with our current code printer is that the expressions use the symbols y0, y1, ..., y13 instead of accessing the values directly from the arrays with state_vals[0], state_vals[1], ..., state_vals[13]. We could go back and rename our SymPy symbols to use brackets, but another way would be to override the _print_Symbol() method to print these symbols as we desire. Modify the code printer so that it prints with the proper array access in the expression.

Double Click For Solution: Subclassing

Double Click For Solution: Exact replacement

In [21]:
# write your answer here

Bonus Exercise

Do this exercise if you finish the previous one quickly.

It turns out that calling pow() for low value integer exponents executes slower than simply expanding the multiplication. For example pow(x, 2) could be printed as x*x. Modify the CCodePrinter ._print_Pow method to expand the multiplication if the exponent is less than or equal to 4. You may want to have a look at the source code with printer._print_Pow??

Note that a Pow expression has an .exp for exponent and .base for the item being raised. For example $x^2$ would have:

expr = x**2
expr.base == x
expr.exp == 2

Double Click for Solution

In [22]:
# write your answer here

Common Subexpression Elimination

If you look carefully at the expressions in the two matrices you'll see repeated expressions. These are not ideal in the sense that the computer has to repeat the exact same calculation multiple times. For large expressions this can be a major issue. Compilers, such as gcc, can often eliminate common subexpressions on their own when different optimization flags are invoked but for complex expressions the algorithms in some compilers do not do a thorough job or compilation can take an extremely long time. SymPy has tools to perform common subexpression elimination which is both thorough and reasonably efficient. In particular if gcc is run with the lowest optimization setting -O0 cse can give large speedups.

For example if you have two expressions:

a = x*y + 5
b = x*y + 6

you can convert this to these three expressions:

z = x*y
a = z + 5
b = z + 6

and x*y only has to be computed once.

The cse() function in SymPy returns the subexpression, z = x*y, and the simplified expressions: a = z + 5, b = z + 6.

Here is how it works:

In [23]:
sm.cse?
Object `sm.cse` not found.
In [24]:
sub_exprs, simplified_rhs = sym.cse(rhs_of_odes)
In [25]:
for var, expr in sub_exprs:
    sym.pprint(sym.Eq(var, expr))
x₀ = 0.0158⋅y₀⋅y₁
x₁ = -x₀
x₂ = 27600000.0⋅y₀⋅y₄
x₃ = -x₂
x₄ = 24400.0⋅y₂⋅y₄
       2
x₅ = y₀ 
x₆ = 14520000.0⋅x₅
x₇ = -x₆
x₈ = 5.83⋅y₄
x₉ = 20900000.0⋅y₀⋅y₁₁
x₁₀ = -x₉
x₁₁ = 35500000.0⋅y₀⋅y₅
x₁₂ = -x₁₁
x₁₃ = 13600000.0⋅y₀⋅y₆
x₁₄ = -x₁₃
x₁₅ = 22900000.0⋅y₀⋅y₇
x₁₆ = -x₁₅
x₁₇ = y₀⋅y₈
x₁₈ = 13000000.0⋅x₁₇
x₁₉ = -x₁₈
x₂₀ = 13000000.0⋅y₀⋅y₉
x₂₁ = -x₂₀
x₂₂ = 4.58e-8⋅y₁⋅y₄
x₂₃ = 39500.0⋅y₃⋅y₅
x₂₄ = -x₂₂ + x₂₃
x₂₅ = 2.12e-5⋅y₁
x₂₆ = 26000000.0⋅x₁₇
x₂₇ = 1270.0⋅y₁
x₂₈ = x₂₇⋅y₁₀
x₂₉ = -x₂₈
x₃₀ = x₂₇⋅y₁₂
x₃₁ = -x₃₀
x₃₂ = 0.000155⋅y₁⋅y₈
x₃₃ = -x₃₂
x₃₄ = 4070000.0⋅y₁₀⋅y₅
x₃₅ = 118000000.0⋅y₁₁⋅y₂
x₃₆ = 13300000.0⋅y₂
x₃₇ = x₃₆⋅y₅
x₃₈ = x₃₆⋅y₆
x₃₉ = x₃₆⋅y₉
x₄₀ = 10900000.0⋅y₄⋅y₅
x₄₁ = 36500.0⋅y₄⋅y₆
x₄₂ = 29200.0⋅y₅⋅y₆
x₄₃ = 8840000.0⋅y₅⋅y₉
x₄₄ = 101000.0⋅y₈⋅y₉
x₄₅ = -x₄₄
x₄₆ = x₁₃ + x₃₀ - x₃₇
x₄₇ = 786000.0⋅y₁₀⋅y₁₂
x₄₈ = -x₃₅
x₄₉ = 128000.0⋅y₁₂⋅y₃
x₅₀ = 4070000.0⋅y₁₂⋅y₆
x₅₁ = -x₄
x₅₂ = -x₃₈
x₅₃ = -x₃₉
x₅₄ = 11000000.0⋅y₅⋅y₈
x₅₅ = -x₂₃
x₅₆ = -x₄₉
        2
x₅₇ = y₄ 
x₅₈ = -x₄₀
x₅₉ = -x₄₁
x₆₀ = 13100000.0⋅y₄⋅y₇
x₆₁ = -x₆₀
x₆₂ = 11400000.0⋅y₄
x₆₃ = x₆₂⋅y₈
x₆₄ = -x₆₃
x₆₅ = -x₆₂⋅y₉
x₆₆ = 0.0943⋅y₅
x₆₇ = -x₃₄
x₆₈ = 50200000.0⋅y₁₁
x₆₉ = x₆₈⋅y₁₂
x₇₀ = -x₄₂
x₇₁ = -x₅₄
x₇₂ = -x₄₃
        2
x₇₃ = y₅ 
x₇₄ = x₆₈⋅y₁₀
x₇₅ = -x₅₀
        2
x₇₆ = y₉ 
x₇₇ = 840.0⋅x₇₆
x₇₈ = 2620.0⋅y₁₃
x₇₉ = 3750000.0⋅y₁₂⋅y₇
x₈₀ = -x₇₉
x₈₁ = 773000.0⋅y₉
x₈₂ = x₆₈⋅y₈
x₈₃ = -x₈₂
x₈₄ = 0.0943⋅y₆
x₈₅ = -x₇₄
x₈₆ = -x₄₇
x₈₇ = -x₆₉

cse() can return a number of simplified expressions and to do this it returns a list. In our case we have 1 simplified expression that can be accessed as the first item of the list.

In [26]:
type(simplified_rhs)
Out[26]:
list
In [27]:
len(simplified_rhs)
Out[27]:
$$1$$
In [28]:
simplified_rhs[0]
Out[28]:
$$\left[\begin{matrix}x_{1} + x_{10} + x_{12} + x_{14} + x_{16} + x_{19} + x_{21} + x_{3} + x_{4} + x_{7} + x_{8} + 0.0854\\x_{1} + x_{24} - x_{25} - x_{26} + x_{29} + x_{3} + x_{31} + x_{33} + x_{34} + x_{35} + x_{37} + x_{38} + x_{39} + x_{4} + x_{40} + x_{41} + x_{42} + x_{43} + x_{45} + x_{7} + 2.6 \cdot 10^{-7} y_{6} - 0.131\\x_{0} + x_{11} + x_{2} + x_{25} + x_{26} + x_{28} + x_{32} + x_{44} + x_{46} + x_{47} + x_{48} + x_{49} + x_{50} + x_{51} + x_{52} + x_{53} + x_{54} + x_{6}\\x_{2} + x_{22} + 7260000.0 x_{5} + x_{55} + x_{56} + 5140000.0 x_{57} + 0.0136\\x_{0} + x_{24} + x_{3} + x_{49} + x_{51} - 10280000.0 x_{57} + x_{58} + x_{59} + x_{61} + x_{64} + x_{65} - x_{8} + x_{9} + 0.0188\\x_{12} + x_{22} + x_{41} + x_{46} + x_{55} + x_{58} - x_{66} + x_{67} + x_{69} + x_{70} + x_{71} + x_{72} - 9620000.0 x_{73} + 22800000.0 y_{4} y_{9} + 0.0871\\x_{14} + x_{18} + x_{28} + x_{44} + x_{52} + x_{59} + x_{70} + 4810000.0 x_{73} + x_{74} + x_{75} + x_{77} - 0.09430026 y_{6} + 0.0221\\x_{16} + x_{43} + x_{44} + x_{54} + x_{61} + x_{77} + x_{78} + x_{80} + 1.3 \cdot 10^{-7} y_{6} + 8.19 \cdot 10^{-7}\\x_{15} + x_{19} + x_{33} + x_{34} + x_{39} + x_{45} + x_{47} + x_{64} + x_{71} + x_{81} + x_{83}\\x_{21} + x_{32} + x_{42} + x_{45} + x_{50} + x_{53} + x_{60} + x_{65} + x_{72} - 1680.0 x_{76} - x_{81} + x_{82}\\x_{20} + x_{29} + x_{38} + x_{63} + x_{67} + x_{84} + x_{85} + x_{86}\\x_{10} + x_{25} + x_{48} + x_{66} + x_{8} + x_{81} + x_{83} + x_{84} + x_{85} + x_{87} + 0.0854\\x_{31} + x_{37} + x_{56} + x_{66} + x_{75} + x_{78} + x_{80} + x_{86} + x_{87}\\- x_{78} + x_{79}\end{matrix}\right]$$

You can find common subexpressions among multiple objects also:

In [29]:
jac_of_odes = rhs_of_odes.jacobian(states)

sub_exprs, simplified_exprs = sym.cse((rhs_of_odes, jac_of_odes))
In [30]:
for var, expr in sub_exprs:
    sym.pprint(sym.Eq(var, expr))
x₀ = 0.0158⋅y₀
x₁ = x₀⋅y₁
x₂ = -x₁
x₃ = 27600000.0⋅y₀
x₄ = x₃⋅y₄
x₅ = -x₄
x₆ = 24400.0⋅y₂
x₇ = x₆⋅y₄
       2
x₈ = y₀ 
x₉ = 14520000.0⋅x₈
x₁₀ = -x₉
x₁₁ = 5.83⋅y₄
x₁₂ = 20900000.0⋅y₀
x₁₃ = x₁₂⋅y₁₁
x₁₄ = -x₁₃
x₁₅ = 35500000.0⋅y₀
x₁₆ = x₁₅⋅y₅
x₁₇ = -x₁₆
x₁₈ = 13600000.0⋅y₀
x₁₉ = x₁₈⋅y₆
x₂₀ = -x₁₉
x₂₁ = 22900000.0⋅y₀
x₂₂ = x₂₁⋅y₇
x₂₃ = -x₂₂
x₂₄ = 13000000.0⋅y₀
x₂₅ = x₂₄⋅y₈
x₂₆ = -x₂₅
x₂₇ = x₂₄⋅y₉
x₂₈ = -x₂₇
x₂₉ = 4.58e-8⋅y₁
x₃₀ = x₂₉⋅y₄
x₃₁ = 39500.0⋅y₃
x₃₂ = x₃₁⋅y₅
x₃₃ = -x₃₀ + x₃₂
x₃₄ = 2.12e-5⋅y₁
x₃₅ = 26000000.0⋅y₀
x₃₆ = x₃₅⋅y₈
x₃₇ = 1270.0⋅y₁
x₃₈ = x₃₇⋅y₁₀
x₃₉ = -x₃₈
x₄₀ = x₃₇⋅y₁₂
x₄₁ = -x₄₀
x₄₂ = 0.000155⋅y₁
x₄₃ = x₄₂⋅y₈
x₄₄ = -x₄₃
x₄₅ = 4070000.0⋅y₁₀
x₄₆ = x₄₅⋅y₅
x₄₇ = 118000000.0⋅y₁₁
x₄₈ = x₄₇⋅y₂
x₄₉ = 13300000.0⋅y₂
x₅₀ = x₄₉⋅y₅
x₅₁ = x₄₉⋅y₆
x₅₂ = x₄₉⋅y₉
x₅₃ = 10900000.0⋅y₄
x₅₄ = x₅₃⋅y₅
x₅₅ = 36500.0⋅y₄
x₅₆ = x₅₅⋅y₆
x₅₇ = 29200.0⋅y₅
x₅₈ = x₅₇⋅y₆
x₅₉ = 8840000.0⋅y₅
x₆₀ = x₅₉⋅y₉
x₆₁ = 101000.0⋅y₈
x₆₂ = x₆₁⋅y₉
x₆₃ = -x₆₂
x₆₄ = x₁₉ + x₄₀ - x₅₀
x₆₅ = 786000.0⋅y₁₀
x₆₆ = x₆₅⋅y₁₂
x₆₇ = -x₄₈
x₆₈ = 128000.0⋅y₁₂
x₆₉ = x₆₈⋅y₃
x₇₀ = 4070000.0⋅y₁₂
x₇₁ = x₇₀⋅y₆
x₇₂ = -x₇
x₇₃ = -x₅₁
x₇₄ = -x₅₂
x₇₅ = 11000000.0⋅y₅
x₇₆ = x₇₅⋅y₈
x₇₇ = -x₃₂
x₇₈ = -x₆₉
        2
x₇₉ = y₄ 
x₈₀ = -x₅₄
x₈₁ = -x₅₆
x₈₂ = 13100000.0⋅y₄
x₈₃ = x₈₂⋅y₇
x₈₄ = -x₈₃
x₈₅ = 11400000.0⋅y₄
x₈₆ = x₈₅⋅y₈
x₈₇ = -x₈₆
x₈₈ = -x₈₅⋅y₉
x₈₉ = 0.0943⋅y₅
x₉₀ = -x₄₆
x₉₁ = 50200000.0⋅y₁₁
x₉₂ = x₉₁⋅y₁₂
x₉₃ = 22800000.0⋅y₄
x₉₄ = -x₅₈
x₉₅ = -x₇₆
x₉₆ = -x₆₀
        2
x₉₇ = y₅ 
x₉₈ = 50200000.0⋅y₁₀
x₉₉ = x₉₈⋅y₁₁
x₁₀₀ = -x₇₁
         2
x₁₀₁ = y₉ 
x₁₀₂ = 840.0⋅x₁₀₁
x₁₀₃ = 2620.0⋅y₁₃
x₁₀₄ = 3750000.0⋅y₁₂
x₁₀₅ = x₁₀₄⋅y₇
x₁₀₆ = -x₁₀₅
x₁₀₇ = 773000.0⋅y₉
x₁₀₈ = x₉₁⋅y₈
x₁₀₉ = -x₁₀₈
x₁₁₀ = 0.0943⋅y₆
x₁₁₁ = -x₉₉
x₁₁₂ = -x₆₆
x₁₁₃ = -x₉₂
x₁₁₄ = 29040000.0⋅y₀
x₁₁₅ = 0.0158⋅y₁
x₁₁₆ = 27600000.0⋅y₄
x₁₁₇ = -x₁₁₆
x₁₁₈ = -x₁₁₄ - x₁₁₅ + x₁₁₇
x₁₁₉ = 20900000.0⋅y₁₁
x₁₂₀ = -x₁₁₉
x₁₂₁ = 35500000.0⋅y₅
x₁₂₂ = -x₁₂₁
x₁₂₃ = 13600000.0⋅y₆
x₁₂₄ = -x₁₂₃
x₁₂₅ = 22900000.0⋅y₇
x₁₂₆ = -x₁₂₅
x₁₂₇ = 13000000.0⋅y₈
x₁₂₈ = -x₁₂₇
x₁₂₉ = 13000000.0⋅y₉
x₁₃₀ = -x₁₂₉
x₁₃₁ = -x₀
x₁₃₂ = 24400.0⋅y₄
x₁₃₃ = -x₃
x₁₃₄ = x₁₃₃ + x₆
x₁₃₅ = -x₁₅
x₁₃₆ = -x₁₈
x₁₃₇ = -x₂₁
x₁₃₈ = -x₂₄
x₁₃₉ = -x₁₂
x₁₄₀ = 26000000.0⋅y₈
x₁₄₁ = 1270.0⋅y₁₀
x₁₄₂ = -x₁₄₁
x₁₄₃ = 1270.0⋅y₁₂
x₁₄₄ = -x₁₄₃
x₁₄₅ = 4.58e-8⋅y₄
x₁₄₆ = -x₁₄₅
x₁₄₇ = 0.000155⋅y₈
x₁₄₈ = -x₁₄₇
x₁₄₉ = 13300000.0⋅y₅
x₁₅₀ = 13300000.0⋅y₆
x₁₅₁ = 13300000.0⋅y₉
x₁₅₂ = 39500.0⋅y₅
x₁₅₃ = -x₂₉
x₁₅₄ = 10900000.0⋅y₅
x₁₅₅ = 36500.0⋅y₆
x₁₅₆ = 29200.0⋅y₆
x₁₅₇ = 8840000.0⋅y₉
x₁₅₈ = -x₄₂
x₁₅₉ = 101000.0⋅y₉
x₁₆₀ = -x₁₅₉
x₁₆₁ = -x₆₁
x₁₆₂ = x₁₆₁ + x₄₉
x₁₆₃ = -x₃₇
x₁₆₄ = 4070000.0⋅y₅
x₁₆₅ = 118000000.0⋅y₂
x₁₆₆ = -x₄₇
x₁₆₇ = -x₁₃₂
x₁₆₈ = -x₁₄₉
x₁₆₉ = -x₁₅₀
x₁₇₀ = -x₁₅₁
x₁₇₁ = -x₆
x₁₇₂ = -x₄₉
x₁₇₃ = 11000000.0⋅y₈
x₁₇₄ = x₁₅₉ + x₇₅
x₁₇₅ = 786000.0⋅y₁₂
x₁₇₆ = -x₁₆₅
x₁₇₇ = 128000.0⋅y₃
x₁₇₈ = 4070000.0⋅y₆
x₁₇₉ = -x₆₈
x₁₈₀ = -x₁₅₂
x₁₈₁ = -x₃₁
x₁₈₂ = -x₁₇₇
x₁₈₃ = -x₁₅₄
x₁₈₄ = -x₁₅₅
x₁₈₅ = 13100000.0⋅y₇
x₁₈₆ = -x₁₈₅
x₁₈₇ = 11400000.0⋅y₈
x₁₈₈ = -x₁₈₇
x₁₈₉ = -11400000.0⋅y₉
x₁₉₀ = -x₅₃
x₁₉₁ = -x₅₅
x₁₉₂ = -x₈₂
x₁₉₃ = -x₈₅
x₁₉₄ = -x₄₅
x₁₉₅ = -x₁₅₆
x₁₉₆ = -x₁₇₃
x₁₉₇ = -x₁₅₇
x₁₉₈ = -x₅₇
x₁₉₉ = -x₇₅
x₂₀₀ = -x₅₉
x₂₀₁ = -x₁₆₄
x₂₀₂ = 50200000.0⋅y₁₂
x₂₀₃ = x₃₇ + x₉₁
x₂₀₄ = -x₇₀
x₂₀₅ = x₆₁ + 1680.0⋅y₉
x₂₀₆ = -x₁₇₈
x₂₀₇ = -x₁₀₄
x₂₀₈ = 3750000.0⋅y₇
x₂₀₉ = -x₂₀₈
x₂₁₀ = x₁₃₈ + x₁₉₃
x₂₁₁ = -x₉₁
x₂₁₂ = 50200000.0⋅y₈
x₂₁₃ = -x₂₁₂
x₂₁₄ = x₄₉ + 0.0943
x₂₁₅ = x₁₆₃ + x₂₁₁
x₂₁₆ = -x₁₇₅
x₂₁₇ = -x₉₈
x₂₁₈ = -x₆₅
x₂₁₉ = -x₂₀₂
In [31]:
simplified_exprs[0]
Out[31]:
$$\left[\begin{matrix}x_{10} + x_{11} + x_{14} + x_{17} + x_{2} + x_{20} + x_{23} + x_{26} + x_{28} + x_{5} + x_{7} + 0.0854\\x_{10} + x_{2} + x_{33} - x_{34} - x_{36} + x_{39} + x_{41} + x_{44} + x_{46} + x_{48} + x_{5} + x_{50} + x_{51} + x_{52} + x_{54} + x_{56} + x_{58} + x_{60} + x_{63} + x_{7} + 2.6 \cdot 10^{-7} y_{6} - 0.131\\x_{1} + x_{16} + x_{34} + x_{36} + x_{38} + x_{4} + x_{43} + x_{62} + x_{64} + x_{66} + x_{67} + x_{69} + x_{71} + x_{72} + x_{73} + x_{74} + x_{76} + x_{9}\\x_{30} + x_{4} + x_{77} + x_{78} + 5140000.0 x_{79} + 7260000.0 x_{8} + 0.0136\\x_{1} - x_{11} + x_{13} + x_{33} + x_{5} + x_{69} + x_{72} - 10280000.0 x_{79} + x_{80} + x_{81} + x_{84} + x_{87} + x_{88} + 0.0188\\x_{17} + x_{30} + x_{56} + x_{64} + x_{77} + x_{80} - x_{89} + x_{90} + x_{92} + x_{93} y_{9} + x_{94} + x_{95} + x_{96} - 9620000.0 x_{97} + 0.0871\\x_{100} + x_{102} + x_{20} + x_{25} + x_{38} + x_{62} + x_{73} + x_{81} + x_{94} + 4810000.0 x_{97} + x_{99} - 0.09430026 y_{6} + 0.0221\\x_{102} + x_{103} + x_{106} + x_{23} + x_{60} + x_{62} + x_{76} + x_{84} + 1.3 \cdot 10^{-7} y_{6} + 8.19 \cdot 10^{-7}\\x_{107} + x_{109} + x_{22} + x_{26} + x_{44} + x_{46} + x_{52} + x_{63} + x_{66} + x_{87} + x_{95}\\- 1680.0 x_{101} - x_{107} + x_{108} + x_{28} + x_{43} + x_{58} + x_{63} + x_{71} + x_{74} + x_{83} + x_{88} + x_{96}\\x_{110} + x_{111} + x_{112} + x_{27} + x_{39} + x_{51} + x_{86} + x_{90}\\x_{107} + x_{109} + x_{11} + x_{110} + x_{111} + x_{113} + x_{14} + x_{34} + x_{67} + x_{89} + 0.0854\\x_{100} + x_{103} + x_{106} + x_{112} + x_{113} + x_{41} + x_{50} + x_{78} + x_{89}\\- x_{103} + x_{105}\end{matrix}\right]$$
In [32]:
simplified_exprs[1]
Out[32]:
$$\left[\begin{array}{cccccccccccccc}x_{118} + x_{120} + x_{122} + x_{124} + x_{126} + x_{128} + x_{130} & x_{131} & x_{132} & 0 & x_{134} + 5.83 & x_{135} & x_{136} & x_{137} & x_{138} & x_{138} & 0 & x_{139} & 0 & 0\\x_{118} - x_{140} & x_{131} + x_{142} + x_{144} + x_{146} + x_{148} - 2.12 \cdot 10^{-5} & x_{132} + x_{149} + x_{150} + x_{151} + x_{47} & x_{152} & x_{134} + x_{153} + x_{154} + x_{155} & x_{156} + x_{157} + x_{31} + x_{45} + x_{49} + x_{53} & x_{49} + x_{55} + x_{57} + 2.6 \cdot 10^{-7} & 0 & x_{158} + x_{160} - x_{35} & x_{162} + x_{59} & x_{163} + x_{164} & x_{165} & x_{163} & 0\\x_{114} + x_{115} + x_{116} + x_{121} + x_{123} + x_{140} & x_{0} + x_{141} + x_{143} + x_{147} + 2.12 \cdot 10^{-5} & x_{166} + x_{167} + x_{168} + x_{169} + x_{170} & x_{68} & x_{171} + x_{3} & x_{15} + x_{172} + x_{173} & x_{172} + x_{18} + x_{70} & 0 & x_{174} + x_{35} + x_{42} & x_{172} + x_{61} & x_{175} + x_{37} & x_{176} & x_{177} + x_{178} + x_{37} + x_{65} & 0\\x_{116} + 14520000.0 y_{0} & x_{145} & 0 & x_{179} + x_{180} & x_{29} + x_{3} + 10280000.0 y_{4} & x_{181} & 0 & 0 & 0 & 0 & 0 & 0 & x_{182} & 0\\x_{115} + x_{117} + x_{119} & x_{0} + x_{146} & x_{167} & x_{152} + x_{68} & x_{133} + x_{153} + x_{171} + x_{183} + x_{184} + x_{186} + x_{188} + x_{189} - 20560000.0 y_{4} - 5.83 & x_{190} + x_{31} & x_{191} & x_{192} & x_{193} & x_{193} & 0 & x_{12} & x_{177} & 0\\x_{122} + x_{123} & x_{143} + x_{145} & x_{168} & x_{180} & x_{155} + x_{183} + x_{29} + 22800000.0 y_{9} & x_{135} + x_{172} + x_{181} + x_{190} + x_{194} + x_{195} + x_{196} + x_{197} - 19240000.0 y_{5} - 0.0943 & x_{18} + x_{198} + x_{55} & 0 & x_{199} & x_{200} + x_{93} & x_{201} & x_{202} & x_{203} & 0\\x_{124} + x_{127} & x_{141} & x_{169} & 0 & x_{184} & x_{195} + 9620000.0 y_{5} & x_{136} + x_{172} + x_{191} + x_{198} + x_{204} - 0.09430026 & 0 & x_{159} + x_{24} & x_{205} & x_{203} & x_{98} & x_{206} & 0\\x_{126} & 0 & 0 & 0 & x_{186} & x_{157} + x_{173} & 1.3 \cdot 10^{-7} & x_{137} + x_{192} + x_{207} & x_{174} & x_{205} + x_{59} & 0 & 0 & x_{209} & 2620.0\\x_{125} + x_{128} & x_{148} & x_{151} & 0 & x_{188} & x_{196} + x_{45} & 0 & x_{21} & x_{158} + x_{160} + x_{199} + x_{210} + x_{211} & x_{162} + 773000.0 & x_{164} + x_{175} & x_{213} & x_{65} & 0\\x_{130} & x_{147} & x_{170} & 0 & x_{185} + x_{189} & x_{156} + x_{197} & x_{57} + x_{70} & x_{82} & x_{160} + x_{42} + x_{91} & x_{161} + x_{172} + x_{200} + x_{210} - 3360.0 y_{9} - 773000.0 & 0 & x_{212} & x_{178} & 0\\x_{129} & x_{142} & x_{150} & 0 & x_{187} & x_{194} & x_{214} & 0 & x_{85} & x_{24} & x_{201} + x_{215} + x_{216} & x_{217} & x_{218} & 0\\x_{120} & 2.12 \cdot 10^{-5} & x_{166} & 0 & 5.83 & 0.0943 & 0.0943 & 0 & x_{211} & 773000.0 & x_{211} & x_{139} + x_{176} + x_{213} + x_{217} + x_{219} & x_{211} & 0\\0 & x_{144} & x_{149} & x_{179} & 0 & x_{214} & x_{204} & x_{207} & 0 & 0 & x_{216} & x_{219} & x_{182} + x_{206} + x_{209} + x_{215} + x_{218} & 2620.0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & x_{104} & 0 & 0 & 0 & 0 & x_{208} & -2620.0\end{array}\right]$$

Exercise [15min]

Use common subexpression elimination to print out C code for your two arrays such that:

double x0 = first_sub_expression;
...
double xN = last_sub_expression;

rhs_result[0] = expressions_containing_the_subexpressions;
...
rhs_result[13] = ...;

jac_result[0] = ...;
...
jac_result[195] = ...;

The code you create can be copied and pasted into the provided template above to make a C program. Refer back to the introduction to C code printing above.

To give you a bit of help we will first introduce the Assignment class. The printers know how to print variable assignments that are defined by an Assignment instance.

In [33]:
from sympy.printing.codeprinter import Assignment

print(printer.doprint(Assignment(theta, 5)))
theta = 5;

The following code demonstrates a way to use cse() to simplify single matrix objects. Note that we use ImmutableDenseMatrix because all dense matrics are internally converted to this type in the printers. Check the type of your matrices to see.

In [34]:
class CMatrixPrinter(C99CodePrinter):
    def _print_ImmutableDenseMatrix(self, expr):
        sub_exprs, simplified = sym.cse(expr)
        lines = []
        for var, sub_expr in sub_exprs:
            lines.append('double ' + self._print(Assignment(var, sub_expr)))
        M = sym.MatrixSymbol('M', *expr.shape)
        return '\n'.join(lines) + '\n' + self._print(Assignment(M, expr))
In [35]:
p = CMatrixPrinter()
print(p.doprint(jac_of_odes))
double x0 = 29040000.0*y0;
double x1 = 0.0158*y1;
double x2 = 27600000.0*y4;
double x3 = -x2;
double x4 = -x0 - x1 + x3;
double x5 = 20900000.0*y11;
double x6 = -x5;
double x7 = 35500000.0*y5;
double x8 = -x7;
double x9 = 13600000.0*y6;
double x10 = -x9;
double x11 = 22900000.0*y7;
double x12 = -x11;
double x13 = 13000000.0*y8;
double x14 = -x13;
double x15 = 13000000.0*y9;
double x16 = -x15;
double x17 = 0.0158*y0;
double x18 = -x17;
double x19 = 24400.0*y4;
double x20 = 27600000.0*y0;
double x21 = -x20;
double x22 = 24400.0*y2;
double x23 = x21 + x22;
double x24 = 35500000.0*y0;
double x25 = -x24;
double x26 = 13600000.0*y0;
double x27 = -x26;
double x28 = 22900000.0*y0;
double x29 = -x28;
double x30 = 13000000.0*y0;
double x31 = -x30;
double x32 = 20900000.0*y0;
double x33 = -x32;
double x34 = 26000000.0*y8;
double x35 = 1270.0*y10;
double x36 = -x35;
double x37 = 1270.0*y12;
double x38 = -x37;
double x39 = 4.58e-8*y4;
double x40 = -x39;
double x41 = 0.000155*y8;
double x42 = -x41;
double x43 = 118000000.0*y11;
double x44 = 13300000.0*y5;
double x45 = 13300000.0*y6;
double x46 = 13300000.0*y9;
double x47 = 39500.0*y5;
double x48 = 4.58e-8*y1;
double x49 = -x48;
double x50 = 10900000.0*y5;
double x51 = 36500.0*y6;
double x52 = 4070000.0*y10;
double x53 = 13300000.0*y2;
double x54 = 39500.0*y3;
double x55 = 10900000.0*y4;
double x56 = 29200.0*y6;
double x57 = 8840000.0*y9;
double x58 = 36500.0*y4;
double x59 = 29200.0*y5;
double x60 = 26000000.0*y0;
double x61 = 0.000155*y1;
double x62 = -x61;
double x63 = 101000.0*y9;
double x64 = -x63;
double x65 = 8840000.0*y5;
double x66 = 101000.0*y8;
double x67 = -x66;
double x68 = x53 + x67;
double x69 = 1270.0*y1;
double x70 = -x69;
double x71 = 4070000.0*y5;
double x72 = 118000000.0*y2;
double x73 = -x43;
double x74 = -x19;
double x75 = -x44;
double x76 = -x45;
double x77 = -x46;
double x78 = 128000.0*y12;
double x79 = -x22;
double x80 = -x53;
double x81 = 11000000.0*y8;
double x82 = 4070000.0*y12;
double x83 = 11000000.0*y5;
double x84 = x63 + x83;
double x85 = 786000.0*y12;
double x86 = -x72;
double x87 = 786000.0*y10;
double x88 = 128000.0*y3;
double x89 = 4070000.0*y6;
double x90 = -x78;
double x91 = -x47;
double x92 = -x54;
double x93 = -x88;
double x94 = -x50;
double x95 = -x51;
double x96 = 13100000.0*y7;
double x97 = -x96;
double x98 = 11400000.0*y8;
double x99 = -x98;
double x100 = -11400000.0*y9;
double x101 = -x55;
double x102 = -x58;
double x103 = 13100000.0*y4;
double x104 = -x103;
double x105 = 11400000.0*y4;
double x106 = -x105;
double x107 = -x52;
double x108 = -x56;
double x109 = -x81;
double x110 = -x57;
double x111 = -x59;
double x112 = -x83;
double x113 = -x65;
double x114 = -x71;
double x115 = 50200000.0*y12;
double x116 = 50200000.0*y11;
double x117 = x116 + x69;
double x118 = -x82;
double x119 = x66 + 1680.0*y9;
double x120 = 50200000.0*y10;
double x121 = -x89;
double x122 = 3750000.0*y12;
double x123 = -x122;
double x124 = 3750000.0*y7;
double x125 = -x124;
double x126 = x106 + x31;
double x127 = -x116;
double x128 = 50200000.0*y8;
double x129 = -x128;
double x130 = x53 + 0.0943;
double x131 = x127 + x70;
double x132 = -x85;
double x133 = -x120;
double x134 = -x87;
double x135 = -x115;
M[0] = -29040000.0*y0 - 0.0158*y1 - 20900000.0*y11 - 27600000.0*y4 - 35500000.0*y5 - 13600000.0*y6 - 22900000.0*y7 - 13000000.0*y8 - 13000000.0*y9;
M[1] = -0.0158*y0;
M[2] = 24400.0*y4;
M[3] = 0;
M[4] = -27600000.0*y0 + 24400.0*y2 + 5.83;
M[5] = -35500000.0*y0;
M[6] = -13600000.0*y0;
M[7] = -22900000.0*y0;
M[8] = -13000000.0*y0;
M[9] = -13000000.0*y0;
M[10] = 0;
M[11] = -20900000.0*y0;
M[12] = 0;
M[13] = 0;
M[14] = -29040000.0*y0 - 0.0158*y1 - 27600000.0*y4 - 26000000.0*y8;
M[15] = -0.0158*y0 - 1270.0*y10 - 1270.0*y12 - 4.58e-8*y4 - 0.000155*y8 - 2.12e-5;
M[16] = 118000000.0*y11 + 24400.0*y4 + 13300000.0*y5 + 13300000.0*y6 + 13300000.0*y9;
M[17] = 39500.0*y5;
M[18] = -27600000.0*y0 - 4.58e-8*y1 + 24400.0*y2 + 10900000.0*y5 + 36500.0*y6;
M[19] = 4070000.0*y10 + 13300000.0*y2 + 39500.0*y3 + 10900000.0*y4 + 29200.0*y6 + 8840000.0*y9;
M[20] = 13300000.0*y2 + 36500.0*y4 + 29200.0*y5 + 2.6e-7;
M[21] = 0;
M[22] = -26000000.0*y0 - 0.000155*y1 - 101000.0*y9;
M[23] = 13300000.0*y2 + 8840000.0*y5 - 101000.0*y8;
M[24] = -1270.0*y1 + 4070000.0*y5;
M[25] = 118000000.0*y2;
M[26] = -1270.0*y1;
M[27] = 0;
M[28] = 29040000.0*y0 + 0.0158*y1 + 27600000.0*y4 + 35500000.0*y5 + 13600000.0*y6 + 26000000.0*y8;
M[29] = 0.0158*y0 + 1270.0*y10 + 1270.0*y12 + 0.000155*y8 + 2.12e-5;
M[30] = -118000000.0*y11 - 24400.0*y4 - 13300000.0*y5 - 13300000.0*y6 - 13300000.0*y9;
M[31] = 128000.0*y12;
M[32] = 27600000.0*y0 - 24400.0*y2;
M[33] = 35500000.0*y0 - 13300000.0*y2 + 11000000.0*y8;
M[34] = 13600000.0*y0 + 4070000.0*y12 - 13300000.0*y2;
M[35] = 0;
M[36] = 26000000.0*y0 + 0.000155*y1 + 11000000.0*y5 + 101000.0*y9;
M[37] = -13300000.0*y2 + 101000.0*y8;
M[38] = 1270.0*y1 + 786000.0*y12;
M[39] = -118000000.0*y2;
M[40] = 1270.0*y1 + 786000.0*y10 + 128000.0*y3 + 4070000.0*y6;
M[41] = 0;
M[42] = 14520000.0*y0 + 27600000.0*y4;
M[43] = 4.58e-8*y4;
M[44] = 0;
M[45] = -128000.0*y12 - 39500.0*y5;
M[46] = 27600000.0*y0 + 4.58e-8*y1 + 10280000.0*y4;
M[47] = -39500.0*y3;
M[48] = 0;
M[49] = 0;
M[50] = 0;
M[51] = 0;
M[52] = 0;
M[53] = 0;
M[54] = -128000.0*y3;
M[55] = 0;
M[56] = 0.0158*y1 + 20900000.0*y11 - 27600000.0*y4;
M[57] = 0.0158*y0 - 4.58e-8*y4;
M[58] = -24400.0*y4;
M[59] = 128000.0*y12 + 39500.0*y5;
M[60] = -27600000.0*y0 - 4.58e-8*y1 - 24400.0*y2 - 20560000.0*y4 - 10900000.0*y5 - 36500.0*y6 - 13100000.0*y7 - 11400000.0*y8 - 11400000.0*y9 - 5.83;
M[61] = 39500.0*y3 - 10900000.0*y4;
M[62] = -36500.0*y4;
M[63] = -13100000.0*y4;
M[64] = -11400000.0*y4;
M[65] = -11400000.0*y4;
M[66] = 0;
M[67] = 20900000.0*y0;
M[68] = 128000.0*y3;
M[69] = 0;
M[70] = -35500000.0*y5 + 13600000.0*y6;
M[71] = 1270.0*y12 + 4.58e-8*y4;
M[72] = -13300000.0*y5;
M[73] = -39500.0*y5;
M[74] = 4.58e-8*y1 - 10900000.0*y5 + 36500.0*y6 + 22800000.0*y9;
M[75] = -35500000.0*y0 - 4070000.0*y10 - 13300000.0*y2 - 39500.0*y3 - 10900000.0*y4 - 19240000.0*y5 - 29200.0*y6 - 11000000.0*y8 - 8840000.0*y9 - 0.0943;
M[76] = 13600000.0*y0 + 36500.0*y4 - 29200.0*y5;
M[77] = 0;
M[78] = -11000000.0*y5;
M[79] = 22800000.0*y4 - 8840000.0*y5;
M[80] = -4070000.0*y5;
M[81] = 50200000.0*y12;
M[82] = 1270.0*y1 + 50200000.0*y11;
M[83] = 0;
M[84] = -13600000.0*y6 + 13000000.0*y8;
M[85] = 1270.0*y10;
M[86] = -13300000.0*y6;
M[87] = 0;
M[88] = -36500.0*y6;
M[89] = 9620000.0*y5 - 29200.0*y6;
M[90] = -13600000.0*y0 - 4070000.0*y12 - 13300000.0*y2 - 36500.0*y4 - 29200.0*y5 - 0.09430026;
M[91] = 0;
M[92] = 13000000.0*y0 + 101000.0*y9;
M[93] = 101000.0*y8 + 1680.0*y9;
M[94] = 1270.0*y1 + 50200000.0*y11;
M[95] = 50200000.0*y10;
M[96] = -4070000.0*y6;
M[97] = 0;
M[98] = -22900000.0*y7;
M[99] = 0;
M[100] = 0;
M[101] = 0;
M[102] = -13100000.0*y7;
M[103] = 11000000.0*y8 + 8840000.0*y9;
M[104] = 1.3e-7;
M[105] = -22900000.0*y0 - 3750000.0*y12 - 13100000.0*y4;
M[106] = 11000000.0*y5 + 101000.0*y9;
M[107] = 8840000.0*y5 + 101000.0*y8 + 1680.0*y9;
M[108] = 0;
M[109] = 0;
M[110] = -3750000.0*y7;
M[111] = 2620.0;
M[112] = 22900000.0*y7 - 13000000.0*y8;
M[113] = -0.000155*y8;
M[114] = 13300000.0*y9;
M[115] = 0;
M[116] = -11400000.0*y8;
M[117] = 4070000.0*y10 - 11000000.0*y8;
M[118] = 0;
M[119] = 22900000.0*y0;
M[120] = -13000000.0*y0 - 0.000155*y1 - 50200000.0*y11 - 11400000.0*y4 - 11000000.0*y5 - 101000.0*y9;
M[121] = 13300000.0*y2 - 101000.0*y8 + 773000.0;
M[122] = 786000.0*y12 + 4070000.0*y5;
M[123] = -50200000.0*y8;
M[124] = 786000.0*y10;
M[125] = 0;
M[126] = -13000000.0*y9;
M[127] = 0.000155*y8;
M[128] = -13300000.0*y9;
M[129] = 0;
M[130] = 13100000.0*y7 - 11400000.0*y9;
M[131] = 29200.0*y6 - 8840000.0*y9;
M[132] = 4070000.0*y12 + 29200.0*y5;
M[133] = 13100000.0*y4;
M[134] = 0.000155*y1 + 50200000.0*y11 - 101000.0*y9;
M[135] = -13000000.0*y0 - 13300000.0*y2 - 11400000.0*y4 - 8840000.0*y5 - 101000.0*y8 - 3360.0*y9 - 773000.0;
M[136] = 0;
M[137] = 50200000.0*y8;
M[138] = 4070000.0*y6;
M[139] = 0;
M[140] = 13000000.0*y9;
M[141] = -1270.0*y10;
M[142] = 13300000.0*y6;
M[143] = 0;
M[144] = 11400000.0*y8;
M[145] = -4070000.0*y10;
M[146] = 13300000.0*y2 + 0.0943;
M[147] = 0;
M[148] = 11400000.0*y4;
M[149] = 13000000.0*y0;
M[150] = -1270.0*y1 - 50200000.0*y11 - 786000.0*y12 - 4070000.0*y5;
M[151] = -50200000.0*y10;
M[152] = -786000.0*y10;
M[153] = 0;
M[154] = -20900000.0*y11;
M[155] = 2.12e-5;
M[156] = -118000000.0*y11;
M[157] = 0;
M[158] = 5.83;
M[159] = 0.0943;
M[160] = 0.0943;
M[161] = 0;
M[162] = -50200000.0*y11;
M[163] = 773000.0;
M[164] = -50200000.0*y11;
M[165] = -20900000.0*y0 - 50200000.0*y10 - 50200000.0*y12 - 118000000.0*y2 - 50200000.0*y8;
M[166] = -50200000.0*y11;
M[167] = 0;
M[168] = 0;
M[169] = -1270.0*y12;
M[170] = 13300000.0*y5;
M[171] = -128000.0*y12;
M[172] = 0;
M[173] = 13300000.0*y2 + 0.0943;
M[174] = -4070000.0*y12;
M[175] = -3750000.0*y12;
M[176] = 0;
M[177] = 0;
M[178] = -786000.0*y12;
M[179] = -50200000.0*y12;
M[180] = -1270.0*y1 - 786000.0*y10 - 50200000.0*y11 - 128000.0*y3 - 4070000.0*y6 - 3750000.0*y7;
M[181] = 2620.0;
M[182] = 0;
M[183] = 0;
M[184] = 0;
M[185] = 0;
M[186] = 0;
M[187] = 0;
M[188] = 0;
M[189] = 3750000.0*y12;
M[190] = 0;
M[191] = 0;
M[192] = 0;
M[193] = 0;
M[194] = 3750000.0*y7;
M[195] = -2620.0;

Now create a custom printer that uses cse() on the two matrices simulatneously so that subexpressions are not repeated. Hint: think about how the list printer method, _print_list(self, list_of_exprs), might help here.

Double Click For Solution

In [36]:
# write your answer here

Bonus Exercise: Compile and Run the C Program

Below we provide you with a template for the C program described above. You can use it by passing in a string like:

c_template.format(code='the holy grail')

Use this template and your code printer to create a file called run.c in the working directory.

To compile the code there are several options. The first is gcc (the GNU C Compiler). If you have Linux, Mac, or Windows (w/ mingw installed) you can use the Jupyter notebook ! command to send your command to the terminal. For example:

!gcc run.c -lm -o run

This will compile run.c, link against the C math library with -lm and output, -o, to a file run (Mac/Linux) or run.exe (Windows).

On Mac and Linux the program can be executed with:

!./run

and on Windows:

!run.exe

Other options are using the clang compiler or Windows cl compiler command:

!clang run.c -lm -o run
!cl run.c -lm

Double Click For Solution

In [37]:
c_template = """\
#include <math.h>
#include <stdio.h>

void evaluate_odes(const double state_vals[14], double rhs_result[14], double jac_result[196])
{{
    // We need to fill in the code here using SymPy.
{code}
}}

int main() {{

    // initialize the state vector with some values
    double state_vals[14] = {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14}};
    // create "empty" 1D arrays to hold the results of the computation
    double rhs_result[14];
    double jac_result[196];

    // call the function
    evaluate_odes(state_vals, rhs_result, jac_result);

    // print the computed values to the terminal
    int i;
    printf("The right hand side of the equations evaluates to:\\n");
    for (i=0; i < 14; i++) {{
        printf("%lf\\n", rhs_result[i]);
    }}
    printf("\\nThe Jacobian evaluates to:\\n");
    for (i=0; i < 196; i++) {{
        printf("%lf\\n", jac_result[i]);
    }}

    return 0;
}}\
"""
In [38]:
# write your answer here