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).
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).
These are simply string substitutions.
The name of the variable (tau) is changed in the string to taum.
The name of the variable (tau) is changed in the string to a unique identifier.
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.
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).
After initialisation, an Equation object contains:
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.
The Equation object is finalised by an explicit call to the prepare() method.
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).
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()).
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:
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.
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).
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.
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.
Functions are then recompiled so that differential equations are now independent of static variables.
Finally, the list of undefined identifiers is checked (free_variables()) and a warning is issued if any is found.
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 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.