Compare two sorters

This example show how to compare the result of two sorters.

Import

import numpy as np
import matplotlib.pyplot as plt

import spikeinterface as si
import spikeinterface.extractors as se
import spikeinterface.sorters as ss
import spikeinterface.comparison as sc
import spikeinterface.widgets as sw

First, let’s download a simulated dataset

from the repo ‘https://gin.g-node.org/NeuralEnsemble/ephy_testing_data’

local_path = si.download_dataset(remote_path='mearec/mearec_test_10s.h5')
recording, sorting = se.read_mearec(local_path)
print(recording)
print(sorting)

Out:

MEArecRecordingExtractor: 32 channels - 1 segments - 32.0kHz - 10.000s
  file_path: /home/docs/spikeinterface_datasets/ephy_testing_data/mearec/mearec_test_10s.h5
MEArecSortingExtractor: 10 units - 1 segments - 32.0kHz
  file_path: /home/docs/spikeinterface_datasets/ephy_testing_data/mearec/mearec_test_10s.h5

Then run two spike sorters and compare their output.

sorting_HS = ss.run_herdingspikes(recording)
sorting_TDC = ss.run_tridesclous(recording)

Traceback (most recent call last):
  File "/home/docs/checkouts/readthedocs.org/user_builds/spikeinterface/checkouts/0.91.0/examples/modules/comparison/plot_1_compare_two_sorters.py", line 37, in <module>
    sorting_TDC = ss.run_tridesclous(recording)
  File "/home/docs/checkouts/readthedocs.org/user_builds/spikeinterface/conda/0.91.0/lib/python3.8/site-packages/spikeinterface/sorters/runsorter.py", line 320, in run_tridesclous
    return run_sorter('tridesclous', *args, **kwargs)
  File "/home/docs/checkouts/readthedocs.org/user_builds/spikeinterface/conda/0.91.0/lib/python3.8/site-packages/spikeinterface/sorters/runsorter.py", line 58, in run_sorter
    sorting = run_sorter_local(sorter_name, recording, output_folder=output_folder,
  File "/home/docs/checkouts/readthedocs.org/user_builds/spikeinterface/conda/0.91.0/lib/python3.8/site-packages/spikeinterface/sorters/runsorter.py", line 83, in run_sorter_local
    SorterClass.run_from_folder(output_folder, raise_error, verbose)
  File "/home/docs/checkouts/readthedocs.org/user_builds/spikeinterface/conda/0.91.0/lib/python3.8/site-packages/spikeinterface/sorters/basesorter.py", line 234, in run_from_folder
    raise SpikeSortingError(
spikeinterface.sorters.utils.misc.SpikeSortingError: Spike sorting failed. You can inspect the runtime trace in spikeinterface_log.json

The compare_two_sorters function allows us to compare the spike  sorting output. It returns a SortingComparison object, with methods  to inspect the comparison output easily. The comparison matches the  units by comparing the agreement between unit spike trains.

 Let’s see how to inspect and access this matching.

cmp_HS_TDC = sc.compare_two_sorters(sorting1=sorting_HS, sorting2=sorting_TDC,
                                               sorting1_name='HS', sorting2_name='TDC')

We can check the agreement matrix to inspect the matching.

sw.plot_agreement_matrix(cmp_HS_TDC)

Some useful internal dataframes help to check the match and count   like match_event_count or agreement_scores

print(cmp_HS_TDC.match_event_count)
print(cmp_HS_TDC.agreement_scores)

In order to check which units were matched, the get_matching  methods can be used. If units are not matched they are listed as -1.

sc_to_tdc, tdc_to_sc = cmp_HS_TDC.get_matching()

print('matching HS to TDC')
print(sc_to_tdc)
print('matching TDC to HS')
print(tdc_to_sc)

The :code:get_unit_spike_train` returns the mapped spike train. We can use  it to check the spike times.

matched_ids = sc_to_tdc[sc_to_tdc != -1]

unit_id_HS = matched_ids.index[0]
unit_id_TDC = matched_ids[unit_id_HS]



# check that matched spike trains correspond
st1 = sorting_HS.get_unit_spike_train(unit_id_HS)
st2 = sorting_TDC.get_unit_spike_train(unit_id_TDC)
fig, ax = plt.subplots()
ax.plot(st1, np.zeros(st1.size), '|')
ax.plot(st2, np.ones(st2.size), '|')

plt.show()