Merge branch 'main' into sa-by-dict

Author: chrishalcrow

doc/api.rst

@@ -176,13 +176,16 @@ spikeinterface.preprocessing
     .. autofunction:: clip
     .. autofunction:: common_reference
     .. autofunction:: correct_lsb
+    .. autofunction:: compute_motion
     .. autofunction:: correct_motion
     .. autofunction:: get_motion_presets
     .. autofunction:: get_motion_parameters_preset
     .. autofunction:: load_motion_info
     .. autofunction:: save_motion_info
     .. autofunction:: depth_order
     .. autofunction:: detect_bad_channels
+    .. autofunction:: detect_and_interpolate_bad_channels
+    .. autofunction:: detect_and_remove_bad_channels
     .. autofunction:: directional_derivative
     .. autofunction:: filter
     .. autofunction:: gaussian_filter

doc/how_to/index.rst

@@ -18,4 +18,5 @@ Guides on how to solve specific, short problems in SpikeInterface. Learn how to.
     auto_curation_training
     auto_curation_prediction
     physical_units
+    unsigned_to_signed
     customize_a_plot

doc/how_to/physical_units.rst

@@ -1,3 +1,5 @@
+.. _physical_units:
+
 Working with physical units in SpikeInterface recordings
 ========================================================
 

doc/how_to/unsigned_to_signed.rst

@@ -0,0 +1,102 @@
+Unsigned to Signed Data types
+=============================
+
+As of version 0.103.0 SpikeInterface has changed one of its defaults for interacting with
+:code:`Recording` objects. We no longer autocast unsigned dtypes to signed implicitly. This
+means that some users of SpikeInterface will need to add one additional line of code to their scripts
+to explicitly handle this conversion.
+
+
+Why this matters?
+-----------------
+
+For those that want a deeper understanding of dtypes `NumPy provides a great explanation <https://numpy.org/doc/stable/reference/arrays.dtypes.html>`_.
+For our purposes it is important to know that many pieces of recording equipment opt to store their electrophysiological data as unsigned integers
+(e.g., Intan, Maxwell Biosystems, 3Brain Biocam).
+Similarly to signed integers, in order to convert to real units these file formats only need to store a :code:`gain`
+and an :code:`offset`. Our :code:`RecordingExtractor`'s maintain the dtype that the file format utilizes, which means that some of our
+:code:`RecordingExtractor`'s will have unsigned dtypes.
+
+The problem with using unsigned dtypes is that many types of functions (including the ones we use from :code:`SciPy`) perform poorly with unsigned integers.
+This is made worse by the fact that these failures are silent (i.e. no error is triggered but the operation leads to nonsensical data). So the
+solution required is to convernt unsigned integers into signed integers. Previously we did this under the hood, automatically for users that had
+a :code:`Recording` object with an unsigned dtype.
+
+We decided, however, that implicitly performing this action was not the best course of action, since:
+
+1) *explicit* is always better than *implicit*
+2) some functions would *magically* change the dtype of the :code:`Recording` object, which can cause confusion
+
+So from version 0.103.0, users will now explicitly have to perform this transformation of their data. This will help users better understand how they are
+processing their data during an analysis pipeline as well as better understand the provenance of their pipeline.
+
+
+Using :code:`unsigned_to_signed`
+--------------------------------
+
+For users that receive an error because their :code:`Recording` is unsigned, their is one additional step that must be done:
+
+.. code:: python
+
+    import spikeinterface.extractors as se
+    import spikeinterface.preprocessing as spre
+
+    # Intan is an example of unsigned data
+    recording = se.read_intan('path/to/my/file.rhd', stream_id='0')
+    # to get a signed version of our Recording we use the following function
+    recording_signed = spre.unsigned_to_signed(recording)
+    # we can now apply any preprocessing functions like normal, e.g.
+    recording_filtered = spre.bandpass_filter(recording_signed)
+
+
+Now with the signed dtype of the :code:`Recording` one can use a SpikeInterface pipeline as usual.
+
+
+If you are curious if your :code:`Recording` is unsigned you can simply check the repr or use :code:`get_dtype()`
+
+.. code:: python
+
+    # the repr automatically displays the dtype
+    print(recording)
+    # use method on the Recording object
+    print(recording.get_dtype())
+
+In either case, if the dtype displayed has a :code:`u` at the beginning (e.g. :code:`uint16`) then your recording is
+unsigned. If it doesn't have the :code:`u` (e.g. :code:`int16`) then it is signed and would not need this preprocessing step.
+
+
+Bit depth
+---------
+
+One final important piece of information for some users is the concept of bit depth, which is the number of bits used to
+sample the data. The :code:`bit_depth` argument that can be fed into the :code:`unsigned_to_signed` function.
+This should be used in cases where the ADC bit depth does not match the bit depth of the data type (e.g., if the data is
+stored as :code:`uint16` but the ADC is 12 bits).
+Let's make a concrete example: the Biocam acquisition system from 3Brain uses a 12-bit ADC and stores the data as
+:code:`uint16`. This means that the data is stored in a 16-bit unsigned integer format, but the actual data
+only covers a 12-bit range. Therefore, that the "zero" of the data is not at 0, nor at half of the :code:`uint16` range (i.e. 2^15),
+but rather at 2048 (i.e., 2^12).
+In this case, setting the :code:`bit_depth` argument to 12 will allow the :code:`unsigned_to_signed` function to
+correctly convert the unsigned data to signed data and offset the data to be centered around 0, by subtracting 2048
+while converting the data from unsigned to signed.
+
+.. code:: python
+
+    recording_unsigned = se.read_biocam('path/to/my/file.brw')
+    # we can now convert to signed with the correct bit depth
+    recording_signed = spre.unsigned_to_signed(recording_unsigned, bit_depth=12)
+
+
+Additional Notes
+----------------
+
+1) Some sorters make use of SpikeInterface preprocessing either
+   within their wrappers or within their own code base. So remember to use the "signed" version of
+   your recording for the rest of your pipeline.
+
+2) Using :code:`unsigned_to_signed` in versions less than 0.103.0 does not hurt your scripts. This
+   option was available previously along with the implicit option. Adding this into scripts with old
+   versions of SpikeInterface will still work and will "future-proof" your scripts for when you
+   update to a version greater than or equal to 0.103.0.
+
+3) For additional information on units and scaling in SpikeInterface see :ref:`physical_units`.

