oups

Author: samuelgarcia

src/spikeinterface/benchmark/__init__.py

@@ -1,7 +1,9 @@
 """
 Module to benchmark:
-  * sorters
+  * sorters with or without ground-truth
   * some sorting components (clustering, motion, template matching)
 """
 
+from .residual_analysis import analyse_residual, make_residual_recording
 from .benchmark_sorter import SorterStudy
+from .benchmark_sorter_without_gt import SorterStudyWithoutGroundTruth

src/spikeinterface/benchmark/benchmark_base.py

@@ -9,7 +9,6 @@
 import time
 
 
-
 from spikeinterface.core import SortingAnalyzer, NumpySorting
 from spikeinterface.core.job_tools import fix_job_kwargs, split_job_kwargs
 from spikeinterface import load, create_sorting_analyzer, load_sorting_analyzer
@@ -130,7 +129,9 @@ def create(cls, study_folder, datasets={}, cases={}, levels=None):
                 else:
                     # some study/benchmark has no GT sorting
                     # in that case we still need an analyzer for internal API
-                    gt_sorting = NumpySorting.from_samples_and_labels([np.array([])], [np.array([])], rec.sampling_frequency, unit_ids=None)
+                    gt_sorting = NumpySorting.from_samples_and_labels(
+                        [np.array([])], [np.array([])], rec.sampling_frequency, unit_ids=None
+                    )
                     analyzer = create_sorting_analyzer(
                         gt_sorting, rec, sparse=False, format="binary_folder", folder=local_analyzer_folder
                     )

src/spikeinterface/benchmark/benchmark_sorter_without_gt.py

@@ -29,30 +29,24 @@ def run(self):
         self.result = {"sorting": sorting}
 
     def compute_result(self, residulal_peak_threshold=6, **job_kwargs):
-        
-        sorting = self.result['sorting']
-        analyzer = create_sorting_analyzer(
-            sorting, self.recording, sparse=True, format="memory", **job_kwargs
-        )
+
+        sorting = self.result["sorting"]
+        analyzer = create_sorting_analyzer(sorting, self.recording, sparse=True, format="memory", **job_kwargs)
         analyzer.compute("random_spikes")
         analyzer.compute("templates")
         analyzer.compute("noise_levels")
-        analyzer.compute(
-            {
-                "spike_amplitudes" : {},
-                "amplitude_scalings" : {"handle_collisions": False}
-            },
-            **job_kwargs)
+        analyzer.compute({"spike_amplitudes": {}, "amplitude_scalings": {"handle_collisions": False}}, **job_kwargs)
 
         analyzer.compute("quality_metrics", **job_kwargs)
 
         residual, peaks = analyse_residual(
-            analyzer, detect_peaks_kwargs=dict(
+            analyzer,
+            detect_peaks_kwargs=dict(
                 method="locally_exclusive",
                 peak_sign="neg",
                 detect_threshold=residulal_peak_threshold,
             ),
-            **job_kwargs
+            **job_kwargs,
         )
 
         self.result["sorter_analyzer"] = analyzer
@@ -66,11 +60,9 @@ def compute_result(self, residulal_peak_threshold=6, **job_kwargs):
         ("multi_comp", "pickle"),
         ("sorter_analyzer", "sorting_analyzer"),
         ("peaks_from_residual", "npy"),
-        
     ]
 
 
