Missed time step(s) for noisy / AC current injection in Brian implementation

[appukuttan-shailesh]

It was seen that current injection (NoisyCurrentSource and ACSource) failed to invoke the method '_compute()' for a few time intervals. During these intervals (1 or 2 dts) the current injected continued to be that evaluated for the previous time step.

Splitting the if statement in this line: https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/brian/standardmodels/electrodes.py#L74
and adding a print statement inbetween the two conditions, helps identify this problem. Something like this (continuing from the fix for issue #437):

def _update_current(self):
    if self.running:
        print "s.s.t = ",  simulator.state.t, "\tself.i = ", self.i, "\tself.time[i]*1e3 = ", self.times[self.i]*1e3, "\tdiff = ", (simulator.state.t - (self.times[self.i] * 1e3))
        if simulator.state.t  >= self.times[self.i] * 1e3:
            for cell, idx in zip(self.cell_list, self.indices):
                if not self._is_playable:
                    cell.parent.brian_group.i_inj[idx] = self.amplitudes[self.i] * ampere
                else:
                    cell.parent.brian_group.i_inj[idx] = self._compute(self.times[self.i]) * ampere
            self.i += 1
            if self.i >= len(self.times):
                self.running = False
                if self._is_playable:
                    for cell, idx in zip(self.cell_list, self.indices):
                        cell.parent.brian_group.i_inj[idx] = 0
        else:
            print "_compute() not invoked"

Snippets from some sample outputs for this code are given below:

Note: s.s.t = simulator.state.t

1. start=50.0, stop=125.0, inj_dt = 0.1, setup(timestep = 0.1), run(200.0)

s.s.t =  49.9 	self.i =  0 	self.time[i]*1e3 =  50.0 	diff =  -0.1
_compute() not invoked
s.s.t =  50.0 	self.i =  0 	self.time[i]*1e3 =  50.0 	diff =  0.0
s.s.t =  50.1 	self.i =  1 	self.time[i]*1e3 =  50.1 	diff =  -7.1054273576e-15
_compute() not invoked
s.s.t =  50.2 	self.i =  1 	self.time[i]*1e3 =  50.1 	diff =  0.1
s.s.t =  50.3 	self.i =  2 	self.time[i]*1e3 =  50.2 	diff =  0.1

2. start=40.0, stop=125.0, inj_dt = 0.1, setup(timestep = 0.1), run(200.0)

s.s.t =  39.9 	self.i =  0 	self.time[i]*1e3 =  40.0 	diff =  -0.1
_compute() not invoked
s.s.t =  40.0 	self.i =  0 	self.time[i]*1e3 =  40.0 	diff =  0.0
s.s.t =  40.1 	self.i =  1 	self.time[i]*1e3 =  40.1 	diff =  0.0
s.s.t =  40.2 	self.i =  2 	self.time[i]*1e3 =  40.2 	diff =  -1.42108547152e-14
_compute() not invoked
s.s.t =  40.3 	self.i =  2 	self.time[i]*1e3 =  40.2 	diff =  0.1

3. start=50.0, stop=125.0, inj_dt = 0.01, setup(timestep = 0.01), run(200.0)

s.s.t =  49.99 	self.i =  0 	self.time[i]*1e3 =  50.0 	diff =  -0.00999999999999
_compute() not invoked
s.s.t =  50.0 	self.i =  0 	self.time[i]*1e3 =  50.0 	diff =  0.0
s.s.t =  50.01 	self.i =  1 	self.time[i]*1e3 =  50.01 	diff =  0.0
s.s.t =  50.02 	self.i =  2 	self.time[i]*1e3 =  50.02 	diff =  -7.1054273576e-15
_compute() not invoked
s.s.t =  50.03 	self.i =  2 	self.time[i]*1e3 =  50.02 	diff =  0.00999999999999

4. start=2.0, stop=4.0, inj_dt = 0.01, setup(timestep = 0.01), run(5.0)

s.s.t =  2.1 	self.i =  10 	self.time[i]*1e3 =  2.1 	diff =  0.0
s.s.t =  2.11 	self.i =  11 	self.time[i]*1e3 =  2.11 	diff =  0.0
s.s.t =  2.12 	self.i =  12 	self.time[i]*1e3 =  2.12 	diff =  -4.4408920985e-16
_compute() not invoked
s.s.t =  2.13 	self.i =  12 	self.time[i]*1e3 =  2.12 	diff =  0.01
s.s.t =  2.14 	self.i =  13 	self.time[i]*1e3 =  2.13 	diff =  0.01
s.s.t =  2.15 	self.i =  14 	self.time[i]*1e3 =  2.14 	diff =  0.01

It can be seen that the issue is due to comparison of floating point numbers, and can be overcome using the concept of EPSILON. This would allows us to set an upper bound on the error due to rounding in floating point arithmetic. Will test the same and update here.

[appukuttan-shailesh]

Added a new global method names 'get_time_epsilon()' in file https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/utility/__init__.py
defined as:

def get_time_epsilon():
    return 1e-10    # 100 picoseconds

This is to be imported in the file https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/brian/standardmodels/electrodes.py
by adding:

from pyNN.utility import get_time_epsilon

The test condition for the if statement in _update_current() method here:
https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/brian/standardmodels/electrodes.py#L74
is changed to:

if self.running and abs(simulator.state.t - self.times[self.i] * 1e3) <= get_time_epsilon():

My original intention was to set the value of EPSILON to an even smaller value of 1e-12. But on using that value, the output for case 1 (above) was like this (see here).

Interestingly, the floating point error seemed to progressively increase in magnitude and cross 1e-12 (and subsequently leading to other issues in the above example, as 'self.running' is never set to false). On trying with a larger EPSILON value of 1e-10 things seem to work properly. The output for the earlier examples with this fix are given below:

1. start=50.0, stop=125.0, inj_dt = 0.1, setup(timestep = 0.1), run(200.0)

• see output

2. start=40.0, stop=125.0, inj_dt = 0.1, setup(timestep = 0.1), run(200.0)

• see output

3. start=50.0, stop=125.0, inj_dt = 0.01, setup(timestep = 0.01), run(200.0)

• see output

4. start=2.0, stop=4.0, inj_dt = 0.01, setup(timestep = 0.01), run(5.0)

• see output

But if one takes a closer look at those outputs, the value of difference seems to keep rising and for longer simulations could continue to cross even this value of EPSILON. This process is faster when dt is smaller (e.g. 0.01 ms vs 0.1 ms). So a definite fix for this issue might need to consider the following:

a) Should we go with a fixed value for EPSILON, or have it variable (relative to values being compared). Having a fixed value is more intuitive, but a relative one might alleviate the above issue of gradual rise in the value of difference.