doc/modules/qualitymetrics/example_cutoff.png

@@ -1,30 +1,128 @@
-Noise cutoff (not currently implemented)
-========================================
+Noise cutoff (:code:`noise_cutoff`)
+===================================
 
 Calculation
 -----------
 
 
-Metric describing whether an amplitude distribution is cut off, similar to _amp_cutoff  :ref:`amplitude cutoff <amp_cutoff>` but without a Gaussian assumption.
-A histogram of amplitudes is created and quantifies the distance between the low tail, mean number of spikes and high tail in terms of standard deviations.
+Metric describing whether an amplitude distribution is cut off as it approaches zero, similar to  :ref:`amplitude cutoff <amp_cutoff>` but without a Gaussian assumption.
 
-A SpikeInterface implementation is not yet available.
+The **noise cutoff** metric assesses whether a unit’s spike‐amplitude distribution is truncated
+at the low-end, which may be due to the high amplitude detection threshold in the deconvolution step,
+i.e., if low‐amplitude spikes were missed.  It does not assume a Gaussian shape;
+instead, it directly compares counts in the low‐amplitude bins to counts in high‐amplitude bins.
+
+1. **Build a histogram**
+
+   For each unit, divide all amplitudes into ``n_bins`` equally spaced bins over the range of the amplitude.
+   If the number of spikes is large, you may consider using a larger ``n_bins``. For a small number of spikes, consider a smaller ``n_bins``.
+   Let :math:`n_i` denote the count in the :math:`i`-th bin.
+
+2. **Identify the “low” region**
+    - Compute the amplitude value at the specified ``low_quantile`` (for example, 0.10 = 10th percentile), denoted as :math:`\text{amp}_{low}`.
+    - Find all histogram bins whose upper edge is below that quantile value.  These bins form the "low‐quantile region".
+    - Compute
+
+    .. math::
+        L_{\mathrm{bins}} = \bigl\{i : \text{upper_edge}_i \le \text{amp}_{low}\bigr\}, \quad
+        \mu_{\mathrm{low}} = \frac{1}{|L_{\mathrm{bins}}|}\sum_{i\in L_{\mathrm{bins}}} n_i.
+
+3. **Identify the “high” region**
+
+   - Compute the amplitude value at the specified ``high_quantile`` (for example, 0.25 = top 25th percentile), denoted as :math:`\text{amp}_{high}`.
+   - Find all histogram bins whose lower edge is greater than that quantile value.  These bins form the "high‐quantile region".
+   - Compute
+
+    .. math::
+        H_{\mathrm{bins}} &= \bigl\{i : \text{lower_edge}_i \ge \text{amp}_{high}\bigr\}, \\
+        \mu_{\mathrm{high}} &= \frac{1}{|H_{\mathrm{bins}}|}\sum_{i\in H_{\mathrm{bins}}} n_i, \quad
+        \sigma_{\mathrm{high}} = \sqrt{\frac{1}{|H_{\mathrm{bins}}|-1}\sum_{i\in H_{\mathrm{bins}}}\bigl(n_i-\mu_{\mathrm{high}} \bigr)^2}.
+
+4. **Compute cutoff**
+
+   The *cutoff* is given by how many standard deviations away the low-amplitude bins are from the high-amplitude bins, defined as
+
+    .. math::
+     \mathrm{cutoff} = \frac{\mu_{\mathrm{low}} - \mu_{\mathrm{high}}}{\sigma_{\mathrm{high}}}.
+
+
+   - If no low‐quantile bins exist, a warning is issued and ``cutoff = NaN``.
+   - If no high‐quantile bins exist or :math:`\sigma_{\mathrm{high}} = 0`, a warning is issued and ``cutoff = NaN``.
+
+5. **Compute the low-to-peak ratio**
+
+   - Let :math:`M = \max_i\,n_i` be the height of the largest bin in the histogram.
+   - Define
+
+   .. math::
+        \mathrm{ratio} = \frac{\mu_{\mathrm{low}}}{M}.
+
+
+   - If there are no low bins, :math:`\mathrm{ratio} = NaN`.
+
+
+Together, ``(cutoff, ratio)`` quantify how suppressed the low‐end of the amplitude distribution is relative to the top quantile and to the peak.
 
 Expectation and use
 -------------------
 
 Noise cutoff attempts to describe whether an amplitude distribution is cut off.
