UserDocumentation

  1. PyDSTool basics
    1. Modularity and embedding
    2. Design philosophies
    3. Numerical classes
    4. Symbolic classes
    5. Other classes
  2. Starting your own modeling project
    1. Filespace preparation
    2. Mathematical preparation
  3. Basic I/O and object manipulation
    1. Loading and saving data and objects
    2. Deletion
    3. Sessions
    4. Memory management
    5. Copying of PyDSTool data structures
  4. Visualization
  5. System specification
    1. Two routes to specification
    2. Elements of a system specification
    3. Mathematical expression syntax
      1. Syntax checking
      2. Standard math names
      3. Multiple definitions using the for loop macro
      4. Auxiliary functions
      5. Noteworthy syntax quirk: powers
      6. Sub-expression substitutions in target language
  6. Solving for trajectories
    1. Generators
    2. Trajectories
  7. "Models" and hybrid dynamical systems
    1. Hybrid trajectories
    2. Important model methods
    3. Absolute and relative time
    4. Setting up a Model object
      1. Model objects as non-hybrid systems
      2. Model objects as hybrid systems
  8. Bounds and constraints
  9. Continuation and bifurcation analysis
  10. Porting models to other packages

1. PyDSTool basics

Python is a high level, interpreted, language with powerful object orientation. PyDSTool is primarily a set of Python classes for use within a regular, interactive, Python session. These classes themselves rely on several SciPy and numpy classes, particularly the array class.

The user can build objects from the PyDSTool classes interactively, and can manipulate them using PyDSTool utility functions, or using scripts or modules that the user has written. Visualization tools are provided via SciPy and Pylab in order to study the results of computations, in much the same way as is done in Matlab. Being built to work with common packages such as SciPy, the user benefits from being able to apply well-known algorithms from other packages to PyDSTool objects.

For a full introduction please see the page ProjectOverview.

1.1. Modularity and embedding

As well as providing simulation and analysis routines of its own, PyDSTool is also intended to provide an environment for embedding numerical calculations within eachother, in a hierarchical fashion (a.k.a. "nesting"). A design goal was to make the embedding as straightforward and as general-purpose as possible, with the least amount of effort required by the user. PyDSTool provides core Python classes from which a user can take existing numerical computation packages and apply them to high level objects. This is often done by "wrapping" the external packages, providing them with a common API for use within Python. The packages may take the form of simple executables, or as shared libraries that can be loaded into Python after they are wrapped (e.g. using SWIG). Embedding is also made more transparent through Python's dynamic typing system, whereby different classes can be treated in common ways by providing common interface methods.

1.2. Design philosophies

The core elements of PyDSTool are the classes that represent or relate to trajectories of dynamical systems, which are discussed next.

1.3. Numerical classes

From the bottom up:

Note that in discrete spaces, a "curve" refers to an ordered set of discrete points (a "discretized curve").

1.4. Symbolic classes

From the bottom up:

1.5. Other classes

Beyond the numerical classes are those that create them as part of solving systems of dynamical equations. Generator objects create Trajectories, for instance. A Model object allows one or more Generators to be combined into (possibly hybrid) dynamical systems models, also providing many additional user utilities for interacting with the model.

2. Starting your own modeling project

2.1. Filespace preparation

Create a directory somewhere using the name of your project. Make sure the path to the PyDSTool package is in Python's path. In a new Python script that will execute your computations, put from PyDSTool import * in the header. Any project-specific compilation of program components (e.g. using numerical integrators implemented in C) will require temporary files to be created. These will appear in automatically generated sub-directories of your project directory.

2.2. Mathematical preparation

It is recommended that, to the fullest possible extent, variables and parameters of the dynamical systems studied using PyDSTool are scaled so that they are "order 1" in magnitude. This helps the software accurately deal with occasional rounding errors in the resolution of intervals (see the section on Bounds and Constraints on this page).

3. Basic I/O and object manipulation

3.1. Loading and saving data and objects

