SpikeInterface / spikeinterface

A Python-based module for creating flexible and robust spike sorting pipelines.
https://spikeinterface.readthedocs.io
MIT License
455 stars 174 forks source link

Error in postprocessing.compute_unit_locations(method = "center_of_mass" #3074

Open jjb314 opened 1 week ago

jjb314 commented 1 week ago

Hello, spikeinterface successfully handles 2D probe information from a MaxWell MEA recording up until I try to find unit locations.

Because of the large data size(~5GB), I computed the SortingAnalyzer and templates with sparse=True.

The array returned by compute_unit_locations is completely flat on the second dimension (all values are 122.5), despite the get_channel_locations() function returning proper geometry.

Additionally, all other methods for unit location return values that don't make sense.

Do you know why I would have this issue?

chrishalcrow commented 1 week ago

Hi @jjb314 , you're right that something seems suspicious here! But I can't think of anything obvious: if you could provide some more info on the questions below, that'd be really useful.

Could you print out the result of your_sorting_analyzer.get_channel_locations() please?

Does your sparsity look reasonable? You can check the sparisty mask by running sa_mask = sorting_analyzer.sparsity.mask. This gives the mask for each unit. You can then look at which channel_locations are non-zero on the mask for the e.g. 3rd unit using sorting_analyzer.get_channel_locations()[sa_mask[3]]. Are these only on the channel locations with y-value 122.5?

Do your templates look sensible? You can plot them using plot_unit_templates(), or whatever you're using to plot later. (More info on widgets: https://spikeinterface.readthedocs.io/en/latest/modules/widgets.html )

What happens when you compute the spike locations, using compute_spike_locations()? Do these numbers reflect the proper geometry?

Thanks!

samuelgarcia commented 1 week ago

The spikeinterface-gui can be usefull to debug with some plots: si.plot_sorting_summary(analyzer, backend="spikeinterface_gui") no ?

jjb314 commented 1 week ago

Thanks for the reply. In short, everything has bad geometry except for the channel locations.

Instead of printing out the arrays of channel locations, unit locations, and spiek lcoations, I made a pyplot of each array for quick visualization. here is a pyplot of the values for channel locations in Red, and computed unit locations in Blue. As you can see, beyond the 2nd axis being flat, the first axis also has wildly different geometry.

Red-Channel_Loc__Blue-Unit_Loc png

While I don't fully understand the sparsity masks aspect of the preprocessing, when I run sorting_analyzer.get_channel_locations()[sa_mask[i]] and print every unit, all have identical values. Here is the output for just the first 2 units, but they're all the same. [[367.5 122.5] [332.5 122.5] [402.5 157.5] [402.5 122.5] [262.5 192.5]] [[367.5 122.5] [332.5 122.5] [402.5 157.5] [402.5 122.5] [262.5 192.5]]

Plotting the templates shows something similar, every unit is on the same electrodes and has 0 amplitude Figure 5

The spike locations are messed up in the same way as the unit locations, and occupy the same general space as the blue dots in the earlier figure.

These kinds of errors don't show up when using phy to visualize the sorter output, so I don't think the issue is with the source recording and sorting files.

I'm at a loss for where I messed up. In the meantime I'm going to recompute it all to see if I mixed up steps.

chrishalcrow commented 1 week ago

Hello, Thanks - interesting! So phy produces nice unit templates and stuff like that? You're right that that suggests it's not a problem with the sorting or recording. We can check if it's something to do with the sparsity (which is trying to select which channels the units are localised on) by creating the sorting analyzer with sparisty switched off. This is done using create_sorting_analyzer(sorting=sorting, recording=recording, sparse=False) This will then use all the probe's channels the templates on. Then you can see if there are some channels with "normal" looking units on them.