Merge branch 'main' of github.com:SpikeInterface/spikeinterface into rtd-icons

Author: alejoe91

doc/api.rst

@@ -277,8 +277,6 @@ spikeinterface.comparison
 
     .. autoclass:: CollisionGTComparison
     .. autoclass:: CorrelogramGTComparison
-    .. autoclass:: CollisionGTStudy
-    .. autoclass:: CorrelogramGTStudy
 
 
 

src/spikeinterface/comparison/basecomparison.py

@@ -283,9 +283,10 @@ class MixinSpikeTrainComparison:
        * n_jobs
     """
 
-    def __init__(self, delta_time=0.4, n_jobs=-1):
+    def __init__(self, delta_time=0.4, agreement_method="count", n_jobs=-1):
         self.delta_time = delta_time
         self.n_jobs = n_jobs
+        self.agreement_method = agreement_method
         self.sampling_frequency = None
         self.delta_frames = None
 

src/spikeinterface/comparison/multicomparisons.py

@@ -35,6 +35,9 @@ class MultiSortingComparison(BaseMultiComparison, MixinSpikeTrainComparison):
         Minimum agreement score to match units
     chance_score : float, default: 0.1
         Minimum agreement score to for a possible match
+    agreement_method : "count" | "distance", default: "count"
+        The method to compute agreement scores. The "count" method computes agreement scores from spike counts.
+        The "distance" method computes agreement scores from spike time distance functions.
     n_jobs : int, default: -1
        Number of cores to use in parallel. Uses all available if -1
     spiketrain_mode : "union" | "intersection", default: "union"
@@ -60,6 +63,7 @@ def __init__(
         delta_time=0.4,  # sampling_frequency=None,
         match_score=0.5,
         chance_score=0.1,
+        agreement_method="count",
         n_jobs=-1,
         spiketrain_mode="union",
         verbose=False,
@@ -75,7 +79,9 @@ def __init__(
             chance_score=chance_score,
             verbose=verbose,
         )
-        MixinSpikeTrainComparison.__init__(self, delta_time=delta_time, n_jobs=n_jobs)
+        MixinSpikeTrainComparison.__init__(
+            self, delta_time=delta_time, agreement_method=agreement_method, n_jobs=n_jobs
+        )
         self.set_frames_and_frequency(self.object_list)
         self._spiketrain_mode = spiketrain_mode
         self._spiketrains = None
@@ -93,6 +99,8 @@ def _compare_ij(self, i, j):
             sorting2_name=self.name_list[j],
             delta_time=self.delta_time,
             match_score=self.match_score,
+            chance_score=self.chance_score,
+            agreement_method=self.agreement_method,
             n_jobs=self.n_jobs,
             verbose=False,
         )

src/spikeinterface/comparison/tests/test_multisortingcomparison.py

@@ -41,12 +41,15 @@ def test_compare_multiple_sorters(setup_module):
     )
     msc = compare_multiple_sorters([sorting1, sorting2, sorting3], verbose=True)
     msc_shuffle = compare_multiple_sorters([sorting3, sorting1, sorting2])
+    msc_dist = compare_multiple_sorters([sorting3, sorting1, sorting2], agreement_method="distance")
 
     agr = msc._do_agreement_matrix()
     agr_shuffle = msc_shuffle._do_agreement_matrix()
+    agr_dist = msc_dist._do_agreement_matrix()
 
     print(agr)
     print(agr_shuffle)
+    print(agr_dist)
 
     assert len(msc.get_agreement_sorting(minimum_agreement_count=3).get_unit_ids()) == 3
     assert len(msc.get_agreement_sorting(minimum_agreement_count=2).get_unit_ids()) == 5
@@ -57,7 +60,14 @@ def test_compare_multiple_sorters(setup_module):
     assert len(msc.get_agreement_sorting(minimum_agreement_count=2).get_unit_ids()) == len(
         msc_shuffle.get_agreement_sorting(minimum_agreement_count=2).get_unit_ids()
     )
+    assert len(msc.get_agreement_sorting(minimum_agreement_count=3).get_unit_ids()) == len(
+        msc_dist.get_agreement_sorting(minimum_agreement_count=3).get_unit_ids()
+    )
+    assert len(msc.get_agreement_sorting(minimum_agreement_count=2).get_unit_ids()) == len(
+        msc_dist.get_agreement_sorting(minimum_agreement_count=2).get_unit_ids()
+    )
     assert len(msc.get_agreement_sorting().get_unit_ids()) == len(msc_shuffle.get_agreement_sorting().get_unit_ids())
+    assert len(msc.get_agreement_sorting().get_unit_ids()) == len(msc_dist.get_agreement_sorting().get_unit_ids())
     agreement_2 = msc.get_agreement_sorting(minimum_agreement_count=2, minimum_agreement_count_only=True)
     assert np.all([agreement_2.get_unit_property(u, "agreement_number")] == 2 for u in agreement_2.get_unit_ids())
 

src/spikeinterface/core/baserecordingsnippets.py

@@ -560,7 +560,7 @@ def split_by(self, property="group", outputs="dict"):
             recordings = []
         elif outputs == "dict":
             recordings = {}
-        for value in np.unique(values):
+        for value in np.unique(values).tolist():
             (inds,) = np.nonzero(values == value)
             new_channel_ids = self.get_channel_ids()[inds]
             subrec = self.select_channels(new_channel_ids)

src/spikeinterface/core/sparsity.py

@@ -15,6 +15,8 @@
     method : str
         * "best_channels" : N best channels with the largest amplitude. Use the "num_channels" argument to specify the
                            number of channels.
+        * "closest_channels" : N closest channels according to the distance. Use the "num_channels" argument to specify the
+                           number of channels.
         * "radius" : radius around the best channel. Use the "radius_um" argument to specify the radius in um.
         * "snr" : threshold based on template signal-to-noise ratio. Use the "threshold" argument
                  to specify the SNR threshold (in units of noise levels) and the "amplitude_mode" argument
@@ -320,6 +322,39 @@ def from_best_channels(
             mask[unit_ind, chan_inds] = True
         return cls(mask, templates_or_sorting_analyzer.unit_ids, templates_or_sorting_analyzer.channel_ids)
 
+    ## Some convinient function to compute sparsity from several strategy
+    @classmethod
+    def from_closest_channels(cls, templates_or_sorting_analyzer, num_channels):
+        """
+        Construct sparsity from N closest channels
+        Use the "num_channels" argument to specify the number of channels.
+
+        Parameters
+        ----------
+        templates_or_sorting_analyzer : Templates | SortingAnalyzer
+            A Templates or a SortingAnalyzer object.
+        num_channels : int
+            Number of channels for "best_channels" method.
+
+        Returns
+        -------
+        sparsity : ChannelSparsity
+            The estimated sparsity
+        """
+        from .template_tools import get_template_amplitudes
+
+        mask = np.zeros(
+            (templates_or_sorting_analyzer.unit_ids.size, templates_or_sorting_analyzer.channel_ids.size), dtype="bool"
+        )
+        channel_locations = templates_or_sorting_analyzer.get_channel_locations()
+        distances = np.linalg.norm(channel_locations[:, np.newaxis] - channel_locations[np.newaxis, :], axis=2)
+
+        for unit_ind, unit_id in enumerate(templates_or_sorting_analyzer.unit_ids):
+            chan_inds = np.argsort(distances[unit_ind])
+            chan_inds = chan_inds[:num_channels]
+            mask[unit_ind, chan_inds] = True
+        return cls(mask, templates_or_sorting_analyzer.unit_ids, templates_or_sorting_analyzer.channel_ids)
+
     @classmethod
     def from_radius(cls, templates_or_sorting_analyzer, radius_um, peak_sign="neg"):
         """
