Repository: spikeinterface/spikeinterface
Issue: Use spike_amplitudes extension for amplitude_cutoff metric

The amplitude_cutoff metric relies on spike amplitudes to estimate the pdf of the amplitude distribution. Currently the implementation uses the amplitudes of the sampled waveforms (by default 500 samples), which is too few for a proper estimation.

A proposed solution is to use the `spike_amplitudes` extension, which contains computed amplitudes for all spikes in the recording. If that is not available, the metric should not be calculated and a warning should be issued.

@TomBugnon  @samuelgarcia 

Changed files:
spikeinterface/core/__init__.py, spikeinterface/curation/tests/test_sortingview_curation.py, spikeinterface/qualitymetrics/tests/test_metrics_functions.py, spikeinterface/core/segmentutils.py, spikeinterface/curation/tests/test_auto_merge.py, spikeinterface/qualitymetrics/tests/test_quality_metric_calculator.py, spikeinterface/qualitymetrics/misc_metrics.py, spikeinterface/core/generate.py, spikeinterface/curation/tests/test_remove_redundant.py, spikeinterface/widgets/tests/test_widgets.py, spikeinterface/curation/auto_merge.py

diff --git a/spikeinterface/core/__init__.py b/spikeinterface/core/__init__.py
index 9906c518..775c2df5 100644
--- a/spikeinterface/core/__init__.py
+++ b/spikeinterface/core/__init__.py
@@ -35,7 +35,7 @@ from .unitsaggregationsorting import UnitsAggregationSorting, aggregate_units
 # generator of simple object for testing or examples
 from .generate import (generate_recording, generate_sorting,
   create_sorting_npz, generate_snippets,
-  synthesize_random_firings,  inject_some_duplicat_units,
+  synthesize_random_firings,  inject_some_duplicate_units,
   inject_some_split_units, synthetize_spike_train_bad_isi)
 
 # utils to append and concatenate segment (equivalent to OLD MultiRecordingTimeExtractor)
@@ -51,6 +51,8 @@ from .segmentutils import (
     AppendSegmentSorting,
     split_sorting,
     SplitSegmentSorting,
+    select_segment_sorting,
+    SelectSegmentSorting
 )
 
 # default folder
diff --git a/spikeinterface/core/generate.py b/spikeinterface/core/generate.py
index cade3695..c8aad413 100644
--- a/spikeinterface/core/generate.py
+++ b/spikeinterface/core/generate.py
@@ -223,7 +223,7 @@ def clean_refractory_period(times, refractory_period):
 
     return times
 
-def inject_some_duplicat_units(sorting, num=4, max_shift=5, ratio=None, seed=None):
+def inject_some_duplicate_units(sorting, num=4, max_shift=5, ratio=None, seed=None):
     """
     Inject some duplicate units in a sorting.
     The peak shift can be control in a range.
@@ -248,6 +248,7 @@ def inject_some_duplicat_units(sorting, num=4, max_shift=5, ratio=None, seed=Non
     """
     other_ids = np.arange(np.max(sorting.unit_ids) +1 , np.max(sorting.unit_ids) + num + 1)
     shifts = np.random.RandomState(seed).randint(low=-max_shift, high=max_shift, size=num)
+    shifts[shifts == 0] += max_shift
     unit_peak_shifts = dict(zip(other_ids, shifts))
 
     spiketrains = []
diff --git a/spikeinterface/core/segmentutils.py b/spikeinterface/core/segmentutils.py
index f95aefa9..9db7eee3 100644
--- a/spikeinterface/core/segmentutils.py
+++ b/spikeinterface/core/segmentutils.py
@@ -268,7 +268,8 @@ def split_recording(recording: BaseRecording):
     return recording_list
 
 
-select_segment_recording = define_function_from_class(source_class=SelectSegmentRecording, name='select_segment_recording')
+select_segment_recording = define_function_from_class(source_class=SelectSegmentRecording,
+                                                      name='select_segment_recording')
 
 
 class AppendSegmentSorting(BaseSorting):
@@ -322,6 +323,7 @@ class ProxyAppendSortingSegment(BaseSortingSegment):
 
 append_sortings = define_function_from_class(source_class=AppendSegmentSorting, name='append_sortings')
 