b) It is interesting to note that the instant at which the floating point errors arise is not dependent upon specific timestamps (e.g. 40.1 ms or 50.2 ms), but rather seems to be soon after the initiation of current injection. Could be worth exploring this origin.

[appukuttan-shailesh]

Two factors contributed to the above problem:

1. The accumulation of floating point errors when using numpy.arange() with a floating point time-step.

Example:

import numpy
a = numpy.arange(5, 15, 1)
b = numpy.arange(0.5, 1.5, 0.1) * 10

for i,j in zip(a,b):
    print i, "\t", j, "\t", i-j

prints:

5 	5.0 	0.0
6 	6.0 	0.0
7 	7.0 	0.0
8 	8.0 	8.881784197e-16
9 	9.0 	0.0
10 	10.0 	1.7763568394e-15
11 	11.0 	1.7763568394e-15
12 	12.0 	3.5527136788e-15
13 	13.0 	1.7763568394e-15
14 	14.0 	0.0

numpy.arange() basically adds the step-size to the previous value to determine the next value in the sequence. Addition/subtraction of floats leads to such accumulation of errors, whereas multiplication/division do not. Thus, a better approach would be to construct the array by multiplication as shown below.

import numpy
a = numpy.arange(5, 15, 1)
c = numpy.array([0.5+(i*0.1) for i in range(10)]) * 10
for i,j in zip(a,c):
    print i, "\t", j, "\t", i-j

giving:

5 	5.0 	0.0
6 	6.0 	0.0
7 	7.0 	0.0
8 	8.0 	0.0
9 	9.0 	0.0
10 	10.0 	0.0
11 	11.0 	0.0
12 	12.0 	-1.7763568394e-15
13 	13.0 	0.0
14 	14.0 	0.0

Note: Floating point values (almost) always involve some rounding-off errors (as seen above for 12), and this discussion is only in context of the accumulation of such (small) errors in leading to a larger error.

2. There is also an error involved due to the extent of accuracy offered by double, as seen in the below example. Code:

a = numpy.arange(1, 1000+1, 1)
b = numpy.arange(0.001, 1+0.001, 0.001)

for i in range(len(a)):
    print ("%d\t%0.20f\t%0.20f\t%g" % (i, a[i], b[i]*1000, a[i]-b[i]*1000))

Taking two sample values of the above output corresponding to {a[i]=9, b[i]=0.009} and {a[i]=960, b[i]=0.960}

  9.00000000000000000000	  9.00000000000000177636	-1.77636e-15
960.00000000000000000000	960.00000000000011368684	-1.13687e-13

Quoting from here: http://softwareengineering.stackexchange.com/questions/202843/solutions-for-floating-point-rounding-errors

Due to the nature of double, the larger the whole integer you lose relative accuracy.

[appukuttan-shailesh]

The fix for this matter could be via following changes in https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/brian/standardmodels/electrodes.py

Generate self.times via multiplication (not numpy.arange):

This would help eliminate accumulation of floating-point errors.

method _generate() for ACSource changed from:

def _generate(self):
    self.times = numpy.arange(self.start, self.stop, simulator.state.dt)

to
(Note: this change takes into account issue #442)

def _generate(self):
    temp_num_t = len(numpy.arange(self.start, self.stop + simulator.state.dt * 1e-3, simulator.state.dt * 1e-3))
    self.times = numpy.array([self.start+(i*simulator.state.dt*1e-3) for i in range(temp_num_t)])

and method _generate() for NoisyCurrentSource changed from:

def _generate(self):
    self.times = numpy.arange(self.start, self.stop, simulator.state.dt)

to
(Note: this change takes into account issue #437)

def _generate(self):
    temp_num_t = len(numpy.arange(self.start, self.stop, max(self.dt, simulator.state.dt * 1e-3)))
    self.times = numpy.array([self.start+(i*max(self.dt, simulator.state.dt * 1e-3)) for i in range(temp_num_t)])
    self.times = numpy.append(self.times, self.stop)

Change comparison of timestamps:

(Note: the fix suggested here does not require EPSILON, as was suggested in an earlier post)
The comparison of timestamps on this line: https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/brian/standardmodels/electrodes.py#L74
should be changed from:

if self.running and simulator.state.t >= self.times[self.i] * 1e3:

to

if self.running and abs(simulator.state.t - self.times[self.i] * 1e3) < (simulator.state.dt/2.0):

The above fix suffices as the maximum time resolution we have during the simulation is simulator.state.dt, and hence any timestamp falling within half this value (on either side) could be considered equal.

[apdavison]

Fixed in #467