Solving differential equations with the GNU Scientific Library
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Conventionally, Brian generates its own code performing :doc:`../user/numerical_integration`
according to the chosen algorithm (see the section on :doc:`codegen`).
Another option is to let the differential equation solvers defined in the
`GNU Scientific Library (GSL) <https://www.gnu.org/software/gsl/doc/html/ode-initval.html>`_
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 :doc:`../advanced/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):

.. code-block:: none

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

We would describe the equations to GSL as follows:

.. code-block:: C

   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:

.. code-block:: C

    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