+
 class SplitSegmentSorting(BaseSorting):
     """Splits a sorting with a single segment to multiple segments
     based on the given list of recordings (must be in order)
@@ -373,3 +375,39 @@ class SplitSegmentSorting(BaseSorting):
                         'recording_list': [recording.to_dict() for recording in recording_list]}
 
 split_sorting = define_function_from_class(source_class=SplitSegmentSorting, name='split_sorting')
+
+
+class SelectSegmentSorting(BaseSorting):
+    """
+    Return a new sorting with a single segment from a multi-segment sorting.
+
+    Parameters
+    ----------
+    sorting : BaseSorting
+        The multi-segment sorting
+    segment_indices : list of int
+        The segment indices to select
+    """
+
+    def __init__(self, sorting: BaseSorting, segment_indices: Union[int, List[int]]):
+        BaseSorting.__init__(self, sorting.get_sampling_frequency(), 
+                             sorting.unit_ids)
+        sorting.copy_metadata(self)
+        
+        if isinstance(segment_indices, int):
+            segment_indices = [segment_indices]
+        
+        num_segments = sorting.get_num_segments()
+        assert all(0 <= s < num_segments for s in segment_indices), \\
+            f"'segment_index' must be between 0 and {num_segments - 1}"
+
+        for segment_index in segment_indices:
+            sort_seg = sorting._sorting_segments[segment_index]
+            self.add_sorting_segment(sort_seg)
+
+        self._kwargs = {'sorting': sorting.to_dict(),
+                        'segment_indices': [int(s) for s in segment_indices]}
+
+
+select_segment_sorting = define_function_from_class(source_class=SelectSegmentSorting,
+                                                    name='select_segment_sorting')
\\ No newline at end of file
diff --git a/spikeinterface/curation/auto_merge.py b/spikeinterface/curation/auto_merge.py
index ceae3a0f..141897f7 100644
--- a/spikeinterface/curation/auto_merge.py
+++ b/spikeinterface/curation/auto_merge.py
@@ -135,7 +135,7 @@ def get_potential_auto_merge(
 
     # STEP 2 : remove contaminated auto corr
     if 'remove_contaminated' in steps:
-        nb_violations, contaminations = compute_refrac_period_violations(we, refractory_period_ms=refractory_period_ms,
+        contaminations, nb_violations = compute_refrac_period_violations(we, 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()))
@@ -434,7 +434,7 @@ def check_improve_contaminations_score(we, pair_mask, contaminations,
         # make a lazy fake WaveformExtractor to compute contamination and firing rate
         we_new = MockWaveformExtractor(recording, sorting_merged)
 
-        _, new_contaminations = compute_refrac_period_violations(we_new, refractory_period_ms=refractory_period_ms,
+        new_contaminations, _ = compute_refrac_period_violations(we_new, refractory_period_ms=refractory_period_ms,
                                     censored_period_ms=censored_period_ms)
         c_new = new_contaminations[unit_id1]
         f_new = compute_firing_rate(we_new)[unit_id1]
diff --git a/spikeinterface/curation/tests/test_auto_merge.py b/spikeinterface/curation/tests/test_auto_merge.py
index 15b843d9..e8c310d8 100644
--- a/spikeinterface/curation/tests/test_auto_merge.py
+++ b/spikeinterface/curation/tests/test_auto_merge.py
@@ -22,7 +22,7 @@ set_global_tmp_folder(cache_folder)
 def test_get_auto_merge_list():
 
 
-    rec, sorting = toy_example(num_segments=1, num_units=5, duration=[300.], firing_rate=20.)
+    rec, sorting = toy_example(num_segments=1, num_units=5, duration=[300.], firing_rate=20., seed=0)
     
     num_unit_splited = 1
     num_split = 2
diff --git a/spikeinterface/curation/tests/test_remove_redundant.py b/spikeinterface/curation/tests/test_remove_redundant.py
index add7bdd3..39dc17c9 100644
--- a/spikeinterface/curation/tests/test_remove_redundant.py
+++ b/spikeinterface/curation/tests/test_remove_redundant.py
@@ -7,7 +7,7 @@ import pytest
 import numpy as np
 
 from spikeinterface import WaveformExtractor, load_extractor, extract_waveforms, NumpySorting, set_global_tmp_folder
-from spikeinterface.core.generate import inject_some_duplicat_units
+from spikeinterface.core.generate import inject_some_duplicate_units
 
 from spikeinterface.extractors import toy_example
 
@@ -22,9 +22,9 @@ else:
 set_global_tmp_folder(cache_folder)
 
 def test_remove_redundant_units():
-    rec, sorting = toy_example(num_segments=1, duration=[10.])
+    rec, sorting = toy_example(num_segments=1, duration=[10.], seed=0)
     
-    sorting_with_dup = inject_some_duplicat_units(sorting, ratio=0.8, num=4)
+    sorting_with_dup = inject_some_duplicate_units(sorting, ratio=0.8, num=4, seed=1)
     
     rec = rec.save()
     sorting_with_dup = sorting_with_dup.save()
@@ -33,8 +33,8 @@ def test_remove_redundant_units():
         shutil.rmtree(wf_folder)
     we = extract_waveforms(rec, sorting_with_dup, folder=wf_folder)
     print(we)
-    
-    for remove_strategy in ('minimum_shift', 'highest_amplitude'):
+
+    for remove_strategy in ('max_spikes', 'minimum_shift', 'highest_amplitude'):
         sorting_clean = remove_redundant_units(we, remove_strategy=remove_strategy)
         # print(sorting_clean)
         # print(sorting_clean.unit_ids)
diff --git a/spikeinterface/curation/tests/test_sortingview_curation.py b/spikeinterface/curation/tests/test_sortingview_curation.py
index 87628b7b..0b63511a 100644
--- a/spikeinterface/curation/tests/test_sortingview_curation.py
+++ b/spikeinterface/curation/tests/test_sortingview_curation.py
@@ -47,7 +47,7 @@ def test_sortingview_curation():
     jot_uri = "jot://cENcosqWTArO"
     # curation link: 
     # https://figurl.org/f?v=gs://figurl/spikesortingview-10&d=sha1://1ba03f81e62ec7cb2e3e46898830f92cdf5e026f&label=SpikeInterface%20-%20Sorting%20Summary&s={%22sortingCuration%22:%22jot://cENcosqWTArO%22}
-    sorting_curated_jot = apply_sortingview_curation(sorting, uri=jot_uri, verbose=True)
+    sorting_curated_jot = apply_sortingview_curation(sorting, uri_or_json=jot_uri, verbose=True)
 
     assert len(sorting_curated_jot.unit_ids) == 7
     assert "#0-#1" in sorting_curated_jot.unit_ids
@@ -56,9 +56,10 @@ def test_sortingview_curation():
     assert "noise" in sorting_curated_jot.get_property_keys()
     assert "reject" in sorting_curated_jot.get_property_keys()
 
-    sorting_curated_jot_accepted = apply_sortingview_curation(sorting, uri=jot_uri, include_labels=["accept"])
-    sorting_curated_jot_rejected = apply_sortingview_curation(sorting, uri=jot_uri, exclude_labels=["reject"])
-    sorting_curated_jot_rejected1 = apply_sortingview_curation(sorting, uri=jot_uri, exclude_labels=["noise", "reject"])
+    sorting_curated_jot_accepted = apply_sortingview_curation(sorting, uri_or_json=jot_uri, include_labels=["accept"])
+    sorting_curated_jot_rejected = apply_sortingview_curation(sorting, uri_or_json=jot_uri, exclude_labels=["reject"])
+    sorting_curated_jot_rejected1 = apply_sortingview_curation(sorting, uri_or_json=jot_uri,
+                                                               exclude_labels=["noise", "reject"])
     assert len(sorting_curated_jot_accepted.unit_ids) == 4
     assert len(sorting_curated_jot_rejected.unit_ids) == 5
     assert len(sorting_curated_jot_rejected1.unit_ids) == 4
@@ -66,7 +67,7 @@ def test_sortingview_curation():
     # curation_link: 
     # https://figurl.org/f?v=gs://figurl/spikesortingview-10&d=sha1://1ba03f81e62ec7cb2e3e46898830f92cdf5e026f&label=SpikeInterface%20-%20Sorting%20Summary&s={%22sortingCuration%22:%22sha1://59feb326204cf61356f1a2eb31f04d8e0177c4f1%22}
     sha_uri = "sha1://59feb326204cf61356f1a2eb31f04d8e0177c4f1"
-    sorting_curated_sha = apply_sortingview_curation(sorting, uri=sha_uri, verbose=True)
+    sorting_curated_sha = apply_sortingview_curation(sorting, uri_or_json=sha_uri, verbose=True)
 
     assert len(sorting_curated_sha.unit_ids) == 9
     assert "#8-#9" in sorting_curated_sha.unit_ids
@@ -74,9 +75,10 @@ def test_sortingview_curation():
     assert "mua" in sorting_curated_sha.get_property_keys()
     assert "artifact" in sorting_curated_sha.get_property_keys()
 
-    sorting_curated_sha_accepted = apply_sortingview_curation(sorting, uri=sha_uri, include_labels=["accept"])
-    sorting_curated_sha_mua = apply_sortingview_curation(sorting, uri=sha_uri, exclude_labels=["mua"])
-    sorting_curated_sha_mua1 = apply_sortingview_curation(sorting, uri=sha_uri, exclude_labels=["artifact", "mua"])
+    sorting_curated_sha_accepted = apply_sortingview_curation(sorting, uri_or_json=sha_uri, include_labels=["accept"])
+    sorting_curated_sha_mua = apply_sortingview_curation(sorting, uri_or_json=sha_uri, exclude_labels=["mua"])
+    sorting_curated_sha_mua1 = apply_sortingview_curation(sorting, uri_or_json=sha_uri,
+                                                          exclude_labels=["artifact", "mua"])
     assert len(sorting_curated_sha_accepted.unit_ids) == 3
     assert len(sorting_curated_sha_mua.unit_ids) == 6
     assert len(sorting_curated_sha_mua1.unit_ids) == 5
diff --git a/spikeinterface/qualitymetrics/misc_metrics.py b/spikeinterface/qualitymetrics/misc_metrics.py
index bf28a973..c7137c52 100644
--- a/spikeinterface/qualitymetrics/misc_metrics.py
+++ b/spikeinterface/qualitymetrics/misc_metrics.py
@@ -95,15 +95,16 @@ def compute_firing_rate(waveform_extractor):
 _default_params["firing_rate"] = dict()
 
 
-def compute_presence_ratio(waveform_extractor, num_bin_edges=101):
+def compute_presence_ratio(waveform_extractor, bin_duration_s=60):
     """Calculate the presence ratio, representing the fraction of time the unit is firing.
 
     Parameters
     ----------
     waveform_extractor : WaveformExtractor
         The waveform extractor object.
