Add predefined timestep() function (#1100)

doc/nestml_language/nestml_language_concepts.rst

@@ -561,7 +561,10 @@ The following functions are predefined in NESTML and can be used out of the box.
      - Convert a time into a number of simulation steps. See the section :ref:`Handling of time` for more information.
    * - ``resolution``
      -
-     - Returns the current resolution of the simulation in ms. See the section :ref:`Handling of time` for more information.
+     - In the ``update`` block, or in initialising expressions, returns the current timestep taken in milliseconds. See the section :ref:`Handling of time` for more information.
+   * - ``timestep``
+     -
+     - In the ``update`` block, returns the current timestep taken in milliseconds. See the section :ref:`Handling of time` for more information.
 
 
 Predefined variables and constants
@@ -825,7 +828,7 @@ Block types
 -  ``output`` *``<event_type>``* - Defines which type of event the model can send. Currently, only ``spike`` is supported.
 -  ``update`` - Contains statements that are executed once every simulation timestep (on a fixed grid or from event to event).
 - ``onReceive`` - Can be defined for each spiking input port; contains statements that are executed whenever an incoming spike event arrives. Optional event parameters, such as the weight, can be accessed by referencing the input port name. Priorities can optionally be defined for each ``onReceive`` block; these resolve ambiguity in the model specification of which event handler should be called after which, in case multiple events occur at the exact same moment in time on several input ports, triggering multiple event handlers.
-- ``onCondition`` - Contains statements that are executed when a particular condition holds. The condition is expressed as a (boolean typed) expression. The advantage of having conditions separate from the ``update`` block is that a root-finding algorithm can be used to find the precise time at which a condition holds, within each (fixed resolution) simulation timestep. This makes the model more generic with respect to the simulator that is used.
+- ``onCondition`` - Contains statements that are executed when a particular condition holds. The condition is expressed as a (boolean typed) expression. The advantage of having conditions separate from the ``update`` block is that a root-finding algorithm can be used to find the precise time at which a condition holds (with a higher resolution than the simulation timestep). This makes the model more generic with respect to the simulator that is used.
 
 
 Input
@@ -902,7 +905,7 @@ Each model can only send a single type of event. The type of the event has to be
    output:
        spike
 
-Calling the ``emit_spike()`` function in the ``update`` block results in firing a spike to all target neurons and devices time stamped with the simulation time at the end of the time interval ``t + resolution()``.
+Calling the ``emit_spike()`` function in the ``update`` block results in firing a spike to all target neurons and devices time stamped with the simulation time at the end of the time interval ``t + timestep()``.
 
 
 Equations
@@ -1084,7 +1087,7 @@ A Dirac delta impulse kernel can be defined by using the predefined function ``d
 Handling of time
 ----------------
 
-Inside the ``update`` block, the current time can be retrieved via the predefined, global variable ``t``. The statements executed in the block are reponsible for updating the state of the model between timesteps or events. The statements in this block update the state of the model from the "current" time ``t``, to the next simulation timestep or time of next event ``t + resolution()``. The update step typically involves integration of the ODEs and corresponds to the "free-flight" or "subthreshold" integration; the events themselves are typically handled elsewhere, namely as a convolution with a kernel, or as an ``onReceive`` block.
+Inside the ``update`` block, the current time can be retrieved via the predefined, global variable ``t``. The statements executed in the block are reponsible for updating the state of the model between timesteps or events. The statements in this block update the state of the model from the "current" time ``t``, to the next simulation timestep or time of next event ``t + timestep()``. The update step typically involves integration of the ODEs and corresponds to the "free-flight" or "subthreshold" integration; the events themselves are typically handled elsewhere, namely as a convolution with a kernel, or as an ``onReceive`` block.
 
 
 Integrating the ODEs
@@ -1104,22 +1107,23 @@ Retrieving simulation timing parameters
 
 To retrieve timing parameters from the simulator kernel, two special functions are built into NESTML:
 
--  ``resolution`` returns the current resolution of the simulation in ms.
--  ``steps`` takes one parameter of type ``ms`` and returns the number of simulation steps in the current simulation resolution.
-
-These functions can be used to implement custom buffer lookup logic but should be used with care. In particular, when a non-constant simulation timestep is used, ``steps()`` should be avoided.
+- ``resolution`` returns the current timestep taken. Can be used only inside the ``update`` block and in intialising expressions. The use of this function assumes that the simulator uses fixed resolution steps, therefore it is recommended to use ``timestep()`` instead in order to make the models more generic.
+- ``timestep`` returns the current timestep taken. Can be used only inside the ``update`` block.
+- ``steps`` takes one parameter of type ``ms`` and returns the number of simulation steps in the current simulation resolution. This only makes sense in case of a fixed simulation resolution (such as in NEST); hence, use of this function is not recommended, because it precludes the models from being compatible with other simulation platforms where a non-constant simulation timestep is used.
 
-When using ``resolution()``, it is recommended to use the function call directly in the code, rather than defining it as a parameter. This makes the model more robust in case of non-constant timestep. In some cases, as in the synapse ``update`` block, a step is made between spike events, a timestep which is not constrained by the simulation timestep. For example:
+When using ``resolution()``, it is recommended to use the function call directly in the code, rather than defining it as a parameter. This makes the model more robust in case the resolution is changed during the simulation. In some cases, as in the synapse ``update`` block, a step is made between spike events, unconstrained by the simulation resolution. For example:
 
 .. code-block:: nestml
 
    parameters:
-       h ms = resolution()   # !! NOT RECOMMENDED.
+       h ms = resolution()   # !! NOT RECOMMENDED
 
    update:
-       # update from t to t + resolution()
-       x *= exp(-resolution() / tau)   # let x' = -x / tau
-                                       # evolve the state of x one timestep
+       # update from t to t + timestep()
+       # let x' = -x / tau
+       # x *= exp(-h / tau)  # !! NOT RECOMMENDED
+       # x *= exp(-resolution() / tau)  # !! better but NOT RECOMMENDED
+       x *= exp(-timestep() / tau)  # recommended (supports any timestep)
 
 
 Integration order
@@ -1138,7 +1142,7 @@ The numeric results of a typical simulation run are shown below. Consider a leak
 .. figure:: https://raw.githubusercontent.com/clinssen/nestml/integrate_specific_odes/doc/fig/integration_order_example.png
    :alt: Numerical example for two different integration sequences.
 
-On the left, both pre-synaptic spikes are only processed at the end of the interval in which they occur. The statements in the ``update`` block are run every timestep for a fixed resolution of :math:`1~\text{ms}`, alternating with the statements in the ``onReceive`` handler for the spiking input port. Note that this means that the effect of the spikes becomes visible at the end of the timestep in :math:`I_\text{syn}`, but it takes another timestep before ``integrate_odes()`` is called again and consequently for the effect of the spikes to become visible in the membrane potential. This results in a threshold crossing and the neuron firing a spike. On the right half of the figure, the same presynaptic spike timing is used, but because events are processed at their exact time of occurrence. In this case, the ``update`` statements are called once to update the neuron from time 0 to :math:`1~\text{ms}`, then again to update from :math:`1~\text{ms}` to the time of the first spike, then the spike is processed by running the statements in its ``onReceive`` block, then ``update`` is called to update from the time of the first spike to the second spike, and so on. The time courses of :math:`I_\text{syn}` and :math:`V_\text{m}` are such that the threshold is not reached and the neuron does not fire, illustrating the numerical differences that can occur when the same model is simulated using different strategies.
+On the left, both pre-synaptic spikes are only processed at the end of the interval in which they occur. The statements in the ``update`` block are run every timestep for a fixed timestep of :math:`1~\text{ms}`, alternating with the statements in the ``onReceive`` handler for the spiking input port. Note that this means that the effect of the spikes becomes visible at the end of the timestep in :math:`I_\text{syn}`, but it takes another timestep before ``integrate_odes()`` is called again and consequently for the effect of the spikes to become visible in the membrane potential. This results in a threshold crossing and the neuron firing a spike. On the right half of the figure, the same presynaptic spike timing is used, but because events are processed at their exact time of occurrence. In this case, the ``update`` statements are called once to update the neuron from time 0 to :math:`1~\text{ms}`, then again to update from :math:`1~\text{ms}` to the time of the first spike, then the spike is processed by running the statements in its ``onReceive`` block, then ``update`` is called to update from the time of the first spike to the second spike, and so on. The time courses of :math:`I_\text{syn}` and :math:`V_\text{m}` are such that the threshold is not reached and the neuron does not fire, illustrating the numerical differences that can occur when the same model is simulated using different strategies.
 
 
 Guards

doc/nestml_language/neurons_in_nestml.rst

@@ -242,31 +242,54 @@ Output
 Implementing refractoriness
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-In order to model an absolute refractory state, in which the neuron cannot fire action potentials, an extra parameter (say, ``refr_T``) can be introduced, that defines the duration of the refractory period. A new state variable, ``refr_t``, then specifies the time of the refractory period that has already elapsed, and a second boolean state variable ``is_refactory`` identifies whether or not we are in the refractory state. In the initial state, the neuron is not refractory and the timer is set to zero. When a spike is emitted, the boolean flag is set to true and the timer is set to ``refr_T``. Using a separate flag allows us to freely formulate a condition on ending the timer without having to worry about special (for instance, negative) values representing a non-refactory condition. This is hazardous because of an imprecise floating point representation of real numbers. The check against :math:`\Delta t/2` instead of checking against 0 additionally guards against accumulated discretization errors. Integrating the ODE for :math:`V_\text{m}` is disabled while the flag is set to true. When the timer reaches zero, the flag is set to false. In the ``update`` block, the timer is decremented each timestep. An ``onCondition`` is formulated on ending the refractory period, which allows the time at which the condition becomes true to be determined precisely (whereas it would be aliased to the nearest simulation timestep interval if the condition had been checked in ``update``).
+In order to model an absolute refractory state, in which the neuron cannot fire action potentials, different approaches can be used. In general, an extra parameter (say, ``refr_T``) is introduced, that defines the duration of the refractory period. A new state variable (say, ``refr_t``) can then act as a timer, counting the time of the refractory period that has already elapsed. The dynamics of ``refr_t`` could be specified in the ``update`` block, as follows:
 
 .. code-block:: nestml
 
-   parameters:
-       refr_T ms = 5 ms
+   update:
+       refr_t -= resolution()
 
-   state:
-       refr_t ms = 0 ms    # Refractory period timer
-       is_refractory boolean = false
+The test for refractoriness can then be added in the ``onCondition`` block as follows:
+
+.. code-block:: nestml
+
+   # if not refractory and threshold is crossed...
+   onCondition(refr_t <= 0 ms and V_m > V_th):
+       V_m = E_L    # Reset the membrane potential
+       refr_t = refr_T    # Start the refractoriness timer
+       emit_spike()
+
+The disadvantage of this method is that it requires a call to the ``resolution()`` function, which is only supported by fixed-timestep simulators. To write the model in a more generic way, the refractoriness timer can alternatively be expressed as an ODE:
+
+.. code-block:: nestml
 
    equations:
+       refr_t' = -1 / s    # a timer counting back down to zero
+
+Typically, the membrane potential should remain clamped to the reset or leak potential during the refractory period. It depends on the intended behavior of the model whether the synaptic currents and conductances also continue to be integrated or whether they are reset, and whether incoming spikes during the refractory period are taken into account or ignored.
+
+In order to hold the membrane potential at the reset voltage during refractoriness, it can be simply excluded from the integration call:
+
+.. code-block:: nestml
+
        I_syn' = ...
        V_m' = ...
+       refr_t' = -1 / s    # Count down towards zero
 
    update:
-       if is_refractory:
+       if refr_t > 0 ms:
            # neuron is absolute refractory, do not evolve V_m
-           refr_t -= resolution()
-           integrate_odes(I_syn)
+           integrate_odes(I_syn, refr_t)
        else:
-           # neuron not refractory, so evolve all ODEs
-           integrate_odes(V_m, I_syn)
+           # neuron not refractory
+           integrate_odes(I_syn, V_m)
+
+Note that in some cases, the finite resolution by which real numbers are expressed (as floating point numbers) in computers, can cause unexpected behaviors. If the simulation resolution is not exactly representable as a float (say, Δt = 0.1 ms) then it could be the case that after 20 simulation steps, the timer has not reached zero, but a very small value very close to zero (say, 0.00000001 ms), causing the refractory period to end only in the next timestep. If this kind of behavior is undesired, the simulation resolution and refractory period can be chosen as powers of two (which can be represented exactly as floating points), or a small "epsilon" value can be included in the comparison in the model:
+
+.. code-block:: nestml
+
+   parameters:
+       float_epsilon ms = 1E-9 ms
 
-   onCondition(is_refractory and refr_t <= resolution() / 2):
-       # end of refractory period
-       refr_t = 0 ms
-       is_refractory = false
+   onCondition(refr_t <= float_epsilon ...):
+       # ...