Solving differential equations with the GNU Scientific Library

Conventionally, Brian generates its own code performing Numerical integration according to the chosen algorithm (see the section on Code generation). Another option is to let the differential equation solvers defined in the GNU Scientific Library (GSL) solve the given equations. In addition to offering a few extra integration methods, the GSL integrator comes with the option of having an adaptable timestep. The latter functionality can have benefits for the speed with which large simulations can be run. This is because it allows the use of larger timesteps for the overhead loops in Python, without losing the accuracy of the numerical integration at points where small timesteps are necessary. In addition, a major benefit of using the ODE solvers from GSL is that an estimation is performed on how wrong the current solution is, so that simulations can be performed with some confidence on accuracy. (Note however that the confidence of accuracy is based on estimation!)

StateUpdateMethod

Translation of equations to abstract code

The first part of Brian’s code generation is the translation of equations to what we call ‘abstract code’. In the case of Brian’s stateupdaters so far, this abstract code describes the calculations that need to be done to update differential variables depending on their equations as is explained in the section on State update. In the case of preparing the equations for GSL integration this is a bit different. Instead of writing down the computations that have to be done to reach the new value of the variable after a time step, the equations have to be described in a way that GSL understands. The differential equations have to be defined in a function and the function is given to GSL. This is best explained with an example. If we have the following equations (taken from the adaptive threshold example):

dv/dt = -v/(10*ms) : volt
dvt/dt = (10*mV - vt)/(15*ms) : volt

We would describe the equations to GSL as follows:

v = y[0]
vt = y[1]
f[0] = -v/(10e-3)
f[1] = (10e-3 - vt)

Each differential variable gets an index. Its value at any time is saved in the y-array and the derivatives are saved in the f-array. However, doing this translation in the stateupdater would mean that Brian has to deal with variable descriptions that contain array accessing: something that for example sympy doesn’t do. Because we still want to use Brian’s existing parsing and checking mechanisms, we needed to find a way to describe the abstract code with only ‘normal’ variable names. Our solution is to replace the y[0], f[0], etc. with a ‘normal’ variable name that is later replaced just before the final code generation (in the GSLCodeGenerator). It has a tag and all the information needed to write the final code. As an example, the GSL abstract code for the above equations would be:

v = _gsl_y0
vt = _gsl_y1
_gsl_f0 = -v/(10e-3)
_gsl_f1 = (10e-3 - vt)

In the GSLCodeGenerator these tags get replaced by the actual accessing of the arrays.

Return value of the StateUpdateMethod

So far, for each each code generation language (numpy, cython) there was just one set of rules of how to translate abstract code to real code, described in its respective CodeObject and CodeGenerator. If the target language is set to Cython, the stateupdater will use the CythonCodeObject, just like other objects such as the StateMonitor. However, to achieve the above decribed translations of the abstract code generated by the StateUpdateMethod, we need a special CythonCodeObject for the stateupdater alone (which at its turn can contain the special CodeGenerator), and this CodeObject should be selected based on the chosen StateUpdateMethod.

In order to achieve CodeObject selection based on the chosen stateupdater, the StateUpdateMethod returns a class that can be called with an object, and the appropriate CodeObject is added as an attribute to the given object. The return value of this callable is the abstract code describing the equations in a language that makes sense to the GSLCodeGenerator.

GSLCodeObject

Each target language has its own GSLCodeObject that is derived from the already existing code object of its language. There are only minimal changes to the already existing code object:

  • Overwrite stateupate template: a new version of the stateupdate template is given (stateupdate.cpp for C++ standalone and stateupdate.pyx for cython).

  • Have a GSL specific generator_class: GSLCythonCodeGenerator

  • Add the attribute original_generator_class: the conventional target-language generator is used to do the bulk of the translation to get from abstract code to language-specific code.

This defining of GSL-specific code objects also allowed us to catch compilation errors so we can give the user some information on that it might be GSL-related (overwriting the compile() method in the case of cython). In the case of the C++ CodeObject such overriding wasn’t really possible so compilation errors in this case might be quite undescriptive.

GSLCodeGenerator

This is where the magic happens. Roughly 1000 lines of code define the translation of abstract code to code that uses the GNU Scientific Library’s ODE solvers to achieve state updates.

