<jupyter_start><jupyter_text>Sorter comparison with ground truth recordings@author : <NAME>  Spampinato from Institut de la Vision, Paris, have publish data from mice retina recorded invitro with from dense array.The data contain one ground truth to benchmark spike sorting tools. The ground truth is recorded with juxta celullar pipet.Here the official publication of this open dataset:https://zenodo.org/record/1205233.W9mq1HWLTIFThis datasets was used by Pierre Yger publish spyking circus:https://elifesciences.org/articles/34518In the original dataset, 19 files are provided.here we volontary keep only 9 files.See carrefully the notebook : **prepare_and_check_ground_truth.py.ipynb**  How to reproduce this notebook  1. to create a folder **basedir** somwhere.  2. download all files in zenodo link  3. put them all in a subfolder **basedir/original_files** (20160415_patch2.tar.gz, ...)  4. and then execute cells by cell this notebook<jupyter_code>import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

import spikeextractors as se
import spiketoolkit as st
from spiketoolkit.study import GroundTruthStudy

%matplotlib notebook<jupyter_output>23:00:15 [I] klustakwik KlustaKwik2 version 0.2.6
/home/samuel/.virtualenvs/py36/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject
  return f(*args, **kwds)
/home/samuel/.virtualenvs/py36/lib/python3.6/site-packages/sklearn/externals/six.py:31: DeprecationWarning: The module is deprecated in version 0.21 and will be removed in version 0.23 since we've dropped support for Python 2.7. Please rely on the official version of six (https://pypi.org/project/six/).
  "(https://pypi.org/project/six/).", DeprecationWarning)
/home/samuel/.virtualenvs/py36/lib/python3.6/site-packages/sklearn/externals/joblib/__init__.py:15: DeprecationWarning: sklearn.externals.joblib is deprecated in 0.21 and will be removed in 0.23. Please import this functionality directly from joblib, which can be installed with: pip install joblib. If this warning is raised when loading pickled models, you[...]<jupyter_text>global variables<jupyter_code># my working path
basedir = '/media/samuel/SamCNRS/DataSpikeSorting/pierre/zenodo/'
# basedir = '/home/samuel/DataSpikeSorting/Pierre/zenodo/'
# basedir = '/mnt/data/sam/DataSpikeSorting/pierre_zenodo/'
# basedir = '/media/samuel/dataspikesorting/DataSpikeSortingHD2/Pierre/zenodo/'

# input file
recording_folder = basedir + 'original_files/'

# ground truth information
ground_truth_folder = basedir + 'ground_truth/'

# where output will be
study_folder = basedir + 'study_gt252/'

# sorter list
sorter_list = ['tridesclous']

