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

[Guy Parsey]

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]

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.

Jason
moorepants.info
+01 530-601-9791

--
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]

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]

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
moorepants.info
+01 530-601-9791

[Jason Moore]

You'll need this patch: https://github.com/sympy/sympy/pull/2358 to run my code and you need theano latest dev version.

Jason
moorepants.info
+01 530-601-9791

[Guy Parsey]

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]

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]

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.

Jason
moorepants.info
+01 530-601-9791

[Frédéric Bastien]

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]

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).

Jason
moorepants.info
+01 530-601-9791

[Frédéric Bastien]

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]

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]

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]

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]

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]

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]

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]

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]

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

Try expr.xreplace(subs_dict)

[Jason Moore]

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?

Jason
moorepants.info
+01 530-601-9791

[Matthew Rocklin]

When I asked this long ago [1] the solution came back as shared variables [2]

[1] https://groups.google.com/forum/#!topic/theano-dev/8sOtgZl7L_I

[2] http://deeplearning.net/software/theano/library/compile/shared.html

[Jason Moore]

Cool, I'll look into it.

Jason
moorepants.info
+01 530-601-9791

[Ondřej Čertík]

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]

After looking around I think I'm looking for the codegen routine.

[Frédéric Bastien]

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]

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]

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]

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]

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]

PR for Theano-Piecewise https://github.com/sympy/sympy/pull/2366

[Aaron Meurer]

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]

@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]

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]

A draft for a SymPy CCode op is here https://github.com/Theano/Theano/pull/1486

[Matthew Rocklin]

Hi Guy,  

Thanks again for raising these issues.  This is the status as I see it

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?

This is in https://github.com/sympy/sympy/pull/2366 and should hopefully soon be merged.

 

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? 

I'm not sure I fully understand this piece.  Is it possible to define your interpolation scheme with lower-level operators such as with Piecewise and polynomials?  What is your interpolation scheme? 

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?

Jason answered this I think.  In the dev version theano_function should be able to take SymPy matrices as inputs.  See the description of this PR for an example https://github.com/sympy/sympy/pull/2013

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?

I find your engagement of the SymPy-Theano connection valuable and am happy to support it.  I think you'll find autowrap/ufuncify possibly more mature but less well supported.

 

Stupid questions:

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

Answered: yes

 

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.

Answered: In PR https://github.com/sympy/sympy/pull/2358

> 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)

Answered: use theano_code

> 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?

I don't think we can do this.  If you have C implementations of operations within SymPy that aren't yet supported we'd love help with this.  See https://github.com/sympy/sympy/pull/2369

In general let us know how you're doing with this project and what is keeping you from making progress.

-Matt

[Guy Parsey]

Hey Matt et al.,

Between the responses (and code updates) from Fred, Jason and yourself, all of my present questions have been answered. Now I just need to get more familiar with Theano. I was quick to assume that I would understand all of the printer code, but at least now it is all beginning to make sense. Thank you all for all of your help.

I'm not sure I fully understand this piece.  Is it possible to define your interpolation scheme with lower-level operators such as with Piecewise and polynomials?  What is your interpolation scheme? 

In most cases, the interpolated functions that we have are pre-computed convolutions of an arbitrary energy distribution function (analytic) and cross section (experimental/simulation data defined). In principle we could try and do some type of a polynomial fit, but our focus has been continuity of the ODE system and system Jacobian, so we have been using splrep and splrev from scipy.interpolate (also to avoid oscillations around the infection points of each reaction rate). Once we are integrating the ODE system, the source data for the interpolations is non-changing, so perhaps our choice of python splines is not very efficient. My mention of interpolating the piecewise functions was with regards to forcing continuity of the piecewise derivative at a transition point.

I am completely humbled by the amount of information provided in this topic, thank you again.

Cheers,

Guy

[Guy Parsey]

Hello again everyone,

I thought I understood everything I needed to implement a theano Op wrapping a scipy spline function through theanocode, but have been promptly proven wrong (and was on a small family vacation). 

Firstly, many thanks for the modifications done to theanocode to allow for Piecewise and Undefined functions. I feel as though everything is in place for me to solve my problem, but I am still either lacking or mis-undertsanding something with regards to mapping a custom theano Op through the theanocode.theano_function.

I have created a quick test case to show what I have understood to date which in my mind should have all the pieces necessary to function correctly. I have made a small git repository on GitHub in order to share this example and because I became fed up trying to figure out how to publicly share a BitBucket repository (academic license-where I am hosting my thesis project-which will be made public once functioning correctly).

https://github.com/gparsey/KGMindep_SympyTheanoOp

Running:

>>> ipython SymPy_Theano_KGM_indep.py 

Evaluates three test cases: f0) simple arithmetic operation, f1) sympy piecewise into theano and f2) sympy undefined function wrapped spline into theano using a custom theano Op

Third test case crashes with:

<<<MissingInputError: ('An input of the graph, used to compute TheanoInterpWrapOp.theanointerp(y), was not provided and not given a value', y)

I apologize in advance for: the verbosity of the test cases (trying to figure out what is happening within theanocode) using loc_theanocode (modified sympy.printing.theanocode), the novice nature of my code and whether I included correct references to the SymPy community. I am pretty sure that I am either missing something crucial or doing something silly. Any and all help/comments would be greatly appreciated.

Cheers,

Guy

[Matthew Rocklin]

Have you tested your interpolation op in isolation from SymPy?

A quick glance at the error (quick glance means I can easily be wrong) leads me to think that this particular issue is localized within the domain of Theano.  If this is the case then I recommend asking about your spline op on the thean...@googlegroups.com mailing list.

