Theano + Sympy for system of ODEs: Mapping dictionary issue and a few questions

1,085 views
Skip to first unread message

Guy Parsey

unread,
Aug 5, 2013, 3:13:59 PM8/5/13
to sy...@googlegroups.com
Hello Everyone,
Thank you in advance for reading through my problem and for any input you may have. I must say that I still feel like a novice programmer and my problems may be easily solvable from a different mindset. My present project entails time-integration of extremely stiff and non-linear ODEs with regards to chemical kinetics (one derivative equation for each variable species) and energy equations. Initially we were planning on using the sympy.lambdify function to create callable functions for the main function along with the jacobian and passing said functions to scipy.integrate.odeint, however this method only works for easier test cases (fewer species and/or no energy equations) before being limited by either the list recursion limit or segfaulting due to the limited stack size. I know that both of these limits can be edited, but that fact that I am reaching them makes me feel as though I am doing something extremely inefficiently. Outside of the documentation of SymPy and Theano, I have also been heavily using the BlogPost by Matthew Rocklin http://matthewrocklin.com/blog/work/2013/03/19/SymPy-Theano-part-1/

Presently I am trying to use the mapping between Theano and SymPy (sympy.printing.theanocode theano_function) to make my callable functions and take advantage of the optimization routines. I have two major problems and a few questions:

