diff --git a/+io/+backend/+hdf5/HDF5Reader.m b/+io/+backend/+hdf5/HDF5Reader.m index 9381841b..06def35f 100644 --- a/+io/+backend/+hdf5/HDF5Reader.m +++ b/+io/+backend/+hdf5/HDF5Reader.m @@ -3,6 +3,21 @@ % % This reader is intentionally thin and delegates to the existing HDF5 % utility functions used by matnwb today. + % + % A reader keeps the tree returned by readRootInfo and the object + % addresses used to resolve references (see + % io.backend.hdf5.ReferenceTargetResolver) for as long as it exists. A + % reader should therefore be used for a single read, and not across + % writes to the file. + + properties (Access = private) + % RootInfo - h5info tree of the whole file, set by readRootInfo. + RootInfo = [] + + % ReferenceResolver - Resolves reference targets for this file. + % Created on the first reference read. + ReferenceResolver = [] + end methods function obj = HDF5Reader(filename) @@ -19,6 +34,7 @@ function node = readRootInfo(obj) node = h5info(obj.Filename); + obj.RootInfo = node; end function tf = isReferenceDataset(~, datasetInfo) @@ -98,7 +114,8 @@ fid = H5F.open(obj.Filename, 'H5F_ACC_RDONLY', 'H5P_DEFAULT'); aid = H5A.open_by_name(fid, context, attributeInfo.Name); tid = H5A.get_type(aid); - attributeValue = io.parseReference(aid, tid, attributeInfo.Value); + attributeValue = io.parseReference(aid, tid, attributeInfo.Value, ... + obj.getReferenceResolver()); H5T.close(tid); H5A.close(aid); H5F.close(fid); @@ -134,7 +151,8 @@ % Load all H5T references. This is required, unfortunately also a % bottleneck tid = H5D.get_type(did); - datasetValue = io.parseReference(did, tid, H5D.read(did)); + datasetValue = io.parseReference(did, tid, H5D.read(did), ... + obj.getReferenceResolver()); H5T.close(tid); elseif strcmp(dataspace.Type, 'scalar') datasetValue = H5D.read(did); @@ -163,7 +181,8 @@ end case 'H5T_COMPOUND' isScalar = true; - datasetValue = io.parseCompound(did, datasetValue, isScalar); + datasetValue = io.parseCompound(did, datasetValue, isScalar, ... + obj.getReferenceResolver()); end else % non scalar sid = H5D.get_space(did); @@ -196,4 +215,22 @@ end end end + + methods (Access = private) + function resolver = getReferenceResolver(obj) + if isempty(obj.ReferenceResolver) + obj.ReferenceResolver = io.backend.hdf5.ReferenceTargetResolver( ... + @() obj.listObjectPaths()); + end + resolver = obj.ReferenceResolver; + end + + function objectPaths = listObjectPaths(obj) + rootInfo = obj.RootInfo; + if isempty(rootInfo) + rootInfo = h5info(obj.Filename); + end + objectPaths = io.backend.hdf5.ReferenceTargetResolver.listObjectPaths(rootInfo); + end + end end diff --git a/+io/+backend/+hdf5/ReferenceTargetResolver.m b/+io/+backend/+hdf5/ReferenceTargetResolver.m new file mode 100644 index 00000000..3681a86b --- /dev/null +++ b/+io/+backend/+hdf5/ReferenceTargetResolver.m @@ -0,0 +1,204 @@ +classdef ReferenceTargetResolver < handle +% ReferenceTargetResolver - Resolve HDF5 references to the paths of their targets. +% +% H5R.get_name has no path to work from when an object was reached through +% a reference, so the HDF5 library searches the whole file for it on every +% call. Resolving each reference that way makes reading a file take time +% proportional to the number of references times the number of objects, +% which dominates reading files with many references (for example one +% table per ROI, each with a VectorIndex referencing its target). +% +% This class records the address of every object in the file once, on the +% first lookup, and resolves a reference by dereferencing it and looking up +% the address of the object it points to. A reference that cannot be +% resolved this way (a null reference, a target that is not in the recorded +% paths, or an HDF5 library that does not report object addresses) is +% resolved with H5R.get_name. +% +% Usage: +% resolver = io.backend.hdf5.ReferenceTargetResolver(objectPaths); +% resolver = io.backend.hdf5.ReferenceTargetResolver(@() listPaths()); +% targetPath = resolver.resolve(locationId, referenceType, rawReference); +% +% The recorded addresses belong to the file at the time of the first +% lookup, so a resolver should only be used while reading one file, and +% not across writes to it. + + properties (Access = private) + % ObjectPaths - Paths of all objects in the file, in the order HDF5 + % visits them (see listObjectPaths), or a function returning them. + ObjectPaths + + % AddressToPath - Map from an object address (as a character key) + % to the path the object was first found at. + AddressToPath + + IsAddressMapBuilt (1,1) logical = false + IsLookupAvailable (1,1) logical = true + end + + methods + function obj = ReferenceTargetResolver(objectPaths) + % ReferenceTargetResolver - Create a resolver for a list of object paths. + % + % Input Arguments: + % - objectPaths (cell | function_handle) - Paths of the objects + % that references may point to, typically from listObjectPaths. + % A function returning the paths defers gathering them until a + % reference is first resolved. + arguments + objectPaths {matnwb.common.compatibility.mustBeA(objectPaths, ["cell", "function_handle"])} + end + obj.ObjectPaths = objectPaths; + end + + function targetPath = resolve(obj, locationId, referenceType, rawReference) + % resolve - Return the path of the object a reference points to. + % + % Input Arguments: + % - locationId - HDF5 identifier of the dataset or attribute that + % holds the reference. + % - referenceType - H5R_OBJECT or H5R_DATASET_REGION. + % - rawReference - Raw reference buffer, as read from the file. + % + % Output Arguments: + % - targetPath (char) - The same path H5R.get_name returns. + targetPath = ''; + % A null reference is all zeros and has no target to look up. + if obj.IsLookupAvailable && any(rawReference(:)) + targetPath = obj.lookup(locationId, referenceType, rawReference); + end + if isempty(targetPath) + targetPath = H5R.get_name(locationId, referenceType, rawReference); + end + end + end + + methods (Static) + function objectPaths = listObjectPaths(groupInfo) + % listObjectPaths - List the paths of all groups and datasets in an h5info tree. + % + % Paths are listed depth first with siblings in name order, which + % is the order H5R.get_name searches in. An object reachable by + % more than one hard link is therefore resolved to the same path + % H5R.get_name would return. + % + % Input Arguments: + % - groupInfo (struct) - Group information as returned by h5info. + % + % Output Arguments: + % - objectPaths (cell) - Row of absolute object paths, starting + % with the path of groupInfo itself. + objectPaths = {groupInfo.Name}; + + numDatasets = numel(groupInfo.Datasets); + numGroups = numel(groupInfo.Groups); + childPaths = cell(1, numDatasets + numGroups); + for iDataset = 1:numDatasets + childPaths{iDataset} = joinPath( ... + groupInfo.Name, groupInfo.Datasets(iDataset).Name); + end + for iGroup = 1:numGroups + % h5info names groups by their full path. + childPaths{numDatasets + iGroup} = groupInfo.Groups(iGroup).Name; + end + % Siblings share the parent path, so sorting full paths sorts + % them by name. + [~, childOrder] = sort(childPaths); + + for iChild = childOrder + if iChild <= numDatasets + objectPaths{end+1} = childPaths{iChild}; %#ok + else + subgroupPaths = io.backend.hdf5.ReferenceTargetResolver.listObjectPaths( ... + groupInfo.Groups(iChild - numDatasets)); + objectPaths = [objectPaths, subgroupPaths]; %#ok + end + end + end + end + + methods (Access = private) + function targetPath = lookup(obj, locationId, referenceType, rawReference) + targetPath = ''; + if ~obj.IsAddressMapBuilt + obj.buildAddressMap(locationId); + end + if ~obj.IsLookupAvailable + return + end + + try + objectId = H5R.dereference(locationId, referenceType, rawReference); + objectCleanup = onCleanup(@() H5O.close(objectId)); %#ok + key = objectKey(H5O.get_info(objectId)); + catch + return + end + + if isKey(obj.AddressToPath, key) + targetPath = obj.AddressToPath(key); + end + end + + function buildAddressMap(obj, locationId) + obj.IsAddressMapBuilt = true; + obj.AddressToPath = containers.Map('KeyType', 'char', 'ValueType', 'char'); + + try + objectPaths = obj.ObjectPaths; + if isa(objectPaths, 'function_handle') + objectPaths = objectPaths(); + end + + fileId = H5I.get_file_id(locationId); + fileCleanup = onCleanup(@() H5F.close(fileId)); %#ok + + for iPath = 1:numel(objectPaths) + objectPath = objectPaths{iPath}; + objectId = H5O.open(fileId, objectPath, 'H5P_DEFAULT'); + try + key = objectKey(H5O.get_info(objectId)); + catch ME + % An object left open would keep the file open. + H5O.close(objectId); + rethrow(ME) + end + H5O.close(objectId); + + if ~isKey(obj.AddressToPath, key) + obj.AddressToPath(key) = objectPath; + end + end + catch + % Any failure while recording addresses, including object + % information without an address or token, leaves every + % reference to be resolved by H5R.get_name. + obj.IsLookupAvailable = false; + end + end + end +end + +function key = objectKey(objectInfo) +% objectKey - Key identifying an object within its file. +% +% HDF5 1.10 reports an object's address, and HDF5 1.12 and later an opaque +% token in its place. Both identify the object uniquely within the file. + if isfield(objectInfo, 'addr') + key = sprintf('%d', objectInfo.addr); + elseif isfield(objectInfo, 'token') && isnumeric(objectInfo.token) + key = sprintf('%d,', objectInfo.token); + else + error('NWB:ReferenceTargetResolver:NoObjectKey', ... + 'The HDF5 library does not report object addresses or tokens.') + end +end + +function fullPath = joinPath(parentPath, name) + if strcmp(parentPath, '/') + fullPath = ['/' name]; + else + fullPath = [parentPath '/' name]; + end +end diff --git a/+io/parseCompound.m b/+io/parseCompound.m index 6fb6427d..52ce957a 100644 --- a/+io/parseCompound.m +++ b/+io/parseCompound.m @@ -1,7 +1,9 @@ -function data = parseCompound(datasetId, data, isScalar) +function data = parseCompound(datasetId, data, isScalar, targetResolver) %did is the dataset_id for the containing dataset %data should be a scalar struct with fields as columns + %targetResolver (optional) is passed on to io.parseReference if nargin < 3; isScalar = false; end + if nargin < 4; targetResolver = []; end typeId = H5D.get_type(datasetId); if isempty(data) % A dataset holding no rows is read back as a 0x0 struct without any @@ -48,7 +50,7 @@ name = referenceFieldName{iFieldName}; rawReference = data.(name); rawTypeId = referenceTypeId{iFieldName}; - data.(name) = io.parseReference(datasetId, rawTypeId, rawReference); + data.(name) = io.parseReference(datasetId, rawTypeId, rawReference, targetResolver); end % Close type ids diff --git a/+io/parseReference.m b/+io/parseReference.m index 2043450e..95d529b3 100644 --- a/+io/parseReference.m +++ b/+io/parseReference.m @@ -1,4 +1,10 @@ -function Reference = parseReference(datasetId, typeId, data) +function Reference = parseReference(datasetId, typeId, data, targetResolver) + % targetResolver (optional) resolves the path each reference points to. + % Without one, every reference is resolved with H5R.get_name, which + % searches the whole file. See io.backend.hdf5.ReferenceTargetResolver. + if nargin < 4 + targetResolver = []; + end referenceSize = size(data); %first dimension is always the raw buffer size referenceSize = referenceSize(2:end); @@ -12,13 +18,18 @@ referenceType = H5ML.get_constant_value('H5R_DATASET_REGION'); end for iReference = 1:totalNumReferences - Reference(iReference) = parseSingleReference(datasetId, referenceType, data(:,iReference)); + Reference(iReference) = parseSingleReference( ... + datasetId, referenceType, data(:,iReference), targetResolver); end Reference = reshape(Reference, referenceSize); end -function Reference = parseSingleReference(datasetId, referenceType, data) - target = H5R.get_name(datasetId, referenceType, data); +function Reference = parseSingleReference(datasetId, referenceType, data, targetResolver) + if isempty(targetResolver) + target = H5R.get_name(datasetId, referenceType, data); + else + target = targetResolver.resolve(datasetId, referenceType, data); + end %% H5R_OBJECT if referenceType == H5ML.get_constant_value('H5R_OBJECT') diff --git a/+tests/+unit/+io/+backend/ReferenceTargetResolverTest.m b/+tests/+unit/+io/+backend/ReferenceTargetResolverTest.m new file mode 100644 index 00000000..0f4d6cd6 --- /dev/null +++ b/+tests/+unit/+io/+backend/ReferenceTargetResolverTest.m @@ -0,0 +1,208 @@ +classdef ReferenceTargetResolverTest < matlab.unittest.TestCase +% ReferenceTargetResolverTest - Tests for io.backend.hdf5.ReferenceTargetResolver. +% +% The resolver replaces a per-reference H5R.get_name call with an address +% lookup, so these tests check that it returns exactly what H5R.get_name +% returns, for dataset and attribute references alike. + + properties (Constant) + FileName = "reference-resolver-test.nwb" + NumPlaneSegmentations = 3 + ElectrodesGroupPath = '/general/extracellular_ephys/electrodes/group' + end + + methods (TestClassSetup) + function createTestFile(testCase) + import matlab.unittest.fixtures.WorkingFolderFixture + testCase.applyFixture(WorkingFolderFixture); + + nwb = tests.factory.NWBFile(); + + % A "group" column of object references in a dataset. + tests.factory.ElectrodeTable(nwb); + + % One VectorIndex per plane segmentation, each referencing its + % target in an attribute, and a DynamicTableRegion referencing + % a plane segmentation. + device = types.core.Device(); + nwb.general_devices.set('Microscope', device); + imagingPlane = tests.factory.ImagingPlane(device); + nwb.general_optophysiology.set('ImagingPlane', imagingPlane); + + imageSegmentation = types.core.ImageSegmentation(); + for iTable = 1:testCase.NumPlaneSegmentations + planeSegmentation = tests.factory.PlaneSegmentation(imagingPlane, ... + 'RoiType', 'pixel_mask', 'NumRois', 2, 'ImageShape', [20, 20]); + imageSegmentation.planesegmentation.set( ... + sprintf('PlaneSegmentation%d', iTable), planeSegmentation); + end + fluorescence = types.core.Fluorescence(); + fluorescence.roiresponseseries.set('RoiResponseSeries', ... + tests.factory.RoiResponseSeries(planeSegmentation, 'NumTimepoints', 5)); + + ophysModule = types.core.ProcessingModule('description', 'ophys'); + ophysModule.nwbdatainterface.set('ImageSegmentation', imageSegmentation); + ophysModule.nwbdatainterface.set('Fluorescence', fluorescence); + nwb.processing.set('ophys', ophysModule); + + nwbExport(nwb, testCase.FileName); + end + end + + methods (Test) + function resolveMatchesGetNameForDatasetReferences(testCase) + fileId = H5F.open(testCase.FileName, 'H5F_ACC_RDONLY', 'H5P_DEFAULT'); + fileCleanup = onCleanup(@() H5F.close(fileId)); %#ok + datasetId = H5D.open(fileId, testCase.ElectrodesGroupPath); + datasetCleanup = onCleanup(@() H5D.close(datasetId)); %#ok + rawReferences = H5D.read(datasetId); + + resolver = testCase.createResolver(); + referenceType = H5ML.get_constant_value('H5R_OBJECT'); + for iReference = 1:size(rawReferences, 2) + rawReference = rawReferences(:, iReference); + testCase.verifyEqual( ... + resolver.resolve(datasetId, referenceType, rawReference), ... + H5R.get_name(datasetId, referenceType, rawReference)); + end + end + + function resolveMatchesGetNameForAttributeReferences(testCase) + fileId = H5F.open(testCase.FileName, 'H5F_ACC_RDONLY', 'H5P_DEFAULT'); + fileCleanup = onCleanup(@() H5F.close(fileId)); %#ok + + resolver = testCase.createResolver(); + referenceType = H5ML.get_constant_value('H5R_OBJECT'); + for iTable = 1:testCase.NumPlaneSegmentations + indexPath = sprintf( ... + '/processing/ophys/ImageSegmentation/PlaneSegmentation%d/pixel_mask_index', iTable); + attributeId = H5A.open_by_name(fileId, indexPath, 'target'); + rawReference = H5A.read(attributeId, 'H5ML_DEFAULT'); + + expectedPath = H5R.get_name(attributeId, referenceType, rawReference); + actualPath = resolver.resolve(attributeId, referenceType, rawReference); + H5A.close(attributeId); + + testCase.verifyEqual(actualPath, expectedPath); + testCase.verifyEqual(actualPath, sprintf( ... + '/processing/ophys/ImageSegmentation/PlaneSegmentation%d/pixel_mask', iTable)); + end + end + + function resolveFallsBackForTargetsWithoutRecordedPath(testCase) + % A resolver that knows no paths has to give the same answer + % by falling back to H5R.get_name. + fileId = H5F.open(testCase.FileName, 'H5F_ACC_RDONLY', 'H5P_DEFAULT'); + fileCleanup = onCleanup(@() H5F.close(fileId)); %#ok + datasetId = H5D.open(fileId, testCase.ElectrodesGroupPath); + datasetCleanup = onCleanup(@() H5D.close(datasetId)); %#ok + rawReferences = H5D.read(datasetId); + + resolver = io.backend.hdf5.ReferenceTargetResolver({}); + referenceType = H5ML.get_constant_value('H5R_OBJECT'); + testCase.verifyEqual( ... + resolver.resolve(datasetId, referenceType, rawReferences(:, 1)), ... + H5R.get_name(datasetId, referenceType, rawReferences(:, 1))); + end + + function resolveUsesRecordedPathsInsteadOfGetName(testCase) + % The target is reachable as both /a_alias and /z_target. + % H5R.get_name returns /a_alias, the first path in name order, + % while the address lookup returns the only path it recorded. + % Returning /z_target therefore proves the lookup ran instead + % of the H5R.get_name fallback. + fileName = 'hard-linked-reference.h5'; + createHardLinkedReferenceFile(fileName); + + fileId = H5F.open(fileName, 'H5F_ACC_RDONLY', 'H5P_DEFAULT'); + fileCleanup = onCleanup(@() H5F.close(fileId)); %#ok + datasetId = H5D.open(fileId, '/refs'); + datasetCleanup = onCleanup(@() H5D.close(datasetId)); %#ok + rawReference = H5D.read(datasetId); + + referenceType = H5ML.get_constant_value('H5R_OBJECT'); + testCase.assumeEqual( ... + H5R.get_name(datasetId, referenceType, rawReference), '/a_alias', ... + 'H5R.get_name must return the alias for this test to tell the two paths apart.'); + + resolver = io.backend.hdf5.ReferenceTargetResolver({'/', '/z_target'}); + testCase.verifyEqual( ... + resolver.resolve(datasetId, referenceType, rawReference), '/z_target'); + end + + function objectPathsAreGatheredOnFirstResolve(testCase) + fileId = H5F.open(testCase.FileName, 'H5F_ACC_RDONLY', 'H5P_DEFAULT'); + fileCleanup = onCleanup(@() H5F.close(fileId)); %#ok + datasetId = H5D.open(fileId, testCase.ElectrodesGroupPath); + datasetCleanup = onCleanup(@() H5D.close(datasetId)); %#ok + rawReferences = H5D.read(datasetId); + + callCounter = containers.Map({'count'}, {0}); + filename = char(testCase.FileName); + resolver = io.backend.hdf5.ReferenceTargetResolver( ... + @() countedListObjectPaths(filename, callCounter)); + testCase.verifyEqual(callCounter('count'), 0) + + referenceType = H5ML.get_constant_value('H5R_OBJECT'); + resolver.resolve(datasetId, referenceType, rawReferences(:, 1)); + resolver.resolve(datasetId, referenceType, rawReferences(:, 1)); + testCase.verifyEqual(callCounter('count'), 1) + end + + function listObjectPathsIsDepthFirstInNameOrder(testCase) + datasetB = struct('Name', 'b'); + groupA = struct('Name', '/a', 'Groups', [], 'Datasets', struct('Name', 'x')); + groupC = struct('Name', '/c', 'Groups', [], 'Datasets', []); + rootInfo = struct('Name', '/', 'Groups', [groupC, groupA], 'Datasets', datasetB); + + objectPaths = io.backend.hdf5.ReferenceTargetResolver.listObjectPaths(rootInfo); + + testCase.verifyEqual(objectPaths, {'/', '/a', '/a/x', '/b', '/c'}); + end + + function readerResolvesReferencesOnRead(testCase) + nwb = nwbRead(testCase.FileName, 'ignorecache'); + + planeSegmentation = nwb.processing.get('ophys') ... + .nwbdatainterface.get('ImageSegmentation') ... + .planesegmentation.get('PlaneSegmentation1'); + testCase.verifyEqual(planeSegmentation.pixel_mask_index.target.path, ... + '/processing/ophys/ImageSegmentation/PlaneSegmentation1/pixel_mask'); + end + end + + methods (Access = private) + function resolver = createResolver(testCase) + objectPaths = io.backend.hdf5.ReferenceTargetResolver.listObjectPaths( ... + h5info(testCase.FileName)); + resolver = io.backend.hdf5.ReferenceTargetResolver(objectPaths); + end + end +end + +function objectPaths = countedListObjectPaths(filename, callCounter) + callCounter('count') = callCounter('count') + 1; + objectPaths = io.backend.hdf5.ReferenceTargetResolver.listObjectPaths(h5info(filename)); +end + +function createHardLinkedReferenceFile(fileName) +% createHardLinkedReferenceFile - Write a dataset linked as /z_target and +% /a_alias, and a dataset /refs holding one object reference to /z_target. + fileId = H5F.create(fileName, 'H5F_ACC_TRUNC', 'H5P_DEFAULT', 'H5P_DEFAULT'); + fileCleanup = onCleanup(@() H5F.close(fileId)); %#ok + + targetSpaceId = H5S.create_simple(1, 3, []); + targetSpaceCleanup = onCleanup(@() H5S.close(targetSpaceId)); %#ok + targetId = H5D.create(fileId, '/z_target', 'H5T_NATIVE_DOUBLE', targetSpaceId, 'H5P_DEFAULT'); + targetCleanup = onCleanup(@() H5D.close(targetId)); %#ok + H5D.write(targetId, 'H5ML_DEFAULT', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', [1, 2, 3]); + + H5L.create_hard(fileId, '/z_target', fileId, '/a_alias', 'H5P_DEFAULT', 'H5P_DEFAULT'); + + referenceSpaceId = H5S.create_simple(1, 1, []); + referenceSpaceCleanup = onCleanup(@() H5S.close(referenceSpaceId)); %#ok + referencesId = H5D.create(fileId, '/refs', 'H5T_STD_REF_OBJ', referenceSpaceId, 'H5P_DEFAULT'); + referencesCleanup = onCleanup(@() H5D.close(referencesId)); %#ok + rawReference = H5R.create(fileId, '/z_target', 'H5R_OBJECT', -1); + H5D.write(referencesId, 'H5ML_DEFAULT', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', rawReference); +end