Repository: SpikeInterface/spikeinterface
Pull request #3509 (2024-10-28 14:38:02)

Add `shift start time` function.
This PR adds a 'shift start time' function that a user can use to set a new start time for their recording.

This follows from discussions on the time API #3117 and #3157. Briefly, there are two ways to handle time in SI, a `t_start` recording attribute that is used to set the start time. Under this regime, the (regularly) spaced timeseries is regenerated on the fly from `t_start`, num samples and sampling frequency. Alternatively, a (more memory heavy) `time_vector` can be set to explicitly set the time array for a recording. As I understand it this is useful in the case there are gaps in the recording or there is some clock drift. From the user perspective, they always get a time vector back from `get_times()`.

This PR implements a `shift_start_time` as suggested by @h-mayorquin as a convenient way for users to shift the start time of their recording, abstracting away the technical details of the implementation under the hood (`t_start` vs `time_vector`).

I have updated the PR #3072  containing documentation on times, the relevant proposed section can be found [here](https://spikeinterface--3072.org.readthedocs.build/en/3072/tutorials/core/plot_6_handle_times_ver2.html#shifting-the-start-time).


Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 4c36d6fb8..1e133694b 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -6,7 +6,7 @@ repos:
     -   id: end-of-file-fixer
     -   id: trailing-whitespace
 -   repo: https://github.com/psf/black
-    rev: 24.10.0
+    rev: 24.8.0
     hooks:
     -   id: black
         files: ^src/
diff --git a/doc/modules/core.rst b/doc/modules/core.rst
index 5df9a7e6b..8aa1815a5 100644
--- a/doc/modules/core.rst
+++ b/doc/modules/core.rst
@@ -385,7 +385,7 @@ and merging unit groups.
 
     sorting_analyzer_select = sorting_analyzer.select_units(unit_ids=[0, 1, 2, 3])
     sorting_analyzer_remove = sorting_analyzer.remove_units(remove_unit_ids=[0])
-    sorting_analyzer_merge = sorting_analyzer.merge_units([[0, 1], [2, 3]])
+    sorting_analyzer_merge = sorting_analyzer.merge_units([0, 1], [2, 3])
 
 All computed extensions will be automatically propagated or merged when curating. Please refer to the
 :ref:`modules/curation:Curation module` documentation for more information.
diff --git a/doc/modules/curation.rst b/doc/modules/curation.rst
index 37de99280..d115b33e4 100644
--- a/doc/modules/curation.rst
+++ b/doc/modules/curation.rst
@@ -88,7 +88,7 @@ The ``censored_period_ms`` parameter is the time window in milliseconds to consi
 The :py:func:`~spikeinterface.curation.remove_redundand_units` function removes
 redundant units from the sorting output. Redundant units are units that share over
 a certain percentage of spikes, by default 80%.
-The function can act both on a ``BaseSorting`` or a ``SortingAnalyzer`` object.
+The function can acto both on a ``BaseSorting`` or a ``SortingAnalyzer`` object.
 
 .. code-block:: python
 
@@ -102,18 +102,13 @@ The function can act both on a ``BaseSorting`` or a ``SortingAnalyzer`` object.
     )
 
     # remove redundant units from SortingAnalyzer object
