Rename ``resolution()`` to ``timestep()``/floating point epsilon

[clinssen]

Depends on #879.

The update block can be used to integrate the subthreshold dynamics of the neuron across one timestep. However, the timestep might not be constant, for instance when integrating from spike to spike in a synapse model. The word "resolution" implies a fixed timestep. Hence I would suggest to rename it to "timestep".

[clinssen]

An issue with this is the use of resolution() in an onCondition expression, for instance in refractoriness:

onCondition(is_refractory and refr_t <= resolution() / 2):

What should the timestep refer to here in case of non-constant simulation resolution? Possibly, NESTML should have a predefined variable "epsilon", which could for instance correspond to the C++ std::numeric_limits<float>::epsilon().

[heplesser]

An issue with this is the use of resolution() in an onCondition expression, for instance in refractoriness:

onCondition(is_refractory and refr_t <= resolution() / 2):

What should the timestep refer to here in case of non-constant simulation resolution? Possibly, NESTML should have a predefined variable "epsilon", which could for instance correspond to the C++ std::numeric_limits<float>::epsilon().

I wonder what mathematical concept you want to express with refr_t < resolution() / 2 here, as you in the remainder of the comment mention some epsilon?

[clinssen]

@heplesser: you came up with this solution initially, as we wanted to refactor the refractoriness mechanism from an integer-based countdown method to a floating point countdown method, to make models more generic in supporting simulation platforms that have a non-constant timestep. Because of floating point rounding errors, we introduced an epsilon equal to resolution/2. We could possibly replace this with a numeric_floats::epsilon. Do you think that would be an adequate solution?

[heplesser]

resolution is clearly a concept that makes sense only for fixed-time-step simulation schemes. Since NESTML aims to be generic and describe models, not implementations, resolution is indeed not a suitable term.

At the same time, I think that timestep is not a suitable term either, as it is also an implementation detail. In the model specification, strictly speaking neither should appear, although one might consider implementation hints, e.g., for fixed-time-step models.

The underlying challenge here is how to reliably implement discrete events such as return from refractoriness. Let us assume that $t_r$ is the precise, model-defined time at which the refractory period of the neuron ends and that we have a sequence of update time points $t_1 < t_2 < t3$ with no further update events in between them (except possibly return from refractoriness). These time points might be times at which incoming spikes arrive, outgoing spikes are generated or times on a fixed update grid. Updates are defined over the left-open, right-closed intervals $(t_1, t_2]$ and $(t_2, t_3]$. If $t_r\in (t_1, t_2]$, we can distinguish two cases (I am not sure how to specify in NESTML which approach to use)

• If we have precise-timing style updates (see Morrison et al (2007)), we need to break the first update interval into $(t_1, t_r]$ and $(t_r, t_2]$ with the neuron refractory only in the first part. In this case, we have no rounding problems, we treat $t_r$ simply as a double.
• If we switch refractoriness only at otherwise existing update times, then small numerical errors in calculating t_r and t_2 could lead to delayed return from refractoriness. But this is an implementation rather than a model specification problem:
  • In a step-based simulation scheme, all times should be expressed as discrete steps (as NEST does it classically) and the rounding problem disappears.
  • For continuous-time schemes (including the precise-timing case above) one can consider to allow for some numerical jitter to avoid confusing results, although such measures may again lead to confusing results if a user injects carefully crafter spike trains. But an approach here could be to allow a relative error of a few bits of precision.

Now for the specification of refractoriness in general, I wonder if the onCondition() expression above is sensible at all, or if it would not be better to express refractoriness in a declarative way as follows

1. At any time, $t_s$ is the time of the last spike and $t_r = t_s+\tau_r$ is the end of the refractory period following that spike.
2. The neuron is thus refractory in the period $(t_s, t_r]$ and this is expressed in NESTML as is_refractory = t_s < t <= t_s + tau_r

I am not entirely happy with writing ... < ... <= ... since the user might deviate from the left-open, right-closed logic, so

is_refractory = t in Interval(t_s, t_s + tau_r)

could be a better solution, where Interval() always is left-open, right-closed.

[clinssen]

The discussion of whether a floating point epsilon has also come up in a discussion about the ignore_and_fire neuron.

When the simulation resolution is 1 ms and the rate is set to 100/s, the actual number of spikes achieved over a 1 second interval is not 100, but 91. This is because it takes 10 iterations to update the phase to 0.99999999999999989. Only on the next iteration (with a phase value close to 1.1) does the threshold check get triggered and a spike get fired, so we are only firing a spike every 11 iterations rather than every 10.

Clearly, a floating point epsilon primitive in the NESTML language would help here.

@tomtetzlaff argues against the inclusion of a floating point epsilon in NESTML: "At the level of the model description, this would be confusing because, mathematically, i.e. by solving the ODE for the phase, the firing rate of the defined model is correct." I would tend to agree as we define numbers as "real" in NESTML, not as "float" or "double".

On the other hand, "the solution requires the user to think about potential problems of a numerical implementation of the model -- something any modeler should be aware of." (@tomtetzlaff)

@diesmann: "This problems of the floating point representation are known to every scientist and it is better to rather expose than to hide them."

Automatically adjusting time comparisons