--

[Guy Parsey]

Hey Matt,

I am pretty sure that I have tested the theano Op separately from Sympy, but again, I am probably missing something silly.

After running the test cases, or evaluating 

  k = TestInterpOp()

  k.CreateSuite()

from the SymPy_Theano_KGM_indep.py file, one can run the following lines

In [2]: s,K,Kp = k.SymInterp()

In [3]: op = k.TheanoInterpOp()

In [4]: x = theano.tensor.dvector()

In [5]: f = theano.function([x],op(x))

In [6]: s(4.0)

Out[6]: array(16.0)

In [7]: f([4.0])

Out[7]: array([ 16.])

where 's' is the sympy wrapped spline function (undefined function) and 'f' is the theano.function of the theano op created around the spline. 

Is this what you mean?

Cheers,

Guy

[Guy Parsey]

Correction (independent of testing the theano op), 's' is a simple wrapper to the spline function, it is not the sympy wrapper. the sympy wrapper is K and can be evaluated as

In [17]: K._imp_(4.0)

Out[17]: array(16.0)

[Matthew Rocklin]

Here is a printout of the error message: 

MissingInputError: ('An input of the graph, used to compute TheanoInterpWrapOp.theanointerp(y), was not provided and not given a value', y)

What this says is that some nodes in your graph (TheanoInterpWrapOp.theanointerp(y),) weren't given access to all of the inputs that they needed.  This would happen in Theano if, for example, 

x = theano.tensor.vector('x')

y = theano.tensor.vector('y')

z = x + y

f = theano.function([x], [z])  

Notice that we're only giving function x and asking it to compute z.  This code will produce a similar error to what you're receiving.

[Guy Parsey]

That is precisely what I don't understand. In your example we are neglecting to give the function all of its inputs and the error message:

MissingInputError: ('An input of the graph, used to compute Elemwise{add,no_inplace}(x, y), was not provided and not given a value', y) 

is saying that we have forgotten the y input. In the test case I am doing, which yields the message:

MissingInputError: ('An input of the graph, used to compute TheanoInterpWrapOp.theanointerp(y), was not provided and not given a value', y)

Though the message is practically identical, when running the test case with my verbose print statements, the following is printed before being passed to the Op pulled from the TheanoPrinter.cache

children:  [y]

child types:  [<class 'theano.tensor.basic.TensorVariable'>]

followed by the return line:

self.cache[newkey](*children)

where self.cache[newkey] is the theano op. Doesn't this mean that the theano variable y is being passed to the theano Op or does this y not carry a value? Instead of passing the theano variable y to the op, should I be passing it to a theano.function of the Op?

Cheers,

Guy

[Matthew Rocklin]

Looking at your code it's difficult for me to see what you're doing with the cache.  It looks like you're trying to fill it with a particular value so that SymPy's theano_function call latches onto something you've already built.  

Instead, I recommend that you use theano_code, to transform individual sympy expressions into Theano variables like so

theano_inps = [theano_code(inp) for inp in inplist]

theano_outs = [theano_code(expr) for expr in exprlist]

Then, if you want to build more Theano expressions with your custom Theano op you can do so in standard Theano using these standard Theano variables.  I think that this will be simpler than reverse engineering a caching system in order to inject a pre-built theano variable into the graph.

In short, only use theano_function if you're only using pure SymPy.  If you're doing more complex things then translate your SymPy expressions to Theano expressions and work in pure Theano from then on.

[Frédéric Bastien]

Hi,

I was in vacation. Have you been able to make this work?

Fred

[Guy Parsey]

Hello,

Due to an up coming conference and a need to get my code in a working state, I have taken a brute force approach which seems to be working pretty well (I am in the middle of testing it). I am very much interested in working with Matt's advice on building up the ODE function from Theano variables, but due to time constraints and an unfamiliarity with Theano, my short-term solution is using Sympy theano_function while removing complicated parts (spline interpolations wrapped as function_symbols+argument_symbols) and replacing them as other symbols. A python wrapper takes in an argument list (in the format for scipy.integrate.odeint for now), a standard lambda function uses mapped arguments to evaluate each function+argument pair only once per call, the spline results are then passed to the theano_function function (this also works for C codegen) as extra arguments. This method has worked with all of my test cases so far (including a few I was having trouble with before), but I have some larger simulations soon to come. 

I use this method on my system of ODEs along with the jacobian (sympy haddles the creation of this really well). Though sympy.lambdify so far seems to be quite a bit faster to generate a callable function, this wrapper method has been about 5x faster with function calls, 25x faster with jacobian calls (A good number of the jacobian elements are zero, could this be what is causing this difference in timing?) and does not encounter the limits that initially drove me here. It definitely feels like a brute force solution, but broken code doesn't return results to present :) 

With regards to my test example I posted, I feel as though I am missing something rather simple but it also makes complete sense that creating the system of equations from theano variables would be a lot more efficient if Theano evaluation is the final goal. Once the time crunch is not as aggressive I will return to trying to find a cleaner solution.

Thank you all very much for your help and advice. I have learnt quite lot with regards to both SymPy and Theano.

Cheers,

Guy

[Frédéric Bastien]

Hi,

When you tell 5x and 25x faster, what do you time?

Fred

[Guy Parsey]

Hey Fred,
It is probably an over simplified version of timing. I am using the iPython timeit routine while calling the numerically implemented functions (either the sympy.lambdify created function or my wrapped function) with some example numerical input. I am planning on using the created function in an integration process, so the calling time is more important than the set up time.
Cheers,
Guy