-The metric is loosely based on [Hill]_'s amplitude cutoff, but is here adapted (originally by [IBL]_) to avoid making the Gaussianity assumption on spike distributions.
-Noise cutoff provides an estimate of false negative rate, so a lower value indicates fewer missed spikes (a more complete unit).
+Larger values of ``cutoff`` and ``ratio`` suggest that the distribution is cut-off.
+IBL uses the default value of 1 (equivalent to e.g. ``low_quantile=0.01, n_bins=100``) to choose the number of
+lower bins, with a suggested threshold of 5 for ``cutoff`` to determine whether a unit is cut off or not.
+In practice, the IBL threshold is quite conservative, and a lower threshold might work better for your data.
+We suggest plotting the data using the  :py:func:`~spikeinterface.widgets.plot_amplitudes` widget to view your data when choosing your threshold.
+It is suggested to use this metric when the amplitude histogram is **unimodal**.
+
+The metric is loosely based on [Hill]_'s amplitude cutoff, but is here adapted (originally by [IBL2024]_) to avoid making the Gaussian assumption on spike distributions.
+
+Example code
+------------
+
+.. code-block:: python
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from spikeinterface.full as si
+
+    # Suppose `sorting_analyzer` has been computed with spike amplitudes:
+    # Select units you are interested in visualizing
+    unit_ids = ...
+
+    # Compute noise_cutoff:
+    summary_dict = compute_noise_cutoff(
+        sorting_analyzer=sorting_analyzer
+        high_quantile=0.25,
+        low_quantile=0.10,
+        n_bins=100,
+        unit_ids=unit_ids
+    )
+
+Reference
+---------
+
+.. autofunction:: spikeinterface.qualitymetrics.misc_metrics.compute_noise_cutoff
+
+Examples with plots
+-------------------
+
+Here is shown the histogram of two units, with the vertical lines separating low- and high-amplitude regions.
+
+- On the left, we have a unit with no truncation at the left end, and the cutoff and ratio are small.
+- On the right, we have a unit with truncation at -1, and the cutoff and ratio are much larger.
 