Alternatively, we should evaluate whether floating point comparisons can automatically be made safe in a generic way, by including the epsilon when generating code in a fully automated fashion. For instance, to mark time-like variables during the definition of state variables or parameters as "time-like variables", like this:

duration    ms = 100. ms   is_time

This would inform the compiler to perform the following replacements:

`timer == duration`  ->  `abs(timer-duration) < epsilon`
`timer <  duration`  ->  `timer < duration -  epsilon`
`timer >  duration`  ->  `timer > duration +  epsilon`
`timer <= duration`  ->  `timer < duration +  epsilon`
`timer >= duration`  ->  `timer < duration -  epsilon`

Use of global time variable

In order to implement an interval check, we have different options:

• Decrease a timer by timestep() inside the update block
• Implement the timer as an ODE
• Write a start time t_start to a local variable (e.g. start of the refractory period), then compare the current time t with respect to t >= t_start + interval_duration. @diesmann suggests that we should do benchmarking, as well as that this could have consequences for neuromorphic hardware, as we might now need access to the global time t at each neuron or synapse model.

[clinssen]

@tomtetzlaff provided the following minimal reproducer script.

To illustrate the floating-point-precision issue with a minimal model, I wrote the attached nestml model and test script. The model implements a timer, similar to the timer we use for refractoriness, but without all the overhead we have in normal neuron models. Executing the test script should produce the figure attached to this email.

The timer starts counting at the specified starting time tstart (here 1ms) by means of the simple ODE timer' = 1. After each time step (vertical gray lines in the figure), this ODE increases the timer by the length of a time step dt (here 0.1ms). When the timer reaches the desired duration D (here, D=1ms), a spike is emitted and the timer is instantly reset to zero. With the parameters chosen in the attached example, the first spike should hence occur at time t=2ms. However, as you can see in the figure, the spike is generated at t=2.1ms, one time step too late, even though the timer seems to hit the threshold D at the right time t=2ms. Only after close inspection, we find that the timer is not exactly identical to the duration at this point in time (see printouts). As expected, the problem disappears if we choose dt as a power of 2, e.g., dt=0.125ms.

timer.zip

[abigailm]

I feel like this discussion is panicking at the wrong disco. Firstly, I absolutely disagree with @diesmann 's statement "This problems of the floating point representation are known to every scientist and it is better to rather expose than to hide them." This is not taught in science degrees and so either you are lucky and someone tells you, or you discover it painfully for yourself, or you stay ignorant. There is no telling which of these states a NESTML user is currently in! Equally I disagree with @tomtetzlaff 's statement: "the solution requires the user to think about potential problems of a numerical implementation of the model -- something any modeler should be aware of." Indeed, a major goal of NESTML is to remove as many as possible of these potential problems in the numerics from the user, who typically doesn't know or doesn't care.

Thus, we don't expect the user to write down how the numerics are to be solved - the user simply writes down the dynamics and NESTML figures out the best way to solve them. The user can indeed have an influence on this at the code-generation stage, or can even write a solver in NESTML - it is not forbidden - but the default behaviour is that NESTML takes care of this in a smart way. Similarly, NESTML should allow the user to state how long the refractory period is, and provide functionality to query whether a neuron is currently refractory, and take care of a graceful exit from refractoriness behind the scenes (in a well documented way of course). Constructs such as integer or (even worse) float decrements should not typically be part of the model description (whilst not being forbidden).

[clinssen]

Proposed three-pronged solution:

1. we introduce a floating point epsilon constant in NESTML, that can be rendered in a default way or its value can be overridden by the code generator options, allowing the user to write, for example, r <= -eps ms;
2. addition of an "include" statement to the language, so different implementations of refractotiness can be easily swapped out (Include statment in NESTML #1121), for example ODE-based vs. (t_last - t) < tau;
3. detailed documentation including numerical examples (covering at least the cases where the refractory period is a multiple of the simulation resolution, and where it isn't).

[tomtetzlaff]

I'm totally in favour of this three-pronged solution. Just two comments:

1. I agree with @abigailm that a nice solution for the refractoriness dynamics would be something such as

if t < t_last_spike - t_ref:
    be_refractory()

This is how we describe models in math. No counters, no ODEs. The downside may be that getting access to the global time t may be expensive, in particular in (asynchronous) neuromorphic systems. Biology indeed uses a local, cell-intrinsic dynamics to implement refractoriness, something we usually describe via ODEs.

2. As mentioned earlier, introducing a floating point epsilon in NESTML means polluting the model description with numerics. But I agree that it would be nice to have this option, so let's be practical. In the documentation, this should be described as a general feature of NESTML. It should become clear that the issues arising from a finite floating point precision are not specific to time variables, such as refractoriness states. There is absolutely no conceptual difference between

if V > VTh:
   do_stuff()

and

if r <= 0 ms:
  do_stuff()

Both r and V (and many other state variables in NESTML/NEST) are floats that are incremented by floats, and in both cases, the threshold may be missed due to the finite floating-point precision. For variables such as the voltage, we are used to accept this imprecision. With introducing the epsilon, the user can now decide to continue this old tradition, or to make it "better":

if V > VTh + eps:
   do_stuff()

or

if r <= 0 ms - eps:
  do_stuff()