Repository: SpikeInterface/spikeinterface (python)

Issue:
Error : Compressing Maxwell hdf5 dataformat to Zarr using SpikeInterface
Trying to compress maxwell hdf5 file to zarr format as discussed in the paper "Compression strategies for large-scale electrophysiology data".

Tried the following snippet :

`
from flac_numcodecs import Flac

compressor = Flac(level =8)

local_path= '/mnt/disk15tb/mmpatil/Spikesorting/Data/May16_analysison1024Amp/16848/Network/000041/data.raw.h5'
recording1 = se.read_maxwell(local_path)

recording_zarr = recording1.save(format = "zarr",folder="/mnt/disk15tb/mmpatil/Spikesorting/Data/May16_analysison1024Amp/16848/Network/000041/compressed.zarr",compressor=compressor,
                                channel_chunk_size =2,n_jobs=64,chunk_duration="1s")
`


But get the folowing Assertion error:


_RemoteTraceback                          Traceback (most recent call last)
_RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/usr/lib/python3.10/concurrent/futures/process.py", line 246, in _process_worker
    r = call_item.fn(*call_item.args, **call_item.kwargs)
  File "/usr/lib/python3.10/concurrent/futures/process.py", line 205, in _process_chunk
    return [fn(*args) for args in chunk]
  File "/usr/lib/python3.10/concurrent/futures/process.py", line 205, in 
    return [fn(*args) for args in chunk]
  File "/home/mmp/.local/lib/python3.10/site-packages/spikeinterface/core/job_tools.py", line 399, in function_wrapper
    return _func(segment_index, start_frame, end_frame, _worker_ctx)
  File "/home/mmp/.local/lib/python3.10/site-packages/spikeinterface/core/core_tools.py", line 702, in _write_zarr_chunk
    zarr_dataset[start_frame:end_frame, :] = traces
  File "/home/mmp/.local/lib/python3.10/site-packages/zarr/core.py", line 1391, in __setitem__
    self.set_basic_selection(pure_selection, value, fields=fields)
  File "/home/mmp/.local/lib/python3.10/site-packages/zarr/core.py", line 1486, in set_basic_selection
    return self._set_basic_selection_nd(selection, value, fields=fields)
  File "/home/mmp/.local/lib/python3.10/site-packages/zarr/core.py", line 1790, in _set_basic_selection_nd
    self._set_selection(indexer, value, fields=fields)
  File "/home/mmp/.local/lib/python3.10/site-packages/zarr/core.py", line 1842, in _set_selection
    self._chunk_setitem(chunk_coords, chunk_selection, chunk_value, fields=fields)
  File "/home/mmp/.local/lib/python3.10/site-packages/zarr/core.py", line 2137, in _chunk_setitem
    self._chunk_setitem_nosync(chunk_coords, chunk_selection, value,
...
    404     finally:
    405         # Break a reference cycle with the exception in self._exception
    406         self = None

AssertionError: Data type not supported. Only int16 is supported.

any guidance /support will be greatly appreciated.

Pull request:
Deal with unsigned int to int conversion
Fixes #1702 so you can do the following for a recording of typoe unsigned:

```
rec_int16 = rec_uint16.astype("int16")
```

Internally, it takes care of upcasting and removing the `2**(nbits-1)` offset before adding the sign (to avoid overflow errors).

Also added some tests ;)

Patch:
diff --git a/doc/api.rst b/doc/api.rst
index 4dbf1e10e..bc95281b9 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -141,13 +141,17 @@ spikeinterface.preprocessing
 
 .. automodule:: spikeinterface.preprocessing
 
+    .. autofunction:: astype
+    .. autofunction:: average_across_direction
     .. autofunction:: bandpass_filter
     .. autofunction:: blank_staturation
     .. autofunction:: center
     .. autofunction:: clip
     .. autofunction:: common_reference
     .. autofunction:: correct_lsb
+    .. autofunction:: depth_order
     .. autofunction:: detect_bad_channels
+    .. autofunction:: directional_derivative
     .. autofunction:: filter
     .. autofunction:: gaussian_bandpass_filter
     .. autofunction:: highpass_filter
@@ -158,7 +162,10 @@ spikeinterface.preprocessing
     .. autofunction:: phase_shift
     .. autofunction:: rectify
     .. autofunction:: remove_artifacts
