diff --git a/elephant/spike_train_correlation.py b/elephant/spike_train_correlation.py index 783860114..43124f528 100644 --- a/elephant/spike_train_correlation.py +++ b/elephant/spike_train_correlation.py @@ -900,10 +900,14 @@ def run_p(spiketrain_j: neo.core.SpikeTrain, """ # Create a boolean array where each element represents whether a spike # in spiketrain_j lies within +- dt of any spike in spiketrain_i. - tiled_spikes_j = np.isclose( - spiketrain_j.times.magnitude[:, np.newaxis], - spiketrain_i.times.magnitude, - atol=dt.item()) + # The comparison is done explicitly rather than with `np.isclose`: + # `np.isclose` also applies its default relative tolerance rtol=1e-5, + # which widens the synchronicity window by 1e-5 * |spike time|. Because + # spike times are absolute, that made the window grow with the distance + # of the recording from t = 0 instead of staying at the requested dt. + tiled_spikes_j = np.abs( + spiketrain_j.times.magnitude[:, np.newaxis] + - spiketrain_i.times.magnitude) <= dt.item() # Determine which spikes in spiketrain_j satisfy the time window # condition. tiled_spike_indices = np.any(tiled_spikes_j, axis=1) diff --git a/elephant/test/test_spike_train_correlation.py b/elephant/test/test_spike_train_correlation.py index d1cad9547..d0e4a5ef6 100644 --- a/elephant/test/test_spike_train_correlation.py +++ b/elephant/test/test_spike_train_correlation.py @@ -839,6 +839,43 @@ def test_sttc_unsorted_spiketimes(self): spiketrain_B3, dt=0.10 * pq.s) self.assertAlmostEqual(sttc_unsorted_E8_B3, sttc_sorted_E8_B3) + def test_sttc_synchronicity_window_is_exactly_dt(self): + # Two spikes 12 ms apart are not synchronous for dt = 5 ms, no matter + # how far the recording is from t = 0. Previously the window was built + # with np.isclose, whose default rtol=1e-5 widened it by + # 1e-5 * |spike time|, so the same pair of spikes was counted as + # synchronous once the recording started late enough. + for t_offset in (0., 1000., 100000.): + spiketrain_i = neo.SpikeTrain( + [10. + t_offset], units='s', + t_start=t_offset * pq.s, t_stop=(20. + t_offset) * pq.s) + spiketrain_j = neo.SpikeTrain( + [10.012 + t_offset], units='s', + t_start=t_offset * pq.s, t_stop=(20. + t_offset) * pq.s) + sttc = sc.sttc(spiketrain_i, spiketrain_j, dt=0.005 * pq.s) + self.assertAlmostEqual(sttc, -0.0005, places=6) + + def test_sttc_invariant_under_time_shift(self): + # The STTC of a pair of spike trains must not change when the whole + # recording is shifted in time. + spiketrain_i = neo.SpikeTrain( + [1.3, 7.56, 15.87, 28.23, 30.9, 34.2, 38.2, 43.2], + units='s', t_stop=50. * pq.s) + spiketrain_j = neo.SpikeTrain( + [1.02, 2.71, 18.82, 28.46, 28.79, 43.6], + units='s', t_stop=50. * pq.s) + target = sc.sttc(spiketrain_i, spiketrain_j, dt=0.005 * pq.s) + + for t_offset in (100., 3600., 86400.): + shifted_i = neo.SpikeTrain( + spiketrain_i.times.magnitude + t_offset, units='s', + t_start=t_offset * pq.s, t_stop=(50. + t_offset) * pq.s) + shifted_j = neo.SpikeTrain( + spiketrain_j.times.magnitude + t_offset, units='s', + t_start=t_offset * pq.s, t_stop=(50. + t_offset) * pq.s) + self.assertAlmostEqual( + sc.sttc(shifted_i, shifted_j, dt=0.005 * pq.s), target) + def test_sttc_validation_test(self): """This test checks the results of elephants implementation of the spike time tiling coefficient against the results of the