perf(io): read gapped multi-dimensional selections in one call by ehennestad · Pull Request #923 · NeurodataWithoutBorders/matnwb

ehennestad · 2026-10-02T11:19:57Z

Motivation

Background — A DataStub read with two or more subscripts, such as stub(:, columns), goes through io.space.findShapes, which turns the indices of each dimension into hyperslab shapes, and io.space.getReadSpace, which ORs one hyperslab per shape into the file selection for a single H5D.read. Region references resolve through the same two functions. The rows of a ragged DynamicTable column are always a union of index ranges: row r spans elements index(r-1)+1 through index(r).

Problem — A user selects a subset of a multi-dimensional dataset with gaps between the selected indices, for example the waveforms of every other unit, or an irregular subset of channels or samples of a large ElectricalSeries. The read takes seconds to minutes while the same read of a contiguous range takes milliseconds. In the first example below, the 11 samples around each of 2,500 spikes, 2,296 runs of a 16-channel LFP series, take 41 s to read in one call, while reading each window in a call of its own takes 2.1 s. The cause is findShapes: it searches for the longest regularly strided block, removes it and repeats, so a selection of R irregular runs costs R passes and the time grows faster than linearly. A single call with 10,000 runs takes 113 s. #909 works around this in getRow by reading each run of a multi-dimensional column in a call of its own, which costs one file open per run: the 2,490 correct trials of a 5,000-trial table, 1,213 runs, take 1.25 s to get.

Solution — findShapes now splits the sorted indices into runs of consecutive indices with diff, in linear time, and only searches for strided blocks when most indices are isolated, as in 1:2:N. The same windows read in 0.28 s, 10,000 runs are split in 0.06 s, and getRow reads the elements of all requested rows of a column in one call, so the correct trials of a 5,000-trial table in the second example take 0.30 s instead of 1.25 s. HDF5 itself merges each hyperslab into the selection built so far, which costs 1.4 s at 10,000 hyperslabs and 166 s at 100,000, so a selection with more than 1,000 hyperslabs is read in groups along the dimension with the most runs and joined in memory.

What changed

io.space.findShapes emits one Block with step 1 per run of consecutive indices and one Point per isolated index, in increasing order. A pure stride such as 1:2:N still becomes one strided Block: when more than half of the indices are isolated, the strided block search runs first, and each pass is accepted only when its block covers a quarter of the remaining indices. That bounds the passes by a logarithm of the number of indices, and irregular isolated indices fall back to the run split after one pass.
HDF5LazyArray.load_mat_style reads a selection with more than MaxHyperslabsPerRead (1,000) hyperslabs in groups of whole runs along the dimension with the most shapes, each group with its own file selection and H5D.read, and joins the groups along that dimension. Memory use stays that of the selected elements.
The io.space.shape.Block constructor validates with an arguments block instead of inputParser, which brings construction from 45 µs to 12 µs per block. Thousands of blocks are built for a selection with thousands of runs.
getRow reads the elements of all requested rows of a file-backed ragged column in one call for every rank. The per-run loop that fix(dynamictable): read ragged columns once per index level in getRow #909 added for multi-dimensional columns is removed.
A DataStub read now closes the copied file dataspace it reads through. It was left open before.

Implementation notes

findShapes keeps findOptimalBlock as the strided search. The loop around it runs only while the remaining indices are mostly isolated (numel(runs) > numel(indices)/2) and the best block covers at least a quarter of them (minimumStrideCoverage). The remainder is split into runs by findRuns and createRunShapes. [1:2:99, 4] becomes Block(1:2:99) and Point(4); 1,000 random points in 1:50000 become one shape per run with no strided block.
readInGroups in load_mat_style splits the sorted unique indices of the split dimension at run boundaries into groups of floor(MaxHyperslabsPerRead / hyperslabs per shape of that dimension) runs. Each group is read and reshaped like a whole selection with readSelection and reshapeLoadedData, so the other dimensions come out in the requested order, and the groups are joined with cat. The split dimension is reordered afterwards when its subscript is unsorted or has repeats. The groups never recurse into readInGroups.
MaxHyperslabsPerRead is 1,000 because the time to OR N hyperslabs is about 6 µs × N + 13 ns × N² on this machine: at 1,000 hyperslabs the quadratic part is a third of the linear part. Each group is read from the dataset that load_mat_style already has open, so a group adds one H5D.read call and no file or dataset open. Grouping is linear in the number of runs (see Performance).
readDataRows in getRow reads the sorted unique elements of all rows with one readRows call and takes each row from that block by position. A column whose rows are all empty is not read; readEmptyRow builds the empty rows.
load_mat_style opens the file and dataset once and closes them after its last read. It is split into readSelection (read from the open dataset, convert types, close the dataspaces) and reshapeLoadedData (reshape and reorder) so the grouped path reuses them. The reorder step is unchanged, including its handling of ':' subscripts.