1st major problem: Though piecewise functions exist in SymPy (sympy.functions.elementary.piecewise) there is no counterpart in Theano. Looking at the source of the inspiration for theanocode (https://github.com/nouiz/theano_sympy/    graph_translation.py) I see that some of the SymPy equivalents were defined as lambda functions. Is there an equivalent way to add Theano conditional expressions wrapped into a function to add to the mapping dictionary in theanocode.py?

2nd major problem: Similar to the problem above in that I am not sure that the Theano counterpart is; some of the terms that I use are interpolated functions (with one ODE variable as input) that we have wrapped symbolically while providing a numerical implementation (so that symbolic derivatives can be made, resulting in their own interpolations). Is it possible to recreate the interpolation function as a Theano operation for use within the system of ODEs? 

Remain questions:
I presently have to flatten my input to theano_function to a list of expressions and then wrap to return to a form (Jacobian is a matrix not a vector); is it possible to have a matrix of different expressions as an input to theano_function with a vector output?

I know that a huge amount of Theano speed up is due to parallelization of matrix operations (which I do not have), should I be focusing on SymPy Autowrap/Ufuncify or my own code generation instead of trying to get Theano to play nicely?

Stupid questions:
Does sympy.printing.theanocode.theano_function automatically optimize the compiled graph?

Minor comment:
Perhaps unnecessary for most uses of the theano_function, but I needed to modify function inputs so as to be able to use the keyword argument 'on_unused_input=ignore' as opposed to 'raise' so that I did not need to have all symbols in all equations. This may be avoided by having the unused symbols somehow (I don't know how) included in each expression.

Thank you again for your time in reading my problems and any potential help you may think of. I can attach code if necessary, I just didn't want to make my post more confusing.
Have an excellent day.
Sincerely,
Guy Parsey

Jason Moore

unread,
Aug 5, 2013, 9:03:24 PM8/5/13
to sy...@googlegroups.com
Guy,

We're working on the same problem for sympy.physics.mechanics. Matthew Rocklin added support for matrix conversions in the theano code that is in SymPy and I used that, but found that theano was slower that lambdify for most of my cases (I only have two cores, so I'm not taking advantage of the Theano parallel stuff). I think writing specific code gen for ode integration is going to be the best bet. I'm happy to collaborate on this with you.

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

Frédéric Bastien

unread,
Aug 6, 2013, 9:25:50 AM8/6/13
to sympy
Hi,

I don't know what is sympy.functions.elementary.piecewise. Do Jason answered that part? If not, I'll look into it to know how to make Theano reproduce it. About converting any Sympy symbol to Theano symbol, when there  isn't a one to one matching, you can create a one to a full Theano graph conversion. When this is possible, it is probably the best, as if you make a new Theano op, it work work on the GPU. But if you make a Theano graph, there is good change that the graph will already work on the GPU.

If it is not possible to make a Theano graph for a sympy symbol, it is possible to make a new Theano op that just wrap the sympy implementation. Also, if this is a bottleneck, I recently added an example that show how to use numba with Theano so speed up the python code in a Theano op.


Theano do not parallelize on the CPU, except for the call to BLAS, when the BLAS library is parallel. On the GPU, it is parallel.

Having on_unused_input=ignore is normal for complicated generated code. That is why it was added. But when the code is simpler and not generated, but all user manually coded, most of time it mean the user didn't do what he wanted. If you know it is normal that you have unused input, there is no problem to use that flag.


Jason, about the case where Theano is slower, can you send me the Theano code? I would like to look at it. I'm very surprised that Theano is slower then Sympy in this case and would like to know why it is like this.

Fred

Jason Moore

unread,
Aug 6, 2013, 9:37:41 AM8/6/13
to sy...@googlegroups.com
Fred,

I think on_used_input=ignore should be a default in Matthew's theano printing code or that arg needs to be pushed up to his layer. I hit that issue too.

The code I have is here: https://github.com/PythonDynamics/pydy-code-gen

See the results.txt file for basic speed comparisons. I'm generating the ODEs for an n-link pendulum with mechanics and then see how fast it generates and simulates with scipy.odeint.

The code doesn't work at the moment. I haven't touched it in a month and looks like some things have changed in sympy and/or theano. I'll work on the bugs now.

But "python benchmark.py" should run it with sympy master and ?some? version of Theano.

Jason Moore

unread,
Aug 6, 2013, 10:39:57 AM8/6/13
to sy...@googlegroups.com
You'll need this patch: https://github.com/sympy/sympy/pull/2358 to run my code and you need theano latest dev version.

Guy Parsey

unread,
Aug 6, 2013, 11:04:33 AM8/6/13
to sy...@googlegroups.com
Hello,

Theano/Sympy questions
Fred, in terms of replacing sympy.Pieceiwse with a Theano equivalent, since sympy.Piecewise attempts each condition (from ExprCondPair) until one is valid, I would think that the closest equivalent would be a recursive theano ifelse (ie. sympy.Piecewise((expr1,cond1),(expr2,cond2)) ~ theano.ifelse(cond1,expr1,theano.ifelse(cond2,expr2,None))) statement so that the expressions are not evaluated until the condition is achieved. This might be too much of a patch though and I am not sure how to implement it with the same argument structure as sympy.Pieceiwse.

With regards to using an externally defined theano graph or op (theano wrapped sympy implementation), how would one  pass them to theano_function? It feels as though it would be analogous to the autowrap helpers argument, or would this be an issue of merging graphs before calling theano.function? If I can understand this step, I feel as though both of my problems would be solvable (a wrapping of a piecewise along with a wrapping of an interpolated function)

In all of my simplified test cases (no pieceiwse or interpolations, 9 heavily-linked non-linear ODEs), I also show that the lambdify function is faster to evaluate than my created theano function. When passing a list of expressions, if there is a commons subexpression that exists in two expressions (not twice in the same expression), is it treated as such during theano compilation or is each expression handled separately?

Codegen questions (maybe a different topic):
Jason, I have started playing with the sympy.printing.codegen (C for now, Fortran later) as an alternative to theano (I would very much like to get multiple methods working, and GPU acceleration makes me want to keep theano implementation). There is a size-able speedup relative to lambdify or theano_function, which is to be expected, but I arrive at the same problem as above with regards to interpolated or piecewise functions. The initial crash from using Piecewise functions comes from Routine calling sympy.tensor.index_methods get_contraction_structure which explicitly states no support for Piecewise function types. In principle, if one made a symbolic implemented function (preferably without an analytical representation), could one link the python object to the compiled code? What attributes would said implemented_function need to have?

Thank you for your responses, I am glad to hear that at least I am not completely missing something.
Cheers,
Guy

Guy Parsey

unread,
Aug 6, 2013, 11:30:22 AM8/6/13
to sy...@googlegroups.com
Addon to the Codegen section:
With the example above of having an implemented function in some expression in sympy, assuming one can write the c-code equivalent, how would one link the c-code to the codegen process as opposed to having codegen generate it?
Cheers,
Guy

Jason Moore

unread,
Aug 6, 2013, 11:34:01 AM8/6/13
to sy...@googlegroups.com
Guy,

I don't know enough about the details of the codegen module, but theoretically there should be a mapping from sympy functions to a c-code equivalent that you can add and/or override. It is just a matter of getting real familiar with the code gen modules.

You may want to check out https://github.com/IgnitionProject/ignition. Andy showed me this at SciPy and setting up your own mappings for various things seemed more clean than what is in SymPy.

Frédéric Bastien

unread,
Aug 6, 2013, 11:53:58 AM8/6/13
to sympy
About the benchmark, the problem is a known one, Theano isn't very fast to parse Theano input. In your case, it is this parsing that take the most time in theano when n is small. But when n increase, Theano is faster then lambdify:

n = 1
The derivation took 0.154258966446 seconds.
Running with theano method.
It took 0.515482187271 seconds to compute M and F with theano 1000 times at an average of 0.000515482187271 seconds per computation.
Running with lambdify method.
It took 0.0773530006409 seconds to compute M and F with lambdify 1000 times at an average of 7.73530006409e-05 seconds per computation.
n = 6
The derivation took 6.48748397827 seconds.
Running with theano method.
It took 1.83713889122 seconds to compute M and F with theano 1000 times at an average of 0.00183713889122 seconds per computation.
Running with lambdify method.
It took 0.807645082474 seconds to compute M and F with lambdify 1000 times at an average of 0.000807645082474 seconds per computation.
n = 11
The derivation took 33.0718569756 seconds.
Running with theano method.
It took 3.446434021 seconds to compute M and F with theano 1000 times at an average of 0.003446434021 seconds per computation.
Running with lambdify method.
It took 3.89024996758 seconds to compute M and F with lambdify 1000 times at an average of 0.00389024996758 seconds per computation.
n = 16
The derivation took 104.534656048 seconds.
Running with theano method.
It took 5.41041302681 seconds to compute M and F with theano 1000 times at an average of 0.00541041302681 seconds per computation.
Running with lambdify method.
It took 11.5178751945 seconds to compute M and F with lambdify 1000 times at an average of 0.0115178751945 seconds per computation.


The option to make it faster would be to give less inputs to the theano function. The first step would be to make the constants as theano constant instead of input to the theano function. That will make put the constant in the c code, so it could also speed up even the computation itself.

Jason Moore

unread,
Aug 6, 2013, 12:00:27 PM8/6/13
to sy...@googlegroups.com
How do you specify which values are constants and which aren't?

I do pass in a lot of constants and they should be hard coded in C if possible (or some other way to make them available).

Frédéric Bastien

unread,
Aug 6, 2013, 12:21:09 PM8/6/13
to sympy
You can just use the corresponding ndarray with its value. Theano will convert it to a constant. To explicitly make a constant, you can call theano.tensor.constant(a_ndarray_or_python_object).

I don't know enought about theano_function(), but can you make a SymPy constant and maybe it will convert it to a Theano constant?

Frédéric Bastien

unread,
Aug 6, 2013, 1:28:16 PM8/6/13
to sympy
Hi,

I used 2 "tricks" to speed up the Theano version. I made a PR to your repo with this. First Theano is strongly typed. You where passed numpy scalar (numpy.float64) as input, but the theano variable where representing numpy.ndarray. Changing the input to pass numpy.ndarray (of 0 dimensions) changed the Theano run time from 0.5s to 0.15s. Then using a special attribute a_theano_function.trust_input=True, this lower to 0.08s for the run time. This is with n=1. Here is the new benchmark results:

n = 1
The derivation took 0.142597913742 seconds.
Running with lambdify method.
It took 0.0880799293518 seconds to compute M and F with lambdify 1000 times at an average of 8.80799293518e-05 seconds per computation.
Running with theano method.
It took 0.0840051174164 seconds to compute M and F with theano 1000 times at an average of 8.40051174164e-05 seconds per computation.
n = 4
The derivation took 2.22830605507 seconds.
Running with lambdify method.
It took 0.745465040207 seconds to compute M and F with lambdify 1000 times at an average of 0.000745465040207 seconds per computation.
Running with theano method.
It took 0.180059194565 seconds to compute M and F with theano 1000 times at an average of 0.000180059194565 seconds per computation.
n = 7
The derivation took 8.44271087646 seconds.
Running with lambdify method.
It took 2.92468690872 seconds to compute M and F with lambdify 1000 times at an average of 0.00292468690872 seconds per computation.
Running with theano method.
It took 0.314508914948 seconds to compute M and F with theano 1000 times at an average of 0.000314508914948 seconds per computation.
n = 10
The derivation took 22.0363950729 seconds.
Running with lambdify method.
It took 7.42660617828 seconds to compute M and F with lambdify 1000 times at an average of 0.00742660617828 seconds per computation.
Running with theano method.
It took 0.471493005753 seconds to compute M and F with theano 1000 times at an average of 0.000471493005753 seconds per computation.

With thoses 2 tricks, it solve the Theano overhead for the inputs as now Theano is always faster then lambdify. But yes, ideally Theano should do this automatically or warn about the bad input type that get automatically converted.

Frédéric Bastien

unread,
Aug 6, 2013, 1:36:29 PM8/6/13
to sympy
On Tue, Aug 6, 2013 at 11:04 AM, Guy Parsey <guy.p...@gmail.com> wrote:
Hello,

Theano/Sympy questions
Fred, in terms of replacing sympy.Pieceiwse with a Theano equivalent, since sympy.Piecewise attempts each condition (from ExprCondPair) until one is valid, I would think that the closest equivalent would be a recursive theano ifelse (ie. sympy.Piecewise((expr1,cond1),(expr2,cond2)) ~ theano.ifelse(cond1,expr1,theano.ifelse(cond2,expr2,None))) statement so that the expressions are not evaluated until the condition is achieved. This might be too much of a patch though and I am not sure how to implement it with the same argument structure as sympy.Pieceiwse.

It would work. But for the timming, I would also try with the switch op instead of ifelse. The ifelse currently have much bigger oveahead then the switch op, and even if the switch op evaluate both branch, it happen relatively frequently that it is faster then the ifelse. Especially if each expression have small size, like a scalar or a small vector.
 
With regards to using an externally defined theano graph or op (theano wrapped sympy implementation), how would one  pass them to theano_function? It feels as though it would be analogous to the autowrap helpers argument, or would this be an issue of merging graphs before calling theano.function? If I can understand this step, I feel as though both of my problems would be solvable (a wrapping of a piecewise along with a wrapping of an interpolated function)

There is 2 steps. 1) make a theano op like descrided here: http://deeplearning.net/software/theano/tutorial/extending_theano.html

