Generators

Overview
Generator basics
1. Exposing the RHS specification, Jacobian, and user-defined auxiliary functions
The ODE / DAE integrators
Difference equations and maps
Time-series trajectories and lookup tables
Trajectories or array data as external inputs to Generators

1. Overview

The Generator classes represent the basic simulation objects of PyDSTool, implementing ODE/DAE integration, and the computation of discrete maps. More generally, they are also sources of curves, which may or may not be interpreted as trajectories. For instance, explicit or implicit functions (including those that are parameterized) can be used to specify Generators.

In programming terms, Generators are class instances that produce Trajectory objects (also viewed as curves, especially if they are not parameterized) on demand (see Trajectories). They take many forms, and may or may not use user-specified parameters and initial condition specifications in order to determine the exact Trajectory object that they produce.

In the very clear dynamical systems and application-oriented setting of PyDSTool, we believe there is no inappropriate clash of terminology in naming these classes Generators. The core of the Python language also contains a special class known as a generator, and we trust that no-one finds this confusing.

2. Generator basics

Generators compute multi-dimensional trajectories for the primary system state variables. However, additional auxiliary variables may be specified that can be explicit functions of the system state. The values of auxiliary variables are also part of a computed Trajectory object.

Generators come in two principal sub-classes, depending on the representation of the independent variable. This variable may either be represented in a continuous domain or a discrete domain. All Generator instances are sub-typed from one of these super classes, ctsGen or discGen.

Generators can also contain a specification for events that can be detected during the computation of a trajectory. These events may be specified as zero crossings of functions in the state variables or of the independent variable.

Depending on the type of Generator, the user may need to specify parameter values and initial conditions that uniquely determine the Trajectory that will be created when the Generator's compute method is called.

Some forms of Generator also permit autonomous external inputs to be supplied in the form of a Trajectory object.

After computation of a Trajectory, any errors or warnings produced during the computation are recorded in the Generator object, and can be queried or reset.

The "right-hand side" function specifications for Generators are not provided by hand-written code, which would be inconvenient and would burden the user with too much book-keeping with array indices, etc. Instead, the user specifies a dictionary of definitions for these functions, keyed by the variable names. In turn, these dictionaries can be created automatically using the ModelSpec and ModelConstructor class utilities for PyDSTool symbolic expressions, described elsewhere. The ModelSpec Quantity classes can also be used in place of string definitions (see the example tests/vode_withJac_Symbolic_test.py).

Preset sub-domains of the real numbers can be specified for each state variable, for the independent variable of the system ("time"), and for the parameters. These bounds can help with error checking and validation of a model, and can be crucial for setting up some types of hybrid system. The Generator initialization keys for these bounds are 'xdomain', 'tdomain', and 'pdomain', respectively.

The arguments to the Generator's initialization call can be made most easily by first assembling them in an argument structure. See the GettingStarted page for an example, under sub-heading "A simple ODE".

2.1. Exposing the RHS specification, Jacobian, and user-defined auxiliary functions

The RHS vector field function, any Jacobians, and all user-defined auxiliary functions are exposed to the user of the generator as python functions, and are called using the same calling sequence as they were defined. For a generator g, these are respectively defined as the methods Rhs, Jacobian, JacobianP (for Jac with respect to parameters), and auxfns.<aux_func_name>. Internal references to parameter names are mapped dynamically to the current values of the parameters in the generator's pars dictionary. If external inputs are defined for the generator then the auxiliary function wrappers have to be contained with the final __parsinps__ dictionary argument supplied. This dictionary must contain the values of the parameters and external inputs to evaluate the function.

The various forms of Generator and their usage are now described in detail.

3. The ODE / DAE integrators

3.1. Integrators using Python code vector fields

3.1.1. VODE integrator

The VODE integrator is used by the Generator.Vode_ODEsystem class. This class is a wrapper around the SciPy wrapper of the CVODE solver, and the specification of vector fields is done using Python functions. These are defined dynamically from the user specifications at the Generator's initialization, and are called from the C code to evaluate the RHS at each time step (hence they are known as call-back functions).

The SciPy wrapper is an elementary interface to CVODE. There is no in-built event detection to CVODE, which we therefore must provide at the Python level. This is not an efficient situation by any means, but it is a simple and somewhat reliable method for both stiff and non-stiff systems with events, and there is a small setup time. This Generator class also does not require 'external' configuration and installation -- it only relies on SciPy. For instance, there is no need for an external compiler that would include installation of SWIG, gcc, etc.. On a fast computer, non-stiff and low-dimensional vector fields integrate surprisingly quickly with this Generator, given the overheads involved!

Note that the only control over the time mesh that will form the basis of the computed Trajectory object is the parameter init_step. This makes this Generator's output resemble that of a fixed time-step integrator, but remember that internally CVODE is an adaptive time-step solver that is quite efficient for stiff systems. The limitation of the SciPy interface to CVODE is such that at the Python level only state values at every init_step time points are available for later use in a Trajectory. The max_step and min_step parameters will still be used internally by CVODE to control its adaptive step-size between the output steps that are of size init_step. The initial internal step size of CVODE will also be set at init_step.

