Compute missing spikes on the right tail of amplitude distributions (#4353)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: ecobost

src/spikeinterface/metrics/quality/misc_metrics.py

@@ -873,21 +873,16 @@ def compute_amplitude_cutoffs(
 
     all_fraction_missing = {}
 
-    if sorting_analyzer.has_extension("spike_amplitudes"):
-        extension = sorting_analyzer.get_extension("spike_amplitudes")
-        all_amplitudes = extension.get_data()
-        invert_amplitudes = np.median(all_amplitudes) > 0
-    elif sorting_analyzer.has_extension("amplitude_scalings"):
-        # amplitude scalings are positive, we need to invert them
-        invert_amplitudes = True
-        extension = sorting_analyzer.get_extension("amplitude_scalings")
-
+    available_extension = (
+        "spike_amplitudes" if sorting_analyzer.has_extension("spike_amplitudes") else "amplitude_scalings"
+    )
+    extension = sorting_analyzer.get_extension(available_extension)
     amplitudes_by_units = extension.get_data(outputs="by_unit", concatenated=True, periods=periods)
 
     for unit_id in unit_ids:
         amplitudes = amplitudes_by_units[unit_id]
 
-        if invert_amplitudes:
+        if np.median(amplitudes) < 0:  # amplitude_cutoff expects positive amplitudes
             amplitudes = -amplitudes
         all_fraction_missing[unit_id] = amplitude_cutoff(
             amplitudes, num_histogram_bins, histogram_smoothing_value, amplitudes_bins_min_ratio
@@ -903,7 +898,7 @@ class AmplitudeCutoff(BaseMetric):
     metric_name = "amplitude_cutoff"
     metric_function = compute_amplitude_cutoffs
     metric_params = {
-        "num_histogram_bins": 100,
+        "num_histogram_bins": 200,
         "histogram_smoothing_value": 3,
         "amplitudes_bins_min_ratio": 5,
     }
@@ -1012,15 +1007,10 @@ def compute_noise_cutoffs(
     noise_cutoff_dict = {}
     noise_ratio_dict = {}
 
-    if sorting_analyzer.has_extension("spike_amplitudes"):
-        extension = sorting_analyzer.get_extension("spike_amplitudes")
-        all_amplitudes = extension.get_data()
-        invert_amplitudes = np.median(all_amplitudes) > 0
-    elif sorting_analyzer.has_extension("amplitude_scalings"):
-        # amplitude scalings are positive, we need to invert them
-        invert_amplitudes = True
-        extension = sorting_analyzer.get_extension("amplitude_scalings")
-
+    available_extension = (
+        "spike_amplitudes" if sorting_analyzer.has_extension("spike_amplitudes") else "amplitude_scalings"
+    )
+    extension = sorting_analyzer.get_extension(available_extension)
     amplitudes_by_units = extension.get_data(outputs="by_unit", concatenated=True, periods=periods)
 
     for unit_id in unit_ids:
@@ -1029,7 +1019,7 @@ def compute_noise_cutoffs(
             cutoff, ratio = np.nan, np.nan
             continue
 
-        if invert_amplitudes:
+        if np.median(amplitudes) < 0:  # _noise_cutoff expects positive amplitudes
             amplitudes = -amplitudes
 
         cutoff, ratio = _noise_cutoff(amplitudes, high_quantile=high_quantile, low_quantile=low_quantile, n_bins=n_bins)
@@ -1520,6 +1510,9 @@ def amplitude_cutoff(amplitudes, num_histogram_bins=500, histogram_smoothing_val
     """
     Calculate approximate fraction of spikes missing from a distribution of amplitudes.
 
+    Find the missing spikes from the left tail of the distribution. Assumes cutoff happens at spikes
+    with lower amplitudes.
+
     See compute_amplitude_cutoffs for additional documentation
 
     Parameters
@@ -1544,27 +1537,20 @@ def amplitude_cutoff(amplitudes, num_histogram_bins=500, histogram_smoothing_val
     if len(amplitudes) / num_histogram_bins < amplitudes_bins_min_ratio:
         return np.nan
     else:
-        h, b = np.histogram(amplitudes, num_histogram_bins, density=True)
-
-        # TODO : use something better than scipy.ndimage.gaussian_filter1d
         from scipy.ndimage import gaussian_filter1d
 
-        pdf = gaussian_filter1d(h, histogram_smoothing_value)
-        support = b[:-1]
-        bin_size = np.mean(np.diff(support))
-        peak_index = np.argmax(pdf)
+        # Approximate amplitude pdf with np.histogram
+        h = np.histogram(amplitudes, num_histogram_bins)[0]
+        pdf = gaussian_filter1d(h, histogram_smoothing_value, mode="nearest")
 
-        pdf_above = np.abs(pdf[peak_index:] - pdf[0])
-
-        if len(np.where(pdf_above == pdf_above.min())[0]) > 1:
-            warnings.warn(
-                "Amplitude PDF does not have a unique minimum! More spikes might be required for a correct "
-                "amplitude_cutoff computation!"
-            )
+        # Find number of missed spikes
+        cutoff_point = pdf[0]  # >> pdf[-1] if spikes were cutoff (at lower amplitudes)
+        G = np.where(pdf >= cutoff_point)[0][-1]  # last occurence where pdf was greater than cutoff
+        num_missed_spikes = np.sum(pdf[G + 1 :])  # theoretically missing spikes on the left side
 
-        G = np.argmin(pdf_above) + peak_index
-        fraction_missing = np.sum(pdf[G:]) * bin_size
-        fraction_missing = np.min([fraction_missing, 0.5])
+        # Compute fraction of missed spikes
+        fraction_missing = num_missed_spikes / (len(amplitudes) + num_missed_spikes)
+        fraction_missing = min(fraction_missing, 0.5)
 
         return fraction_missing
 

src/spikeinterface/metrics/quality/tests/test_metrics_functions.py

@@ -30,6 +30,7 @@
     compute_sd_ratio,
     _noise_cutoff,
     _get_synchrony_counts,
+    amplitude_cutoff,
 )
 
 from spikeinterface.metrics.quality.pca_metrics import (
@@ -130,6 +131,20 @@ def test_noise_cutoff():
     assert ratio1 <= ratio2
 
 
+def test_amplitude_cutoff():
+    """
+    Generate two artificial gaussians, one truncated and one not. Check the metrics are higher for the truncated one.
+    """
+    np.random.seed(1)
+    amps = np.random.normal(0, 1, 1000)
+    amps_trunc = amps[amps > -1]
+
+    fraction_missing_1 = amplitude_cutoff(amplitudes=amps, num_histogram_bins=20)
+    fraction_missing_2 = amplitude_cutoff(amplitudes=amps_trunc, num_histogram_bins=20)
+
+    assert fraction_missing_1 < fraction_missing_2
+
+
 def test_synchrony_counts_no_sync():
 
     spike_times, spike_units = synthesize_random_firings(num_units=1, duration=1, firing_rates=1.0)