Valid unit periods extension (#4299)

Co-authored-by: Alessio Buccino <alejoe9187@gmail.com>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Chris Halcrow <57948917+chrishalcrow@users.noreply.github.com>
Co-authored-by: Samuel Garcia <sam.garcia.die@gmail.com>

Author: 5 people

doc/index.rst

@@ -30,7 +30,7 @@ SpikeInterface is made of several modules to deal with different aspects of the
 - visualize recordings and spike sorting outputs in several ways (matplotlib, sortingview, jupyter, ephyviewer)
 - export a report and/or export to phy
 - curate your sorting with several strategies (ml-based, metrics based, manual, ...)
-- offer a powerful Qt-based or we-based viewer in a separate package `spikeinterface-gui <https://github.com/SpikeInterface/spikeinterface-gui>`_ for manual curation that replace phy.
+- offer a powerful desktop or web viewer in a separate package `spikeinterface-gui <https://github.com/SpikeInterface/spikeinterface-gui>`_ for manual curation that replace phy.
 - have powerful sorting components to build your own sorter.
 - have a full motion/drift correction framework (See :ref:`motion_correction`)
 

doc/modules/core.rst

@@ -93,10 +93,11 @@ with 16 channels:
     timestamps = np.arange(num_samples) / sampling_frequency + 300
     recording.set_times(times=timestamps, segment_index=0)
 
-**Note**:
-Raw data formats often store data as integer values for memory efficiency. To give these integers meaningful physical units (uV), you can apply a gain and an offset.
-Many devices have their own gains and offsets necessary to convert their data and these values are handled by SpikeInterface for its extractors. This
-is triggered by the :code:`return_in_uV` parameter in :code:`get_traces()`, (see above example), which will return the traces in uV. Read more in our how to guide, :ref:`physical_units`.
+.. note::
+
+    Raw data formats often store data as integer values for memory efficiency. To give these integers meaningful physical units (uV), you can apply a gain and an offset.
+    Many devices have their own gains and offsets necessary to convert their data and these values are handled by SpikeInterface for its extractors. This
+    is triggered by the :code:`return_in_uV` parameter in :code:`get_traces()`, (see above example), which will return the traces in uV. Read more in our how to guide, :ref:`physical_units`.
 
 
 Sorting
@@ -180,8 +181,9 @@ a numpy.array with dtype `[("sample_index", "int64"), ("unit_index", "int64"), (
 For computations which are done unit-by-unit, like computing isi-violations per unit, it is better that
 spikes from a single unit are concurrent in memory. For these other cases, we can re-order the
 `spike_vector` in different ways:
-  * order by unit, then segment, then sample
-  * order by segment, then unit, then sample
+
+* order by unit, then segment, then sample
+* order by segment, then unit, then sample
 
 This is done using `sorting.to_reordered_spike_vector()`. The first time a reordering is done, the
 reordered spiketrain is cached in memory by default. Users should rarely have to worry about these
@@ -458,9 +460,11 @@ It represents unsorted waveform cutouts. Some acquisition systems, in fact, allo
 threshold and only record the times at which a peak was detected and the waveform cut out around
 the peak.
 
-**NOTE**: while we support this class (mainly for legacy formats), this approach is a bad practice
-and is highly discouraged! Most modern spike sorters, in fact, require the raw traces to perform
-template matching to recover spikes!
+.. note::
+
+    While we support this class (mainly for legacy formats), this approach is a bad practice
+    and is highly discouraged! Most modern spike sorters, in fact, require the raw traces to perform
+    template matching to recover spikes!
 
 Here we assume :code:`snippets` is a :py:class:`~spikeinterface.core.BaseSnippets` object
 with 16 channels:
@@ -548,9 +552,11 @@ Sparsity is defined as the subset of channels on which waveforms (and related in
 sparsity is not global, but it is unit-specific. Importantly, saving sparse waveforms, especially for high-density probes,
 dramatically reduces the size of the waveforms extension if computed.
 
-**NOTE** As of :code:`0.101.0` all :code:`SortingAnalyzer`'s have a default of :code:`sparse=True`. This was first
-introduced in :code:`0.99.0` for :code:`WaveformExtractor`'s and will be the default going forward. To obtain dense
-waveforms you will need to set :code:`sparse=False` at the creation of the :code:`SortingAnalyzer`.
+.. note::
+
+    As of :code:`0.101.0` all :code:`SortingAnalyzer`'s have a default of :code:`sparse=True`. This was first
+    introduced in :code:`0.99.0` for :code:`WaveformExtractor`'s and will be the default going forward. To obtain dense
+    waveforms you will need to set :code:`sparse=False` at the creation of the :code:`SortingAnalyzer`.
 
 
 Sparsity can be computed from a :py:class:`~spikeinterface.core.SortingAnalyzer` object with the
@@ -854,10 +860,12 @@ The same functions are also available for
 :py:func:`~spikeinterface.core.select_segment_sorting`).
 
 
-**Note** :py:func:`~spikeinterface.core.append_recordings` and:py:func:`~spikeinterface.core.concatenate_recordings`
-have the same goal, aggregate recording pieces on the time axis but with 2 different strategies! One is keeping the
-multi segments concept, the other one is breaking it!
-See this example for more detail :ref:`example_segments`.
+.. note::
+
+    :py:func:`~spikeinterface.core.append_recordings` and:py:func:`~spikeinterface.core.concatenate_recordings`
+    have the same goal, aggregate recording pieces on the time axis but with 2 different strategies! One is keeping the
+    multi segments concept, the other one is breaking it!
+    See this example for more detail :ref:`example_segments`.
 
 
 

doc/modules/exporters.rst

@@ -12,9 +12,11 @@ and behavioral data. It can be used to decode behavior, make tuning curves, comp
 The :py:func:`~spikeinterface.exporters.to_pynapple_tsgroup` function allows you to convert a
 SortingAnalyzer to Pynapple's ``TsGroup`` object on the fly.
 
-**Note** : When creating the ``TsGroup``, we will use the underlying time support of the SortingAnalyzer.
-How this works depends on your acquisition system. You can use the ``get_times`` method on a recording
-(``my_recording.get_times()``) to find the time support of your recording.
+.. note::
+
+    When creating the ``TsGroup``, we will use the underlying time support of the SortingAnalyzer.
+    How this works depends on your acquisition system. You can use the ``get_times`` method on a recording
+    (``my_recording.get_times()``) to find the time support of your recording.
 
 When constructed, if ``attach_unit_metadata`` is set to ``True``, any relevant unit information
 is propagated to the ``TsGroup``. The ``to_pynapple_tsgroup`` checks if unit locations, quality
@@ -54,13 +56,15 @@ The :py:func:`~spikeinterface.exporters.export_to_phy` function allows you to us
 `Phy template GUI <https://github.com/cortex-lab/phy>`_ for visual inspection and manual curation of spike sorting
 results.
 
-**Note** : :py:func:`~spikeinterface.exporters.export_to_phy` speed and the size of the folder will highly depend
-on the sparsity of the :code:`SortingAnalyzer` itself or the external specified sparsity.
-The Phy viewer enables one to explore PCA projections, spike amplitudes, waveforms and quality of spike sorting results.
-So if these pieces of information have already been computed as extensions (see :ref:`modules/postprocessing:Extensions as AnalyzerExtensions`),
-then exporting to Phy should be fast (and the user has better control of the parameters for the extensions).
-If not pre-computed, then the required extensions (e.g., :code:`spike_amplitudes`, :code:`principal_components`)
-can be computed directly at export time.
+.. note::
+
+    :py:func:`~spikeinterface.exporters.export_to_phy` speed and the size of the folder will highly depend
+    on the sparsity of the :code:`SortingAnalyzer` itself or the external specified sparsity.
+    The Phy viewer enables one to explore PCA projections, spike amplitudes, waveforms and quality of spike sorting results.
+    So if these pieces of information have already been computed as extensions (see :ref:`modules/postprocessing:Extensions as AnalyzerExtensions`),
+    then exporting to Phy should be fast (and the user has better control of the parameters for the extensions).
+    If not pre-computed, then the required extensions (e.g., :code:`spike_amplitudes`, :code:`principal_components`)
+    can be computed directly at export time.
 
 The input of the :py:func:`~spikeinterface.exporters.export_to_phy` is a :code:`SortingAnalyzer` object.
 
@@ -131,12 +135,14 @@ The report includes summary figures of the spike sorting output (e.g. amplitude
 depth VS amplitude) as well as unit-specific reports, that include waveforms, templates, template maps,
 ISI distributions, and more.
 
-**Note** : similarly to :py:func:`~spikeinterface.exporters.export_to_phy` the
-:py:func:`~spikeinterface.exporters.export_report` depends on the sparsity of the :code:`SortingAnalyzer` itself and
-on which extensions have been computed. For example, :code:`spike_amplitudes` and :code:`correlograms` related plots
-will be automatically included in the report if the associated extensions are computed in advance.
-The function can perform these computations as well, but it is a better practice to compute everything that's needed
-beforehand.
+.. note::
+
+    Similarly to :py:func:`~spikeinterface.exporters.export_to_phy` the
+    :py:func:`~spikeinterface.exporters.export_report` depends on the sparsity of the :code:`SortingAnalyzer` itself and
+    on which extensions have been computed. For example, :code:`spike_amplitudes` and :code:`correlograms` related plots
+    will be automatically included in the report if the associated extensions are computed in advance.
+    The function can perform these computations as well, but it is a better practice to compute everything that's needed
+    beforehand.
 
 Note that every unit will generate a summary unit figure, so the export process can be slow for spike sorting outputs
 with many units!

doc/modules/postprocessing.rst

@@ -163,8 +163,10 @@ Extensions are generally saved in two ways, suitable for two workflows:
    :code:`sorting_analyzer.compute('waveforms', save=False)`).
 
 