Upon a call to run(), the code objects necessary for the simulation get made. The code for this is described in the device. Part of making the code objects is generating the code that describes the code objects. This starts with a call to translate, which in the case of GSL brings us to the GSLCodeGenerator.translate(). This method is built up as follows:

  • Some GSL-specific preparatory work:

    • Check whether the equations contain variable names that are reserved for the GSL code.

    • Add the ‘gsl tags’ (see section on StateUpdateMethod) to the variables known to Brian as non-scalars. This is necessary to ensure that all equations containing ‘gsl tags’ are considered vector equations, and thus added to Brian’s vector code.

    • Add GSL integrator meta variables as official Brian variables, so these are also taken into account upon translation. The meta variables that are possible are described in the user manual (e.g. GSL’s step taken in a single overhead step ‘_step_count’).

    • Save function names. The original generators delete the function names from the variables dictionary once they are processed. However, we need to know later in the GSL part of the code generation whether a certain encountered variable name refers to a function or not.

  • Brian’s general preparatory work. This piece of code is directly copied from the base CodeGenerator and is thus similar to what is done normally.

  • A call to original_generator.translate() to get the abstract code translated into code that is target-language specific.

  • A lot of statements to translate the target-language specific code to GSL-target-language specific code, described in more detail below.

The biggest difference between conventional Brian code and GSL code is that the stateupdate-decribing lines are contained directly in the main() or in a separate function, respectively. In both cases, the equations describing the system refer to parameters that are in the Brian namespace (e.g. “dv/dt = -v/tau” needs access to “tau”). How can we access Brian’s namespace in this separate function that is needed with GSL?

To explain the solution we first need some background information on this ‘separate function’ that is given to the GSL integrators: _GSL_func. This function always gets three arguments:

  • double t: the current time. This is relevant when the equations are dependent on time.

  • const double _GSL_y[]’: an array containing the current values of the differential variables (const because the cannot be changed by _GSL_func itself).

  • double f[]: an array containing the derivatives of the differential variables (i.e. the equations describing the differential system).

  • void * params: a pointer.

The pointer can be a pointer to whatever you want, and can thus point to a data structure containing the system parameters (such as tau). To achieve a structure containing all the parameters of the system, a considerable amount of code has to be added/changed to that generated by conventional Brian:

  • The data structure, _GSL_dataholder, has to be defined with all variables needed in the vector code. For this reason, also the datatype of each variable is required.

    • This is done in the method GSLCodeGenerator.write_dataholder

  • Instead of referring to the variables by their name only (e.g. dv/dt = -v/tau), the variables have to be accessed as part of the data structure (e.g. dv/dt = -v/_GSL_dataholder->tau in the case of cpp). Also, as mentioned earlier, we want to translate the ‘gsl tags’ to what they should be in the final code (e.g. _gsl_f0 to f[0]).

    • This is done in the method GSLCodeGenerator.translate_vector_code. It works based on the to_replace dictionary (generated in the methods GSLCodeGenerator.diff_var_to_replace and GSLCodeGenerator.to_replace_vector_vars) that simply contains the old variables as keys and new variables as values, and is given to the word_replace function.

  • The values of the variables in the data structure have to be set to the values of the variables in the Brian namespace.

    • This is done in the method GSLCodeGenerator.unpack_namespace, and for the ‘scalar’ variables that require calculation first it is done in the method GSLCodeGenerator.translate_scalar_code.

In addition, a few more ‘support’ functions are generated for the GSL script:

  • int _set_dimension(size_t * dimension): sets the dimension of the system. Required for GSL.

  • double* _assign_memory_y(): allocates the right amount of memory for the y array (also according to the dimension of the system).

  • int _fill_y_vector(_dataholder* _GSL_dataholder, double* _GSL_y, int _idx): pulls out the values for each differential variable out of the ‘Brian’ array into the y-vector. This happens in the vector loop (e.g. y[0] = _GSL_dataholder->_ptr_array_neurongroup_v[_idx]; for C++).

  • int _empty_y_vector(_dataholder* _GSL_dataholder, double* _GSL_y, int _idx): the opposite of _fill_y_vector. Pulls final numerical solutions from the y array and gives it back to Brian’s namespace.

  • double* _set_GSL_scale_array(): sets the array bound for each differential variable, for which the values are based on method_options['absolute_error'] and method_options['absolute_error_per_variable'].

All of this is written in support functions so that the vector code in the main() can stay almost constant for any simulation.

Stateupdate templates

There is many extra things that need to be done for each simulation when using GSL compared to conventional Brian stateupdaters. These are summarized in this section.

Things that need to be done for every type of simulation (either before, in or after main()):

  • Cython-only: define the structs and functions that we will be using in cython language.

  • Prepare the gsl_odeiv2_system: give function pointer, set dimension, give pointer to _GSL_dataholder as params.

  • Allocate the driver (name for the struct that contains the info necessary to perform GSL integration)

  • Define dt.

Things that need to be done every loop iteration for every type of simulation:

  • Define t and t1 (t + dt).

  • Transfer the values in the Brian arrays to the y-array that will be given to GSL.

  • Set _GSL_dataholder._idx (in case we need to access array variables in _GSL_func).

  • Initialize the driver (reset counters, set dt_start).

  • Apply driver (either with adaptable- or fixed time step).

  • Optionally save certain meta-variables

  • Transfer values from GSL’s y-vector to Brian arrays