@@ -600,7 +635,9 @@ def create_dense(cls, sorting_analyzer):
 def compute_sparsity(
     templates_or_sorting_analyzer: "Templates | SortingAnalyzer",
     noise_levels: np.ndarray | None = None,
-    method: "radius" | "best_channels" | "snr" | "amplitude" | "energy" | "by_property" | "ptp" = "radius",
+    method: (
+        "radius" | "best_channels" | "closest_channels" | "snr" | "amplitude" | "energy" | "by_property" | "ptp"
+    ) = "radius",
     peak_sign: "neg" | "pos" | "both" = "neg",
     num_channels: int | None = 5,
     radius_um: float | None = 100.0,
@@ -635,7 +672,7 @@ def compute_sparsity(
         # to keep backward compatibility
         templates_or_sorting_analyzer = templates_or_sorting_analyzer.sorting_analyzer
 
-    if method in ("best_channels", "radius", "snr", "amplitude", "ptp"):
+    if method in ("best_channels", "closest_channels", "radius", "snr", "amplitude", "ptp"):
         assert isinstance(
             templates_or_sorting_analyzer, (Templates, SortingAnalyzer)
         ), f"compute_sparsity(method='{method}') need Templates or SortingAnalyzer"
@@ -647,6 +684,9 @@ def compute_sparsity(
     if method == "best_channels":
         assert num_channels is not None, "For the 'best_channels' method, 'num_channels' needs to be given"
         sparsity = ChannelSparsity.from_best_channels(templates_or_sorting_analyzer, num_channels, peak_sign=peak_sign)
+    elif method == "closest_channels":
+        assert num_channels is not None, "For the 'closest_channels' method, 'num_channels' needs to be given"
+        sparsity = ChannelSparsity.from_closest_channels(templates_or_sorting_analyzer, num_channels)
     elif method == "radius":
         assert radius_um is not None, "For the 'radius' method, 'radius_um' needs to be given"
         sparsity = ChannelSparsity.from_radius(templates_or_sorting_analyzer, radius_um, peak_sign=peak_sign)
@@ -698,7 +738,7 @@ def estimate_sparsity(
     num_spikes_for_sparsity: int = 100,
     ms_before: float = 1.0,
     ms_after: float = 2.5,
-    method: "radius" | "best_channels" | "amplitude" | "snr" | "by_property" | "ptp" = "radius",
+    method: "radius" | "best_channels" | "closest_channels" | "amplitude" | "snr" | "by_property" | "ptp" = "radius",
     peak_sign: "neg" | "pos" | "both" = "neg",
     radius_um: float = 100.0,
     num_channels: int = 5,
@@ -747,7 +787,7 @@ def estimate_sparsity(
     # Can't be done at module because this is a cyclic import, too bad
     from .template import Templates
 
-    assert method in ("radius", "best_channels", "snr", "amplitude", "by_property", "ptp"), (
+    assert method in ("radius", "best_channels", "closest_channels", "snr", "amplitude", "by_property", "ptp"), (
         f"method={method} is not available for `estimate_sparsity()`. "
         "Available methods are 'radius', 'best_channels', 'snr', 'amplitude', 'by_property', 'ptp' (deprecated)"
     )
@@ -802,6 +842,9 @@ def estimate_sparsity(
             sparsity = ChannelSparsity.from_best_channels(
                 templates, num_channels, peak_sign=peak_sign, amplitude_mode=amplitude_mode
             )
+        elif method == "closest_channels":
+            assert num_channels is not None, "For the 'closest_channels' method, 'num_channels' needs to be given"
+            sparsity = ChannelSparsity.from_closest_channels(templates, num_channels)
         elif method == "radius":
             assert radius_um is not None, "For the 'radius' method, 'radius_um' needs to be given"
             sparsity = ChannelSparsity.from_radius(templates, radius_um, peak_sign=peak_sign)

src/spikeinterface/core/tests/test_sparsity.py

@@ -195,6 +195,21 @@ def test_estimate_sparsity():
     )
     assert np.array_equal(np.sum(sparsity.mask, axis=1), np.ones(num_units) * 3)
 
+    # closest_channels : the mask should exactly 3 channels per units
+    sparsity = estimate_sparsity(
+        sorting,
+        recording,
+        num_spikes_for_sparsity=50,
+        ms_before=1.0,
+        ms_after=2.0,
+        method="closest_channels",
+        num_channels=3,
+        chunk_duration="1s",
+        progress_bar=True,
+        n_jobs=1,
+    )
+    assert np.array_equal(np.sum(sparsity.mask, axis=1), np.ones(num_units) * 3)
+
     # by_property
     sparsity = estimate_sparsity(
         sorting,
@@ -287,6 +302,7 @@ def test_compute_sparsity():
     # using object SortingAnalyzer
     sparsity = compute_sparsity(sorting_analyzer, method="best_channels", num_channels=2, peak_sign="neg")
     sparsity = compute_sparsity(sorting_analyzer, method="radius", radius_um=50.0, peak_sign="neg")
+    sparsity = compute_sparsity(sorting_analyzer, method="closest_channels", num_channels=2)
     sparsity = compute_sparsity(sorting_analyzer, method="snr", threshold=5, peak_sign="neg")
     sparsity = compute_sparsity(
         sorting_analyzer, method="snr", threshold=5, peak_sign="neg", amplitude_mode="peak_to_peak"
@@ -304,6 +320,7 @@ def test_compute_sparsity():
     sparsity = compute_sparsity(templates, method="radius", radius_um=50.0, peak_sign="neg")
     sparsity = compute_sparsity(templates, method="snr", noise_levels=noise_levels, threshold=5, peak_sign="neg")
     sparsity = compute_sparsity(templates, method="amplitude", threshold=5, amplitude_mode="peak_to_peak")
+    sparsity = compute_sparsity(templates, method="closest_channels", num_channels=2)
 
     with pytest.warns(DeprecationWarning):
         sparsity = compute_sparsity(templates, method="ptp", noise_levels=noise_levels, threshold=5)

src/spikeinterface/extractors/neoextractors/openephys.py

@@ -164,7 +164,7 @@ def __init__(
             **neo_kwargs,
         )
         # get streams to find correct probe
-        stream_names, stream_ids = self.get_streams(folder_path, experiment_names)
+        stream_names, stream_ids = self.get_streams(folder_path, load_sync_channel, experiment_names)
         if stream_name is None and stream_id is None:
             stream_name = stream_names[0]
         elif stream_name is None:

src/spikeinterface/postprocessing/localization_tools.py

@@ -363,7 +363,7 @@ def solve_monopolar_triangulation(wf_data, local_contact_locations, max_distance
             output = scipy.optimize.least_squares(estimate_distance_error, x0=x0, bounds=bounds, args=args)
             return tuple(output["x"])
         except Exception as e:
-            print(f"scipy.optimize.least_squares error: {e}")
+            warnings.warn(f"scipy.optimize.least_squares error: {e}")
             return (np.nan, np.nan, np.nan, np.nan)
 
     if optimizer == "minimize_with_log_penality":
@@ -378,7 +378,7 @@ def solve_monopolar_triangulation(wf_data, local_contact_locations, max_distance
             alpha = (wf_data * q).sum() / np.square(q).sum()
             return (*output["x"], alpha)
         except Exception as e:
-            print(f"scipy.optimize.minimize error: {e}")
+            warnings.warn(f"scipy.optimize.minimize error: {e}")
             return (np.nan, np.nan, np.nan, np.nan)