Fix Maxwell routed channel mapping

doc/source/authors.rst

@@ -106,6 +106,7 @@ and may not be the current affiliation of a contributor.
 * Aditya Singh (github)
 * Stanley C.
 * Axel Coiffard-Donaty
+* Arthur031221 (github)
 
 1. Centre de Recherche en Neuroscience de Lyon, CNRS UMR5292 - INSERM U1028 - Université Claude Bernard Lyon 1
 2. Unité de Neuroscience, Information et Complexité, CNRS UPR 3293, Gif-sur-Yvette, France

neo/rawio/baserawio.py

@@ -1695,9 +1695,6 @@ def _get_analogsignal_chunk(
             else:
                 raise RuntimeError("Should never happen")
 
-            if buffer_slice is not None:
-                raw_sigs = raw_sigs[:, buffer_slice]
-
         else:
             raise NotImplementedError()
 

neo/rawio/maxwellrawio.py

@@ -66,6 +66,33 @@ def __init__(self, filename="", rec_name=None):
     def _source_name(self):
         return self.filename
 
+    def _get_ids_and_electrodes(self, version, stream_id, h5file, mapping):
+        """Return channel IDs, electrode IDs and their signal-buffer indices."""
+        mapped_channel_ids = np.array(mapping["channel"])
+        electrode_ids = np.array(mapping["electrode"])
+        # Use the legacy format version already handled by _parse_header.
+        if int(version) == 20160704:
+            mask = mapped_channel_ids >= 0
+            return mapped_channel_ids[mask], electrode_ids[mask], slice(None)
+
+        routed = h5file["wells"][stream_id][self.rec_name]["groups"]["routed"]
+        channel_ids = np.array(routed["channels"])
+        routed_channel_ids_mask = np.isin(mapped_channel_ids, channel_ids[channel_ids >= 0])
+        unique_channel_ids, first_indices, counts = np.unique(
+            mapped_channel_ids[routed_channel_ids_mask], return_index=True, return_counts=True
+        )
+        duplicate_channel_ids = unique_channel_ids[counts > 1]
+        if duplicate_channel_ids.size:
+            warnings.warn(
+                f"Ignoring routed channels with multiple electrode mappings: {duplicate_channel_ids.tolist()}"
+            )
+        valid_channels = (channel_ids >= 0) & np.isin(channel_ids, unique_channel_ids[counts == 1])
+        buffer_indices = np.flatnonzero(valid_channels)
+        channel_ids = channel_ids[valid_channels]
+        channel_electrodes = dict(zip(unique_channel_ids, electrode_ids[routed_channel_ids_mask][first_indices]))
+        electrode_ids = np.array([channel_electrodes[channel_id] for channel_id in channel_ids])
+        return channel_ids, electrode_ids, buffer_indices
+
     def _parse_header(self):
         import h5py
 
@@ -139,9 +166,6 @@ def _parse_header(self):
                     gain_uV = 3.3 / (1024 * gain) * 1e6
                 hdf5_path = "sig"
                 mapping = h5file["mapping"]
-                ids = np.array(mapping["channel"])
-                ids = ids[ids >= 0]
-                self._stream_buffer_slice[stream_id] = ids
             elif int(version) > 20160704:
                 settings = h5file["wells"][stream_id][self.rec_name]["settings"]
                 sr = settings["sampling"][0]
@@ -162,24 +186,23 @@ def _parse_header(self):
                     well_indices_to_remove.append(stream_index)
                     continue
 
-                self._stream_buffer_slice[stream_id] = None
-
             buffer_id = stream_id
             shape = h5file[hdf5_path].shape
+            max_sig_length = max(max_sig_length, shape[1])
+            channel_ids, electrode_ids, buffer_indices = self._get_ids_and_electrodes(
+                version, stream_id, h5file, mapping
+            )
+            if channel_ids.size == 0:
+                well_indices_to_remove.append(stream_index)
+                continue
             self._buffer_descriptions[0][0][buffer_id] = {
                 "type": "hdf5",
                 "file_path": str(self.filename),
                 "hdf5_path": hdf5_path,
                 "shape": shape,
                 "time_axis": 1,
             }
-            self._stream_buffer_slice[stream_id] = slice(None)
-
-            channel_ids = np.array(mapping["channel"])
-            electrode_ids = np.array(mapping["electrode"])
-            mask = channel_ids >= 0
-            channel_ids = channel_ids[mask]
-            electrode_ids = electrode_ids[mask]
+            self._stream_buffer_slice[stream_id] = buffer_indices
 
             for i, chan_id in enumerate(channel_ids):
                 elec_id = electrode_ids[i]
@@ -190,12 +213,11 @@ def _parse_header(self):
                     (ch_name, str(chan_id), sr, "uint16", "uV", gain_uV, offset_uV, stream_id, buffer_id)
                 )
 
