.. currentmodule:: brian

.. index:: equations

.. _moreonequations:

More on equations
=================

The :class:`Equations` class is a central part of Brian, since
models are generally specified with an :class:`Equations` object.
Here we explain advanced aspects of this class.

.. index::
	pair: equations; external variables
	pair: equations; namespaces

External variables
------------------
Equations may contain external variables. When an :class:`Equations` object is initialised,
a dictionary is built with the values of all external variables. These values are taken
from the namespace where the :class:`Equations` object was defined. It is possible to go one or
several levels up in the namespaces by specifying the keyword ``level`` (default=0).
The value of these parameters can in general be changed during the simulation and the
modifications are taken into account, except in two situations: when the equations are
frozen (see below) or when the integration is exact (linear equations). In those cases,
the values of the parameters are the ones at initialisation time.

Alternatively, the string defining the equations can be evaluated within a given namespace
by providing keywords at initialisation time, e.g.::

  eqs=Equations('dx/dt=-x/tau : volt',tau=10*ms)
  
In that case, the values of all external variables are taken from the specified
dictionary (given by the keyword arguments), even if variables with the same name
exist in the namespace where the string was defined. The two methods for passing the
values of external variables are mutually exclusive, that is, either all external variables
are explicitly specified with keywords (if not, they are left unspecified even if there
are variables with the same names in the namespace where the string was defined), or all
values are taken from the calling namespace.

More can be done with keyword arguments. If the value is a string, then the name of the
variable is replaced, e.g.::

  eqs=Equations('dx/dt=-x/tau : volt',tau=10*ms,x='Vm')

changes the variable name x to Vm. This is useful for writing functions which return
equations where the variable name is provided by the user.

Finally, if the value is ``None`` then the name of the variable is replaced by a unique
name, e.g.::

  eqs=Equations('dx/dt=-x/tau : volt',tau=10*ms,x=None)

This is useful to avoid conflicts in the names of hidden variables.

Issues
^^^^^^
* There can be problems if a variable with the same name as the variable of a
  differential equation exists in the namespace where the :class:`Equations` object was defined.

.. index::
	pair: equations; combining

Combining equations
-------------------
:class:`Equations` can be combined using the sum operator. For example::

  eqs=Equations('dx/dt=(y-x)/tau : volt')
  eqs+=Equations('dy/dt=-y/tau: volt')
  
Note that some variables may be undefined when defining the first equation. No error is
raised when variables are undefined and absent from the calling namespace.
When two :class:`Equations` objects are added, the consistency is checked. For example it is not
possible to add two :class:`Equations` objects which define the same variable.

.. index::
	pair: equations; membrane potential

Which variable is the membrane potential?
-----------------------------------------
Several objects, such as :class:`Threshold` or :class:`Reset` objects
can be initialised without specifying which variable is the membrane potential, in which
case it is assumed that it is the first variable.
Internally, the variables of an :class:`Equations` object are reorderered so that the first one
is most likely to be the membrane potential (using :meth:`Equations.get_Vm`).
The first variable is, with decreasing priority :

* v
* V
* vm
* Vm
* the first defined variable.

.. index::
	pair: equations; numerical integration

.. _numerical-integration:

Numerical integration
---------------------
The currently available integration methods are:

* Exact integration when the equations are linear.
* Euler integration (explicit, first order).
* Runge-Kutta integration (explicit, second order).
* Exponential Euler integration (implicit, first order).

The method is selected when a :class:`NeuronGroup` is initialized.
If the equations are linear, exact integration is automatically selected.
Otherwise, Euler integration is selected by default, unless the keyword
``implicit=True`` is passed, which selects the exponential Euler method. A second-order method
can be selected using the keyword ``order=2`` (explicit Runge-Kutta method, midpoint estimation).
It is possible to override this behaviour with the ``method`` keyword when initialising
a :class:`NeuronGroup`. Possible values are ``linear``, ``nonlinear``,
``Euler``, ``RK``, ``exponential_Euler``.

.. index::
	pair: equations; linear
	pair: numerical integration; exact
	pair: numerical integration; semi-exact