2) Modify theano_function to use the new Theano op for the equivalent SymPy op. It could be possible to make an autowrap for this, but I don't know enough of the detail to be sure.
 
In all of my simplified test cases (no pieceiwse or interpolations, 9 heavily-linked non-linear ODEs), I also show that the lambdify function is faster to evaluate than my created theano function. When passing a list of expressions, if there is a commons subexpression that exists in two expressions (not twice in the same expression), is it treated as such during theano compilation or is each expression handled separately?

Theano will merge duplicated subexpression automatically by default. What is your function? If it is too "simple", for example just working on scalar, Theano have by default much overhead. As said for the other benchmark, there is tricks to remove a big part or even all of it, but it ask for some work on the developer.

Matthew Rocklin

unread,
Aug 6, 2013, 5:48:13 PM8/6/13
to sy...@googlegroups.com

Thanks Guy and Jason for starting this conversation and thanks Fred for jumping in with Theano expertise.

I'm out camping in Wisconsin this week but should be available starting Thursday.  I'm excited to help contribute to this machinery.

Guy Parsey

unread,
Aug 7, 2013, 11:37:16 AM8/7/13
to sy...@googlegroups.com, no...@nouiz.org
Hey Fred,
I am starting to play with making a theano.Op, and in principle, linking a custom theano.Op to an implemented sympy function will solve all of my problems (at least in the sense of getting my system running, not necessarily speed). My initial post spoke of both mapping Piecewise and interpolated functions into theano; after speaking with my advisor, due to wanting a smoothing criterion over our piecewise datasets, I think we will end up treating them as interpolated functions over smoothed data sets. My idea for how to use an interpolation routine will be a wrapped theano.Op, so I am still rather curious as to how implementing Piecewise would work in such way that theano can manipulate it in graph form (built out of native theano operations).