Related pull requests

fix(io): keep the requested order of unsorted DataStub subscripts #924 (fix-datastub-unsorted-subscript-order) is merged, and this PR is rebased on it. Its fixes to the reorder step, the ':' stand-in loop and isSelectionNormal(i), sit inside reshapeLoadedData, and its reordered* tests in HDF5LazyArrayTest pass here.
perf(io): open the HDF5 file once per DataStub read #910 (hdf5-lazy-array-open-once-per-read) is merged, and this PR is rebased on it. perf(io): open the HDF5 file once per DataStub read #910 opens the file and dataset once at the top of load_mat_style and reads a contiguous range of a 1-D dataset as one hyperslab. This PR keeps both: readSelection, readInGroups and getReadSpace use the dataset that load_mat_style opened, so a read in groups opens the file once. The tests of both PRs are in HDF5LazyArrayTest.
perf(io): read several DataStub selections in one pass #912 (read-selections-in-one-pass) is stacked on this PR, in stack #925: fix(dynamictable): read ragged columns once per index level in getRow #909 → perf(io): read gapped multi-dimensional selections in one call #923 → perf(io): read several DataStub selections in one pass #912. It adds loadSelections to LazyArray, DataStub and DataPipe. Its getRow change, which batched the per-run reads of fix(dynamictable): read ragged columns once per index level in getRow #909 through loadSelections, is dropped because this PR reads all rows in one call; getRow.m is identical in perf(io): read several DataStub selections in one pass #912 and this PR.

Examples

Spike-triggered windows of an LFP series

The snippet writes a 16-channel, 5-minute LFP series at 1 kHz, reads it back and selects the 11 samples around each of 2,500 spikes, once as a single DataStub read and once with one read per spike.

rng(923)
numChannels = 16;
numSamples = 300000; % 5 minutes at 1 kHz
lfp = types.core.TimeSeries('description', 'LFP, one row per channel', ...
    'data', rand(numChannels, numSamples), 'data_unit', 'volts', ...
    'starting_time', 0, 'starting_time_rate', 1000);
nwb = NwbFile('session_description', 'spike-triggered LFP', 'identifier', 'sta', ...
    'session_start_time', datetime(2018, 4, 25, 'TimeZone', 'local'));
nwb.acquisition.set('lfp', lfp);
nwbExport(nwb, 'sta.nwb');
nwbIn = nwbRead('sta.nwb', 'ignorecache');
stub = nwbIn.acquisition.get('lfp').data;

