1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328

.. 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=(yx)/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
.. _numericalintegration:
Numerical integration

The currently available integration methods are:
* Exact integration when the equations are linear.
* Euler integration (explicit, first order).
* RungeKutta 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 secondorder method
can be selected using the keyword ``order=2`` (explicit RungeKutta 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; semiexact
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(XB) 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(XB),
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; HodginHuxley type equations
Exponential Euler integration
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The exponential Euler method is used for HodgkinHuxley 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 nontrivial, 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 userdefined 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; nonautonomous
pair: differential equations; timedependent
pair: equations; timedependent
pair: equations; nonautonomous
Nonautonomous 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 speedup.
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 speedup (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=(yx)/(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 righthand 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.
