diff --git a/+tests/+unit/+io/+backend/+doubles/HDF5LazyArraySpy.m b/+tests/+unit/+io/+backend/+doubles/HDF5LazyArraySpy.m new file mode 100644 index 000000000..a0bc7affa --- /dev/null +++ b/+tests/+unit/+io/+backend/+doubles/HDF5LazyArraySpy.m @@ -0,0 +1,46 @@ +classdef HDF5LazyArraySpy < io.backend.hdf5.HDF5LazyArray +% HDF5LazyArraySpy - Test spy that records reads from an HDF5 dataset. +% Delegates every read to HDF5LazyArray and records the calls. Pass it to +% the DataStub constructor to see how often, and for which elements, code +% under test reads the dataset. +% +% An empty load_mat_style selection is recorded twice, because +% HDF5LazyArray reads the first element through load_mat_style to get the +% data type. + + properties (SetAccess = private) + % LoadCount - Number of calls to load_mat_style and load_h5_style. + LoadCount (1,1) double = 0 + + % MatStyleSelections - Subscripts passed to each load_mat_style + % call, one cell array per call. + MatStyleSelections (1,:) cell = {} + end + + methods + function obj = HDF5LazyArraySpy(filename, path) + arguments + filename (1,1) string + path (1,1) string + end + obj@io.backend.hdf5.HDF5LazyArray(filename, path); + end + + function reset(obj) + % reset - Forget the calls recorded so far. + obj.LoadCount = 0; + obj.MatStyleSelections = {}; + end + + function data = load_mat_style(obj, varargin) + obj.LoadCount = obj.LoadCount + 1; + obj.MatStyleSelections{end+1} = varargin; + data = load_mat_style@io.backend.hdf5.HDF5LazyArray(obj, varargin{:}); + end + + function data = load_h5_style(obj, varargin) + obj.LoadCount = obj.LoadCount + 1; + data = load_h5_style@io.backend.hdf5.HDF5LazyArray(obj, varargin{:}); + end + end +end diff --git a/+tests/+unit/DynamicTableRaggedReadTest.m b/+tests/+unit/DynamicTableRaggedReadTest.m new file mode 100644 index 000000000..5c5a02f19 --- /dev/null +++ b/+tests/+unit/DynamicTableRaggedReadTest.m @@ -0,0 +1,346 @@ +classdef DynamicTableRaggedReadTest < tests.abstract.NwbTestCase +% DynamicTableRaggedReadTest - Tests for reading rows of ragged DynamicTable columns. +% +% The shared fixture is a Units table with four units and three ragged columns: +% waveforms doubly ragged. The units have 3, 0, 2 and 5 spikes, and each +% spike is a [4 x 2] block: 4 samples on 2 electrodes. +% labels ragged text. +% peaks ragged compound, stored as a table. + + properties (TestParameter) + RowIndices = struct( ... + 'firstRow', 1, ... + 'rowWithoutSpikes', 2, ... + 'allRows', 1:4, ... + 'rowsWithGap', [1 4], ... + 'unsortedRows', [4 2], ... + 'lastRow', 4) + end + + properties (Access = private) + FileName (1,1) string + MemoryUnits + FileUnits + end + + methods (TestClassSetup) + function exportUnitsTable(testCase) + testCase.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); + + nwb = tests.factory.NWBFile(); + nwb.units = testCase.createUnitsTable(); + testCase.FileName = "ragged_read.nwb"; + nwbExport(nwb, testCase.FileName); + + testCase.MemoryUnits = nwb.units; + nwbIn = nwbRead(testCase.FileName, 'ignorecache'); + testCase.FileUnits = nwbIn.units; + end + end + + methods (Test) + function testGetRowFromFileMatchesMemory(testCase, RowIndices) + testCase.verifyEqual( ... + testCase.FileUnits.getRow(RowIndices), ... + testCase.MemoryUnits.getRow(RowIndices)); + end + + function testToTableFromFileMatchesMemory(testCase) + testCase.verifyEqual( ... + testCase.FileUnits.toTable(), ... + testCase.MemoryUnits.toTable()); + end + + function testGetRowFromBoundDataPipeMatchesUnbound(testCase, RowIndices) + nwb = tests.factory.NWBFile(); + nwb.units = testCase.createUnitsTable(UseDataPipes=true); + expectedRows = nwb.units.getRow(RowIndices, 'columns', {'waveforms'}); + + % Exporting binds each DataPipe to its dataset in the file. + nwbExport(nwb, testCase.getRandomFilename()); + testCase.assertTrue(nwb.units.waveforms.data.isBound) + + testCase.verifyEqual( ... + nwb.units.getRow(RowIndices, 'columns', {'waveforms'}), ... + expectedRows); + end + + function testEmptyRowDoesNotDependOnOtherRequestedRows(testCase) + % The second unit has no spike times. + [spikeTimes, spikeTimesIndex] = util.create_indexed_column( ... + {[1 2 3], zeros(1, 0), [4 5]}); + nwb = tests.factory.NWBFile(); + nwb.units = types.core.Units( ... + 'description', 'units with an empty row', ... + 'colnames', {'spike_times'}, ... + 'id', types.hdmf_common.ElementIdentifiers('data', int64(0:2)'), ... + 'spike_times', spikeTimes, ... + 'spike_times_index', spikeTimesIndex); + fileName = testCase.getRandomFilename(); + nwbExport(nwb, fileName); + nwbIn = nwbRead(fileName, 'ignorecache'); + + testCase.verifyEmptyRowIsSameAloneAndWithOtherRows(nwb.units) + testCase.verifyEmptyRowIsSameAloneAndWithOtherRows(nwbIn.units) + end + + function testGetRowWhenColumnHasNoElements(testCase) + % Every row of the column is empty, so its dataset has no + % elements. An extendable dataset with nothing appended is + % written that way. + amplitudes = types.hdmf_common.VectorData( ... + 'description', 'amplitudes', ... + 'data', types.untyped.DataPipe('maxSize', Inf, 'dataType', 'double')); + amplitudesIndex = types.hdmf_common.VectorIndex( ... + 'description', 'index into amplitudes', ... + 'data', uint64([0; 0]), ... + 'target', types.untyped.ObjectView(amplitudes)); + nwb = tests.factory.NWBFile(); + nwb.acquisition.set('rows_without_amplitudes', types.hdmf_common.DynamicTable( ... + 'description', 'table whose rows have no amplitudes', ... + 'colnames', {'amplitudes'}, ... + 'amplitudes', amplitudes, ... + 'amplitudes_index', amplitudesIndex, ... + 'id', types.hdmf_common.ElementIdentifiers('data', int64([0; 1])))); + fileName = testCase.getRandomFilename(); + nwbExport(nwb, fileName); + nwbIn = nwbRead(fileName, 'ignorecache'); + tableIn = nwbIn.acquisition.get('rows_without_amplitudes'); + + vector = tableIn.vectordata.get('amplitudes'); + testCase.assertClass(vector.data, 'types.untyped.DataPipe') + testCase.verifyRowsHaveNoAmplitudes(tableIn) + + vector.data = types.untyped.DataStub(fileName, '/acquisition/rows_without_amplitudes/amplitudes'); + testCase.verifyRowsHaveNoAmplitudes(tableIn) + end + + function testGetRowWhenCompoundColumnHasNoElements(testCase) + % Every row of the compound column is empty, so its dataset has + % no elements. Export writes a compound dataset only with rows, + % so the rows are dropped from the file afterwards, which leaves + % the dataset as another NWB writer would write it. + peaks = types.hdmf_common.VectorData( ... + 'description', 'peaks', ... + 'data', table(uint32(1), {'first'}, true, ... + 'VariableNames', {'index', 'label', 'flag'})); + peaksIndex = types.hdmf_common.VectorIndex( ... + 'description', 'index into peaks', ... + 'data', uint64([1; 1]), ... + 'target', types.untyped.ObjectView(peaks)); + nwb = tests.factory.NWBFile(); + nwb.acquisition.set('rows_without_peaks', types.hdmf_common.DynamicTable( ... + 'description', 'table whose rows have no peaks', ... + 'colnames', {'peaks'}, ... + 'peaks', peaks, ... + 'peaks_index', peaksIndex, ... + 'id', types.hdmf_common.ElementIdentifiers('data', int64([0; 1])))); + fileName = testCase.getRandomFilename(); + nwbExport(nwb, fileName); + dropDatasetRows(fileName, '/acquisition/rows_without_peaks/peaks'); + h5write(fileName, '/acquisition/rows_without_peaks/peaks_index', uint64([0; 0])); + nwbIn = nwbRead(fileName, 'ignorecache'); + tableIn = nwbIn.acquisition.get('rows_without_peaks'); + + rows = tableIn.getRow(1:2); + + % Each empty row keeps the members and types of the compound. + expectedRow = table(uint32.empty(0, 1), cell(0, 1), false(0, 1), ... + 'VariableNames', {'index', 'label', 'flag'}); + testCase.verifyEqual(rows.peaks, {expectedRow; expectedRow}); + testCase.verifyEqual(height(tableIn.toTable()), 2); + end + + function testGetRowOfCompoundColumnStoredAsStructOfColumns(testCase) + % A compound column built in memory may hold a scalar struct whose + % fields are the columns of the compound type, instead of a table. + peaks = types.hdmf_common.VectorData( ... + 'description', 'peaks as a struct of columns', ... + 'data', struct('amplitude', (1:6)', 'time', (11:16)')); + peaksIndex = types.hdmf_common.VectorIndex( ... + 'description', 'index into peaks', ... + 'data', uint64([1; 3; 4; 6]), ... + 'target', types.untyped.ObjectView(peaks)); + dynamicTable = types.hdmf_common.DynamicTable( ... + 'description', 'table with a ragged compound column', ... + 'colnames', {'peaks'}, ... + 'peaks', peaks, ... + 'peaks_index', peaksIndex, ... + 'id', types.hdmf_common.ElementIdentifiers('data', int64(0:3)')); + + rows = dynamicTable.getRow([1 4]); + + expectedPeaks = { ... + struct('amplitude', 1, 'time', 11); ... + struct('amplitude', (5:6)', 'time', (15:16)')}; + testCase.verifyEqual(rows.peaks, expectedPeaks); + end + + function testGetRowReadsEachLevelOnce(testCase) + [units, spies] = testCase.readUnitsWithSpies(); + + units.getRow(1:4, 'columns', {'waveforms'}); + + for columnName = string(fieldnames(spies))' + testCase.verifyEqual(spies.(columnName).LoadCount, 1, ... + sprintf('Expected one read of "%s".', columnName)); + end + end + + function testGetRowReadsOnlyRequestedElements(testCase) + % Units 1 and 4 hold spikes 1-3 and 6-10, whose waveforms are + % columns 1-6 and 11-20 of the data. Unit 3 lies between them. + [units, spies] = testCase.readUnitsWithSpies(); + + units.getRow([1 4], 'columns', {'waveforms'}); + + testCase.verifyEqual(getReadElements(spies.waveforms_index_index), [1, 3, 4]); + testCase.verifyEqual(getReadElements(spies.waveforms_index), [1:3, 5:10]); + testCase.verifyEqual(getReadElements(spies.waveforms), [1:6, 11:20]); + end + end + + methods (Access = private) + function [units, spies] = readUnitsWithSpies(testCase) + % readUnitsWithSpies - Read the units table with a spy on each dataset of the waveforms column. + nwbIn = nwbRead(testCase.FileName, 'ignorecache'); + units = nwbIn.units; + spies = struct(); + for columnName = ["waveforms", "waveforms_index", "waveforms_index_index"] + vector = units.(columnName); + spy = tests.unit.io.backend.doubles.HDF5LazyArraySpy( ... + vector.data.filename, vector.data.path); + vector.data = types.untyped.DataStub( ... + vector.data.filename, vector.data.path, [], [], spy); + % Assigning data runs validation, which also reads the dataset. + spy.reset() + spies.(columnName) = spy; + end + end + + function verifyRowsHaveNoAmplitudes(testCase, dynamicTable) + rows = dynamicTable.getRow(1:2); + testCase.verifyEqual(rows.amplitudes, {zeros(0, 0); zeros(0, 0)}); + testCase.verifyEqual(height(dynamicTable.toTable()), 2); + end + + function verifyEmptyRowIsSameAloneAndWithOtherRows(testCase, units) + emptyRowAlone = units.getRow(2).spike_times{1}; + allRows = units.getRow(1:3); + testCase.verifyEqual(allRows.spike_times{2}, emptyRowAlone); + end + end + + methods (Static, Access = private) + function units = createUnitsTable(options) + arguments + options.UseDataPipes (1,1) logical = false + end + + spikesPerUnit = [3 0 2 5]; + numSamples = 4; + numElectrodes = 2; + numSpikes = sum(spikesPerUnit); + + % Each column of waveformData is the waveform of one spike on one electrode. + waveformData = reshape(1:(numSamples*numElectrodes*numSpikes), numSamples, []); + waveformsIndexData = uint64(numElectrodes:numElectrodes:(numElectrodes*numSpikes))'; + waveformsIndexIndexData = uint64(cumsum(spikesPerUnit))'; + if options.UseDataPipes + waveformData = types.untyped.DataPipe( ... + 'data', waveformData, 'maxSize', [numSamples, Inf], 'axis', 2); + waveformsIndexData = types.untyped.DataPipe( ... + 'data', waveformsIndexData, 'maxSize', Inf); + waveformsIndexIndexData = types.untyped.DataPipe( ... + 'data', waveformsIndexIndexData, 'maxSize', Inf); + end + + waveforms = types.hdmf_common.VectorData( ... + 'description', 'waveforms', ... + 'data', waveformData); + waveformsIndex = types.hdmf_common.VectorIndex( ... + 'description', 'index into waveforms, one row per spike', ... + 'data', waveformsIndexData, ... + 'target', types.untyped.ObjectView(waveforms)); + waveformsIndexIndex = types.hdmf_common.VectorIndex( ... + 'description', 'index into waveforms_index, one row per unit', ... + 'data', waveformsIndexIndexData, ... + 'target', types.untyped.ObjectView(waveformsIndex)); + + units = types.core.Units( ... + 'description', 'units with ragged columns', ... + 'colnames', {'waveforms'}, ... + 'id', types.hdmf_common.ElementIdentifiers('data', int64(0:3)'), ... + 'waveforms', waveforms, ... + 'waveforms_index', waveformsIndex, ... + 'waveforms_index_index', waveformsIndexIndex); + + labels = types.hdmf_common.VectorData( ... + 'description', 'labels', ... + 'data', {'a'; 'b'; 'c'; 'd'; 'e'; 'f'; 'g'}); + labelsIndex = types.hdmf_common.VectorIndex( ... + 'description', 'index into labels', ... + 'data', uint64([2; 3; 4; 7]), ... + 'target', types.untyped.ObjectView(labels)); + units.addColumn('labels', labels, 'labels_index', labelsIndex); + + peaks = types.hdmf_common.VectorData( ... + 'description', 'peaks', ... + 'data', table((1:6)', (11:16)', 'VariableNames', {'amplitude', 'time'})); + peaksIndex = types.hdmf_common.VectorIndex( ... + 'description', 'index into peaks', ... + 'data', uint64([1; 3; 4; 6]), ... + 'target', types.untyped.ObjectView(peaks)); + units.addColumn('peaks', peaks, 'peaks_index', peaksIndex); + end + end +end + +function elements = getReadElements(spy) +% getReadElements - Sorted elements that a spy's load_mat_style calls selected in their last subscript. +lastSubscripts = cellfun(@(selection) reshape(selection{end}, 1, []), ... + spy.MatStyleSelections, 'UniformOutput', false); +elements = sort([lastSubscripts{:}]); +end + +function dropDatasetRows(fileName, datasetPath) +% dropDatasetRows - Replace a dataset with one of the same type and attributes that has no rows. +attributeNames = {h5info(fileName, datasetPath).Attributes.Name}; + +fileId = H5F.open(fileName, 'H5F_ACC_RDWR', 'H5P_DEFAULT'); +fileCleanup = onCleanup(@() H5F.close(fileId)); %#ok + +datasetId = H5D.open(fileId, datasetPath); +typeId = H5D.get_type(datasetId); +typeCleanup = onCleanup(@() H5T.close(typeId)); %#ok +attributes = cellfun(@(name) readAttribute(datasetId, name), attributeNames, ... + 'UniformOutput', false); +H5D.close(datasetId); + +H5L.delete(fileId, datasetPath, 'H5P_DEFAULT'); + +spaceId = H5S.create_simple(1, 0, []); +spaceCleanup = onCleanup(@() H5S.close(spaceId)); %#ok +datasetId = H5D.create(fileId, datasetPath, typeId, spaceId, 'H5P_DEFAULT'); +% The attributes identify the dataset as a VectorData when the file is read. +cellfun(@(attribute) writeAttribute(datasetId, attribute), attributes); +H5D.close(datasetId); +end + +function attribute = readAttribute(objectId, name) +attributeId = H5A.open(objectId, name); +attribute = struct('name', name, ... + 'typeId', H5A.get_type(attributeId), ... + 'spaceId', H5A.get_space(attributeId), ... + 'value', {H5A.read(attributeId)}); +H5A.close(attributeId); +end + +function writeAttribute(objectId, attribute) +attributeId = H5A.create(objectId, attribute.name, attribute.typeId, attribute.spaceId, ... + 'H5P_DEFAULT'); +H5A.write(attributeId, attribute.typeId, attribute.value); +H5A.close(attributeId); +H5S.close(attribute.spaceId); +H5T.close(attribute.typeId); +end diff --git a/+types/+util/+dynamictable/getRow.m b/+types/+util/+dynamictable/getRow.m index 0ab8bad74..9af6fa97a 100644 --- a/+types/+util/+dynamictable/getRow.m +++ b/+types/+util/+dynamictable/getRow.m @@ -152,105 +152,341 @@ function mustBeCellstrOrEmpty(value) % VectorIndex columns after it, the outermost index last. A column without % an index returns its rows as an array, or as a struct or table with one % array per member for a compound column. A ragged column returns a cell -% array with one cell per row, nested once per index level: the outermost -% index gives the element range of each row, and the function recurses into -% the chain below it for the elements of that range. -column = vectorChainNames{end}; -if isprop(DynamicTable, column) - Vector = DynamicTable.(column); -elseif isprop(DynamicTable, 'vectorindex') && DynamicTable.vectorindex.isKey(column) % Schema version < 2.3.0 - Vector = DynamicTable.vectorindex.get(column); +% array with one cell per row, nested once per index level. +vectorChain = cell(size(vectorChainNames)); +for iVector = 1:numel(vectorChainNames) + vectorChain{iVector} = getColumnVector(DynamicTable, vectorChainNames{iVector}); +end + +if isscalar(vectorChain) + columnRows = orientRows(vectorChain{1}, readRows(vectorChain{1}, rowIndices)); else - Vector = DynamicTable.vectordata.get(column); + columnRows = getRaggedRows(vectorChain, rowIndices); +end end -if isscalar(vectorChainNames) - if isa(Vector.data, 'types.untyped.DataStub') || ... - isa(Vector.data,'types.untyped.DataPipe') - if isa(Vector.data, 'types.untyped.DataStub') - dataSize = Vector.data.dims; - else - dataSize = Vector.data.internal.maxSize; - end - if length(dataSize) == 2 && dataSize(2) == 1 - % A column vector is indexed with one subscript. - numSubscripts = 1; - else - numSubscripts = length(dataSize); - end +function vector = getColumnVector(DynamicTable, columnName) +% getColumnVector - Get the VectorData or VectorIndex object of a column by name. +if isprop(DynamicTable, columnName) + vector = DynamicTable.(columnName); +elseif isprop(DynamicTable, 'vectorindex') && DynamicTable.vectorindex.isKey(columnName) % Schema version < 2.3.0 + vector = DynamicTable.vectorindex.get(columnName); +else + vector = DynamicTable.vectordata.get(columnName); +end +end + +function selected = getRaggedRows(vectorChain, rowIndices) +% getRaggedRows - Get rows of a ragged column, reading only the requested elements. +% +% vectorChain lists the column's VectorData first and its VectorIndex objects +% after it, the outermost index last. Each index level is read in one call, +% for exactly the index elements the requested rows need. Data in a file is +% read for exactly the elements of the requested rows, in one call or in one +% call per stretch of consecutive elements, and the rows are then taken from +% those reads in memory. Data in memory is indexed once per row (see +% readDataRows). +% +% selected is a cell array with one cell per requested row. For a doubly +% ragged column, each cell holds a cell array with one cell per element. + +% Row indices are shifted by subtraction, which saturates for unsigned integers. +rowIndices = double(reshape(rowIndices, 1, [])); + +numLevels = numel(vectorChain); +% rowStarts{iLevel} and rowStops{iLevel} give the element range, in the +% level below, of each row requested from level iLevel. +rowStarts = cell(1, numLevels); +rowStops = cell(1, numLevels); + +% Read the index levels from the outermost down. levelRows lists the rows +% requested from the current level in output order, and after the loop the +% data elements requested from the column. +levelRows = rowIndices; +for iLevel = numLevels:-1:2 + indexVector = vectorChain{iLevel}; + assert(isa(indexVector, 'types.hdmf_common.VectorIndex') || isa(indexVector, 'types.core.VectorIndex'), ... + 'NWB:DynamicTable:GetRow:InternalError', ... + 'Internal VectorIndex Stack is not using VectorIndex objects!'); + + [rowStarts{iLevel}, rowStops{iLevel}] = readElementRanges(indexVector, levelRows); + levelRows = expandRanges(rowStarts{iLevel}, rowStops{iLevel}); +end + +rowValues = readDataRows(vectorChain{1}, rowStarts{2}, rowStops{2}, levelRows); + +% Nest the rows of each level under their rows in the level above. An empty +% row has a stop below its start. +for iLevel = 3:numLevels + rowLengths = max(rowStops{iLevel} - rowStarts{iLevel} + 1, 0); + rowValues = mat2cell(rowValues, rowLengths(:), 1); +end +selected = rowValues; +end + +function [starts, stops] = readElementRanges(indexVector, rowIndices) +% readElementRanges - Element range of VectorIndex rows in the level below. +% +% Row r spans elements index(r-1)+1 through index(r) of the level below, +% where index(0) is 0. The index values of all rows are read in one call. +starts = ones(size(rowIndices)); +stops = zeros(size(rowIndices)); +if isempty(rowIndices) + return +end + +% Each row needs index(r) and index(r-1), so neighbouring rows share +% elements. The selection is read sorted and without duplicates, which +% keeps the HDF5 reader on its fast path. index(0) is not stored. +rowsToRead = unique([rowIndices, rowIndices - 1]); +rowsToRead(rowsToRead == 0) = []; +indexValues = readIndexValues(indexVector, rowsToRead); + +% Look up each row's stop, and its previous row's stop, in the single read. +% Row 1 has no previous row and keeps its default start of 1. +[~, stopPositions] = ismember(rowIndices, rowsToRead); +stops = reshape(indexValues(stopPositions), size(rowIndices)); +hasPreviousRow = rowIndices > 1; +[~, previousPositions] = ismember(rowIndices(hasPreviousRow) - 1, rowsToRead); +starts(hasPreviousRow) = indexValues(previousPositions) + 1; +end + +function elements = expandRanges(starts, stops) +% expandRanges - Concatenate the ranges starts(i):stops(i) into one row vector. +hasElements = starts <= stops; +starts = starts(hasElements); +stops = stops(hasElements); +if isempty(starts) + elements = zeros(1, 0); + return +end + +% Consecutive elements differ by 1 within a range. The first element of a +% range differs from the last element of the range before it by the gap. +steps = ones(1, sum(stops - starts + 1)); +rangeFirst = cumsum([1, stops(1:end-1) - starts(1:end-1) + 1]); +steps(rangeFirst) = [starts(1), starts(2:end) - stops(1:end-1)]; +elements = cumsum(steps); +end + +function rowValues = readDataRows(dataVector, starts, stops, elements) +% readDataRows - Get the data of ragged rows, one cell per row. +% +% Row i holds data elements starts(i) through stops(i), and elements lists +% the elements of all rows. +% +% A column in memory is indexed once per row, which copies each element +% once. A column in a file is read in as few calls as possible, and the rows +% are then taken from those reads in memory. A column indexed with one +% subscript reads any set of elements in one call. With more subscripts, a +% selection with gaps becomes one hyperslab per stretch of consecutive +% elements, and building it takes time that grows faster than linearly with +% the number of stretches, so each stretch is read in a call of its own. +% A run is the set of elements read by one call. +[numSubscripts, rowAxis] = getRowDimension(dataVector); +% A DataPipe is indexed through its subsref, which reads from the file once +% the pipe is bound to it and from the pipe's own data before that. Both are +% read in blocks, like a DataStub. +isFileBacked = isa(dataVector.data, 'types.untyped.DataStub') ... + || isa(dataVector.data, 'types.untyped.DataPipe'); + +if isFileBacked + elements = unique(elements); + numElements = numel(elements); + % runFirst(j) and runLast(j) are the positions in elements of the first + % and last element of run j. + if numElements == 0 + [runFirst, runLast] = deal(zeros(1, 0)); + elseif numSubscripts == 1 + [runFirst, runLast] = deal(1, numElements); else - if iscolumn(Vector.data) - % A column vector is indexed with one subscript. - numSubscripts = 1; - elseif istable(Vector.data) - % A compound column held as a table is indexed by row only. - numSubscripts = 1; - else - numSubscripts = ndims(Vector.data); - end + runFirst = find([true, diff(elements) > 1]); + runLast = [runFirst(2:end) - 1, numElements]; end - - subscripts = repmat({':'}, 1, numSubscripts); - if isa(Vector.data, 'types.untyped.DataPipe') - subscripts{Vector.data.axis} = rowIndices; - else - subscripts{end} = rowIndices; + + blocks = cell(size(runFirst)); + for iRun = 1:numel(runFirst) + blocks{iRun} = readRows(dataVector, elements(runFirst(iRun):runLast(iRun))); end - - if (isstruct(Vector.data) && isscalar(Vector.data)) || istable(Vector.data) - if istable(Vector.data) - columnRows = table(); - fields = Vector.data.Properties.VariableNames; - else - columnRows = struct(); - fields = fieldnames(Vector.data); - end - - for iField = 1:length(fields) - fieldName = fields{iField}; - memberData = Vector.data.(fieldName); - columnRows.(fieldName) = memberData(subscripts{:}); - end + % runOfPosition(k) is the run that holds elements(k). + runOfPosition = zeros(1, numElements); + runOfPosition(runFirst) = 1; + runOfPosition = cumsum(runOfPosition); + % The elements of a row are consecutive in elements and lie in one run. + [~, startPositions] = ismember(starts, elements); +end + +isEmptyRow = starts > stops; +emptyRow = []; +if any(isEmptyRow) + emptyRow = readEmptyRow(dataVector); +end + +rowValues = cell(numel(starts), 1); +for iRow = 1:numel(starts) + if isEmptyRow(iRow) + rowValues{iRow} = emptyRow; + continue + end + if isFileBacked + iRun = runOfPosition(startPositions(iRow)); + % The row's elements as positions in the block of its run. + blockRows = startPositions(iRow) - runFirst(iRun) + 1 + (0:(stops(iRow) - starts(iRow))); + block = selectRows(blocks{iRun}, blockRows, numSubscripts, rowAxis); else - columnRows = Vector.data(subscripts{:}); + block = selectRows(dataVector.data, starts(iRow):stops(iRow), numSubscripts, rowAxis); end + rowValues{iRow} = orientRows(dataVector, block); +end +end + +function emptyRow = readEmptyRow(dataVector) +% readEmptyRow - Value of an empty row of a ragged column. +% +% The empty row is read as an empty selection from the column instead of +% taken from a block of rows, so it has the type and shape of a direct read. +% A file-backed column reads its first element to learn the type of an +% empty selection, so a dataset without elements gets an empty of its data +% type instead. +data = dataVector.data; +if isa(data, 'types.untyped.DataStub') + isEmptyDataset = any(data.dims == 0); +elseif isa(data, 'types.untyped.DataPipe') && data.isBound + isEmptyDataset = any(size(data) == 0); +else + % Indexing in-memory data with an empty selection needs no elements. + isEmptyDataset = false; +end + +if isEmptyDataset + emptyRow = orientRows(dataVector, createEmptyValue(data.dataType)); +else + emptyRow = orientRows(dataVector, readRows(dataVector, zeros(1, 0))); +end +end - % A DataPipe can hold its rows along any axis. The table needs them along - % the first axis, or its variables would have unequal heights. - if isa(Vector.data, 'types.untyped.DataPipe') - columnRows = permute(columnRows, ... - circshift(1:ndims(columnRows), -(Vector.data.axis-1))); +function value = createEmptyValue(dataType) +% createEmptyValue - Empty value of the MATLAB type a dataset is read as. +% +% A compound dataset is read as a table with one variable per member, and +% every other dataset as an array. The empty array is 0x0, the shape an +% empty selection of a file-backed column has. +if isstruct(dataType) + compoundMemberNames = fieldnames(dataType); + compoundMembers = struct(); + for iMember = 1:numel(compoundMemberNames) + compoundMembers.(compoundMemberNames{iMember}) = createEmptyArray(dataType.(compoundMemberNames{iMember}), 1); end + value = struct2table(compoundMembers); else - assert(isa(Vector, 'types.hdmf_common.VectorIndex') || isa(Vector, 'types.core.VectorIndex'),... - 'NWB:DynamicTable:GetRow:InternalError',... - 'Internal VectorIndex Stack is not using VectorIndex objects!'); - if isa(Vector.data, 'types.untyped.DataStub') || isa(Vector.data, 'types.untyped.DataPipe') - stopInds = uint64(Vector.data.load(rowIndices)); + value = createEmptyArray(dataType, 0); +end +end + +function value = createEmptyArray(matlabType, numColumns) +% createEmptyArray - Empty array with no rows of the MATLAB type a dataset is read as. +% +% The types follow io.parseCompound, which builds the columns of a +% compound dataset without rows the same way. +switch matlabType + case {'char', 'cell'} + % Text and non-boolean enums are read as cell arrays. + value = cell(0, numColumns); + case 'logical' + value = false(0, numColumns); + otherwise + % Numeric types and the reference classes construct an empty + % instance from their class name. + value = feval([matlabType '.empty'], 0, numColumns); +end +end + +function values = readIndexValues(indexVector, rowsToRead) +% readIndexValues - Read elements of a VectorIndex as a double column vector. +if isa(indexVector.data, 'types.untyped.DataStub') || isa(indexVector.data, 'types.untyped.DataPipe') + values = indexVector.data.load(rowsToRead); +else + values = indexVector.data(rowsToRead); +end +values = double(values(:)); +end + +function block = readRows(vector, rowIndices) +% readRows - Read rows of a column, leaving the rows along the column's row axis. +[numSubscripts, rowAxis] = getRowDimension(vector); +block = selectRows(vector.data, rowIndices, numSubscripts, rowAxis); +end + +function selected = orientRows(vector, block) +% orientRows - Put the row axis first for a DataPipe column. +% +% A DataPipe can hold its rows along any axis. The table needs them along +% the first axis, or its variables would have unequal heights. +if isa(vector.data, 'types.untyped.DataPipe') + selected = permute(block, circshift(1:ndims(block), -(vector.data.axis-1))); +else + selected = block; +end +end + +function [numSubscripts, rowAxis] = getRowDimension(vector) +% getRowDimension - Number of subscripts used to index a column, and the one that selects rows. +data = vector.data; +if isa(data, 'types.untyped.DataStub') || isa(data, 'types.untyped.DataPipe') + if isa(data, 'types.untyped.DataStub') + dataSize = data.dims; + else + dataSize = data.internal.maxSize; + end + if length(dataSize) == 2 && dataSize(2) == 1 + % A column vector is indexed with one subscript. + numSubscripts = 1; + else + numSubscripts = length(dataSize); + end +else + if iscolumn(data) + % A column vector is indexed with one subscript. + numSubscripts = 1; + elseif istable(data) + % A compound column held as a table is indexed by row only. + numSubscripts = 1; else - stopInds = uint64(Vector.data(rowIndices)); + numSubscripts = ndims(data); end +end - startIndInd = rowIndices - 1; - zeroMask = startIndInd == 0; - startInds = zeros(size(startIndInd)); - if ~isempty(startIndInd(~zeroMask)) - if isa(Vector.data, 'types.untyped.DataStub') || isa(Vector.data, 'types.untyped.DataPipe') - startInds(~zeroMask) = Vector.data.load(startIndInd(~zeroMask)); - else - startInds(~zeroMask) = Vector.data(startIndInd(~zeroMask)); - end +if isa(data, 'types.untyped.DataPipe') + rowAxis = data.axis; +else + rowAxis = numSubscripts; +end +end + +function selected = selectRows(data, rowIndices, numSubscripts, rowAxis) +% selectRows - Select rows of in-memory or file-backed column data. +% +% The same indexing reads rows from the column and takes rows from a block +% already read from it. +subscripts = repmat({':'}, 1, numSubscripts); +subscripts{rowAxis} = rowIndices; + +if (isstruct(data) && isscalar(data)) || istable(data) + if istable(data) + selected = table(); + fields = data.Properties.VariableNames; + else + selected = struct(); + fields = fieldnames(data); end - startInds = startInds + 1; - - columnRows = cell(length(rowIndices), 1); - for iRange = 1:length(rowIndices) - startInd = startInds(iRange); - stopInd = stopInds(iRange); - columnRows{iRange} = getColumnRows(DynamicTable,... - vectorChainNames(1:(end-1)),... - startInd:stopInd); + + for iField = 1:length(fields) + fieldName = fields{iField}; + memberData = data.(fieldName); + selected.(fieldName) = memberData(subscripts{:}); end +else + selected = data(subscripts{:}); end end