-    num_bin_edges : int, optional, default: 101
-        The number of bins edges to use to compute the presence ratio.
+    bin_duration_s : float, optional, default: 60
+        The duration of each bin in seconds. If the duration is less than this value, 
+        presence_ratio is set to NaN
 
     Returns
     -------
@@ -123,23 +124,31 @@ def compute_presence_ratio(waveform_extractor, num_bin_edges=101):
 
     seg_length = [recording.get_num_samples(i) for i in range(num_segs)]
     total_length = np.sum(seg_length)
-
-    presence_ratio = {}
-    for unit_id in unit_ids:
-        spiketrain = []
-        for segment_index in range(num_segs):
-            st = sorting.get_unit_spike_train(unit_id=unit_id, segment_index=segment_index)
-            st = st + np.sum(seg_length[:segment_index])
-            spiketrain.append(st)
-        spiketrain = np.concatenate(spiketrain)
-        h, b = np.histogram(spiketrain, np.linspace(0, total_length, num_bin_edges))
-        presence_ratio[unit_id] = np.sum(h > 0) / (num_bin_edges - 1)
+    bin_duration_samples = int((bin_duration_s * recording.sampling_frequency))
+    num_bin_edges = total_length // bin_duration_samples + 1
+    bins = np.arange(num_bin_edges) * bin_duration_samples
+
+    if total_length < bin_duration_samples:
+        warnings.warn(f"Bin duration of {bin_duration_s}s is larger than recording duration. "
+                      f"Presence ratios are set to NaN.")
+        presence_ratio = {unit_id: np.nan for unit_id in sorting.unit_ids}
+    else:
+        presence_ratio = {}
+        for unit_id in unit_ids:
+            spiketrain = []
+            for segment_index in range(num_segs):
+                st = sorting.get_unit_spike_train(unit_id=unit_id, segment_index=segment_index)
+                st = st + np.sum(seg_length[:segment_index])
+                spiketrain.append(st)
+            spiketrain = np.concatenate(spiketrain)
+            h, _ = np.histogram(spiketrain, bins=bins)
+            presence_ratio[unit_id] = np.sum(h > 0) / (num_bin_edges - 1)
 
     return presence_ratio
 
 
 _default_params["presence_ratio"] = dict(
-    num_bin_edges=101
+    bin_duration_s=60
 )
 
 
