0 length spike trains cause error in waveform extractors

[colehurwitz]

Hello! One of the researchers at Columbia ran into this issue with waveform extractors that I thought I would post here. If there are 0 length spike trains, the waveform extractor will throw an error. Is this the intended behavior? She can overcome this by getting rid of 0 length spike trains in Phy or by constructing a new sorting extractor (which excludes 0 length spike trains), but seems like this should be handled internally? Curious your thoughts.

we = si.extract_waveforms(recording_f, sorting, folder,
                          load_if_exists=False,
                          overwrite=True)         

IndexError                                Traceback (most recent call last)
Input In [24], in <cell line: 2>()
      1 folder = '/home/xiaji/Downloads/SpikeInterface/Official_Tutorial_SI_0.94_May2022/waveform'
----> 2 we = si.extract_waveforms(recording_f, sorting, folder,
      3                           progress_bar=True,
      4                           n_jobs=4, total_memory="500M", overwrite=True, load_if_exists=False)
      5 print(we)
File ~/miniconda3/envs/si_env/lib/python3.9/site-packages/spikeinterface/core/waveform_extractor.py:798, in extract_waveforms(recording, sorting, folder, load_if_exists, precompute_template, ms_before, ms_after, max_spikes_per_unit, unselect_spike_on_border, overwrite, return_scaled, dtype, use_relative_path, **job_kwargs)
    795 we = WaveformExtractor.create(recording, sorting, folder, use_relative_path=use_relative_path)
    796 we.set_params(ms_before=ms_before, ms_after=ms_after, max_spikes_per_unit=max_spikes_per_unit, dtype=dtype,
    797               return_scaled=return_scaled)
--> 798 we.run_extract_waveforms(**job_kwargs)
    800 if precompute_template is not None:
    801     we.precompute_templates(modes=precompute_template)
File ~/miniconda3/envs/si_env/lib/python3.9/site-packages/spikeinterface/core/waveform_extractor.py:652, in WaveformExtractor.run_extract_waveforms(self, **job_kwargs)
    649 return_scaled = self.return_scaled
    650 unit_ids = self.sorting.unit_ids
--> 652 selected_spikes = self.sample_spikes()
    654 selected_spike_times = {}
    655 for unit_id in self.sorting.unit_ids:
File ~/miniconda3/envs/si_env/lib/python3.9/site-packages/spikeinterface/core/waveform_extractor.py:624, in WaveformExtractor.sample_spikes(self)
    621 nbefore = self.nbefore
    622 nafter = self.nafter
--> 624 selected_spikes = select_random_spikes_uniformly(self.recording, self.sorting,
    625                                                  self._params['max_spikes_per_unit'], nbefore, nafter)
    627 # store in a 2 columns (spike_index, segment_index) in a npy file
    628 for unit_id in self.sorting.unit_ids:
File ~/miniconda3/envs/si_env/lib/python3.9/site-packages/spikeinterface/core/waveform_extractor.py:720, in select_random_spikes_uniformly(recording, sorting, max_spikes_per_unit, nbefore, nafter)
    718 assert nafter is not None
    719 spike_times = sorting.get_unit_spike_train(unit_id=unit_id, segment_index=segment_index)
--> 720 sampled_spike_times = spike_times[inds]
    721 num_samples = recording.get_num_samples(segment_index=segment_index)
    722 mask = (sampled_spike_times >= nbefore) & (sampled_spike_times < (num_samples - nafter))
IndexError: too many indices for array: array is 0-dimensional, but 1 were indexed

[alejoe91]

Hi @colehurwitz

We have a function sorting.remove_empty_units() which returns the new "full" sorting.

I agree that this should be done automatically at waveform extraction. Something like:

if sorting.has_empth_units(): # to be implemented
    self.sorting = sorting.remove_empty_units()

I wouldn't apply the unit removal after spike sorting though. Empty units could mean that a template was found which didn't match any spikes, so it could be informative. Waht do you think?

[colehurwitz]

Yep! Great solution. I agree that the template could be informative, but one could remove these units temporarily for waveform extraction.

[tinaxia2016]

@alejoe91 When I try to use the function remove_empty_units(), I got the following error message:

sorting = se.read_phy('folder_path/mec0/')
sorting = sorting.remove_empty_units()

error message:

File ~/miniconda3/envs/si_env/lib/python3.9/site-packages/spikeinterface/core/basesorting.py:165, in BaseSorting.remove_empty_units(self)
    163 for segment_index in range(self.get_num_segments()):
    164     for unit in self.get_unit_ids():
--> 165         if len(self.get_unit_spike_train(unit, segment_index=segment_index)) > 0:
    166             units_to_keep.append(unit)
    167 units_to_keep = np.unique(units_to_keep)

TypeError: len() of unsized object

[alejoe91]

Hi @tinaxia2016

Sorry for the delay! Could you share the phy folder with me? I'll take a look!