[SciPy-user] PyDSTool and SciPy (integrators)

[Rob Clewley]

Dear colleagues,

Picking up on our previous email:

One of the features of PyDSTool that seems most broadly useful to the
scipy community is the facility for fast integration of ODEs using C
and Fortran based codes from E. Hairer (dopri853 and radau5). Our
implementation dynamically creates C-based right hand sides which are
compiled and linked against the dopri/radau .o files to create a DLL
file, which is then loaded as a python module. Integration of vector
fields is very fast, as everything is handled at the C-level, rather
than through callbacks to python functions as in scipy.odeint.  (This
is a *very* important distinction regarding performance, vs. merely
having a C-based integrator but retaining the user's vector field
function at the Python level). Though we have not done rigorous
testing, we have seen speed-ups of 10-50+ times over scipy.odeint.

In addition, we are able to generate 'event' functions and 'auxiliary'
functions that are calculated alongside and 'online' with the
integration of the ODEs. These (nearly) arbitrary functions of the
phase space variables and parameters may be used to detect, e.g. when
the computed trajectory crosses a boundary, achieves a local maximum,
or when some function of the phase space variables and parameters
reaches a zero-crossing. Event functions may be used to halt
integration when a particular condition is achieved (critical for use
in hybrid systems, for example), while auxiliary functions may be used
to calculate quantities of interest that are complementary to the ODE
calculations (e.g. ionic currents across membranes in neuronal model
ODEs that describe changes in potential as ion channels open and
close).

We are also able to use autonomous external inputs as part of the ODE
vector fields. For instance, we can use real data from experimentally
obtained voltage traces from a neuron as input to an ODE model of
neuronal membrane voltage dynamics.

One catch to our system is in the generation of C-based right hand
sides, which must be compiled. We use distutils to accomplish this in
a simple call that is naturally platform independent. This is a great
benefit to us and saves us relying on a third party library or
application, but this solution has some frustrating quirks (we very
much understand that distutils was not *designed* to be used this
way!). Also, we have extensive machinery to generate appropriate C
vector field descriptions from Python vector field descriptions, and
our event/auxiliary function mechanisms currently depend on
this. Model description, compilation, and integration are rather
intertwined at present.

As it now stands, a model (ODE right-hand side, event/auxiliary
functions, inputs) is defined as a set of strings, which are converted
to a representation as Variable, Parameter, and Expression classes (or
the string step is skipped and the model is just given in terms of
classes). The class representation is checked for consistency and is
used to generate a text file containing C code, including event
functions, etc., which is then compiled and linked against the
integrator and event detection code. SWIG is used to generate an
interface to python; the resulting module is then loaded into the
(interactive) python session.

When the vector field is integrated, the results (trajectories,
auxiliary function values, event values) are returned numpy arrays,
packaged as Trajectory classes (which have added functionality over
numpy arrays).

If our integration facilities were to be incorporated into scipy, the
underlying code for our integrators could be adapted to used python
callbacks, like odeint, though this would negate much of the speed
advantage. If the C compilation route were maintained, then we would
need to figure out a better way to use distutils, or find an
alternative compilation mechanism. In either case, we would have to
determine if and how to keep the event/auxiliary function and
autonomous external input facilities, which are closely bound to the
model building routines. The results of integration could be returned
as plain numpy arrays, but even doing this requires quite some
disentanglement.

Since integration is a very common task, and our integrators perform
well, this appears to us to be the best place to start moving PyDSTool
capabilities to scipy, either directly or in scikit form. However,
such a move would require significant effort on our part. We welcome
comments/advice/suggestions regarding community interest in and the
planning/design of such an effort. In particular, the issues of model
building, C compilation vs. python callbacks, and the use of distutils
are points where we are looking for input. We have already begun
seeking help from a consultant on fixing some of the quirks of using
distutils this way, in case we can continue to viably use it. But, for
instance, we'd like to get an update from Scipy developers about the
intended/expected future of distutils in scipy vs. the regular python
version (we use scipy's for Fortran compatilibity). Maybe it would be
reasonable to consider a modification of distutils (at least in its
API) so that it could be viewed as a platform-indepedent way to
compile DLLs on the fly in the way we need. Maybe you can see a much
better solution to solve our problem, given the ever-changing python
technology available (Spyke looks like it might be relevant, for
instance).


Regards,

Erik Sherwood
Rob Clewley

-- 
Robert H. Clewley, Ph. D.
Assistant Professor
Department of Mathematics and Statistics
Georgia State University
720 COE, 30 Pryor St
Atlanta, GA 30303, USA

tel: 404-413-6420 fax: 404-651-2246
http://www.mathstat.gsu.edu/~matrhc
http://brainsbehavior.gsu.edu/