The reading and writing of text-format data files is supported. This can be used to import or export trajectory data, for instance. The commands for these are importPointset and exportPointset. To export arrays more easily than using scipy.io, run exportPointset(arrayToPointset(a), {'fname.dat': []}) to save the array in text format to file 'fname.dat'.

It is also possible to load and save full PyDSTool objects, for re-using objects in different sessions. These functions are loadObjects and saveObjects. Multiple objects can be saved in a single file, and named objects (such as trajectories or generators) can be loaded individually from a file containing multiple objects.

N.B. For Windows users, saving and loading PyDSTool objects may be a little slow. This is because there is a bug in Python's "pickler" in Windows implementations that concerns IEEE 754 special values such as Inf and NaN. As a result, the workaround is a non-binary data format which is less compact.

3.2. Deletion

The regular Python command del can be used on any PyDSTool object in order to delete it from memory. The whole session can be cleared using the restart command (see below, under Memory Management).

3.3. Sessions

The current PyDSTool session can be saved and reloaded using the commands saveSession and loadSession. By default, saving the session will only save objects that belong to PyDSTool classes. Setting the argument deepSearch=True in a call to saveSession will cause the system to search all top-level lists, dictionaries and tuples for PyDSTool objects, and save those top-level lists, etc. if there is a match.

Non-PyDSTool objects can be saved separately using the saveObjects and loadObjects commands, or the scipy.io.save function can be called to save every variable in memory, to be re-imported.

3.4. Memory management

PyDSTool objects presently resident in memory can be identified using the who command. who can take an optional argument, either an individual type name or a list of type names, that filter the output to only include PyDSTool objects of those types.

A more specialized option is to add the returnlevel=N option, where N > 0. In the argument's absence, N defaults to 0. When N is 1, the call to who does not print anything to the screen. Instead, it returns a list of the objects found. This can be useful when you want to pass to a function all objects of a certain type that have been defined in the default namespace. N = 2 causes who to return a dictionary of object name keys mapping to the objects.

The default namespace for the search is the globals() dictionary of the calling scope. To specify an alternative namespace, add the option objdict=<scope_name_dictionary> to the call.

