Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

[Bug] IndexError in spike_train_synchrony.py, annotate_synchrofacts #493

Open
Moritz-Alexander-Kern opened this issue Jun 8, 2022 · 6 comments · May be fixed by #637
Open

[Bug] IndexError in spike_train_synchrony.py, annotate_synchrofacts #493

Moritz-Alexander-Kern opened this issue Jun 8, 2022 · 6 comments · May be fixed by #637
Labels
bug Indicates an unexpected problem or unintended behavior

Comments

@Moritz-Alexander-Kern
Copy link
Member

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 Moritz-Alexander-Kern added the bug Indicates an unexpected problem or unintended behavior label Jun 8, 2022
@Moritz-Alexander-Kern
Copy link
Member Author

Moritz-Alexander-Kern commented Jun 8, 2022

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:

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
Copy link
Contributor

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
Copy link
Contributor

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
Copy link
Member Author

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.

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
Copy link

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
Copy link
Member Author

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

# for free to join this conversation on GitHub. Already have an account? # to comment
Labels
bug Indicates an unexpected problem or unintended behavior
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants