Check implementations of NoisyCurrentSource

[pgleeson]

This file runs differently on Neuron, Nest, Brian (NoisyCurrentSource with dt = 1ms, 5ms, 10ms).

NEST (probably correct)

NEURON:

Brian:

@pedroernesto will look into this...

[apdavison]

This line is incorrect in the pyNN.neuron implementation:
https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/neuron/standardmodels/electrodes.py#L190
It should be something like:

self.times = numpy.arange(self.start, self.stop + simulator.state.dt, max(self.dt, simulator.state.dt))

not

self.times = numpy.arange(self.start, self.stop + simulator.state.dt, simulator.state.dt)

and then similarly in line 191

I'm not sure what went wrong with the Brian implementation.

[appukuttan-shailesh]

In the above figure, there are two main problems with the Brian output:

1. current injection not stopping (at stop = 450 ms)
2. the lack of noise

Fix for problem 1:

Th method '_update_current()' here:
https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/brian/standardmodels/electrodes.py#L73
needs to be updated to set the injected current to be zero after t_stop. This problem also exists for ACSource (see Issue #422).

The fix is to change '_update_current()'
from:

    def _update_current(self):
        if self.running and 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

to:

    def _update_current(self):
        if self.running and 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
        else:
            if self._is_playable and self.i == len(self.times) and simulator.state.t >= (self.stop * 1e3):
                for cell, idx in zip(self.cell_list, self.indices):
                    cell.parent.brian_group.i_inj[idx] = 0
                    self.i += 1

The (new) if statement has three clauses:
i) self._is_playable (true)
ensuring it applies only to NoisyCurrentSource and ACSource
ii) self.i == len(self.times)
ensures the condition evaluates to true only once
iii) simulator.state.t >= (self.stop * 1e3)
ensures evaluation only at t_stop

Will check if there is a more optimal way of doing this.

Fix for problem 2:

Brian requires time to be specified with units of seconds. Most of the required conversion is done, except for simulator.state.dt which is still specified in ms. Thus the following lines to be changed:

i) https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/brian/standardmodels/electrodes.py#L163
from:

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

to (in accordance with the fix for NEURON suggested above by @apdavison):

self.times = numpy.arange(self.start, self.stop, max(self.dt, simulator.state.dt * 1e-3))

ii) https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/brian/standardmodels/electrodes.py#L166
from:

return self.mean + (self.stdev * self.dt) * numpy.random.randn()

to (in accordance with the fix for NEURON suggested above by @apdavison):

return self.mean + self.stdev * numpy.random.randn()

This seems to be similar to the way NoisyCurrentSource is implemented in NEST. To ensure consistency across simulators, this could be adopted for now.

Regarding the NEURON fix:

Just adding a bit to the fix for NEURON suggested above by @apdavison:
Changing the method _generate() in line https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/neuron/standardmodels/electrodes.py#L187
from:

def _generate(self):
        ## Not efficient at all... Is there a way to have those vectors computed on the fly ?
        ## Otherwise should have a buffer mechanism
        self.times = numpy.arange(self.start, self.stop + simulator.state.dt, simulator.state.dt)
        tmp = numpy.arange(0, self.stop - self.start, simulator.state.dt)
        self.amplitudes = self.mean + (self.stdev * self.dt) * numpy.random.randn(len(tmp))
        self.amplitudes[-1] = 0.0

to:

def _generate(self):
        ## Not efficient at all... Is there a way to have those vectors computed on the fly ?
        ## Otherwise should have a buffer mechanism
        self.times = numpy.arange(self.start, self.stop, max(self.dt, simulator.state.dt))
        self.times = numpy.append(self.times, self.stop)        
        self.amplitudes = self.mean + self.stdev * numpy.random.randn(len(self.times))
        self.amplitudes[-1] = 0.0

This updates the way the amplitudes are calculated (in accordance with Brian fix and NEST). Also, explicitly adds t_stop to the list of times (with corresponding amplitude zero). This ensures current injection in the last time interval.

Output with above fixes

NEST:

Brian:

NEURON:

[appukuttan-shailesh]

A neater fix for problem 1 (from the previous post) might be to add t_stop to the array of times, and have a new check within the existing if statement.

This would involve modifying the method _generate() in line https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/brian/standardmodels/electrodes.py#L162
from:

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

to

def _generate(self):
        self.times = numpy.arange(self.start, self.stop, max(self.dt, simulator.state.dt * 1e-3))
        self.times = numpy.append(self.times, self.stop)

And the method _update_current() in line https://github.com/NeuralEnsemble/PyNN/blob/master/pyNN/brian/standardmodels/electrodes.py#L73
from:

def _update_current(self):
        if self.running and 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

to

def _update_current(self):
        if self.running and 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

This would potentially reduce the number of computations involved than in the earlier fix (for problem 1)

[apdavison]

Fixed in #467