diff --git a/doc/source/authors.rst b/doc/source/authors.rst index 031523c60..f850324f4 100644 --- a/doc/source/authors.rst +++ b/doc/source/authors.rst @@ -104,6 +104,7 @@ and may not be the current affiliation of a contributor. * Sai Asish Yamani [51] * Kevin Doran [52] * Aditya Singh (github) +* 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 diff --git a/neo/rawio/baserawio.py b/neo/rawio/baserawio.py index 4207f3c03..ae0991c8d 100644 --- a/neo/rawio/baserawio.py +++ b/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() diff --git a/neo/rawio/maxwellrawio.py b/neo/rawio/maxwellrawio.py index 614f1e741..f9af17716 100644 --- a/neo/rawio/maxwellrawio.py +++ b/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,10 +186,15 @@ 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), @@ -173,13 +202,7 @@ def _parse_header(self): "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] @@ -187,15 +210,24 @@ def _parse_header(self): offset_uV = 0 buffer_id = stream_id sig_channels.append( - (ch_name, str(chan_id), sr, "uint16", "uV", gain_uV, offset_uV, stream_id, buffer_id) + ( + 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) diff --git a/neo/test/rawiotest/test_baserawio.py b/neo/test/rawiotest/test_baserawio.py new file mode 100644 index 000000000..c678a7d5c --- /dev/null +++ b/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() diff --git a/neo/test/rawiotest/test_maxwellrawio.py b/neo/test/rawiotest/test_maxwellrawio.py index aa4a30c7d..d5602165c 100644 --- a/neo/test/rawiotest/test_maxwellrawio.py +++ b/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()