error from amplitude_median metric computation

[ckwalters]

Hi SI team! I am running into an error when I try to compute amplitude median. It has happened both sortings I have tried-- two different kilosort4 sortings of the same session (one split into groups and run within SI, one run outside of SI and extracted). I was actually wanting to do something more simple and get mean amplitude of waveforms in each cluster in microvolts-- is this what the templates extension + the average operator is returning? Thanks!

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[80], [line 2](vscode-notebook-cell:?execution_count=80&line=2)
      [1](vscode-notebook-cell:?execution_count=80&line=1) metric_names=['amplitude_median']
----> [2](vscode-notebook-cell:?execution_count=80&line=2) metrics = si.compute_quality_metrics(sorting_analyzer, metric_names=metric_names)
      [3](vscode-notebook-cell:?execution_count=80&line=3) metrics

File [c:\Users\CKW\anaconda3\envs\kilosort4\lib\site-packages\spikeinterface\core\sortinganalyzer.py:1194](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1194), in AnalyzerExtension.function_factory.<locals>.FuncWrapper.__call__(self, sorting_analyzer, load_if_exists, *args, **kwargs)
   [1191](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1191)         ext = sorting_analyzer.get_extension(self.extension_name)
   [1192](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1192)         return ext
-> [1194](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1194) ext = sorting_analyzer.compute(cls.extension_name, *args, **kwargs)
   [1195](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1195) return ext.get_data()

File [c:\Users\CKW\anaconda3\envs\kilosort4\lib\site-packages\spikeinterface\core\sortinganalyzer.py:763](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:763), in SortingAnalyzer.compute(self, input, save, **kwargs)
    [752](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:752) """
    [753](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:753) Compute one extension or several extensiosn.
    [754](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:754) Internally calls compute_one_extension() or compute_several_extensions() depending on the input type.
   (...)
    [760](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:760)     If the input is a dict then compute several extensions with compute_several_extensions(extensions=input)
    [761](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:761) """
    [762](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:762) if isinstance(input, str):
--> [763](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:763)     return self.compute_one_extension(extension_name=input, save=save, **kwargs)
    [764](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:764) elif isinstance(input, dict):
    [765](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:765)     params_, job_kwargs = split_job_kwargs(kwargs)

File [c:\Users\CKW\anaconda3\envs\kilosort4\lib\site-packages\spikeinterface\core\sortinganalyzer.py:831](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:831), in SortingAnalyzer.compute_one_extension(self, extension_name, save, **kwargs)
    [829](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:829) extension_instance = extension_class(self)
    [830](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:830) extension_instance.set_params(save=save, **params)
--> [831](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:831) extension_instance.run(save=save, **job_kwargs)
    [833](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:833) self.extensions[extension_name] = extension_instance
    [835](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:835) return extension_instance

File [c:\Users\CKW\anaconda3\envs\kilosort4\lib\site-packages\spikeinterface\core\sortinganalyzer.py:1321](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1321), in AnalyzerExtension.run(self, save, **kwargs)
   [1317](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1317) if save and not self.sorting_analyzer.is_read_only():
   [1318](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1318)     # this also reset the folder or zarr group
   [1319](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1319)     self._save_params()
-> [1321](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1321) self._run(**kwargs)
   [1323](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1323) if save and not self.sorting_analyzer.is_read_only():
   [1324](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/core/sortinganalyzer.py:1324)     self._save_data(**kwargs)

File [c:\Users\CKW\anaconda3\envs\kilosort4\lib\site-packages\spikeinterface\qualitymetrics\quality_metric_calculator.py:129](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/quality_metric_calculator.py:129), in ComputeQualityMetrics._run(self, verbose, **job_kwargs)
    [126](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/quality_metric_calculator.py:126) func = _misc_metric_name_to_func[metric_name]
    [128](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/quality_metric_calculator.py:128) params = qm_params[metric_name] if metric_name in qm_params else {}
--> [129](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/quality_metric_calculator.py:129) res = func(self.sorting_analyzer, unit_ids=non_empty_unit_ids, **params)
    [130](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/quality_metric_calculator.py:130) # QM with uninstall dependencies might return None
    [131](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/quality_metric_calculator.py:131) if res is not None:

