Novice scipy user: having some efficiency issues with using scipy.integrate.odeint for my research project, and not sure how to proceed...
Just to give you a bit regarding my background: I am a math undergraduate working on biology modelling problems. I have a strong background in Python, but it is only over the last few months that I have begun to dabble in scientific computation using scipy, numpy and most recently, Cython. This semester, I am taking my first proper course on numerical methods. Currently, in my research work, I am writing a program that solves a large number of non-linear ODEs. I use `scipy.integrate.odeint` as my workhorse, although I have dabbled with using `pyDSTool` as well, but found that the limited capacity I had in understanding advanced methods didn't allow me to use `pyDSTool` so that I had better performance than vanilla `scipy.integrate.odeint`. My supervisor recently pointed out a fairly major error in my previous implementation of the program: I was "breaking up" `scipy.integrate.odeint`'s flow. So if I had to solve a problem between times 0 and 10, I would make `odeint` work between `t in [0, 1]`, take the result at `t=1`, use it to update (non-linear) parameters for the next interval (`[1, 2]`) and then finally run the next interval using `t=1` results as the initial conditions. I thought it would be a good idea to do this in order to temporarily reduce the complexity of the problem by having some of the coefficients be constant over each sub-interval (something like a "quasi-equilibrium" condition?), and it also allowed me to straightforwardly implement non-DE rules to update parameters, if I needed to. However, problems with doing this include introducing errors, and causing loss of speed due to the breaking/restarting of `odeint`'s internals. So, I changed the implementation so that everything was done in "one-shot" (luckily the only non-DE rule I had so far was easily converted to a DE rule, although there is one coming up which I don't see any easy conversion for...[that is the subject of this question]( http://scicomp.stackexchange.com/questions/14765/scipy-integrate-odeint-how-...) ). However, I am finding that the "one-shot" method is *slower* than the "broken up" method! For a single cell (the number of ODEs increases linearly with increasing number of cells), solving for the interval between `[0, 1]` takes about 3 times longer than the "broken up method" (where I was breaking up the problem into intervals `[0, 0.5]` and `[0.5, 1]` for the same `rtol` and `atol`. Solving for the behaviour of a single cell over `[0, 10]` takes roughly 3 times longer than `[0, 1]` (hey! not bad!). Solving between `[0, 1000]` takes roughly 1500 times longer (weird :( ). Solving 2 cells (and solving between `[0,1]`) takes about 3 times longer than the single cell case, and solving 10 cells over the same time interval takes 60 times longer! So the times aren't growing linearly...whether I am increasing the length of the time interval, or increasing the number of cells. Alright, part of the problem can be explained by the fact that even though the number of ODEs grows linearly with the number of cells, the work done to calculate the Jacobian doesn't grow quite grow linearly -- still, in all cases, the "one shot" method is significantly slower than the "broken up" method (and the greater than linear effort required in calculating the Jacobian also affected the "broken up" method). Also, the Jacobian doesn't become harder to calculate when increasing the time interval...so shouldn't the time taken to solve the ODE grow linearly there at least? Alright, that's about all the certain information I have. My guess as to the number one culprit is that in the "one shot" method, the function used to calculate the Jacobian is doing a *lot* more work: recall, non-linear ODEs, so some internal iteration is required, and basically every iteration, all the coefficients are calculated again, while in the "broken up" method, I could choose which ones were held constant for a small time period, before recalculation. I thought reimplementing by making significant use of Cython (I used "annotation profiling" (`cython -a my_file.pyx`) to make sure I am cutting down on C lines wherever possible...so I am pretty sure I have not done a horribly ugly Python to C conversion. Still, even Cythonizing doesn't help with the fact that a lot more work is being done. I'm at a bit of a loss right now as to how to proceed - what questions should I begin asking myself in order to get a handle on things?
Hi Brian, I know we've talked about this before, but here are some more general suggestions about your problem.
So, I changed the implementation so that everything was done in "one-shot" (luckily the only non-DE rule I had so far was easily converted to a DE rule, although there is one coming up which I don't see any easy conversion for...[that is the subject of this question](http://scicomp.stackexchange.com/questions/14765/scipy-integrate-odeint-how-...)).
Having looked at that, your problem is not so bad. You will need to solve a piecewise (a.k.a. "hybrid") system, but you shouldn't be discretizing it based on time but based on events -- i.e. when your particle crosses into a new discrete area you are calling p. You then have a discrete state update for your parameters and then you can restart your smooth integrator until the next transition. This is a well-studied problem in this form and is numerically soluble with a hybrid solver such as that provided by PyDSTool. Between positional transitions, your dynamical system is smooth and can be solved with a regular ODE solver. The Jacobian within each domain should be simple enough to pre-calculate (as a function of p or <p>). There are a few hybrid model examples with PyDSTool on the tutorial and I am willing to help a bit with the setup once you've given a new script for your problem a shot. Take a copy of a helpful example (e.g. the SLIP pogo stick dynamics) and adapt it to set up what you can. Put in some comments etc. of what you need to happen. PyDSTool will be able to convert the ODEs to C code automatically but not the transition rules, which will still happen in python. This will not be the *fastest* way to solve it but, more importantly, this way will give you an accurate solution in a form that you can understand and manipulate, and you can worry about optimizing speed later, IMO. -Rob
Hi Rob, Thanks for your response! I didn't realize I had one here, until today. Since I wrote that question, I have been done two important things: 1) learn more math: last semester was numerical methods, this semester is dynamical systems 2) learned how to use Cython -- already, this has helped me speed up my simulations to the point where making things like bifurcation diagrams is a bigger deal for me (another thing I know PyDSTool can do too) This semester, in my dynamical systems course, we are going to be using lots of XPPAUT, but I plan on using PyDSTool. I have already had a look at tutorials like: http://www.ni.gsu.edu/~rclewley/PyDSTool/Tutorial/Tutorial_VdP.html I am excited! Again, sorry for not seeing this message earlier, and thanks for offering to help! I will post on the discussion board if I find something specific I could use some help with. Kind regards, Brian On Wed, Oct 8, 2014 at 8:26 AM, Rob Clewley <rob.clewley@gmail.com> wrote:
Hi Brian,
I know we've talked about this before, but here are some more general suggestions about your problem.
So, I changed the implementation so that everything was done in
"one-shot"
(luckily the only non-DE rule I had so far was easily converted to a DE rule, although there is one coming up which I don't see any easy conversion for...[that is the subject of this question]( http://scicomp.stackexchange.com/questions/14765/scipy-integrate-odeint-how-...) ).
Having looked at that, your problem is not so bad. You will need to solve a piecewise (a.k.a. "hybrid") system, but you shouldn't be discretizing it based on time but based on events -- i.e. when your particle crosses into a new discrete area you are calling p. You then have a discrete state update for your parameters and then you can restart your smooth integrator until the next transition. This is a well-studied problem in this form and is numerically soluble with a hybrid solver such as that provided by PyDSTool. Between positional transitions, your dynamical system is smooth and can be solved with a regular ODE solver. The Jacobian within each domain should be simple enough to pre-calculate (as a function of p or <p>). There are a few hybrid model examples with PyDSTool on the tutorial and I am willing to help a bit with the setup once you've given a new script for your problem a shot. Take a copy of a helpful example (e.g. the SLIP pogo stick dynamics) and adapt it to set up what you can. Put in some comments etc. of what you need to happen.
PyDSTool will be able to convert the ODEs to C code automatically but not the transition rules, which will still happen in python. This will not be the *fastest* way to solve it but, more importantly, this way will give you an accurate solution in a form that you can understand and manipulate, and you can worry about optimizing speed later, IMO.
-Rob _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user
participants (2)
-
Brian Merchant -
Rob Clewley