The VODE Generator class can handle terminal and non-terminal events, which are known as "high-level" events (a sub-class of the Event class) because they are specified directly as Python functions. The event checking for this integrator takes place after the calculation of each time step at the Python level. The Generator adapts the time step used by the C code when events are found in order to resolve their position accurately (if this option is selected).

Important: Unfortunately, non-terminal event points are not resolved to the same accuracy as for the C-based integrators. They are only linearly interpolated to the required 'precision' between data points separated by the chosen time-step. The resulting event points are still added to the underlying time mesh. This limitation does not affect the treatment of terminal events. Thus we suggest running the VODE integrator at a sufficiently small time-step to ensure non-terminal event accuracy.

The Vode_ODEsystem Generator does not support backward integration. Continued integration is available by the 'c' option to the compute call.

3.1.2. Making changes to the Python code through the API

Arbitrary Python code can also be added to the very beginning and end of the vector field specification function, using the Generator initialization option keys vfcodeinsert_start and vfcodeinsert_end. Note that these string literals must contain appropriate indentation ("top level" in the function is a four space indentation). Care must be taken to make any temporary variable declarations have names that will not clash with other, automatically-generated, local names. For the most part this will just be a case of not using the state, parameter, or input names. No warnings will be given if you mess this up!

The pseudo-parser that converts the specifications into actual Python functions makes various changes associated with proper library naming for math calls, and renames certain internal function calls to avoid name clashes that would occur with multiple instances of the same Generator. However, the parser does not currently process any code inserted using the above options. Thus, references to auxiliary functions will fail because their unique name has not yet been created. Also, math references in the inserted code should follow the usual Python naming scheme, i.e. sin(t) should be entered as math.sin(t).

3.2. Integrators using C code vector fields

PyDSTool supports ODE integrators that compile C code specifications of the vector field, in such a way that it can dynamically load a library containing the integration solver after compilation. Presently Hairer and Wanner's Dopri and Radau methods are implemented (see below), but a Taylor series method (using ADOLC) will appear soon.

On Win32 platforms, MinGW's gcc is assumed to be the default compiler. If you use Cygwin's instead, you must add the keyword option compiler="cygwin" to all initializations of C-based Generators.

3.2.1. Compiler messages

In old versions of PyDSTool the user may observe that the compiler displays a series of messages about the compilation. This is because compilation is done using utilities imported from the distutils package. In future versions of PyDSTool we hope to be able to control the unwanted output so that only error messages will appear (this is partly a problem with the distutils package). Compiler messages may be invisible on Win32 platforms using IDLE, which does not echo from stderr. This means that important compiler error messages may not be seen unless python is run from a different shell system (e.g. Windows command prompt or IPython).

Note: You may get an error message saying that SWIG failed during compilation of your vector field if your PyDSTool directory structure contains directory names that include spaces.

3.2.2. When is the DLL imported?

Once compilation is done, the resulting Dynamic Link Library is not imported into PyDSTool until the compute method is invoked. At that point, the integrator becomes available to the class.

3.2.3. Memory vs. speed efficiency

The implementation of the interface to the C code requires that the user sets a maximum on the number of points that can be computed by the integrator at any one time, before it is reset. The algorithmic parameter max_pts defaults to 10000, but may be set to up to over 100 times this amount, if long runs for large systems are necessary. The reason that this parameter should be set as small as possible is to minimize dynamic memory allocation, which in turn will keep the memory allocation time overheads to a minimum (they are substantial for large allocations).

If memory runs out during an integration, and you have not reached the "hard" upper limit imposed by the system, then you will get an error report that will include an estimate the new value of max_pts necessary to complete the current integration. After changing this value using one of the parameter-setting methods, you options for recovery depend on your situation. If you are running the integrator directly from a Generator, the Generator method compute can be used with the "continue integration" option 'c' in order to pick up from where the memory ran out. However, if memory ran out from an integration made via a Model object (for hybrid systems) then there is presently no alternative but to restart the entire calculation afresh.

For integration runs that are expected to be large and take a long time, it helps to have an estimate of the average step-size. From this the expected number of points for storage can be estimated before integration begins.

3.2.4. Making changes to the C code through the API

Arbitrary C code can be added to the very beginning and end of the vector field specification function, prior to DLL building, using the Generator initialization option keys vfcodeinsert_start and vfcodeinsert_end. Note that unless triple-quotes are used, these string literals must use \\n rather than \n in order to put the symbols '\n' into the C code, rather than an actual newline.

User-supplied C code is not checked for syntactic correctness until compile time.

By default, only the standard libraries stdio, stdlib, string, math are imported to the vector field file for use by code added through the API. For additional library access, use the include keyword option when calling makeLib or makeLibSource to specify a list of local .h files to include, or see the next section on adding code by hand. To use a library held either in the local project directory or in a standard library, simplify specify file's name as usual, i.e.

myDopriGen.makeLib(include=['mysourcelib.h','anotherlib.h'])

3.2.5. Making changes to the C code by hand