% 2,500 spikes and the 11 samples around each one.
spikeSamples = sort(randperm(numSamples - 20, 2500)) + 10;
windowSamples = unique(spikeSamples + (-5:5)')';

tic; windows = stub(:, windowSamples); oneCallSeconds = toc;
tic
perSpike = cell(1, numel(spikeSamples));
for iSpike = 1:numel(spikeSamples)
    perSpike{iSpike} = stub(:, spikeSamples(iSpike) + (-5:5));
end
perSpikeSeconds = toc;
fprintf('samples selected: %d in %d runs\n', numel(windowSamples), sum(diff(windowSamples) > 1) + 1);
fprintf('one call: %.3f s\n', oneCallSeconds);
fprintf('one call per spike: %.3f s\n', perSpikeSeconds);
fprintf('equal: %d\n', isequal(windows, lfp.data(:, windowSamples)));

Before — One call takes 41 s, twenty times longer than 2,500 separate reads.

samples selected: 26478 in 2296 runs
one call: 41.185 s
one call per spike: 2.117 s
equal: 1

After — One call takes 0.28 s and returns the same data.

samples selected: 26478 in 2296 runs
one call: 0.282 s
one call per spike: 1.706 s
equal: 1

The correct trials of a trials table

The snippet writes a 5,000-trial table with a correct column and a ragged lick_positions column holding the x and y of 1 to 20 licks per trial, reads it back and gets the rows of the correct trials.

rng(923)
numTrials = 5000;
licksPerTrial = randi([1 20], 1, numTrials);
lickPositions = types.hdmf_common.VectorData('description', 'x and y of each lick', ...
    'data', rand(2, sum(licksPerTrial)));
lickPositionsIndex = types.hdmf_common.VectorIndex('description', 'index into lick_positions, one row per trial', ...
    'data', uint64(cumsum(licksPerTrial))', 'target', types.untyped.ObjectView(lickPositions));
trialStarts = 10*(0:numTrials-1)';
nwb = NwbFile('session_description', 'trials with lick positions', 'identifier', 'trials', ...
    'session_start_time', datetime(2018, 4, 25, 'TimeZone', 'local'));
nwb.intervals_trials = types.core.TimeIntervals('description', 'trials', ...
    'colnames', {'start_time', 'stop_time', 'correct', 'lick_positions'}, ...
    'start_time', types.hdmf_common.VectorData('description', 'trial start', 'data', trialStarts), ...
    'stop_time', types.hdmf_common.VectorData('description', 'trial end', 'data', trialStarts + 5), ...
    'correct', types.hdmf_common.VectorData('description', 'whether the trial was correct', 'data', rand(numTrials, 1) < 0.5), ...
    'lick_positions', lickPositions, 'lick_positions_index', lickPositionsIndex, ...
    'id', types.hdmf_common.ElementIdentifiers('data', int64(0:numTrials-1)'));
nwbExport(nwb, 'trials.nwb');
nwbIn = nwbRead('trials.nwb', 'ignorecache');
trials = nwbIn.intervals_trials;

correctTrials = find(trials.correct.data.load());
tic; rows = trials.getRow(correctTrials); seconds = toc;
fprintf('correct trials: %d in %d runs\n', numel(correctTrials), sum(diff(correctTrials) > 1) + 1);
fprintf('getRow: %.3f s\n', seconds);
fprintf('equal: %d\n', isequal(rows, nwb.intervals_trials.getRow(correctTrials)));

Before — With #909, the correct trials are read one run at a time.

correct trials: 2490 in 1213 runs
getRow: 1.249 s
equal: 1

After — The rows are read in one call per column.

correct trials: 2490 in 1213 runs
getRow: 0.298 s
equal: 1

Performance

All timings are medians of three runs after a warm-up call, on one machine: MATLAB 26.1.0.3203278 (R2026a), HDF5 1.14.4, macOS. Before is the #909 branch this PR builds on.

findShapes alone. Selections of R runs with random lengths 1 to 8 and gaps 1 to 20.

R	Before	After
100	0.018 s	0.002 s
1,000	1.21 s	0.008 s
10,000	113 s (single call)	0.063 s

DataStub reads. The gapped selection has 2,500 runs of 1 to 8 columns with gaps of 1 to 20 between them (11,338 elements, span 5 to 37,189). The span read is stub(:, 5:37189).

Read	`8 x 50000` before	`8 x 50000` after	`64 x 211722` before	`64 x 211722` after
gapped, one call	7.46 s	0.127 s	7.57 s	0.113 s
gapped, one call per run	1.55 s	1.32 s	1.62 s	1.36 s
span	0.003 s	0.002 s	0.006 s	0.006 s
pure stride `1:2:end`	0.002 s	0.006 s	0.020 s	0.022 s
1,000 random points	0.121 s	0.035 s	0.168 s	0.034 s

getRow end to end. The Units tables are the #787 script from #909 with 20 units (26,583 spikes) and the full size (172 units, 211,722 spikes). The 5,000-row table has an 8 x n ragged column with 1 to 20 elements per row.

Call	Before	After
20 units, `toTable` after `nwbRead`	0.184 s	0.129 s
20 units, `getRow(1:2:20)`	0.111 s	0.069 s
20 units, `getRow([1 20])`	0.036 s	0.027 s
172 units, `toTable` after `nwbRead`	1.17 s	0.75 s
172 units, `getRow(1:2:172)`	0.669 s	0.360 s
172 units, `getRow([1 172])`	0.024 s	0.016 s
5,000 rows, `getRow(1:2:5000)`	2.01 s	0.153 s
5,000 rows, `toTable`	0.059 s	0.044 s

Where HDF5 becomes the limit. H5S.select_hyperslab with H5S_SELECT_OR for N blocks of 4 columns on a 2-D space, one call per block, single run.

N	Time
1,000	0.028 s
10,000	1.43 s
100,000	166 s

Reading in groups of at most 1,000 hyperslabs keeps a one-call gapped read linear in the number of runs (8 x n dataset, one call, median of three):

Runs	Elements	One call
2,500	11,372	0.120 s
10,000	44,803	0.447 s
50,000	224,487	2.18 s

Without grouping, 50,000 runs would be one selection of 50,000 hyperslabs, which the table above puts between 1.4 s and 166 s to build.

How to test

Run the two snippets above on the base branch and here and compare the one call and getRow lines. Then run the tests. The SpaceTest cases for runs, strides and many runs describe the new findShapes output; testManyRunsSelectExactlyTheInput takes minutes on the base branch. selectionWithManyRunsIsReadInGroups in HDF5LazyArrayTest reads 2,001 runs of a 3-D dataset, more than one selection holds, with sorted, reversed and repeated indices. testGetRowWithGapsReadsEachLevelOnce in DynamicTableRaggedReadTest fails on the base branch, where the waveforms of units 1 and 4 are read in two calls.

results = runtests({'tests.unit.io.SpaceTest', 'tests.unit.io.backend.HDF5LazyArrayTest', ...
    'tests.unit.DynamicTableRaggedReadTest', 'tests.unit.dataStubTest'});
disp(table(results))

Checklist

Have you ensured the PR description clearly describes the problem and solutions?

🤖 Generated with Claude Code

codecov · 2026-10-02T11:35:42Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 95.28%. Comparing base (8aac84e) to head (9ed6562).

Additional details and impacted files

@@                           Coverage Diff                           @@
##           fix-get-row-ragged-read-performance     #923      +/-   ##
=======================================================================
+ Coverage                                95.25%   95.28%   +0.02%     
=======================================================================
  Files                                      239      239              
  Lines                                     8880     8915      +35     
=======================================================================
+ Hits                                      8459     8495      +36     
+ Misses                                     421      420       -1

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

io.space.findShapes searched for the longest regularly strided block, removed it and repeated, so a selection of R irregular runs cost R passes and the time grew faster than linearly: a single call with 10,000 runs took 113 s. Every multi-subscript DataStub read and every region reference goes through it, so a gapped selection of a multi-dimensional dataset took seconds while the covering span took milliseconds. findShapes now splits the sorted indices into runs of consecutive indices with diff, one Block per run and one Point per isolated index. The strided search runs only when most indices are isolated, and each pass must cover a quarter of the remaining indices, so 1:2:N is still one strided Block while irregular points fall back to the run split after one pass. HDF5 merges each hyperslab into the selection built so far, which costs 1.4 s at 10,000 hyperslabs and 166 s at 100,000, so load_mat_style reads a selection with more than 1,000 hyperslabs in groups of whole runs along the dimension with the most shapes and joins them in memory. The Block constructor uses an arguments block instead of inputParser, which matters when thousands of blocks are built. getRow reads the elements of all requested rows of a file-backed ragged column in one call for every rank; the per-run loop for multi-dimensional columns is no longer needed. Every other row of a 5,000-row table with an 8 x n column takes 0.15 s instead of 2.0 s. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

ehennestad force-pushed the perf-find-shapes-contiguous-runs branch from efb0019 to 9ed6562 Compare October 2, 2026 12:34

ehennestad marked this pull request as draft October 2, 2026 13:22

ehennestad mentioned this pull request Oct 2, 2026

fix(io): keep the requested order of unsorted DataStub subscripts #924

Merged

1 task done

ehennestad added this pull request to stack #925 October 2, 2026 13:35

ehennestad mentioned this pull request Oct 2, 2026

perf(io): read several DataStub selections in one pass #912

Draft

2 of 3 tasks

ehennestad force-pushed the perf-find-shapes-contiguous-runs branch from 9ed6562 to 9bc4495 Compare October 6, 2026 05:38

ehennestad mentioned this pull request Oct 6, 2026

perf(io): open the HDF5 file once per DataStub read #910

Merged

2 of 3 tasks

ehennestad force-pushed the perf-find-shapes-contiguous-runs branch from 9bc4495 to c5329aa Compare October 6, 2026 05:47

ehennestad force-pushed the perf-find-shapes-contiguous-runs branch from c5329aa to 7c4f63b Compare October 6, 2026 18:09

ehennestad mentioned this pull request Oct 6, 2026

refactor(dynamictable): rename getRow locals and use an arguments block #940

Open

2 of 3 tasks

ehennestad removed this pull request from stack #925 October 6, 2026 20:44

ehennestad force-pushed the perf-find-shapes-contiguous-runs branch from 7c4f63b to a2ca013 Compare October 6, 2026 20:44

ehennestad added this pull request to stack #941 October 6, 2026 20:45

ehennestad added this to the v2.12.0 milestone Oct 7, 2026

Navigation Menu

Repository navigation

perf(io): read gapped multi-dimensional selections in one call - #923

perf(io): read gapped multi-dimensional selections in one call#923
ehennestad wants to merge 1 commit into
fix-get-row-ragged-read-performancefrom
perf-find-shapes-contiguous-runs

ehennestad commented Oct 2, 2026 •

edited

Loading

Uh oh!

codecov Bot commented Oct 2, 2026 •

edited

Loading

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

ehennestad commented Oct 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Motivation

What changed

Related pull requests

Examples

Spike-triggered windows of an LFP series

The correct trials of a trials table

Performance

How to test

Checklist

Uh oh!

codecov Bot commented Oct 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

ehennestad commented Oct 2, 2026 •

edited

Loading

codecov Bot commented Oct 2, 2026 •

edited

Loading