Theano will merge duplicated subexpression automatically by default. What is your function? If it is too "simple", for example just working on scalar, Theano have by default much overhead. As said for the other benchmark, there is tricks to remove a big part or even all of it, but it ask for some work on the developer.

The basic system/functions can be described as follows; we have some number of species (n_0 ...n_k) and we are creating a function for the vector of time derivatives. The different types of expressions come about from the reaction rates (K_0 ..... K_l) The whole system is described symbolically, for symbolic determination of the jacobian (spline interpolated functions are wrapped with a sympy implemented function so that taking a symbolic derivative also yields the derivative of the numerical implementation). Dummy example below with three types of base expressions
d(n_i)/dt = K_0*n_2*n_0 + K_1*n_3 + K_2*n_4*n_2*n_4/(sum(n_0....n_k))
Each of the K terms can be described differently:
--constant
--analytical parameterized of some n_j    (Kf ~ a0*(n_j**a1)*exp(a2/n_j)  , a0..a2 are constants)
--spline interpolation function of some n_j  (K_interp(n_j))
--analytical weighted with piecewise/interpolated (new_Kf ~ Kf*exp(-K_interp(n_j))
I also have two energy equations, but these incorporate all of the same types of terms, just not in a nice summation as in the species equations. Any given base expression (eg. K_0*n_2*n_0) will appear in multiple time derivative equations, so if the interpolation is slow, it would be nice to get theano to recognize that it is a 'common expression.' I can make test cases with a few species and a few reactions (Ks), but the true systems I would like to be simulating will have 40-100 species and 300-2000 reaction rates. Are these functions too 'simple' to be sped by using theano? There is a lot of optimization that can be done due to the high number of shared expressions, so I was hoping to take advantage of the theano graph optimizations.

My initial test cases were using the numpy.float64 dtype, so I will switch it over to ndarrays if possible.
I hope to have questions or a working spline interpolation theano.Op in the next few hours.
Cheers,
Guy

Quick Codegen question: In principle, can one create a python callback within a compiled c code to a python function that is not compiled (interpolation function)?

Matthew Rocklin

unread,
Aug 8, 2013, 1:19:26 PM8/8/13
to sy...@googlegroups.com
Ok, Just went through this thread again. 

It should be simple to translate SymPy.Piecewise to a recursive Theano.switch (after translating SymPy.LT to theano.lt, etc.)  I'll get on this soon.  Does this sound reasonable to you Fred?

Jason's PR allowing pass through of keyword arguments looks close to me.  I suspect that the `on_unused_input=ignore` issue will soon be possible.  Question, should this be default?  Are there other defaults that should be used?  I'm curious why you both are making functions with unused inputs. 

Jason, was Fred's note on directly passing constants useful?  If not why not?  It sounds like you should subs your constant symbols for numeric values prior to creating the theano graph.

With regards to using an externally defined theano graph or op (theano wrapped sympy implementation), how would one  pass them to theano_function?
You wouldn't.  You would use the other common interface function, theano_code to create theano variables instead.  You would then use raw theano operations (like theano.function) to compile your graph from theano variables.  Theano_function is a convenient all-in-one interface function for the common case.  If you want to do anything custom I recommend using SymPy and Theano natively, translating from one to the other using theano_code instead.  It does nothing except translate.

Does sympy.printing.theanocode.theano_function automatically optimize the compiled graph?
Yes.  theano_function translates sympy expressions into theano expressions using theano_code.  It then calls theano.function on these expressions.  theano.function optimizes by default. 

> SymPy C Codegen and Theano

@Fred, how hard would it be to leverage SymPy's C codegen in Theano?  This might be a lot cleaner than wrapping raw SymPy operations and might substantially extend Theano's support of scalar expressions.  Do you have a performant Bessel function op?  I'll bet SymPy could be made to do this quite well.

@Aaron / @Ondrej, if you're reading this thread could you point us to the best place to start looking at C codegen in SymPy?  Alternatively can you point to an active community member who would be able to do so?


> ODE Integration

What are people's thoughts on scipy.integrate.odeint?  Earlier this year I was looking at nodepy, a purely mathematical treatment/catalog of time stepping schemes.  It should be easy to generate such schemes in Theano using simple abstractions like Butcher matrices.  The code I was playing around with was https://github.com/mrocklin/tsgen.  The goal was that we could easily experiment with and rapidly benchmark different time stepping schemes for specific problems.


Aaron Meurer

unread,
Aug 8, 2013, 1:24:00 PM8/8/13
to sy...@googlegroups.com
Maybe Ondrej knows, but I don't. I would just look at the code. I
believe the meat of it is in the printers, while the rest is just nice
wrappers.

Aaron Meurer

Jason Moore

unread,
Aug 8, 2013, 1:37:34 PM8/8/13
to sy...@googlegroups.com
On Thu, Aug 8, 2013 at 1:19 PM, Matthew Rocklin <mroc...@gmail.com> wrote:
Ok, Just went through this thread again. 

It should be simple to translate SymPy.Piecewise to a recursive Theano.switch (after translating SymPy.LT to theano.lt, etc.)  I'll get on this soon.  Does this sound reasonable to you Fred?

Jason's PR allowing pass through of keyword arguments looks close to me.  I suspect that the `on_unused_input=ignore` issue will soon be possible.  Question, should this be default?  Are there other defaults that should be used?  I'm curious why you both are making functions with unused inputs.

If I generate tons of long ode's I don't really care to search through them to find all of the unique arguments for each equation. We setup a multibody problem with lots of constants and lots of time varying state variables (+ others too). Every ode generated will only use a subset of these variables for the whole system. So it is just easier to pass in all the numerical values it could possibly need instead of trying to figure out which ones it does need. But we could surely be more precise about this, it's just easier to pass them all in.
 
 

Jason, was Fred's note on directly passing constants useful?  If not why not?  It sounds like you should subs your constant symbols for numeric values prior to creating the theano graph.

I'll try that, but subs can be time consuming for our bigger problems.

Matthew Rocklin

unread,
Aug 8, 2013, 1:44:05 PM8/8/13
to sy...@googlegroups.com
> I'll try that, but subs can be time consuming for our bigger problems.

Try expr.xreplace(subs_dict)

Jason Moore

unread,
Aug 8, 2013, 2:03:57 PM8/8/13
to sy...@googlegroups.com
Actually subbing in the constants before the code gen is a problem. They are constant in the sense that they don't change with respect to calling the theano generated function for a given simulation. But if I wanted to change the mass constant of one of the objects in my system it would suck to have to regenerate a theano function just because I want to change the numerical value of one of these constants.

So I can have the constant values fixed with respect to one odeint call but want to be able to change them for the next simulation without having to regenerate a theano function. I think Luke solves this in generated C++ code with a static data member. Is there there a similar concept for theano generated code?

Matthew Rocklin

unread,
Aug 8, 2013, 2:07:36 PM8/8/13
to sy...@googlegroups.com
When I asked this long ago [1] the solution came back as shared variables [2]

Jason Moore

unread,
Aug 8, 2013, 2:23:47 PM8/8/13
to sy...@googlegroups.com
Cool, I'll look into it.

Ondřej Čertík

unread,
Aug 8, 2013, 3:28:18 PM8/8/13
to sympy
C code gen is here:

https://github.com/sympy/sympy/blob/master/sympy/utilities/codegen.py

Is this what you were looking for?

Ondrej

Matthew Rocklin

unread,
Aug 8, 2013, 4:04:32 PM8/8/13
to sy...@googlegroups.com
After looking around I think I'm looking for the codegen routine.


Frédéric Bastien

unread,
Aug 8, 2013, 4:15:37 PM8/8/13
to sympy
On Wed, Aug 7, 2013 at 11:37 AM, Guy Parsey <guy.p...@gmail.com> wrote:
Hey Fred,
I am starting to play with making a theano.Op, and in principle, linking a custom theano.Op to an implemented sympy function will solve all of my problems (at least in the sense of getting my system running, not necessarily speed). My initial post spoke of both mapping Piecewise and interpolated functions into theano; after speaking with my advisor, due to wanting a smoothing criterion over our piecewise datasets, I think we will end up treating them as interpolated functions over smoothed data sets. My idea for how to use an interpolation routine will be a wrapped theano.Op, so I am still rather curious as to how implementing Piecewise would work in such way that theano can manipulate it in graph form (built out of native theano operations).

Theano will merge duplicated subexpression automatically by default. What is your function? If it is too "simple", for example just working on scalar, Theano have by default much overhead. As said for the other benchmark, there is tricks to remove a big part or even all of it, but it ask for some work on the developer.

The basic system/functions can be described as follows; we have some number of species (n_0 ...n_k) and we are creating a function for the vector of time derivatives. The different types of expressions come about from the reaction rates (K_0 ..... K_l) The whole system is described symbolically, for symbolic determination of the jacobian (spline interpolated functions are wrapped with a sympy implemented function so that taking a symbolic derivative also yields the derivative of the numerical implementation). Dummy example below with three types of base expressions
d(n_i)/dt = K_0*n_2*n_0 + K_1*n_3 + K_2*n_4*n_2*n_4/(sum(n_0....n_k))
Each of the K terms can be described differently:
--constant
--analytical parameterized of some n_j    (Kf ~ a0*(n_j**a1)*exp(a2/n_j)  , a0..a2 are constants)
--spline interpolation function of some n_j  (K_interp(n_j))
--analytical weighted with piecewise/interpolated (new_Kf ~ Kf*exp(-K_interp(n_j))
I also have two energy equations, but these incorporate all of the same types of terms, just not in a nice summation as in the species equations. Any given base expression (eg. K_0*n_2*n_0) will appear in multiple time derivative equations, so if the interpolation is slow, it would be nice to get theano to recognize that it is a 'common expression.' I can make test cases with a few species and a few reactions (Ks), but the true systems I would like to be simulating will have 40-100 species and 300-2000 reaction rates. Are these functions too 'simple' to be sped by using theano? There is a lot of optimization that can be done due to the high number of shared expressions, so I was hoping to take advantage of the theano graph optimizations.

Theano will automatically execute only once each common exact sub expression. This is something sure. What is not sure is not sure is common "equivalent" but not exact sub expression. My guess with the result of the benchmark here is that Theano should be faster then lambdify. But only benchmark can confirm this.
 
My initial test cases were using the numpy.float64 dtype, so I will switch it over to ndarrays if possible.
I hope to have questions or a working spline interpolation theano.Op in the next few hours.
Cheers,
Guy

Quick Codegen question: In principle, can one create a python callback within a compiled c code to a python function that is not compiled (interpolation function)?

From C code, you can execute a pure python function. Theano do this for op that have only a pure python code without any equivalent c code.

Frédéric Bastien

unread,
Aug 8, 2013, 4:30:11 PM8/8/13
to sympy
On Thu, Aug 8, 2013 at 1:19 PM, Matthew Rocklin <mroc...@gmail.com> wrote:
Ok, Just went through this thread again. 

It should be simple to translate SymPy.Piecewise to a recursive Theano.switch (after translating SymPy.LT to theano.lt, etc.)  I'll get on this soon.  Does this sound reasonable to you Fred?

It sound reasonable and is the first thing I suggest to try.
 
Jason's PR allowing pass through of keyword arguments looks close to me.  I suspect that the `on_unused_input=ignore` issue will soon be possible.  Question, should this be default?  Are there other defaults that should be used?  I'm curious why you both are making functions with unused inputs. 

Jason, was Fred's note on directly passing constants useful?  If not why not?  It sounds like you should subs your constant symbols for numeric values prior to creating the theano graph.

I don't know how frequently it would happen that this option is needed. It was discussed that pylearn2 always pass it. I don't know if he do it or not, but if the "framework" generate unseless input and that this is independent of the user, my guess is that this should be the default. But if this is dependent of the user, I think it should not be the default.
 
> SymPy C Codegen and Theano

@Fred, how hard would it be to leverage SymPy's C codegen in Theano?  This might be a lot cleaner than wrapping raw SymPy operations and might substantially extend Theano's support of scalar expressions.  Do you have a performant Bessel function op?  I'll bet SymPy could be made to do this quite well.

@Aaron / @Ondrej, if you're reading this thread could you point us to the best place to start looking at C codegen in SymPy?  Alternatively can you point to an active community member who would be able to do so?


@Matt, you already did a new Theano op with C code. I think it is the only "easy" way to wrap other people c code in Theano. If the person already know this C code AND a little of Python AND NumPy C-API, it isn't very hard to a new Theano op with C code. Otherwise, doing the first such op ask to learn a few think and could ask a few days. You already did this, so you have a good idea of the work it need.

Now the questions is how is done the SymPy code gen? Is just just string template that is filled with dtype and other stuff? If we can just call one SymPy function with the information of what we want and it return a string with the C code it could be relatively easy. The only questions is about how to handle the variable name to pass the information around. At worst, we wrap the sympy c code in a c function, then make a small wrapper c code that take the Theano c variable name and call this function. So not very hard as Theano provide what is needed.

 
> ODE Integration

What are people's thoughts on scipy.integrate.odeint?  Earlier this year I was looking at nodepy, a purely mathematical treatment/catalog of time stepping schemes.  It should be easy to generate such schemes in Theano using simple abstractions like Butcher matrices.  The code I was playing around with was https://github.com/mrocklin/tsgen.  The goal was that we could easily experiment with and rapidly benchmark different time stepping schemes for specific problems.

Honestly, I don't know what ode is and I didn't understand most of this paragraph due to missing knowledge. But I can tell that there isn't any "butcher" op in Theano. Maybe it could be done with a graph of the existing op, but I can't tell. At worst, we can make a new Theano op for that.

Frédéric Bastien

unread,
Aug 8, 2013, 4:38:11 PM8/8/13
to sympy
Yes and not. Using shared variable or explicit input to a Theano function is equivalent. We can't make more optimization with shared variable then with a Theano inputs. But shared variable allow more optimization when they are used with the "updates" mechanism: we specify the new value of a theano shared variable after the execution of the theano function. Theano can do the update inplace on the shared variable memory region, so this can save memory alloc and copy.

If you use a given set of constant value for only 1 theano function call, it is not a constant in theano mind set. Theano compilation is relatively slow, so to amortize this, you need to call the theano function many times. But if you use many times the same set of constant value, it could be faster to recompile a new theano function for each set of constant value. But this depend on many factor including the execution time of the function, the compilation time for that function and the number of time you reuse the function. So to know if this is usefull, you need to benchmark it. But if you don't spend at least a few minutes of run time execution per compilation, my guess is that it isn't worth to recompile. Or at least, it is not worth to spend time doing all this calculation/benchmark if you won't spend more time waiting for your execution then the time you spend testing this.

Frédéric Bastien

unread,
Aug 8, 2013, 4:40:24 PM8/8/13
to sympy
One last thing, theano.function() have a parameter givens, that allow to do substitution. That could be an easy way to don't change your Theano/SymPy graph, but just specify at compile time what is constant. But you will probably need to don't put those constant as input to the Theano function, which should be easy to do, just remove the constant value from the list.

Matthew Rocklin

unread,
Aug 8, 2013, 4:54:44 PM8/8/13
to sy...@googlegroups.com
It should be simple to translate SymPy.Piecewise to a recursive Theano.switch (after translating SymPy.LT to theano.lt, etc.)  I'll get on this soon.  Does this sound reasonable to you Fred?

It sound reasonable and is the first thing I suggest to try.

Working on this now.

> SymPy C Codegen and Theano

@Fred, how hard would it be to leverage SymPy's C codegen in Theano?  This might be a lot cleaner than wrapping raw SymPy operations and might substantially extend Theano's support of scalar expressions.  Do you have a performant Bessel function op?  I'll bet SymPy could be made to do this quite well.

@Aaron / @Ondrej, if you're reading this thread could you point us to the best place to start looking at C codegen in SymPy?  Alternatively can you point to an active community member who would be able to do so?


@Matt, you already did a new Theano op with C code. I think it is the only "easy" way to wrap other people c code in Theano. If the person already know this C code AND a little of Python AND NumPy C-API, it isn't very hard to a new Theano op with C code. Otherwise, doing the first such op ask to learn a few think and could ask a few days. You already did this, so you have a good idea of the work it need.

Now the questions is how is done the SymPy code gen? Is just just string template that is filled with dtype and other stuff? If we can just call one SymPy function with the information of what we want and it return a string with the C code it could be relatively easy. The only questions is about how to handle the variable name to pass the information around. At worst, we wrap the sympy c code in a c function, then make a small wrapper c code that take the Theano c variable name and call this function. So not very hard as Theano provide what is needed.

It looks like codegen is the relevant high-level api call 

In [1]: from sympy.utilities.codegen import
 
In [2]: expr = sin(x)**2

In [3]: [(c_name, c_code), (h_name, c_header)] = codegen(("f", expr), 'C', 'test', header=False)

In [4]: print c_code
#include "test.h"
#include <math.h>

double f(double x) {

   return pow(sin(x), 2);

}

The work I did rarely dealt with making and using functions.  I'll go over past work and see what I can do.  Expect some calls for help though!

Matthew Rocklin

unread,
Aug 8, 2013, 5:00:49 PM8/8/13
to sy...@googlegroups.com

Aaron Meurer

unread,
Aug 8, 2013, 5:05:08 PM8/8/13
to sy...@googlegroups.com
If you figure out how to do something that's under-documented, please document it. Our code generation needs a lot more documentation, especially high-level narrative documentation. 

Aaron Meurer

Matthew Rocklin

unread,
Aug 8, 2013, 7:27:56 PM8/8/13
to sy...@googlegroups.com
@Fred, what easy things can we do on the Theano side to support performant scalar codes?  Your two tricks above were impressive.  How can we make this behavior default?  Should this happen on the SymPy side or the Theano side (Theano side seems better if there is a place for it).

Frédéric Bastien

unread,
Aug 8, 2013, 11:20:14 PM8/8/13
to sympy
For your codegen.py example, it is easy to reuse it manually in a Theano Op. A good example is to look at the GpuEye op (and just pretend the the GPU kernel is a CPU one). I don't know of a simpler CPU example that demonstrate that:

https://github.com/Theano/Theano/blob/master/theano/sandbox/cuda/basic_ops.py#L3421

You put the c function in c_support_code() function and call it with code returned by the c_code() fct.

If we want to automate that, we need a way to know the order of the parameter passed to the c function if there is more then 1 parameter. Is that possible? Ideally, we need to separate the include from the c function, but I think it work even if we don't.

About the Theano speed up, the problem is located in the parsing of the Theano fct inputs. This is done in pure python code. I think we can always use the trust_input=True argument. In the case there is an error, you will get an error message that is harder to understands, as it is an op that will complain that it have bad input. But you won't be easily able to know to witch input it correspond.

It could also happen that you don't pass a good input and there is no error. For example if it is used in an op with only python code. But in that case, the computation is probably done correctly, but with the input object received, not the one specified when compiling the theano function. I can't certify this always happen, but I would be surprised to see something else being done.

So if you want, just always add trust_input=True, but be ready to have uglier error message and remove the flag to have a better error message in case you have a problem.

Ideally, Theano should move that to C code, but I don't think it will happen shortly.

Matthew Rocklin

unread,
Aug 9, 2013, 10:21:04 AM8/9/13
to sy...@googlegroups.com
A draft for a SymPy CCode op is here https://github.com/Theano/Theano/pull/1486