diff --git a/+matnwb/+neurodata/DynamicTableBase.m b/+matnwb/+neurodata/DynamicTableBase.m index d638191ec..f5035390b 100644 --- a/+matnwb/+neurodata/DynamicTableBase.m +++ b/+matnwb/+neurodata/DynamicTableBase.m @@ -9,6 +9,7 @@ properties (Abstract) id colnames + vectordata end methods @@ -112,6 +113,145 @@ function addColumn(obj, columnName, columnVector) types.util.dynamictable.addVarargColumn(obj, columnVectorPairs{:}); end + function addRaggedArray(obj, columnName, data, options) + % addRaggedArray - Add a ragged-array column to the DynamicTable. + % + % A ragged array stores a variable number of elements per row. The + % values are held in a single VectorData column, and a companion + % VectorIndex ('_index') marks each row's boundary. See + % the "Tables and ragged arrays" section of the NWB format + % specification. + % + % Syntax: + % dynamicTable.addRaggedArray(columnName, data) build and add a + % ragged column named columnName, plus its VectorIndex. + % + % dynamicTable.addRaggedArray(__, Name, Value) provide optional + % arguments (see below). + % + % Input Arguments: + % - columnName (string) - + % Name of the new column. + % + % - data (cell) - + % A cell array with one cell per row; each cell holds that row's + % elements (e.g. {[1 2 3], [4 5]} for a 2-row table). + % + % Name-Value Arguments: + % - description (string) - + % Description stored on the VectorData column. + % + % - table (DynamicTable) - + % If provided, the column is created as a DynamicTableRegion that + % references this table (row indices) instead of a VectorData. + % + % See also util.create_indexed_column, addColumn, addDoublyRaggedArray + + arguments + obj (1,1) {matnwb.common.validation.mustBeDynamicTable} + columnName (1,1) string + data cell + options.description (1,1) string = "no description" + options.table = [] + end + + if isempty(options.table) + [vector, index] = util.create_indexed_column( ... + data, char(options.description)); + else + [vector, index] = util.create_indexed_column( ... + data, char(options.description), options.table); + end + obj.addColumn(columnName, vector, columnName + "_index", index); + end + + function addDoublyRaggedArray(obj, columnName, data, options) + % addDoublyRaggedArray - Add a doubly-ragged-array column to the DynamicTable. + % + % A doubly-ragged array stores, for each row, a variable number of + % sub-groups, each holding a variable number of fixed-length elements + % (e.g. the Units table 'waveforms' column: per unit, a variable number + % of spike events, each with a waveform per electrode). It is backed by + % a VectorData column and two VectorIndex levels + % ('_index' over sub-groups and '_index_index' + % over rows). See the "Doubly ragged arrays" section of the NWB format + % specification. + % + % Syntax: + % dynamicTable.addDoublyRaggedArray(columnName, data) build and add a + % doubly-ragged column named columnName, plus its two VectorIndex + % levels. + % + % Input Arguments: + % - columnName (string) - + % Name of the new column. + % + % - data (cell) - + % A cell array with one cell per row. Each cell is either a numeric + % matrix [nSubGroups x nSamples] (one element per sub-group), or a + % cell array whose j-th entry is a [nElements x nSamples] matrix for + % sub-group j. See util.create_doubly_indexed_column. + % + % Name-Value Arguments: + % - description (string) - + % Description stored on the VectorData column. + % + % See also util.create_doubly_indexed_column, addColumn, addRaggedArray + + arguments + obj (1,1) {matnwb.common.validation.mustBeDynamicTable} + columnName (1,1) string + data cell + options.description (1,1) string = "no description" + end + + [vector, index, indexIndex] = ... + util.create_doubly_indexed_column(data, options.description); + + % The number of rows equals the length of the outermost index. + rowCount = numel(indexIndex.data); + + % Initialize id for a new table before checking editability (which + % inspects the id column). + if isempty(obj.id) || isempty(obj.id.data) + types.util.dynamictable.internal.initDynamicTableId(obj, rowCount); + end + + obj.assertIsEditable('NWB:DynamicTable:AddDoublyRaggedArray:Uneditable') + + tableHeight = types.util.dynamictable.internal.getColumnHeight(obj.id); + assert(rowCount == tableHeight, ... + 'NWB:DynamicTable:AddDoublyRaggedArray:MissingRows', ... + 'Column `%s` has %d rows, but the table height is %d.', ... + columnName, rowCount, tableHeight) + + % Assign the data column and both index levels to the appropriate + % storage (a typed property for schema-defined columns such as + % Units.waveforms, otherwise the generic vectordata set). addColumn + % is not used here because its height check only follows a single + % index level. + names = [columnName, columnName + "_index", columnName + "_index_index"]; + values = {vector, index, indexIndex}; + for iName = 1:numel(names) + thisName = char(names(iName)); + storageTarget = types.util.dynamictable.resolveColumnStorage( ... + obj, thisName, values{iName}); + switch storageTarget + case 'property' + obj.(thisName) = values{iName}; + case 'vectordata' + obj.vectordata.set(thisName, values{iName}); + end + end + + % Only the data column is listed in colnames; the index levels are + % implicit. Schema-defined columns are added by a property post-set + % hook, so guard against duplicates. + if ~any(strcmp(obj.colnames, char(columnName))) + obj.colnames{end+1} = char(columnName); + end + end + function row = getRow(obj, rowIndices, options) % getRow - Return one or more DynamicTable rows. % diff --git a/+tests/+system/PyNWBIOTest.py b/+tests/+system/PyNWBIOTest.py index 477845155..439404a1c 100644 --- a/+tests/+system/PyNWBIOTest.py +++ b/+tests/+system/PyNWBIOTest.py @@ -214,10 +214,30 @@ def getContainer(self, file): return file -class UnitTimesIOTest(PyNWBIOTest): +class UnitsTableIOTest(PyNWBIOTest): def addContainer(self, file): - self.file.units = Units('units', waveform_rate=1., resolution=3.) - self.file.units.add_unit(waveform_mean=[5], waveform_sd=[7], waveforms=np.full((1, 1), 9), - spike_times=[11]) + # Build a Units table with multi-electrode, doubly-ragged waveforms + # using the idiomatic add_unit method. This mirrors the MatNWB + # addRaggedArray / addDoublyRaggedArray construction in UnitsTableIOTest.m. + # + # Each unit's waveforms are passed as a 3-D array with shape + # (num_spikes, num_electrodes, num_samples). Note the LAST dimension is + # samples, not electrodes: the `waveforms` column is doubly indexed, so + # add_unit stores a 2-D [num_waveforms, num_samples] dataset with the + # electrode dimension carried by waveforms_index. (The add_unit docstring + # line "the third dimension shows ... different electrodes" describes a + # non-indexed 3-D dataset and does not match this input order.) + # + # Unit 1 has 2 spikes, unit 2 has 3 spikes; each spike has 2 electrodes + # and 3 samples. + self.file.units = Units(name='units', waveform_rate=1., resolution=3., + description='data on spiking units') + self.file.units.add_unit(spike_times=[1., 2.], + waveform_mean=[1., 2., 3.], waveform_sd=[7., 8., 9.], + waveforms=np.arange(1, 13, dtype=np.int32).reshape(2, 2, 3)) + self.file.units.add_unit(spike_times=[3., 4., 5.], + waveform_mean=[4., 5., 6.], waveform_sd=[10., 11., 12.], + waveforms=np.arange(13, 31, dtype=np.int32).reshape(3, 2, 3)) + def getContainer(self, file): return file.units diff --git a/+tests/+system/UnitTimesIOTest.m b/+tests/+system/UnitTimesIOTest.m deleted file mode 100644 index ab415e066..000000000 --- a/+tests/+system/UnitTimesIOTest.m +++ /dev/null @@ -1,57 +0,0 @@ -classdef UnitTimesIOTest < tests.system.PyNWBIOTest - methods - function addContainer(~, file) - file.units = types.core.Units( ... - 'colnames', {'spike_times'; 'waveform_mean'; 'waveform_sd'; 'waveforms'} ... - , 'description', 'data on spiking units' ... - ); - - file.units.addRow( ... - 'spike_times', 11 ... - , 'waveforms', int32(9) ... - , 'waveform_sd', 7 ... - , 'waveform_mean', 5 ... - ); - % the following is solely for a2a comparison with files. - file.units.spike_times.description = 'the spike times for each unit in seconds'; - file.units.waveforms.description = ['Individual waveforms for each spike. ' ... - 'If the dataset is three-dimensional, the third dimension shows the response ' ... - 'from different electrodes that all observe this unit simultaneously. ' ... - 'In this case, the `electrodes` column of this Units table should be used to ' ... - 'indicate which electrodes are associated with this unit, ' ... - 'and the electrodes dimension here should be in the same order as the ' ... - 'electrodes referenced in the `electrodes` column of this table.']; - file.units.waveform_sd.description = ['the spike waveform standard deviation for each ' ... - 'spike unit']; - file.units.waveform_mean.description = 'the spike waveform mean for each spike unit'; - file.units.spike_times_index = types.hdmf_common.VectorIndex( ... - 'target', types.untyped.ObjectView(file.units.spike_times) ... - , 'description', 'Index for VectorData ''spike_times''' ... - , 'data', 1 ... - ); - file.units.waveforms_index = types.hdmf_common.VectorIndex( ... - 'target', types.untyped.ObjectView(file.units.waveforms) ... - , 'description', 'Index for VectorData ''waveforms''' ... - , 'data', 1 ... - ); - file.units.waveforms_index_index = types.hdmf_common.VectorIndex( ... - 'target', types.untyped.ObjectView(file.units.waveforms_index) ... - , 'description', 'Index for VectorData ''waveforms_index''' ... - , 'data', 1 ... - ); - - % Set optional Units table dataset attributes via promoted container - % API - file.units.spike_times_resolution = 3; - file.units.waveform_mean_sampling_rate = 1; - file.units.waveform_sd_sampling_rate = 1; - - % Skip waveforms_sampling_rate because PyNWB does not export it. - file.units.waveforms_sampling_rate = 1; - end - - function c = getContainer(~, file) - c = file.units; - end - end -end diff --git a/+tests/+system/UnitsTableIOTest.m b/+tests/+system/UnitsTableIOTest.m new file mode 100644 index 000000000..a914da053 --- /dev/null +++ b/+tests/+system/UnitsTableIOTest.m @@ -0,0 +1,63 @@ +classdef UnitsTableIOTest < tests.system.PyNWBIOTest + methods + function addContainer(~, file) + % Build a Units table with multi-electrode, doubly-ragged waveforms + % using the idiomatic addRaggedArray / addDoublyRaggedArray methods. + % This mirrors PyNWB's Units.add_unit in PyNWBIOTest.py so the two + % round-tripped containers compare equal. + file.units = types.core.Units('description', 'data on spiking units'); + + % spike_times: a ragged column, one list of times per unit. + file.units.addRaggedArray('spike_times', {[1 2], [3 4 5]}); + + % waveforms: a doubly-ragged column. data{unit}{spike} is a + % [numElectrodes x numSamples] matrix (one row per electrode's + % waveform). This matches PyNWB add_unit's per-unit 3-D input + % (num_spikes, num_electrodes, num_samples): both store a 2-D + % [num_waveforms, num_samples] dataset with the electrode dimension + % carried by waveforms_index. + % Unit 1 has 2 spikes, unit 2 has 3 spikes; each spike has 2 + % electrodes and 3 samples. + waveforms = { ... + {int32([1 2 3; 4 5 6]), int32([7 8 9; 10 11 12])}, ... + {int32([13 14 15; 16 17 18]), int32([19 20 21; 22 23 24]), int32([25 26 27; 28 29 30])} ... + }; + file.units.addDoublyRaggedArray('waveforms', waveforms); + + % waveform_mean / waveform_sd: fixed 2-D columns [numSamples x numUnits]. + file.units.waveform_mean = types.hdmf_common.VectorData( ... + 'description', 'the spike waveform mean for each spike unit', ... + 'data', [1 4; 2 5; 3 6]); + file.units.waveform_sd = types.hdmf_common.VectorData( ... + 'description', 'the spike waveform standard deviation for each spike unit', ... + 'data', [7 10; 8 11; 9 12]); + + % Match PyNWB's column order and auto-generated descriptions so the + % round-tripped containers compare equal. + file.units.colnames = {'spike_times'; 'waveform_mean'; 'waveform_sd'; 'waveforms'}; + file.units.spike_times.description = 'the spike times for each unit in seconds'; + file.units.spike_times_index.description = 'Index for VectorData ''spike_times'''; + file.units.waveforms.description = ['Individual waveforms for each spike. ' ... + 'If the dataset is three-dimensional, the third dimension shows the response ' ... + 'from different electrodes that all observe this unit simultaneously. ' ... + 'In this case, the `electrodes` column of this Units table should be used to ' ... + 'indicate which electrodes are associated with this unit, ' ... + 'and the electrodes dimension here should be in the same order as the ' ... + 'electrodes referenced in the `electrodes` column of this table.']; + file.units.waveforms_index.description = 'Index for VectorData ''waveforms'''; + file.units.waveforms_index_index.description = 'Index for VectorData ''waveforms_index'''; + + % Optional Units dataset attributes (match PyNWB waveform_rate / resolution). + file.units.spike_times_resolution = 3; + file.units.waveform_mean_sampling_rate = 1; + file.units.waveform_sd_sampling_rate = 1; + + % Skip waveforms_sampling_rate because PyNWB does not export it. + file.units.waveforms_sampling_rate = 1; + end + + function c = getContainer(~, file) + c = file.units; + end + end +end diff --git a/+tests/+unit/createDoublyIndexedColumnTest.m b/+tests/+unit/createDoublyIndexedColumnTest.m new file mode 100644 index 000000000..056d9ec1c --- /dev/null +++ b/+tests/+unit/createDoublyIndexedColumnTest.m @@ -0,0 +1,172 @@ +classdef createDoublyIndexedColumnTest < tests.abstract.NwbTestCase +% createDoublyIndexedColumnTest - Unit tests for util.create_doubly_indexed_column + + methods (TestClassSetup) + function setupClass(testCase) + testCase.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); + end + end + + methods (Test) + function testSingleElectrodeShortcut(testCase) + % Matrix-per-row form: each row's [nSpikes x nSamples] matrix, one + % waveform per spike (single electrode). + nSamples = 4; + unit1 = reshape(1:(3*nSamples), 3, nSamples); % 3 spikes + unit2 = 100 + reshape(1:(4*nSamples), 4, nSamples); % 4 spikes + + [vector, index, indexIndex] = ... + util.create_doubly_indexed_column({unit1, unit2}, 'waveforms'); + + % Data: samples along dim 1, one column per waveform (7 total). + testCase.verifyEqual(size(vector.data), [nSamples, 7]); + testCase.verifyEqual(vector.data(:, 1), unit1(1, :).'); + testCase.verifyEqual(vector.data(:, 4), unit2(1, :).'); + + % Inner index: one entry per spike event (single waveform each). + testCase.verifyEqual(index.data, uint64((1:7)')); + % Outer index: one entry per unit (cumulative spike counts). + testCase.verifyEqual(indexIndex.data, uint64([3; 7])); + + % Targets are wired: index -> vector, indexIndex -> index. + testCase.verifySameHandle(index.target.target, vector); + testCase.verifySameHandle(indexIndex.target.target, index); + end + + function testNDArrayShortcut(testCase) + % The first dimension is the subgroup axis; trailing dimensions are + % preserved as the fixed element payload. + unit1 = reshape(1:(2*3*4), 2, 3, 4); + unit2 = 100 + reshape(1:(1*3*4), 1, 3, 4); + + [vector, index, indexIndex] = ... + util.create_doubly_indexed_column({unit1, unit2}, 'volumes'); + + testCase.verifyEqual(size(vector.data), [4, 3, 3]); + testCase.verifyEqual(vector.data(:, :, 1), reshape(unit1(1, :, :), 3, 4).'); + testCase.verifyEqual(vector.data(:, :, 3), reshape(unit2(1, :, :), 3, 4).'); + testCase.verifyEqual(index.data, uint64((1:3)')); + testCase.verifyEqual(indexIndex.data, uint64([2; 3])); + end + + function testGeneralNestedCell(testCase) + % Cell-per-row form: DATA{i}{j} is [nElements x nSamples]. + nSamples = 4; + % unit1: event1 has 2 electrodes, event2 has 1 electrode + % unit2: event1 has 3 electrodes + unit1 = {ones(2, nSamples), 2 * ones(1, nSamples)}; + unit2 = {3 * ones(3, nSamples)}; + + [vector, index, indexIndex] = ... + util.create_doubly_indexed_column({unit1, unit2}); + + testCase.verifyEqual(size(vector.data), [nSamples, 6]); + % Inner index counts elements (electrodes) per sub-group: 2, 1, 3. + testCase.verifyEqual(index.data, uint64([2; 3; 6])); + % Outer index counts sub-groups per row: 2, 1. + testCase.verifyEqual(indexIndex.data, uint64([2; 3])); + end + + function testNDArrayNestedCell(testCase) + % Cell-per-row form supports arbitrary trailing element dimensions. + unit1 = {reshape(1:(2*3*4), 2, 3, 4), ... + 100 + reshape(1:(1*3*4), 1, 3, 4)}; + unit2 = {200 + reshape(1:(3*3*4), 3, 3, 4)}; + + [vector, index, indexIndex] = ... + util.create_doubly_indexed_column({unit1, unit2}); + + testCase.verifyEqual(size(vector.data), [4, 3, 6]); + testCase.verifyEqual(vector.data(:, :, 1), reshape(unit1{1}(1, :, :), 3, 4).'); + testCase.verifyEqual(vector.data(:, :, 3), reshape(unit1{2}(1, :, :), 3, 4).'); + testCase.verifyEqual(vector.data(:, :, 6), reshape(unit2{1}(3, :, :), 3, 4).'); + testCase.verifyEqual(index.data, uint64([2; 3; 6])); + testCase.verifyEqual(indexIndex.data, uint64([2; 3])); + end + + function testNDArrayExportKeepsSchemaDimensionOrder(testCase) + unit1 = {reshape(1:(2*3*4), 2, 3, 4), ... + 100 + reshape(1:(1*3*4), 1, 3, 4)}; + unit2 = {200 + reshape(1:(3*3*4), 3, 3, 4)}; + + vector = util.create_doubly_indexed_column({unit1, unit2}); + + fileName = testCase.getRandomFilename(); + fileId = H5F.create(fileName); + cleanup = onCleanup(@() H5F.close(fileId)); + io.writeDataset(fileId, '/wf', vector.data); + delete(cleanup) + + fileId = H5F.open(fileName, 'H5F_ACC_RDONLY', 'H5P_DEFAULT'); + fileCleanup = onCleanup(@() H5F.close(fileId)); + datasetId = H5D.open(fileId, '/wf'); + datasetCleanup = onCleanup(@() H5D.close(datasetId)); + spaceId = H5D.get_space(datasetId); + spaceCleanup = onCleanup(@() H5S.close(spaceId)); + [~, h5Dimensions, ~] = H5S.get_simple_extent_dims(spaceId); + testCase.verifyEqual(h5Dimensions, [6, 3, 4]); + delete(spaceCleanup) + delete(datasetCleanup) + delete(fileCleanup) + end + + function testEmptyRow(testCase) + % A row with no sub-groups is represented by [] (or {}). + nSamples = 4; + unit1 = reshape(1:(2*nSamples), 2, nSamples); + unit3 = reshape(1:(3*nSamples), 3, nSamples); + + [vector, index, indexIndex] = ... + util.create_doubly_indexed_column({unit1, [], unit3}); + + testCase.verifyEqual(size(vector.data), [nSamples, 5]); + testCase.verifyEqual(index.data, uint64((1:5)')); + % Middle (empty) unit repeats the previous cumulative value. + testCase.verifyEqual(indexIndex.data, uint64([2; 2; 5])); + end + + function testInconsistentSampleLengthErrors(testCase) + testCase.verifyError( ... + @() util.create_doubly_indexed_column({ones(3, 4), ones(2, 5)}), ... + 'NWB:CreateDoublyIndexedColumn:InconsistentSampleLength'); + end + + function testInvalidRowTypeErrors(testCase) + testCase.verifyError( ... + @() util.create_doubly_indexed_column({ones(3, 4), "not numeric"}), ... + 'NWB:CreateDoublyIndexedColumn:InvalidRow'); + end + + function testRoundTripInUnitsTable(testCase) + % Build a Units table with the doubly-indexed waveforms column, + % export, read back, and verify the on-disk index structure. + nSamples = 4; + unit1 = reshape(1:(3*nSamples), 3, nSamples); + unit2 = 100 + reshape(1:(4*nSamples), 4, nSamples); + + [vector, index, indexIndex] = ... + util.create_doubly_indexed_column({unit1, unit2}, 'spike waveforms'); + + units = types.core.Units( ... + 'colnames', {'waveforms'}, ... + 'description', 'units', ... + 'waveforms', vector, ... + 'waveforms_index', index, ... + 'waveforms_index_index', indexIndex, ... + 'id', types.hdmf_common.ElementIdentifiers('data', [0; 1])); + + nwb = NwbFile( ... + 'identifier', 'doubly_indexed_test', ... + 'session_description', 'test', ... + 'session_start_time', datetime(2024, 1, 1, 'TimeZone', 'local')); + nwb.units = units; + + fileName = testCase.getRandomFilename(); + nwbExport(nwb, fileName); + + back = nwbRead(fileName, 'ignorecache'); + testCase.verifyEqual(back.units.waveforms_index.data.load(), uint64((1:7)')); + testCase.verifyEqual(back.units.waveforms_index_index.data.load(), uint64([3; 7])); + end + end +end diff --git a/+tests/+unit/dynamicTableRaggedArrayTest.m b/+tests/+unit/dynamicTableRaggedArrayTest.m new file mode 100644 index 000000000..6f5b4c7b5 --- /dev/null +++ b/+tests/+unit/dynamicTableRaggedArrayTest.m @@ -0,0 +1,106 @@ +classdef dynamicTableRaggedArrayTest < tests.abstract.NwbTestCase +% dynamicTableRaggedArrayTest - Tests for the DynamicTable addRaggedArray and +% addDoublyRaggedArray convenience methods. + + methods (TestClassSetup) + function setupClass(testCase) + testCase.applyFixture(matlab.unittest.fixtures.WorkingFolderFixture); + end + end + + methods (Test) + function testAddRaggedArray(testCase) + dt = types.hdmf_common.DynamicTable('description', 'test'); + dt.addRaggedArray('spikes', {[1 2 3], [4 5]}, 'description', 'spike times'); + + testCase.verifyTrue(any(strcmp(dt.colnames, 'spikes'))); + testCase.verifyEqual(numel(dt.id.data), 2); + % Cumulative row boundaries for a 3- and 2-element row. + testCase.verifyEqual(dt.vectordata.get('spikes_index').data, uint64([3; 5])); + testCase.verifyEqual(dt.vectordata.get('spikes').description, 'spike times'); + end + + function testAddRaggedArrayWithTableRegion(testCase) + target = types.hdmf_common.DynamicTable( ... + 'description', 'target', ... + 'colnames', {'x'}, ... + 'x', types.hdmf_common.VectorData('description', 'x', 'data', (1:5)'), ... + 'id', types.hdmf_common.ElementIdentifiers('data', (0:4)')); + + dt = types.hdmf_common.DynamicTable('description', 'test'); + dt.addRaggedArray('regions', {[0 1], [2 3 4]}, ... + 'description', 'regions', 'table', target); + + column = dt.vectordata.get('regions'); + testCase.verifyClass(column, 'types.hdmf_common.DynamicTableRegion'); + testCase.verifyEqual(dt.vectordata.get('regions_index').data, uint64([2; 5])); + end + + function testAddDoublyRaggedArrayGeneric(testCase) + nSamples = 4; + unit1 = reshape(1:(3*nSamples), 3, nSamples); + unit2 = 100 + reshape(1:(4*nSamples), 4, nSamples); + + dt = types.hdmf_common.DynamicTable('description', 'test'); + dt.addDoublyRaggedArray('wf', {unit1, unit2}, 'description', 'waveforms'); + + testCase.verifyEqual(dt.colnames, {'wf'}); + testCase.verifyEqual(numel(dt.id.data), 2); + testCase.verifyEqual(size(dt.vectordata.get('wf').data), [nSamples, 7]); + testCase.verifyEqual(dt.vectordata.get('wf_index').data, uint64((1:7)')); + testCase.verifyEqual(dt.vectordata.get('wf_index_index').data, uint64([3; 7])); + end + + function testAddDoublyRaggedArrayWithNDArrayElements(testCase) + unit1 = {reshape(1:(2*3*4), 2, 3, 4), ... + 100 + reshape(1:(1*3*4), 1, 3, 4)}; + unit2 = {200 + reshape(1:(3*3*4), 3, 3, 4)}; + + dt = types.hdmf_common.DynamicTable('description', 'test'); + dt.addDoublyRaggedArray('wf', {unit1, unit2}, 'description', 'volumes'); + + testCase.verifyEqual(size(dt.vectordata.get('wf').data), [4, 3, 6]); + testCase.verifyEqual(dt.vectordata.get('wf').data(:, :, 1), ... + reshape(unit1{1}(1, :, :), 3, 4).'); + testCase.verifyEqual(dt.vectordata.get('wf').data(:, :, 6), ... + reshape(unit2{1}(3, :, :), 3, 4).'); + testCase.verifyEqual(dt.vectordata.get('wf_index').data, uint64([2; 3; 6])); + testCase.verifyEqual(dt.vectordata.get('wf_index_index').data, uint64([2; 3])); + end + + function testAddDoublyRaggedArrayOnUnitsRoundTrip(testCase) + nSamples = 4; + unit1 = reshape(1:(3*nSamples), 3, nSamples); + unit2 = 100 + reshape(1:(4*nSamples), 4, nSamples); + + units = types.core.Units('colnames', {}, 'description', 'units'); + units.addDoublyRaggedArray('waveforms', {unit1, unit2}, ... + 'description', 'spike waveforms'); + + testCase.verifyTrue(any(strcmp(units.colnames, 'waveforms'))); + + nwb = NwbFile( ... + 'identifier', 'ragged_method_test', ... + 'session_description', 'test', ... + 'session_start_time', datetime(2024, 1, 1, 'TimeZone', 'local')); + nwb.units = units; + + fileName = testCase.getRandomFilename(); + nwbExport(nwb, fileName); + + back = nwbRead(fileName, 'ignorecache'); + testCase.verifyEqual(back.units.waveforms_index.data.load(), uint64((1:7)')); + testCase.verifyEqual(back.units.waveforms_index_index.data.load(), uint64([3; 7])); + end + + function testAddDoublyRaggedArrayHeightMismatchErrors(testCase) + dt = types.hdmf_common.DynamicTable('description', 'test'); + dt.addColumn('a', types.hdmf_common.VectorData( ... + 'description', 'a', 'data', [1 2 3]')); % 3 rows + + testCase.verifyError( ... + @() dt.addDoublyRaggedArray('wf', {reshape(1:8, 2, 4), reshape(1:8, 2, 4)}), ... + 'NWB:DynamicTable:AddDoublyRaggedArray:MissingRows'); + end + end +end diff --git a/+tests/+unit/howToRaggedArrayExamplesTest.m b/+tests/+unit/howToRaggedArrayExamplesTest.m new file mode 100644 index 000000000..9868ec0da --- /dev/null +++ b/+tests/+unit/howToRaggedArrayExamplesTest.m @@ -0,0 +1,18 @@ +classdef howToRaggedArrayExamplesTest < tests.abstract.NwbTestCase +% howToRaggedArrayExamplesTest - Executes the runnable code embedded in the +% "Storing Ragged and Doubly-Ragged Array Columns" how-to guide, keeping the +% documented snippets correct against the API. The example function contains +% its own assertions. + + methods (Test) + function testRaggedArrayExamplesRun(testCase) + exampleDir = fullfile(misc.getMatnwbDir(), ... + 'docs', 'source', 'pages', 'how_to', 'examples'); + testCase.applyFixture(matlab.unittest.fixtures.PathFixture(exampleDir)); + + % Runs to completion only if every embedded snippet and its + % assertions pass. + ragged_arrays_examples(); + end + end +end diff --git a/+util/create_doubly_indexed_column.m b/+util/create_doubly_indexed_column.m new file mode 100644 index 000000000..6135c14e2 --- /dev/null +++ b/+util/create_doubly_indexed_column.m @@ -0,0 +1,149 @@ +function [vector, index, index_index] = create_doubly_indexed_column(data, description) +%CREATE_DOUBLY_INDEXED_COLUMN creates the objects for a doubly-ragged column +%in an NWB DynamicTable (a VectorData indexed by two levels of VectorIndex). +% +% Some NWB columns are doubly ragged: each table row contains a variable +% number of sub-groups, and each sub-group contains a variable number of +% fixed-length elements. The Units table 'waveforms' column is the canonical +% example: each unit (row) has a number of spike events (sub-groups), and +% each spike event has a waveform for each electrode (elements), where every +% waveform has the same number of samples. +% +% [VECTOR, INDEX, INDEX_INDEX] = CREATE_DOUBLY_INDEXED_COLUMN(DATA) returns a +% VectorData (VECTOR) and two VectorIndex objects: INDEX has one entry per +% sub-group and targets VECTOR; INDEX_INDEX has one entry per row and targets +% INDEX. Assign them to a column and its '_index'/'_index_index' +% properties (e.g. waveforms, waveforms_index, waveforms_index_index). +% +% DATA is a cell array with one cell per table row. Each cell is either: +% - a numeric array [nSubGroups x trailingDimensions] (single-element +% shortcut): each slice along dimension 1 is one sub-group containing +% exactly one element. This matches spike-sorted units on a single +% electrode, where DATA{i} is that unit's [nSpikes x nSamples] waveform +% matrix. +% - a cell array where DATA{i}{j} is a numeric array +% [nElements x trailingDimensions] holding the elements of sub-group j +% (general, multi-element case). +% +% All elements across all rows must have the same trailing dimensions. A row +% may be empty ([] or {}) to represent a row with no sub-groups. +% +% [VECTOR, INDEX, INDEX_INDEX] = CREATE_DOUBLY_INDEXED_COLUMN(DATA, DESCRIPTION) +% sets the string DESCRIPTION on the returned VectorData. +% +% Example (single electrode, 2 units with 3 and 4 spikes, 40 samples): +% waveforms = {unit1Waveforms, unit2Waveforms}; % each [nSpikes x 40] +% [wf, wfIdx, wfIdxIdx] = util.create_doubly_indexed_column(waveforms, 'spike waveforms'); +% units.waveforms = wf; +% units.waveforms_index = wfIdx; +% units.waveforms_index_index = wfIdxIdx; +% +% See also util.create_indexed_column + + arguments + data cell + description (1,1) string = "no description" + end + + numRows = numel(data); + outerCounts = zeros(numRows, 1); % number of sub-groups per row + rowChunks = cell(1, numRows); % data chunks with ragged axis last + rowInnerCounts = cell(1, numRows); % element count per sub-group, per row + payloadSize = []; % fixed trailing dimensions, once known + + for iRow = 1:numRows + rowData = data{iRow}; + if isnumeric(rowData) || islogical(rowData) + [chunk, counts, payloadSize] = fromArray(rowData, payloadSize, iRow); + elseif iscell(rowData) + [chunk, counts, payloadSize] = fromCell(rowData, payloadSize, iRow); + else + error("NWB:CreateDoublyIndexedColumn:InvalidRow", ... + "Each element of DATA must be a numeric array or a cell array. " + ... + "Element %d is a %s.", iRow, class(rowData)); + end + rowChunks{iRow} = chunk; + rowInnerCounts{iRow} = counts; + outerCounts(iRow) = numel(counts); + end + + allData = concatenateChunks(rowChunks, payloadSize); + innerCounts = vertcat(rowInnerCounts{:}); % one entry per sub-group + + vector = types.hdmf_common.VectorData( ... + 'description', char(description), ... + 'data', allData); + + index = types.hdmf_common.VectorIndex( ... + 'description', 'Index into the data column, one entry per sub-group', ... + 'target', types.untyped.ObjectView(vector), ... + 'data', uint64(cumsum(innerCounts))); + + index_index = types.hdmf_common.VectorIndex( ... + 'description', 'Index into the index column, one entry per row', ... + 'target', types.untyped.ObjectView(index), ... + 'data', uint64(cumsum(outerCounts))); +end + +function [chunk, counts, payloadSize] = fromArray(rowData, payloadSize, iRow) + % Shortcut form: [nSubGroups x trailingDimensions], one element per sub-group. + if isempty(rowData) + chunk = []; + counts = zeros(0, 1); + return + end + payloadSize = checkPayloadSize(getPayloadSize(rowData), payloadSize, iRow); + chunk = moveFirstDimensionToLast(rowData); + counts = ones(size(rowData, 1), 1); % one element per sub-group +end + +function [chunk, counts, payloadSize] = fromCell(rowData, payloadSize, iRow) + % General form: rowData{j} = [nElements x trailingDimensions] for sub-group j. + numGroups = numel(rowData); + chunks = cell(1, numGroups); + counts = zeros(numGroups, 1); + for j = 1:numGroups + element = rowData{j}; + if ~(isnumeric(element) || islogical(element)) + error("NWB:CreateDoublyIndexedColumn:InvalidElement", ... + "DATA{%d}{%d} must be a numeric array.", iRow, j); + end + if ~isempty(element) + payloadSize = checkPayloadSize(getPayloadSize(element), payloadSize, iRow); + chunks{j} = moveFirstDimensionToLast(element); + else + chunks{j} = []; + end + counts(j) = size(element, 1); + end + chunk = concatenateChunks(chunks, payloadSize); +end + +function payloadSize = checkPayloadSize(thisSize, payloadSize, iRow) + if isempty(payloadSize) + payloadSize = thisSize; + elseif ~isequal(thisSize, payloadSize) + error("NWB:CreateDoublyIndexedColumn:InconsistentSampleLength", ... + "All elements must have the same trailing dimensions. " + ... + "Expected [%s], but an element in row %d has [%s].", ... + join(string(payloadSize), " "), iRow, join(string(thisSize), " ")); + end +end + +function payloadSize = getPayloadSize(array) + arraySize = size(array); + payloadSize = arraySize(2:end); +end + +function data = moveFirstDimensionToLast(data) + data = permute(data, ndims(data):-1:1); +end + +function data = concatenateChunks(chunks, payloadSize) + nonEmptyChunks = chunks(~cellfun(@isempty, chunks)); + if isempty(nonEmptyChunks) + data = []; + return + end + data = cat(numel(payloadSize) + 1, nonEmptyChunks{:}); +end diff --git a/docs/source/pages/how_to/examples/ragged_arrays_examples.m b/docs/source/pages/how_to/examples/ragged_arrays_examples.m new file mode 100644 index 000000000..f230fdb1f --- /dev/null +++ b/docs/source/pages/how_to/examples/ragged_arrays_examples.m @@ -0,0 +1,88 @@ +function ragged_arrays_examples() +% ragged_arrays_examples - Runnable snippets for the "Storing Ragged and +% Doubly-Ragged Array Columns" how-to guide. +% +% Each region between "% snippet: " and "% end snippet" markers is +% embedded into docs/source/pages/how_to/ragged_arrays.rst via literalinclude. +% The assertions keep the guide's code correct: this function is executed by +% tests/unit/howToRaggedArrayExamplesTest. + + electrodesTable = local_electrodesTable(); + + % snippet: ragged-spike-times + units = types.core.Units('colnames', {}, 'description', 'units'); + units.addRaggedArray('spike_times', {[0.1 0.2 0.3], [0.5 0.6]}, ... + 'description', 'spike times'); + % end snippet + assert(isequal(units.spike_times_index.data(:).', uint64([3 5]))) + + % snippet: ragged-electrodes-region + units.addRaggedArray('electrodes', {[0 1 2], [0 1 2]}, 'table', electrodesTable); + % end snippet + assert(isa(units.electrodes, 'types.hdmf_common.DynamicTableRegion')) + + % ---- Doubly-ragged, single electrode per unit ---- + % Each unit's matrix is [numWaveforms x numSamples]; with one electrode, + % numWaveforms is the number of spikes. + numSamples = 40; + unit1 = rand(3, numSamples); % 3 spikes + unit2 = rand(4, numSamples); % 4 spikes + + % snippet: doubly-single-electrode + units = types.core.Units('colnames', {}, 'description', 'units'); + units.addDoublyRaggedArray('waveforms', {unit1, unit2}, ... + 'description', 'spike waveforms'); + % end snippet + assert(isequal(units.waveforms_index.data(:).', uint64(1:7))) + assert(isequal(units.waveforms_index_index.data(:).', uint64([3 7]))) + + % ---- Doubly-ragged, multiple channels per spike ---- + % Each spike's matrix is [numWaveforms x numSamples]; with several + % electrodes, numWaveforms is the number of electrodes. + numElectrodes = 3; + m1 = { rand(numElectrodes, numSamples), rand(numElectrodes, numSamples) }; + m2 = { rand(numElectrodes, numSamples), rand(numElectrodes, numSamples), rand(numElectrodes, numSamples) }; + + % snippet: doubly-multi-channel + units = types.core.Units('colnames', {}, 'description', 'units'); + units.addDoublyRaggedArray('waveforms', {m1, m2}, ... + 'description', 'multi-channel spike waveforms'); + units.addRaggedArray('electrodes', {[0 1 2], [0 1 2]}, 'table', electrodesTable); + % end snippet + assert(isequal(size(units.waveforms.data), [numSamples, 15])) + assert(isequal(units.waveforms_index.data(:).', uint64([3 6 9 12 15]))) + assert(isequal(units.waveforms_index_index.data(:).', uint64([2 5]))) + + % snippet: helper-functions + [waveforms, waveformsIndex, waveformsIndexIndex] = ... + util.create_doubly_indexed_column({unit1, unit2}, 'spike waveforms'); + + units = types.core.Units( ... + 'colnames', {'waveforms'}, ... + 'description', 'units', ... + 'waveforms', waveforms, ... + 'waveforms_index', waveformsIndex, ... + 'waveforms_index_index', waveformsIndexIndex, ... + 'id', types.hdmf_common.ElementIdentifiers('data', [0; 1])); + % end snippet + assert(isequal(units.waveforms_index_index.data(:).', uint64([3 7]))) +end + +function electrodesTable = local_electrodesTable() + % Minimal electrodes table for the DynamicTableRegion examples. + device = types.core.Device(); + group = types.core.ElectrodeGroup( ... + 'description', 'example group', 'location', 'cortex', ... + 'device', types.untyped.SoftLink(device)); + groupView = types.untyped.ObjectView(group); + electrodesTable = types.core.ElectrodesTable( ... + 'colnames', {'location', 'group', 'group_name'}, ... + 'description', 'electrodes', ... + 'location', types.hdmf_common.VectorData('description', 'location', ... + 'data', {'cortex'; 'cortex'; 'cortex'}), ... + 'group', types.hdmf_common.VectorData('description', 'group', ... + 'data', [groupView; groupView; groupView]), ... + 'group_name', types.hdmf_common.VectorData('description', 'group name', ... + 'data', {'g'; 'g'; 'g'}), ... + 'id', types.hdmf_common.ElementIdentifiers('data', (0:2)')); +end diff --git a/docs/source/pages/how_to/index.rst b/docs/source/pages/how_to/index.rst index 97ca668cc..7f06d7a62 100644 --- a/docs/source/pages/how_to/index.rst +++ b/docs/source/pages/how_to/index.rst @@ -4,6 +4,7 @@ Work with Data :maxdepth: 1 working_with_containers + ragged_arrays Use Extensions ============== diff --git a/docs/source/pages/how_to/ragged_arrays.rst b/docs/source/pages/how_to/ragged_arrays.rst new file mode 100644 index 000000000..0916a49e0 --- /dev/null +++ b/docs/source/pages/how_to/ragged_arrays.rst @@ -0,0 +1,164 @@ +.. _how_to_ragged_arrays: + +Storing Ragged and Doubly-Ragged Array Columns +=============================================== + +This guide shows you how to store columns of a :class:`types.hdmf_common.DynamicTable` +that hold a *variable* number of items per row. + +.. contents:: + :local: + :depth: 2 + +Overview +-------- +Most :class:`types.hdmf_common.DynamicTable` columns have exactly one value per row. +Some columns instead need a *variable* number of values per row. NWB stores these as +**ragged arrays**: + +- A **ragged array** stores a variable number of elements per row (for example, the + spike times of each unit). It is backed by a :class:`types.hdmf_common.VectorData` + column plus a companion :class:`types.hdmf_common.VectorIndex` + (named ``_index``) that marks each row's boundary. +- A **doubly-ragged array** adds a second level of variability: each row holds a + variable number of sub-groups, and each sub-group holds a variable number of + fixed-length elements. The canonical example is the :class:`types.core.Units` + ``waveforms`` column: per unit, a variable number of spike events, each with one + waveform per recording electrode. It is backed by a + :class:`types.hdmf_common.VectorData` column plus two + :class:`types.hdmf_common.VectorIndex` levels (``_index`` over sub-groups + and ``_index_index`` over rows). + +For a full description of how NWB represents these on disk — including a diagram of the +doubly-ragged layout — see the "Tables and ragged arrays" and "Doubly ragged arrays" +sections of the `NWB format specification +`_. + +MatNWB provides two :class:`types.hdmf_common.DynamicTable` methods that build and wire +these objects for you in a single call: + +- ``addRaggedArray`` for ragged columns. +- ``addDoublyRaggedArray`` for doubly-ragged columns. + +.. note:: + + ``addRaggedArray`` and ``addDoublyRaggedArray`` build a whole column in one call. A + ragged (single-index) column can also be filled row by row with ``addRow``, but a + **doubly-ragged** column cannot be built reliably that way — ``addRow`` infers each + value's structure from its array shape, so it may index the data incorrectly. Use + ``addDoublyRaggedArray`` for those. + +Ragged arrays +------------- +Pass a cell array with one cell per row; each cell holds that row's elements. For +example, to store the spike times of two units in the :class:`types.core.Units` table: + +.. literalinclude:: examples/ragged_arrays_examples.m + :language: matlab + :start-after: % snippet: ragged-spike-times + :end-before: % end snippet + :dedent: + +The ``spike_times`` column now holds all five values, and ``spike_times_index`` is +``[3 5]`` — marking that the first unit owns values 1-3 and the second owns values 4-5. + +Referencing another table +~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Provide the ``table`` argument to store a +:class:`types.hdmf_common.DynamicTableRegion` (a column of row references) instead of a +plain :class:`types.hdmf_common.VectorData`. For example, to record which electrodes +each unit was detected on (``electrodesTable`` is an existing electrodes table): + +.. literalinclude:: examples/ragged_arrays_examples.m + :language: matlab + :start-after: % snippet: ragged-electrodes-region + :end-before: % end snippet + :dedent: + +The values are 0-based row indices into the referenced table. + +Doubly-ragged arrays +-------------------- +Use ``addDoublyRaggedArray`` for columns such as :class:`types.core.Units` +``waveforms``. There are two input forms depending on how many electrodes contribute a +waveform per spike. + +.. note:: + + Provide each waveform in natural "one row per waveform, columns are samples" order — + the same order the schema uses (``waveforms`` has dimensions + ``[num_waveforms, num_samples]``). The method transposes and stores the data in the + layout NWB expects, so you do **not** apply MatNWB's usual reversed-dimension + convention here. + +Single electrode (shortcut form) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +When each spike has a single waveform (one electrode per unit), pass one +``[numWaveforms x numSamples]`` matrix per unit — one row per waveform, which for a +single electrode is one row per spike (here ``unit1`` is ``[3 x 40]`` and ``unit2`` is +``[4 x 40]``): + +.. literalinclude:: examples/ragged_arrays_examples.m + :language: matlab + :start-after: % snippet: doubly-single-electrode + :end-before: % end snippet + :dedent: + +This yields ``waveforms_index = [1 2 3 4 5 6 7]`` (one waveform per spike) and +``waveforms_index_index = [3 7]`` (3 spikes, then 4 spikes). + +Multiple channels (nested form) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +When each spike is recorded across several electrodes, each spike has one waveform *per +electrode*. Use a nested cell array where ``data{unit}{spike}`` is a +``[numWaveforms x numSamples]`` matrix — one row per waveform, which here is one per +electrode. Here unit 1 has 2 spikes and unit 2 has 3 spikes, each recorded on 3 +electrodes (``m1`` and ``m2`` are the two units' cell arrays of ``[3 x 40]`` per-spike +matrices): + +.. literalinclude:: examples/ragged_arrays_examples.m + :language: matlab + :start-after: % snippet: doubly-multi-channel + :end-before: % end snippet + :dedent: + +This yields ``waveforms.data`` of size ``[40 15]`` (``numSamples`` × ``numWaveforms``, +where ``numWaveforms`` = 5 spikes × 3 electrodes = 15), ``waveforms_index = [3 6 9 12 +15]`` (3 electrodes per spike), and ``waveforms_index_index = [2 5]`` (2 spikes, then 3 +spikes). The ``electrodes`` column is paired in the same order as the waveform rows within +each spike. + +.. warning:: + + For a multi-channel unit, the order of the waveform rows within each spike must match + the order of the electrodes listed in that unit's ``electrodes`` row, and each spike + of a given unit should have the same number of electrodes. + +Understanding the two index levels +---------------------------------- +For the multi-channel example above: + +- ``waveforms.data`` (``[numSamples x numWaveforms]`` = ``[40 x 15]``) holds all 15 + individual waveforms concatenated, with samples down the rows. +- ``waveforms_index`` (``[3 6 9 12 15]``) has one entry per spike event and marks where + each spike's waveforms end, so spike 1 owns waveforms 1-3, spike 2 owns 4-6, and so on. +- ``waveforms_index_index`` (``[2 5]``) has one entry per unit and marks where each + unit's spike events end, so unit 1 owns spikes 1-2 and unit 2 owns spikes 3-5. + +Building columns without adding them to a table +----------------------------------------------- +If you need the underlying objects (for example, to pass them to a constructor), use the +helper functions that ``addRaggedArray`` and ``addDoublyRaggedArray`` build on: + +- ``util.create_indexed_column`` returns a + :class:`types.hdmf_common.VectorData` (or :class:`types.hdmf_common.DynamicTableRegion`) + and its :class:`types.hdmf_common.VectorIndex`. +- ``util.create_doubly_indexed_column`` returns a + :class:`types.hdmf_common.VectorData` and two :class:`types.hdmf_common.VectorIndex` + levels. + +.. literalinclude:: examples/ragged_arrays_examples.m + :language: matlab + :start-after: % snippet: helper-functions + :end-before: % end snippet + :dedent: