Merge b863024 into d036012

pyNN/brian/standardmodels/electrodes.py

@@ -95,30 +95,37 @@ def inject_into(self, cell_list):
             if not cell.celltype.injectable:
                 raise TypeError("Can't inject current into a spike source.")
         self.cell_list.extend(cell_list)
+        self.prev_amp_dict = {}
         for cell in cell_list:
-            self.indices.extend([cell.parent.id_to_index(cell)])
+            cell_idx = cell.parent.id_to_index(cell)
+            self.indices.extend([cell_idx])
+            self.prev_amp_dict[cell_idx] = 0.0
 
     def _update_current(self):
         # check if current timestamp is within dt/2 of target time; Brian uses seconds as unit of time
         if self.running and abs(simulator.state.t - self.times[self.i] * 1e3) < (simulator.state.dt/2.0):
             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
+                    cell.parent.brian_group.i_inj[idx] += self.amplitudes[self.i] - self.prev_amp_dict[idx] #* ampere
+                    self.prev_amp_dict[idx] = self.amplitudes[self.i]
                 else:
-                    cell.parent.brian_group.i_inj[idx] = self._compute(self.times[self.i]) * ampere
+                    amp_val = self._compute(self.times[self.i])
+                    self.amplitudes = numpy.append(self.amplitudes, amp_val)
+                    cell.parent.brian_group.i_inj[idx] += amp_val - self.prev_amp_dict[idx]
+                    self.prev_amp_dict[idx] = amp_val #* ampere
             self.i += 1
             if self.i >= len(self.times):
                 self.running = False
                 if self._is_playable:
                     # ensure that currents are set to 0 after t_stop
                     for cell, idx in zip(self.cell_list, self.indices):
-                        cell.parent.brian_group.i_inj[idx] = 0
+                        cell.parent.brian_group.i_inj[idx] -= self.prev_amp_dict[idx]
 
-    def record(self):
+    def _record_old(self):
         self.i_state_monitor = brian.StateMonitor(self.cell_list[0].parent.brian_group[self.indices[0]], 'i_inj', record=0, when='start')
         simulator.state.network.add(self.i_state_monitor)
 
-    def _get_data(self):
+    def _get_data_old(self):
         # code based on brian/recording.py:_get_all_signals()
         # because we use `when='start'`, we need to add the
         # value at the end of the final time step.
@@ -129,6 +136,32 @@ def _get_data(self):
         i_all_values = numpy.append(device._values, current_i_value)
         return (t_all_values / ms, i_all_values / nA)
 
+    def record(self):
+        pass
+
+    def _get_data(self):
+        def find_nearest(array, value):
+            array = numpy.asarray(array)
+            return (numpy.abs(array - value)).argmin()
+
+        len_t = int(round((simulator.state.t * 1e-3) / (simulator.state.dt * 1e-3))) + 1
+        times = numpy.array([(i * simulator.state.dt * 1e-3) for i in range(len_t)])
+        amps = numpy.array([0.0] * len_t)
+
+        for idx, [t1, t2] in enumerate(zip(self.times, self.times[1:])):
+            if t2 < simulator.state.t * 1e-3:
+                idx1 = find_nearest(times, t1)
+                idx2 = find_nearest(times, t2)
+                amps[idx1:idx2] = [self.amplitudes[idx]] * len(amps[idx1:idx2])
+                if idx == len(self.times)-2:
+                    if not self._is_playable and not self._is_computed:
+                        amps[idx2:] = [self.amplitudes[idx+1]] * len(amps[idx2:])
+            else:
+                if t1 < simulator.state.t * 1e-3:
+                    idx1 = find_nearest(times, t1)
+                    amps[idx1:] = [self.amplitudes[idx]] * len(amps[idx1:])
+                break
+        return (times / ms, amps / nA)
 
 class StepCurrentSource(BrianCurrentSource, electrodes.StepCurrentSource):
     __doc__ = electrodes.StepCurrentSource.__doc__
@@ -161,8 +194,9 @@ def __init__(self, **parameters):
 
     def _generate(self):
         # Note: Brian uses seconds as unit of time
-        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)])
+        temp_num_t = int(round(((self.stop + simulator.state.dt * 1e-3) - self.start) / (simulator.state.dt * 1e-3)))
+        self.times = numpy.array([self.start + (i * simulator.state.dt * 1e-3) for i in range(temp_num_t)])
+        self.amplitudes = numpy.zeros(0)
 
     def _compute(self, time):
         # Note: Brian uses seconds as unit of time; frequency is specified in Hz; thus no conversion required
@@ -217,9 +251,10 @@ def __init__(self, **parameters):
         self._generate()
 
     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)])
+        temp_num_t = int(round((self.stop - self.start) / 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)
+        self.amplitudes = numpy.zeros(0)
 
     def _compute(self, time):
         return self.mean + self.stdev * numpy.random.randn()

pyNN/neuron/standardmodels/electrodes.py

@@ -210,8 +210,9 @@ def __init__(self, **parameters):
     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)
+        temp_num_t = int(round(((self.stop + simulator.state.dt) - self.start) / simulator.state.dt))
+        self.times = numpy.array([self.start + (i * simulator.state.dt) for i in range(temp_num_t)])
+        tmp = numpy.array([i * simulator.state.dt for i in range(temp_num_t)])
         self.amplitudes = self.offset + self.amplitude * numpy.sin(tmp * 2 * numpy.pi * self.frequency / 1000. + 2 * numpy.pi * self.phase / 360)
         self.amplitudes[-1] = 0.0
 
@@ -238,7 +239,8 @@ def __init__(self, **parameters):
     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))
+        temp_num_t = int(round((self.stop - self.start) / max(self.dt, simulator.state.dt)))
+        self.times = numpy.array([self.start + (i * max(self.dt, simulator.state.dt)) for i in range(temp_num_t)])
         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

test/system/scenarios/test_electrodes.py

@@ -369,8 +369,8 @@ def issue487(sim):
     assert_true (numpy.isclose(v_step_2_arr[0:int(step_2.times[0]/dt)], v_rest).all())
 
 
-@register(exclude=["brian"])  # todo: fix for Brian
-def issue_465_474(sim):
+@register()
+def issue_465_474_630(sim):
     """
     Checks the current traces recorded for each of the four types of
     electrodes in pyNN, and verifies that:
@@ -441,12 +441,18 @@ def issue_465_474(sim):
     assert_true (i_noise.t_start == 0.0 * pq.ms and numpy.isclose(float(i_noise.times[-1]), simtime))
     assert_true (i_step.t_start == 0.0 * pq.ms and numpy.isclose(float(i_step.times[-1]), simtime))
 
-    # test to check current changes at the expected time instant
+    # test to check current changes at start time instant
     assert_true (i_ac[(int(start / sim_dt)) - 1, 0] == 0 * pq.nA and i_ac[int(start / sim_dt), 0] != 0 * pq.nA)
     assert_true (i_dc[int(start / sim_dt) - 1, 0] == 0 * pq.nA and i_dc[int(start / sim_dt), 0] != 0 * pq.nA)
     assert_true (i_noise[int(start / sim_dt) - 1, 0] == 0 * pq.nA and i_noise[int(start / sim_dt), 0] != 0 * pq.nA)
     assert_true (i_step[int(start / sim_dt) - 1, 0] == 0 * pq.nA and i_step[int(start / sim_dt), 0] != 0 * pq.nA)
 
+    # test to check current changes appropriately at stop time instant - issue #630
+    assert_true (i_ac[(int(stop / sim_dt)) - 1, 0] != 0.0 * pq.nA and i_ac[(int(stop / sim_dt)), 0] == 0.0 * pq.nA )
+    assert_true (i_dc[(int(stop / sim_dt)) - 1, 0] != 0.0 * pq.nA and i_ac[(int(stop / sim_dt)), 0] == 0.0 * pq.nA )
+    assert_true (i_noise[(int(stop / sim_dt)) - 1, 0] != 0.0 * pq.nA and i_ac[(int(stop / sim_dt)), 0] == 0.0 * pq.nA )
+    assert_true (i_step[(int(stop / sim_dt)) - 1, 0] != 0.2 * pq.nA and i_ac[(int(stop / sim_dt)), 0] == 0.0 * pq.nA )
+
     # test to check vm changes at the time step following current initiation
     assert_true (numpy.isclose(float(v_ac[int(start / sim_dt), 0].item()), v_rest) and v_ac[int(start / sim_dt) + 1] != v_rest * pq.mV)
     assert_true (numpy.isclose(float(v_dc[int(start / sim_dt), 0].item()), v_rest) and v_dc[int(start / sim_dt) + 1] != v_rest * pq.mV)
@@ -601,6 +607,44 @@ def get_len(data):
             assert_true (abs(step.times[1]-0.86) < 1e-9)
 
 
+@register()
+def issue631(sim):
+    """
+    Test to ensure that recording of multiple electrode currents do not
+    interfere with one another.
+    """
+    sim_dt = 0.1
+    sim.setup(timestep=sim_dt, min_delay=sim_dt)
+
+    cells = sim.Population(1, sim.IF_curr_exp(v_rest = -65.0, v_thresh=-55.0, tau_refrac=5.0))#, i_offset=-1.0*amp))
+    dc_source = sim.DCSource(amplitude=0.5, start=25, stop=50)
+    ac_source = sim.ACSource(start=75, stop=125, amplitude=0.5, offset=0.25, frequency=100.0, phase=0.0)
+    noisy_source = sim.NoisyCurrentSource(mean=0.5, stdev=0.05, start=150, stop=175, dt=1.0)
+    step_source = sim.StepCurrentSource(times=[200, 225, 250], amplitudes=[0.4, 0.6, 0.2])
+
+    cells[0].inject(dc_source)
+    cells[0].inject(ac_source)
+    cells[0].inject(noisy_source)
+    cells[0].inject(step_source)
+
+    dc_source.record()
+    ac_source.record()
+    noisy_source.record()
+    step_source.record()
+
+    sim.run(275.0)
+
+    i_dc = dc_source.get_data()
+    i_ac = ac_source.get_data()
+    i_noisy = noisy_source.get_data()
+    i_step = step_source.get_data()
+
+    assert_true (numpy.all(i_dc.magnitude[:int(25.0 / sim_dt) - 1:] == 0) and numpy.all(i_dc.magnitude[int(50.0 / sim_dt):] == 0))
+    assert_true (numpy.all(i_ac.magnitude[:int(75.0 / sim_dt) - 1:] == 0) and numpy.all(i_ac.magnitude[int(125.0 / sim_dt):] == 0))
+    assert_true (numpy.all(i_noisy.magnitude[:int(150.0 / sim_dt) - 1:] == 0) and numpy.all(i_noisy.magnitude[int(175.0 / sim_dt):] == 0))
+    assert_true (numpy.all(i_step.magnitude[:int(200.0 / sim_dt) - 1:] == 0))
+
+
 if __name__ == '__main__':
     from pyNN.utility import get_simulator
     sim, args = get_simulator()
@@ -614,6 +658,7 @@ def get_len(data):
     issue451(sim)
     issue483(sim)
     issue487(sim)
-    issue_465_474(sim)
+    issue_465_474_630(sim)
     issue497(sim)
     issue512(sim)
+    issue631(sim)

test/system/test_brian.py

@@ -78,3 +78,23 @@ def test_tsodyks_markram_synapse():
     sim.run(100.0)
     tau_psc = prj._brian_synapses[0][0].tau_syn.data * 1e3  # s --> ms
     assert_array_equal(tau_psc, numpy.arange(0.2, 0.7, 0.1))
+
+
+def test_issue_634():
+    if not have_brian:
+        raise SkipTest
+    sim = pyNN.brian
+    sim.setup()
+    cells = sim.Population(1, sim.IF_curr_exp(v_thresh=-55.0, tau_refrac=5.0, v_rest=-60.0))
+    dc1_source = sim.DCSource(amplitude=0.5, start=25, stop=50)
+    dc2_source = sim.DCSource(amplitude=-0.5, start=40, stop=80)
+    cells[0].inject(dc1_source)
+    cells[0].inject(dc2_source)
+    dc1_source.record()
+    dc2_source.record()
+    dc1_source._record_old()
+    sim.run(100.0)
+    i_dc1 = dc1_source.get_data()
+    i_dc2 = dc2_source.get_data()
+    i_dc_total_times, i_dc_total = dc1_source._get_data_old()
+    assert_array_equal(numpy.squeeze(i_dc1.magnitude+i_dc2.magnitude), i_dc_total)