# selected recordings (8/19)
rec_names = [
    '20160415_patch2',
    '20170803_patch1',
    '20160426_patch3', 
    '20170725_patch1',
    '20170621_patch1',
    '20160426_patch2', 
    '20170728_patch2',
    '20170713_patch1',
]<jupyter_output><empty_output><jupyter_text>Step 1 : detect ground truth on the juxta cellularThis external script:    * unzip all    * detect juxta peak    * make figures to check everything<jupyter_code>%run detect_ground_truth_spike_on_juxta.py<jupyter_output><empty_output><jupyter_text>Step 2 : check groud truth carcteristicsThis is done in a separate script to keep this doc simple.<jupyter_code>gt_info = pd.read_excel(ground_truth_folder+'gt_info.xlsx', index_col=0)
gt_info.loc[rec_names, :]<jupyter_output><empty_output><jupyter_text>Step 3: setup studyIn this step:  * we create a dict of (recordings+sorting)  * and call the function ```setup_comparison_study(study_folder, gt_dict)```Internanly spiketoolkit copy all recordings and ground truth sorting into a n organised folder.Important note :  * the file have 256 channels but only 252 are usefull.    The PRB file contain all channels needed so we need to explicit **grouping_property='group'**    to be sure to only take in account the channel in the unique group.  * This step have to run only once and take some times because of copy.<jupyter_code>gt_dict = {}
for rec_name in rec_names:

    # find raw file
    dirname = recording_folder + rec_name + '/'
    for f in os.listdir(dirname):
        if f.endswith('.raw') and not f.endswith('juxta.raw'):
            mea_filename = dirname + f

    # raw files have an internal offset that depend on the channel count
    # a simple built header can be parsed to get it
    with open(mea_filename.replace('.raw', '.txt'), mode='r') as f:
        offset = int(re.findall('padding = (\d+)', f.read())[0])

    # recording
    rec = se.BinDatRecordingExtractor(mea_filename, 20000., 256, 'uint16', offset=offset, frames_first=True)

    # this reduce channel count to 252
    rec = se.load_probe_file(rec, basedir + 'mea_256.prb')

    # gt sorting
    gt_indexes = np.fromfile(ground_truth_folder + rec_name + '/juxta_peak_indexes.raw', dtype='int64')
    sorting_gt = se.NumpySortingExtractor()
    sorting_gt.set_times_labels(gt_indexes, np.zeros(gt_indexes.size, dtype='int64'))
    sorting_gt.set_sampling_frequency(20000.0)

    gt_dict[rec_name] = (rec, sorting_gt)


study = GroundTruthStudy.setup(study_folder, gt_dict)<jupyter_output><empty_output><jupyter_text>Step 4: run all sortersspiketoolkit have convinient function for that loop over all recordings and all sorter, compute and save the output sorting into a subfolder.Note:  * this can be run several times with other sorter  * the ```engine='loop'``` prevent to recompute something alreadu computed<jupyter_code>study.run_sorters(sorter_list, mode='keep', engine='loop')<jupyter_output><empty_output><jupyter_text>Step 4 : collect resultsThis collect in the folder all output sorting and compare then the ground truth.Internally it is just a loop and a use of the class ```GroundTruthComparison```.  * comparisons: contains all GroundTruthComparison  * dataframes: contains some synthetic dataframes that summary results<jupyter_code># avoid re run sorter if NB is reopen
study = GroundTruthStudy(study_folder)
%%capture --no-display
study.run_comparisons(exhaustive_gt=False)

comparisons = study.comparisons
dataframes = study.aggregate_dataframes()