+    .. autofunction:: resample
     .. autofunction:: scale
+    .. autofunction:: silence_periods
+    .. autofunction:: unsigned_to_signed
     .. autofunction:: whiten
     .. autofunction:: zero_channel_pad
     .. autofunction:: zscore
diff --git a/doc/modules/preprocessing.rst b/doc/modules/preprocessing.rst
index ed4f06bf1..bc851ff08 100644
--- a/doc/modules/preprocessing.rst
+++ b/doc/modules/preprocessing.rst
@@ -269,6 +269,29 @@ strategies:
 * :py:func:`~spikeinterface.preprocessing.remove_artifacts()`
 
 
+astype() / unsigned_to_signed()
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Similarly to :code:`numpy.astype()`, the :code:`astype()` casts the traces to the desired :code:`dtype`:
+
+.. code-block:: python
+
+    rec_int16 = astype(rec_float, "int16")
+
+
+For recordings whose traces are unsigned (e.g. Maxwell Biosystems), the :code:`unsigned_to_signed()` function makes them
+signed by removing the unsigned "offset". For example, :code:`uint16` traces will be first upcast to :code:`uint32`, 2**15
+is subtracted, and the traces are finally cast to :code:`int16`:
+
+
+.. code-block:: python
+
+    rec_int16 = unsigned_to_signed(rec_uint16)
+
+* :py:func:`~spikeinterface.preprocessing.astype()`
+* :py:func:`~spikeinterface.preprocessing.unsigned_to_signed()`
+
+
 zero_channel_pad()
 ^^^^^^^^^^^^^^^^^^
 
diff --git a/src/spikeinterface/preprocessing/astype.py b/src/spikeinterface/preprocessing/astype.py
index f79aabaab..2b6cff47d 100644
--- a/src/spikeinterface/preprocessing/astype.py
+++ b/src/spikeinterface/preprocessing/astype.py
@@ -7,6 +7,8 @@ class AstypeRecording(BasePreprocessor):
     """The spikeinterface analog of numpy.astype
 
     Converts a recording to another dtype on the fly.
+
+    For recording with an unsigned dtype, please use the `unsigned_to_signed` preprocessing function.
     """
 
     name = "astype"
diff --git a/src/spikeinterface/preprocessing/preprocessinglist.py b/src/spikeinterface/preprocessing/preprocessinglist.py
index a23aada52..de83d509a 100644
--- a/src/spikeinterface/preprocessing/preprocessinglist.py
+++ b/src/spikeinterface/preprocessing/preprocessinglist.py
@@ -36,6 +36,7 @@ from .average_across_direction import AverageAcrossDirectionRecording, average_a
 from .directional_derivative import DirectionalDerivativeRecording, directional_derivative
 from .depth_order import DepthOrderRecording, depth_order
 from .astype import AstypeRecording, astype
