From 8075b6a2a592279c8e8fda5de82af6b3a36be43a Mon Sep 17 00:00:00 2001 From: Aditya Singh Date: Wed, 5 Aug 2026 08:07:10 -0700 Subject: [PATCH] Use an exact +- dt window in spike_time_tiling_coefficient run_p built the synchronicity window with np.isclose(t_j, t_i, atol=dt). np.isclose also applies its default relative tolerance rtol=1e-5, so the test is |t_j - t_i| <= dt + 1e-5 * |t_i| rather than |t_j - t_i| <= dt. Spike times are absolute, so the window grew with the distance of the recording from t = 0, and it was asymmetric between the two trains because rtol multiplies only the second argument. run_t computes the T terms with dt exactly, so P and T were being derived from different windows. The effect is a silently wrong coefficient rather than an error. Two spikes 12 ms apart with dt = 5 ms score -0.0005 when the recording starts at t = 0 and +1.0 when the same pair starts at t = 1000 s. For two independent 5 Hz Poisson trains over 600 s with dt = 5 ms, the STTC is +0.027 at t_start = 0 and +0.9998 once the recording carries wall-clock style timestamps, where the correct value is -0.003 throughout. Compare the times directly instead. This matches the documented [-dt, +dt] window, matches what run_t already does, and makes the result invariant under a shift of the whole recording in time. The comparison against C. Cutts' original C implementation in test_sttc_validation_test still passes. The relative tolerance came in with the vectorization in 8bac14c (#564); it looks like an artifact of moving to np.isclose rather than a deliberate tolerance. --- elephant/spike_train_correlation.py | 12 ++++-- elephant/test/test_spike_train_correlation.py | 37 +++++++++++++++++++ 2 files changed, 45 insertions(+), 4 deletions(-) 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