Exact integration
^^^^^^^^^^^^^^^^^
If the differential equations are linear, then the update phase
X(t)->X(t+dt) can be calculated exactly with a matrix product.
First, the equations are examined to determine whether they are linear
with the method :meth:`~equations.Equations.islinear` and the function
:func:`~inspection.is_affine` (this is currently done using dynamic typing).
Second, the matrix M and the vector B such that dX/dt=M(X-B) are calculated with
the function :func:`~stateupdater.get_linear_equations` [1]_.
Third, the matrix A such that X(t+dt)=A*(X(t)-B)+B is calculated at initialisation
of a specific state updater object, :class:`LinearStateUpdater`,
as A=expm(M*dt), where expm is the matrix exponential.

**Important remark**: since the update matrix and vector are precalculated,
the values of all external variables in the equations are frozen at
initialisation. If external variables are modified after initialisation,
those modifications are *not* taken into account during the simulation.

**Inexact exact integration**: If the equation cannot be put into the form dX/dt=M(X-B),
for example if the equation is dX/dt=MX+A where M is not invertible, then the equations
are not integrated exactly, but using a system equivalent to Euler integration but with
dt 100 times smaller than specified. Updates are of the form X(t+dt)=A*X(t)+C where the
matrix A and vector C are computed by applying Euler integration 100 times to the
differential equations.

.. index::
	pair: numerical integration; Euler

Euler integration
^^^^^^^^^^^^^^^^^
The Euler is a first order explicit integration method. It is the default one for
nonlinear equations. It is simply implemented as X(t+dt)=X(t)+f(X)*dt.

.. index::
	pair: numerical integration; exponential Euler
	pair: numerical integration; Hodgin-Huxley type equations

Exponential Euler integration
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The exponential Euler method is used for Hodgkin-Huxley type
equations, are which stiff.
Equations of that type are conditionally linear, that is, the differential equation
for each variable is linear in that variable (i.e., linear if all other variables
are considered constant).
The idea is thus to solve the differential equation for each variable over one time step,
assuming that all other variables are constant over that time step. The numerical scheme
is still first order, but it is more stable than the forward Euler method.
Each equation can be written as dx/dt=a*x+b, where a and b depend on the other variables
and thus change after each time step.
The values of a and b are obtained during the update phase by calculating a*x+b for x=0 and x=1
(note that these values are different for every neuron, thus we calculate vectors A and B).
Then x(t+dt) is calculated in the same way as for the exact integration method above.

.. index::
	single: equations; stochastic
	pair: differential equations; stochastic

Stochastic differential equations
---------------------------------
Noise is introduced in differential equations with the keyword
``xi``, which means normalised gaussian noise (the derivative of the Brownian term).
Currently, this is implemented simply by adding a normal random number to the variable
at the end of the integration step (independently for each neuron).
The unit of white noise is non-trivial, it is ``second**(-.5)``. Thus, a typical stochastic
equation reads::

  dx/dt=-x/tau+sigma*xi/tau**.5

where ``sigma`` is in the same units as ``x``. We note the following two facts:

* The noise term is independent between neurons. Thus, one cannot use this method to analyse
  the response to frozen noise (where all neurons receive the same input noise). One would need
  to use an external variable representing the input, updated by a user-defined operation.
* The noise term is independent between equations. This can however be solved by the following
  trick::
  
    dx/dt=-x/tau+sigmax*u/tau**.5 : volt
    dy/dt=-y/tau+sigmay*u/tau**.5 : volt
    u=xi : second**(-.5)

Important notice
^^^^^^^^^^^^^^^^
It is not possible to modulate the noise term with a variable (e.g.
``v*xi/tau**.5``). One reason is that, with multiplicative noise, there is
an ambiguity between the Ito and the Stratonovich interpretation.
Unfortunately, this limitation also applies to parameters, i.e.,
``sigma*xi/tau**.5`` is not possible if sigma is a parameter, as in the
following example::

  eqs=Equations('''dx/dt=-x/tau + sigma*xi/tau**.5 : volt
                   sigma : volt''')
  group = NeuronGroup(1, eqs, threshold='x>vt')