+from .unsigned_to_signed import UnsignedToSignedRecording, unsigned_to_signed
 
 
 preprocessers_full_list = [
@@ -70,6 +71,7 @@ preprocessers_full_list = [
     AverageAcrossDirectionRecording,
     DirectionalDerivativeRecording,
     AstypeRecording,
+    UnsignedToSignedRecording,
 ]
 
 installed_preprocessers_list = [pp for pp in preprocessers_full_list if pp.installed]
diff --git a/src/spikeinterface/preprocessing/unsigned_to_signed.py b/src/spikeinterface/preprocessing/unsigned_to_signed.py
new file mode 100644
index 000000000..bf032945f
--- /dev/null
+++ b/src/spikeinterface/preprocessing/unsigned_to_signed.py
@@ -0,0 +1,56 @@
+import numpy as np
+
+from ..core.core_tools import define_function_from_class
+from .basepreprocessor import BasePreprocessor, BasePreprocessorSegment
+from .filter import fix_dtype
+
+
+class UnsignedToSignedRecording(BasePreprocessor):
+    """
+    Converts a recording with unsigned traces to a signed one.
+    """
+
+    name = "unsigned_to_signed"
+
+    def __init__(
+        self,
+        recording,
+    ):
+        dtype = np.dtype(recording.dtype)
+        assert dtype.kind == "u", "Recording is not unsigned!"
+        itemsize = dtype.itemsize
+        assert itemsize < 8, "Cannot convert uint64 to int64."
+        dtype_signed = dtype.str.replace("uint", "int")
+
+        BasePreprocessor.__init__(self, recording, dtype=dtype_signed)
+
+        for parent_segment in recording._recording_segments:
+            rec_segment = UnsignedToSignedRecordingSegment(parent_segment, dtype_signed)
+            self.add_recording_segment(rec_segment)
+
+        self._kwargs = dict(
+            recording=recording,
+        )
+
+
+class UnsignedToSignedRecordingSegment(BasePreprocessorSegment):
+    def __init__(self, parent_recording_segment, dtype_signed):
+        BasePreprocessorSegment.__init__(self, parent_recording_segment)
+        self.dtype_signed = dtype_signed
+
+    def get_traces(self, start_frame, end_frame, channel_indices):
+        if channel_indices is None:
+            channel_indices = slice(None)
+        traces = self.parent_recording_segment.get_traces(start_frame, end_frame, channel_indices)
+        # if uint --> take care of offset
+        traces_dtype = traces.dtype
+        nbits = traces_dtype.itemsize * 8
+        signed_dtype = f"int{2 * (traces_dtype.itemsize) * 8}"
+        offset = 2 ** (nbits - 1)
+        # upcast to int with double itemsize
+        traces = traces.astype(signed_dtype, copy=False) - offset
+        return traces.astype(self.dtype_signed, copy=False)
+
+
+# function for API
+unsigned_to_signed = define_function_from_class(source_class=UnsignedToSignedRecording, name="unsigned_to_signed")


Test patch:
diff --git a/src/spikeinterface/preprocessing/tests/test_astype.py b/src/spikeinterface/preprocessing/tests/test_astype.py
new file mode 100644
index 000000000..f4ef2fe2c
--- /dev/null
+++ b/src/spikeinterface/preprocessing/tests/test_astype.py
@@ -0,0 +1,24 @@
+import pytest
+from pathlib import Path
+import shutil
+
+from spikeinterface import set_global_tmp_folder, NumpyRecording
+from spikeinterface.core import generate_recording
+
+from spikeinterface.preprocessing import astype
+
+import numpy as np
+
+
+def test_astype():
+    rng = np.random.RandomState(0)
+    traces = (rng.randn(10000, 4) * 100).astype("float32")
+    rec_float32 = NumpyRecording(traces, sampling_frequency=30000)
+    traces_int16 = traces.astype("int16")
+    np.testing.assert_array_equal(traces_int16, astype(rec_float32, "int16").get_traces())
+    traces_float64 = traces.astype("float64")
+    np.testing.assert_array_equal(traces_float64, astype(rec_float32, "float64").get_traces())
+
+
+if __name__ == "__main__":
+    test_astype()
diff --git a/src/spikeinterface/preprocessing/tests/test_unsigned_to_signed.py b/src/spikeinterface/preprocessing/tests/test_unsigned_to_signed.py
new file mode 100644
index 000000000..94d58cce9
--- /dev/null
+++ b/src/spikeinterface/preprocessing/tests/test_unsigned_to_signed.py
@@ -0,0 +1,29 @@
+import pytest
+from pathlib import Path
+import shutil
+
+from spikeinterface import set_global_tmp_folder, NumpyRecording
+from spikeinterface.core import generate_recording
+
+from spikeinterface.preprocessing import unsigned_to_signed
+
+import numpy as np
+
+
+def test_unsigned_to_signed():
+    rng = np.random.RandomState(0)
+    traces = rng.rand(10000, 4) * 100 + 2**15
+    traces_uint16 = traces.astype("uint16")
+    traces = rng.rand(10000, 4) * 100 + 2**31
+    traces_uint32 = traces.astype("uint32")
+    rec_uint16 = NumpyRecording(traces_uint16, sampling_frequency=30000)
+    rec_uint32 = NumpyRecording(traces_uint32, sampling_frequency=30000)
+
+    traces_int16 = (traces_uint16.astype("int32") - 2**15).astype("int16")
+    np.testing.assert_array_equal(traces_int16, unsigned_to_signed(rec_uint16).get_traces())
+    traces_int32 = (traces_uint32.astype("int64") - 2**31).astype("int32")
+    np.testing.assert_array_equal(traces_int32, unsigned_to_signed(rec_uint32).get_traces())
+
+
+if __name__ == "__main__":
+    test_unsigned_to_signed()