Equations

An Equation is a set of single lines in a string:
  1. dx/dt = f : unit (differential equation)
  2. x = f : unit (equation)
  3. x = y (alias)
  4. x : unit (parameter)

The equations may be defined on multiple lines with the character . Comments using # may also be included.

Two special variables are defined: t (time) and xi (white noise). Ultimately, it should be possible (using Sympy) to define equations implicitly, e.g.: ‘tau*dv/dt=-v : unit’ (although it makes unit specification ambiguous).

An equation can be seen as a set of functions or code and a namespace to evaluate them. A key part of object construction is the construction of the namespace (after parsing).

Namespace construction

The namespaces are stored in eq._namespace. Each equation (string) has a specific namespace.

Proposition for a simplification: there could be just one namespace per Equation object rather than per string. Possible conflicts would be dealt with when equations are added (with prefix as when inserting static variables, see below).

Variable substitution

These are simply string substitutions.

  • Equation(‘dv/dt=-v/tau:volt’,tau=’taum’)

The name of the variable (tau) is changed in the string to taum.

  • Equation(‘dv/dt=-v/tau:volt’,tau=None)

The name of the variable (tau) is changed in the string to a unique identifier.

Explicit namespace

  • Equation(‘dv/dt=-v/tau:volt’,tau=2*ms)

The namespace is explicitly given: {‘tau’:2*ms}. In this case, Brian does not try to build a namespace “magically”, so the namespace must be exhaustive. Units need not be passed.

Implicit namespace

  • Equation(‘dv/dt=-v/tau:volt’)

The namespace is built from the globals and locals in the caller’s frame. For each identifier in the string, the name is looked up in: 1) locals of caller, 2) globals of caller, 3) globals of equations.py module (typically units). Identifiers can be any Python object (for example functions).

Issues

  • Special variables (xi and t) are not taken into account at this stage, i.e., they are integrated in the namespace if present. This should probably be fixed and a warning should be raised. A warning is raised for t at the preparation stage (see below).
  • If identifiers are not found, then no error is raised. This is to allow equations to be built in several pieces, which is useful in particular for library objects.
  • If an identifier is found whose name is the same as the name of a variable, then no error is raised here and it is included in the namespace. This is difficult to avoid in the case when equations are built in several pieces (e.g. the conflict appears only when the pieces are put together). A warning is issued at the preparation stage (see below).

Attributes after initialisation

After initialisation, an Equation object contains:

  • a namespace (_namespace)
  • a dictionary of units for all variables (_units)
  • a dictionary of strings corresponding to each variable (right hand side of each equation), including parameters and aliases (_string). Parameters are defined as differential equations with RHS 0*unit/second. All comments are removed and multiline strings are concatenated.
  • a list of variables of non-differential equations (_eq_names)
  • a list of variables of differential equations, including parameters (_diffeq_names)
  • a list of variables of differential equations, excluding parameters (_diffeq_names_nonzero)
  • a dictionary of aliases (_alias), mapping a variable name to its alias

There is no explicit list of parameters, maybe it should be added. Nothing more is done at initialisation time (no units checking, etc). The reason is that the equation set might not be complete at this time, in the case when equations are built in several pieces. Various checks are done using the prepare() method.

Finalisation (prepare())

The Equation object is finalised by an explicit call to the prepare() method.

Finding Vm

The first step is to find the name of the membrane potential variable (getVm()). This is useful when the variable name for threshold or reset is not given (e.g. threshold=10*mV). The method looks for one these names: ‘v’,’V’,’vm’,’Vm’. If one is present, it is the membrane potential variable. If none or more than one is present, no variable is found. If it is found, the corresponding variable is swapped with the first variable in the _diffeq_names list (note: not in the _diffeq_names_nonzero list). Otherwise, nothing happens. This way, the first variable in the list is the membrane potential. Possibly, a warning could be issued if it is not found. The problem it might issue warnings too often. A better way would be to issue warnings only when threshold and reset are used ambiguously (i.e., no Vm found and more than 1 variable).

Cleaning the namespace

Then variables and t are removed from the namespace if present (N.B.: xi does not appear to be removed), and warnings are issued using log_warn (method clean_namespace()).

Compiling functions

This is done by the compile_functions() method. Python functions are created from the string definition of equations. For each equation/differential equation, the list of identifiers is obtained from the string definition, then only those referring to variables are kept. A Python lambda function of these remaining identifiers is then compiled (using eval) and put in the _function dictionary.

Compiled functions are used for:

  • checking units
  • obtaining the list of arguments (this could be done independently)
  • state updates

This step might be avoided and replaced by eval calls. It might actually be a little simpler because arguments would be replaced by namespace. It seems to be faster with the current implementation, but the string could be compiled with compile() (then evaluated in the relevant namespace). Besides, with the way it is currently evaluated in the Euler update: f(*[S[var] for var in f.func_code.co_varnames]), it is not faster than direct evaluation in the namespace.

Checking units

This is done by the check_units() method. First, the static equations are ordered (see next section).

To check the units of a static equation, one calls the associated function (giving the RHS) where the arguments are units (e.g., 1*volt for v, etc.) and adds the units of the LHS. A dimension error is raised if it is not homogeneous. Currently, the message states “The differential equation is not homogeneous” but it should be adapted to non-differential equations. One problem with this way of checking units is that the RHS function may not be defined at the point it is checked.

Differential equations are checked in the same way, with two specificities: the units of RHS should be the units of the variable divided by second (dx/dt), and noise (xi) has units of second**-.5 (this is put in the globals of the function, which might not be a very clean way to do it).

Ordering static variables

It seems that this method (set_eq_order()) is already called by check_units() and therefore it is probably not necessary to call it here. This method computes the dependency graph of (static) equations on other static variables, which must have no cycle (otherwise an error is raised). From that graph, an update list is built and put in _eq_names. Then for each variable (static or differential), the list of dependent static variables is built and sorted in update order. The result is put in the _dependencies dictionary.

This is a necessary step to calculate the RHS of any equation: it gives the ordered list of static variables to calculate first before calculating the RHS.

Inserting static variables into differential equations

The value of static variables are then replaced by their string value (RHS) in all differential equations (substitute_eq()). The previous step (ordering) ensures that the result if correct and does not depend on static variables anymore. To avoid namespace conflicts, all identifiers in the namespace of a static variable is augmented by a prefix: name+’_’ (e.g. ‘x_y’ for identifier y in equation ‘x=2*y’). Then namespaces are merged.

It might not be optimal to do it in this way, because some of calculations will be done several times in an update step. It might be better to keep the static variables separate.

Recompiling functions

Functions are then recompiled so that differential equations are now independent of static variables.

Checking free variables

Finally, the list of undefined identifiers is checked (free_variables()) and a warning is issued if any is found.

Freezing

Freezing is done by calling compile_functions(freeze=True). Each string expression is then frozen with optimiser.freeze(), which replaces identifiers by their float value. This step does not necessarily succeed, in which case a warning (not an error) is issued.

Adding Equation objects

Adding equations consists simply in merging the lists/dictionaries of variables, namespaces, strings, units and functions. Conflicts raise an error. This step must precede preparation of the object.