Alterations can be made to the C code specifying the vector field, etc., before compilation, if the nobuild option is set to True in the keywords of the Generator's initialization call. The partially-defined Generator object will be returned from the initialization, but the DLL will need to be built, either using the makeLib method or by the individual makeLibSource and compileLib methods.

User-supplied C code is not checked for syntactic correctness until compile time.

3.2.6. Making changes to the C code after initial trajectory computed

Once compute has been called, the present inability of Python to properly un-load DLL modules means that the user will not be able to make changes to the vector field specification, or the auxiliary variable or function definitions, without restarting the session and deleting the old .so (.pyd on Win32 platforms) library in the local dopri853_temp directory.

Interactive changes can be made to parameters specified through the API to the Generator object, however.

3.2.7. Dopri integrator

The Dopri integrator is used by the Generator.Dopri_ODEsystem class. It is based in C, as a wrapper around the original C code implementation by Hairer and Wanner of a Dormand-Prince algorithm. Events, auxiliary functions / variables, and so forth are handled in the C interface to this code. The Dopri Generator class sets up the integrator from the same specifications as the VODE integrator, but converts them into C routines. These routines contribute to a C file containing code for the vector field, auxiliary functions / variables, and so forth, which is compiled and linked to the integrator automatically (by default) by the Dopri_ODEsystem class at runtime.

The Dopri integrator does support backward integration -- simply call compute with the second argument 'b'. Continuing a previous integration is also supported using the 'c' option.

3.2.8. Radau integrator

The Radau stiff integrator is available in the class Generator.Radau_ODEsystem. The integrator uses embedded C code in much the same way as Dopri: only a few algorithmic options are different. See the /tests directory for examples.

At present, the automatic code generation may leave a messy trail of temporary directories starting in the source code directory, and will almost certainly produce an uncontrolled stream of log messages on stdio. If the building of the extension DLL fails please let me know the error and I'll try to help patch it up.

The control of target file placement in the directory structure through options to numpy_distutils calls appears to contain bugs. You may find that on your platform the compilation process for Radau results in an error whereby PyDSTool cannot find a compiled library target. In this case you can try to locate the .pyd, .so or .py file yourself deep inside the temporary directory structure, and copy it back to the directory containing your original PyDSTool model script. You should then be able to continue executing your script. We aim to change the manner by which we generate DLLs in version 0.90.

3.2.8.1. Mass-matrix systems and DAEs

Differential-algebraic systems are supported by Radau, when the mass-matrix for the system is made singular. The mass matrix is assumed to be the identity unless provided explicitly when a Radau generator is created. Two examples are provided in the /tests/ directory: DAE_example.py and freefinger_noforce_radau.py.

An important extension to Radau was made by Erik Sherwood in allowing the underlying Fortran code to support non-constant mass matrices, without resorting to an entirely separate (and much more sophisticated) integrator.

3.2.9. Preset integration time-steps

The Dopri and Radau integrators have been interfaced to PyDSTool in a way that supports the specification of explicit time points at which integration is forced to occur. This feature can be of use when repeating trajectory computations for different parameters using an adaptive time-step solver, but when comparison of numerically-computed trajectories is only truly fair when the time meshes for the trajectories are identical. (For instance, this feature is useful to systems biology simulations such as SloppyCell.)

Selection of preset times is done by specifying the specialtimes array and the use_special Boolean switch when creating a Dopri or Radau Generator. The VODE Generator currently supports this feature only in the version from the subversion repository (June 2007 -- see CodeTopics). Note that if event detection is set to be 'precise' then event times will appear in the time mesh of the output trajectory regardless of the setting of specialtimes.

3.2.10. Backwards integration

The Dopri and Radau integrators both support backwards integration. When calling compute, an optional argument 'b' can be provided. Erik Sherwood has tweaked the initial step size selection algorithm for backwards integration, to prevent it from choosing initial steps larger than the user-specified maximum time step.

3.3. Taylor-series integration

This is not yet implemented as we have not been able to get an automatic differentiation package to work with SciPy data structures. We hope to resolve this soon. In the meantime, Erik Sherwood has implemented a model export utility from PyDSTool to the Matlab-based specification necessary to use ADMC++ as part of the ODETools package that he helped to develop.

4. Difference equations and maps

These are implemented in pure Python in the MapSystem sub-class of Generator.

5. Time-series trajectories and lookup tables

Discrete data points from time-series can be treated as "discretized curves", or as discretely sampled trajectories from a dynamical system source. The InterpolateTable class provides simple linear interpolation over an ordered set of datapoints in order to be able to treat a time-series as a Trajectory in the context of other PyDSTool algorithms. Non-interpolated trajectories are supported using the LookupTable class.

6. Trajectories or array data as external inputs to Generators

Both the Python-based and C-based Generators support the inclusion of 'external inputs'. This means that a previously-computed trajectory (including from time-series data generated by a physical experiment) can be used as an autonomous input to an ODE/DAE/difference equation right-hand side. This is a powerful tool for data-driven modeling -- see the page DataDrivenModels.

Wiki

Page

User