The optional argument deepSearch=True will search for PyDSTool objects contained up to one-level deep in lists, tuples and dictionaries defined in the supplied namespace objdict (defaults to the calling namespace's globals dictionary).

The current session can be restarted using the restart command. By default this merely clears the global registries of declared symbolic and ModelSpec names, but the additional argument delall=X where X=1 or 2 will delete all PyDSTool objects (X=2 causes a "deep search" to find them contained in lists, tuples, and dictionaries).

Important: At this time, restart cannot clear from memory any dynamically loaded C-based ODE integrators such as Dopri and Radau. This is a problem with the underlying behaviour of Python itself, and we hope to fix this in the future.

3.5. Copying of PyDSTool data structures

The larger data structures in PyDSTool generally contain several levels of objects, which themselves include arrays, IEEE754 Floating Point special values, and sometimes dynamically created method functions, etc. This includes classes such as Model, Generator, Trajectory, Variable, Pointset, and the symbolic and model specification classes QuantSpec, Quantity, ModelSpec, and so on.

The standard Python copy and deepcopy functions will both behave identically for all these major PyDSTool classes, and will always create full copies without references. This has been implemented for the sake of safety, as shallow copy-by-reference is rarely a desirable outcome for these large data structures.

Copying may however be a little slow for very large data structures, as the only typesafe way to deep copy complex classes is through a tweaked version of the standard Python "pickler". A web search will show that there are problems with using deepcopy on complex data types, and even the C pickler cannot be relied upon to deal with IEEE754 special values properly (at least not on Windows platforms). Hence, we use a patched version of the standard pickler, which on Windows platforms uses a non-binary data format and is therefore less efficient. This issue is discussed further, in relation to handling dynamically-created functions, on the TechDocumentation page.

4. Visualization

Pylab's MatPlotLib is the graphics library currently recommended the most for 2D graphics, although no core PyDSTool functionality relies on any particular graphics library. The command interface of MatPlotLib is intentionally similar to that of Matlab. For 3D plots, you are free to choose your own as there is no obvious standard yet. One option is "old" Scipy's own gplt interface to GnuPlot, but the output quality is lower than that of MatPlotLib. A 3D toolkit emerged for MatPlotLib, that is called [WWW] mplot3D, but it is unclear whether it is still in active development.

The most important thing to remember when plotting trajectories is that a sample of the trajectory must be extracted from the trajectory object. This is partly because trajectory objects hide their underlying time mesh, but also because the plotted mesh resolution often has to be chosen differently to the time mesh of the original data source anyway. In an update of version 0.83, trajectories that can be sampled with default options (by calling sample with no arguments) will be acceptable plot call arguments, without first needing to be converted to Pointsets or arrays.

MatPlotLib is automatically imported when PyDSTool is imported (or ignored if it is not present). The PyDSTool interface to MatPlotLib has been slightly altered in two ways, for convenience.

(1) The plot function is wrapped, and should be called without its original package prefix pylab.. The wrapped version permits singleton numeric types to be used in the first and second arguments, and more importantly will automatically sample trajectory objects where possible.

(2) The save_fig function provides an easy way to save a figure in multiple formats with one call.

5. System specification

5.1. Two routes to specification

There are two routes to the specification of a dynamical system.

In addition, you can import a vector field definition from [WWW] VFGEN using an XML format, or from SloppyCell using an included conversion tool.

5.2. Elements of a system specification

A system specification is made up of some or all of the following types of definition:

Auxiliary variables are meant for the output of useful quantities derived from state variables, parameters etc. Thus, a specification of a state variable may not refer to an auxiliary variable. Specification of valid domains for these variables, parameters, and the independent variable may be optional, depending on the target Generator. These are supplied as two-element lists containing the interval range endpoints that can include Inf. Also see BoundsSafety for uses of the domain information. If you want to avoid costly re-calculation of terms, either define auxiliary ("helper") functions or utilize the reuseterms option in your system definition. Examples of use are provided in the /tests/ directory.

User-defined events (see Events) and autonomous external inputs (see DataDrivenModels) can also be provided to the system specification.

5.3. Mathematical expression syntax

PyDSTool's API provides a simple syntax for the specification of expressions and functions that define Generator objects -- i.e. differential or difference equation right-hand sides, and auxiliary ('helper') or explicit functions. We have not undertaken the time-consuming effort of building a true parser to deal with this, and there are some restrictions on what can be specified through the API, and what syntax or naming errors can be detected.

The expression syntax is independent of the target language, and is "compiled" into target language specifications by the FuncSpec class (see FunctionalSpec). Symbolic manipulation of the abstract syntax prior to "compilation" is described on the page Symbolic. The specification of structured models, that exhibit hierarchical modularity and other such symmetries, is described on the page ModelSpec, along with examples.

5.3.1. Syntax checking

There is a small amount of syntax checking provided at 'compile' time, mainly to check that all string tokens reference known variables, parameters, inputs, auxiliary functions or math objects. No semantic checking is done, and errors not caught at the time of initializing a Generator may lead to cryptic Python errors during trajectory generation.

5.3.2. Standard math names

Standard math names such as sin, pow, and pi do not need prefixing with a library reference such as in math.sin (for a Python target). Time must be referenced as t, and references to variables, parameters, and external inputs is by the corresponding declared name. This is known as "index-free" notation. Internally, The FuncSpec.py module will translate these names to array index references once and for all. The user does not have to know or use these indices. As of version 0.85, some support for special math names (such as the error function provided by SciPy) is provided. This support will grow in future releases.

5.3.3. Multiple definitions using the for loop macro

A 'for' loop macro is provided to ease repetitive definitions in all specification strings. This is pre-processed to `unroll' the loop into a flat list of definitions before the parser processes the specification string. The syntax is:

for(i, <ilo>, <ihi>, <expr_in_i>)

where <ilo> and <ihi> are integers, and the expression in i (or any other alphabet character) has all occurrences of i in square brackets replaced with the appropriate integer. Parameter declarations for such expressions are strings of the form x[i] for a variable x.

5.3.4. Auxiliary functions

Auxiliary functions are accessible to variable specifications and event definitions by name. As well as the option for user-defined functions, several in-built auxiliary functions are provided. Currently, these are:

The first two are present for dealing with discrete-time systems (i.e., MapSystem), or legacy ODE code (e.g., from XPP) when the discrete events they involve are not computed beyond the accuracy of the integrator's step size. For accurate discrete event handling in continuous-time systems, the hybrid system Model class should be used.

The user may specify other, user-defined, auxiliary functions. These can help to simplify and organize the evaluation of right-hand sides, and to improve computational efficiency. User-defined functions are exposed in the generator as python functions (as of March 20th 2007 in the subversion repository). For instance, a function 'h' of one argument for a generator g is exposed as g.auxfns.h and has the same calling sequence as the internal function. Any parameter values used in the definition of the function are also mapped dynamically to the current values of g.pars.

5.3.5. Noteworthy syntax quirk: powers

For Python language targets you should always use the syntax pow(x,p) rather than x^p or x**p in order to express powers involving floating point values. Python does not behave well with the carat operator ^ for floating point numbers. In order to make it easier to convert legacy code to work with PyDSTool, the utility convertPowers(<string>, <target>) is available, where <target> is "^", "**", and "pow".

If you use the symbolic expression classes (see Symbolic), you can build expressions using the ** operator. Internally, the represenation for powers are converted to using the pow function. (Later, this perhaps should be changed to SciPy's power function.)

5.3.6. Sub-expression substitutions in target language

There is a helpful option to use a two-pass sub-expression parser and term reuser in expressions appearing in the target language vector field function definition, and in auxiliary functions / variable defintions. The user can list pairs of (subexpression, symbolname) strings in a dictionary at the initialization of the Generator, through the keyword reuseterms. If the subexpressions are actually found in the function code, this causes 

double symbolname = subexpression;

in C-code Generators, and 

symbolname = subexpression

in Python-code Generators. to appear in the preamble of the appropriate function, and all occurrences of subexpression will be replaced in the RHS. The two-pass aspect allows 'second-order subexpressions'. For instance, consider the specification

2*sin(afunc(p))+cos(afunc(p))*sin(afunc(p))-afunc(p)/5.

A reuseterms dictionary  {'afunc(p)': 'afp', 'sin(afp)': 'sa', 'other(x)': 'ox'}  will result in the final specification:

2*sa+cos(afp)*sa-afp/5.

preceded by the declarations: 

double afp = afunc(p);
double sa = sin(afp);

(in C) 

afp = afunc(p)
sa = sin(afp)

(in Python)

Note that, because other(x) did not appear in the code specification, no declarations for ox are made. This avoids some compiler warnings about unused symbols.

The two-pass system ensures that second-order references are not processed until the second pass, and that the declarations appear in the correct order.

6. Solving for trajectories

6.1. Generators

Trajectories of dynamical systems are determined by invoking the compute method of a Generator object. This may cause an ODE to be solved, for instance if the Generator is a sub-class of ODEsystem.

The abstract and target-language independent specification of dynamical systems is converted into target simulation objects that perform trajectory computations on demand. The target class is Generator, and has many sub-types corresponding to the different types of dynamical system supported.

The Generator classes are discussed in detail on the page Generators.

6.2. Trajectories

The trajectory class can represent continuous curves or discrete sets of ordered points. In the former case, it provides the illusion of dense output for trajectories that are defined only by an underlying discrete mesh. This is done by linear interpolation between mesh points (see Trajectories for details about this). Trajectories defined by explicit functions (including Taylor series, which will be implemented at a later date), are intrinsically dense.

In essence, trajectories are simply multiple Variable objects brought together in one class, under a common independent variable.

Calling a trajectory can be done in two ways. Either a single independent variable value can be passed as the first argument, or a list (or tuple or array) of such values, which will result in the return of a Point or Pointset object (respectively). Optionally, a list of coordinate names may be specified in the second argument, the order and content of which determine the fields of the Point of Pointset returned (and their order). Passing multiple values of the independent variable in the call ("vectorizing" the call) is much more efficient than repeatedly calling a Trajectory object at single values.

The implementation of Trajectory objects is described on the page TechDocumentation.

7. "Models" and hybrid dynamical systems

The PyDSTool Model class serves two purposes. The first is to generalize the forms of dynamical system that can be simulated in PyDSTool, by supporting hybrid dynamical systems (see the page HybridSystems). The second is to provide additional utilities for introspection and manipulation of a dynamical system, and organization of associated data, above what the Generator class already provides.

A model instance holds information that determines which Generator is used to calculate a trajectory segment, given the initial state. When a terminal event occurs during the calculation, the model uses a rule to determine which Generator to use next.

7.1. Hybrid trajectories

Hybrid trajectories are stored inside the Model object they are associated with. They can be accessed much the same way as regular trajectory objects when treated as continuous curves. However, hybrid trajectories can also be treated as discrete mappings, from one "terminal event" to the next. To access a hybrid trajectory this way, we call the Model instance that contains it with the additional option 'asmap=True': my_model(<trajname>, <integer>, asmap=True), where the integer represents the number of the partitions to access. These integers may range from 0 to the number of partitions minus one. The number of partitions can be found by calling my_model.numPartitions(<trajname>).

In a Model object, alongside the hybrid trajectory is information about the event parameters used during its calculation. In fact, there is an object containing various types of information about the trajectory, referenced as my_model.trajectories[<traj_name>]. The event parameters are stored in a copy of the EventStruct object associated with each Generator used to create the hybrid trajectory in a dictionary attribute 'genEvtStructs', accessed in the following way: my_model.trajectories[<traj_name>].genEvtStructs[<gen_name>].

7.2. Important model methods

When you have built a model, use help(<methodname>) to read the documentation string for any of these methods.

7.3. Absolute and relative time

When used as a complete model in its own right, all references to time in a Generator object are to the absolute ("global") time. However, when embedded in a Model object as a hybrid system, the concept of absolute and relative ("local") time may become distinct for certain purposes. A Generator object has an attribute globalt0 for use by a Model object for these purposes. The primary reason that a Generator object is set up this way is when it is re-used as a trajectory segment as part of a hybrid model. More information is provided on the page HybridSystems.

The built-in auxiliary function globalindepvar provides access to a hybrid system's absolute time from within an individual Generator (see earlier).

7.4. Setting up a Model object

7.4.1. Model objects as non-hybrid systems

To take advantage of the introspection and book-keeping utilities built into the Model class, you may wish to convert a single Generator corresponding to a non-hybrid system into a Model object. This, most trivial, use of the Model class is easily set up. Simply call the embed(<generator>, [<model_name>]) function, which returns a Model containing the single Generator. If the model name is not specified, the Generator's name plus the suffix '_model' is applied.

Because this form of Model object is not a true hybrid system, the function embed is provided to wrap individual Generators into a Model object.

7.4.2. Model objects as hybrid systems

This is discussed further on the HybridSystems page.

8. Bounds and constraints

PyDSTool supports the optional specification of bounds on variables and parameters. This is most useful when performing parameter and model estimation, but also provides a means to check the consistency of iterative processes as they run. Presently, variable domains may be singleton values (i.e., not an interval), but parameter domains must be an interval.

Bounds on variables are specified to a generator using the xdomain key, while those on parameters use the pdomain key.

These issues are described further on the BoundsSafety page.

9. Continuation and bifurcation analysis

As of v0.80, PyDSTool incorporates the sub-package PyCont, written by Drew Lamar as part of the PyDSTool group. PyCont provides a suite of continuation and bifurcation tools that can be applied to existing PyDSTool Model objects. Further information about the package can be found on the PyCont wiki page.

10. Porting models to other packages

Access to automatic differentiation via the ODETools package is provided by creating a model as an instance of the ADMC_ODEsystem class of Generators.

In the future we expect to implement utilities to export model definitions directly into formats that can be read by the original DsTool package and XPP, among others.

last edited 2008-04-30 03:18:46 by RobClewley