From 989423851488cb91a5cc639d32ba50bc801d077b Mon Sep 17 00:00:00 2001 From: Kevin <1439017+kevindoran@users.noreply.github.com> Date: Thu, 23 Jul 2026 17:59:22 +0100 Subject: [PATCH 1/2] Fix #3303, and the performance regression the previous fix attempt introduced. Previous comments (starting around 1814bc1e) attempted to fix issue #3303, but in doing so, made a performance regression by forcing an array copy. It causes ~25x slowdown on my machine. Here we remove the regression, and properly chunk-read in the event of smaller subsets of large arrays being requested. This also fixed #3303. --- neo/rawio/biocamrawio.py | 77 ++++++++------ neo/test/rawiotest/test_biocamrawio.py | 142 +++++++++++++++++++++++++ 2 files changed, 186 insertions(+), 33 deletions(-) diff --git a/neo/rawio/biocamrawio.py b/neo/rawio/biocamrawio.py index 3aac8d2fc..43e58dcf6 100644 --- a/neo/rawio/biocamrawio.py +++ b/neo/rawio/biocamrawio.py @@ -25,6 +25,10 @@ from neo.core import NeoReadWriteError +# Maximum bytes to read at once. +_MAX_READ_BYTES = 64 * 1024 * 1024 # 64 MB + + class BiocamRawIO(BaseRawIO): """ Class for reading data from a Biocam h5 file. @@ -136,12 +140,8 @@ def _get_signal_t_start(self, block_index, seg_index, stream_index): raise ValueError("`stream_index must be 0") return self._segment_t_start(block_index, seg_index) - def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, stream_index, channel_indexes): - if i_start is None: - i_start = 0 - if i_stop is None: - i_stop = self._num_frames - + def _read_block(self, i_start, i_stop): + """Read frames [i_start, i_stop) for all channels, as (n_frames, n_channels).""" # read functions are different based on the version of biocam if self._read_function is readHDF5t_brw4_sparse: if self._fill_gaps_strategy is None: @@ -163,34 +163,45 @@ def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, strea else: data = self._read_function(self._filehandle, i_start, i_stop, self._num_channels) - # 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) + # Newer style data is read as a flat (n_samples * n_channels) array. The h5py read is + # C-contiguous, so this reshape is a free view: it copies nothing and allocates nothing. + # Older style data is already (n_samples, n_channels). + if data.ndim == 1: + data = data.reshape(i_stop - i_start, self._num_channels) + return data + + def _get_analogsignal_chunk(self, block_index, seg_index, i_start, i_stop, stream_index, channel_indexes): + if i_start is None: + i_start = 0 + if i_stop is None: + i_stop = self._num_frames + n_frames = i_stop - i_start + # Read a single sample to determine the number of selected channels. + if channel_indexes is None: + channel_indexes = slice(None) + n_ch = self._read_block(i_start, i_start + 1)[:, channel_indexes].shape[1] + dtype = self.header["signal_channels"]["dtype"][0] + itemsize = np.dtype(dtype).itemsize + read_bytes = n_frames * self._num_channels * itemsize + out_bytes = n_frames * n_ch * itemsize + + # Reading a subset of channels from a very large recording. + # The BioCam hdf5 layout is such that all channels are chunked together, so they cannot be + # read independently. In the event that only a subset of channels of a very large recording + # is needed, we cannot slice a single channel. Instead, we must loop and read all channels, + # slicing as we go. We will do this chunking if the read is large and the output to be + # returned is smaller. + read_by_chunks = read_bytes > _MAX_READ_BYTES and out_bytes < read_bytes // 2 + if read_by_chunks: + # Read in chunks, copying desired channels into a preallocated array. + out = np.empty((n_frames, n_ch), dtype=dtype) + max_frames = max(1, _MAX_READ_BYTES // (itemsize * self._num_channels)) + for start in range(i_start, i_stop, max_frames): + stop = min(start + max_frames, i_stop) + out[start - i_start : stop - i_start] = self._read_block(start, stop)[:, channel_indexes] 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 + out = self._read_block(i_start, i_stop)[:, channel_indexes] + return out def open_biocam_file_header(filename) -> dict: diff --git a/neo/test/rawiotest/test_biocamrawio.py b/neo/test/rawiotest/test_biocamrawio.py index b7dd0c906..6cb50c5d9 100644 --- a/neo/test/rawiotest/test_biocamrawio.py +++ b/neo/test/rawiotest/test_biocamrawio.py @@ -4,6 +4,11 @@ import unittest +import numpy as np +import pytest +import h5py + +from neo.rawio import biocamrawio from neo.rawio.biocamrawio import BiocamRawIO from neo.test.rawiotest.common_rawio_test import BaseTestRawIO @@ -23,5 +28,142 @@ class TestBiocamRawIO( ] +def write_minimal_brw(path, n_ch, n_frames, version=102): + """Write a minimal Biocam file whose samples are all distinct. + + Sample `s` of channel `c` holds the value `s * n_ch + c`, so an error in reshaping, + channel selection or block stitching produces wrong values rather than plausible ones. + Version 102 stores the raw data as a flat array, version 101 stores it 2-D. + + BioCam 3.x stores the 3BData group's format in its "Version" attribute. + v100 lays Raw out as a (T, W) matrix; v101 as a flat (T*W, 1) interleaved array. + + Info on 100-102 file formats: + + https://resources.3brain.com/f1/downloads/BrainWave_Documentation_for_BRWv.%203.x_and_BXRv.%202.x_FileFormats_v1.2.0.pdf + + Returns the expected (n_frames, n_ch) contents. + + TODO: this can be shared with the parallel pull request for issue #1883. + """ + values = np.arange(n_frames * n_ch, dtype=np.uint16) + with h5py.File(path, "w") as f: + rec_vars = f.create_group("3BRecInfo/3BRecVars") + rec_vars.create_dataset("BitDepth", data=np.array([12], dtype=np.uint8)) + rec_vars.create_dataset("MaxVolt", data=np.array([4125.0])) + rec_vars.create_dataset("MinVolt", data=np.array([-4125.0])) + rec_vars.create_dataset("NRecFrames", data=np.array([n_frames], dtype=np.int64)) + rec_vars.create_dataset("SamplingRate", data=np.array([17852.77])) + rec_vars.create_dataset("SignalInversion", data=np.array([1], dtype=np.int32)) + f.create_dataset("3BRecInfo/3BMeaStreams/Raw/Chs", data=np.arange(2 * n_ch, dtype=np.int32).reshape(n_ch, 2)) + f.create_dataset("3BData/Raw", data=values.reshape(n_frames, n_ch) if version == 100 else values) + f["3BData"].attrs["Version"] = version + frames = values.reshape(n_frames, n_ch) + return frames + + +class TestGetAnalogsignalChunk: + """Tests for BiocamRawIO._get_analogsignal_chunk. + + It's a little tricky to test this function, as what we care about is not input or output but performance. + A reasonably effective way to inspect the behaviour is to count how many times the inner _read block is called, + which is the function that actually reads the raw file. This is is done for the tests of this class. To do it, we + replace the private _read_block method of BiocamRawIO with a custom method that tracks its calls. + """ + + n_ch = 8 + n_frames = 20 + # Some ways neo and SpikeInterface are known to ask for channels: + channel_queries = [None, slice(None), slice(2, 6), slice(0, n_ch, 2), [0, 1, 2], [3, 0, 2], [5], np.array([1, 4])] + # None means "the whole segment"; the rest include a non-zero start and a single frame. + frame_ranges = [(None, None), (0, 5), (7, 13), (n_frames - 1, n_frames)] + # 64 bytes is 4 frames of 8 uint16 channels, so requests split into several blocks, + # the last of which is partial for most of the ranges above. + small_budget = 64 + small_budget_frames = 4 + + def make_reader(self, tmp_path, version=101): + """Return a BiocamRawIO , and the expected frames.""" + path = tmp_path / "minimal.brw" + frames = write_minimal_brw(path, self.n_ch, self.n_frames, version) + reader = BiocamRawIO(filename=path) + reader.parse_header() + return reader, frames + + @pytest.fixture + def frames_per_read(self, monkeypatch): + """Record how many frames each raw read covers.""" + reads = [] + original = BiocamRawIO._read_block + + def spy(self, i_start, i_stop): + reads.append(i_stop - i_start) + return original(self, i_start, i_stop) + + monkeypatch.setattr(BiocamRawIO, "_read_block", spy) + return reads + + @pytest.mark.parametrize("version", [100, 101]) + @pytest.mark.parametrize("channel_indexes", channel_queries) + @pytest.mark.parametrize("i_start, i_stop", frame_ranges) + def test_values(self, tmp_path, version, channel_indexes, i_start, i_stop): + """Values, shape and dtype are correct for every way of selecting channels.""" + reader, expected_all = self.make_reader(tmp_path, version) + + sig = reader.get_analogsignal_chunk( + block_index=0, seg_index=0, i_start=i_start, i_stop=i_stop, channel_indexes=channel_indexes + ) + + start = 0 if i_start is None else i_start + stop = self.n_frames if i_stop is None else i_stop + expected = expected_all[start:stop][:, slice(None) if channel_indexes is None else channel_indexes] + assert sig.dtype == expected.dtype + np.testing.assert_array_equal(sig, expected) + + @pytest.mark.parametrize("channel_indexes", channel_queries) + @pytest.mark.parametrize("i_start, i_stop", frame_ranges) + def test_chunked_read_matches_direct_read(self, tmp_path, monkeypatch, channel_indexes, i_start, i_stop): + """Shrinking the read budget must not change the data that comes back.""" + reader, _ = self.make_reader(tmp_path) + kwargs = dict(block_index=0, seg_index=0, i_start=i_start, i_stop=i_stop, channel_indexes=channel_indexes) + + direct = reader.get_analogsignal_chunk(**kwargs) + monkeypatch.setattr(biocamrawio, "_MAX_READ_BYTES", self.small_budget) + chunked = reader.get_analogsignal_chunk(**kwargs) + + assert chunked.dtype == direct.dtype + np.testing.assert_array_equal(chunked, direct) + + def test_read_budget(self, tmp_path, monkeypatch, frames_per_read): + """A long request for few channels is split, so the whole recording is never resident. + + See https://github.com/SpikeInterface/spikeinterface/issues/3303 + """ + reader, expected_all = self.make_reader(tmp_path) + # Another testing hack here, to make sure multiple reads are done. + monkeypatch.setattr(biocamrawio, "_MAX_READ_BYTES", self.small_budget) + + sig = reader.get_analogsignal_chunk( + block_index=0, seg_index=0, i_start=0, i_stop=self.n_frames, channel_indexes=[2] + ) + + np.testing.assert_array_equal(sig, expected_all[:, [2]]) + largest = max(frames_per_read) + assert largest <= self.small_budget_frames, f"a single read covered {largest} frames" + assert len(frames_per_read) > 2, "the request was not split, so the chunked path is untested" + + def test_requesting_all_channels(self, tmp_path, monkeypatch, frames_per_read): + """Requesting every channel shouldn't cause multiple reads (nothing to be gained).""" + reader, _ = self.make_reader(tmp_path) + monkeypatch.setattr(biocamrawio, "_MAX_READ_BYTES", self.small_budget) + + sig = reader.get_analogsignal_chunk( + block_index=0, seg_index=0, i_start=0, i_stop=self.n_frames, channel_indexes=None + ) + + assert max(frames_per_read) == self.n_frames, "the full-channel request was split, which only adds a copy" + assert sig.base is not None, "the full-channel read should be a view on the raw read, not a copy" + + if __name__ == "__main__": unittest.main() From 86a61e883101834ed372bab0d3bb238c76d68597 Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Fri, 31 Jul 2026 12:17:13 -0600 Subject: [PATCH 2/2] improve tests --- neo/test/rawiotest/test_biocamrawio.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/neo/test/rawiotest/test_biocamrawio.py b/neo/test/rawiotest/test_biocamrawio.py index db9e12ed2..7a59d4821 100644 --- a/neo/test/rawiotest/test_biocamrawio.py +++ b/neo/test/rawiotest/test_biocamrawio.py @@ -94,19 +94,6 @@ def test_chunked_read_matches_direct_read(self): assert chunked.dtype == direct.dtype np.testing.assert_array_equal(chunked, direct) - def test_all_channel_read_is_a_view(self): - """Asking for every channel returns the raw read itself rather than a copy of it. - - Copying there doubles the memory of a full read and costs a pass over the data, which is - the regression that pull request #1885 fixed. - """ - for entity in self.shape_variants: - with self.subTest(entity=entity): - reader = BiocamRawIO(filename=self.get_local_path(entity)) - reader.parse_header() - sig = reader.get_analogsignal_chunk(block_index=0, seg_index=0, channel_indexes=None) - assert sig.base is not None - def test_biocamrawio_gain(tmp_path): """Test that BiocamRawIO correctly reads the gain from a Biocam HDF5 file.