Merge branch 'master' into fix-spikegadgets

Author: zm711

doc/source/releases.rst

@@ -6,6 +6,7 @@ Release notes
 .. toctree::
    :maxdepth: 1
 
+   releases/0.13.3.rst
    releases/0.13.2.rst
    releases/0.13.1.rst
    releases/0.13.0.rst

doc/source/releases/0.13.3.rst

@@ -0,0 +1,33 @@
+========================
+Neo 0.13.3 release notes
+========================
+
+28 August 2024
+
+This release of Neo contains bug fixes, still with a focus on the planned 1.0 release,
+and will be the last release not to support NumPy 2.0.
+
+See all `pull requests`_ included in this release and the `list of closed issues`_.
+
+
+Updated dependencies
+--------------------
+
+Neo has a limit of NumPy >= 1.19.5, < 2.0.0 and Quantities >= 14.0.1, < 0.16.0
+
+
+Bug fixes and improvements in IO modules
+----------------------------------------
+
+Bug fixes and/or improvements have been made to :class:`PlexonIO`, :class:`SpikeGLXIO`,
+:class:`BiocamIO`.
+
+Acknowledgements
+----------------
+
+Thanks to Zach McKenzie, Heberto Mayorquin, and Alessio Buccino for their contributions to this release.
+
+
+.. _`pull requests`: https://github.com/NeuralEnsemble/python-neo/pulls?q=is%3Apr+is%3Aclosed+milestone%3A0.13.3
+
+.. _`list of closed issues`: https://github.com/NeuralEnsemble/python-neo/issues?q=is%3Aissue+is%3Aclosed+milestone%3A0.13.3

neo/rawio/biocamrawio.py

@@ -17,6 +17,7 @@
     _spike_channel_dtype,
     _event_channel_dtype,
 )
+from neo.core import NeoReadWriteError
 
 
 class BiocamRawIO(BaseRawIO):
@@ -122,15 +123,53 @@ def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, strea
             i_start = 0
         if i_stop is None:
             i_stop = self._num_frames
-        if channel_indexes is None:
-            channel_indexes = slice(None)
+
+        # read functions are different based on the version of biocam
         data = self._read_function(self._filehandle, i_start, i_stop, self._num_channels)
-        return data[:, channel_indexes]
 
+        # older style data returns array of (n_samples, n_channels), should be a view
+        # but if memory issues come up we should doublecheck out how the file is being stored
+        if data.ndim > 1:
+            if channel_indexes is None:
+                channel_indexes = slice(None)
+            sig_chunk = data[:, channel_indexes]
+
+        # newer style data returns an initial flat array (n_samples * n_channels)
+        # we iterate through channels rather than slicing
+        # Due to the fact that Neo and SpikeInterface tend to prefer slices we need to add
+        # some careful checks around slicing of None in the case we need to iterate through
+        # channels. First check if None. Then check if slice and only if slice check that it is slice(None)
+        else:
+            if channel_indexes is None:
+                channel_indexes = [ch for ch in range(self._num_channels)]
+            elif isinstance(channel_indexes, slice):
+                start = channel_indexes.start or 0
+                stop = channel_indexes.stop or self._num_channels
+                step = channel_indexes.step or 1
+                channel_indexes = [ch for ch in range(start, stop, step)]
+
+            sig_chunk = np.zeros((i_stop - i_start, len(channel_indexes)), dtype=data.dtype)
+            # iterate through channels to prevent loading all channels into memory which can cause
+            # memory exhaustion. See https://github.com/SpikeInterface/spikeinterface/issues/3303
+            for index, channel_index in enumerate(channel_indexes):
+                sig_chunk[:, index] = data[channel_index :: self._num_channels]
+
+        return sig_chunk
 