-
 class SorterStudyWithoutGroundTruth(BenchmarkStudy):
     """
     This class is an alternative to SorterStudy when the dataset do not have groundtruth.
@@ -87,19 +79,21 @@ def create_benchmark(self, key):
         sorter_folder = self.folder / "sorters" / self.key_to_str(key)
         benchmark = SorterBenchmarkWithoutGroundTruth(recording, gt_sorting, params, sorter_folder)
         return benchmark
-    
+
     def _get_comparison_groups(self):
         # multicomparison are done on all cases sharing the same dataset key.
         case_keys = list(self.cases.keys())
         groups = {}
         for case_key in case_keys:
-            data_key = self.cases[case_key]['dataset']
+            data_key = self.cases[case_key]["dataset"]
             if data_key not in groups:
                 groups[data_key] = []
             groups[data_key].append(case_key)
         return groups
 
-    def compute_results(self, case_keys=None, verbose=False, delta_time=0.4, match_score=0.5, chance_score=0.1, **result_params):
+    def compute_results(
+        self, case_keys=None, verbose=False, delta_time=0.4, match_score=0.5, chance_score=0.1, **result_params
+    ):
         # Here we need a hack because the results is not computed case by case but all at once
 
         assert case_keys is None, "SorterStudyWithoutGroundTruth do not permit compute_results for sub cases"
@@ -115,7 +109,7 @@ def compute_results(self, case_keys=None, verbose=False, delta_time=0.4, match_s
 
         for data_key, group in groups.items():
 
-            sorting_list = [self.get_result(key)['sorting'] for key in group]
+            sorting_list = [self.get_result(key)["sorting"] for key in group]
             name_list = [key for key in group]
             multi_comp = compare_multiple_sorters(
                 sorting_list,
@@ -132,7 +126,7 @@ def compute_results(self, case_keys=None, verbose=False, delta_time=0.4, match_s
             # and then the same multi comp is stored for each case_key
             for key in case_keys:
                 benchmark = self.benchmarks[key]
-                benchmark.result['multi_comp'] = multi_comp
+                benchmark.result["multi_comp"] = multi_comp
                 benchmark.save_result(self.folder / "results" / self.key_to_str(key))
 
     def plot_residual_peak_amplitudes(self, figsize=None):
@@ -143,7 +137,7 @@ def plot_residual_peak_amplitudes(self, figsize=None):
 
         for data_key, group in groups.items():
             fig, ax = plt.subplots(figsize=figsize)
-            
+
             lim0, lim1 = np.inf, -np.inf
 
             for key in group:
@@ -158,7 +152,6 @@ def plot_residual_peak_amplitudes(self, figsize=None):
             if lim0 > 0:
                 lim0 = 0
 
-
             for key in group:
                 peaks = self.get_result(key)["peaks_from_residual"]
                 print(peaks.size)
@@ -167,7 +160,7 @@ def plot_residual_peak_amplitudes(self, figsize=None):
                 ax.plot(bins[:-1], count, color=colors[key], label=self.cases[key]["label"])
 
             ax.legend()
-    
+
     # def plot_quality_metrics_comparison_on_agreement(self, qm_name='rp_contamination', figsize=None):
     #     import matplotlib.pyplot as plt
 
@@ -189,13 +182,13 @@ def plot_residual_peak_amplitudes(self, figsize=None):
 
     #                     multi_comp = self.get_result(key1)['multi_comp']
     #                     comp = multi_comp.comparisons[key1, key2]
-                        
+
     #                     match_12 = comp.hungarian_match_12
     #                     if match_12.dtype.kind =='i':
     #                         mask = match_12.values != -1
     #                     if match_12.dtype.kind =='U':
     #                         mask = match_12.values != ''
-                        
+
     #                     common_unit1_ids = match_12[mask].index
     #                     common_unit2_ids = match_12[mask].values
     #                     metrics1 = self.get_result(key1)["sorter_analyzer"].get_extension("quality_metrics").get_data()
@@ -215,5 +208,3 @@ def plot_residual_peak_amplitudes(self, figsize=None):
     #                         ax.set_yticks([])
     #                         ax.set_xticklabels([])
     #                         ax.set_yticklabels([])
-
-

src/spikeinterface/benchmark/residual_analysis.py

@@ -1,17 +1,15 @@
 from spikeinterface.core.generate import InjectTemplatesRecording
 
 
-
 def analyse_residual(
-        analyzer,
-        detect_peaks_kwargs=dict(
-            method="locally_exclusive",
-            peak_sign="both",
-            detect_threshold=6.,
-
-        ),
-        **job_kwargs
-    ):
+    analyzer,
+    detect_peaks_kwargs=dict(
+        method="locally_exclusive",
+        peak_sign="both",
+        detect_threshold=6.0,
+    ),
+    **job_kwargs,
+):
     """
     This create the residual by removing each spike from the recording.
     This take in account the spike amplitude scaling, analyzer need "amplitude_scalings" extensions.
@@ -41,8 +39,6 @@ def analyse_residual(
     return residual, peaks
 
 
-
-
 def make_residual_recording(analyzer):
     """
     This make a lazy recording residual from an anlyzer.
@@ -56,7 +52,7 @@ def make_residual_recording(analyzer):
     residual : Recording
         The resdiual
     """
-    
+
     templates = analyzer.get_extension("templates").get_templates(outputs="Templates")
     neg_templates_array = templates.templates_array.copy()
     neg_templates_array *= -1

src/spikeinterface/benchmark/tests/test_benchmark_sorter_without_gt.py

@@ -40,14 +40,12 @@ def _create_simple_study_no_gt(study_folder):
                 "sorter_name": "spykingcircus2",
             },
         },
