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