-def open_biocam_file_header(filename):
+
+def open_biocam_file_header(filename) -> dict:
     """Open a Biocam hdf5 file, read and return the recording info, pick the correct method to access raw data,
-    and return this to the caller."""
+    and return this to the caller
+
+    Parameters
+    ----------
+    filename: str
+        The file to be parsed
+
+    Returns
+    -------
+    dict
+        The information necessary to read a biocam file (gain, n_samples, n_channels, etc)."""
     import h5py
 
     rf = h5py.File(filename, "r")
@@ -154,9 +193,9 @@ def open_biocam_file_header(filename):
         elif file_format in (101, 102) or file_format is None:
             num_channels = int(rf["3BData/Raw"].shape[0] / num_frames)
         else:
-            raise Exception("Unknown data file format.")
+            raise NeoReadWriteError("Unknown data file format.")
 
-        # # get channels
+        # get channels
         channels = rf["3BRecInfo/3BMeaStreams/Raw/Chs"][:]
 
         # determine correct function to read data
@@ -166,14 +205,14 @@ def open_biocam_file_header(filename):
             elif signal_inv == -1:
                 read_function = readHDF5t_100_i
             else:
-                raise Exception("Unknown signal inversion")
+                raise NeoReadWriteError("Unknown signal inversion")
         else:
             if signal_inv == 1:
                 read_function = readHDF5t_101
             elif signal_inv == -1:
                 read_function = readHDF5t_101_i
             else:
-                raise Exception("Unknown signal inversion")
+                raise NeoReadWriteError("Unknown signal inversion")
 
         gain = (max_uv - min_uv) / (2**bit_depth)
         offset = min_uv
@@ -200,19 +239,22 @@ def open_biocam_file_header(filename):
         scale_factor = experiment_settings["ValueConverter"]["ScaleFactor"]
         sampling_rate = experiment_settings["TimeConverter"]["FrameRate"]
 
+        num_channels = None
         for key in rf:
             if key[:5] == "Well_":
                 num_channels = len(rf[key]["StoredChIdxs"])
                 if len(rf[key]["Raw"]) % num_channels:
-                    raise RuntimeError(f"Length of raw data array is not multiple of channel number in {key}")
+                    raise NeoReadWriteError(f"Length of raw data array is not multiple of channel number in {key}")
                 num_frames = len(rf[key]["Raw"]) // num_channels
                 break
-        try:
+
+        if num_channels is not None:
             num_channels_x = num_channels_y = int(np.sqrt(num_channels))
-        except NameError:
-            raise RuntimeError("No Well found in the file")
+        else:
+            raise NeoReadWriteError("No Well found in the file")
+
         if num_channels_x * num_channels_y != num_channels:
-            raise RuntimeError(f"Cannot determine structure of the MEA plate with {num_channels} channels")
+            raise NeoReadWriteError(f"Cannot determine structure of the MEA plate with {num_channels} channels")
         channels = 1 + np.concatenate(np.transpose(np.meshgrid(range(num_channels_x), range(num_channels_y))))
 
         gain = scale_factor * (max_uv - min_uv) / (max_digital - min_digital)
@@ -231,6 +273,11 @@ def open_biocam_file_header(filename):
         )
 
 
+######################################################################
+# Helper functions to obtain the raw data split by Biocam version.
+
+
+# return the full array for the old datasets
 def readHDF5t_100(rf, t0, t1, nch):
     return rf["3BData/Raw"][t0:t1]
 
@@ -239,15 +286,16 @@ def readHDF5t_100_i(rf, t0, t1, nch):
     return 4096 - rf["3BData/Raw"][t0:t1]
 
 
+# return flat array that we will iterate through
 def readHDF5t_101(rf, t0, t1, nch):
-    return rf["3BData/Raw"][nch * t0 : nch * t1].reshape((t1 - t0, nch), order="C")
+    return rf["3BData/Raw"][nch * t0 : nch * t1]
 
 
 def readHDF5t_101_i(rf, t0, t1, nch):
-    return 4096 - rf["3BData/Raw"][nch * t0 : nch * t1].reshape((t1 - t0, nch), order="C")
+    return 4096 - rf["3BData/Raw"][nch * t0 : nch * t1]
 
 
 def readHDF5t_brw4(rf, t0, t1, nch):
     for key in rf:
         if key[:5] == "Well_":
-            return rf[key]["Raw"][nch * t0 : nch * t1].reshape((t1 - t0, nch), order="C")
+            return rf[key]["Raw"][nch * t0 : nch * t1]

neo/rawio/blackrockrawio.py

@@ -986,7 +986,8 @@ def __read_nsx_dataheader_variant_c(
         index = 0
 
         if offset is None:
-            offset = self.__nsx_basic_header[nsx_nb]["bytes_in_headers"]
+            # This is read as an uint32 numpy scalar from the header so we transform it to python int
+            offset = int(self.__nsx_basic_header[nsx_nb]["bytes_in_headers"])
 
         ptp_dt = [
             ("reserved", "uint8"),

neo/rawio/medrawio.py

@@ -240,9 +240,18 @@ def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, strea
             self.sess.set_channel_active(self._stream_info[stream_index]["raw_chans"])
             num_channels = len(self._stream_info[stream_index]["raw_chans"])
             self.sess.set_reference_channel(self._stream_info[stream_index]["raw_chans"][0])
+        
+        # in the case we have a slice or we give an ArrayLike we need to iterate through the channels
+        # in order to activate them.
         else:
-            if any(channel_indexes < 0):
-                raise IndexError(f"Can not index negative channels: {channel_indexes}")
+            if isinstance(channel_indexes, slice):
+                start = channel_indexes.start or 0
+                stop = channel_indexes.stop or len(self._stream_info[stream_index]["raw_chans"])
+                step = channel_indexes.step or 1
+                channel_indexes = [ch for ch in range(start, stop, step)]
+            else:
+                if any(channel_indexes < 0):
+                    raise IndexError(f"Can not index negative channels: {channel_indexes}")
             # Set all channels to be inactive, then selectively set some of them to be active
             self.sess.set_channel_inactive("all")
             for i, channel_idx in enumerate(channel_indexes):

neo/rawio/plexonrawio.py

@@ -24,6 +24,7 @@
 
 import datetime
 from collections import OrderedDict
+import re
 
 import numpy as np
 
@@ -146,7 +147,7 @@ def _parse_header(self):
 
             # Update tqdm with the number of bytes processed in this iteration
             if self.progress_bar:
-                progress_bar.update(length)
+                progress_bar.update(length)  # This was clever, Sam : )
 
         if self.progress_bar:
             progress_bar.close()
@@ -231,6 +232,7 @@ def _parse_header(self):
         # signals channels
         sig_channels = []
         all_sig_length = []
+        source_id = []
         if self.progress_bar:
             chan_loop = trange(nb_sig_chan, desc="Parsing signal channels", leave=True)
         else:
@@ -242,6 +244,7 @@ def _parse_header(self):
             length = self._data_blocks[5][chan_id]["size"].sum() // 2
             if length == 0:
                 continue  # channel not added
+            source_id.append(h["SrcId"])
             all_sig_length.append(length)
             sampling_rate = float(h["ADFreq"])
             sig_dtype = "int16"
@@ -255,7 +258,7 @@ def _parse_header(self):
                     0.5 * (2 ** global_header["BitsPerSpikeSample"]) * h["Gain"] * h["PreampGain"]
                 )
             offset = 0.0
-            stream_id = "0"
+            stream_id = "0"  # This is overwritten later
             sig_channels.append((name, str(chan_id), sampling_rate, sig_dtype, units, gain, offset, stream_id))
 
         sig_channels = np.array(sig_channels, dtype=_signal_channel_dtype)
@@ -264,22 +267,49 @@ def _parse_header(self):
             signal_streams = np.array([], dtype=_signal_stream_dtype)
 
         else:
-            # detect groups (aka streams)
-            all_sig_length = all_sig_length = np.array(all_sig_length)
-            groups = set(zip(sig_channels["sampling_rate"], all_sig_length))
+            # Detect streams
+            all_sig_length = np.asarray(all_sig_length)
+
+            # names are WB{number}, FPX{number}, SPKCX{number}, AI{number}, etc
+            pattern = r"^\D+"  # Match any non-digit character at the beginning of the string
+            channels_prefixes = np.asarray([re.match(pattern, name).group(0) for name in sig_channels["name"]])
+            buffer_stream_groups = set(zip(channels_prefixes, sig_channels["sampling_rate"], all_sig_length))
+
+            # There are explanations of the streams based on channel names
+            # provided by a Plexon Engineer, see here:
+            # https://github.com/NeuralEnsemble/python-neo/pull/1495#issuecomment-2184256894
+            channel_prefix_to_stream_name = {
+                "WB": "WB-Wideband",
+                "FP": "FPl-Low Pass Filtered ",
+                "SP": "SPKC-High Pass Filtered",
+                "AI": "AI-Auxiliary Input",
+            }
+
+            # Using a mapping to ensure consistent order of stream_index
+            channel_prefix_to_stream_id = {
+                "WB": "0",
+                "FP": "1",
+                "SP": "2",
+                "AI": "3",
+            }
 
             signal_streams = []
             self._signal_length = {}
             self._sig_sampling_rate = {}
-            for stream_index, (sr, length) in enumerate(groups):
-                stream_id = str(stream_index)
+
+            for stream_index, (channel_prefix, sr, length) in enumerate(buffer_stream_groups):
+                # The users of plexon can modify the prefix of the channel names (e.g. `my_prefix` instead of `WB`). This is not common but in that case
+                # We assign the channel_prefix both as stream_name and stream_id
+                stream_name = channel_prefix_to_stream_name.get(channel_prefix, channel_prefix)
+                stream_id = channel_prefix_to_stream_id.get(channel_prefix, channel_prefix)
+
                 mask = (sig_channels["sampling_rate"] == sr) & (all_sig_length == length)
                 sig_channels["stream_id"][mask] = stream_id
 
                 self._sig_sampling_rate[stream_index] = sr
                 self._signal_length[stream_index] = length
 
-                signal_streams.append(("Signals " + stream_id, stream_id))
+                signal_streams.append((stream_name, stream_id))
 
             signal_streams = np.array(signal_streams, dtype=_signal_stream_dtype)
 
@@ -363,8 +393,8 @@ def _segment_t_start(self, block_index, seg_index):
     def _segment_t_stop(self, block_index, seg_index):
         t_stop = float(self._last_timestamps) / self._global_ssampling_rate
         if hasattr(self, "_signal_length"):
-            for stream_id in self._signal_length:
-                t_stop_sig = self._signal_length[stream_id] / self._sig_sampling_rate[stream_id]
+            for stream_index in self._signal_length.keys():
+                t_stop_sig = self._signal_length[stream_index] / self._sig_sampling_rate[stream_index]
                 t_stop = max(t_stop, t_stop_sig)
         return t_stop
 

neo/test/iotest/test_plexonio.py

@@ -18,6 +18,7 @@ class TestPlexonIO(
         "plexon/File_plexon_1.plx",
         "plexon/File_plexon_2.plx",
         "plexon/File_plexon_3.plx",
+        "plexon/4chDemoPLX.plx"
     ]
 
 

neo/test/rawiotest/rawio_compliance.py

@@ -172,7 +172,7 @@ def read_analogsignals(reader):
         channel_names = signal_channels["name"][mask]
         channel_ids = signal_channels["id"][mask]
 
-        # acces by channel inde/ids/names should give the same chunk
+        # acces by channel index/ids/names should give the same chunk
         channel_indexes2 = channel_indexes[::2]
         channel_names2 = channel_names[::2]
         channel_ids2 = channel_ids[::2]
@@ -214,6 +214,29 @@ def read_analogsignals(reader):
             )
             np.testing.assert_array_equal(raw_chunk0, raw_chunk1)
 
+        # test slice(None). This should return the same array as giving
+        # all channel indexes or using None as an argument in `get_analogsignal_chunk`
+        # see https://github.com/NeuralEnsemble/python-neo/issues/1533
+
+        raw_chunk_slice_none = reader.get_analogsignal_chunk(
+            block_index=block_index,
+            seg_index=seg_index,
+            i_start=i_start,
+            i_stop=i_stop,
+            stream_index=stream_index,
+            channel_indexes=slice(None)
+        )
+        raw_chunk_channel_indexes = reader.get_analogsignal_chunk(
+            block_index=block_index,
+            seg_index=seg_index,
+            i_start=i_start,
+            i_stop=i_stop,
+            stream_index=stream_index,
+            channel_indexes=channel_indexes
+        )
+
+        np.testing.assert_array_equal(raw_chunk_slice_none, raw_chunk_channel_indexes)
+
         # test prefer_slice=True/False
         if nb_chan >= 3:
             for prefer_slice in (True, False):

neo/test/rawiotest/test_plexonrawio.py

@@ -15,6 +15,7 @@ class TestPlexonRawIO(
         "plexon/File_plexon_1.plx",
         "plexon/File_plexon_2.plx",
         "plexon/File_plexon_3.plx",
+        "plexon/4chDemoPLX.plx"
     ]