-
     }
 
-    study = SorterStudyWithoutGroundTruth.create(
-        study_folder, datasets=datasets, cases=cases
-    )
+    study = SorterStudyWithoutGroundTruth.create(study_folder, datasets=datasets, cases=cases)
     # print(study)
 
+
 @pytest.mark.skip()
 def test_SorterStudyWithoutGroundTruth(create_simple_study):
     # job_kwargs = dict(n_jobs=2, chunk_duration="1s")
@@ -64,9 +62,10 @@ def test_SorterStudyWithoutGroundTruth(create_simple_study):
     print(study)
 
 
-
 if __name__ == "__main__":
-    study_folder_simple = Path(__file__).resolve().parents[4] / "cache_folder" / "benchmarks" / "test_SorterStudyWithoutGroundTruth"
+    study_folder_simple = (
+        Path(__file__).resolve().parents[4] / "cache_folder" / "benchmarks" / "test_SorterStudyWithoutGroundTruth"
+    )
     print(study_folder_simple)
     if study_folder_simple.exists():
         shutil.rmtree(study_folder_simple)

src/spikeinterface/benchmark/tests/test_residual_analysis.py

@@ -10,9 +10,10 @@
 
 job_kwargs = dict(n_jobs=-1, progress_bar=True)
 
+
 @pytest.mark.skip()
 def test_analyse_residual():
-    _, _ , analyzer = make_dataset()
+    _, _, analyzer = make_dataset()
     if not analyzer.has_extension("amplitude_scalings"):
         analyzer.compute("amplitude_scalings", **job_kwargs)
     print(analyzer)

src/spikeinterface/comparison/comparisontools.py

@@ -461,19 +461,19 @@ def make_possible_match(agreement_scores, min_score):
 def _empty_match_series(unit1_ids, unit2_ids):
     # construct empty series of match with the correct dtype for best match and hugarian match
     import pandas as pd
+
     match_12 = pd.Series(data=np.zeros(unit1_ids.size, dtype=unit2_ids.dtype), index=unit1_ids)
-    if unit2_ids.dtype.kind == 'i':
+    if unit2_ids.dtype.kind == "i":
         match_12[:] = -1
-    elif unit2_ids.dtype.kind == 'U':
-        match_12[:] = ''
-    elif unit2_ids.dtype.kind == 'O':
-        match_12[:] = ''
+    elif unit2_ids.dtype.kind == "U":
+        match_12[:] = ""
+    elif unit2_ids.dtype.kind == "O":
+        match_12[:] = ""
     else:
         raise ValueError("make_best_match or make_hungarian_match has unit_ids dtype wich are not  'i' or 'U'")
     return match_12
 
 
-
 def make_best_match(agreement_scores, min_score) -> "tuple[pd.Series, pd.Series]":
     """
     Given an agreement matrix and a min_score threshold.

src/spikeinterface/comparison/paircomparisons.py

@@ -203,7 +203,7 @@ def get_matching_unit_list2(self, unit2):
         return self.possible_match_21[unit2]
 
     def get_agreement_fraction(self, unit1=None, unit2=None):
-        if unit1 is None or unit1 == -1 or unit1 == '' or unit2 is None or unit2 == -1 or unit2 == '':
+        if unit1 is None or unit1 == -1 or unit1 == "" or unit2 is None or unit2 == -1 or unit2 == "":
             return 0
         else:
             return self.agreement_scores.at[unit1, unit2]
@@ -571,7 +571,7 @@ def get_false_positive_units(self, redundant_score=None):
         false_positive_ids = []
         for u2 in self.unit2_ids:
             if u2 not in matched_units2:
-                if self.best_match_21[u2] == -1 or self.best_match_21[u2] == '':
+                if self.best_match_21[u2] == -1 or self.best_match_21[u2] == "":
                     false_positive_ids.append(u2)
                 else:
                     u1 = self.best_match_21[u2]
@@ -615,7 +615,7 @@ def get_redundant_units(self, redundant_score=None):
         matched_units2 = list(self.hungarian_match_12.values)
         redundant_ids = []
         for u2 in self.unit2_ids:
-            if u2 not in matched_units2 and self.best_match_21[u2] != -1 and self.best_match_21[u2] != '':
+            if u2 not in matched_units2 and self.best_match_21[u2] != -1 and self.best_match_21[u2] != "":
                 u1 = self.best_match_21[u2]
                 if u2 != self.best_match_12[u1]:
                     score = self.agreement_scores.at[u1, u2]