@@ -207,7 +216,6 @@ def compute_isi_violations(waveform_extractor, isi_threshold_ms=1.5, min_isi_ms=
     It computes several metrics related to isi violations:
         * isi_violations_ratio: the relative firing rate of the hypothetical neurons that are
                                 generating the ISI violations. Described in [1]. See Notes.
-        * isi_violation_rate: number of ISI violations divided by total rate
         * isi_violation_count: number of ISI violations
 
     Parameters
@@ -372,7 +380,8 @@ _default_params["rp_violations"] = dict(
 
 
 def compute_amplitudes_cutoff(waveform_extractor, peak_sign='neg',
-                              num_histogram_bins=500, histogram_smoothing_value=3):
+                              num_histogram_bins=100, histogram_smoothing_value=3,
+                              amplitudes_bins_min_ratio=5):
     """Calculate approximate fraction of spikes missing from a distribution of amplitudes.
 
     Parameters
@@ -381,10 +390,14 @@ def compute_amplitudes_cutoff(waveform_extractor, peak_sign='neg',
         The waveform extractor object.
     peak_sign : {'neg', 'pos', 'both'}
         The sign of the template to compute best channels.
-    num_histogram_bins : int, optional, default: 500
+    num_histogram_bins : int, optional, default: 100
         The number of bins to use to compute the amplitude histogram.
     histogram_smoothing_value : int, optional, default: 3
         Controls the smoothing applied to the amplitude histogram.
+    amplitudes_bins_min_ratio : int, optional, default: 5
+        The minimum ratio between number of amplitudes for a unit and the number of bins.
+        If the ratio is less than this threshold, the amplitude_cutoff for the unit is set 
+        to NaN
 
     Returns
     -------
@@ -399,14 +412,11 @@ def compute_amplitudes_cutoff(waveform_extractor, peak_sign='neg',
 
     Notes
     -----
-    This approach assumes the amplitude histogram is symmetric (not valid in the presence of drift)
-
-    Important note: here the amplitudes are extracted from the waveform extractor.
-    It means that the number of spike to estimate amplitude is low
-    See: WaveformExtractor.set_params(max_spikes_per_unit=500)
-
+    This approach assumes the amplitude histogram is symmetric (not valid in the presence of drift).
+    If available, amplitudes are extracted from the "spike_amplitude" extension (recommended). 
+    If the "spike_amplitude" extension is not available, the amplitudes are extracted from the waveform extractor,
+    which usually has waveforms for a small subset of spikes (500 by default).
     """
-
     recording = waveform_extractor.recording
     sorting = waveform_extractor.sorting
     unit_ids = sorting.unit_ids
@@ -415,17 +425,36 @@ def compute_amplitudes_cutoff(waveform_extractor, peak_sign='neg',
 
     extremum_channels_ids = get_template_extremum_channel(waveform_extractor, peak_sign=peak_sign)
 
+    spike_amplitudes = None
+    invert_amplitudes = False
+    if waveform_extractor.is_extension("spike_amplitudes"):
+        amp_calculator = waveform_extractor.load_extension("spike_amplitudes")
+        spike_amplitudes = amp_calculator.get_data(outputs="by_unit")
+        if amp_calculator._params["peak_sign"] == "pos":
+            invert_amplitudes = True
+    else:
+        if peak_sign == "pos":
+            invert_amplitudes = True
+
     all_fraction_missing = {}
+    nan_units = []
     for unit_id in unit_ids:
-        waveforms = waveform_extractor.get_waveforms(unit_id)
-        chan_id = extremum_channels_ids[unit_id]
-        chan_ind = recording.id_to_index(chan_id)
-        amplitudes = waveforms[:, before, chan_ind]
+        if spike_amplitudes is None:
+            waveforms = waveform_extractor.get_waveforms(unit_id)
+            chan_id = extremum_channels_ids[unit_id]
+            chan_ind = recording.id_to_index(chan_id)
+            amplitudes = waveforms[:, before, chan_ind]
+        else:
+            amplitudes = np.concatenate([spike_amps[unit_id] for spike_amps in spike_amplitudes])
+
+        if len(amplitudes) / num_histogram_bins < amplitudes_bins_min_ratio:
+            nan_units.append(unit_id)
+            all_fraction_missing[unit_id] = np.nan
+            continue
 
         # change amplitudes signs in case peak_sign is pos
-        if peak_sign == "pos":
+        if invert_amplitudes:
             amplitudes = -amplitudes
-
         h, b = np.histogram(amplitudes, num_histogram_bins, density=True)
 
         # TODO : change with something better scipy.ndimage.gaussian_filter1d
@@ -445,12 +474,16 @@ def compute_amplitudes_cutoff(waveform_extractor, peak_sign='neg',
         fraction_missing = np.min([fraction_missing, 0.5])
         all_fraction_missing[unit_id] = fraction_missing
 
+    if len(nan_units) > 0:
+        warnings.warn(f"Units {nan_units} have too few spikes and "
+                       "amplitude_cutoff is set to NaN")
+
     return all_fraction_missing
 
 
 _default_params["amplitude_cutoff"] = dict(
     peak_sign='neg',
-    num_histogram_bins=500,
+    num_histogram_bins=100,
     histogram_smoothing_value=3
 )
 
diff --git a/spikeinterface/qualitymetrics/tests/test_metrics_functions.py b/spikeinterface/qualitymetrics/tests/test_metrics_functions.py
index 8f70c00d..e1032b3d 100644
--- a/spikeinterface/qualitymetrics/tests/test_metrics_functions.py
+++ b/spikeinterface/qualitymetrics/tests/test_metrics_functions.py
@@ -20,6 +20,27 @@ else:
     cache_folder = Path("cache_folder") / "qualitymetrics"
 
 
+def _simulated_data():
+    max_time = 100.0
+
+    trains = [synthetize_spike_train_bad_isi(max_time, 10, 2),
+              synthetize_spike_train_bad_isi(max_time, 5, 4),
+              synthetize_spike_train_bad_isi(max_time, 5, 10)]
+
+    labels = [np.ones((len(trains[i]),), dtype='int') * i for i in range(len(trains))]
+
+    spike_times = np.concatenate(trains)
+    spike_clusters = np.concatenate(labels)
+
+    order = np.argsort(spike_times)
+
+    indexes = np.arange(0, max_time + 1, 1 / 30000)
+    spike_times = np.searchsorted(indexes, spike_times[order], side="left")
+    spike_clusters = spike_clusters[order]
+
+    return {"duration": max_time, "times": spike_times, "labels": spike_clusters }
+
+
 def setup_module():
     for folder_name in ('toy_rec', 'toy_sorting', 'toy_waveforms'):
         if (cache_folder / folder_name).is_dir():
@@ -84,24 +105,7 @@ def test_nearest_neighbors_metrics():
 
 @pytest.fixture
 def simulated_data():
-    max_time = 100.0
-
-    trains = [synthetize_spike_train_bad_isi(max_time, 10, 2),
-              synthetize_spike_train_bad_isi(max_time, 5, 4),
-              synthetize_spike_train_bad_isi(max_time, 5, 10)]
-
-    labels = [np.ones((len(trains[i]),), dtype='int') * i for i in range(len(trains))]
-
-    spike_times = np.concatenate(trains)
-    spike_clusters = np.concatenate(labels)
-
-    order = np.argsort(spike_times)
-
-    indexes = np.arange(0, max_time + 1, 1 / 30000)
-    spike_times = np.searchsorted(indexes, spike_times[order], side="left")
-    spike_clusters = spike_clusters[order]
-
-    return {"duration": max_time, "times": spike_times, "labels": spike_clusters }
+    return _simulated_data()
 
 
 def setup_dataset(spike_data, score_detection=1):
@@ -120,7 +124,6 @@ def setup_dataset(spike_data, score_detection=1):
 
 
 def test_calculate_firing_rate_num_spikes(simulated_data):
-
     we = setup_dataset(simulated_data)
     firing_rates = compute_firing_rate(we)
     num_spikes = compute_num_spikes(we)
@@ -131,26 +134,23 @@ def test_calculate_firing_rate_num_spikes(simulated_data):
 
 def test_calculate_amplitude_cutoff(simulated_data):
     we = setup_dataset(simulated_data, score_detection=0.5)
-    amp_cuts = compute_amplitudes_cutoff(we)
-    assert amp_cuts == {0: 0.021888823653840105, 1: 0.28298787894332295, 2: 0.28298787894332295}
+    amp_cuts = compute_amplitudes_cutoff(we, num_histogram_bins=10)
+    assert amp_cuts == {0: 0.3307144004373338, 1: 0.43482247296942045, 2: 0.43482247296942045}
 
 
 def test_calculate_snrs(simulated_data):
-
     we = setup_dataset(simulated_data, score_detection=0.5)
     snr = compute_snrs(we)
     assert np.allclose(np.array(list(snr.values())), np.array([12.92, 12.99, 12.99]), atol=0.05)
 
 
 def test_calculate_presence_ratio(simulated_data):
-
     we = setup_dataset(simulated_data)
-    ratios = compute_presence_ratio(we)
+    ratios = compute_presence_ratio(we, bin_duration_s=10)
     assert ratios == {0: 1.0, 1: 1.0, 2: 1.0}
 
 
 def test_calculate_isi_violations(simulated_data):
-
     we = setup_dataset(simulated_data)
     ratio, count = compute_isi_violations(we, 1, 0.0)
 
@@ -158,7 +158,6 @@ def test_calculate_isi_violations(simulated_data):
     assert count == {0: 2, 1: 4, 2: 10}
 
 def test_calculate_rp_violations(simulated_data):
-
     we = setup_dataset(simulated_data)
     rp_contamination, rp_violations = compute_refrac_period_violations(we, 1, 0.0)
 
@@ -168,5 +167,6 @@ def test_calculate_rp_violations(simulated_data):
 
 if __name__ == '__main__':
     setup_module()
-
-    test_calculate_pc_metrics()
+    sim_data = _simulated_data()
+    test_calculate_amplitude_cutoff(sim_data)
+    test_calculate_presence_ratio(sim_data)
diff --git a/spikeinterface/qualitymetrics/tests/test_quality_metric_calculator.py b/spikeinterface/qualitymetrics/tests/test_quality_metric_calculator.py
index 0c658e92..57b5227d 100644
--- a/spikeinterface/qualitymetrics/tests/test_quality_metric_calculator.py
+++ b/spikeinterface/qualitymetrics/tests/test_quality_metric_calculator.py
@@ -1,13 +1,14 @@
 import unittest
 import pytest
+import warnings
 from pathlib import Path
 import numpy as np
 
-from spikeinterface import WaveformExtractor, load_extractor, extract_waveforms
+from spikeinterface import WaveformExtractor, load_extractor, extract_waveforms, split_recording, select_segment_sorting
 from spikeinterface.extractors import toy_example
 from spikeinterface.core import get_template_channel_sparsity
 
-from spikeinterface.postprocessing import WaveformPrincipalComponent
+from spikeinterface.postprocessing import compute_principal_components, compute_spike_amplitudes
 from spikeinterface.preprocessing import scale
 from spikeinterface.qualitymetrics import QualityMetricCalculator, get_default_qm_params
 
@@ -32,7 +33,7 @@ class QualityMetricsExtensionTest(WaveformExtensionCommonTestSuite, unittest.Tes
     def setUp(self):
         super().setUp()
         self.cache_folder = cache_folder
-        recording, sorting = toy_example(num_segments=2, num_units=10, duration=60)
+        recording, sorting = toy_example(num_segments=2, num_units=10, duration=120)
         if (cache_folder / 'toy_rec_long').is_dir():
             recording = load_extractor(self.cache_folder / 'toy_rec_long')
         else:
@@ -43,11 +44,27 @@ class QualityMetricsExtensionTest(WaveformExtensionCommonTestSuite, unittest.Tes
             sorting = sorting.save(folder=self.cache_folder / 'toy_sorting_long')
         we_long = extract_waveforms(recording, sorting,
                                     self.cache_folder / 'toy_waveforms_long',
-                                    max_spikes_per_unit=None,
-                                    overwrite=True)
+                                    max_spikes_per_unit=500,
+                                    overwrite=True,
+                                    seed=0)
+        # make a short we for testing amp cutoff
+        recording_one = split_recording(recording)[0]
+        sorting_one = select_segment_sorting(sorting, [0])
+
+        nsec_short = 30
+        recording_short = recording_one.frame_slice(start_frame=0,
+                                                    end_frame=int(nsec_short * recording.sampling_frequency))
+        sorting_short = sorting_one.frame_slice(start_frame=0,
+                                                end_frame=int(nsec_short * recording.sampling_frequency))
+        we_short = extract_waveforms(recording_short, sorting_short,
+                                     self.cache_folder / 'toy_waveforms_short',
+                                     max_spikes_per_unit=500,
+                                     overwrite=True,
+                                     seed=0)
         self.sparsity_long = get_template_channel_sparsity(we_long, method="radius",
                                                            radius_um=50)
         self.we_long = we_long
+        self.we_short = we_short
 
     def test_metrics(self):
         we = self.we_long
@@ -66,9 +83,7 @@ class QualityMetricsExtensionTest(WaveformExtensionCommonTestSuite, unittest.Tes
         # print(metrics)
 
         # with PCs
-        pca = WaveformPrincipalComponent(we)
-        pca.set_params(n_components=5, mode='by_channel_local')
-        pca.run()
+        _ = compute_principal_components(we, n_components=5, mode='by_channel_local')
         metrics = self.extension_class.get_extension_function()(we, seed=0)
         assert 'isolation_distance' in metrics.columns
 
@@ -89,6 +104,44 @@ class QualityMetricsExtensionTest(WaveformExtensionCommonTestSuite, unittest.Tes
         #     assert np.allclose(metrics[metric_name], metrics_par[metric_name])
         # print(metrics_sparse)
 
+    def test_amplitude_cutoff(self):
+        we = self.we_short
+        _ = compute_spike_amplitudes(we, peak_sign="neg")
+
+        # If too few spikes, should raise a warning and set amplitude cutoffs to nans
+        with pytest.warns(UserWarning) as w:
+            metrics = self.extension_class.get_extension_function()(
+                we, metric_names=['amplitude_cutoff'], peak_sign="neg")
+        assert all(np.isnan(cutoff) for cutoff in metrics["amplitude_cutoff"].values)
+
+        # now we decrease the number of bins and check that amplitude cutoffs are correctly computed
+        qm_params=dict(amplitude_cutoff=dict(num_histogram_bins=5))
+        with warnings.catch_warnings():
+            warnings.simplefilter("error")
+            metrics = self.extension_class.get_extension_function()(
+                we, metric_names=['amplitude_cutoff'], peak_sign="neg", qm_params=qm_params)
+        assert all(not np.isnan(cutoff) for cutoff in metrics["amplitude_cutoff"].values)
+
+
+    def test_presence_ratio(self):
+        we = self.we_long
+
+        total_duration = we.recording.get_total_duration()
+        # If bin_duration_s is larger than total duration, should raise a warning and set presence ratios to nans
+        qm_params=dict(presence_ratio=dict(bin_duration_s=total_duration+1))
+        with pytest.warns(UserWarning) as w:
+            metrics = self.extension_class.get_extension_function()(
+                we, metric_names=['presence_ratio'], qm_params=qm_params)
+        assert all(np.isnan(ratio) for ratio in metrics["presence_ratio"].values)
+
+        # now we decrease the bin_duration_s and check that presenc ratios are correctly computed
+        qm_params=dict(presence_ratio=dict(bin_duration_s=total_duration // 10))
+        with warnings.catch_warnings():
+            warnings.simplefilter("error")
+            metrics = self.extension_class.get_extension_function()(
+                we, metric_names=['presence_ratio'], qm_params=qm_params)
+        assert all(not np.isnan(ratio) for ratio in metrics["presence_ratio"].values)
+
     def test_peak_sign(self):
         we = self.we_long
         rec = we.recording
@@ -102,26 +155,29 @@ class QualityMetricsExtensionTest(WaveformExtensionCommonTestSuite, unittest.Tes
         we_inv.set_params(ms_before=3., ms_after=4., max_spikes_per_unit=None)
         we_inv.run_extract_waveforms(n_jobs=1, chunk_size=30000)
 
-        # neg_qm_params = get_default_qm_params()
-        # neg_qm_params["snr"]["peak_sign"] = "neg"
-        # neg_qm_params["amplitude_cutoff"]["peak_sign"] = "neg"
-        # pos_qm_params = get_default_qm_params()
-        # pos_qm_params["snr"]["peak_sign"] = "pos"
-        # pos_qm_params["amplitude_cutoff"]["peak_sign"] = "pos"
+        # compute amplitudes
+        _ = compute_spike_amplitudes(we, peak_sign="neg")
+        _ = compute_spike_amplitudes(we_inv, peak_sign="pos")
+
 
         # without PC
         metrics = self.extension_class.get_extension_function()(
             we, metric_names=['snr', 'amplitude_cutoff'], peak_sign="neg")
         metrics_inv = self.extension_class.get_extension_function()(
             we_inv, metric_names=['snr', 'amplitude_cutoff'], peak_sign="pos")
+        # print(metrics)
+        # print(metrics_inv)
+        # for SNR we allow a 5% tollerance because of waveform sub-sampling
         assert np.allclose(metrics["snr"].values,
-                           metrics_inv["snr"].values, atol=1e-4)
+                           metrics_inv["snr"].values, rtol=0.05)
+        # for amplitude_cutoff, since spike amplitudes are computed, values should be exactly the same
         assert np.allclose(metrics["amplitude_cutoff"].values,
-                           metrics_inv["amplitude_cutoff"].values, atol=1e-4)
+                           metrics_inv["amplitude_cutoff"].values, atol=1e-5)
+
 
 if __name__ == '__main__':
     test = QualityMetricsExtensionTest()
     test.setUp()
     # test.test_extension()
-    test.test_metrics()
+    test.test_presence_ratio()
     # test.test_peak_sign()
diff --git a/spikeinterface/widgets/tests/test_widgets.py b/spikeinterface/widgets/tests/test_widgets.py
index 0252a7ef..e6fff8a3 100644
--- a/spikeinterface/widgets/tests/test_widgets.py
+++ b/spikeinterface/widgets/tests/test_widgets.py
@@ -50,10 +50,11 @@ class TestWidgets(unittest.TestCase):
         # at the moment only mpl will be tested on github actions
         sw.set_default_plotter_backend('matplotlib')
         
+        metric_names = ["snr", "isi_violations", "num_spikes"]
         _ = compute_spike_amplitudes(self.we)
         _ = compute_unit_locations(self.we)
         _ = compute_spike_locations(self.we)
-        _ = compute_quality_metrics(self.we)
+        _ = compute_quality_metrics(self.we, metric_names=metric_names)
         _ = compute_template_metrics(self.we)
         _ = compute_correlograms(self.we)
         _ = compute_template_similarity(self.we)
@@ -223,6 +224,8 @@ if __name__ == '__main__':
     # mytest.test_plot_unit_templates()
     # mytest.test_plot_unit_depths()
     # mytest.test_plot_unit_templates()
-    mytest.test_plot_unit_summary()
+    # mytest.test_plot_unit_summary()
+    mytest.test_quality_metrics()
+    mytest.test_sorting_summary()
 
     plt.show()