-    # note this returns a cleaned sorting
-    clean_sorting = remove_redundant_units(
+    clean_sorting_analyzer = remove_redundant_units(
         sorting_analyzer,
         duplicate_threshold=0.9,
         remove_strategy="min_shift"
     )
-    # in order to have a SortingAnalyer with only the non-redundant units one must
-    # select the designed units remembering to give format and folder if one wants
-    # a persistent SortingAnalyzer.
-    clean_sorting_analyzer = sorting_analyzer.select_units(clean_sorting.unit_ids)
 
-We recommend using the ``SortingAnalyzer`` approach, since the ``min_shift`` strategy keeps
+We recommend usinf the ``SortingAnalyzer`` approach, since the ``min_shift`` strategy keeps
 the unit (among the redundant ones), with a better template alignment.
 
 
diff --git a/pyproject.toml b/pyproject.toml
index fc09ad919..a43ab63c8 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -7,7 +7,7 @@ authors = [
 ]
 description = "Python toolkit for analysis, visualization, and comparison of spike sorting output"
 readme = "README.md"
-requires-python = ">=3.9,<3.13"  # Only numpy 2.1 supported on python 3.13 for windows. We need to wait for fix on neo
+requires-python = ">=3.9,<4.0"
 classifiers = [
     "Programming Language :: Python :: 3 :: Only",
     "License :: OSI Approved :: MIT License",
diff --git a/src/spikeinterface/benchmark/benchmark_matching.py b/src/spikeinterface/benchmark/benchmark_matching.py
index 1934b65ef..c53567f46 100644
--- a/src/spikeinterface/benchmark/benchmark_matching.py
+++ b/src/spikeinterface/benchmark/benchmark_matching.py
@@ -33,7 +33,7 @@ class MatchingBenchmark(Benchmark):
         sorting["unit_index"] = spikes["cluster_index"]
         sorting["segment_index"] = spikes["segment_index"]
         sorting = NumpySorting(sorting, self.recording.sampling_frequency, unit_ids)
-        self.result = {"sorting": sorting, "spikes": spikes}
+        self.result = {"sorting": sorting}
         self.result["templates"] = self.templates
 
     def compute_result(self, with_collision=False, **result_params):
@@ -45,7 +45,6 @@ class MatchingBenchmark(Benchmark):
 
     _run_key_saved = [
         ("sorting", "sorting"),
-        ("spikes", "npy"),
         ("templates", "zarr_templates"),
     ]
     _result_key_saved = [("gt_collision", "pickle"), ("gt_comparison", "pickle")]
@@ -72,15 +71,9 @@ class MatchingStudy(BenchmarkStudy):
 
         return plot_performances_vs_snr(self, **kwargs)
 
-    def plot_performances_comparison(self, **kwargs):
-        from .benchmark_plot_tools import plot_performances_comparison
-
-        return plot_performances_comparison(self, **kwargs)
-
     def plot_collisions(self, case_keys=None, figsize=None):
         if case_keys is None:
             case_keys = list(self.cases.keys())
-        import matplotlib.pyplot as plt
 
         fig, axs = plt.subplots(ncols=len(case_keys), nrows=1, figsize=figsize, squeeze=False)
 
@@ -97,6 +90,70 @@ class MatchingStudy(BenchmarkStudy):
 
         return fig
 
+    def plot_comparison_matching(
+        self,
+        case_keys=None,
+        performance_names=["accuracy", "recall", "precision"],
+        colors=["g", "b", "r"],
+        ylim=(-0.1, 1.1),
+        figsize=None,
+    ):
+
+        if case_keys is None:
+            case_keys = list(self.cases.keys())
+
+        num_methods = len(case_keys)
+        import pylab as plt
+
+        fig, axs = plt.subplots(ncols=num_methods, nrows=num_methods, figsize=(10, 10))
+        for i, key1 in enumerate(case_keys):
+            for j, key2 in enumerate(case_keys):
+                if len(axs.shape) > 1:
+                    ax = axs[i, j]
+                else:
+                    ax = axs[j]
+                comp1 = self.get_result(key1)["gt_comparison"]
+                comp2 = self.get_result(key2)["gt_comparison"]
+                if i <= j:
+                    for performance, color in zip(performance_names, colors):
+                        perf1 = comp1.get_performance()[performance]
+                        perf2 = comp2.get_performance()[performance]
+                        ax.plot(perf2, perf1, ".", label=performance, color=color)
+
+                    ax.plot([0, 1], [0, 1], "k--", alpha=0.5)
+                    ax.set_ylim(ylim)
+                    ax.set_xlim(ylim)
+                    ax.spines[["right", "top"]].set_visible(False)
+                    ax.set_aspect("equal")
+
+                    label1 = self.cases[key1]["label"]
+                    label2 = self.cases[key2]["label"]
+                    if j == i:
+                        ax.set_ylabel(f"{label1}")
+                    else:
+                        ax.set_yticks([])
+                    if i == j:
+                        ax.set_xlabel(f"{label2}")
+                    else:
+                        ax.set_xticks([])
+                    if i == num_methods - 1 and j == num_methods - 1:
+                        patches = []
+                        import matplotlib.patches as mpatches
+
+                        for color, name in zip(colors, performance_names):
+                            patches.append(mpatches.Patch(color=color, label=name))
+                        ax.legend(handles=patches, bbox_to_anchor=(1.05, 1), loc="upper left", borderaxespad=0.0)
+                else:
+                    ax.spines["bottom"].set_visible(False)
+                    ax.spines["left"].set_visible(False)
+                    ax.spines["top"].set_visible(False)
+                    ax.spines["right"].set_visible(False)
+                    ax.set_xticks([])
+                    ax.set_yticks([])
+        plt.tight_layout(h_pad=0, w_pad=0)
+
+        return fig
+
     def get_count_units(self, case_keys=None, well_detected_score=None, redundant_score=None, overmerged_score=None):
         import pandas as pd
 
@@ -139,7 +196,6 @@ class MatchingStudy(BenchmarkStudy):
         plot_study_unit_counts(self, case_keys, figsize=figsize)
 
     def plot_unit_losses(self, before, after, metric=["precision"], figsize=None):
-        import matplotlib.pyplot as plt
 
         fig, axs = plt.subplots(ncols=1, nrows=len(metric), figsize=figsize, squeeze=False)
 
diff --git a/src/spikeinterface/benchmark/benchmark_plot_tools.py b/src/spikeinterface/benchmark/benchmark_plot_tools.py
index e15636eba..a6e9b6dac 100644
--- a/src/spikeinterface/benchmark/benchmark_plot_tools.py
+++ b/src/spikeinterface/benchmark/benchmark_plot_tools.py
@@ -235,71 +235,9 @@ def plot_performances_vs_snr(study, case_keys=None, figsize=None, metrics=["accu
             ax.scatter(x, y, marker=".", label=label)
             ax.set_title(k)
 
-        ax.set_ylim(-0.05, 1.05)
+        ax.set_ylim(0, 1.05)
 
         if count == 2:
             ax.legend()
 
     return fig
-
-
-def plot_performances_comparison(
-    study,
-    case_keys=None,
-    figsize=None,
-    metrics=["accuracy", "recall", "precision"],
-    colors=["g", "b", "r"],
-    ylim=(-0.1, 1.1),
-):
-    import matplotlib.pyplot as plt
-
-    if case_keys is None:
-        case_keys = list(study.cases.keys())
-
-    num_methods = len(case_keys)
-    assert num_methods >= 2, "plot_performances_comparison need at least 2 cases!"
-
-    fig, axs = plt.subplots(ncols=num_methods - 1, nrows=num_methods - 1, figsize=(10, 10), squeeze=False)
-    for i, key1 in enumerate(case_keys):
-        for j, key2 in enumerate(case_keys):
-
-            if i < j:
-                ax = axs[i, j - 1]
-
-                comp1 = study.get_result(key1)["gt_comparison"]
-                comp2 = study.get_result(key2)["gt_comparison"]
-
-                for performance, color in zip(metrics, colors):
-                    perf1 = comp1.get_performance()[performance]
-                    perf2 = comp2.get_performance()[performance]
-                    ax.scatter(perf2, perf1, marker=".", label=performance, color=color)
-
-                ax.plot([0, 1], [0, 1], "k--", alpha=0.5)
-                ax.set_ylim(ylim)
-                ax.set_xlim(ylim)
-                ax.spines[["right", "top"]].set_visible(False)
-                ax.set_aspect("equal")
-
-                label1 = study.cases[key1]["label"]
-                label2 = study.cases[key2]["label"]
-
-                if i == j - 1:
-                    ax.set_xlabel(label2)
-                    ax.set_ylabel(label1)
-
-            else:
-                if j >= 1 and i < num_methods - 1:
-                    ax = axs[i, j - 1]
-                    ax.spines[["right", "top", "left", "bottom"]].set_visible(False)
-                    ax.set_xticks([])
-                    ax.set_yticks([])
-
-    ax = axs[num_methods - 2, 0]
-    patches = []
-    from matplotlib.patches import Patch
-
-    for color, name in zip(colors, metrics):
-        patches.append(Patch(color=color, label=name))
-    ax.legend(handles=patches)
-    fig.tight_layout()
-    return fig
diff --git a/src/spikeinterface/core/analyzer_extension_core.py b/src/spikeinterface/core/analyzer_extension_core.py
index 447bbe562..bc5de63d0 100644
--- a/src/spikeinterface/core/analyzer_extension_core.py
+++ b/src/spikeinterface/core/analyzer_extension_core.py
@@ -22,23 +22,21 @@ from .sorting_tools import random_spikes_selection
 
 class ComputeRandomSpikes(AnalyzerExtension):
     """
-    AnalyzerExtension that select somes random spikes.
-    This allows for a subsampling of spikes for further calculations and is important
-    for managing that amount of memory and speed of computation in the analyzer.
+    AnalyzerExtension that select some random spikes.
 
     This will be used by the `waveforms`/`templates` extensions.
 
-    This internally uses `random_spikes_selection()` parameters.
+    This internally use `random_spikes_selection()` parameters are the same.
 
     Parameters
     ----------
-    method : "uniform" | "all", default: "uniform"
+    method: "uniform" | "all", default: "uniform"
         The method to select the spikes
-    max_spikes_per_unit : int, default: 500
+    max_spikes_per_unit: int, default: 500
         The maximum number of spikes per unit, ignored if method="all"
-    margin_size : int, default: None
+    margin_size: int, default: None
         A margin on each border of segments to avoid border spikes, ignored if method="all"
-    seed : int or None, default: None
+    seed: int or None, default: None
         A seed for the random generator, ignored if method="all"
 
     Returns
@@ -106,7 +104,7 @@ class ComputeRandomSpikes(AnalyzerExtension):
         return self._some_spikes
 
     def get_selected_indices_in_spike_train(self, unit_id, segment_index):
-        # useful for WaveformExtractor backwards compatibility
+        # usefull for Waveforms extractor backwars compatibility
         # In Waveforms extractor "selected_spikes" was a dict (key: unit_id) of list (segment_index) of indices of spikes in spiketrain
         sorting = self.sorting_analyzer.sorting
         random_spikes_indices = self.data["random_spikes_indices"]
@@ -135,16 +133,16 @@ class ComputeWaveforms(AnalyzerExtension):
 
     Parameters
     ----------
-    ms_before : float, default: 1.0
+    ms_before: float, default: 1.0
         The number of ms to extract before the spike events
-    ms_after : float, default: 2.0
+    ms_after: float, default: 2.0
         The number of ms to extract after the spike events
-    dtype : None | dtype, default: None
+    dtype: None | dtype, default: None
         The dtype of the waveforms. If None, the dtype of the recording is used.
 
     Returns
     -------
-    waveforms : np.ndarray
+    waveforms: np.ndarray
         Array with computed waveforms with shape (num_random_spikes, num_samples, num_channels)
     """
 
@@ -382,12 +380,7 @@ class ComputeTemplates(AnalyzerExtension):
         assert isinstance(operators, list)
         for operator in operators:
             if isinstance(operator, str):
-                if operator not in ("average", "std", "median", "mad"):
-                    error_msg = (
-                        f"You have entered an operator {operator} in your `operators` argument which is "
-                        f"not supported. Please use any of ['average', 'std', 'median', 'mad'] instead."
-                    )
-                    raise ValueError(error_msg)
+                assert operator in ("average", "std", "median", "mad")
             else:
                 assert isinstance(operator, (list, tuple))
                 assert len(operator) == 2
@@ -412,13 +405,9 @@ class ComputeTemplates(AnalyzerExtension):
             self._compute_and_append_from_waveforms(self.params["operators"])
 
         else:
-            bad_operator_list = [
-                operator for operator in self.params["operators"] if operator not in ("average", "std")
-            ]
-            if len(bad_operator_list) > 0:
-                raise ValueError(
-                    f"Computing templates with operators {bad_operator_list} requires the 'waveforms' extension"
-                )
+            for operator in self.params["operators"]:
+                if operator not in ("average", "std"):
+                    raise ValueError(f"Computing templates with operators {operator} needs the 'waveforms' extension")
 
             recording = self.sorting_analyzer.recording
             sorting = self.sorting_analyzer.sorting
@@ -452,7 +441,7 @@ class ComputeTemplates(AnalyzerExtension):
 
     def _compute_and_append_from_waveforms(self, operators):
         if not self.sorting_analyzer.has_extension("waveforms"):
-            raise ValueError(f"Computing templates with operators {operators} requires the 'waveforms' extension")
+            raise ValueError(f"Computing templates with operators {operators} needs the 'waveforms' extension")
 
         unit_ids = self.sorting_analyzer.unit_ids
         channel_ids = self.sorting_analyzer.channel_ids
@@ -477,7 +466,7 @@ class ComputeTemplates(AnalyzerExtension):
 
         assert self.sorting_analyzer.has_extension(
             "random_spikes"
-        ), "compute 'templates' requires the random_spikes extension. You can run sorting_analyzer.compute('random_spikes')"
+        ), "compute templates requires the random_spikes extension. You can run sorting_analyzer.get_random_spikes()"
         some_spikes = self.sorting_analyzer.get_extension("random_spikes").get_random_spikes()
         for unit_index, unit_id in enumerate(unit_ids):
             spike_mask = some_spikes["unit_index"] == unit_index
@@ -560,17 +549,9 @@ class ComputeTemplates(AnalyzerExtension):
         if operator != "percentile":
             key = operator
         else:
-            assert percentile is not None, "You must provide percentile=... if `operator=percentile`"
+            assert percentile is not None, "You must provide percentile=..."
             key = f"percentile_{percentile}"
 
-        if key not in self.data.keys():
-            error_msg = (
-                f"You have entered `operator={key}`, but the only operators calculated are "
-                f"{list(self.data.keys())}. Please use one of these as your `operator` in the "
-                f"`get_data` function."
-            )
-            raise ValueError(error_msg)
-
         templates_array = self.data[key]
 
         if outputs == "numpy":
@@ -585,7 +566,7 @@ class ComputeTemplates(AnalyzerExtension):
                 probe=self.sorting_analyzer.get_probe(),
             )
         else:
-            raise ValueError("outputs must be `numpy` or `Templates`")
+            raise ValueError("outputs must be numpy or Templates")
 
     def get_templates(self, unit_ids=None, operator="average", percentile=None, save=True, outputs="numpy"):
         """
@@ -595,26 +576,26 @@ class ComputeTemplates(AnalyzerExtension):
 
         Parameters
         ----------
-        unit_ids : list or None
+        unit_ids: list or None
             Unit ids to retrieve waveforms for
-        operator : "average" | "median" | "std" | "percentile", default: "average"
+        operator: "average" | "median" | "std" | "percentile", default: "average"
             The operator to compute the templates
-        percentile : float, default: None
+        percentile: float, default: None
             Percentile to use for operator="percentile"
-        save : bool, default: True
+        save: bool, default True
             In case, the operator is not computed yet it can be saved to folder or zarr
-        outputs : "numpy" | "Templates", default: "numpy"
+        outputs: "numpy" | "Templates"
             Whether to return a numpy array or a Templates object
 
         Returns
         -------
-        templates : np.array | Templates
+        templates: np.array
             The returned templates (num_units, num_samples, num_channels)
         """
         if operator != "percentile":
             key = operator
         else:
-            assert percentile is not None, "You must provide percentile=... if `operator='percentile'`"
+            assert percentile is not None, "You must provide percentile=..."
             key = f"pencentile_{percentile}"
 
         if key in self.data:
@@ -651,7 +632,7 @@ class ComputeTemplates(AnalyzerExtension):
                 is_scaled=self.sorting_analyzer.return_scaled,
             )
         else:
-            raise ValueError("`outputs` must be 'numpy' or 'Templates'")
+            raise ValueError("outputs must be numpy or Templates")
 
     def get_unit_template(self, unit_id, operator="average"):
         """
@@ -661,7 +642,7 @@ class ComputeTemplates(AnalyzerExtension):
         ----------
         unit_id: str | int
             Unit id to retrieve waveforms for
-        operator: str, default: "average"
+        operator: str
              The operator to compute the templates
 
         Returns
@@ -710,23 +691,22 @@ class ComputeNoiseLevels(AnalyzerExtension):
     need_recording = True
     use_nodepipeline = False
     need_job_kwargs = False
-    need_backward_compatibility_on_load = True
 
     def __init__(self, sorting_analyzer):
         AnalyzerExtension.__init__(self, sorting_analyzer)
 
-    def _set_params(self, **noise_level_params):
-        params = noise_level_params.copy()
+    def _set_params(self, num_chunks_per_segment=20, chunk_size=10000, seed=None):
+        params = dict(num_chunks_per_segment=num_chunks_per_segment, chunk_size=chunk_size, seed=seed)
         return params
 
     def _select_extension_data(self, unit_ids):
-        # this does not depend on units
+        # this do not depend on units
         return self.data
 
     def _merge_extension_data(
         self, merge_unit_groups, new_unit_ids, new_sorting_analyzer, keep_mask=None, verbose=False, **job_kwargs
     ):
-        # this does not depend on units
+        # this do not depend on units
         return self.data.copy()
 
     def _run(self, verbose=False):
@@ -737,15 +717,6 @@ class ComputeNoiseLevels(AnalyzerExtension):
     def _get_data(self):
         return self.data["noise_levels"]
 
-    def _handle_backward_compatibility_on_load(self):
-        # The old parameters used to be params=dict(num_chunks_per_segment=20, chunk_size=10000, seed=None)
-        # now it is handle more explicitly using random_slices_kwargs=dict()
-        for key in ("num_chunks_per_segment", "chunk_size", "seed"):
-            if key in self.params:
-                if "random_slices_kwargs" not in self.params:
-                    self.params["random_slices_kwargs"] = dict()
-                self.params["random_slices_kwargs"][key] = self.params.pop(key)
-
 
 register_result_extension(ComputeNoiseLevels)
 compute_noise_levels = ComputeNoiseLevels.function_factory()
diff --git a/src/spikeinterface/core/baserecording.py b/src/spikeinterface/core/baserecording.py
index 5e2e9e401..6d9d2a827 100644
--- a/src/spikeinterface/core/baserecording.py
+++ b/src/spikeinterface/core/baserecording.py
@@ -509,6 +509,35 @@ class BaseRecording(BaseRecordingSnippets):
             rs.t_start = None
             rs.sampling_frequency = self.sampling_frequency
 
+    def shift_times(self, shift: int | float, segment_index: int | None = None) -> None:
+        """
+        Shift all times by a scalar value.
+
+        Parameters
+        ----------
+        shift : int | float
+            The shift to apply. If positive, times will be increased by `shift`.
+            e.g. shifting by 1 will be like the recording started 1 second later.
+            If negative, the start time will be decreased i.e. as if the recording
+            started earlier.
+
+        segment_index : int | None
+            The segment on which to shift the times.
+            If `None`, all segments will be shifted.
+        """
+        if segment_index is None:
+            segments_to_shift = range(self.get_num_segments())
+        else:
+            segments_to_shift = (segment_index,)
+
+        for idx in segments_to_shift:
+            rs = self._recording_segments[idx]
+
+            if self.has_time_vector(segment_index=idx):
+                rs.time_vector += shift
+            else:
+                rs.t_start += shift
+
     def sample_index_to_time(self, sample_ind, segment_index=None):
         """
         Transform sample index into time in seconds
diff --git a/src/spikeinterface/core/baserecordingsnippets.py b/src/spikeinterface/core/baserecordingsnippets.py
index 2ec3664a4..310533c96 100644
--- a/src/spikeinterface/core/baserecordingsnippets.py
+++ b/src/spikeinterface/core/baserecordingsnippets.py
@@ -172,10 +172,8 @@ class BaseRecordingSnippets(BaseExtractor):
         number_of_device_channel_indices = np.max(list(device_channel_indices) + [0])
         if number_of_device_channel_indices >= self.get_num_channels():
             error_msg = (
-                f"The given Probe either has 'device_channel_indices' that does not match channel count \n"
-                f"{len(device_channel_indices)} vs {self.get_num_channels()} \n"
-                f"or it's max index {number_of_device_channel_indices} is the same as the number of channels {self.get_num_channels()} \n"
-                f"If using all channels remember that python is 0-indexed so max device_channel_index should be {self.get_num_channels() - 1} \n"
+                f"The given Probe have 'device_channel_indices' that do not match channel count \n"
+                f"{number_of_device_channel_indices} vs {self.get_num_channels()} \n"
                 f"device_channel_indices are the following: {device_channel_indices} \n"
                 f"recording channels are the following: {self.get_channel_ids()} \n"
             )
diff --git a/src/spikeinterface/core/job_tools.py b/src/spikeinterface/core/job_tools.py
index 7a6172369..5240edcee 100644
--- a/src/spikeinterface/core/job_tools.py
+++ b/src/spikeinterface/core/job_tools.py
@@ -149,12 +149,12 @@ def divide_segment_into_chunks(num_frames, chunk_size):
 
 
 def divide_recording_into_chunks(recording, chunk_size):
-    recording_slices = []
+    all_chunks = []
     for segment_index in range(recording.get_num_segments()):
         num_frames = recording.get_num_samples(segment_index)
         chunks = divide_segment_into_chunks(num_frames, chunk_size)
-        recording_slices.extend([(segment_index, frame_start, frame_stop) for frame_start, frame_stop in chunks])
-    return recording_slices
+        all_chunks.extend([(segment_index, frame_start, frame_stop) for frame_start, frame_stop in chunks])
+    return all_chunks
 
 
 def ensure_n_jobs(recording, n_jobs=1):
@@ -185,22 +185,6 @@ def ensure_n_jobs(recording, n_jobs=1):
     return n_jobs
 
 
-def chunk_duration_to_chunk_size(chunk_duration, recording):
-    if isinstance(chunk_duration, float):
-        chunk_size = int(chunk_duration * recording.get_sampling_frequency())
-    elif isinstance(chunk_duration, str):
-        if chunk_duration.endswith("ms"):
-            chunk_duration = float(chunk_duration.replace("ms", "")) / 1000.0
-        elif chunk_duration.endswith("s"):
-            chunk_duration = float(chunk_duration.replace("s", ""))
-        else:
-            raise ValueError("chunk_duration must ends with s or ms")
-        chunk_size = int(chunk_duration * recording.get_sampling_frequency())
-    else:
-        raise ValueError("chunk_duration must be str or float")
-    return chunk_size
-
-
 def ensure_chunk_size(
     recording, total_memory=None, chunk_size=None, chunk_memory=None, chunk_duration=None, n_jobs=1, **other_kwargs
 ):
@@ -247,7 +231,18 @@ def ensure_chunk_size(
         num_channels = recording.get_num_channels()
         chunk_size = int(total_memory / (num_channels * n_bytes * n_jobs))
     elif chunk_duration is not None:
-        chunk_size = chunk_duration_to_chunk_size(chunk_duration, recording)
+        if isinstance(chunk_duration, float):
+            chunk_size = int(chunk_duration * recording.get_sampling_frequency())
+        elif isinstance(chunk_duration, str):
+            if chunk_duration.endswith("ms"):
+                chunk_duration = float(chunk_duration.replace("ms", "")) / 1000.0
+            elif chunk_duration.endswith("s"):
+                chunk_duration = float(chunk_duration.replace("s", ""))
+            else:
+                raise ValueError("chunk_duration must ends with s or ms")
+            chunk_size = int(chunk_duration * recording.get_sampling_frequency())
+        else:
+            raise ValueError("chunk_duration must be str or float")
     else:
         # Edge case to define single chunk per segment for n_jobs=1.
         # All chunking parameters equal None mean single chunk per segment
@@ -387,13 +382,11 @@ class ChunkRecordingExecutor:
                 f"chunk_duration={chunk_duration_str}",
             )
 
-    def run(self, recording_slices=None):
+    def run(self):
         """
         Runs the defined jobs.
         """
-
-        if recording_slices is None:
-            recording_slices = divide_recording_into_chunks(self.recording, self.chunk_size)
+        all_chunks = divide_recording_into_chunks(self.recording, self.chunk_size)
 
         if self.handle_returns:
             returns = []
@@ -402,17 +395,17 @@ class ChunkRecordingExecutor:
 
         if self.n_jobs == 1:
             if self.progress_bar:
-                recording_slices = tqdm(recording_slices, ascii=True, desc=self.job_name)
+                all_chunks = tqdm(all_chunks, ascii=True, desc=self.job_name)
 
             worker_ctx = self.init_func(*self.init_args)
-            for segment_index, frame_start, frame_stop in recording_slices:
+            for segment_index, frame_start, frame_stop in all_chunks:
                 res = self.func(segment_index, frame_start, frame_stop, worker_ctx)
                 if self.handle_returns:
                     returns.append(res)
                 if self.gather_func is not None:
                     self.gather_func(res)
         else:
-            n_jobs = min(self.n_jobs, len(recording_slices))
+            n_jobs = min(self.n_jobs, len(all_chunks))
 
             # parallel
             with ProcessPoolExecutor(
@@ -421,10 +414,10 @@ class ChunkRecordingExecutor:
                 mp_context=mp.get_context(self.mp_context),
                 initargs=(self.func, self.init_func, self.init_args, self.max_threads_per_process),
             ) as executor:
-                results = executor.map(function_wrapper, recording_slices)
+                results = executor.map(function_wrapper, all_chunks)
 
                 if self.progress_bar:
-                    results = tqdm(results, desc=self.job_name, total=len(recording_slices))
+                    results = tqdm(results, desc=self.job_name, total=len(all_chunks))
 
                 for res in results:
                     if self.handle_returns:
diff --git a/src/spikeinterface/core/node_pipeline.py b/src/spikeinterface/core/node_pipeline.py
index 53c2445c7..d90a20902 100644
--- a/src/spikeinterface/core/node_pipeline.py
+++ b/src/spikeinterface/core/node_pipeline.py
@@ -489,7 +489,6 @@ def run_node_pipeline(
     names=None,
     verbose=False,
     skip_after_n_peaks=None,
-    recording_slices=None,
 ):
     """
     Machinery to compute in parallel operations on peaks and traces.
@@ -541,10 +540,6 @@ def run_node_pipeline(
     skip_after_n_peaks : None | int
         Skip the computation after n_peaks.
         This is not an exact because internally this skip is done per worker in average.
-    recording_slices : None | list[tuple]
-        Optionaly give a list of slices to run the pipeline only on some chunks of the recording.
-        It must be a list of (segment_index, frame_start, frame_stop).
-        If None (default), the function iterates over the entire duration of the recording.
 
     Returns
     -------
@@ -583,7 +578,7 @@ def run_node_pipeline(
         **job_kwargs,
     )
 
-    processor.run(recording_slices=recording_slices)
+    processor.run()
 
     outs = gather_func.finalize_buffers(squeeze_output=squeeze_output)
     return outs
diff --git a/src/spikeinterface/core/recording_tools.py b/src/spikeinterface/core/recording_tools.py
index 4aabbfd58..77d427bc8 100644
--- a/src/spikeinterface/core/recording_tools.py
+++ b/src/spikeinterface/core/recording_tools.py
@@ -18,8 +18,6 @@ from .job_tools import (
     fix_job_kwargs,
     ChunkRecordingExecutor,
     _shared_job_kwargs_doc,
-    chunk_duration_to_chunk_size,
-    split_job_kwargs,
 )
 
 
@@ -514,38 +512,33 @@ def determine_cast_unsigned(recording, dtype):
     return cast_unsigned
 
 
-def get_random_recording_slices(
+def get_random_data_chunks(
     recording,
-    method="full_random",
+    return_scaled=False,
     num_chunks_per_segment=20,
-    chunk_duration="500ms",
-    chunk_size=None,
+    chunk_size=10000,
+    concatenated=True,
+    seed=0,
     margin_frames=0,
-    seed=None,
 ):
     """
-    Get random slice of a recording across segments.
+    Extract random chunks across segments
 
-    This is used for instance in get_noise_levels() and get_random_data_chunks() to estimate noise on traces.
+    This is used for instance in get_noise_levels() to estimate noise on traces.
 
     Parameters
     ----------
     recording : BaseRecording
         The recording to get random chunks from
-    method : "full_random"
-        The method used to get random slices.
-          * "full_random" : legacy method,  used until version 0.101.0, there is no constrain on slices
-            and they can overlap.
+    return_scaled : bool, default: False
+        If True, returned chunks are scaled to uV
     num_chunks_per_segment : int, default: 20
         Number of chunks per segment
-    chunk_duration : str | float | None, default "500ms"
-        The duration of each chunk in 's' or 'ms'
-    chunk_size : int | None
-        Size of a chunk in number of frames. This is ued only if chunk_duration is None.
-        This is kept for backward compatibility, you should prefer 'chunk_duration=500ms' instead.
+    chunk_size : int, default: 10000
+        Size of a chunk in number of frames
     concatenated : bool, default: True
         If True chunk are concatenated along time axis
-    seed : int, default: None
+    seed : int, default: 0
         Random seed
     margin_frames : int, default: 0
         Margin in number of frames to avoid edge effects
@@ -554,89 +547,42 @@ def get_random_recording_slices(
     -------
     chunk_list : np.array
         Array of concatenate chunks per segment
-
-
     """
     # TODO: if segment have differents length make another sampling that dependant on the length of the segment
     # Should be done by changing kwargs with total_num_chunks=XXX and total_duration=YYYY
     # And randomize the number of chunk per segment weighted by segment duration
 
-    if method == "full_random":
-        if chunk_size is None:
-            if chunk_duration is not None:
-                chunk_size = chunk_duration_to_chunk_size(chunk_duration, recording)
-            else:
-                raise ValueError("get_random_recording_slices need chunk_size or chunk_duration")
-
-        # check chunk size
-        num_segments = recording.get_num_segments()
-        for segment_index in range(num_segments):
-            chunk_size_limit = recording.get_num_frames(segment_index) - 2 * margin_frames
-            if chunk_size > chunk_size_limit:
-                chunk_size = chunk_size_limit - 1
-                warnings.warn(
-                    f"chunk_size is greater than the number "
-                    f"of samples for segment index {segment_index}. "
-                    f"Using {chunk_size}."
-                )
-        rng = np.random.default_rng(seed)
-        recording_slices = []
-        low = margin_frames
-        size = num_chunks_per_segment
-        for segment_index in range(num_segments):
-            num_frames = recording.get_num_frames(segment_index)
-            high = num_frames - chunk_size - margin_frames
-            random_starts = rng.integers(low=low, high=high, size=size)
-            random_starts = np.sort(random_starts)
-            recording_slices += [
-                (segment_index, start_frame, (start_frame + chunk_size)) for start_frame in random_starts
-            ]
-    else:
-        raise ValueError(f"get_random_recording_slices : wrong method {method}")
-
-    return recording_slices
-
-
-def get_random_data_chunks(recording, return_scaled=False, concatenated=True, **random_slices_kwargs):
-    """
-    Extract random chunks across segments.
-
-    Internally, it uses `get_random_recording_slices()` and retrieves the traces chunk as a list
-    or a concatenated unique array.
-
-    Please read `get_random_recording_slices()` for more details on parameters.
-
-
-    Parameters
-    ----------
-    recording : BaseRecording
-        The recording to get random chunks from
-    return_scaled : bool, default: False
-        If True, returned chunks are scaled to uV
-    num_chunks_per_segment : int, default: 20
-        Number of chunks per segment
-    concatenated : bool, default: True
-        If True chunk are concatenated along time axis
-    **random_slices_kwargs : dict
-        Options transmited to  get_random_recording_slices(), please read documentation from this
-        function for more details.
-
-    Returns
-    -------
-    chunk_list : np.array | list of np.array
-        Array of concatenate chunks per segment
-    """
-    recording_slices = get_random_recording_slices(recording, **random_slices_kwargs)
+    # check chunk size
+    num_segments = recording.get_num_segments()
+    for segment_index in range(num_segments):
+        chunk_size_limit = recording.get_num_frames(segment_index) - 2 * margin_frames
+        if chunk_size > chunk_size_limit:
+            chunk_size = chunk_size_limit - 1
+            warnings.warn(
+                f"chunk_size is greater than the number "
+                f"of samples for segment index {segment_index}. "
+                f"Using {chunk_size}."
+            )
 
+    rng = np.random.default_rng(seed)
     chunk_list = []
-    for segment_index, start_frame, end_frame in recording_slices:
-        traces_chunk = recording.get_traces(
-            start_frame=start_frame,
-            end_frame=end_frame,
-            segment_index=segment_index,
-            return_scaled=return_scaled,
-        )
-        chunk_list.append(traces_chunk)
+    low = margin_frames
+    size = num_chunks_per_segment
+    for segment_index in range(num_segments):
+        num_frames = recording.get_num_frames(segment_index)
+        high = num_frames - chunk_size - margin_frames
+        random_starts = rng.integers(low=low, high=high, size=size)
+        segment_trace_chunk = [
+            recording.get_traces(
+                start_frame=start_frame,
+                end_frame=(start_frame + chunk_size),
+                segment_index=segment_index,
+                return_scaled=return_scaled,
+            )
+            for start_frame in random_starts
+        ]
+
+        chunk_list.extend(segment_trace_chunk)
 
     if concatenated:
         return np.concatenate(chunk_list, axis=0)
@@ -691,52 +637,19 @@ def get_closest_channels(recording, channel_ids=None, num_channels=None):
     return np.array(closest_channels_inds), np.array(dists)
 
 
-def _noise_level_chunk(segment_index, start_frame, end_frame, worker_ctx):
-    recording = worker_ctx["recording"]
-
-    one_chunk = recording.get_traces(
-        start_frame=start_frame,
-        end_frame=end_frame,
-        segment_index=segment_index,
-        return_scaled=worker_ctx["return_scaled"],
-    )
-
-    if worker_ctx["method"] == "mad":
-        med = np.median(one_chunk, axis=0, keepdims=True)
-        # hard-coded so that core doesn't depend on scipy
-        noise_levels = np.median(np.abs(one_chunk - med), axis=0) / 0.6744897501960817
-    elif worker_ctx["method"] == "std":
-        noise_levels = np.std(one_chunk, axis=0)
-
-    return noise_levels
-
-
-def _noise_level_chunk_init(recording, return_scaled, method):
-    worker_ctx = {}
-    worker_ctx["recording"] = recording
-    worker_ctx["return_scaled"] = return_scaled
-    worker_ctx["method"] = method
-    return worker_ctx
-
-
 def get_noise_levels(
     recording: "BaseRecording",
     return_scaled: bool = True,
     method: Literal["mad", "std"] = "mad",
     force_recompute: bool = False,
-    random_slices_kwargs: dict = {},
-    **kwargs,
+    **random_chunk_kwargs,
 ) -> np.ndarray:
     """
     Estimate noise for each channel using MAD methods.
     You can use standard deviation with `method="std"`
 
     Internally it samples some chunk across segment.
-    And then, it uses the MAD estimator (more robust than STD) or the STD on each chunk.
-    Finally the average of all MAD/STD values is performed.
-
-    The result is cached in a property of the recording, so that the next call on the same
-    recording will use the cached result unless `force_recompute=True`.
+    And then, it use MAD estimator (more robust than STD)
 
     Parameters
     ----------
@@ -749,11 +662,8 @@ def get_noise_levels(
         The method to use to estimate noise levels
     force_recompute : bool
         If True, noise levels are recomputed even if they are already stored in the recording extractor
-    random_slices_kwargs : dict
-        Options transmited to  get_random_recording_slices(), please read documentation from this
-        function for more details.
-
-    {}
+    random_chunk_kwargs : dict
+        Kwargs for get_random_data_chunks
 
     Returns
     -------
@@ -769,56 +679,19 @@ def get_noise_levels(
     if key in recording.get_property_keys() and not force_recompute:
         noise_levels = recording.get_property(key=key)
     else:
-        # This is to keep backward compatibility
-        # lets keep for a while and remove this maybe in 0.103.0
-        # chunk_size used to be in the signature and now is ambiguous
-        random_slices_kwargs_, job_kwargs = split_job_kwargs(kwargs)
-        if len(random_slices_kwargs_) > 0 or "chunk_size" in job_kwargs:
-            msg = (
-                "get_noise_levels(recording, num_chunks_per_segment=20) is deprecated\n"
-                "Now, you need to use get_noise_levels(recording, random_slices_kwargs=dict(num_chunks_per_segment=20, chunk_size=1000))\n"
-                "Please read get_random_recording_slices() documentation for more options."
-            )
-            # if the user use both the old and the new behavior then an error is raised
-            assert len(random_slices_kwargs) == 0, msg
-            warnings.warn(msg)
-            random_slices_kwargs = random_slices_kwargs_
-            if "chunk_size" in job_kwargs:
-                random_slices_kwargs["chunk_size"] = job_kwargs["chunk_size"]
-
-        recording_slices = get_random_recording_slices(recording, **random_slices_kwargs)
-
-        noise_levels_chunks = []
-
-        def append_noise_chunk(res):
-            noise_levels_chunks.append(res)
-
-        func = _noise_level_chunk
-        init_func = _noise_level_chunk_init
-        init_args = (recording, return_scaled, method)
-        executor = ChunkRecordingExecutor(
-            recording,
-            func,
-            init_func,
-            init_args,
-            job_name="noise_level",
-            verbose=False,
-            gather_func=append_noise_chunk,
-            **job_kwargs,
-        )
-        executor.run(recording_slices=recording_slices)
-        noise_levels_chunks = np.stack(noise_levels_chunks)
-        noise_levels = np.mean(noise_levels_chunks, axis=0)
-
-        # set property
+        random_chunks = get_random_data_chunks(recording, return_scaled=return_scaled, **random_chunk_kwargs)
+
+        if method == "mad":
+            med = np.median(random_chunks, axis=0, keepdims=True)
+            # hard-coded so that core doesn't depend on scipy
+            noise_levels = np.median(np.abs(random_chunks - med), axis=0) / 0.6744897501960817
+        elif method == "std":
+            noise_levels = np.std(random_chunks, axis=0)
         recording.set_property(key, noise_levels)
 
     return noise_levels
 
 
-get_noise_levels.__doc__ = get_noise_levels.__doc__.format(_shared_job_kwargs_doc)
-
-
 def get_chunk_with_margin(
     rec_segment,
     start_frame,
diff --git a/src/spikeinterface/core/template_tools.py b/src/spikeinterface/core/template_tools.py
index 3c8663df7..934b18ed4 100644
--- a/src/spikeinterface/core/template_tools.py
+++ b/src/spikeinterface/core/template_tools.py
@@ -31,12 +31,7 @@ def get_dense_templates_array(one_object: Templates | SortingAnalyzer, return_sc
             )
         ext = one_object.get_extension("templates")
         if ext is not None:
-            if "average" in ext.data:
-                templates_array = ext.data.get("average")
-            elif "median" in ext.data:
-                templates_array = ext.data.get("median")
-            else:
-                raise ValueError("Average or median templates have not been computed.")
+            templates_array = ext.data["average"]
         else:
             raise ValueError("SortingAnalyzer need extension 'templates' to be computed to retrieve templates")
     else:
diff --git a/src/spikeinterface/curation/__init__.py b/src/spikeinterface/curation/__init__.py
index 0302ffe5b..657b936fb 100644
--- a/src/spikeinterface/curation/__init__.py
+++ b/src/spikeinterface/curation/__init__.py
@@ -3,7 +3,7 @@ from .curation_tools import find_duplicated_spikes
 from .remove_redundant import remove_redundant_units, find_redundant_units
 from .remove_duplicated_spikes import remove_duplicated_spikes
 from .remove_excess_spikes import remove_excess_spikes
-from .auto_merge import compute_merge_unit_groups, auto_merge_units, get_potential_auto_merge
+from .auto_merge import get_potential_auto_merge
 
 
 # manual sorting,
diff --git a/src/spikeinterface/curation/auto_merge.py b/src/spikeinterface/curation/auto_merge.py
index 4f4cff144..19336e594 100644
--- a/src/spikeinterface/curation/auto_merge.py
+++ b/src/spikeinterface/curation/auto_merge.py
@@ -1,7 +1,5 @@
 from __future__ import annotations
 
-import warnings
-
 from typing import Tuple
 import numpy as np
 import math
@@ -14,90 +12,49 @@ except ImportError:
     HAVE_NUMBA = False
 
 from ..core import SortingAnalyzer, Templates
+from ..core.template_tools import get_template_extremum_channel
+from ..postprocessing import compute_correlograms
 from ..qualitymetrics import compute_refrac_period_violations, compute_firing_rates
 
 from .mergeunitssorting import MergeUnitsSorting
 from .curation_tools import resolve_merging_graph
 
-_compute_merge_presets = {
-    "similarity_correlograms": [
-        "num_spikes",
-        "remove_contaminated",
-        "unit_locations",
-        "template_similarity",
-        "correlogram",
-        "quality_score",
-    ],
-    "temporal_splits": [
-        "num_spikes",
-        "remove_contaminated",
-        "unit_locations",
-        "template_similarity",
-        "presence_distance",
-        "quality_score",
-    ],
-    "x_contaminations": [
-        "num_spikes",
-        "remove_contaminated",
-        "unit_locations",
-        "template_similarity",
-        "cross_contamination",
-        "quality_score",
-    ],
-    "feature_neighbors": [
-        "num_spikes",
-        "snr",
-        "remove_contaminated",
-        "unit_locations",
-        "knn",
-        "quality_score",
-    ],
-}
+
+_possible_presets = ["similarity_correlograms", "x_contaminations", "temporal_splits", "feature_neighbors"]
 
 _required_extensions = {
-    "unit_locations": ["templates", "unit_locations"],
+    "unit_locations": ["unit_locations"],
     "correlogram": ["correlograms"],
-    "snr": ["templates", "noise_levels"],
-    "template_similarity": ["templates", "template_similarity"],
-    "knn": ["templates", "spike_locations", "spike_amplitudes"],
-}
-
-
-_default_step_params = {
-    "num_spikes": {"min_spikes": 100},
-    "snr": {"min_snr": 2},
-    "remove_contaminated": {"contamination_thresh": 0.2, "refractory_period_ms": 1.0, "censored_period_ms": 0.3},
-    "unit_locations": {"max_distance_um": 150},
-    "correlogram": {
-        "corr_diff_thresh": 0.16,
-        "censor_correlograms_ms": 0.15,
-        "sigma_smooth_ms": 0.6,
-        "adaptative_window_thresh": 0.5,
-    },
-    "template_similarity": {"template_diff_thresh": 0.25},
-    "presence_distance": {"presence_distance_thresh": 100},
-    "knn": {"k_nn": 10},
-    "cross_contamination": {
-        "cc_thresh": 0.1,
-        "p_value": 0.2,
-        "refractory_period_ms": 1.0,
-        "censored_period_ms": 0.3,
-    },
-    "quality_score": {"firing_contamination_balance": 1.5, "refractory_period_ms": 1.0, "censored_period_ms": 0.3},
+    "template_similarity": ["template_similarity"],
+    "knn": ["spike_locations", "spike_amplitudes"],
 }
 
 
-def compute_merge_unit_groups(
+def get_potential_auto_merge(
     sorting_analyzer: SortingAnalyzer,
     preset: str | None = "similarity_correlograms",
-    resolve_graph: bool = True,
-    steps_params: dict = None,
-    compute_needed_extensions: bool = True,
+    resolve_graph: bool = False,
+    min_spikes: int = 100,
+    min_snr: float = 2,
+    max_distance_um: float = 150.0,
+    corr_diff_thresh: float = 0.16,
+    template_diff_thresh: float = 0.25,
+    contamination_thresh: float = 0.2,
+    presence_distance_thresh: float = 100,
+    p_value: float = 0.2,
+    cc_thresh: float = 0.1,
+    censored_period_ms: float = 0.3,
+    refractory_period_ms: float = 1.0,
+    sigma_smooth_ms: float = 0.6,
+    adaptative_window_thresh: float = 0.5,
+    censor_correlograms_ms: float = 0.15,
+    firing_contamination_balance: float = 2.5,
+    k_nn: int = 10,
+    knn_kwargs: dict | None = None,
+    presence_distance_kwargs: dict | None = None,
     extra_outputs: bool = False,
     steps: list[str] | None = None,
-    force_copy: bool = True,
-    **job_kwargs,
-) -> list[tuple[int | str, int | str]] | Tuple[list[tuple[int | str, int | str]], dict]:
+) -> list[tuple[int | str, int | str]] | Tuple[tuple[int | str, int | str], dict]:
     """
     Algorithm to find and check potential merges between units.
 
@@ -121,9 +78,6 @@ def compute_merge_unit_groups(
 
         Q = f(1 - (k + 1)C)
 
-    IMPORTANT: internally, all computations are relying on extensions of the analyzer, that are computed
-    with default parameters if not present (i.e. correlograms, template_similarity, ...) If you want to
-    have a finer control on these values, please precompute the extensions before applying the auto_merge
 
     Parameters
     ----------
@@ -144,11 +98,47 @@ def compute_merge_unit_groups(
         * | "feature_neighbors": focused on finding unit pairs whose spikes are close in the feature space using kNN.
           | It uses the following steps: "num_spikes", "snr", "remove_contaminated", "unit_locations",
           | "knn", "quality_score"
+
         If `preset` is None, you can specify the steps manually with the `steps` parameter.
-    resolve_graph : bool, default: True
+    resolve_graph : bool, default: False
         If True, the function resolves the potential unit pairs to be merged into multiple-unit merges.
-    compute_needed_extensions : bool, default : True
-        Should we force the computation of needed extensions, if not already computed?
+    min_spikes : int, default: 100
+        Minimum number of spikes for each unit to consider a potential merge.
+        Enough spikes are needed to estimate the correlogram
+    min_snr : float, default 2
+        Minimum Signal to Noise ratio for templates to be considered while merging
+    max_distance_um : float, default: 150
+        Maximum distance between units for considering a merge
+    corr_diff_thresh : float, default: 0.16
+        The threshold on the "correlogram distance metric" for considering a merge.
+        It needs to be between 0 and 1
+    template_diff_thresh : float, default: 0.25
+        The threshold on the "template distance metric" for considering a merge.
+        It needs to be between 0 and 1
+    contamination_thresh : float, default: 0.2
+        Threshold for not taking in account a unit when it is too contaminated.
+    presence_distance_thresh : float, default: 100
+        Parameter to control how present two units should be simultaneously.
+    p_value : float, default: 0.2
+        The p-value threshold for the cross-contamination test.
+    cc_thresh : float, default: 0.1
+        The threshold on the cross-contamination for considering a merge.
+    censored_period_ms : float, default: 0.3
+        Used to compute the refractory period violations aka "contamination".
+    refractory_period_ms : float, default: 1
+        Used to compute the refractory period violations aka "contamination".
+    sigma_smooth_ms : float, default: 0.6
+        Parameters to smooth the correlogram estimation.
+    adaptative_window_thresh : float, default: 0.5
+        Parameter to detect the window size in correlogram estimation.
+    censor_correlograms_ms : float, default: 0.15
+        The period to censor on the auto and cross-correlograms.
+    firing_contamination_balance : float, default: 2.5
+        Parameter to control the balance between firing rate and contamination in computing unit "quality score".
+    k_nn : int, default 5
+        The number of neighbors to consider for every spike in the recording.
+    knn_kwargs : dict, default None
+        The dict of extra params to be passed to knn.
     extra_outputs : bool, default: False
         If True, an additional dictionary (`outs`) with processed data is returned.
     steps : None or list of str, default: None
@@ -156,141 +146,157 @@ def compute_merge_unit_groups(
         Pontential steps : "num_spikes", "snr", "remove_contaminated", "unit_locations", "correlogram",
         "template_similarity", "presence_distance", "cross_contamination", "knn", "quality_score"
         Please check steps explanations above!
-    steps_params : dict
-        A dictionary whose keys are the steps, and keys are steps parameters.
-    force_copy : boolean, default: True
-        When new extensions are computed, the default is to make a copy of the analyzer, to avoid overwriting
-        already computed extensions. False if you want to overwrite
+    presence_distance_kwargs : None|dict, default: None
+        A dictionary of kwargs to be passed to compute_presence_distance().
 
     Returns
     -------
-    merge_unit_groups:
-        List of groups that need to be merge.
-        When `resolve_graph` is true (default) a list of tuples of 2+ elements
-        If `resolve_graph` is false then a list of tuple of 2 elements is returned instead.
+    potential_merges:
+        A list of tuples of 2 elements (if `resolve_graph`if false) or 2+ elements (if `resolve_graph` is true).
+        List of pairs that could be merged.
     outs:
         Returned only when extra_outputs=True
         A dictionary that contains data for debugging and plotting.
 
     References
     ----------
-    This function used to be inspired and built upon similar functions from Lussac [Llobet]_,
+    This function is inspired and built upon similar functions from Lussac [Llobet]_,
     done by Aurelien Wyngaard and Victor Llobet.
     https://github.com/BarbourLab/lussac/blob/v1.0.0/postprocessing/merge_units.py
-
-    However, it has been greatly consolidated and refined depending on the presets.
     """
     import scipy
 
     sorting = sorting_analyzer.sorting
     unit_ids = sorting.unit_ids
 
-    if preset is None and steps is None:
-        raise ValueError("You need to specify a preset or steps for the auto-merge function")
-    elif steps is not None:
-        # steps has precedence on presets
-        pass
-    elif preset is not None:
-        if preset not in _compute_merge_presets:
-            raise ValueError(f"preset must be one of {list(_compute_merge_presets.keys())}")
-        steps = _compute_merge_presets[preset]
-
-    # check at least one extension is needed
-    at_least_one_extension_to_compute = False
+    # to get fast computation we will not analyse pairs when:
+    #    * not enough spikes for one of theses
+    #    * auto correlogram is contaminated
+    #    * to far away one from each other
+
+    all_steps = [
+        "num_spikes",
+        "snr",
+        "remove_contaminated",
+        "unit_locations",
+        "correlogram",
+        "template_similarity",
+        "presence_distance",
+        "knn",
+        "cross_contamination",
+        "quality_score",
+    ]
+
+    if preset is not None and preset not in _possible_presets:
+        raise ValueError(f"preset must be one of {_possible_presets}")
+
+    if steps is None:
+        if preset is None:
+            if steps is None:
+                raise ValueError("You need to specify a preset or steps for the auto-merge function")
+        elif preset == "similarity_correlograms":
+            steps = [
+                "num_spikes",
+                "remove_contaminated",
+                "unit_locations",
+                "template_similarity",
+                "correlogram",
+                "quality_score",
+            ]
+        elif preset == "temporal_splits":
+            steps = [
+                "num_spikes",
+                "remove_contaminated",
+                "unit_locations",
+                "template_similarity",
+                "presence_distance",
+                "quality_score",
+            ]
+        elif preset == "x_contaminations":
+            steps = [
+                "num_spikes",
+                "remove_contaminated",
+                "unit_locations",
+                "template_similarity",
+                "cross_contamination",
+                "quality_score",
+            ]
+        elif preset == "feature_neighbors":
+            steps = [
+                "num_spikes",
+                "snr",
+                "remove_contaminated",
+                "unit_locations",
+                "knn",
+                "quality_score",
+            ]
+
     for step in steps:
-        assert step in _default_step_params, f"{step} is not a valid step"
         if step in _required_extensions:
             for ext in _required_extensions[step]:
-                if sorting_analyzer.has_extension(ext):
-                    continue
-                if not compute_needed_extensions:
+                if not sorting_analyzer.has_extension(ext):
                     raise ValueError(f"{step} requires {ext} extension")
-                at_least_one_extension_to_compute = True
-
-    if force_copy and at_least_one_extension_to_compute:
-        # To avoid erasing the extensions of the user
-        sorting_analyzer = sorting_analyzer.copy()
 
     n = unit_ids.size
-    pair_mask = np.triu(np.arange(n), 1) > 0
+    pair_mask = np.triu(np.arange(n)) > 0
     outs = dict()
 
     for step in steps:
 
-        if step in _required_extensions:
-            for ext in _required_extensions[step]:
-                if sorting_analyzer.has_extension(ext):
-                    continue
-
-                # special case for templates
-                if ext == "templates" and not sorting_analyzer.has_extension("random_spikes"):
-                    sorting_analyzer.compute(["random_spikes", "templates"], **job_kwargs)
-                else:
-                    sorting_analyzer.compute(ext, **job_kwargs)
-
-        params = _default_step_params.get(step).copy()
-        if steps_params is not None and step in steps_params:
-            params.update(steps_params[step])
+        assert step in all_steps, f"{step} is not a valid step"
 
         # STEP : remove units with too few spikes
         if step == "num_spikes":
-
             num_spikes = sorting.count_num_spikes_per_unit(outputs="array")
-            to_remove = num_spikes < params["min_spikes"]
+            to_remove = num_spikes < min_spikes
             pair_mask[to_remove, :] = False
             pair_mask[:, to_remove] = False
-            outs["num_spikes"] = to_remove
 
         # STEP : remove units with too small SNR
         elif step == "snr":
             qm_ext = sorting_analyzer.get_extension("quality_metrics")
             if qm_ext is None:
-                sorting_analyzer.compute("quality_metrics", metric_names=["snr"], **job_kwargs)
+                sorting_analyzer.compute("noise_levels")
+                sorting_analyzer.compute("quality_metrics", metric_names=["snr"])
                 qm_ext = sorting_analyzer.get_extension("quality_metrics")
 
             snrs = qm_ext.get_data()["snr"].values
-            to_remove = snrs < params["min_snr"]
+            to_remove = snrs < min_snr
             pair_mask[to_remove, :] = False
             pair_mask[:, to_remove] = False
-            outs["snr"] = to_remove
 
         # STEP : remove contaminated auto corr
         elif step == "remove_contaminated":
             contaminations, nb_violations = compute_refrac_period_violations(
-                sorting_analyzer,
-                refractory_period_ms=params["refractory_period_ms"],
-                censored_period_ms=params["censored_period_ms"],
+                sorting_analyzer, refractory_period_ms=refractory_period_ms, censored_period_ms=censored_period_ms
             )
             nb_violations = np.array(list(nb_violations.values()))
             contaminations = np.array(list(contaminations.values()))
-            to_remove = contaminations > params["contamination_thresh"]
+            to_remove = contaminations > contamination_thresh
             pair_mask[to_remove, :] = False
             pair_mask[:, to_remove] = False
-            outs["remove_contaminated"] = to_remove
 
         # STEP : unit positions are estimated roughly with channel
-        elif step == "unit_locations":
+        elif step == "unit_locations" in steps:
             location_ext = sorting_analyzer.get_extension("unit_locations")
             unit_locations = location_ext.get_data()[:, :2]
 
             unit_distances = scipy.spatial.distance.cdist(unit_locations, unit_locations, metric="euclidean")
-            pair_mask = pair_mask & (unit_distances <= params["max_distance_um"])
+            pair_mask = pair_mask & (unit_distances <= max_distance_um)
             outs["unit_distances"] = unit_distances
 
         # STEP : potential auto merge by correlogram
-        elif step == "correlogram":
+        elif step == "correlogram" in steps:
             correlograms_ext = sorting_analyzer.get_extension("correlograms")
             correlograms, bins = correlograms_ext.get_data()
-            censor_ms = params["censor_correlograms_ms"]
-            sigma_smooth_ms = params["sigma_smooth_ms"]
-            mask = (bins[:-1] >= -censor_ms) & (bins[:-1] < censor_ms)
+            mask = (bins[:-1] >= -censor_correlograms_ms) & (bins[:-1] < censor_correlograms_ms)
             correlograms[:, :, mask] = 0
             correlograms_smoothed = smooth_correlogram(correlograms, bins, sigma_smooth_ms=sigma_smooth_ms)
             # find correlogram window for each units
             win_sizes = np.zeros(n, dtype=int)
             for unit_ind in range(n):
                 auto_corr = correlograms_smoothed[unit_ind, unit_ind, :]
-                thresh = np.max(auto_corr) * params["adaptative_window_thresh"]
+                thresh = np.max(auto_corr) * adaptative_window_thresh
                 win_size = get_unit_adaptive_window(auto_corr, thresh)
                 win_sizes[unit_ind] = win_size
             correlogram_diff = compute_correlogram_diff(
@@ -300,7 +306,7 @@ def compute_merge_unit_groups(
                 pair_mask=pair_mask,
             )
             # print(correlogram_diff)
-            pair_mask = pair_mask & (correlogram_diff < params["corr_diff_thresh"])
+            pair_mask = pair_mask & (correlogram_diff < corr_diff_thresh)
             outs["correlograms"] = correlograms
             outs["bins"] = bins
             outs["correlograms_smoothed"] = correlograms_smoothed
@@ -308,21 +314,22 @@ def compute_merge_unit_groups(
             outs["win_sizes"] = win_sizes
 
         # STEP : check if potential merge with CC also have template similarity
-        elif step == "template_similarity":
+        elif step == "template_similarity" in steps:
             template_similarity_ext = sorting_analyzer.get_extension("template_similarity")
             templates_similarity = template_similarity_ext.get_data()
             templates_diff = 1 - templates_similarity
-            pair_mask = pair_mask & (templates_diff < params["template_diff_thresh"])
+            pair_mask = pair_mask & (templates_diff < template_diff_thresh)
             outs["templates_diff"] = templates_diff
 
         # STEP : check the vicinity of the spikes
-        elif step == "knn":
-            pair_mask = get_pairs_via_nntree(sorting_analyzer, **params, pair_mask=pair_mask)
+        elif step == "knn" in steps:
+            if knn_kwargs is None:
+                knn_kwargs = dict()
+            pair_mask = get_pairs_via_nntree(sorting_analyzer, k_nn, pair_mask, **knn_kwargs)
 
         # STEP : check how the rates overlap in times
-        elif step == "presence_distance":
-            presence_distance_kwargs = params.copy()
-            presence_distance_thresh = presence_distance_kwargs.pop("presence_distance_thresh")
+        elif step == "presence_distance" in steps:
+            presence_distance_kwargs = presence_distance_kwargs or dict()
             num_samples = [
                 sorting_analyzer.get_num_samples(segment_index) for segment_index in range(sorting.get_num_segments())
             ]
@@ -333,243 +340,40 @@ def compute_merge_unit_groups(
             outs["presence_distances"] = presence_distances
 
         # STEP : check if the cross contamination is significant
-        elif step == "cross_contamination":
-            refractory = (
-                params["censored_period_ms"],
-                params["refractory_period_ms"],
-            )
+        elif step == "cross_contamination" in steps:
+            refractory = (censored_period_ms, refractory_period_ms)
             CC, p_values = compute_cross_contaminations(
-                sorting_analyzer, pair_mask, params["cc_thresh"], refractory, contaminations
+                sorting_analyzer, pair_mask, cc_thresh, refractory, contaminations
             )
-            pair_mask = pair_mask & (p_values > params["p_value"])
+            pair_mask = pair_mask & (p_values > p_value)
             outs["cross_contaminations"] = CC, p_values
 
         # STEP : validate the potential merges with CC increase the contamination quality metrics
-        elif step == "quality_score":
+        elif step == "quality_score" in steps:
             pair_mask, pairs_decreased_score = check_improve_contaminations_score(
                 sorting_analyzer,
                 pair_mask,
                 contaminations,
-                params["firing_contamination_balance"],
-                params["refractory_period_ms"],
-                params["censored_period_ms"],
+                firing_contamination_balance,
+                refractory_period_ms,
+                censored_period_ms,
             )
             outs["pairs_decreased_score"] = pairs_decreased_score
 
     # FINAL STEP : create the final list from pair_mask boolean matrix
     ind1, ind2 = np.nonzero(pair_mask)
-    merge_unit_groups = list(zip(unit_ids[ind1], unit_ids[ind2]))
+    potential_merges = list(zip(unit_ids[ind1], unit_ids[ind2]))
+
+    # some methods return identities ie (1,1) which we can cleanup first.
+    potential_merges = [(ids[0], ids[1]) for ids in potential_merges if ids[0] != ids[1]]
 
     if resolve_graph:
-        merge_unit_groups = resolve_merging_graph(sorting, merge_unit_groups)
+        potential_merges = resolve_merging_graph(sorting, potential_merges)
 
     if extra_outputs:
-        return merge_unit_groups, outs
+        return potential_merges, outs
     else:
-        return merge_unit_groups
-
-
-def auto_merge_units(
-    sorting_analyzer: SortingAnalyzer, compute_merge_kwargs: dict = {}, apply_merge_kwargs: dict = {}, **job_kwargs
-) -> SortingAnalyzer:
-    """
-    Compute merge unit groups and apply it on a SortingAnalyzer.
-    Internally uses `compute_merge_unit_groups()`
-    """
-    merge_unit_groups = compute_merge_unit_groups(
-        sorting_analyzer, extra_outputs=False, **compute_merge_kwargs, **job_kwargs
-    )
-
-    merged_analyzer = sorting_analyzer.merge_units(merge_unit_groups, **apply_merge_kwargs, **job_kwargs)
-    return merged_analyzer
-
-
-def get_potential_auto_merge(
-    sorting_analyzer: SortingAnalyzer,
-    preset: str | None = "similarity_correlograms",
-    resolve_graph: bool = False,
-    min_spikes: int = 100,
-    min_snr: float = 2,
-    max_distance_um: float = 150.0,
-    corr_diff_thresh: float = 0.16,
-    template_diff_thresh: float = 0.25,
-    contamination_thresh: float = 0.2,
-    presence_distance_thresh: float = 100,
-    p_value: float = 0.2,
-    cc_thresh: float = 0.1,
-    censored_period_ms: float = 0.3,
-    refractory_period_ms: float = 1.0,
-    sigma_smooth_ms: float = 0.6,
-    adaptative_window_thresh: float = 0.5,
-    censor_correlograms_ms: float = 0.15,
-    firing_contamination_balance: float = 1.5,
-    k_nn: int = 10,
-    knn_kwargs: dict | None = None,
-    presence_distance_kwargs: dict | None = None,
-    extra_outputs: bool = False,
-    steps: list[str] | None = None,
-) -> list[tuple[int | str, int | str]] | Tuple[tuple[int | str, int | str], dict]:
-    """
-    This function is deprecated. Use compute_merge_unit_groups() instead.
-    This will be removed in 0.103.0
-
-    Algorithm to find and check potential merges between units.
-
-    The merges are proposed based on a series of steps with different criteria:
-
-        * "num_spikes": enough spikes are found in each unit for computing the correlogram (`min_spikes`)
-        * "snr": the SNR of the units is above a threshold (`min_snr`)
-        * "remove_contaminated": each unit is not contaminated (by checking auto-correlogram - `contamination_thresh`)
-        * "unit_locations": estimated unit locations are close enough (`max_distance_um`)
-        * "correlogram": the cross-correlograms of the two units are similar to each auto-corrleogram (`corr_diff_thresh`)
-        * "template_similarity": the templates of the two units are similar (`template_diff_thresh`)
-        * "presence_distance": the presence of the units is complementary in time (`presence_distance_thresh`)
-        * "cross_contamination": the cross-contamination is not significant (`cc_thresh` and `p_value`)
-        * "knn": the two units are close in the feature space
-        * "quality_score": the unit "quality score" is increased after the merge
-
-    The "quality score" factors in the increase in firing rate (**f**) due to the merge and a possible increase in
-    contamination (**C**), wheighted by a factor **k** (`firing_contamination_balance`).
-
-    .. math::
-
-        Q = f(1 - (k + 1)C)
-
-    IMPORTANT: internally, all computations are relying on extensions of the analyzer, that are computed
-    with default parameters if not present (i.e. correlograms, template_similarity, ...) If you want to
-    have a finer control on these values, please precompute the extensions before applying the auto_merge
-
-    Parameters
-    ----------
-    sorting_analyzer : SortingAnalyzer
-        The SortingAnalyzer
-    preset : "similarity_correlograms" | "x_contaminations" | "temporal_splits" | "feature_neighbors" | None, default: "similarity_correlograms"
-        The preset to use for the auto-merge. Presets combine different steps into a recipe and focus on:
-
-        * | "similarity_correlograms": mainly focused on template similarity and correlograms.
-          | It uses the following steps: "num_spikes", "remove_contaminated", "unit_locations",
-          | "template_similarity", "correlogram", "quality_score"
-        * | "x_contaminations": similar to "similarity_correlograms", but checks for cross-contamination instead of correlograms.
-          | It uses the following steps: "num_spikes", "remove_contaminated", "unit_locations",
-          | "template_similarity", "cross_contamination", "quality_score"
-        * | "temporal_splits": focused on finding temporal splits using presence distance.
-          | It uses the following steps: "num_spikes", "remove_contaminated", "unit_locations",
-          | "template_similarity", "presence_distance", "quality_score"
-        * | "feature_neighbors": focused on finding unit pairs whose spikes are close in the feature space using kNN.
-          | It uses the following steps: "num_spikes", "snr", "remove_contaminated", "unit_locations",
-          | "knn", "quality_score"
-
-        If `preset` is None, you can specify the steps manually with the `steps` parameter.
-    resolve_graph : bool, default: False
-        If True, the function resolves the potential unit pairs to be merged into multiple-unit merges.
-    min_spikes : int, default: 100
-        Minimum number of spikes for each unit to consider a potential merge.
-        Enough spikes are needed to estimate the correlogram
-    min_snr : float, default 2
-        Minimum Signal to Noise ratio for templates to be considered while merging
-    max_distance_um : float, default: 150
-        Maximum distance between units for considering a merge
-    corr_diff_thresh : float, default: 0.16
-        The threshold on the "correlogram distance metric" for considering a merge.
-        It needs to be between 0 and 1
-    template_diff_thresh : float, default: 0.25
-        The threshold on the "template distance metric" for considering a merge.
-        It needs to be between 0 and 1
-    contamination_thresh : float, default: 0.2
-        Threshold for not taking in account a unit when it is too contaminated.
-    presence_distance_thresh : float, default: 100
-        Parameter to control how present two units should be simultaneously.
-    p_value : float, default: 0.2
-        The p-value threshold for the cross-contamination test.
-    cc_thresh : float, default: 0.1
-        The threshold on the cross-contamination for considering a merge.
-    censored_period_ms : float, default: 0.3
-        Used to compute the refractory period violations aka "contamination".
-    refractory_period_ms : float, default: 1
-        Used to compute the refractory period violations aka "contamination".
-    sigma_smooth_ms : float, default: 0.6
-        Parameters to smooth the correlogram estimation.
-    adaptative_window_thresh : float, default: 0.5
-        Parameter to detect the window size in correlogram estimation.
-    censor_correlograms_ms : float, default: 0.15
-        The period to censor on the auto and cross-correlograms.
-    firing_contamination_balance : float, default: 1.5
-        Parameter to control the balance between firing rate and contamination in computing unit "quality score".
-    k_nn : int, default 5
-        The number of neighbors to consider for every spike in the recording.
-    knn_kwargs : dict, default None
-        The dict of extra params to be passed to knn.
-    extra_outputs : bool, default: False
-        If True, an additional dictionary (`outs`) with processed data is returned.
-    steps : None or list of str, default: None
-        Which steps to run, if no preset is used.
-        Pontential steps : "num_spikes", "snr", "remove_contaminated", "unit_locations", "correlogram",
-        "template_similarity", "presence_distance", "cross_contamination", "knn", "quality_score"
-        Please check steps explanations above!
-    presence_distance_kwargs : None|dict, default: None
-        A dictionary of kwargs to be passed to compute_presence_distance().
-
-    Returns
-    -------
-    potential_merges:
-        A list of tuples of 2 elements (if `resolve_graph`if false) or 2+ elements (if `resolve_graph` is true).
-        List of pairs that could be merged.
-    outs:
-        Returned only when extra_outputs=True
-        A dictionary that contains data for debugging and plotting.
-
-    References
-    ----------
-    This function is inspired and built upon similar functions from Lussac [Llobet]_,
-    done by Aurelien Wyngaard and Victor Llobet.
-    https://github.com/BarbourLab/lussac/blob/v1.0.0/postprocessing/merge_units.py
-    """
-    warnings.warn(
-        "get_potential_auto_merge() is deprecated. Use compute_merge_unit_groups() instead",
-        DeprecationWarning,
-        stacklevel=2,
-    )
-
-    presence_distance_kwargs = presence_distance_kwargs or dict()
-    knn_kwargs = knn_kwargs or dict()
-    return compute_merge_unit_groups(
-        sorting_analyzer,
-        preset,
-        resolve_graph,
-        steps_params={
-            "num_spikes": {"min_spikes": min_spikes},
-            "snr": {"min_snr": min_snr},
-            "remove_contaminated": {
-                "contamination_thresh": contamination_thresh,
-                "refractory_period_ms": refractory_period_ms,
-                "censored_period_ms": censored_period_ms,
-            },
-            "unit_locations": {"max_distance_um": max_distance_um},
-            "correlogram": {
-                "corr_diff_thresh": corr_diff_thresh,
-                "censor_correlograms_ms": censor_correlograms_ms,
-                "sigma_smooth_ms": sigma_smooth_ms,
-                "adaptative_window_thresh": adaptative_window_thresh,
-            },
-            "template_similarity": {"template_diff_thresh": template_diff_thresh},
-            "presence_distance": {"presence_distance_thresh": presence_distance_thresh, **presence_distance_kwargs},
-            "knn": {"k_nn": k_nn, **knn_kwargs},
-            "cross_contamination": {
-                "cc_thresh": cc_thresh,
-                "p_value": p_value,
-                "refractory_period_ms": refractory_period_ms,
-                "censored_period_ms": censored_period_ms,
-            },
-            "quality_score": {
-                "firing_contamination_balance": firing_contamination_balance,
-                "refractory_period_ms": refractory_period_ms,
-                "censored_period_ms": censored_period_ms,
-            },
-        },
-        compute_needed_extensions=True,
-        extra_outputs=extra_outputs,
-        steps=steps,
-    )
+        return potential_merges
 
 
 def get_pairs_via_nntree(sorting_analyzer, k_nn=5, pair_mask=None, **knn_kwargs):
@@ -857,10 +661,10 @@ def check_improve_contaminations_score(
         f_new = compute_firing_rates(sorting_analyzer_new)[unit_id1]
 
         # old and new scores
-        k = 1 + firing_contamination_balance
-        score_1 = f_1 * (1 - k * c_1)
-        score_2 = f_2 * (1 - k * c_2)
-        score_new = f_new * (1 - k * c_new)
+        k = firing_contamination_balance
+        score_1 = f_1 * (1 - (k + 1) * c_1)
+        score_2 = f_2 * (1 - (k + 1) * c_2)
+        score_new = f_new * (1 - (k + 1) * c_new)
 
         if score_new < score_1 or score_new < score_2:
             # the score is not improved
diff --git a/src/spikeinterface/exporters/report.py b/src/spikeinterface/exporters/report.py
index ab0840138..3a4be9213 100644
--- a/src/spikeinterface/exporters/report.py
+++ b/src/spikeinterface/exporters/report.py
@@ -20,10 +20,10 @@ def export_report(
     **job_kwargs,
 ):
     """
-    Exports a SI spike sorting report. The report includes summary figures of the spike sorting output.
-    What is plotted depends on what has been calculated. Unit locations and unit waveforms are always included.
-    Unit waveform densities, correlograms and spike amplitudes are plotted if `waveforms`, `correlograms`,
-    and `spike_amplitudes` have been computed for the given `sorting_analyzer`.
+    Exports a SI spike sorting report. The report includes summary figures of the spike sorting output
+    (e.g. amplitude distributions, unit localization and depth VS amplitude) as well as unit-specific reports,
+    that include waveforms, templates, template maps, ISI distributions, and more.
+
 
     Parameters
     ----------
diff --git a/src/spikeinterface/extractors/neoextractors/plexon2.py b/src/spikeinterface/extractors/neoextractors/plexon2.py
index e0604f749..2f360ed86 100644
--- a/src/spikeinterface/extractors/neoextractors/plexon2.py
+++ b/src/spikeinterface/extractors/neoextractors/plexon2.py
@@ -28,10 +28,6 @@ class Plexon2RecordingExtractor(NeoBaseRecordingExtractor):
             ids: ["source3.1" , "source3.2", "source3.3", "source3.4"]
     all_annotations : bool, default: False
         Load exhaustively all annotations from neo.
-    reading_attempts : int, default: 25
-        Number of attempts to read the file before raising an error
-        This opening process is somewhat unreliable and might fail occasionally. Adjust this higher
-        if you encounter problems in opening the file.
 
     Examples
     --------
@@ -41,16 +37,8 @@ class Plexon2RecordingExtractor(NeoBaseRecordingExtractor):
 
     NeoRawIOClass = "Plexon2RawIO"
 
-    def __init__(
-        self,
-        file_path,
-        stream_id=None,
-        stream_name=None,
-        use_names_as_ids=True,
-        all_annotations=False,
-        reading_attempts: int = 25,
-    ):
-        neo_kwargs = self.map_to_neo_kwargs(file_path, reading_attempts=reading_attempts)
+    def __init__(self, file_path, stream_id=None, stream_name=None, use_names_as_ids=True, all_annotations=False):
+        neo_kwargs = self.map_to_neo_kwargs(file_path)
         NeoBaseRecordingExtractor.__init__(
             self,
             stream_id=stream_id,
@@ -62,18 +50,8 @@ class Plexon2RecordingExtractor(NeoBaseRecordingExtractor):
         self._kwargs.update({"file_path": str(file_path)})
 
     @classmethod
-    def map_to_neo_kwargs(cls, file_path, reading_attempts: int = 25):
-
+    def map_to_neo_kwargs(cls, file_path):
         neo_kwargs = {"filename": str(file_path)}
-
-        from packaging.version import Version
-        import neo
-
-        neo_version = Version(neo.__version__)
-
-        if neo_version > Version("0.13.3"):
-            neo_kwargs["reading_attempts"] = reading_attempts
-
         return neo_kwargs
 
 
diff --git a/src/spikeinterface/postprocessing/localization_tools.py b/src/spikeinterface/postprocessing/localization_tools.py
index 59c12a992..e6278fc59 100644
--- a/src/spikeinterface/postprocessing/localization_tools.py
+++ b/src/spikeinterface/postprocessing/localization_tools.py
@@ -3,8 +3,6 @@ from __future__ import annotations
 import warnings
 
 import numpy as np
-from spikeinterface.core import SortingAnalyzer, Templates, compute_sparsity
-from spikeinterface.core.template_tools import _get_nbefore, get_dense_templates_array, get_template_extremum_channel
 
 try:
     import numba
@@ -14,6 +12,10 @@ except ImportError:
     HAVE_NUMBA = False
 
 
+from spikeinterface.core import compute_sparsity, SortingAnalyzer, Templates
+from spikeinterface.core.template_tools import get_template_extremum_channel, _get_nbefore, get_dense_templates_array
+
+
 def compute_monopolar_triangulation(
     sorting_analyzer_or_templates: SortingAnalyzer | Templates,
     unit_ids=None,
@@ -75,11 +77,7 @@ def compute_monopolar_triangulation(
 
     contact_locations = sorting_analyzer_or_templates.get_channel_locations()
 
-    if sorting_analyzer_or_templates.sparsity is None:
-        sparsity = compute_sparsity(sorting_analyzer_or_templates, method="radius", radius_um=radius_um)
-    else:
-        sparsity = sorting_analyzer_or_templates.sparsity
-
+    sparsity = compute_sparsity(sorting_analyzer_or_templates, method="radius", radius_um=radius_um)
     templates = get_dense_templates_array(
         sorting_analyzer_or_templates, return_scaled=get_return_scaled(sorting_analyzer_or_templates)
     )
@@ -108,7 +106,7 @@ def compute_monopolar_triangulation(
         # wf is (nsample, nchan) - chann is only nieghboor
         wf = templates[i, :, :][:, chan_inds]
         if feature == "ptp":
-            wf_data = np.ptp(wf, axis=0)
+            wf_data = wf.ptp(axis=0)
         elif feature == "energy":
             wf_data = np.linalg.norm(wf, axis=0)
         elif feature == "peak_voltage":
@@ -159,13 +157,9 @@ def compute_center_of_mass(
 
     assert feature in ["ptp", "mean", "energy", "peak_voltage"], f"{feature} is not a valid feature"
 
-    if sorting_analyzer_or_templates.sparsity is None:
-        sparsity = compute_sparsity(
-            sorting_analyzer_or_templates, peak_sign=peak_sign, method="radius", radius_um=radius_um
-        )
-    else:
-        sparsity = sorting_analyzer_or_templates.sparsity
-
+    sparsity = compute_sparsity(
+        sorting_analyzer_or_templates, peak_sign=peak_sign, method="radius", radius_um=radius_um
+    )
     templates = get_dense_templates_array(
         sorting_analyzer_or_templates, return_scaled=get_return_scaled(sorting_analyzer_or_templates)
     )
@@ -186,7 +180,7 @@ def compute_center_of_mass(
         wf = templates[i, :, :]
 
         if feature == "ptp":
-            wf_data = np.ptp(wf[:, chan_inds], axis=0)
+            wf_data = (wf[:, chan_inds]).ptp(axis=0)
         elif feature == "mean":
             wf_data = (wf[:, chan_inds]).mean(axis=0)
         elif feature == "energy":
@@ -656,55 +650,8 @@ if HAVE_NUMBA:
     enforce_decrease_shells = numba.jit(enforce_decrease_shells_data, nopython=True)
 
 
-def compute_location_max_channel(
-    templates_or_sorting_analyzer: SortingAnalyzer | Templates,
-    unit_ids=None,
-    peak_sign: "neg" | "pos" | "both" = "neg",
-    mode: "extremum" | "at_index" | "peak_to_peak" = "extremum",
-) -> np.ndarray:
-    """
-    Localize a unit using max channel.
-
-    This uses internally `get_template_extremum_channel()`
-
-
-    Parameters
-    ----------
-    templates_or_sorting_analyzer : SortingAnalyzer | Templates
-        A SortingAnalyzer or Templates object
-    unit_ids: list[str] | list[int] | None
-        A list of unit_id to restrict the computation
-    peak_sign :  "neg" | "pos" | "both"
-        Sign of the template to find extremum channels
-    mode : "extremum" | "at_index" | "peak_to_peak", default: "at_index"
-        Where the amplitude is computed
-        * "extremum" : take the peak value (max or min depending on `peak_sign`)
-        * "at_index" : take value at `nbefore` index
-        * "peak_to_peak" : take the peak-to-peak amplitude
-
-    Returns
-    -------
-    unit_locations: np.ndarray
-        2d
-    """
-    extremum_channels_index = get_template_extremum_channel(
-        templates_or_sorting_analyzer, peak_sign=peak_sign, mode=mode, outputs="index"
-    )
-    contact_locations = templates_or_sorting_analyzer.get_channel_locations()
-    if unit_ids is None:
-        unit_ids = templates_or_sorting_analyzer.unit_ids
-    else:
-        unit_ids = np.asarray(unit_ids)
-    unit_locations = np.zeros((unit_ids.size, 2), dtype="float32")
-    for i, unit_id in enumerate(unit_ids):
-        unit_locations[i, :] = contact_locations[extremum_channels_index[unit_id]]
-
-    return unit_locations
-
-
 _unit_location_methods = {
     "center_of_mass": compute_center_of_mass,
     "grid_convolution": compute_grid_convolution,
     "monopolar_triangulation": compute_monopolar_triangulation,
-    "max_channel": compute_location_max_channel,
 }
diff --git a/src/spikeinterface/postprocessing/template_metrics.py b/src/spikeinterface/postprocessing/template_metrics.py
index 6e7bcf21b..306e9594b 100644
--- a/src/spikeinterface/postprocessing/template_metrics.py
+++ b/src/spikeinterface/postprocessing/template_metrics.py
@@ -97,7 +97,7 @@ class ComputeTemplateMetrics(AnalyzerExtension):
 
     extension_name = "template_metrics"
     depend_on = ["templates"]
-    need_recording = False
+    need_recording = True
     use_nodepipeline = False
     need_job_kwargs = False
 
diff --git a/src/spikeinterface/postprocessing/template_similarity.py b/src/spikeinterface/postprocessing/template_similarity.py
index 6c30e2730..cfa9d89fe 100644
--- a/src/spikeinterface/postprocessing/template_similarity.py
+++ b/src/spikeinterface/postprocessing/template_similarity.py
@@ -44,7 +44,7 @@ class ComputeTemplateSimilarity(AnalyzerExtension):
 
     extension_name = "template_similarity"
     depend_on = ["templates"]
-    need_recording = False
+    need_recording = True
     use_nodepipeline = False
     need_job_kwargs = False
     need_backward_compatibility_on_load = True
diff --git a/src/spikeinterface/postprocessing/unit_locations.py b/src/spikeinterface/postprocessing/unit_locations.py
index 3f6dd47ee..4029fc88c 100644
--- a/src/spikeinterface/postprocessing/unit_locations.py
+++ b/src/spikeinterface/postprocessing/unit_locations.py
@@ -39,7 +39,7 @@ class ComputeUnitLocations(AnalyzerExtension):
 
     extension_name = "unit_locations"
     depend_on = ["templates"]
-    need_recording = False
+    need_recording = True
     use_nodepipeline = False
     need_job_kwargs = False
     need_backward_compatibility_on_load = True
diff --git a/src/spikeinterface/preprocessing/decimate.py b/src/spikeinterface/preprocessing/decimate.py
index d5fc9d202..334ebb02d 100644
--- a/src/spikeinterface/preprocessing/decimate.py
+++ b/src/spikeinterface/preprocessing/decimate.py
@@ -63,15 +63,18 @@ class DecimateRecording(BasePreprocessor):
                 f"Consider combining DecimateRecording with FrameSliceRecording for fine control on the recording start/end frames."
             )
         self._decimation_offset = decimation_offset
-        decimated_sampling_frequency = self._orig_samp_freq / self._decimation_factor
+        resample_rate = self._orig_samp_freq / self._decimation_factor
 
-        BasePreprocessor.__init__(self, recording, sampling_frequency=decimated_sampling_frequency)
+        BasePreprocessor.__init__(self, recording, sampling_frequency=resample_rate)
 
+        # in case there was a time_vector, it will be dropped for sanity.
+        # This is not necessary but consistent with ResampleRecording
         for parent_segment in recording._recording_segments:
+            parent_segment.time_vector = None
             self.add_recording_segment(
                 DecimateRecordingSegment(
                     parent_segment,
-                    decimated_sampling_frequency,
+                    resample_rate,
                     self._orig_samp_freq,
                     decimation_factor,
                     decimation_offset,
@@ -90,26 +93,22 @@ class DecimateRecordingSegment(BaseRecordingSegment):
     def __init__(
         self,
         parent_recording_segment,
-        decimated_sampling_frequency,
+        resample_rate,
         parent_rate,
         decimation_factor,
         decimation_offset,
         dtype,
     ):
-        if parent_recording_segment.time_vector is not None:
-            time_vector = parent_recording_segment.time_vector[decimation_offset::decimation_factor]
-            decimated_sampling_frequency = None
-            t_start = None
+        if parent_recording_segment.t_start is None:
+            new_t_start = None
         else:
-            time_vector = None
-            if parent_recording_segment.t_start is None:
-                t_start = None
-            else:
-                t_start = parent_recording_segment.t_start + (decimation_offset / parent_rate)
+            new_t_start = parent_recording_segment.t_start + decimation_offset / parent_rate
 
         # Do not use BasePreprocessorSegment bcause we have to reset the sampling rate!
         BaseRecordingSegment.__init__(
-            self, sampling_frequency=decimated_sampling_frequency, t_start=t_start, time_vector=time_vector
+            self,
+            sampling_frequency=resample_rate,
+            t_start=new_t_start,
         )
         self._parent_segment = parent_recording_segment
         self._decimation_factor = decimation_factor
diff --git a/src/spikeinterface/preprocessing/silence_periods.py b/src/spikeinterface/preprocessing/silence_periods.py
index 85169011d..8f38f0146 100644
--- a/src/spikeinterface/preprocessing/silence_periods.py
+++ b/src/spikeinterface/preprocessing/silence_periods.py
@@ -71,10 +71,8 @@ class SilencedPeriodsRecording(BasePreprocessor):
 
         if mode in ["noise"]:
             if noise_levels is None:
-                random_slices_kwargs = random_chunk_kwargs.copy()
-                random_slices_kwargs["seed"] = seed
                 noise_levels = get_noise_levels(
-                    recording, return_scaled=False, random_slices_kwargs=random_slices_kwargs
+                    recording, return_scaled=False, concatenated=True, seed=seed, **random_chunk_kwargs
                 )
             noise_generator = NoiseGeneratorRecording(
                 num_channels=recording.get_num_channels(),
diff --git a/src/spikeinterface/preprocessing/whiten.py b/src/spikeinterface/preprocessing/whiten.py
index 57400c119..195969ff7 100644
--- a/src/spikeinterface/preprocessing/whiten.py
+++ b/src/spikeinterface/preprocessing/whiten.py
@@ -19,8 +19,6 @@ class WhitenRecording(BasePreprocessor):
     recording : RecordingExtractor
         The recording extractor to be whitened.
     dtype : None or dtype, default: None
-        Datatype of the output recording (covariance matrix estimation
-        and whitening are performed in float32).
         If None the the parent dtype is kept.
         For integer dtype a int_scale must be also given.
     mode : "global" | "local", default: "global"
@@ -76,9 +74,7 @@ class WhitenRecording(BasePreprocessor):
         dtype_ = fix_dtype(recording, dtype)
 
         if dtype_.kind == "i":
-            assert (
-                int_scale is not None
-            ), "For recording with dtype=int you must set the output dtype to float OR set a int_scale"
+            assert int_scale is not None, "For recording with dtype=int you must set dtype=float32 OR set a int_scale"
 
         if W is not None:
             W = np.asarray(W)
@@ -128,7 +124,7 @@ class WhitenRecordingSegment(BasePreprocessorSegment):
     def get_traces(self, start_frame, end_frame, channel_indices):
         traces = self.parent_recording_segment.get_traces(start_frame, end_frame, slice(None))
         traces_dtype = traces.dtype
-        # if uint --> force float
+        # if uint --> force int
         if traces_dtype.kind == "u":
             traces = traces.astype("float32")
 
@@ -189,7 +185,6 @@ def compute_whitening_matrix(
 
     """
     random_data = get_random_data_chunks(recording, concatenated=True, return_scaled=False, **random_chunk_kwargs)
-    random_data = random_data.astype(np.float32)
 
     regularize_kwargs = regularize_kwargs if regularize_kwargs is not None else {"method": "GraphicalLassoCV"}
 
diff --git a/src/spikeinterface/sorters/basesorter.py b/src/spikeinterface/sorters/basesorter.py
index c59fa29c0..3502d2754 100644
--- a/src/spikeinterface/sorters/basesorter.py
+++ b/src/spikeinterface/sorters/basesorter.py
@@ -145,10 +145,9 @@ class BaseSorter:
         elif recording.check_serializability("pickle"):
             recording.dump(output_folder / "spikeinterface_recording.pickle", relative_to=output_folder)
         else:
-            raise RuntimeError(
-                "This recording is not serializable and so can not be sorted. Consider `recording.save()` to save a "
-                "compatible binary file."
-            )
+            # TODO: deprecate and finally remove this after 0.100
+            d = {"warning": "The recording is not serializable to json"}
+            rec_file.write_text(json.dumps(d, indent=4), encoding="utf8")
 
         return output_folder
 
diff --git a/src/spikeinterface/sorters/internal/spyking_circus2.py b/src/spikeinterface/sorters/internal/spyking_circus2.py
index eed693b34..211adba99 100644
--- a/src/spikeinterface/sorters/internal/spyking_circus2.py
+++ b/src/spikeinterface/sorters/internal/spyking_circus2.py
@@ -27,7 +27,7 @@ class Spykingcircus2Sorter(ComponentsBasedSorter):
         "general": {"ms_before": 2, "ms_after": 2, "radius_um": 75},
         "sparsity": {"method": "snr", "amplitude_mode": "peak_to_peak", "threshold": 0.25},
         "filtering": {"freq_min": 150, "freq_max": 7000, "ftype": "bessel", "filter_order": 2},
-        "whitening": {"mode": "local", "regularize": False},
+        "whitening": {"mode": "local", "regularize": True},
         "detection": {"peak_sign": "neg", "detect_threshold": 4},
         "selection": {
             "method": "uniform",
@@ -100,12 +100,6 @@ class Spykingcircus2Sorter(ComponentsBasedSorter):
         except:
             HAVE_HDBSCAN = False
 
-        try:
-            import torch
-        except ImportError:
-            HAVE_TORCH = False
-            print("spykingcircus2 could benefit from using torch. Consider installing it")
-
         assert HAVE_HDBSCAN, "spykingcircus2 needs hdbscan to be installed"
 
         # this is importanted only on demand because numba import are too heavy
diff --git a/src/spikeinterface/sortingcomponents/clustering/circus.py b/src/spikeinterface/sortingcomponents/clustering/circus.py
index 99c59f493..b7e71d3b4 100644
--- a/src/spikeinterface/sortingcomponents/clustering/circus.py
+++ b/src/spikeinterface/sortingcomponents/clustering/circus.py
@@ -92,7 +92,7 @@ class CircusClustering:
         # SVD for time compression
         few_peaks = select_peaks(peaks, recording=recording, method="uniform", n_peaks=10000, margin=(nbefore, nafter))
         few_wfs = extract_waveform_at_max_channel(
-            recording, few_peaks, ms_before=ms_before, ms_after=ms_after, **job_kwargs
+            recording, few_peaks, ms_before=ms_before, ms_after=ms_after, **params["job_kwargs"]
         )
 
         wfs = few_wfs[:, :, 0]
@@ -141,7 +141,7 @@ class CircusClustering:
             all_pc_data = run_node_pipeline(
                 recording,
                 pipeline_nodes,
-                job_kwargs,
+                params["job_kwargs"],
                 job_name="extracting features",
             )
 
@@ -176,7 +176,7 @@ class CircusClustering:
             _ = run_node_pipeline(
                 recording,
                 pipeline_nodes,
-                job_kwargs,
+                params["job_kwargs"],
                 job_name="extracting features",
                 gather_mode="npy",
                 gather_kwargs=dict(exist_ok=True),
diff --git a/src/spikeinterface/sortingcomponents/matching/circus.py b/src/spikeinterface/sortingcomponents/matching/circus.py
index 3b97f2dc6..d1b2139c5 100644
--- a/src/spikeinterface/sortingcomponents/matching/circus.py
+++ b/src/spikeinterface/sortingcomponents/matching/circus.py
@@ -18,15 +18,6 @@ spike_dtype = [
     ("segment_index", "int64"),
 ]
 
-try:
-    import torch
-    import torch.nn.functional as F
-
-    HAVE_TORCH = True
-    from torch.nn.functional import conv1d
-except ImportError:
-    HAVE_TORCH = False
-
 from .base import BaseTemplateMatching
 
 
@@ -52,9 +43,9 @@ def compress_templates(
 
     temporal, singular, spatial = np.linalg.svd(templates_array, full_matrices=False)
     # Keep only the strongest components
-    temporal = temporal[:, :, :approx_rank].astype(np.float32)
-    singular = singular[:, :approx_rank].astype(np.float32)
-    spatial = spatial[:, :approx_rank, :].astype(np.float32)
+    temporal = temporal[:, :, :approx_rank]
+    singular = singular[:, :approx_rank]
+    spatial = spatial[:, :approx_rank, :]
 
     if return_new_templates:
         templates_array = np.matmul(temporal * singular[:, np.newaxis, :], spatial)
@@ -116,22 +107,18 @@ class CircusOMPSVDPeeler(BaseTemplateMatching):
 
     Parameters
     ----------
-    amplitude : tuple
+    amplitude: tuple
         (Minimal, Maximal) amplitudes allowed for every template
-    max_failures : int
+    max_failures: int
         Stopping criteria of the OMP algorithm, as number of retry while updating amplitudes
-    sparse_kwargs : dict
+    sparse_kwargs: dict
         Parameters to extract a sparsity mask from the waveform_extractor, if not
         already sparse.
-    rank : int, default: 5
+    rank: int, default: 5
         Number of components used internally by the SVD
-    vicinity : int
+    vicinity: int
         Size of the area surrounding a spike to perform modification (expressed in terms
         of template temporal width)
-    engine : string in ["numpy", "torch", "auto"]. Default "auto"
-        The engine to use for the convolutions
-    torch_device : string in ["cpu", "cuda", None]. Default "cpu"
-        Controls torch device if the torch engine is selected
     -----
     """
 
@@ -161,8 +148,6 @@ class CircusOMPSVDPeeler(BaseTemplateMatching):
         ignore_inds=[],
         vicinity=2,
         precomputed=None,
-        engine="numpy",
-        torch_device="cpu",
     ):
 
         BaseTemplateMatching.__init__(self, recording, templates, return_output=True, parents=None)
@@ -173,19 +158,6 @@ class CircusOMPSVDPeeler(BaseTemplateMatching):
         self.nafter = templates.nafter
         self.sampling_frequency = recording.get_sampling_frequency()
         self.vicinity = vicinity * self.num_samples
-        assert engine in ["numpy", "torch", "auto"], "engine should be numpy, torch or auto"
-        if engine == "auto":
-            if HAVE_TORCH:
-                self.engine = "torch"
-            else:
-                self.engine = "numpy"
-        else:
-            if engine == "torch":
-                assert HAVE_TORCH, "please install torch to use the torch engine"
-            self.engine = engine
-
-        assert torch_device in ["cuda", "cpu", None]
-        self.torch_device = torch_device
 
         self.amplitudes = amplitudes
         self.stop_criteria = stop_criteria
@@ -211,7 +183,6 @@ class CircusOMPSVDPeeler(BaseTemplateMatching):
             self.unit_overlaps_tables[i][self.unit_overlaps_indices[i]] = np.arange(len(self.unit_overlaps_indices[i]))
 
         self.margin = 2 * self.num_samples
-        self.is_pushed = False
 
     def _prepare_templates(self):
 
@@ -283,14 +254,6 @@ class CircusOMPSVDPeeler(BaseTemplateMatching):
         self.temporal = np.moveaxis(self.temporal, [0, 1, 2], [1, 2, 0])
         self.singular = self.singular.T[:, :, np.newaxis]
 
-    def _push_to_torch(self):
-        if self.engine == "torch":
-            self.spatial = torch.as_tensor(self.spatial, device=self.torch_device)
-            self.singular = torch.as_tensor(self.singular, device=self.torch_device)
-            self.temporal = torch.as_tensor(self.temporal.copy(), device=self.torch_device).swapaxes(0, 1)
-            self.temporal = torch.flip(self.temporal, (2,))
-        self.is_pushed = True
-
     def get_extra_outputs(self):
         output = {}
         for key in self._more_output_keys:
@@ -305,15 +268,15 @@ class CircusOMPSVDPeeler(BaseTemplateMatching):
         import scipy
         from scipy import ndimage
 
-        if not self.is_pushed:
-            self._push_to_torch()
-
         (potrs,) = scipy.linalg.get_lapack_funcs(("potrs",), dtype=np.float32)
+
         (nrm2,) = scipy.linalg.get_blas_funcs(("nrm2",), dtype=np.float32)
 
-        omp_tol = np.finfo(np.float32).eps
-        neighbor_window = self.num_samples - 1
+        overlaps_array = self.overlaps
 
+        omp_tol = np.finfo(np.float32).eps
+        num_samples = self.nafter + self.nbefore
+        neighbor_window = num_samples - 1
         if isinstance(self.amplitudes, list):
             min_amplitude, max_amplitude = self.amplitudes
         else:
@@ -321,36 +284,27 @@ class CircusOMPSVDPeeler(BaseTemplateMatching):
             min_amplitude = min_amplitude[:, np.newaxis]
             max_amplitude = max_amplitude[:, np.newaxis]
 
-        if self.engine == "torch":
-            blank = np.zeros((neighbor_window, self.num_channels), dtype=np.float32)
-            traces = np.vstack((blank, traces, blank))
-            num_timesteps = traces.shape[0]
-            torch_traces = torch.as_tensor(traces.T[np.newaxis, :, :], device=self.torch_device)
-            num_templates, num_channels = self.temporal.shape[0], self.temporal.shape[1]
-            spatially_filtered_data = torch.matmul(self.spatial, torch_traces)
-            scaled_filtered_data = (spatially_filtered_data * self.singular).swapaxes(0, 1)
-            scaled_filtered_data_ = scaled_filtered_data.reshape(1, num_templates * num_channels, num_timesteps)
-            scalar_products = conv1d(scaled_filtered_data_, self.temporal, groups=num_templates, padding="valid")
-            scalar_products = scalar_products.cpu().numpy()[0, :, self.num_samples - 1 : -neighbor_window]
-        else:
-            num_timesteps = traces.shape[0]
-            num_peaks = num_timesteps - neighbor_window
-            conv_shape = (self.num_templates, num_peaks)
-            scalar_products = np.zeros(conv_shape, dtype=np.float32)
-            # Filter using overlap-and-add convolution
-            spatially_filtered_data = np.matmul(self.spatial, traces.T[np.newaxis, :, :])
-            scaled_filtered_data = spatially_filtered_data * self.singular
-            from scipy import signal
+        num_timesteps = len(traces)
 
-            objective_by_rank = signal.oaconvolve(scaled_filtered_data, self.temporal, axes=2, mode="valid")
-            scalar_products += np.sum(objective_by_rank, axis=0)
-
-        num_peaks = scalar_products.shape[1]
+        num_peaks = num_timesteps - num_samples + 1
+        conv_shape = (self.num_templates, num_peaks)
+        scalar_products = np.zeros(conv_shape, dtype=np.float32)
 
         # Filter using overlap-and-add convolution
         if len(self.ignore_inds) > 0:
-            scalar_products[self.ignore_inds] = -np.inf
             not_ignored = ~np.isin(np.arange(self.num_templates), self.ignore_inds)
+            spatially_filtered_data = np.matmul(self.spatial[:, not_ignored, :], traces.T[np.newaxis, :, :])
+            scaled_filtered_data = spatially_filtered_data * self.singular[:, not_ignored, :]
+            objective_by_rank = scipy.signal.oaconvolve(
+                scaled_filtered_data, self.temporal[:, not_ignored, :], axes=2, mode="valid"
+            )
+            scalar_products[not_ignored] += np.sum(objective_by_rank, axis=0)
+            scalar_products[self.ignore_inds] = -np.inf
+        else:
+            spatially_filtered_data = np.matmul(self.spatial, traces.T[np.newaxis, :, :])
+            scaled_filtered_data = spatially_filtered_data * self.singular
+            objective_by_rank = scipy.signal.oaconvolve(scaled_filtered_data, self.temporal, axes=2, mode="valid")
+            scalar_products += np.sum(objective_by_rank, axis=0)
 
         num_spikes = 0
 
@@ -368,7 +322,7 @@ class CircusOMPSVDPeeler(BaseTemplateMatching):
         is_in_vicinity = np.zeros(0, dtype=np.int32)
 
         if self.stop_criteria == "omp_min_sps":
-            stop_criteria = self.omp_min_sps * np.maximum(self.norms, np.sqrt(self.num_channels * self.num_samples))
+            stop_criteria = self.omp_min_sps * np.maximum(self.norms, np.sqrt(self.num_channels * num_samples))
         elif self.stop_criteria == "max_failures":
             num_valids = 0
             nb_failures = self.max_failures
@@ -400,11 +354,11 @@ class CircusOMPSVDPeeler(BaseTemplateMatching):
 
                 if num_selection > 0:
                     delta_t = selection[1] - peak_index
-                    idx = np.flatnonzero((delta_t < self.num_samples) & (delta_t > -self.num_samples))
+                    idx = np.flatnonzero((delta_t < num_samples) & (delta_t > -num_samples))
                     myline = neighbor_window + delta_t[idx]
                     myindices = selection[0, idx]
 
-                    local_overlaps = self.overlaps[best_cluster_ind]
+                    local_overlaps = overlaps_array[best_cluster_ind]
                     overlapping_templates = self.unit_overlaps_indices[best_cluster_ind]
                     table = self.unit_overlaps_tables[best_cluster_ind]
 
@@ -482,10 +436,10 @@ class CircusOMPSVDPeeler(BaseTemplateMatching):
                 for i in modified:
                     tmp_best, tmp_peak = sub_selection[:, i]
                     diff_amp = diff_amplitudes[i] * self.norms[tmp_best]
-                    local_overlaps = self.overlaps[tmp_best]
+                    local_overlaps = overlaps_array[tmp_best]
                     overlapping_templates = self.units_overlaps[tmp_best]
                     tmp = tmp_peak - neighbor_window
-                    idx = [max(0, tmp), min(num_peaks, tmp_peak + self.num_samples)]
+                    idx = [max(0, tmp), min(num_peaks, tmp_peak + num_samples)]
                     tdx = [idx[0] - tmp, idx[1] - tmp]
                     to_add = diff_amp * local_overlaps[:, tdx[0] : tdx[1]]
                     scalar_products[overlapping_templates, idx[0] : idx[1]] -= to_add
@@ -546,27 +500,27 @@ class CircusPeeler(BaseTemplateMatching):
 
     Parameters
     ----------
-    peak_sign : str
+    peak_sign: str
         Sign of the peak (neg, pos, or both)
-    exclude_sweep_ms : float
+    exclude_sweep_ms: float
         The number of samples before/after to classify a peak (should be low)
-    jitter : int
+    jitter: int
         The number of samples considered before/after every peak to search for
         matches
-    detect_threshold : int
+    detect_threshold: int
         The detection threshold
-    noise_levels : array
+    noise_levels: array
         The noise levels, for every channels
-    random_chunk_kwargs : dict
+    random_chunk_kwargs: dict
         Parameters for computing noise levels, if not provided (sub optimal)
-    max_amplitude : float
+    max_amplitude: float
         Maximal amplitude allowed for every template
-    min_amplitude : float
+    min_amplitude: float
         Minimal amplitude allowed for every template
-    use_sparse_matrix_threshold : float
+    use_sparse_matrix_threshold: float
         If density of the templates is below a given threshold, sparse matrix
         are used (memory efficient)
-    sparse_kwargs : dict
+    sparse_kwargs: dict
         Parameters to extract a sparsity mask from the waveform_extractor, if not
         already sparse.
     -----
diff --git a/src/spikeinterface/sortingcomponents/matching/tdc.py b/src/spikeinterface/sortingcomponents/matching/tdc.py
index 125baa3bd..56457fe2f 100644
--- a/src/spikeinterface/sortingcomponents/matching/tdc.py
+++ b/src/spikeinterface/sortingcomponents/matching/tdc.py
@@ -2,11 +2,15 @@ from __future__ import annotations
 
 import numpy as np
 from spikeinterface.core import (
+    get_noise_levels,
     get_channel_distances,
+    compute_sparsity,
     get_template_extremum_channel,
 )
 
-from spikeinterface.sortingcomponents.peak_detection import DetectPeakLocallyExclusive, DetectPeakMatchedFiltering
+from spikeinterface.sortingcomponents.peak_detection import DetectPeakLocallyExclusive
+from spikeinterface.core.template import Templates
+
 from .base import BaseTemplateMatching, _base_matching_dtype
 
 
@@ -21,7 +25,7 @@ except ImportError:
 
 class TridesclousPeeler(BaseTemplateMatching):
     """
-    Template-matching used by Tridesclous sorter.
+    Template-matching ported from Tridesclous sorter.
 
     The idea of this peeler is pretty simple.
     1. Find peaks
@@ -30,10 +34,8 @@ class TridesclousPeeler(BaseTemplateMatching):
     4. remove it from traces.
     5. in the residual find peaks again
 
-    Contrary tp circus_peeler or wobble, this template matching is working directly one the waveforms.
-    There is no SVD decomposition
-
-
+    This method is quite fast but don't give exelent results to resolve
+    spike collision when templates have high similarity.
     """
 
     def __init__(
@@ -43,29 +45,26 @@ class TridesclousPeeler(BaseTemplateMatching):
         parents=None,
         templates=None,
         peak_sign="neg",
-        exclude_sweep_ms=0.5,
         peak_shift_ms=0.2,
         detect_threshold=5,
         noise_levels=None,
-        use_fine_detector=True,
-        # TODO optimize theses radius
-        detection_radius_um=80.0,
-        cluster_radius_um=150.0,
-        amplitude_fitting_radius_um=150.0,
-        sample_shift=2,
-        ms_before=0.5,
-        ms_after=0.8,
-        max_peeler_loop=2,
-        amplitude_limits=(0.7, 1.4),
+        radius_um=100.0,
+        num_closest=5,
+        sample_shift=3,
+        ms_before=0.8,
+        ms_after=1.2,
+        num_peeler_loop=2,
+        num_template_try=1,
     ):
 
         BaseTemplateMatching.__init__(self, recording, templates, return_output=True, parents=None)
 
+        # maybe in base?
+        self.templates_array = templates.get_dense_templates()
+
         unit_ids = templates.unit_ids
         channel_ids = recording.channel_ids
 
-        num_templates = unit_ids.size
-
         sr = recording.sampling_frequency
 
         self.nbefore = templates.nbefore
@@ -83,9 +82,8 @@ class TridesclousPeeler(BaseTemplateMatching):
         s1 = -(templates.nafter - nafter_short)
         if s1 == 0:
             s1 = None
-
         # TODO check with out copy
-        self.sparse_templates_array_short = templates.templates_array[:, slice(s0, s1), :].copy()
+        self.templates_short = self.templates_array[:, slice(s0, s1), :].copy()
 
         self.peak_shift = int(peak_shift_ms / 1000 * sr)
 
@@ -94,12 +92,12 @@ class TridesclousPeeler(BaseTemplateMatching):
         self.abs_thresholds = noise_levels * detect_threshold
 
         channel_distance = get_channel_distances(recording)
-        self.neighbours_mask = channel_distance <= detection_radius_um
+        self.neighbours_mask = channel_distance < radius_um
 
         if templates.sparsity is not None:
-            self.sparsity_mask = templates.sparsity.mask
+            self.template_sparsity = templates.sparsity.mask
         else:
-            self.sparsity_mask = np.ones((unit_ids.size, channel_ids.size), dtype=bool)
+            self.template_sparsity = np.ones((unit_ids.size, channel_ids.size), dtype=bool)
 
         extremum_chan = get_template_extremum_channel(templates, peak_sign=peak_sign, outputs="index")
         # as numpy vector
@@ -111,108 +109,72 @@ class TridesclousPeeler(BaseTemplateMatching):
         # distance between units
         import scipy
 
-        # nearby cluster for each channel
+        unit_distances = scipy.spatial.distance.cdist(unit_locations, unit_locations, metric="euclidean")
+
+        # seach for closet units and unitary discriminant vector
+        closest_units = []
+        for unit_ind, unit_id in enumerate(unit_ids):
+            order = np.argsort(unit_distances[unit_ind, :])
+            closest_u = np.arange(unit_ids.size)[order].tolist()
+            closest_u.remove(unit_ind)
+            closest_u = np.array(closest_u[:num_closest])
+
+            # compute unitary discriminent vector
+            (chans,) = np.nonzero(self.template_sparsity[unit_ind, :])
+            template_sparse = self.templates_array[unit_ind, :, :][:, chans]
+            closest_vec = []
+            # against N closets
+            for u in closest_u:
+                vec = self.templates_array[u, :, :][:, chans] - template_sparse
+                vec /= np.sum(vec**2)
+                closest_vec.append((u, vec))
+            # against noise
+            closest_vec.append((None, -template_sparse / np.sum(template_sparse**2)))
+
+            closest_units.append(closest_vec)
+
+        self.closest_units = closest_units
+
+        # distance channel from unit
+        import scipy
+
         distances = scipy.spatial.distance.cdist(channel_locations, unit_locations, metric="euclidean")
-        near_cluster_mask = distances <= cluster_radius_um
+        near_cluster_mask = distances < radius_um
+
+        # nearby cluster for each channel
         self.possible_clusters_by_channel = []
         for channel_index in range(distances.shape[0]):
             (cluster_inds,) = np.nonzero(near_cluster_mask[channel_index, :])
             self.possible_clusters_by_channel.append(cluster_inds)
 
-        # precompute template norms ons sparse channels
-        self.template_norms = np.zeros(num_templates, dtype="float32")
-        for i in range(unit_ids.size):
-            chan_mask = self.sparsity_mask[i, :]
-            n = np.sum(chan_mask)
-            template = templates.templates_array[i, :, :n]
-            self.template_norms[i] = np.sum(template**2)
-
-        #
-        distances = scipy.spatial.distance.cdist(channel_locations, channel_locations, metric="euclidean")
-        self.near_chan_mask = distances <= amplitude_fitting_radius_um
-
         self.possible_shifts = np.arange(-sample_shift, sample_shift + 1, dtype="int64")
 
-        self.max_peeler_loop = max_peeler_loop
-        self.amplitude_limits = amplitude_limits
-
-        self.fast_spike_detector = DetectPeakLocallyExclusive(
-            recording=recording,
-            peak_sign=peak_sign,
-            detect_threshold=detect_threshold,
-            exclude_sweep_ms=exclude_sweep_ms,
-            radius_um=detection_radius_um,
-            noise_levels=noise_levels,
-        )
+        self.num_peeler_loop = num_peeler_loop
+        self.num_template_try = num_template_try
 
-        ##get prototype from best channel of each template
-        prototype = np.zeros(self.nbefore + self.nafter, dtype="float32")
-        for i in range(num_templates):
-            template = templates.templates_array[i, :, :]
-            chan_ind = np.argmax(np.abs(template[self.nbefore, :]))
-            if template[self.nbefore, chan_ind] != 0:
-                prototype += template[:, chan_ind] / np.abs(template[self.nbefore, chan_ind])
-        prototype /= np.abs(prototype[self.nbefore])
-
-        # import matplotlib.pyplot as plt
-        # fig,ax = plt.subplots()
-        # ax.plot(prototype)
-        # plt.show()
-
-        self.use_fine_detector = use_fine_detector
-        if self.use_fine_detector:
-            self.fine_spike_detector = DetectPeakMatchedFiltering(
-                recording=recording,
-                prototype=prototype,
-                ms_before=templates.nbefore / sr * 1000.0,
-                peak_sign="neg",
-                detect_threshold=detect_threshold,
-                exclude_sweep_ms=exclude_sweep_ms,
-                radius_um=detection_radius_um,
-                weight_method=dict(
-                    z_list_um=np.array([50.0]),
-                    sigma_3d=2.5,
-                    mode="exponential_3d",
-                ),
-                noise_levels=None,
-            )
-
-        self.detector_margin0 = self.fast_spike_detector.get_trace_margin()
-        self.detector_margin1 = self.fine_spike_detector.get_trace_margin() if use_fine_detector else 0
-        self.peeler_margin = max(self.nbefore, self.nafter) * 2
-        self.margin = max(self.peeler_margin, self.detector_margin0, self.detector_margin1)
+        self.margin = max(self.nbefore, self.nafter) * 2
 
     def get_trace_margin(self):
         return self.margin
 
     def compute_matching(self, traces, start_frame, end_frame, segment_index):
-
-        # TODO check if this is usefull
-        residuals = traces.copy()
+        traces = traces.copy()
 
         all_spikes = []
         level = 0
-        spikes_prev_loop = np.zeros(0, dtype=_base_matching_dtype)
-        use_fine_detector_level = False
         while True:
-            # print('level', level)
-            spikes = self._find_spikes_one_level(residuals, spikes_prev_loop, use_fine_detector_level, level)
-            if spikes.size > 0:
-                all_spikes.append(spikes)
+            # spikes = _tdc_find_spikes(traces, d, level=level)
+            spikes = self._find_spikes_one_level(traces, level=level)
+            keep = spikes["cluster_index"] >= 0
 
-            level += 1
+            if not np.any(keep):
+                break
+            all_spikes.append(spikes[keep])
 
-            # TODO concatenate all spikes for this instead of prev loop
-            spikes_prev_loop = spikes
+            level += 1
 
-            if (spikes.size == 0) or (level == self.max_peeler_loop):
-                if self.use_fine_detector and not use_fine_detector_level:
-                    # extra loop with fine detector
-                    use_fine_detector_level = True
-                    level = self.max_peeler_loop - 1
-                    continue
-                else:
-                    break
+            if level == self.num_peeler_loop:
+                break
 
         if len(all_spikes) > 0:
             all_spikes = np.concatenate(all_spikes)
@@ -223,34 +185,13 @@ class TridesclousPeeler(BaseTemplateMatching):
 
         return all_spikes
 
-    def _find_spikes_one_level(self, traces, spikes_prev_loop, use_fine_detector, level):
-
-        # print(use_fine_detector, level)
-
-        # TODO change the threhold dynaically depending the level
-        # peak_traces = traces[self.detector_margin : -self.detector_margin, :]
-
-        # peak_sample_ind, peak_chan_ind = DetectPeakLocallyExclusive.detect_peaks(
-        #     peak_traces, self.peak_sign, self.abs_thresholds, self.peak_shift, self.neighbours_mask
-        # )
-
-        if use_fine_detector:
-            peak_detector = self.fine_spike_detector
-        else:
-            peak_detector = self.fast_spike_detector
+    def _find_spikes_one_level(self, traces, level=0):
 
-        detector_margin = peak_detector.get_trace_margin()
-        if self.peeler_margin > detector_margin:
-            margin_shift = self.peeler_margin - detector_margin
-            sl = slice(margin_shift, -margin_shift)
-        else:
-            sl = slice(None)
-            margin_shift = 0
-        peak_traces = traces[sl, :]
-        (peaks,) = peak_detector.compute(peak_traces, None, None, 0, self.margin)
-        peak_sample_ind = peaks["sample_index"]
-        peak_chan_ind = peaks["channel_index"]
-        peak_sample_ind += margin_shift
+        peak_traces = traces[self.margin // 2 : -self.margin // 2, :]
+        peak_sample_ind, peak_chan_ind = DetectPeakLocallyExclusive.detect_peaks(
+            peak_traces, self.peak_sign, self.abs_thresholds, self.peak_shift, self.neighbours_mask
+        )
+        peak_sample_ind += self.margin // 2
 
         peak_amplitude = traces[peak_sample_ind, peak_chan_ind]
         order = np.argsort(np.abs(peak_amplitude))[::-1]
@@ -259,438 +200,153 @@ class TridesclousPeeler(BaseTemplateMatching):
 
         spikes = np.zeros(peak_sample_ind.size, dtype=_base_matching_dtype)
         spikes["sample_index"] = peak_sample_ind
-        spikes["channel_index"] = peak_chan_ind
+        spikes["channel_index"] = peak_chan_ind  # TODO need to put the channel from template
 
-        distances_shift = np.zeros(self.possible_shifts.size)
+        possible_shifts = self.possible_shifts
+        distances_shift = np.zeros(possible_shifts.size)
 
-        delta_sample = max(self.nbefore, self.nafter)  #  TODO check this maybe add margin
-        # neighbors_spikes_inds = get_neighbors_spikes(spikes["sample_index"], spikes["channel_index"], delta_sample, self.near_chan_mask)
-
-        # neighbors in actual and previous level
-        neighbors_spikes_inds = get_neighbors_spikes(
-            np.concatenate([spikes["sample_index"], spikes_prev_loop["sample_index"]]),
-            np.concatenate([spikes["channel_index"], spikes_prev_loop["channel_index"]]),
-            delta_sample,
-            self.near_chan_mask,
-        )
-
-        for i in range(spikes.size):
+        for i in range(peak_sample_ind.size):
             sample_index = peak_sample_ind[i]
 
             chan_ind = peak_chan_ind[i]
             possible_clusters = self.possible_clusters_by_channel[chan_ind]
 
             if possible_clusters.size > 0:
-                cluster_index = get_most_probable_cluster(
-                    traces,
-                    self.sparse_templates_array_short,
-                    possible_clusters,
-                    sample_index,
-                    chan_ind,
-                    self.nbefore_short,
-                    self.nafter_short,
-                    self.sparsity_mask,
-                )
-
-                # import matplotlib.pyplot as plt
-                # fig, ax = plt.subplots()
-                # chans = np.any(self.sparsity_mask[possible_clusters, :], axis=0)
-                # wf = traces[sample_index - self.nbefore : sample_index + self.nafter][:, chans]
-                # ax.plot(wf.T.flatten(), color='k')
-                # dense_templates_array = self.templates.get_dense_templates()
-                # for c_ind in possible_clusters:
-                #     template = dense_templates_array[c_ind, :, :][:, chans]
-                #     ax.plot(template.T.flatten())
-                #     if c_ind == cluster_index:
-                #         ax.plot(template.T.flatten(), color='m', ls='--')
-                #     ax.set_title(f"use_fine_detector{use_fine_detector} level{level}")
-                # plt.show()
-
-                chan_sparsity_mask = self.sparsity_mask[cluster_index, :]
-
-                # find best shift
-                numba_best_shift_sparse(
-                    traces,
-                    self.sparse_templates_array_short[cluster_index, :, :],
-                    sample_index,
-                    self.nbefore_short,
-                    self.possible_shifts,
-                    distances_shift,
-                    chan_sparsity_mask,
-                )
-
-                ind_shift = np.argmin(distances_shift)
-                shift = self.possible_shifts[ind_shift]
-
-                # TODO DEBUG shift later
-                spikes["sample_index"][i] += shift
-
-                spikes["cluster_index"][i] = cluster_index
-
-                # check that the the same cluster is not already detected at same place
-                # this can happen for small template the substract forvever the traces
-                outer_neighbors_inds = [ind for ind in neighbors_spikes_inds[i] if ind > i and ind >= spikes.size]
-                is_valid = True
-                for b in outer_neighbors_inds:
-                    b = b - spikes.size
-                    if (spikes[i]["sample_index"] == spikes_prev_loop[b]["sample_index"]) and (
-                        spikes[i]["cluster_index"] == spikes_prev_loop[b]["cluster_index"]
-                    ):
-                        is_valid = False
-
-                if is_valid:
-                    # temporary assign a cluster to neighbors if not done yet
-                    inner_neighbors_inds = [ind for ind in neighbors_spikes_inds[i] if (ind > i and ind < spikes.size)]
-                    for b in inner_neighbors_inds:
-                        spikes["cluster_index"][b] = get_most_probable_cluster(
-                            traces,
-                            self.sparse_templates_array_short,
-                            possible_clusters,
-                            spikes["sample_index"][b],
-                            spikes["channel_index"][b],
-                            self.nbefore_short,
-                            self.nafter_short,
-                            self.sparsity_mask,
-                        )
-
-                    amp = fit_one_amplitude_with_neighbors(
-                        spikes[i],
-                        spikes[inner_neighbors_inds],
-                        traces,
-                        self.sparsity_mask,
-                        self.templates.templates_array,
-                        self.template_norms,
-                        self.nbefore,
-                        self.nafter,
-                    )
-
-                    low_lim, up_lim = self.amplitude_limits
-                    if low_lim <= amp <= up_lim:
-                        spikes["amplitude"][i] = amp
-                        wanted_channel_mask = np.ones(traces.shape[1], dtype=bool)  # TODO move this before the loop
-                        construct_prediction_sparse(
-                            spikes[i : i + 1],
-                            traces,
-                            self.templates.templates_array,
-                            self.sparsity_mask,
-                            wanted_channel_mask,
-                            self.nbefore,
-                            additive=False,
-                        )
-                    elif low_lim > amp:
-                        # print("bad amp", amp)
-                        spikes["cluster_index"][i] = -1
-
-                        # import matplotlib.pyplot as plt
-                        # fig, ax = plt.subplots()
-                        # sample_ind = spikes["sample_index"][i]
-                        # print(chan_sparsity_mask)
-                        # wf = traces[sample_ind - self.nbefore : sample_ind + self.nafter][:, chan_sparsity_mask]
-                        # dense_templates_array = self.templates.get_dense_templates()
-                        # template = dense_templates_array[cluster_index, :, :][:, chan_sparsity_mask]
-                        # ax.plot(wf.T.flatten())
-                        # ax.plot(template.T.flatten())
-                        # ax.plot(template.T.flatten() * amp)
-                        # ax.set_title(f"amp{amp} use_fine_detector{use_fine_detector} level{level}")
-                        # plt.show()
-                    else:
-                        # amp > up_lim
-                        # TODO should try other cluster for the fit!!
-                        # spikes["cluster_index"][i] = -1
-
-                        # force amplitude to be one and need a fiting at next level
-                        spikes["amplitude"][i] = 1
-
-                        # print(amp)
-                        # import matplotlib.pyplot as plt
-                        # fig, ax = plt.subplots()
-                        # sample_ind = spikes["sample_index"][i]
-                        # wf = traces[sample_ind - self.nbefore : sample_ind + self.nafter][:, chan_sparsity_mask]
-                        # dense_templates_array = self.templates.get_dense_templates()
-                        # template = dense_templates_array[cluster_index, :, :][:, chan_sparsity_mask]
-                        # ax.plot(wf.T.flatten())
-                        # ax.plot(template.T.flatten())
-                        # ax.plot(template.T.flatten() * amp)
-                        # ax.set_title(f"amp{amp} use_fine_detector{use_fine_detector} level{level}")
-                        # plt.show()
-
-                        # import matplotlib.pyplot as plt
-                        # fig, ax = plt.subplots()
-                        # chans = np.any(self.sparsity_mask[possible_clusters, :], axis=0)
-                        # wf = traces[sample_index - self.nbefore : sample_index + self.nafter][:, chans]
-                        # ax.plot(wf.T.flatten(), color='k')
-                        # dense_templates_array = self.templates.get_dense_templates()
-                        # for c_ind in possible_clusters:
-                        #     template = dense_templates_array[c_ind, :, :][:, chans]
-                        #     ax.plot(template.T.flatten())
-                        #     if c_ind == cluster_index:
-                        #         ax.plot(template.T.flatten(), color='m', ls='--')
-                        #     ax.set_title(f"use_fine_detector{use_fine_detector} level{level}")
-                        # plt.show()
-
-                else:
-                    # not valid because already detected
-                    spikes["cluster_index"][i] = -1
-
-            else:
-                # no possible cluster in neighborhood for this channel
-                spikes["cluster_index"][i] = -1
-
-        # delta_sample = self.nbefore + self.nafter
-        # # TODO benchmark this and make this faster
-        # neighbors_spikes_inds = get_neighbors_spikes(spikes["sample_index"], spikes["channel_index"], delta_sample, self.near_chan_mask)
-        # for i in range(spikes.size):
-        #     amp = fit_one_amplitude_with_neighbors(spikes[i], spikes[neighbors_spikes_inds[i]],  traces,
-        #                                      self.sparsity_mask, self.templates.templates_array, self.nbefore, self.nafter)
-        #     spikes["amplitude"][i] = amp
-
-        keep = spikes["cluster_index"] >= 0
-        spikes = spikes[keep]
+                # ~ s0 = sample_index - d['nbefore']
+                # ~ s1 = sample_index + d['nafter']
 
-        # keep = (spikes["amplitude"] >= 0.7) & (spikes["amplitude"] <= 1.4)
-        # spikes = spikes[keep]
+                # ~ wf = traces[s0:s1, :]
 
-        # sparse_templates_array = self.templates.templates_array
-        # wanted_channel_mask = np.ones(traces.shape[1], dtype=bool)
-        # assert np.sum(wanted_channel_mask) == traces.shape[1] # TODO remove this DEBUG later
-        # construct_prediction_sparse(spikes, traces, sparse_templates_array, self.sparsity_mask, wanted_channel_mask, self.nbefore, additive=False)
+                s0 = sample_index - self.nbefore_short
+                s1 = sample_index + self.nafter_short
+                wf_short = traces[s0:s1, :]
 
-        return spikes
-
-
-def get_most_probable_cluster(
-    traces,
-    sparse_templates_array,
-    possible_clusters,
-    sample_index,
-    chan_ind,
-    nbefore_short,
-    nafter_short,
-    template_sparsity_mask,
-):
-    s0 = sample_index - nbefore_short
-    s1 = sample_index + nafter_short
-    wf_short = traces[s0:s1, :]
-
-    ## numba with cluster+channel spasity
-    union_channels = np.any(template_sparsity_mask[possible_clusters, :], axis=0)
-    distances = numba_sparse_distance(
-        wf_short, sparse_templates_array, template_sparsity_mask, union_channels, possible_clusters
-    )
+                ## pure numpy with cluster spasity
+                # distances = np.sum(np.sum((templates[possible_clusters, :, :] - wf[None, : , :])**2, axis=1), axis=1)
 
-    ind = np.argmin(distances)
-    cluster_index = possible_clusters[ind]
+                ## pure numpy with cluster+channel spasity
+                # union_channels, = np.nonzero(np.any(d['template_sparsity'][possible_clusters, :], axis=0))
+                # distances = np.sum(np.sum((templates[possible_clusters][:, :, union_channels] - wf[: , union_channels][None, : :])**2, axis=1), axis=1)
 
-    return cluster_index
+                ## numba with cluster+channel spasity
+                union_channels = np.any(self.template_sparsity[possible_clusters, :], axis=0)
+                # distances = numba_sparse_dist(wf, templates, union_channels, possible_clusters)
+                distances = numba_sparse_dist(wf_short, self.templates_short, union_channels, possible_clusters)
 
+                # DEBUG
+                # ~ ind = np.argmin(distances)
+                # ~ cluster_index = possible_clusters[ind]
 
-def get_neighbors_spikes(sample_inds, chan_inds, delta_sample, near_chan_mask):
+                for ind in np.argsort(distances)[: self.num_template_try]:
+                    cluster_index = possible_clusters[ind]
 
-    neighbors_spikes_inds = []
-    for i in range(sample_inds.size):
+                    chan_sparsity = self.template_sparsity[cluster_index, :]
+                    template_sparse = self.templates_array[cluster_index, :, :][:, chan_sparsity]
 
-        inds = np.flatnonzero(np.abs(sample_inds - sample_inds[i]) < delta_sample)
-        neighb = []
-        for ind in inds:
-            if near_chan_mask[chan_inds[i], chan_inds[ind]] and i != ind:
-                neighb.append(ind)
-        neighbors_spikes_inds.append(neighb)
+                    # find best shift
 
-    return neighbors_spikes_inds
+                    ## pure numpy version
+                    # for s, shift in enumerate(possible_shifts):
+                    #     wf_shift = traces[s0 + shift: s1 + shift, chan_sparsity]
+                    #     distances_shift[s] = np.sum((template_sparse - wf_shift)**2)
+                    # ind_shift = np.argmin(distances_shift)
+                    # shift = possible_shifts[ind_shift]
 
+                    ## numba version
+                    numba_best_shift(
+                        traces,
+                        self.templates_array[cluster_index, :, :],
+                        sample_index,
+                        self.nbefore,
+                        possible_shifts,
+                        distances_shift,
+                        chan_sparsity,
+                    )
+                    ind_shift = np.argmin(distances_shift)
+                    shift = possible_shifts[ind_shift]
+
+                    sample_index = sample_index + shift
+                    s0 = sample_index - self.nbefore
+                    s1 = sample_index + self.nafter
+                    wf_sparse = traces[s0:s1, chan_sparsity]
+
+                    # accept or not
+
+                    centered = wf_sparse - template_sparse
+                    accepted = True
+                    for other_ind, other_vector in self.closest_units[cluster_index]:
+                        v = np.sum(centered * other_vector)
+                        if np.abs(v) > 0.5:
+                            accepted = False
+                            break
+
+                    if accepted:
+                        # ~ if ind != np.argsort(distances)[0]:
+                        # ~ print('not first one', np.argsort(distances), ind)
+                        break
+
+                if accepted:
+                    amplitude = 1.0
+
+                    # remove template
+                    template = self.templates_array[cluster_index, :, :]
+                    s0 = sample_index - self.nbefore
+                    s1 = sample_index + self.nafter
+                    traces[s0:s1, :] -= template * amplitude
 
-def fit_one_amplitude_with_neighbors(
-    spike, neighbors_spikes, traces, template_sparsity_mask, sparse_templates_array, template_norms, nbefore, nafter
-):
-    """
-    Fit amplitude one spike of one spike with/without neighbors
+                else:
+                    cluster_index = -1
+                    amplitude = 0.0
 
-    """
+            else:
+                cluster_index = -1
+                amplitude = 0.0
 
-    import scipy.linalg
-
-    cluster_index = spike["cluster_index"]
-    sample_index = spike["sample_index"]
-    chan_sparsity_mask = template_sparsity_mask[cluster_index, :]
-    num_chans = np.sum(chan_sparsity_mask)
-    if num_chans == 0:
-        # protect against empty template because too sparse
-        return 0.0
-    start, stop = sample_index - nbefore, sample_index + nafter
-    if neighbors_spikes is None or (neighbors_spikes.size == 0):
-        template = sparse_templates_array[cluster_index, :, :num_chans]
-        wf = traces[start:stop, :][:, chan_sparsity_mask]
-        # TODO precompute template norms
-        amplitude = np.sum(template.flatten() * wf.flatten()) / template_norms[cluster_index]
-    else:
-
-        lim0 = min(start, np.min(neighbors_spikes["sample_index"]) - nbefore)
-        lim1 = max(stop, np.max(neighbors_spikes["sample_index"]) + nafter)
-
-        local_traces = traces[lim0:lim1, :][:, chan_sparsity_mask]
-        mask_not_fitted = (neighbors_spikes["amplitude"] == 0.0) & (neighbors_spikes["cluster_index"] >= 0)
-        local_spike = spike.copy()
-        local_spike["sample_index"] -= lim0
-        local_spike["amplitude"] = 1.0
-
-        local_neighbors_spikes = neighbors_spikes.copy()
-        local_neighbors_spikes["sample_index"] -= lim0
-        local_neighbors_spikes["amplitude"][:] = 1.0
-
-        num_spikes_to_fit = 1 + np.sum(mask_not_fitted)
-        x = np.zeros((lim1 - lim0, num_chans, num_spikes_to_fit), dtype="float32")
-        wanted_channel_mask = chan_sparsity_mask
-        construct_prediction_sparse(
-            np.array([local_spike]),
-            x[:, :, 0],
-            sparse_templates_array,
-            template_sparsity_mask,
-            chan_sparsity_mask,
-            nbefore,
-            True,
-        )
+            spikes["cluster_index"][i] = cluster_index
+            spikes["amplitude"][i] = amplitude
 
-        j = 1
-        for i in range(neighbors_spikes.size):
-            if mask_not_fitted[i]:
-                # add to one regressor
-                construct_prediction_sparse(
-                    local_neighbors_spikes[i : i + 1],
-                    x[:, :, j],
-                    sparse_templates_array,
-                    template_sparsity_mask,
-                    chan_sparsity_mask,
-                    nbefore,
-                    True,
-                )
-                j += 1
-            elif local_neighbors_spikes[neighbors_spikes[i]]["sample_index"] >= 0:
-                # remove from traces
-                construct_prediction_sparse(
-                    local_neighbors_spikes[i : i + 1],
-                    local_traces,
-                    sparse_templates_array,
-                    template_sparsity_mask,
-                    chan_sparsity_mask,
-                    nbefore,
-                    False,
-                )
-            # else:
-            #     pass
-
-        x = x.reshape(-1, num_spikes_to_fit)
-        y = local_traces.flatten()
-
-        res = scipy.linalg.lstsq(x, y, cond=None, lapack_driver="gelsd")
-        amplitudes = res[0]
-        amplitude = amplitudes[0]
-
-        # import matplotlib.pyplot as plt
-        # x_plot = x.reshape((lim1 - lim0, num_chans, num_spikes_to_fit)).swapaxes(0, 1).reshape(-1, num_spikes_to_fit)
-        # pred = x @ amplitudes
-        # pred_plot = pred.reshape(-1, num_chans).T.flatten()
-        # y_plot = y.reshape(-1, num_chans).T.flatten()
-        # fig, ax = plt.subplots()
-        # ax.plot(x_plot, color='b')
-        # print(x_plot.shape, y_plot.shape)
-        # ax.plot(y_plot, color='g')
-        # ax.plot(pred_plot , color='r')
-        # ax.set_title(f"{amplitudes}")
-        # # ax.set_title(f"{amplitudes} {amp_dot}")
-        # plt.show()
-
-    return amplitude
+        return spikes
 
 
 if HAVE_NUMBA:
 
     @jit(nopython=True)
-    def construct_prediction_sparse(
-        spikes, traces, sparse_templates_array, template_sparsity_mask, wanted_channel_mask, nbefore, additive
-    ):
-        #  must have np.sum(wanted_channel_mask) == traces.shape[0]
-        total_chans = wanted_channel_mask.shape[0]
-        for spike in spikes:
-            ind0 = spike["sample_index"] - nbefore
-            ind1 = ind0 + sparse_templates_array.shape[1]
-            cluster_index = spike["cluster_index"]
-            amplitude = spike["amplitude"]
-            chan_in_template = 0
-            chan_in_trace = 0
-            for chan in range(total_chans):
-                if wanted_channel_mask[chan]:
-                    if template_sparsity_mask[cluster_index, chan]:
-                        if additive:
-                            traces[ind0:ind1, chan_in_trace] += (
-                                sparse_templates_array[cluster_index, :, chan_in_template] * amplitude
-                            )
-                        else:
-                            traces[ind0:ind1, chan_in_trace] -= (
-                                sparse_templates_array[cluster_index, :, chan_in_template] * amplitude
-                            )
-                        chan_in_template += 1
-                    chan_in_trace += 1
-                else:
-                    if template_sparsity_mask[cluster_index, chan]:
-                        chan_in_template += 1
-
-    @jit(nopython=True)
-    def numba_sparse_distance(
-        wf, sparse_templates_array, template_sparsity_mask, wanted_channel_mask, possible_clusters
-    ):
+    def numba_sparse_dist(wf, templates, union_channels, possible_clusters):
         """
         numba implementation that compute distance from template with sparsity
-
-        wf is dense
-        sparse_templates_array is sparse with the template_sparsity_mask
+        handle by two separate vectors
         """
-        width, total_chans = wf.shape
+        total_cluster, width, num_chan = templates.shape
         num_cluster = possible_clusters.shape[0]
         distances = np.zeros((num_cluster,), dtype=np.float32)
         for i in prange(num_cluster):
             cluster_index = possible_clusters[i]
             sum_dist = 0.0
-            chan_in_template = 0
-            for chan in range(total_chans):
-                if wanted_channel_mask[chan]:
-                    if template_sparsity_mask[cluster_index, chan]:
-                        for s in range(width):
-                            v = wf[s, chan]
-                            t = sparse_templates_array[cluster_index, s, chan_in_template]
-                            sum_dist += (v - t) ** 2
-                        chan_in_template += 1
-                    else:
-                        for s in range(width):
-                            v = wf[s, chan]
-                            t = 0
-                            sum_dist += (v - t) ** 2
-                else:
-                    if template_sparsity_mask[cluster_index, chan]:
-                        chan_in_template += 1
+            for chan_ind in range(num_chan):
+                if union_channels[chan_ind]:
+                    for s in range(width):
+                        v = wf[s, chan_ind]
+                        t = templates[cluster_index, s, chan_ind]
+                        sum_dist += (v - t) ** 2
             distances[i] = sum_dist
         return distances
 
     @jit(nopython=True)
-    def numba_best_shift_sparse(
-        traces, sparse_template, sample_index, nbefore, possible_shifts, distances_shift, chan_sparsity
-    ):
+    def numba_best_shift(traces, template, sample_index, nbefore, possible_shifts, distances_shift, chan_sparsity):
         """
         numba implementation to compute several sample shift before template substraction
         """
-        width = sparse_template.shape[0]
-        total_chans = traces.shape[1]
+        width, num_chan = template.shape
         n_shift = possible_shifts.size
         for i in range(n_shift):
             shift = possible_shifts[i]
             sum_dist = 0.0
-            chan_in_template = 0
-            for chan in range(total_chans):
-                if chan_sparsity[chan]:
+            for chan_ind in range(num_chan):
+                if chan_sparsity[chan_ind]:
                     for s in range(width):
-                        v = traces[sample_index - nbefore + s + shift, chan]
-                        t = sparse_template[s, chan_in_template]
+                        v = traces[sample_index - nbefore + s + shift, chan_ind]
+                        t = template[s, chan_ind]
                         sum_dist += (v - t) ** 2
-                    chan_in_template += 1
             distances_shift[i] = sum_dist
 
         return distances_shift
diff --git a/src/spikeinterface/sortingcomponents/matching/wobble.py b/src/spikeinterface/sortingcomponents/matching/wobble.py
index 59e171fe5..2531a922d 100644
--- a/src/spikeinterface/sortingcomponents/matching/wobble.py
+++ b/src/spikeinterface/sortingcomponents/matching/wobble.py
@@ -8,15 +8,6 @@ from typing import List, Tuple, Optional
 from .base import BaseTemplateMatching, _base_matching_dtype
 from spikeinterface.core.template import Templates
 
-try:
-    import torch
-    import torch.nn.functional as F
-
-    HAVE_TORCH = True
-    from torch.nn.functional import conv1d
-except ImportError:
-    HAVE_TORCH = False
-
 
 @dataclass
 class WobbleParameters:
@@ -50,10 +41,6 @@ class WobbleParameters:
         Maximum value for ampltiude scaling of templates.
     scale_amplitudes : bool
         If True, scale amplitudes of templates to match spikes.
-    engine : string in ["numpy", "torch", "auto"]. Default "auto"
-        The engine to use for the convolutions
-    torch_device : string in ["cpu", "cuda", None]. Default "cpu"
-        Controls torch device if the torch engine is selected
 
     Notes
     -----
@@ -75,8 +62,6 @@ class WobbleParameters:
     scale_min: float = 0
     scale_max: float = np.inf
     scale_amplitudes: bool = False
-    engine: str = "numpy"
-    torch_device: str = "cpu"
 
     def __post_init__(self):
         assert self.amplitude_variance >= 0, "amplitude_variance must be a non-negative scalar"
@@ -359,31 +344,14 @@ class WobbleMatch(BaseTemplateMatching):
         parents=None,
         templates=None,
         parameters={},
-        engine="numpy",
-        torch_device="cpu",
     ):
 
         BaseTemplateMatching.__init__(self, recording, templates, return_output=True, parents=None)
 
-        templates_array = templates.get_dense_templates().astype(np.float32)
+        templates_array = templates.get_dense_templates().astype(np.float32, casting="safe")
 
         # Aggregate useful parameters/variables for handy access in downstream functions
         params = WobbleParameters(**parameters)
-
-        assert engine in ["numpy", "torch", "auto"], "engine should be numpy, torch or auto"
-        if engine == "auto":
-            if HAVE_TORCH:
-                self.engine = "torch"
-            else:
-                self.engine = "numpy"
-        else:
-            if engine == "torch":
-                assert HAVE_TORCH, "please install torch to use the torch engine"
-            self.engine = engine
-
-        assert torch_device in ["cuda", "cpu", None]
-        self.torch_device = torch_device
-
         template_meta = TemplateMetadata.from_parameters_and_templates(params, templates_array)
         if not templates.are_templates_sparse():
             sparsity = WobbleSparsity.from_parameters_and_templates(params, templates_array)
@@ -398,21 +366,13 @@ class WobbleMatch(BaseTemplateMatching):
         pairwise_convolution = convolve_templates(
             compressed_templates, params.jitter_factor, params.approx_rank, template_meta.jittered_indices, sparsity
         )
-
         norm_squared = compute_template_norm(sparsity.visible_channels, templates_array)
-
-        spatial = np.moveaxis(spatial, [0, 1, 2], [1, 0, 2])
-        temporal = np.moveaxis(temporal, [0, 1, 2], [1, 2, 0])
-        singular = singular.T[:, :, np.newaxis]
-
-        compressed_templates = (temporal, singular, spatial, temporal_jittered)
         template_data = TemplateData(
             compressed_templates=compressed_templates,
             pairwise_convolution=pairwise_convolution,
             norm_squared=norm_squared,
         )
 
-        self.is_pushed = False
         self.params = params
         self.template_meta = template_meta
         self.sparsity = sparsity
@@ -424,24 +384,11 @@ class WobbleMatch(BaseTemplateMatching):
         # self.margin = int(buffer_ms*1e-3 * recording.sampling_frequency)
         self.margin = 300  # To ensure equivalence with spike-psvae version of the algorithm
 
-    def _push_to_torch(self):
-        if self.engine == "torch":
-            temporal, singular, spatial, temporal_jittered = self.template_data.compressed_templates
-            spatial = torch.as_tensor(spatial, device=self.torch_device)
-            singular = torch.as_tensor(singular, device=self.torch_device)
-            temporal = torch.as_tensor(temporal.copy(), device=self.torch_device).swapaxes(0, 1)
-            temporal = torch.flip(temporal, (2,))
-            self.template_data.compressed_templates = (temporal, singular, spatial, temporal_jittered)
-        self.is_pushed = True
-
     def get_trace_margin(self):
         return self.margin
 
     def compute_matching(self, traces, start_frame, end_frame, segment_index):
 
-        if not self.is_pushed:
-            self._push_to_torch()
-
         # Unpack method_kwargs
         # nbefore, nafter = method_kwargs["nbefore"], method_kwargs["nafter"]
         # template_meta = method_kwargs["template_meta"]
@@ -453,9 +400,7 @@ class WobbleMatch(BaseTemplateMatching):
         assert traces.dtype == np.float32, "traces must be specified as np.float32"
 
         # Compute objective
-        objective = compute_objective(
-            traces, self.template_data, self.params.approx_rank, self.engine, self.torch_device
-        )
+        objective = compute_objective(traces, self.template_data, self.params.approx_rank)
         objective_normalized = 2 * objective - self.template_data.norm_squared[:, np.newaxis]
 
         # Compute spike train
@@ -841,11 +786,10 @@ def compress_templates(templates, approx_rank) -> tuple[np.ndarray, np.ndarray,
     temporal, singular, spatial = np.linalg.svd(templates, full_matrices=False)
 
     # Keep only the strongest components
-    temporal = temporal[:, :, :approx_rank].astype(np.float32)
+    temporal = temporal[:, :, :approx_rank]
     temporal = np.flip(temporal, axis=1)
-    singular = singular[:, :approx_rank].astype(np.float32)
-    spatial = spatial[:, :approx_rank, :].astype(np.float32)
-
+    singular = singular[:, :approx_rank]
+    spatial = spatial[:, :approx_rank, :]
     return temporal, singular, spatial
 
 
@@ -883,6 +827,7 @@ def upsample_and_jitter(temporal, jitter_factor, num_samples):
 
     shape_temporal_jittered = (-1, num_samples, approx_rank)
     temporal_jittered = np.reshape(temporal_jittered[:, shifted_index, :], shape_temporal_jittered)
+
     temporal_jittered = np.flip(temporal_jittered, axis=1)
     return temporal_jittered
 
@@ -944,7 +889,7 @@ def convolve_templates(compressed_templates, jitter_factor, approx_rank, jittere
     return pairwise_convolution
 
 
-def compute_objective(traces, template_data, approx_rank, engine="numpy", torch_device=None) -> np.ndarray:
+def compute_objective(traces, template_data, approx_rank) -> np.ndarray:
     """Compute objective by convolving templates with voltage traces.
 
     Parameters
@@ -953,39 +898,31 @@ def compute_objective(traces, template_data, approx_rank, engine="numpy", torch_
         Voltage traces for a chunk of the recording.
     template_data : TemplateData
         Dataclass object for aggregating template data together.
+    approx_rank : int
+        Rank of the compressed template matrices.
 
     Returns
     -------
     objective : ndarray (template_meta.num_templates, traces.shape[0]+template_meta.num_samples-1)
             Template matching objective for each template.
     """
-    temporal, singular, spatial, _ = template_data.compressed_templates
-    if engine == "torch":
-        nt = temporal.shape[2] - 1
-        num_channels = traces.shape[1]
-        blank = np.zeros((nt, num_channels), dtype=np.float32)
-        traces = np.vstack((blank, traces, blank))
-        torch_traces = torch.as_tensor(traces.T[None, :, :], device=torch_device)
-        num_templates, num_channels = temporal.shape[0], temporal.shape[1]
-        num_timesteps = torch_traces.shape[2]
-        spatially_filtered_data = torch.matmul(spatial, torch_traces)
-        scaled_filtered_data = (spatially_filtered_data * singular).swapaxes(0, 1)
-        scaled_filtered_data_ = scaled_filtered_data.reshape(1, num_templates * num_channels, num_timesteps)
-        objective = conv1d(scaled_filtered_data_, temporal, groups=num_templates, padding="valid")
-        objective = objective.cpu().numpy()[0, :, :]
-    elif engine == "numpy":
-        num_channels, num_templates = temporal.shape[0], temporal.shape[1]
-        num_timesteps = temporal.shape[2]
-        objective_len = get_convolution_len(traces.shape[0], num_timesteps)
-        conv_shape = (num_templates, objective_len)
-        objective = np.zeros(conv_shape, dtype=np.float32)
-        # Filter using overlap-and-add convolution
-        spatially_filtered_data = np.matmul(spatial, traces.T[np.newaxis, :, :])
-        scaled_filtered_data = spatially_filtered_data * singular
-        from scipy import signal
+    temporal, singular, spatial, temporal_jittered = template_data.compressed_templates
+    num_templates = temporal.shape[0]
+    num_samples = temporal.shape[1]
+    objective_len = get_convolution_len(traces.shape[0], num_samples)
+    conv_shape = (num_templates, objective_len)
+    objective = np.zeros(conv_shape, dtype=np.float32)
+    spatial_filters = np.moveaxis(spatial[:, :approx_rank, :], [0, 1, 2], [1, 0, 2])
+    temporal_filters = np.moveaxis(temporal[:, :, :approx_rank], [0, 1, 2], [1, 2, 0])
+    singular_filters = singular.T[:, :, np.newaxis]
+
+    # Filter using overlap-and-add convolution
+    spatially_filtered_data = np.matmul(spatial_filters, traces.T[np.newaxis, :, :])
+    scaled_filtered_data = spatially_filtered_data * singular_filters
+    from scipy import signal
 
-        objective_by_rank = signal.oaconvolve(scaled_filtered_data, temporal, axes=2, mode="full")
-        objective += np.sum(objective_by_rank, axis=0)
+    objective_by_rank = signal.oaconvolve(scaled_filtered_data, temporal_filters, axes=2, mode="full")
+    objective += np.sum(objective_by_rank, axis=0)
     return objective
 
 
diff --git a/src/spikeinterface/sortingcomponents/motion/dredge.py b/src/spikeinterface/sortingcomponents/motion/dredge.py
index bfedd4e1e..e2b6b1a2b 100644
--- a/src/spikeinterface/sortingcomponents/motion/dredge.py
+++ b/src/spikeinterface/sortingcomponents/motion/dredge.py
@@ -22,19 +22,20 @@ but the original function dredge_ap() and dredge_online_lfp() can be used direct
 
 """
 
-import gc
 import warnings
 
-import numpy as np
 from tqdm.auto import trange
+import numpy as np
+
+import gc
 
 from .motion_utils import (
     Motion,
-    get_spatial_bin_edges,
     get_spatial_windows,
     get_window_domains,
-    make_2d_motion_histogram,
     scipy_conv1d,
+    make_2d_motion_histogram,
+    get_spatial_bin_edges,
 )
 
 
@@ -978,7 +979,7 @@ def xcorr_windows(
 
     if max_disp_um is None:
         if rigid:
-            max_disp_um = int(np.ptp(spatial_bin_edges_um) // 4)
+            max_disp_um = int(spatial_bin_edges_um.ptp() // 4)
         else:
             max_disp_um = int(win_scale_um // 4)
 
diff --git a/src/spikeinterface/sortingcomponents/peak_detection.py b/src/spikeinterface/sortingcomponents/peak_detection.py
index d03744f8f..d608c5d10 100644
--- a/src/spikeinterface/sortingcomponents/peak_detection.py
+++ b/src/spikeinterface/sortingcomponents/peak_detection.py
@@ -57,7 +57,6 @@ def detect_peaks(
     folder=None,
     names=None,
     skip_after_n_peaks=None,
-    recording_slices=None,
     **kwargs,
 ):
     """Peak detection based on threshold crossing in term of k x MAD.
@@ -84,10 +83,6 @@ def detect_peaks(
     skip_after_n_peaks : None | int
         Skip the computation after n_peaks.
         This is not an exact because internally this skip is done per worker in average.
-    recording_slices : None | list[tuple]
-        Optionaly give a list of slices to run the pipeline only on some chunks of the recording.
-        It must be a list of (segment_index, frame_start, frame_stop).
-        If None (default), the function iterates over the entire duration of the recording.
 
     {method_doc}
     {job_doc}
@@ -140,7 +135,6 @@ def detect_peaks(
         folder=folder,
         names=names,
         skip_after_n_peaks=skip_after_n_peaks,
-        recording_slices=recording_slices,
     )
     return outs
 
@@ -609,13 +603,13 @@ class DetectPeakMatchedFiltering(PeakDetector):
     params_doc = (
         DetectPeakByChannel.params_doc
         + """
-    radius_um : float
+    radius_um: float
         The radius to use to select neighbour channels for locally exclusive detection.
-    prototype : array
+    prototype: array
         The canonical waveform of action potentials
-    ms_before : float
-        The time in ms before the maximial value of the absolute prototype
-    weight_method : dict
+    rank : int (default 1)
+        The rank for SVD convolution of spatiotemporal templates with the traces
+    weight_method: dict
         Parameter that should be provided to the get_convolution_weights() function
         in order to know how to estimate the positions. One argument is mode that could
         be either gaussian_2d (KS like) or exponential_3d (default)
@@ -631,12 +625,12 @@ class DetectPeakMatchedFiltering(PeakDetector):
         detect_threshold=5,
         exclude_sweep_ms=0.1,
         radius_um=50,
+        rank=1,
         noise_levels=None,
         random_chunk_kwargs={"num_chunks_per_segment": 5},
         weight_method={},
     ):
         PeakDetector.__init__(self, recording, return_output=True)
-        from scipy.sparse import csr_matrix
 
         if not HAVE_NUMBA:
             raise ModuleNotFoundError('matched_filtering" needs numba which is not installed')
@@ -670,7 +664,7 @@ class DetectPeakMatchedFiltering(PeakDetector):
             self.num_templates *= 2
 
         self.weights = self.weights.reshape(self.num_templates * self.num_z_factors, -1)
-        self.weights = csr_matrix(self.weights)
+
         random_data = get_random_data_chunks(recording, return_scaled=False, **random_chunk_kwargs)
         conv_random_data = self.get_convolved_traces(random_data)
         medians = np.median(conv_random_data, axis=1)
@@ -740,10 +734,10 @@ class DetectPeakMatchedFiltering(PeakDetector):
         return (local_peaks,)
 
     def get_convolved_traces(self, traces):
-        from scipy.signal import oaconvolve
+        import scipy.signal
 
-        tmp = oaconvolve(self.prototype[None, :], traces.T, axes=1, mode="valid")
-        scalar_products = self.weights.dot(tmp)
+        tmp = scipy.signal.oaconvolve(self.prototype[None, :], traces.T, axes=1, mode="valid")
+        scalar_products = np.dot(self.weights, tmp)
         return scalar_products
 
 
diff --git a/src/spikeinterface/widgets/autocorrelograms.py b/src/spikeinterface/widgets/autocorrelograms.py
index c211a277f..c8acd93dc 100644
--- a/src/spikeinterface/widgets/autocorrelograms.py
+++ b/src/spikeinterface/widgets/autocorrelograms.py
@@ -9,13 +9,7 @@ class AutoCorrelogramsWidget(CrossCorrelogramsWidget):
     # the doc is copied form CrossCorrelogramsWidget
 
     def __init__(self, *args, **kargs):
-        _ = kargs.pop("min_similarity_for_correlograms", 0.0)
-        CrossCorrelogramsWidget.__init__(
-            self,
-            *args,
-            **kargs,
-            min_similarity_for_correlograms=None,
-        )
+        CrossCorrelogramsWidget.__init__(self, *args, **kargs)
 
     def plot_matplotlib(self, data_plot, **backend_kwargs):
         import matplotlib.pyplot as plt
diff --git a/src/spikeinterface/widgets/crosscorrelograms.py b/src/spikeinterface/widgets/crosscorrelograms.py
index 88dd80332..cdb2041aa 100644
--- a/src/spikeinterface/widgets/crosscorrelograms.py
+++ b/src/spikeinterface/widgets/crosscorrelograms.py
@@ -21,8 +21,7 @@ class CrossCorrelogramsWidget(BaseWidget):
         List of unit ids
     min_similarity_for_correlograms : float, default: 0.2
         For sortingview backend. Threshold for computing pair-wise cross-correlograms.
-        If template similarity between two units is below this threshold, the cross-correlogram is not displayed.
-        For auto-correlograms plot, this is automatically set to None.
+        If template similarity between two units is below this threshold, the cross-correlogram is not displayed
     window_ms : float, default: 100.0
         Window for CCGs in ms. If correlograms are already computed (e.g. with SortingAnalyzer),
         this argument is ignored
diff --git a/src/spikeinterface/widgets/traces.py b/src/spikeinterface/widgets/traces.py
index f5dadc780..86f2350a8 100644
--- a/src/spikeinterface/widgets/traces.py
+++ b/src/spikeinterface/widgets/traces.py
@@ -52,8 +52,6 @@ class TracesWidget(BaseWidget):
         If dict, keys should be the same as recording keys
     scale : float, default: 1
         Scale factor for the traces
-    vspacing_factor : float, default: 1.5
-        Vertical spacing between channels as a multiple of maximum channel amplitude
     with_colorbar : bool, default: True
         When mode is "map", a colorbar is added
     tile_size : int, default: 1500
@@ -84,7 +82,6 @@ class TracesWidget(BaseWidget):
         tile_size=1500,
         seconds_per_row=0.2,
         scale=1,
-        vspacing_factor=1.5,
         with_colorbar=True,
         add_legend=True,
         backend=None,
@@ -171,7 +168,7 @@ class TracesWidget(BaseWidget):
         traces0 = list_traces[0]
         mean_channel_std = np.mean(np.std(traces0, axis=0))
         max_channel_amp = np.max(np.max(np.abs(traces0), axis=0))
-        vspacing = max_channel_amp * vspacing_factor
+        vspacing = max_channel_amp * 1.5
 
         if rec0.get_channel_groups() is None:
             color_groups = False
diff --git a/src/spikeinterface/widgets/unit_summary.py b/src/spikeinterface/widgets/unit_summary.py
index 946611011..755e60ccb 100644
--- a/src/spikeinterface/widgets/unit_summary.py
+++ b/src/spikeinterface/widgets/unit_summary.py
@@ -107,90 +107,82 @@ class UnitSummaryWidget(BaseWidget):
         # and use custum grid spec
         fig = self.figure
         nrows = 2
-        ncols = 2
-        if sorting_analyzer.has_extension("correlograms"):
-            ncols += 1
-        if sorting_analyzer.has_extension("waveforms"):
+        ncols = 3
+        if sorting_analyzer.has_extension("correlograms") or sorting_analyzer.has_extension("spike_amplitudes"):
             ncols += 1
         if sorting_analyzer.has_extension("spike_amplitudes"):
             nrows += 1
         gs = fig.add_gridspec(nrows, ncols)
-        col_counter = 0
 
-        # Unit locations and unit waveform plots are always generated
-        ax_unit_locations = fig.add_subplot(gs[:2, col_counter])
-        _ = UnitLocationsWidget(
-            sorting_analyzer,
-            unit_ids=[unit_id],
-            unit_colors=unit_colors,
-            plot_legend=False,
-            backend="matplotlib",
-            ax=ax_unit_locations,
-            **unitlocationswidget_kwargs,
-        )
-        col_counter += 1
-
-        unit_locations = sorting_analyzer.get_extension("unit_locations").get_data(outputs="by_unit")
-        unit_location = unit_locations[unit_id]
-        x, y = unit_location[0], unit_location[1]
-        ax_unit_locations.set_xlim(x - 80, x + 80)
-        ax_unit_locations.set_ylim(y - 250, y + 250)
-        ax_unit_locations.set_xticks([])
-        ax_unit_locations.set_xlabel(None)
-        ax_unit_locations.set_ylabel(None)
-
-        ax_unit_waveforms = fig.add_subplot(gs[:2, col_counter])
-        _ = UnitWaveformsWidget(
+        if sorting_analyzer.has_extension("unit_locations"):
+            ax1 = fig.add_subplot(gs[:2, 0])
+            # UnitLocationsPlotter().do_plot(dp.plot_data_unit_locations, ax=ax1)
+            w = UnitLocationsWidget(
+                sorting_analyzer,
+                unit_ids=[unit_id],
+                unit_colors=unit_colors,
+                plot_legend=False,
+                backend="matplotlib",
+                ax=ax1,
+                **unitlocationswidget_kwargs,
+            )
+
+            unit_locations = sorting_analyzer.get_extension("unit_locations").get_data(outputs="by_unit")
+            unit_location = unit_locations[unit_id]
+            x, y = unit_location[0], unit_location[1]
+            ax1.set_xlim(x - 80, x + 80)
+            ax1.set_ylim(y - 250, y + 250)
+            ax1.set_xticks([])
+            ax1.set_xlabel(None)
+            ax1.set_ylabel(None)
+
+        ax2 = fig.add_subplot(gs[:2, 1])
+        w = UnitWaveformsWidget(
             sorting_analyzer,
             unit_ids=[unit_id],
             unit_colors=unit_colors,
             plot_templates=True,
-            plot_waveforms=sorting_analyzer.has_extension("waveforms"),
             same_axis=True,
             plot_legend=False,
             sparsity=sparsity,
             backend="matplotlib",
-            ax=ax_unit_waveforms,
+            ax=ax2,
             **unitwaveformswidget_kwargs,
         )
-        col_counter += 1
 
-        ax_unit_waveforms.set_title(None)
+        ax2.set_title(None)
 
-        if sorting_analyzer.has_extension("waveforms"):
-            ax_waveform_density = fig.add_subplot(gs[:2, col_counter])
-            UnitWaveformDensityMapWidget(
-                sorting_analyzer,
-                unit_ids=[unit_id],
-                unit_colors=unit_colors,
-                use_max_channel=True,
-                same_axis=False,
-                backend="matplotlib",
-                ax=ax_waveform_density,
-                **unitwaveformdensitymapwidget_kwargs,
-            )
-            col_counter += 1
-            ax_waveform_density.set_ylabel(None)
+        ax3 = fig.add_subplot(gs[:2, 2])
+        UnitWaveformDensityMapWidget(
+            sorting_analyzer,
+            unit_ids=[unit_id],
+            unit_colors=unit_colors,
+            use_max_channel=True,
+            same_axis=False,
+            backend="matplotlib",
+            ax=ax3,
+            **unitwaveformdensitymapwidget_kwargs,
+        )
+        ax3.set_ylabel(None)
 
         if sorting_analyzer.has_extension("correlograms"):
-            ax_correlograms = fig.add_subplot(gs[:2, col_counter])
+            ax4 = fig.add_subplot(gs[:2, 3])
             AutoCorrelogramsWidget(
                 sorting_analyzer,
                 unit_ids=[unit_id],
                 unit_colors=unit_colors,
                 backend="matplotlib",
-                ax=ax_correlograms,
+                ax=ax4,
                 **autocorrelogramswidget_kwargs,
             )
-            col_counter += 1
 
-            ax_correlograms.set_title(None)
-            ax_correlograms.set_yticks([])
+            ax4.set_title(None)
+            ax4.set_yticks([])
 
         if sorting_analyzer.has_extension("spike_amplitudes"):
-            ax_spike_amps = fig.add_subplot(gs[2, : col_counter - 1])
-            ax_amps_distribution = fig.add_subplot(gs[2, col_counter - 1])
-            axes = np.array([ax_spike_amps, ax_amps_distribution])
+            ax5 = fig.add_subplot(gs[2, :3])
+            ax6 = fig.add_subplot(gs[2, 3])
+            axes = np.array([ax5, ax6])
             AmplitudesWidget(
                 sorting_analyzer,
                 unit_ids=[unit_id],


Test patch:
diff --git a/src/spikeinterface/core/tests/test_analyzer_extension_core.py b/src/spikeinterface/core/tests/test_analyzer_extension_core.py
index 6f5bef3c6..626899ab6 100644
--- a/src/spikeinterface/core/tests/test_analyzer_extension_core.py
+++ b/src/spikeinterface/core/tests/test_analyzer_extension_core.py
@@ -2,8 +2,6 @@ import pytest
 
 import shutil
 
-from pathlib import Path
-
 from spikeinterface.core import generate_ground_truth_recording
 from spikeinterface.core import create_sorting_analyzer
 from spikeinterface.core import Templates
@@ -252,17 +250,16 @@ def test_compute_several(create_cache_folder):
 
 
 if __name__ == "__main__":
-    cache_folder = Path(__file__).resolve().parents[4] / "cache_folder" / "core"
-    # test_ComputeWaveforms(format="memory", sparse=True, create_cache_folder=cache_folder)
-    # test_ComputeWaveforms(format="memory", sparse=False, create_cache_folder=cache_folder)
-    # test_ComputeWaveforms(format="binary_folder", sparse=True, create_cache_folder=cache_folder)
-    # test_ComputeWaveforms(format="binary_folder", sparse=False, create_cache_folder=cache_folder)
-    # test_ComputeWaveforms(format="zarr", sparse=True, create_cache_folder=cache_folder)
-    # test_ComputeWaveforms(format="zarr", sparse=False, create_cache_folder=cache_folder)
-    # test_ComputeRandomSpikes(format="memory", sparse=True, create_cache_folder=cache_folder)
-    test_ComputeRandomSpikes(format="binary_folder", sparse=False, create_cache_folder=cache_folder)
-    test_ComputeTemplates(format="memory", sparse=True, create_cache_folder=cache_folder)
-    test_ComputeNoiseLevels(format="memory", sparse=False, create_cache_folder=cache_folder)
-
-    # test_get_children_dependencies()
-    # test_delete_on_recompute(cache_folder)
+
+    test_ComputeWaveforms(format="memory", sparse=True)
+    test_ComputeWaveforms(format="memory", sparse=False)
+    test_ComputeWaveforms(format="binary_folder", sparse=True)
+    test_ComputeWaveforms(format="binary_folder", sparse=False)
+    test_ComputeWaveforms(format="zarr", sparse=True)
+    test_ComputeWaveforms(format="zarr", sparse=False)
+    test_ComputeRandomSpikes(format="memory", sparse=True)
+    test_ComputeTemplates(format="memory", sparse=True)
+    test_ComputeNoiseLevels(format="memory", sparse=False)
+
+    test_get_children_dependencies()
+    test_delete_on_recompute()
diff --git a/src/spikeinterface/core/tests/test_node_pipeline.py b/src/spikeinterface/core/tests/test_node_pipeline.py
index 028eaecf1..deef2291c 100644
--- a/src/spikeinterface/core/tests/test_node_pipeline.py
+++ b/src/spikeinterface/core/tests/test_node_pipeline.py
@@ -4,7 +4,7 @@ from pathlib import Path
 import shutil
 
 from spikeinterface import create_sorting_analyzer, get_template_extremum_channel, generate_ground_truth_recording
-from spikeinterface.core.job_tools import divide_recording_into_chunks
+
 
 # from spikeinterface.sortingcomponents.peak_detection import detect_peaks
 from spikeinterface.core.node_pipeline import (
@@ -191,8 +191,8 @@ def test_run_node_pipeline(cache_folder_creation):
             unpickled_node = pickle.loads(pickled_node)
 
 
-def test_skip_after_n_peaks_and_recording_slices():
-    recording, sorting = generate_ground_truth_recording(num_channels=10, num_units=10, durations=[10.0], seed=2205)
+def test_skip_after_n_peaks():
+    recording, sorting = generate_ground_truth_recording(num_channels=10, num_units=10, durations=[10.0])
 
     # job_kwargs = dict(chunk_duration="0.5s", n_jobs=2, progress_bar=False)
     job_kwargs = dict(chunk_duration="0.5s", n_jobs=1, progress_bar=False)
@@ -211,27 +211,18 @@ def test_skip_after_n_peaks_and_recording_slices():
     node1 = AmplitudeExtractionNode(recording, parents=[node0], param0=6.6, return_output=True)
     nodes = [node0, node1]
 
-    # skip
     skip_after_n_peaks = 30
     some_amplitudes = run_node_pipeline(
         recording, nodes, job_kwargs, gather_mode="memory", skip_after_n_peaks=skip_after_n_peaks
     )
+
     assert some_amplitudes.size >= skip_after_n_peaks
     assert some_amplitudes.size < spikes.size
 
-    # slices : 1 every 4
-    recording_slices = divide_recording_into_chunks(recording, 10_000)
-    recording_slices = recording_slices[::4]
-    some_amplitudes = run_node_pipeline(
-        recording, nodes, job_kwargs, gather_mode="memory", recording_slices=recording_slices
-    )
-    tolerance = 1.2
-    assert some_amplitudes.size < (spikes.size // 4) * tolerance
-
 
 # the following is for testing locally with python or ipython. It is not used in ci or with pytest.
 if __name__ == "__main__":
     # folder = Path("./cache_folder/core")
     # test_run_node_pipeline(folder)
 
-    test_skip_after_n_peaks_and_recording_slices()
+    test_skip_after_n_peaks()
diff --git a/src/spikeinterface/core/tests/test_recording_tools.py b/src/spikeinterface/core/tests/test_recording_tools.py
index dad5273f1..23a1574f2 100644
--- a/src/spikeinterface/core/tests/test_recording_tools.py
+++ b/src/spikeinterface/core/tests/test_recording_tools.py
@@ -11,7 +11,6 @@ from spikeinterface.core.generate import NoiseGeneratorRecording
 from spikeinterface.core.recording_tools import (
     write_binary_recording,
     write_memory_recording,
-    get_random_recording_slices,
     get_random_data_chunks,
     get_chunk_with_margin,
     get_closest_channels,
@@ -168,17 +167,6 @@ def test_write_memory_recording():
         shm.unlink()
 
 
-def test_get_random_recording_slices():
-    rec = generate_recording(num_channels=1, sampling_frequency=1000.0, durations=[10.0, 20.0])
-    rec_slices = get_random_recording_slices(
-        rec, method="full_random", num_chunks_per_segment=20, chunk_duration="500ms", margin_frames=0, seed=0
-    )
-    assert len(rec_slices) == 40
-    for seg_ind, start, stop in rec_slices:
-        assert stop - start == 500
-        assert seg_ind in (0, 1)
-
-
 def test_get_random_data_chunks():
     rec = generate_recording(num_channels=1, sampling_frequency=1000.0, durations=[10.0, 20.0])
     chunks = get_random_data_chunks(rec, num_chunks_per_segment=50, chunk_size=500, seed=0)
@@ -194,17 +182,16 @@ def test_get_closest_channels():
 
 
 def test_get_noise_levels():
-    job_kwargs = dict(n_jobs=1, progress_bar=True)
     rec = generate_recording(num_channels=2, sampling_frequency=1000.0, durations=[60.0])
 
-    noise_levels_1 = get_noise_levels(rec, return_scaled=False, **job_kwargs)
-    noise_levels_2 = get_noise_levels(rec, return_scaled=False, **job_kwargs)
+    noise_levels_1 = get_noise_levels(rec, return_scaled=False)
+    noise_levels_2 = get_noise_levels(rec, return_scaled=False)
 
     rec.set_channel_gains(0.1)
     rec.set_channel_offsets(0)
-    noise_levels = get_noise_levels(rec, return_scaled=True, force_recompute=True, **job_kwargs)
+    noise_levels = get_noise_levels(rec, return_scaled=True, force_recompute=True)
 
-    noise_levels = get_noise_levels(rec, return_scaled=True, method="std", **job_kwargs)
+    noise_levels = get_noise_levels(rec, return_scaled=True, method="std")
 
     # Generate a recording following a gaussian distribution to check the result of get_noise.
     std = 6.0
@@ -214,10 +201,8 @@ def test_get_noise_levels():
     recording = NumpyRecording(traces, 30000)
 
     assert np.all(noise_levels_1 == noise_levels_2)
-    assert np.allclose(get_noise_levels(recording, return_scaled=False, **job_kwargs), [std, std], rtol=1e-2, atol=1e-3)
-    assert np.allclose(
-        get_noise_levels(recording, method="std", return_scaled=False, **job_kwargs), [std, std], rtol=1e-2, atol=1e-3
-    )
+    assert np.allclose(get_noise_levels(recording, return_scaled=False), [std, std], rtol=1e-2, atol=1e-3)
+    assert np.allclose(get_noise_levels(recording, method="std", return_scaled=False), [std, std], rtol=1e-2, atol=1e-3)
 
 
 def test_get_noise_levels_output():
@@ -231,21 +216,10 @@ def test_get_noise_levels_output():
     traces = rng.normal(loc=10.0, scale=std, size=(num_samples, num_channels))
     recording = NumpyRecording(traces_list=traces, sampling_frequency=sampling_frequency)
 
-    std_estimated_with_mad = get_noise_levels(
-        recording,
-        method="mad",
-        return_scaled=False,
-        random_slices_kwargs=dict(num_chunks_per_segment=40, chunk_size=1_000, seed=seed),
-    )
-    print(std_estimated_with_mad)
+    std_estimated_with_mad = get_noise_levels(recording, method="mad", return_scaled=False, chunk_size=1_000)
     assert np.allclose(std_estimated_with_mad, [std, std], rtol=1e-2, atol=1e-3)
 
-    std_estimated_with_std = get_noise_levels(
-        recording,
-        method="std",
-        return_scaled=False,
-        random_slices_kwargs=dict(num_chunks_per_segment=40, chunk_size=1_000, seed=seed),
-    )
+    std_estimated_with_std = get_noise_levels(recording, method="std", return_scaled=False, chunk_size=1_000)
     assert np.allclose(std_estimated_with_std, [std, std], rtol=1e-2, atol=1e-3)
 
 
@@ -359,16 +333,14 @@ def test_do_recording_attributes_match():
 
 if __name__ == "__main__":
     # Create a temporary folder using the standard library
-    # import tempfile
-
-    # with tempfile.TemporaryDirectory() as tmpdirname:
-    #     tmp_path = Path(tmpdirname)
-    #     test_write_binary_recording(tmp_path)
-    # test_write_memory_recording()
-
-    test_get_random_recording_slices()
-    # test_get_random_data_chunks()
-    # test_get_closest_channels()
-    # test_get_noise_levels()
-    # test_get_noise_levels_output()
-    # test_order_channels_by_depth()
+    import tempfile
+
+    with tempfile.TemporaryDirectory() as tmpdirname:
+        tmp_path = Path(tmpdirname)
+        test_write_binary_recording(tmp_path)
+    test_write_memory_recording()
+
+    test_get_random_data_chunks()
+    test_get_closest_channels()
+    test_get_noise_levels()
+    test_order_channels_by_depth()
diff --git a/src/spikeinterface/core/tests/test_time_handling.py b/src/spikeinterface/core/tests/test_time_handling.py
index a129316ee..9b7ed11bb 100644
--- a/src/spikeinterface/core/tests/test_time_handling.py
+++ b/src/spikeinterface/core/tests/test_time_handling.py
@@ -15,7 +15,10 @@ class TestTimeHandling:
     is generated on the fly. Both time representations are tested here.
     """
 
-    # Fixtures #####
+    # #########################################################################
+    # Fixtures
+    # #########################################################################
+
     @pytest.fixture(scope="session")
     def time_vector_recording(self):
         """
@@ -95,7 +98,10 @@ class TestTimeHandling:
         raw_recording, times_recording, all_times = time_recording_fixture
         return (raw_recording, times_recording, all_times)
 
-    # Tests #####
+    # #########################################################################
+    # Tests
+    # #########################################################################
+
     def test_has_time_vector(self, time_vector_recording):
         """
         Test the `has_time_vector` function returns `False` before
@@ -305,7 +311,87 @@ class TestTimeHandling:
 
         assert np.array_equal(sorting_analyzer.get_total_duration(), raw_recording.get_total_duration())
 
-    # Helpers ####
+    @pytest.mark.parametrize("fixture_name", ["time_vector_recording", "t_start_recording"])
+    @pytest.mark.parametrize("shift", [-123.456, 123.456])
+    def test_shift_time_all_segments(self, request, fixture_name, shift):
+        """
+        Shift the times in every segment using the `None` default, then
+        check that every segment of the recording is shifted as expected.
+        """
+        _, times_recording, all_times = self._get_fixture_data(request, fixture_name)
+
+        num_segments, orig_seg_data = self._store_all_times(times_recording)
+
+        times_recording.shift_times(shift)  # use default `segment_index=None`
+
+        for idx in range(num_segments):
+            assert np.allclose(
+                orig_seg_data[idx], times_recording.get_times(segment_index=idx) - shift, rtol=0, atol=1e-8
+            )
+
+    @pytest.mark.parametrize("fixture_name", ["time_vector_recording", "t_start_recording"])
+    @pytest.mark.parametrize("shift", [-123.456, 123.456])
+    def test_shift_times_different_segments(self, request, fixture_name, shift):
+        """
+        Shift each segment separately, and check the shifted segment only
+        is shifted as expected.
+        """
+        _, times_recording, all_times = self._get_fixture_data(request, fixture_name)
+
+        num_segments, orig_seg_data = self._store_all_times(times_recording)
+
+        # For each segment, shift the segment only and check the
+        # times are updated as expected.
+        for idx in range(num_segments):
+
+            scaler = idx + 2
+            times_recording.shift_times(shift * scaler, segment_index=idx)
+
+            assert np.allclose(
+                orig_seg_data[idx], times_recording.get_times(segment_index=idx) - shift * scaler, rtol=0, atol=1e-8
+            )
+
+            # Just do a little check that we are not
+            # accidentally changing some other segments,
+            # which should remain unchanged at this point in the loop.
+            if idx != num_segments - 1:
+                assert np.array_equal(orig_seg_data[idx + 1], times_recording.get_times(segment_index=idx + 1))
+
+    @pytest.mark.parametrize("fixture_name", ["time_vector_recording", "t_start_recording"])
+    def test_save_and_load_time_shift(self, request, fixture_name, tmp_path):
+        """
+        Save the shifted data and check the shift is propagated correctly.
+        """
+        _, times_recording, all_times = self._get_fixture_data(request, fixture_name)
+
+        shift = 100
+        times_recording.shift_times(shift=shift)
+
+        times_recording.save(folder=tmp_path / "my_file")
+
+        loaded_recording = si.load_extractor(tmp_path / "my_file")
+
+        for idx in range(times_recording.get_num_segments()):
+            assert np.array_equal(
+                times_recording.get_times(segment_index=idx), loaded_recording.get_times(segment_index=idx)
+            )
+
+    def _store_all_times(self, recording):
+        """
+        Convenience function to store original times of all segments to a dict.
+        """
+        num_segments = recording.get_num_segments()
+        seg_data = {}
+
+        for idx in range(num_segments):
+            seg_data[idx] = copy.deepcopy(recording.get_times(segment_index=idx))
+
+        return num_segments, seg_data
+
+    # #########################################################################
+    # Helpers
+    # #########################################################################
+
     def _check_times_match(self, recording, all_times):
         """
         For every segment in a recording, check the `get_times()`
diff --git a/src/spikeinterface/curation/tests/test_auto_merge.py b/src/spikeinterface/curation/tests/test_auto_merge.py
index 4c05f41a4..33fd06d27 100644
--- a/src/spikeinterface/curation/tests/test_auto_merge.py
+++ b/src/spikeinterface/curation/tests/test_auto_merge.py
@@ -3,16 +3,16 @@ import pytest
 
 from spikeinterface.core import create_sorting_analyzer
 from spikeinterface.core.generate import inject_some_split_units
-from spikeinterface.curation import compute_merge_unit_groups, auto_merge
+from spikeinterface.curation import get_potential_auto_merge
 
 
 from spikeinterface.curation.tests.common import make_sorting_analyzer, sorting_analyzer_for_curation
 
 
 @pytest.mark.parametrize(
-    "preset", ["x_contaminations", "feature_neighbors", "temporal_splits", "similarity_correlograms", None]
+    "preset", ["x_contaminations", "feature_neighbors", "temporal_splits", "similarity_correlograms"]
 )
-def test_compute_merge_unit_groups(sorting_analyzer_for_curation, preset):
+def test_get_auto_merge_list(sorting_analyzer_for_curation, preset):
 
     print(sorting_analyzer_for_curation)
     sorting = sorting_analyzer_for_curation.sorting
@@ -47,38 +47,32 @@ def test_compute_merge_unit_groups(sorting_analyzer_for_curation, preset):
     )
 
     if preset is not None:
-        # do not resolve graph for checking true pairs
-        merge_unit_groups, outs = compute_merge_unit_groups(
+        potential_merges, outs = get_potential_auto_merge(
             sorting_analyzer,
             preset=preset,
-            resolve_graph=False,
-            # min_spikes=1000,
-            # max_distance_um=150.0,
-            # contamination_thresh=0.2,
-            # corr_diff_thresh=0.16,
-            # template_diff_thresh=0.25,
-            # censored_period_ms=0.0,
-            # refractory_period_ms=4.0,
-            # sigma_smooth_ms=0.6,
-            # adaptative_window_thresh=0.5,
-            # firing_contamination_balance=1.5,
+            min_spikes=1000,
+            max_distance_um=150.0,
+            contamination_thresh=0.2,
+            corr_diff_thresh=0.16,
+            template_diff_thresh=0.25,
+            censored_period_ms=0.0,
+            refractory_period_ms=4.0,
+            sigma_smooth_ms=0.6,
+            adaptative_window_thresh=0.5,
+            firing_contamination_balance=1.5,
             extra_outputs=True,
-            **job_kwargs,
         )
         if preset == "x_contaminations":
-            assert len(merge_unit_groups) == num_unit_splited
+            assert len(potential_merges) == num_unit_splited
             for true_pair in other_ids.values():
                 true_pair = tuple(true_pair)
-                assert true_pair in merge_unit_groups
+                assert true_pair in potential_merges
     else:
         # when preset is None you have to specify the steps
         with pytest.raises(ValueError):
-            merge_unit_groups = compute_merge_unit_groups(sorting_analyzer, preset=preset)
-        merge_unit_groups = compute_merge_unit_groups(
-            sorting_analyzer,
-            preset=preset,
-            steps=["num_spikes", "snr", "remove_contaminated", "unit_locations"],
-            **job_kwargs,
+            potential_merges = get_potential_auto_merge(sorting_analyzer, preset=preset)
+        potential_merges = get_potential_auto_merge(
+            sorting_analyzer, preset=preset, steps=["min_spikes", "min_snr", "remove_contaminated", "unit_positions"]
         )
 
     # DEBUG
@@ -99,7 +93,7 @@ def test_compute_merge_unit_groups(sorting_analyzer_for_curation, preset):
 
     # m = correlograms.shape[2] // 2
 
-    # for unit_id1, unit_id2 in merge_unit_groups[:5]:
+    # for unit_id1, unit_id2 in potential_merges[:5]:
     #     unit_ind1 = sorting_with_split.id_to_index(unit_id1)
     #     unit_ind2 = sorting_with_split.id_to_index(unit_id2)
 
@@ -135,6 +129,4 @@ def test_compute_merge_unit_groups(sorting_analyzer_for_curation, preset):
 
 if __name__ == "__main__":
     sorting_analyzer = make_sorting_analyzer(sparse=True)
-    # preset = "x_contaminations"
-    preset = None
-    test_compute_merge_unit_groups(sorting_analyzer, preset=preset)
+    test_get_auto_merge_list(sorting_analyzer)
diff --git a/src/spikeinterface/extractors/tests/test_neoextractors.py b/src/spikeinterface/extractors/tests/test_neoextractors.py
index 3da92331a..fcdd766f4 100644
--- a/src/spikeinterface/extractors/tests/test_neoextractors.py
+++ b/src/spikeinterface/extractors/tests/test_neoextractors.py
@@ -368,7 +368,7 @@ class Plexon2RecordingTest(RecordingCommonTestSuite, unittest.TestCase):
     ExtractorClass = Plexon2RecordingExtractor
     downloads = ["plexon"]
     entities = [
-        ("plexon/4chDemoPL2.pl2", {"stream_name": "WB-Wideband"}),
+        ("plexon/4chDemoPL2.pl2", {"stream_id": "3"}),
     ]
 
 
diff --git a/src/spikeinterface/postprocessing/tests/test_unit_locations.py b/src/spikeinterface/postprocessing/tests/test_unit_locations.py
index 545edb349..c40a917a2 100644
--- a/src/spikeinterface/postprocessing/tests/test_unit_locations.py
+++ b/src/spikeinterface/postprocessing/tests/test_unit_locations.py
@@ -13,7 +13,6 @@ class TestUnitLocationsExtension(AnalyzerExtensionCommonTestSuite):
             dict(method="grid_convolution", radius_um=150, weight_method={"mode": "gaussian_2d"}),
             dict(method="monopolar_triangulation", radius_um=150),
             dict(method="monopolar_triangulation", radius_um=150, optimizer="minimize_with_log_penality"),
-            dict(method="max_channel"),
         ],
     )
     def test_extension(self, params):
diff --git a/src/spikeinterface/preprocessing/tests/test_decimate.py b/src/spikeinterface/preprocessing/tests/test_decimate.py
index aab17560a..100972f76 100644
--- a/src/spikeinterface/preprocessing/tests/test_decimate.py
+++ b/src/spikeinterface/preprocessing/tests/test_decimate.py
@@ -8,14 +8,19 @@ from spikeinterface.preprocessing.decimate import DecimateRecording
 import numpy as np
 
 
-@pytest.mark.parametrize("num_segments", [1, 2])
-@pytest.mark.parametrize("decimation_offset", [0, 1, 5, 21, 101])
-@pytest.mark.parametrize("decimation_factor", [1, 7, 50])
-def test_decimate(num_segments, decimation_offset, decimation_factor):
-    segment_num_samps = [20000, 40000]
+@pytest.mark.parametrize("N_segments", [1, 2])
+@pytest.mark.parametrize("decimation_offset", [0, 1, 9, 10, 11, 100, 101])
+@pytest.mark.parametrize("decimation_factor", [1, 9, 10, 11, 100, 101])
+@pytest.mark.parametrize("start_frame", [0, 1, 5, None, 1000])
+@pytest.mark.parametrize("end_frame", [0, 1, 5, None, 1000])
+def test_decimate(N_segments, decimation_offset, decimation_factor, start_frame, end_frame):
+    rec = generate_recording()
+
+    segment_num_samps = [101 + i for i in range(N_segments)]
+
     rec = NumpyRecording([np.arange(2 * N).reshape(N, 2) for N in segment_num_samps], 1)
 
-    parent_traces = [rec.get_traces(i) for i in range(num_segments)]
+    parent_traces = [rec.get_traces(i) for i in range(N_segments)]
 
     if decimation_offset >= min(segment_num_samps) or decimation_offset >= decimation_factor:
         with pytest.raises(ValueError):
@@ -23,59 +28,19 @@ def test_decimate(num_segments, decimation_offset, decimation_factor):
         return
 
     decimated_rec = DecimateRecording(rec, decimation_factor, decimation_offset=decimation_offset)
-    decimated_parent_traces = [parent_traces[i][decimation_offset::decimation_factor] for i in range(num_segments)]
+    decimated_parent_traces = [parent_traces[i][decimation_offset::decimation_factor] for i in range(N_segments)]
 
-    for start_frame in [0, 1, 5, None, 1000]:
-        for end_frame in [0, 1, 5, None, 1000]:
-            if start_frame is None:
-                start_frame = max(decimated_rec.get_num_samples(i) for i in range(num_segments))
-            if end_frame is None:
-                end_frame = max(decimated_rec.get_num_samples(i) for i in range(num_segments))
+    if start_frame is None:
+        start_frame = max(decimated_rec.get_num_samples(i) for i in range(N_segments))
+    if end_frame is None:
+        end_frame = max(decimated_rec.get_num_samples(i) for i in range(N_segments))
 
-            for i in range(num_segments):
-                assert decimated_rec.get_num_samples(i) == decimated_parent_traces[i].shape[0]
-                assert np.all(
-                    decimated_rec.get_traces(i, start_frame, end_frame)
-                    == decimated_parent_traces[i][start_frame:end_frame]
-                )
-
-    for i in range(num_segments):
+    for i in range(N_segments):
         assert decimated_rec.get_num_samples(i) == decimated_parent_traces[i].shape[0]
         assert np.all(
             decimated_rec.get_traces(i, start_frame, end_frame) == decimated_parent_traces[i][start_frame:end_frame]
         )
 
 
-def test_decimate_with_times():
-    rec = generate_recording(durations=[5, 10])
-
-    # test with times
-    times = [rec.get_times(0) + 10, rec.get_times(1) + 20]
-    for i, t in enumerate(times):
-        rec.set_times(t, i)
-
-    decimation_factor = 2
-    decimation_offset = 1
-    decimated_rec = DecimateRecording(rec, decimation_factor, decimation_offset=decimation_offset)
-
-    for segment_index in range(rec.get_num_segments()):
-        assert np.allclose(
-            decimated_rec.get_times(segment_index),
-            rec.get_times(segment_index)[decimation_offset::decimation_factor],
-        )
-
-    # test with t_start
-    rec = generate_recording(durations=[5, 10])
-    t_starts = [10, 20]
-    for t_start, rec_segment in zip(t_starts, rec._recording_segments):
-        rec_segment.t_start = t_start
-    decimated_rec = DecimateRecording(rec, decimation_factor, decimation_offset=decimation_offset)
-    for segment_index in range(rec.get_num_segments()):
-        assert np.allclose(
-            decimated_rec.get_times(segment_index),
-            rec.get_times(segment_index)[decimation_offset::decimation_factor],
-        )
-
-
 if __name__ == "__main__":
     test_decimate()
diff --git a/src/spikeinterface/preprocessing/tests/test_scaling.py b/src/spikeinterface/preprocessing/tests/test_scaling.py
index e32d96901..321d7c9df 100644
--- a/src/spikeinterface/preprocessing/tests/test_scaling.py
+++ b/src/spikeinterface/preprocessing/tests/test_scaling.py
@@ -55,11 +55,11 @@ def test_scaling_in_preprocessing_chain():
     recording.set_channel_gains(gains)
     recording.set_channel_offsets(offsets)
 
-    centered_recording = CenterRecording(scale_to_uV(recording=recording), seed=2205)
+    centered_recording = CenterRecording(scale_to_uV(recording=recording))
     traces_scaled_with_argument = centered_recording.get_traces(return_scaled=True)
 
     # Chain preprocessors
-    centered_recording_scaled = CenterRecording(scale_to_uV(recording=recording), seed=2205)
+    centered_recording_scaled = CenterRecording(scale_to_uV(recording=recording))
     traces_scaled_with_preprocessor = centered_recording_scaled.get_traces()
 
     np.testing.assert_allclose(traces_scaled_with_argument, traces_scaled_with_preprocessor)
@@ -68,8 +68,3 @@ def test_scaling_in_preprocessing_chain():
     traces_scaled_with_preprocessor_and_argument = centered_recording_scaled.get_traces(return_scaled=True)
 
     np.testing.assert_allclose(traces_scaled_with_preprocessor, traces_scaled_with_preprocessor_and_argument)
-
-
-if __name__ == "__main__":
-    test_scale_to_uV()
-    test_scaling_in_preprocessing_chain()
diff --git a/src/spikeinterface/preprocessing/tests/test_silence.py b/src/spikeinterface/preprocessing/tests/test_silence.py
index 20d4f6dfc..6c2e8ec8b 100644
--- a/src/spikeinterface/preprocessing/tests/test_silence.py
+++ b/src/spikeinterface/preprocessing/tests/test_silence.py
@@ -9,8 +9,6 @@ from spikeinterface.core import get_noise_levels
 
 import numpy as np
 
-from pathlib import Path
-
 
 def test_silence(create_cache_folder):
 
@@ -48,5 +46,4 @@ def test_silence(create_cache_folder):
 
 
 if __name__ == "__main__":
-    cache_folder = Path(__file__).resolve().parents[4] / "cache_folder"
-    test_silence(cache_folder)
+    test_silence()
diff --git a/src/spikeinterface/preprocessing/tests/test_whiten.py b/src/spikeinterface/preprocessing/tests/test_whiten.py
index b40627d83..04b731de4 100644
--- a/src/spikeinterface/preprocessing/tests/test_whiten.py
+++ b/src/spikeinterface/preprocessing/tests/test_whiten.py
@@ -5,15 +5,13 @@ from spikeinterface.core import generate_recording
 
 from spikeinterface.preprocessing import whiten, scale, compute_whitening_matrix
 
-from pathlib import Path
-
 
 def test_whiten(create_cache_folder):
     cache_folder = create_cache_folder
     rec = generate_recording(num_channels=4, seed=2205)
 
     print(rec.get_channel_locations())
-    random_chunk_kwargs = {"seed": 2205}
+    random_chunk_kwargs = {}
     W1, M = compute_whitening_matrix(rec, "global", random_chunk_kwargs, apply_mean=False, radius_um=None)
     # print(W)
     # print(M)
@@ -49,5 +47,4 @@ def test_whiten(create_cache_folder):
 
 
 if __name__ == "__main__":
-    cache_folder = Path(__file__).resolve().parents[4] / "cache_folder"
-    test_whiten(cache_folder)
+    test_whiten()
diff --git a/src/spikeinterface/sortingcomponents/tests/test_template_matching.py b/src/spikeinterface/sortingcomponents/tests/test_template_matching.py
index 7cd899a3b..cbf1d2993 100644
--- a/src/spikeinterface/sortingcomponents/tests/test_template_matching.py
+++ b/src/spikeinterface/sortingcomponents/tests/test_template_matching.py
@@ -9,8 +9,8 @@ from spikeinterface.sortingcomponents.matching import find_spikes_from_templates
 from spikeinterface.sortingcomponents.tests.common import make_dataset
 
 
-# job_kwargs = dict(n_jobs=-1, chunk_duration="500ms", progress_bar=True)
-job_kwargs = dict(n_jobs=1, chunk_duration="500ms", progress_bar=True)
+job_kwargs = dict(n_jobs=-1, chunk_duration="500ms", progress_bar=True)
+# job_kwargs = dict(n_jobs=1, chunk_duration="500ms", progress_bar=True)
 
 
 def get_sorting_analyzer():
@@ -45,7 +45,7 @@ def test_find_spikes_from_templates(method, sorting_analyzer):
         "templates": templates,
     }
     method_kwargs = {}
-    if method in ("naive", "tdc-peeler", "circus", "tdc-peeler2"):
+    if method in ("naive", "tdc-peeler", "circus"):
         method_kwargs["noise_levels"] = noise_levels
 
     # method_kwargs["wobble"] = {
@@ -61,28 +61,26 @@ def test_find_spikes_from_templates(method, sorting_analyzer):
 
     # print(info)
 
-    DEBUG = True
+    # DEBUG = True
 
-    if DEBUG:
-        import matplotlib.pyplot as plt
-        import spikeinterface.full as si
+    # if DEBUG:
+    #     import matplotlib.pyplot as plt
+    #     import spikeinterface.full as si
 
-        sorting_analyzer.compute("waveforms")
-        sorting_analyzer.compute("templates")
+    #     sorting_analyzer.compute("waveforms")
+    #     sorting_analyzer.compute("templates")
 
-        gt_sorting = sorting_analyzer.sorting
+    #     gt_sorting = sorting_analyzer.sorting
 
-        sorting = NumpySorting.from_times_labels(
-            spikes["sample_index"], spikes["cluster_index"], recording.sampling_frequency
-        )
+    #     sorting = NumpySorting.from_times_labels(spikes["sample_index"], spikes["cluster_index"], recording.sampling_frequency)
 
-        ##metrics = si.compute_quality_metrics(sorting_analyzer, metric_names=["snr"])
+    #     ##metrics = si.compute_quality_metrics(sorting_analyzer, metric_names=["snr"])
 
-        # fig, ax = plt.subplots()
-        # comp = si.compare_sorter_to_ground_truth(gt_sorting, sorting)
-        # si.plot_agreement_matrix(comp, ax=ax)
-        # ax.set_title(method)
-        # plt.show()
+    #     fig, ax = plt.subplots()
+    #     comp = si.compare_sorter_to_ground_truth(gt_sorting, sorting)
+    #     si.plot_agreement_matrix(comp, ax=ax)
+    #     ax.set_title(method)
+    # plt.show()
 
 
 if __name__ == "__main__":
@@ -90,6 +88,6 @@ if __name__ == "__main__":
     # method = "naive"
     # method = "tdc-peeler"
     # method =  "circus"
-    # method = "circus-omp-svd"
-    method = "wobble"
+    method = "circus-omp-svd"
+    # method = "wobble"
     test_find_spikes_from_templates(method, sorting_analyzer)
diff --git a/src/spikeinterface/sortingcomponents/tests/test_wobble.py b/src/spikeinterface/sortingcomponents/tests/test_wobble.py
index 0d46b790a..d6d1e1e0b 100644
--- a/src/spikeinterface/sortingcomponents/tests/test_wobble.py
+++ b/src/spikeinterface/sortingcomponents/tests/test_wobble.py
@@ -44,7 +44,7 @@ def test_compress_templates():
         elif test_case == "num_channels == num_samples":
             num_channels = rng.integers(1, 100)
             num_samples = num_channels
-        templates = rng.random((num_templates, num_samples, num_channels), dtype=np.float32)
+        templates = rng.random((num_templates, num_samples, num_channels))
         full_rank = np.minimum(num_samples, num_channels)
         approx_rank = rng.integers(1, full_rank)
 
@@ -66,31 +66,15 @@ def test_compress_templates():
         assert np.all(singular_full >= 0)
         # check that svd matrices are orthonormal if applicable
         if num_channels > num_samples:
-            assert np.allclose(
-                np.matmul(temporal_full, temporal_full.transpose(0, 2, 1)),
-                np.eye(num_samples, dtype=np.float32),
-                atol=1e-3,
-            )
+            assert np.allclose(np.matmul(temporal_full, temporal_full.transpose(0, 2, 1)), np.eye(num_samples))
         elif num_samples > num_channels:
-            assert np.allclose(
-                np.matmul(spatial_full, spatial_full.transpose(0, 2, 1)),
-                np.eye(num_channels, dtype=np.float32),
-                atol=1e-3,
-            )
+            assert np.allclose(np.matmul(spatial_full, spatial_full.transpose(0, 2, 1)), np.eye(num_channels))
         elif num_channels == num_samples:
-            assert np.allclose(
-                np.matmul(temporal_full, temporal_full.transpose(0, 2, 1)),
-                np.eye(num_samples, dtype=np.float32),
-                atol=1e-3,
-            )
-            assert np.allclose(
-                np.matmul(spatial_full, spatial_full.transpose(0, 2, 1)),
-                np.eye(num_channels, dtype=np.float32),
-                atol=1e-3,
-            )
+            assert np.allclose(np.matmul(temporal_full, temporal_full.transpose(0, 2, 1)), np.eye(num_samples))
+            assert np.allclose(np.matmul(spatial_full, spatial_full.transpose(0, 2, 1)), np.eye(num_channels))
             # check that the full rank svd matrices reconstruct the original templates
             reconstructed_templates = np.matmul(temporal_full * singular_full[:, np.newaxis, :], spatial_full)
-            assert np.allclose(reconstructed_templates, templates, atol=1e-3)
+            assert np.allclose(reconstructed_templates, templates)
 
 
 def test_upsample_and_jitter():
@@ -227,33 +211,18 @@ def test_compute_objective():
     approx_rank = rng.integers(1, num_samples)
     num_channels = rng.integers(1, 100)
     chunk_len = rng.integers(num_samples * 2, num_samples * 10)
-    traces = rng.random((chunk_len, num_channels), dtype=np.float32)
+    traces = rng.random((chunk_len, num_channels))
     temporal = rng.random((num_templates, num_samples, approx_rank))
     singular = rng.random((num_templates, approx_rank))
     spatial = rng.random((num_templates, approx_rank, num_channels))
-
-    spatial_transformed = np.moveaxis(spatial, [0, 1, 2], [1, 0, 2])
-    temporal_transformed = np.moveaxis(temporal, [0, 1, 2], [1, 2, 0])
-    singular_transformed = singular.T[:, :, np.newaxis]
-
-    compressed_templates_transformed = (
-        temporal_transformed,
-        singular_transformed,
-        spatial_transformed,
-        temporal_transformed,
-    )
-    norm_squared = np.random.rand(num_templates)
-
-    template_data_transformed = wobble.TemplateData(
-        compressed_templates=compressed_templates_transformed, pairwise_convolution=[], norm_squared=norm_squared
-    )
-    # Act: run compute_objective
-    objective = wobble.compute_objective(traces, template_data_transformed, approx_rank, engine="numpy")
-
     compressed_templates = (temporal, singular, spatial, temporal)
+    norm_squared = np.random.rand(num_templates)
     template_data = wobble.TemplateData(
         compressed_templates=compressed_templates, pairwise_convolution=[], norm_squared=norm_squared
     )
+
+    # Act: run compute_objective
+    objective = wobble.compute_objective(traces, template_data, approx_rank)
     expected_objective = compute_objective_loopy(traces, template_data, approx_rank)
 
     # Assert: check shape and equivalence to expected_objective