[SciPy-user] PyDSTool and SciPy (integrators)

Thu Apr 17 09:17:15 EDT 2008

If it is of any use to get started, I have wrapped dopri5, and dopri853 using
cython, it allows for both c and python callbacks, and I am actively trying to
get your event detection setup working. I would be happy to talk with the
PyDSTool people if you guys are interested in a Cython approach.

Sadly as radau5 is using fortran my method is not yet ready for this system,
though I am hoping that I can learn enough f2py to make this a happen.

My wrapper also supports both range initial conditions and odeint style
specific step points (this is similar to how ode45 works, that is if the user
gives to time points, it returns the ode solvers chosen irregular steps,
otherwise it returns steps at the given intervals).

Gabriel

On Thu, Apr 17, 2008 at 01:24:39AM -0400, Rob Clewley wrote:
> 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/
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user