-            max_sig_length = max(max_sig_length, shape[1])
-
         self._t_stop = max_sig_length / sr
 
         if len(well_indices_to_remove) > 0:
             signal_streams = np.delete(signal_streams, np.array(well_indices_to_remove))
+            signal_buffers = np.delete(signal_buffers, np.array(well_indices_to_remove))
 
         sig_channels = np.array(sig_channels, dtype=_signal_channel_dtype)
 

neo/test/rawiotest/test_baserawio.py

@@ -0,0 +1,57 @@
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+try:
+    import h5py
+except ImportError:
+    h5py = None
+
+from neo.rawio.baserawio import BaseRawWithBufferApiIO, _signal_stream_dtype
+
+
+class BufferReader(BaseRawWithBufferApiIO):
+    def _get_analogsignal_buffer_description(self, block_index, seg_index, buffer_id):
+        return self.description
+
+
+@unittest.skipIf(h5py is None, "requires h5py")
+class TestHDF5BufferSelection(unittest.TestCase):
+    def test_stream_selection_applied_once(self):
+        samples = np.arange(30, dtype="uint16").reshape(6, 5)
+        with tempfile.TemporaryDirectory() as directory:
+            filename = Path(directory) / "signals.h5"
+            with h5py.File(filename, "w") as file:
+                file["time_first"] = samples
+                file["channel_first"] = samples.T
+            for time_axis, dataset in [(0, "time_first"), (1, "channel_first")]:
+                for buffer_slice in [
+                    None,
+                    slice(0, -1),
+                    slice(-1, None),
+                    np.array([0, 2, 4]),
+                    np.array([True, False, True, False, True]),
+                ]:
+                    with self.subTest(time_axis=time_axis, buffer_slice=buffer_slice):
+                        reader = BufferReader()
+                        reader.header = {
+                            "signal_streams": np.array([("signals", "0", "0")], dtype=_signal_stream_dtype)
+                        }
+                        reader.description = {
+                            "type": "hdf5",
+                            "file_path": str(filename),
+                            "hdf5_path": dataset,
+                            "time_axis": time_axis,
+                        }
+                        reader._stream_buffer_slice = {"0": buffer_slice}
+                        expected = samples[1:4]
+                        if buffer_slice is not None:
+                            expected = expected[:, buffer_slice]
+                        for channel_indexes in [None, [0], slice(0, 1), [-1, 0]]:
+                            selected = expected if channel_indexes is None else expected[:, channel_indexes]
+                            np.testing.assert_array_equal(
+                                reader._get_analogsignal_chunk(0, 0, 1, 4, 0, channel_indexes), selected
+                            )
+                        reader._hdf5_analogsignal_buffers[0][0]["0"].close()

neo/test/rawiotest/test_maxwellrawio.py

@@ -1,5 +1,16 @@
 import unittest
+import tempfile
+import warnings
+from pathlib import Path
 
+import numpy as np
+
+try:
+    import h5py
+except ImportError:
+    h5py = None
+
+from neo.io import MaxwellIO
 from neo.rawio.maxwellrawio import MaxwellRawIO, auto_install_maxwell_hdf5_compression_plugin
 from neo.test.rawiotest.common_rawio_test import BaseTestRawIO
 
@@ -20,5 +31,121 @@ def setUp(self):
         auto_install_maxwell_hdf5_compression_plugin(force_download=False)
 
 
