[SciPy-User] Novice scipy user: having some efficiency issues with using scipy.integrate.odeint for my research project, and not sure how to proceed...

[Brian Merchant]

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-can-odeint-access-a-parameter-set-that-is-evolving-i)
).

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?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20141005/5e156208/attachment.html>