[SciPy-user] PyDSTool and SciPy (integrators)

[Rob Clewley]

On Thu, Apr 17, 2008 at 8:58 PM, Anne Archibald
<peridot.faceted at gmail.com> wrote:
>  This sounds like it might have been a problem for me: what was
>  happening was I was integrating an ODE within a domain. Outside the
>  domain my RHS function raised exceptions, so it was not even safe to
>  evaluate it there. I tried hacking up a solution with scipy's
>  one-step-at-a-time integrator modes, but the best I could do was stop
>  "near" the boundary, which failed if the integrator took bigger steps
>  than my notion of "near".

I don't think our event detection routine solves this, but IIRC the
one implemented for the sloppycell package by some colleagues of ours
at Cornell would solve your problem. It's based around the lsoda
integrator I think, but I think the API etc. is possibly even more
entangled with the specifics of their systems biology package than
ours is.

>  I think the right return format depends what people want from the
>  solution. I see two major kinds of problem one might use an ODE solver
>  for:
>
>  * Integrate out to some stopping point and give the values of
>  dependent and independent variable there.
>  * Integrate out to some stopping point and give back a representation
>  of the solution.
>
>  (The first case can in principle be obtained by evaluating the second
>  at the endpoint, but at the cost of replacing O(1) space by O(n)
>  space.)

I personally do not find the former very useful, as I almost always
want to know something about the intermediate points. My understanding
is the dopri and radau integrators inherently return the whole
solution, so maybe an additional calling option supporting the former
usage could just throw everything out but the last point before the
whole solution gets copied into a numpy object.

>  For the second it is clear that some kind of object to represent the
>  solution is a good idea. What exactly it should contain is open for
>  discussion; I envision keeping estimated function values and
>  derivatives at each step, along with convergence information, and
>  providing an interface to view it as either an array of values or a
>  function akin to our interpolating spline objects. Convergence
>  information and auxiliary function values may also be useful. Even for
>  the first, simple case, I think returning an object is a more
>  manageable way to provide convergence information.

This is definitely something that the community should discuss. IMO,
the more that can be stored for re-use, error analysis, etc. after the
fact, the better. Memory is (relatively) cheap for these problems.
Simple (even default) options can always be provided to turn them off.
Your only problem is that all of the integrators we're talking about
are *old* inasmuch as their APIs were not designed with this kind of
en masse export of information out to the user. Any of them (current
scipy integrators and PyDSTool's) would need their interfaces
substantially re-written to accommodate this.

>  Are there other styles of using an ODE solver that people want?
>  Perhaps a "running" mode, where the solver keeps track of whatever
>  state information it has and is asked to advance by a certain amount
>  (say, because more experimental data has become available)?

There are two issues that need separating here:

1) the short term refactoring, reuse and adaptation of existing
PyDSTool code for scipy
2) the long term design goals for supporting ODEs in scipy

Much as I want to keep (2) in mind, I think getting some short-term
progress on improving scipy's ODE support based on existing code would
be good. I foresee it taking a long time to get ODE support truly
"right" for as many users as possible, and in the meantime I think
some our code and other short-term solutions would be an appropriate
stop gap.

-Rob