+@unittest.skipIf(h5py is None, "requires h5py")
+class TestMaxwellRawIOSynthetic(unittest.TestCase):
+    def _check_mapping(self, mapping, routed_ids, channel_ids, electrode_ids, rows, legacy=False):
+        samples = np.array([[10 * (i + 1) + j for j in range(4)] for i in range(len(routed_ids))], dtype="uint16")
+        with tempfile.TemporaryDirectory() as directory:
+            filename = Path(directory) / "recording.h5"
+            with h5py.File(filename, "w") as file:
+                file["version"] = [b"20160704" if legacy else b"20190530"]
+                if legacy:
+                    settings = file.create_group("settings")
+                    file["sig"] = samples
+                    mapping_group = file
+                else:
+                    recording = file.create_group("wells/well000/rec0000")
+                    settings = recording.create_group("settings")
+                    mapping_group = settings
+                    routed = recording.create_group("groups/routed")
+                    routed["channels"] = routed_ids
+                    routed["raw"] = samples
+                settings["sampling"] = [20000.0]
+                settings["lsb"] = [1e-6]
+                mapping_group["mapping"] = np.array(mapping, dtype=[("channel", "i4"), ("electrode", "i4")])
+
+            reader = MaxwellRawIO(filename=str(filename))
+            try:
+                with warnings.catch_warnings(record=True) as caught:
+                    warnings.simplefilter("always")
+                    reader.parse_header()
+                channels = reader.header["signal_channels"]
+                self.assertEqual(channels["id"].tolist(), [str(channel) for channel in channel_ids])
+                self.assertEqual(
+                    channels["name"].tolist(),
+                    [f"ch{channel} elec{electrode}" for channel, electrode in zip(channel_ids, electrode_ids)],
+                )
+                duplicated_ids = {
+                    channel
+                    for channel, _ in mapping
+                    if channel >= 0 and sum(other == channel for other, _ in mapping) > 1
+                }
+                if duplicated_ids.intersection(routed_ids) and not legacy:
+                    self.assertTrue(any("multiple electrode mappings" in str(w.message) for w in caught))
+                if not channel_ids:
+                    self.assertEqual(reader.signal_streams_count(), 0)
+                    self.assertEqual(reader.header["signal_buffers"].size, 0)
+                expected = samples[rows, :].T
+                if channel_ids:
+                    np.testing.assert_array_equal(reader.get_analogsignal_chunk(), expected)
+                    np.testing.assert_array_equal(reader.get_analogsignal_chunk(i_start=1, i_stop=3), expected[1:3])
+                    np.testing.assert_array_equal(
+                        reader.get_analogsignal_chunk(channel_indexes=slice(0, 1)), expected[:, :1]
+                    )
+                    for i, channel in enumerate(channel_ids):
+                        np.testing.assert_array_equal(
+                            reader.get_analogsignal_chunk(channel_ids=[str(channel)]),
+                            expected[:, i : i + 1],
+                        )
+                    if len(channel_ids) > 1:
+                        np.testing.assert_array_equal(
+                            reader.get_analogsignal_chunk(channel_indexes=[len(channel_ids) - 1, 0]),
+                            expected[:, [-1, 0]],
+                        )
+            finally:
+                reader.h5_file.close()
+
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore", UserWarning)
+                io = MaxwellIO(filename=str(filename))
+            try:
+                segment = io.read_segment()
+                if channel_ids:
+                    np.testing.assert_array_equal(segment.analogsignals[0].magnitude, expected)
+                else:
+                    self.assertEqual(len(segment.analogsignals), 0)
+            finally:
+                io.h5_file.close()
+
+    def test_routed_channel_mapping(self):
+        cases = [
+            ([(2, 20), (5, 50)], [2, 5], [2, 5], [20, 50], [0, 1]),
+            ([(2, 20), (5, 50)], [2, 8, 5], [2, 5], [20, 50], [0, 2]),
+            ([(2, 20), (5, 50)], [8, 5, 2], [5, 2], [50, 20], [1, 2]),
+            ([(8, 80)], [2, 5], [], [], []),
+            ([(2, 20), (2, 21), (5, 50)], [2, 5], [5], [50], [1]),
+            ([(2, 20), (2, 21), (5, 50)], [5, 2], [5], [50], [0]),
+            ([(-1, 99), (2, 20), (2, 21), (5, 50), (8, 80)], [2, 5], [5], [50], [1]),
+            ([(8, 80), (5, 50), (5, 51), (2, 20), (-1, 99)], [5, 2], [2], [20], [1]),
+            ([(2, 20), (5, 20)], [5, 2], [5, 2], [20, 20], [0, 1]),
+            ([(2, 20), (2, 21), (5, 50), (5, 51)], [2, 5], [], [], []),
+            ([(2, 20), (5, 50), (-1, 99)], [-1, 5, 2], [5, 2], [50, 20], [1, 2]),
+            ([(2, 20), (5, 50), (8, 80), (8, 81)], [5, 2], [5, 2], [50, 20], [0, 1]),
+            (
+                [
+                    (399, 25732),
+                    (431, 25732),
+                    (399, 13364),
+                    (431, 13364),
+                    (9, 90),
+                    (7, 70),
+                    (5, 50),
+                ],
+                [9, 399, 7, 431, 5],
+                [9, 7, 5],
+                [90, 70, 50],
+                [0, 2, 4],
+            ),
+        ]
+        for mapping, routed_ids, channel_ids, electrode_ids, rows in cases:
+            with self.subTest(mapping=mapping, routed_ids=routed_ids):
+                self._check_mapping(mapping, routed_ids, channel_ids, electrode_ids, rows)
+
+    def test_legacy_channel_mapping(self):
+        for mapping in [[(0, 20), (1, 50)], [(-1, 99), (0, 20), (1, 50)]]:
+            with self.subTest(mapping=mapping):
+                self._check_mapping(mapping, [0, 1], [0, 1], [20, 50], [0, 1], legacy=True)
+
+
 if __name__ == "__main__":
     unittest.main()