However, the problem can usually be solved by some rewriting::

  eqs=Equations('''dy/dt=-y/tau + xi/tau**.5 : 1
                    x=sigma*y : volt
                    sigma : volt''')
  group = NeuronGroup(1, eqs, threshold = 'x>vt')

.. index::
	pair: differential equations; non-autonomous
	pair: differential equations; time-dependent
	pair: equations; time-dependent
	pair: equations; non-autonomous

Non-autonomous equations
------------------------
The time variable ``t`` can be directly inserted into an equation string.
It is replaced at run time by the current value of the time variable for the relevant
neuron group, and also appears as a state variable of the neuron group.

.. index::
	pair: equations; freezing
	pair: differential equations; freezing

Freezing
--------
External variables can be frozen by passing the keyword ``freeze=True``
(default = ``False``) at initialization of a :class:`NeuronGroup` object.
Then when the string defining the equations are compiled into Python functions
(method :meth:`~equations.Equations.compile_functions`),
the external variables are replaced by their float values (units are discarded).
This can result in a significant speed-up.

TODO: more on the implementation.

.. index::
	pair: equations; compilation
	pair: differential equations; compilation

Compilation
-----------
State updates can be compiled into Python code objects by passing the keyword
``compile=True`` at initialization of a a :class:`NeuronGroup`.
Note that this is different from the method :meth:`~equations.Equations.compile_functions`,
which compiles the equation for every variable into a Python function
(not the whole state update process).

When the ``compile`` keyword is set, the method :meth:`~equations.Equations.forward_euler_code`
or :meth:`~equations.Equations.exponential_euler_code` is called. It generates
a string containing the Python code for the update of all state variables (one time step),
then compiles it into Python code object. That compiled object is then called at every time step.
All external variables are frozen in the process (regardless of the value of the ``freeze`` keyword).
This results in a significant speed-up (although the exponential Euler code is not
quite optimised yet). Note that only Python code is generated, thus a
C compiler is not required.

Working with equations
----------------------
:class:`Equations` object can also be used outside simulations.
In the following, we suppose that an :class:`Equations` object is defined as follows::

  eqs=Equations('''
  dx/dt=(y-x)/(10*ms) : volt
  dy/dt=-z/(5*ms) : volt
  z=2*(x+y) : volt
  ''')

.. index::
	pair: equations; applying

Applying an equation
^^^^^^^^^^^^^^^^^^^^
The value of z can be calculated using the :meth:`~equations.Equations.apply` method::

  z=eqs.apply('z',dict(x=3*mV,y=5*mV))
  
The second argument is a dictionary containing the values of all dependent variables
(here the result is ``8*mV``).
The right-hand side of differential equations can also be calculated in the same way::

  x=eqs.apply('x',dict(x=2*mV,y=3*mV))
  y=eqs.apply('y',dict(x=2*mV,y=3*mV))

Note in the second case that only the values of the dynamic variables should be passed.  

.. index::
	pair: equations; fixed points

Calculating a fixed point
^^^^^^^^^^^^^^^^^^^^^^^^^
A fixed point of the equations can be calculated as follows::

  fp=eqs.fixedpoint(x=2*mV,y=3*mV)

where the optional keywords give the initial point (zero if not provided).
Internally, the function ``optimize.fsolve`` from the Scipy package is used to
find a zero of the set of differential equations (thus, convergence is not
guaranteed; in that case, the initial values are returned).
A dictionary with the values of the dynamic variables at the fixed point is returned.

Issues
^^^^^^
* If the equations were previously frozen, then the units disappear from the equations
  and unit consistency problems may arise.
* :class:`Equations` objects need to be "prepared" before use, as follows::

    eqs.prepare()
    
  This is automatically called by the NeuronGroup initialiser.

.. rubric:: Footnotes

.. [1] Note that this approach raises an issue when dX/dt=B. We currently (temporarily)
   solve this problem by adding a small diagonal matrix to M to make it invertible.