print(comparisons.keys())<jupyter_output><empty_output><jupyter_text>Compare run speedThe speed here are totally wrong because sorters have been launch on several computers.This is just to show that it is possible.<jupyter_code>run_times = dataframes['run_times']
run_times
run_times = dataframes['run_times'].reset_index()
fig, ax = plt.subplots()
sns.barplot(data=run_times, x='rec_name', y='run_time', hue='sorter_name', ax=ax)
ax.set_title('Run times')<jupyter_output><empty_output><jupyter_text>Compare global performancesIn our case **dataframes['perf_pooled_with_sum']**,  **dataframes['perf_pooled_with_average']** and **dataframes['perf_by_spiketrain']** because there is only one unique ground truth cells per recordings we will call it perfs.<jupyter_code>perfs = dataframes['perf_by_unit']
perfs.reset_index().sort_values(by=['rec_name', 'sorter_name'])
perfs = dataframes['perf_by_unit']
fig, ax = plt.subplots()
g = sns.swarmplot(data=perfs, x='sorter_name', y='accuracy', hue='sorter_name', ax=ax)
g.legend(loc='center right', bbox_to_anchor=(1.5, 0.5))
fig.suptitle('accuracy')
perfs = dataframes['perf_by_unit']
fig, ax = plt.subplots()
g = sns.swarmplot(data=perfs, x='sorter_name', y='recall', hue='sorter_name', ax=ax)
g.legend(loc='center right', bbox_to_anchor=(1.5, 0.5))
fig.suptitle('recall')
perfs = dataframes['perf_by_unit']
fig, ax = plt.subplots()
g = sns.swarmplot(data=perfs, x='sorter_name', y='precision', hue='sorter_name', ax=ax)
g.legend(loc='center right', bbox_to_anchor=(1.5, 0.5))
fig.suptitle('precision')
perfs = dataframes['perf_by_unit']
fig, ax = plt.subplots()
g = sns.swarmplot(data=perfs, x='sorter_name', y='false_discovery_rate', hue='sorter_name', ax=ax)
g.legend(loc='center right', bbox_to_anchor=(1.5, 0.5))
fig.suptitle('false_discovery_rate')
perfs = dataframes['perf_by_unit']
fig, ax = plt.subplots()
sns.barplot(data=perfs, x='rec_name', y='accuracy', hue='sorter_name', ax=ax)
ax.set_title('accuracy')
ax.set_ylim(0, 1.)
ax.legend(loc='center right', bbox_to_anchor=(1.5, 0.5))
r = plt.setp(ax.xaxis.get_majorticklabels(), rotation=55, ha='right')
fig, ax = plt.subplots()
sns.barplot(data=perfs, x='rec_name', y='recall', hue='sorter_name', ax=ax)
ax.set_title('recall')
ax.set_ylim(0, 1.)
ax.legend(loc='center right', bbox_to_anchor=(1.5, 0.5))
r = plt.setp(ax.xaxis.get_majorticklabels(), rotation=55, ha='right')
fig, ax = plt.subplots()
sns.barplot(data=perfs, x='rec_name', y='precision', hue='sorter_name', ax=ax)
ax.set_title('precision')
ax.set_ylim(0, 1.)
ax.legend(loc='center right', bbox_to_anchor=(1.5, 0.5))
r = plt.setp(ax.xaxis.get_majorticklabels(), rotation=55, ha='right')<jupyter_output><empty_output><jupyter_text>Accuracy vs peak amplitude SNR<jupyter_code># inject in snr in perf tables
perfs_snr = perfs.copy()
perfs_snr['snr'] = None
for rec_name in rec_names:
    peak_snr = gt_info.loc[rec_name, 'peak_snr']
    perfs_snr.loc[perfs_snr['rec_name']==rec_name, 'snr'] = peak_snr

grid = sns.FacetGrid(perfs_snr, hue='sorter_name', col='sorter_name', col_wrap=3)
grid.map(plt.plot, 'snr', 'accuracy', marker="o", ls='None')<jupyter_output><empty_output><jupyter_text>Count how many well detected unitsIn dataframes['count_units'] we can check how many units in ground truth, in the tested sorter.And also *"how many weel detected"* which means settings some threholds and criteria.By default well detected units is for **accuracy>0.95** which is very strict.This can be explore with other thresholds.<jupyter_code># dataframes = aggregate_performances_table(study_folder, accuracy=0.9)
dataframes['count_units'][['num_well_detected']].reset_index().sort_values(by=['sorter_name', 'rec_name'])
df = dataframes['count_units'].group_by('sorter_name')['num_well_detected'].sum()
sns.barplot(df)<jupyter_output><empty_output><jupyter_text>Count how many redundant unitsSome sorter have the following pathology :  detecting a units twice or more!!!This is particulary true for template matching sorters because they have internally a redundant dictionary.We can measure this by counting how many units match the ground truth axcept the best one.<jupyter_code>dataframes['count_units'][['num_redundant']].reset_index().sort_values(by=['sorter_name', 'rec_name'])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11787072243346007, 'fraction_numerical': 0.023098859315589353, 'mean_word_length': 6.1039837947332884, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 50, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '611248', 'license': 'MIT', 'n_tokens_mistral': 3774, 'n_tokens_neox': 3540, 'n_words': 930, 'path': 'Spampinato_mice_retina_mea252ch/sorter_comparison.ipynb', 'repo_name': 'SpikeInterface/spikeresults', 'size': '10533'}