Merge branch 'main' into append-concat

Author: zm711

doc/api.rst

@@ -177,13 +177,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`.

src/spikeinterface/benchmark/tests/common_benchmark_testing.py

@@ -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/core/__init__.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

src/spikeinterface/core/analyzer_extension_core.py

@@ -9,6 +9,7 @@
   * ComputeNoiseLevels which is very convenient to have
 """
 
+import warnings
 import numpy as np
 
 from .sortinganalyzer import AnalyzerExtension, register_result_extension
@@ -92,6 +93,11 @@ def _merge_extension_data(
             new_data["random_spikes_indices"] = np.flatnonzero(selected_mask[keep_mask])
         return new_data
 
+    def _split_extension_data(self, split_units, new_unit_ids, new_sorting_analyzer, verbose=False, **job_kwargs):
+        new_data = dict()
+        new_data["random_spikes_indices"] = self.data["random_spikes_indices"].copy()
+        return new_data
+
     def _get_data(self):
         return self.data["random_spikes_indices"]
 
@@ -243,8 +249,6 @@ def _select_extension_data(self, unit_ids):
     def _merge_extension_data(
         self, merge_unit_groups, new_unit_ids, new_sorting_analyzer, keep_mask=None, verbose=False, **job_kwargs
     ):
-        new_data = dict()
-
         waveforms = self.data["waveforms"]
         some_spikes = self.sorting_analyzer.get_extension("random_spikes").get_random_spikes()
         if keep_mask is not None:
@@ -275,6 +279,11 @@ def _merge_extension_data(
 
         return dict(waveforms=waveforms)
 
+    def _split_extension_data(self, split_units, new_unit_ids, new_sorting_analyzer, verbose=False, **job_kwargs):
+        # splitting only affects random spikes, not waveforms
+        new_data = dict(waveforms=self.data["waveforms"].copy())
+        return new_data
+
     def get_waveforms_one_unit(self, unit_id, force_dense: bool = False):
         """
         Returns the waveforms of a unit id.
@@ -554,6 +563,49 @@ def _merge_extension_data(
 
         return new_data
 
+    def _split_extension_data(self, split_units, new_unit_ids, new_sorting_analyzer, verbose=False, **job_kwargs):
+        if not new_sorting_analyzer.has_extension("waveforms"):
+            warnings.warn(
+                "Splitting templates without the 'waveforms' extension will simply copy the template of the unit that "
+                "was split to the new split units. This is not recommended and may lead to incorrect results. It is "
+                "recommended to compute the 'waveforms' extension before splitting, or to use 'hard' splitting mode.",
+            )
+        new_data = dict()
+        for operator, arr in self.data.items():
+            # we first copy the unsplit units
+            new_array = np.zeros((len(new_sorting_analyzer.unit_ids), arr.shape[1], arr.shape[2]), dtype=arr.dtype)
+            new_analyzer_unit_ids = list(new_sorting_analyzer.unit_ids)
+            unsplit_unit_ids = [unit_id for unit_id in self.sorting_analyzer.unit_ids if unit_id not in split_units]
+            new_indices = np.array([new_analyzer_unit_ids.index(unit_id) for unit_id in unsplit_unit_ids])
+            old_indices = self.sorting_analyzer.sorting.ids_to_indices(unsplit_unit_ids)
+            new_array[new_indices, ...] = arr[old_indices, ...]
+
+            for split_unit_id, new_splits in zip(split_units, new_unit_ids):
+                if new_sorting_analyzer.has_extension("waveforms"):
+                    for new_unit_id in new_splits:
+                        split_unit_index = new_sorting_analyzer.sorting.id_to_index(new_unit_id)
+                        wfs = new_sorting_analyzer.get_extension("waveforms").get_waveforms_one_unit(
+                            new_unit_id, force_dense=True
+                        )
+
+                        if operator == "average":
+                            arr = np.average(wfs, axis=0)
+                        elif operator == "std":
+                            arr = np.std(wfs, axis=0)
+                        elif operator == "median":
+                            arr = np.median(wfs, axis=0)
+                        elif "percentile" in operator:
+                            _, percentile = operator.splot("_")
+                            arr = np.percentile(wfs, float(percentile), axis=0)
+                        new_array[split_unit_index, ...] = arr
+                else:
+                    split_unit_index = self.sorting_analyzer.sorting.id_to_index(split_unit_id)
+                    old_template = arr[split_unit_index, ...]
+                    new_indices = new_sorting_analyzer.sorting.ids_to_indices(new_splits)
+                    new_array[new_indices, ...] = np.tile(old_template, (len(new_splits), 1, 1))
+            new_data[operator] = new_array
+        return new_data
+
     def _get_data(self, operator="average", percentile=None, outputs="numpy"):
         if operator != "percentile":
             key = operator
@@ -646,7 +698,7 @@ def get_templates(self, unit_ids=None, operator="average", percentile=None, save
                 channel_ids=self.sorting_analyzer.channel_ids,
                 unit_ids=unit_ids,
                 probe=self.sorting_analyzer.get_probe(),
-                is_scaled=self.sorting_analyzer.return_in_uV,
+                is_in_uV=self.sorting_analyzer.return_in_uV,
             )
         else:
             raise ValueError("`outputs` must be 'numpy' or 'Templates'")
@@ -727,6 +779,10 @@ def _merge_extension_data(
         # this does not depend on units
         return self.data.copy()
 
+    def _split_extension_data(self, split_units, new_unit_ids, new_sorting_analyzer, verbose=False, **job_kwargs):
+        # this does not depend on units
+        return self.data.copy()
+
     def _run(self, verbose=False, **job_kwargs):
         self.data["noise_levels"] = get_noise_levels(
             self.sorting_analyzer.recording,