[Bug] IndexError in spike_train_synchrony.py, annotate_synchrofacts

[Moritz-Alexander-Kern]

This Bug was originally discovered by @skrausse , thank you for reporting.

Describe the bug
When trying to use the delete_synchrofacts method from Synchrotool class in spike_train_synchrony.py, in specific cases an IndexError is raised.

To Reproduce

1. Install elephant and neo
2. Run the following minimal example:

# fix synchrotool
import elephant, neo
import quantities as pq
from elephant import spike_train_synchrony


import elephant.spike_train_generation as stg
import numpy as np

np.random.seed(1)  # set random seed 2 for no error

SPP = stg.StationaryPoissonProcess(50 * pq.Hz, t_start=0 * pq.ms,
                                   t_stop=1000*pq.ms)
test_sts = SPP.generate_n_spiketrains(2)

sampling_rate = 30000*pq.Hz
test_obj = spike_train_synchrony.Synchrotool(
    test_sts, sampling_rate=sampling_rate, spread=2)
test_obj.delete_synchrofacts(threshold=2, in_place=True)

Gives the following error messages: (if np.random.seed(1), if the random seed is set to 2 , no error is raised)

Traceback (most recent call last):
    test_obj.delete_synchrofacts(threshold=2, in_place=True)
  elephant/spike_train_synchrony.py", line 325, in delete_synchrofacts
    self.annotate_synchrofacts()
  elephant/spike_train_synchrony.py", line 407, in annotate_synchrofacts
    complexity_per_spike = epoch_complexities[spike_to_epoch_idx]
IndexError: index 91 is out of bounds for axis 0 with size 91

Expected behavior
Index should not be out of bounds

Environment

• OS (e.g., Linux):Ubuntu 20.04.4 LTS (64-bit)
• How you installed elephant (conda, pip, source): pip, source
• Python version: 3.8
• neo python package version: 0.10.2
• elephant python package version: 0.11.0
• quantities: 0.13.0
• numpy: 1.22.3

[Moritz-Alexander-Kern]

Update 20.10.2023: the following is no longer up to date, see latest comments
This bug seems to be related to rounding in quantities and could be traced down to:

elephant/elephant/spike_train_synchrony.py

Lines 399 to 407 in 0df4581

	for idx, st in enumerate(self.input_spiketrains):
	# all indices of spikes that are within the half-open intervals
	# defined by the boundaries
	# note that every second entry in boundaries is an upper boundary
	spike_to_epoch_idx = np.searchsorted(
	right_edges,
	st.times.rescale(self.epoch.times.units).magnitude.flatten())
	complexity_per_spike = epoch_complexities[spike_to_epoch_idx]

More specifically line 406, in certain cases the maximum value here can be larger than the maximum value in right_edges. This could be due to rounding errors.

[Kleinjohann]

Hi!
In this minimal example, the problem is actually that the spiketrains do not have a sampling rate. The Synchrotool assumes that the spiketrains are actually discrete and that all spike times can only be multiples of the sampling rate that is passed as a parameter. We have a utility function for discretising such artificial spiketrains, when we apply that no error is raised:

# fix synchrotool
import elephant, neo
import quantities as pq
from elephant import spike_train_synchrony
from elephant.conversion import discretise_spiketimes


import elephant.spike_train_generation as stg
import numpy as np

np.random.seed(1)  # set random seed 2 for no error

SPP = stg.StationaryPoissonProcess(50 * pq.Hz, t_start=0 * pq.ms,
                                   t_stop=1000*pq.ms)
test_sts = SPP.generate_n_spiketrains(2)

sampling_rate = 30000*pq.Hz
test_sts = discretise_spiketimes(test_sts, sampling_rate)

test_obj = spike_train_synchrony.Synchrotool(
    test_sts, sampling_rate=sampling_rate, spread=2)
test_obj.delete_synchrofacts(threshold=2, in_place=True)

[Kleinjohann]

This error can however also be caused by

1. (t_stop - t_start) % sampling_rate != 0 which means that not the entire period of the spike train will be binned since otherwise the last bin would be incomplete
2. st[-1] == st.t_stop which would cause the last spike to be missed

I see multiple options for addressing this:

• add input checks asserting that (t_stop - t_start) % sampling_rate == 0 and that st[-1] < st.t_stop
• catch the error and raise a more meaningful one
• keep the problematic spikes out of the calculation, then array_annotate them with -1 or some other error value and raise a warning

Either way we also need to make this issue more clear in the documentation of the function.

[Moritz-Alexander-Kern]

The issue could be traced to the following:

Thanks @skrausse who posted on INM-6/elephant:

[...] @jo460464 already found the relevant code in the complexity class, where the t_stop bin is excluded and therefore leads to a failure.

elephant/elephant/statistics.py

Lines 1418 to 1446 in 0354d0e

	def _epoch_no_spread(self):
	"""
	Get an epoch object of the complexity distribution with `spread` = 0
	"""
	left_edges = self.time_histogram.times
	durations = self.bin_size * np.ones(self.time_histogram.shape)
	if self.sampling_rate:
	# ensure that spikes are not on the bin edges
	bin_shift = .5 / self.sampling_rate
	left_edges -= bin_shift
	# Ensure that an epoch does not start before the minimum t_start.
	# Note: all spike trains share the same t_start and t_stop.
	if left_edges[0] < self.t_start:
	left_edges[0] = self.t_start
	durations[0] -= bin_shift
	else:
	warnings.warn('No sampling rate specified. '
	'Note that using the complexity epoch to get '
	'precise spike times can lead to rounding errors.')
	complexity = self.time_histogram.magnitude.flatten()
	complexity = complexity.astype(np.uint16)
	epoch = neo.Epoch(left_edges,
	durations=durations,
	array_annotations={'complexity': complexity})
	return epoch

Shifting the bins by half a sampling period leaves events at t_stop to be out of bounds.

[skrausse]

Discussion at elephant meeting (23.10.23):

• binning issue traced back synchrotool -> complexity -> time_histogram -> binned_spiketrains
• if provided with sampling rate / period: Shifting of spiketrains and adding an epoch would solve the problem
  • Issue: binned spiketrain relies on shared t_start and t_stop between spiketrains and provided values to not make assumptions on sampling rate
• proposed fix: condition on binned spiketrain: if sampling rate provided and shift parameter allowed, then shift spiketrains; else warn for spikes being removed if that coincides with t_stop

[Moritz-Alexander-Kern]

We also discussed the suggested changes on branch: https://github.com/jo460464/elephant/tree/enh/iss101