-**NOTE**: We recommend choosing a workflow and sticking with it. Either keep everything on disk or keep everything in memory until
-you'd like to save. A mixture can lead to unexpected behavior. For example, consider the following code
+.. note::
+
+    We recommend choosing a workflow and sticking with it. Either keep everything on disk or keep everything in memory until
+    you'd like to save. A mixture can lead to unexpected behavior. For example, consider the following code
 
 .. code::
 
@@ -257,15 +259,35 @@ spike_amplitudes
 This extension computes the amplitude of each spike as the value of the traces on the extremum channel at the times of
 each spike. The extremum channel is computed from the templates.
 
+
 **NOTE:** computing spike amplitudes is highly recommended before calculating amplitude-based quality metrics, such as
 :ref:`amp_cutoff` and :ref:`amp_median`.
 
 .. code-block:: python
 
-    amplitudes = sorting_analyzer.compute(input="spike_amplitudes", peak_sign="neg")
+    amplitudes = sorting_analyzer.compute(input="spike_amplitudes")
 
 For more information, see :py:func:`~spikeinterface.postprocessing.compute_spike_amplitudes`
 
+
+.. _postprocessing_amplitude_scalings:
+
+amplitude_scalings
+^^^^^^^^^^^^^^^^^^
+
+This extension computes the amplitude scaling of each spike as the value of the linear fit between the template and the
+spike waveform. In case of spatio-temporal collisions, a multi-linear fit is performed using the templates of all units
+involved in the collision.
+
+**NOTE:** computing amplitude scalings is highly recommended before calculating amplitude-based quality metrics, such as
+:ref:`amp_cutoff` and :ref:`amp_median`.
+
+.. code-block:: python
+
+    amplitude_scalings = sorting_analyzer.compute(input="amplitude_scalings")
+
+For more information, see :py:func:`~spikeinterface.postprocessing.compute_amplitude_scalings`
+
 .. _postprocessing_spike_locations:
 
 spike_locations