+.. image:: example_cutoff.png
+    :width: 600
 
 Links to original implementations
 ---------------------------------
 
 * From `IBL implementation <https://github.com/int-brain-lab/ibllib/blob/2e1f91c622ba8dbd04fc53946c185c99451ce5d6/brainbox/metrics/single_units.py>`_
 
+Note: Compared to the original implementation, we have added a comparison between the low-amplitude bins to the largest bin (``noise_ratio``).
+The selection of low-amplitude bins is based on the ``low_quantile`` rather than the number of bins.
 
 Literature
 ----------
 
-Metric introduced by [IBL]_ (adapted from [Hill]_'s amplitude cutoff metric).
+Metric introduced by [IBL2024]_.

doc/modules/qualitymetrics/noise_cutoff.rst

@@ -123,6 +123,8 @@ References
 
 .. [IBL] `Spike sorting pipeline for the International Brain Laboratory. 2022. <https://figshare.com/articles/online_resource/Spike_sorting_pipeline_for_the_International_Brain_Laboratory/19705522/3>`_
 
+.. [IBL2024] `Spike sorting pipeline for the International Brain Laboratory - Version 2. 2024. <https://figshare.com/articles/online_resource/Spike_sorting_pipeline_for_the_International_Brain_Laboratory/19705522?file=49783080>`_
+
 .. [Jackson] `Quantitative assessment of extracellular multichannel recording quality using measures of cluster separation. Society of Neuroscience Abstract. 2005. <https://www.sciencedirect.com/science/article/abs/pii/S0306452204008425>`_
 
 .. [Jain] `UnitRefine: A Community Toolbox for Automated Spike Sorting Curation. 2025 <https://www.biorxiv.org/content/10.1101/2025.03.30.645770v1>`_

doc/references.rst

@@ -202,9 +202,7 @@ docs = [
     "networkx",
     "skops", # For automated curation
     "scikit-learn", # For automated curation
-    # Download data
-    "pooch>=1.8.2",
-    "datalad>=1.0.2",
+    "huggingface_hub", # For automated curation
 
     # for release we need pypi, so this needs to be commented
     "probeinterface @ git+https://github.com/SpikeInterface/probeinterface.git",  # We always build from the latest version

pyproject.toml

@@ -74,7 +74,7 @@ def compute_gt_templates(recording, gt_sorting, ms_before=2.0, ms_after=3.0, ret
         channel_ids=recording.channel_ids,
         unit_ids=gt_sorting.unit_ids,
         probe=recording.get_probe(),
-        is_scaled=return_scaled,
+        is_in_uV=return_scaled,
     )
     return gt_templates
 

src/spikeinterface/benchmark/tests/common_benchmark_testing.py

@@ -109,7 +109,12 @@
     get_chunk_with_margin,
     order_channels_by_depth,
 )
-from .sorting_tools import spike_vector_to_spike_trains, random_spikes_selection, apply_merges_to_sorting
+from .sorting_tools import (
+    spike_vector_to_spike_trains,
+    random_spikes_selection,
+    apply_merges_to_sorting,
+    apply_splits_to_sorting,
+)
 
 from .waveform_tools import extract_waveforms_to_buffers, estimate_templates, estimate_templates_with_accumulator
 from .snippets_tools import snippets_from_sorting