File [c:\Users\CKW\anaconda3\envs\kilosort4\lib\site-packages\spikeinterface\qualitymetrics\misc_metrics.py:890](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:890), in compute_amplitude_medians(sorting_analyzer, peak_sign, unit_ids)
    [888](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:888) all_amplitude_medians = {}
    [889](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:889) if sorting_analyzer.has_extension("spike_amplitudes") or sorting_analyzer.has_extension("waveforms"):
--> [890](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:890)     amplitudes_by_units = _get_amplitudes_by_units(sorting_analyzer, unit_ids, peak_sign)
    [891](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:891)     for unit_id in unit_ids:
    [892](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:892)         all_amplitude_medians[unit_id] = np.median(amplitudes_by_units[unit_id])

File [c:\Users\CKW\anaconda3\envs\kilosort4\lib\site-packages\spikeinterface\qualitymetrics\misc_metrics.py:753](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:753), in _get_amplitudes_by_units(sorting_analyzer, unit_ids, peak_sign)
    [751](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:751)         unit_index = sorting_analyzer.sorting.id_to_index(unit_id)
    [752](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:752)         spike_mask = spikes["unit_index"] == unit_index
--> [753](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:753)         amplitudes_by_units[unit_id] = all_amplitudes[spike_mask]
    [755](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:755) elif sorting_analyzer.has_extension("waveforms"):
    [756](file:///C:/Users/CKW/anaconda3/envs/kilosort4/lib/site-packages/spikeinterface/qualitymetrics/misc_metrics.py:756)     waveforms_ext = sorting_analyzer.get_extension("waveforms")

IndexError: boolean index did not match indexed array along dimension 0; dimension is 3216122 but corresponding boolean dimension is 3216127

[zm711]

Hi Charlie,

Could you put a bit more of the script you ran to generate the sorting_analyzer above? It seems like your spike_mask is longer than your amplitudes vector but it would probably be helpful to see how everything is being generated to figure out if this is a problem with SI or with how kilosort4 is returning the data. If you sort the same data with KS3 or KS2 or (maybe even one of the other sorters), do you get this error?

[ckwalters]

Yup, sorry about that... I just followed the tutorial in the docs, lmk if anything looks off. I'll try loading in some KS2.5 results now! Thanks :)

recording = se.read_spikeglx(recording_path, stream_id=f'imec{probe}.ap')
sorting = se.read_kilosort(folder_path=os.path.join(recording_path,"kilosort4_20240305"))

sorting_analyzer = create_sorting_analyzer(
    sorting=sorting,
    recording=recording,
    sparse=True, # default
    format="memory", # default
)
job_kwargs = dict(n_jobs=40, chunk_duration='1s', progress_bar=True)
sorting_analyzer.compute("random_spikes", method="uniform", max_spikes_per_unit=500)
sorting_analyzer.compute("waveforms",  ms_before=1.5,ms_after=2.,**job_kwargs)
sorting_analyzer.compute("templates", operators=["average", "median", "std"])
sorting_analyzer.compute("noise_levels")
sorting_analyzer.compute("correlograms")
sorting_analyzer.compute("unit_locations")
sorting_analyzer.compute("spike_amplitudes", **job_kwargs)
sorting_analyzer.compute("template_similarity")

import spikeinterface.qualitymetrics as sqm
amplitude_medians = sqm.compute_amplitude_medians(sorting_analyzer)

[ckwalters]

Ah okay, it's working fine for the same session already sorted by KS2.5!

[zm711]

I think they are still working out the bugs in KS4. We can leave this open just so @alejoe91 or @samuelgarcia can double check to see if they want you to try something else, but since this dataset works for KS2.5 my guess is that this is a KS4 thing.

[ckwalters]

Makes sense, thanks, still early days.

[samuelgarcia]

Hi charlie, Maybe running sorting=remove_excess_spikes(sorting, recording) could helps.

[ckwalters]

@samuelgarcia The dimensions are still off by 2 spikes! It's ok though, it's not an urgent function for me. Thanks!