@@ -367,7 +389,29 @@ This extension computes the histograms of inter-spike-intervals. The computed ou
         method="auto"
     )
 
-For more information, see :py:func:`~spikeinterface.postprocessing.compute_isi_histograms`
+valid_unit_periods
+^^^^^^^^^^^^^^^^^^
+
+This extension computes the valid unit periods for each unit based on the estimation of false positive rates
+(using RP violation - see ::doc:`metrics/qualitymetrics/isi_violations`) and false negative rates
+(using amplitude cutoff - see ::doc:`metrics/qualitymetrics/amplitude_cutoff`) computed over chunks of the recording.
+The valid unit periods are the periods where both false positive and false negative rates are below specified
+thresholds. Periods can be either absolute (in seconds), same for all units, or relative, where
+chunks will be unit-specific depending on firing rate (with a target number of spikes per chunk).
+
+.. code-block:: python
+
+    valid_periods = sorting_analyzer.compute(
+        input="valid_unit_periods",
+        period_mode='relative',
+        target_num_spikes=300,
+        fp_threshold=0.1,
+        fn_threshold=0.1,
+    )
+
+For more information, see :py:func:`~spikeinterface.postprocessing.compute_valid_unit_periods`.
+
+
 
 
 Other postprocessing tools

doc/modules/sortingcomponents.rst

@@ -81,7 +81,9 @@ Other variants are also implemented (but less tested or not so useful):
 * **'by_channel_torch'** (requires :code:`torch`): pytorch implementation (GPU-compatible) that uses max pooling for time deduplication
 * **'locally_exclusive_torch'** (requires :code:`torch`): pytorch implementation (GPU-compatible) that uses max pooling for space-time deduplication
 
-**NOTE**: the torch implementations give slightly different results due to a different implementation.
+.. note::
+
+    The torch implementations give slightly different results due to a different implementation.
 
 Peak detection, as many of the other sorting components, can be run in parallel.
 
@@ -274,7 +276,7 @@ handle drift can benefit from drift estimation/correction.
 Especially for acute Neuropixels-like probes, this is a crucial step.
 
 The motion estimation step comes after peak detection and peak localization. Read more about
-it in the :ref:`_motion_correction` modules doc, and a more practical guide in the
+it in the :ref:`motion_correction` modules doc, and a more practical guide in the
 :ref:`handle-drift-in-your-recording` How To.
 
 Here is an example with non-rigid motion estimation:

doc/references.rst

@@ -118,6 +118,8 @@ References
 
 .. [Diggelmann] `Automatic spike sorting for high-density microelectrode arrays. 2018. <https://pubmed.ncbi.nlm.nih.gov/30207864/>`_
 
+.. [Fabre] `Bombcell: automated curation and cell classification of spike-sorted electrophysiology data. 2023. <https://doi.org/10.5281/zenodo.8172822>`
+
 .. [Garcia2024] `A Modular Implementation to Handle and Benchmark Drift Correction for High-Density Extracellular Recordings. 2024. <https://pubmed.ncbi.nlm.nih.gov/38238082/>`_
 
 .. [Garcia2022] `How Do Spike Collisions Affect Spike Sorting Performance? <https://doi.org/10.1523/ENEURO.0105-22.2022>`_