Fix a performance regression (and address #1532)

[kevindoran]

Previous commits (starting around 1814bc1) relating to #1532 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, which is a fix for the problem from #1532.

Edit: I previously got confused with a related issue in SpikeingInterface, #3303, which I think is connected.

[zm711]

Thanks for this @kevindoran. Seems like in trying to fix the original issue I may have kept it the same/made it worse. Appreciate the fix. Could we have you post some sort of benchmark in a comment about the improvement in performance (more like a terminal screenshot or copy of output from the terminal rather than just saying 25x)?

I would propose that you

1. keep the performance improvement in this PR
2. split off the performance regression test into a separate PR

I suggest this because we haven't discussed as maintainers if we want to support performance testing specifically and in general we try to avoid monkey patching as much as we can. So if you open that test as a separate PR that will give us an opportunity to discuss it as a concept for Neo without preventing us from merging this PR (so long as tests pass). How's that sound?

[kevindoran]

The test is a logical test; it doesn't measure execution time or memory usage. The monkey patch just adds a counter: it counts how many read calls there is, which is something we can deduce and assert logically, for example, it shouldn't be called multiple times when the whole stimulus is requested, but should be called multiple times when a small number of channels are requested.

If it was a performance test that measure execution time (and differ between machines), I agree that it wouldn't really belong in the codebase at the moment, as there isn't a dedicated machine they could run on.

[kevindoran]

Regarding the 25x. I was doing preprocessing and a spike sort that required 3 full reads of the recording. With the fix, subsequent reads hdf5 can use read cache in ram (I think). So exact speed-ups will depend on how many reads, and the size of the recording relative to available ram (and disk and cpu).

Phase_01.brw, 8 x 65536-frame chunks = 4.29 GB per pass
timeit -n 1 -r 3

script        time                          throughput   multiple
stock_warm    ~3.6 min (prime + 3 passes)    0.080 GB/s   1x
fixed_warm    ~8 s (prime + 3 passes)        2.03 GB/s    25.4x
stock_cold    ~2.7 min (3 passes)            0.079 GB/s   1x
fixed_cold    ~13 s (3 passes)               0.97 GB/s    12.3x

[zm711]

Thanks for the performance stats that's good enough for my purposes. We just like to have records associated with the PR so we can go back and try to assess things in the future. I think I would still prefer the monkey patch to be a separate proposal. But I'm willing to be overruled if @alejoe91 and @h-mayorquin are cool with it.

[alejoe91]

Suggested change

	if channel_indexes is None:
	channel_indexes = slice(None)
	n_ch = self._read_block(i_start, i_start + 1)[:, channel_indexes].shape[1]
	if channel_indexes is None:
	channel_indexes = slice(None)
	n_ch = self._num_channels
	else:
	n_ch = len(list(channel_indexes)

n_ch is available from the parse header, or is the length of channel indices. The only thing this logic doesn't cover is if we explicitly pass a slice as channel_indexes

[kevindoran]

Yea, if you don't pre-fetch, you need to handle all cases of slice (or don't support slice). Seems a shame to drop support for slice.

[kevindoran]

You will be definitely reading that [start, start+1] bit anyway, in the very next read, so I don't think there is any performance implications of doing the pre-fetch to check shape.

[h-mayorquin]

Unresolving here.

For the first point:

slice.indices(n) is stdlib and does exactly this job. One line covers every form.

For the second point, I don't understand, are we caching the response somehow?

[kevindoran]

Only remaining case is if a bool mask is passed.

[kevindoran]

For the second point, hdf5 has a read cache, and below that the OS cache.

[h-mayorquin]

Bool masks never reach _get_analogsignal_chunk. The base class converts them to integer indices with np.nonzero first, and raises if the length is wrong: baserawio.py#L988-L995. So what arrives is only None, a slice, or an integer array. If you want it robust to direct calls on the private method anyway, np.arange(self._num_channels)[channel_indexes].size covers every case including bool in one expression.

On the caching, the bytes are cached but that is not where the cost is:

• Raw is stored contiguously in every BioCam file (chunks=None on all three in the test data), so HDF5's chunk cache never applies. Only the OS page cache does, and that saves the disk read, not the h5py call, the allocation, or the copy out of the page cache.
• My numbers were already warm. Warmup call, then 11 to 15 reps, median. On the v3.x flat file the probe costs 0.101 ms against 0.101 ms for a full 1000-frame read. On the sparse file it costs a flat 1.72 ms against 2.81 ms for a 100-frame read.
• The sparse cost is Python-level decoding: reading both TOCs and walking every record in the covering block. No cache touches that.

Either way, the probe has to go for a reason that is not about speed. On i_start == i_stop == num_frames it reads one frame past the end, gets an empty array back, and _read_block reshapes it to (1, n_channels):

ValueError: cannot reshape array of size 0 into shape (1, 4096)

Master returns (0, n_ch) there. The real read handles the empty case on its own, since reshape(0, n_ch) is legal, so computing n_ch arithmetically fixes the crash and the overhead together.

[alejoe91]

@kevindoran what is the rationale behinf fixing this to 64MB?
It seems quite small...

[alejoe91]

Maybe this could be a class attribute and we could have a class method set_max_read_bytes to increase it/decrease it?

How is performance affected by it? The smaller this _MAX_READ_BYTES is, the lower the temporary RAM usage (final RAM usage is dictated by start/end frames and number of channels requested), but also the lower the speed. Could you test a couple of lager values to measure this?

[kevindoran]

I tried a sweep, and there was no improvement past 64 MB. I thought of adding as an attribute, but opted to just keep the change small and simple.

[h-mayorquin]

I agree with @alejoe91 that this should be an attribute rather than a module-level constant. Beyond giving users an escape hatch, the sweep below says one value cannot serve both layouts, so the knob ends up doing real work rather than being configurable for its own sake.

I generated a large synthetic file using the same approach as the writer in your tests, 4096 channels by 80000 frames, so a 625 MiB full all-channel read, and swept the budget for a whole-recording request of 3 channels. Same sweep against the event-based sparse file already in the test data.

Two things come out of it.

The budget is too large rather than too small. On the dense layout the minimum sits around 4 MiB at 102 ms, against 381 ms at 64 MiB, so the current value is about 3.7x off. The flat region above 64 MiB is real and matches what you found, but the improvement is below it, not above. Peak allocation is exactly budget + output at every point I measured, so a smaller budget is also a tighter guarantee: 4.5 MiB rather than 64.5 MiB on the same request.

The sparse reader wants the opposite, because its cost is per call rather than per byte. decode_event_based_raw_data re-reads TOC and EventsBasedSparseRawTOC and re-walks the event blob for the covering blocks on every call, so splitting multiplies a fixed cost: 1159 ms at 1 MiB against 183 ms unchunked. It does still need the ceiling, though. readHDF5t_brw4_sparse allocates np.zeros((nch, num_frames)) before it returns, so the sparse path has the same out-of-memory exposure as dense and the compression only helps on disk.

That suggests two defaults keyed on the read function:

class BiocamRawIO(BaseRawIO):
    # Ceiling on the temporary all-channel buffer. The dense layouts read a contiguous
    # slice, so cost scales with bytes and a small buffer stays in cache. The event-based
    # sparse decoder re-reads both tables of contents and re-walks the event blob on every
    # call, so its cost is per call rather than per byte and it wants a buffer large enough
    # that it almost never splits.
    max_read_bytes_dense = 4 * 1024 * 1024
    max_read_bytes_sparse = 512 * 1024 * 1024

budget = (
    self.max_read_bytes_sparse
    if self._read_function is readHDF5t_brw4_sparse
    else self.max_read_bytes_dense
)

[alejoe91]

I see where these tests are coming from, but in general we tend to only have tests on real data.

I would remove these tests and port the test_values and test_chunked_read_matches_direct_read using the existing test files: https://github.com/kevindoran/python-neo/blob/989423851488cb91a5cc639d32ba50bc801d077b/neo/test/rawiotest/test_biocamrawio.py#L26-L27

Could you do that?

[kevindoran]

Following the spec, it's a real brw file. It's much easier to follow unit tests that are self contained and don't rely on the format of some non-human readable file. If the project maintainers wish a different style, go for it. I'll tap out at this point.

[h-mayorquin]

There are three reason for our preference for files. This can admit edge cases but:

• The files on the test library are a public good and can be used for other libraries. This is good and make our job easier.
• Encapsulation. The test carries the format's internals inline, so understanding a failure means reading the writer too. I'd rather reason about "a file with these characteristics fails" than learn a second implementation before I can read the test.
• Independence. A file built in the test code shares the reader's assumptions, so it can't contradict them.

If you want to tap out I am fine with deferring this. I can add the file with these characteristics to gin, do a PR with the requested modifications and credit you.

[alejoe91]

Thanks @kevindoran

I think this is the right way to go!! I just added a couple of comments columns that I think will make the implementaion more performant and easier to maintain (test-wise)