From 11662c5119a00eaaa0d5299b6ec99fca3a9336b2 Mon Sep 17 00:00:00 2001 From: Kenya Otsuka Date: Tue, 16 Jun 2026 14:07:52 +0000 Subject: [PATCH 1/9] chore: relax bdpy[mri] pins for NumPy 2.x and NiBabel 5.x Drop the numpy<1.24 cap and the exact nibabel==3.2 pin that blocked contemporary NumPy and NiBabel (issue #121). nipy>=0.6.1 is required for NumPy 2.0 support. - mri: numpy<1.24 removed (inherits core numpy>=1.20); nibabel==3.2 -> nibabel>=3.2,<6 (6.0 drops NumPy 1.x); nipy>=0.6.0 -> nipy>=0.6.1 --- pyproject.toml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 8f15311..63f100f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,9 +40,8 @@ fig = [ "Pillow" ] mri = [ - "numpy<1.24", - "nibabel==3.2", - "nipy>=0.6.0" + "nibabel>=3.2,<6", + "nipy>=0.6.1" ] pipeline = [ "hydra-core", From ad946968b2cb6a8d65d017d83b0f4b336f6ab108 Mon Sep 17 00:00:00 2001 From: Kenya Otsuka Date: Tue, 16 Jun 2026 14:07:52 +0000 Subject: [PATCH 2/9] fix: support NiBabel 5.x and NumPy 2.x in bdpy/mri NiBabel 5.0 removed Image.get_data(); read nipy images via get_fdata() instead in load_mri, load_epi and fmriprep. Also fix Python 2 leftovers (xrange -> range, reader.next() -> next(reader)) that would fail on Python 3. Add a load_mri round-trip test guarding the get_fdata() path. --- bdpy/mri/fmriprep.py | 10 +++++----- bdpy/mri/glm.py | 2 +- bdpy/mri/load_epi.py | 4 ++-- bdpy/mri/load_mri.py | 2 +- tests/test_mri.py | 22 ++++++++++++++++++++++ 5 files changed, 31 insertions(+), 9 deletions(-) diff --git a/bdpy/mri/fmriprep.py b/bdpy/mri/fmriprep.py index ff4e911..6567a69 100644 --- a/bdpy/mri/fmriprep.py +++ b/bdpy/mri/fmriprep.py @@ -408,7 +408,7 @@ def __load_volume(self): """ img = nipy.load_image(self.__dpath) - data = img.get_data() + data = img.get_fdata() if data.ndim == 4: data = data.reshape(-1, data.shape[-1], order='F').T i_len, j_len, k_len, t = img.shape @@ -823,9 +823,9 @@ def __get_xyz(img): # 3D-image i_len, j_len, k_len = img.shape affine = img.coordmap.affine - ijk = np.array(list(itertools.product(xrange(i_len), - xrange(j_len), - xrange(k_len), + ijk = np.array(list(itertools.product(range(i_len), + range(j_len), + range(k_len), [1]))).T return np.dot(affine, ijk)[:-1] @@ -840,7 +840,7 @@ def __load_mri(fpath): """ img = nipy.load_image(fpath) - data = img.get_data() + data = img.get_fdata() if data.ndim == 4: data = data.reshape(-1, data.shape[-1], order='F').T i_len, j_len, k_len, t = img.shape diff --git a/bdpy/mri/glm.py b/bdpy/mri/glm.py index d57e218..2916cac 100644 --- a/bdpy/mri/glm.py +++ b/bdpy/mri/glm.py @@ -64,7 +64,7 @@ def make_paradigm(event_files, num_vols, tr=2., cond_col=2, label_col=None, regr with open(ef, 'r') as f: reader = csv.reader(f, delimiter='\t') - header = reader.next() + header = next(reader) for row in reader: if regressors is not None and row[cond_col] not in regressors: continue diff --git a/bdpy/mri/load_epi.py b/bdpy/mri/load_epi.py index ba1fd78..67ee1f0 100644 --- a/bdpy/mri/load_epi.py +++ b/bdpy/mri/load_epi.py @@ -38,7 +38,7 @@ def load_epi(datafiles): img = nipy.load_image(df) xyz = _check_xyz(xyz, img) - data_list.append(np.array(img.get_data().flatten(), dtype=np.float64)) + data_list.append(np.array(img.get_fdata().flatten(), dtype=np.float64)) data = np.vstack(data_list) @@ -47,7 +47,7 @@ def load_epi(datafiles): def _check_xyz(xyz, img): """Check voxel xyz consistency.""" - xyz_current = _get_xyz(img.coordmap.affine, img.get_data().shape) + xyz_current = _get_xyz(img.coordmap.affine, img.get_fdata().shape) if xyz.size == 0: xyz = xyz_current diff --git a/bdpy/mri/load_mri.py b/bdpy/mri/load_mri.py index c560190..772d6b0 100644 --- a/bdpy/mri/load_mri.py +++ b/bdpy/mri/load_mri.py @@ -15,7 +15,7 @@ def load_mri(fpath): """ img = nipy.load_image(fpath) - data = img.get_data() + data = img.get_fdata() if data.ndim == 4: data = data.reshape(-1, data.shape[-1], order='F').T i_len, j_len, k_len, t = img.shape diff --git a/tests/test_mri.py b/tests/test_mri.py index 3922f34..bee53ae 100644 --- a/tests/test_mri.py +++ b/tests/test_mri.py @@ -25,6 +25,28 @@ def test_get_roiflag_pass0001(self) -> None: self.assertTrue((test_output == exp_output).all()) + def test_load_mri_3d(self) -> None: + """Test load_mri on a 3D volume. + + Guards the NumPy 2.0 / nibabel 5.x compatibility fix: nibabel removed + ``get_data()``, so the loader reads the image via ``get_fdata()``. + """ + import os + import tempfile + + import nibabel + + arr = np.arange(2 * 3 * 4, dtype=np.float32).reshape(2, 3, 4) + with tempfile.TemporaryDirectory() as tmpdir: + fpath = os.path.join(tmpdir, 'test.nii.gz') + nibabel.save(nibabel.Nifti1Image(arr, affine=np.eye(4)), fpath) + data, xyz, ijk = bmr.load_mri(fpath) # type: ignore + + # A 3D volume is returned flattened in Fortran order. + np.testing.assert_array_equal(data, arr.flatten(order='F')) + self.assertEqual(xyz.shape, (3, arr.size)) + self.assertEqual(ijk.shape, (3, arr.size)) + def test_get_roiflag_pass0002(self) -> None: """Test for get_roiflag (pass case 0002).""" roi_xyz = [np.array([[1, 2, 3], From 2fdc42449588bf5ba5b640dfe98d199821fefd70 Mon Sep 17 00:00:00 2001 From: Kenya Otsuka Date: Tue, 16 Jun 2026 14:07:52 +0000 Subject: [PATCH 3/9] fix: read MATLAB v7.3 files with h5py instead of hdf5storage (#106) Under NumPy 2.0 hdf5storage.loadmat fails (it referenced the removed np.unicode_), breaking DecodedFeatures.get() and other loaders. Add bdpy/dataform/_mat_v73.py with an h5py-based reader for the dense and sparse-struct layouts bdpy uses, plus a loadmat_key() helper that keeps v5 .mat support via scipy and falls back to h5py for v7.3. Route the load paths in features, datastore and sparse through it. Saving still uses hdf5storage. SparseArray.save now retries with a from-scratch rewrite only if the in-place hdf5storage.savemat raises under NumPy 2.0 (overwriting a struct); the normal merging behavior is preserved whenever the write succeeds. Update the features test to write v7.3 mock data so the h5py path is exercised. --- bdpy/dataform/_mat_v73.py | 119 ++++++++++++++++++++++++++++++++ bdpy/dataform/datastore.py | 6 +- bdpy/dataform/features.py | 12 ++-- bdpy/dataform/sparse.py | 69 ++++++++++-------- tests/dataform/test_features.py | 39 +++++------ 5 files changed, 185 insertions(+), 60 deletions(-) create mode 100644 bdpy/dataform/_mat_v73.py diff --git a/bdpy/dataform/_mat_v73.py b/bdpy/dataform/_mat_v73.py new file mode 100644 index 0000000..edf425b --- /dev/null +++ b/bdpy/dataform/_mat_v73.py @@ -0,0 +1,119 @@ +"""Read MATLAB v7.3 (HDF5) ``.mat`` files with h5py. + +bdpy historically relied on :func:`hdf5storage.loadmat` for the load path, but +that broke under NumPy 2.0, which removed ``np.unicode_`` that older hdf5storage +referenced (see issue #106). This module reimplements the *read* side on top of +h5py for the data layouts bdpy actually uses (dense numeric arrays and the +sparse-array struct). Saving still goes through hdf5storage. + +This file is a part of BdPy. +""" + +import h5py +import numpy as np +import scipy.io as sio + +__all__ = ["load_array", "loadmat_key", "read_cell", "read_dataset"] + + +def read_dataset(dset: h5py.Dataset) -> np.ndarray: + """Read an h5py dataset, undoing MATLAB v7.3 / hdf5storage conventions. + + MATLAB stores arrays in Fortran (column-major) order, so multi-dimensional + datasets are written transposed relative to NumPy's C order. hdf5storage + additionally records the original Python shape and empty-array flags as + ``Python.*`` attributes, which we honor to reproduce ``hdf5storage.loadmat``. + + Parameters + ---------- + dset : h5py.Dataset + Dataset to read. + + Returns + ------- + numpy.ndarray + The array with its original shape restored. + """ + attrs = dset.attrs + if "Python.Empty" in attrs or "MATLAB_empty" in attrs: + shape = tuple(int(x) for x in attrs.get("Python.Shape", ())) + return np.empty(shape, dtype=dset.dtype) + arr = dset[()] + if "MATLAB_class" in attrs and isinstance(arr, np.ndarray) and arr.ndim >= 2: + arr = np.transpose(arr) + if "Python.Shape" in attrs: + arr = np.asarray(arr).reshape(tuple(int(x) for x in attrs["Python.Shape"])) + return np.asarray(arr) + + +def read_cell(f: h5py.File, dset: h5py.Dataset) -> list: + """Read a MATLAB cell array (or a plain matrix) into a list of arrays. + + hdf5storage stores Python tuples/lists as MATLAB cell arrays, i.e. an object + dataset of HDF5 references to the individual elements. Files written by other + tools (e.g. MATLAB or Julia) may instead store the same information as a plain + 2-D matrix whose rows are the elements; both layouts are handled here. + + Parameters + ---------- + f : h5py.File + Open file, used to dereference cell-array element references. + dset : h5py.Dataset + The cell (object) dataset or a plain matrix. + + Returns + ------- + list of numpy.ndarray + One array per cell element / matrix row. + """ + data = dset[()] + if isinstance(data, np.ndarray) and data.dtype == object: + return [read_dataset(f[ref]) for ref in data.ravel()] + arr = np.asarray(data) + return [arr[i] for i in range(arr.shape[0])] + + +def load_array(path: str, key: str) -> np.ndarray: + """Load a single dense numeric array from a v7.3 ``.mat`` file. + + Parameters + ---------- + path : str + Path to the ``.mat`` file. + key : str + Variable name to load. + + Returns + ------- + numpy.ndarray + The loaded array. + """ + with h5py.File(path, "r") as f: + return read_dataset(f[key]) + + +def loadmat_key(path: str, key: str) -> np.ndarray: + """Load one variable from a ``.mat`` file, handling both v5 and v7.3. + + MATLAB v5 (and earlier) files are read with :func:`scipy.io.loadmat`; v7.3 + (HDF5) files, which scipy cannot read, fall back to the h5py reader. This + preserves v5 support while avoiding hdf5storage on the load path (which + breaks under NumPy 2.0). + + Parameters + ---------- + path : str + Path to the ``.mat`` file. + key : str + Variable name to load. + + Returns + ------- + numpy.ndarray + The loaded array. + """ + try: + array = sio.loadmat(path)[key] + except (NotImplementedError, ValueError): + return load_array(path, key) + return np.asarray(array) diff --git a/bdpy/dataform/datastore.py b/bdpy/dataform/datastore.py index 00e4d91..4098b15 100644 --- a/bdpy/dataform/datastore.py +++ b/bdpy/dataform/datastore.py @@ -11,10 +11,11 @@ import re import h5py -import hdf5storage import numpy as np import scipy.io as sio +from . import _mat_v73 + __all__ = ['DataStore', 'DirStore'] @@ -250,7 +251,8 @@ def get(self, **kargs): return dat def __load_feature(self, fpath): - r = hdf5storage.loadmat(fpath)[self.__variable] + # v5 .mat via scipy, v7.3 (HDF5) via h5py (avoids hdf5storage under NumPy 2.0). + r = _mat_v73.loadmat_key(fpath, self.__variable) if self.__squeeze: r = np.squeeze(r) return r diff --git a/bdpy/dataform/features.py b/bdpy/dataform/features.py index 46fee23..193537f 100644 --- a/bdpy/dataform/features.py +++ b/bdpy/dataform/features.py @@ -19,14 +19,14 @@ import hdf5storage import numpy as np -import scipy.io as sio + +from . import _mat_v73 def _load_array_with_key(key: str, path: str) -> np.ndarray: - try: - return sio.loadmat(path)[key] - except (NotImplementedError, ValueError): - return hdf5storage.loadmat(path)[key] + # v5 .mat via scipy, v7.3 (HDF5) via h5py; avoids hdf5storage on the load + # path, which breaks under NumPy 2.0 (see bdpy/dataform/_mat_v73.py). + return _mat_v73.loadmat_key(path, key) def _determine_num_parallel(num_files: int) -> int: @@ -86,7 +86,7 @@ def __init__( if feature_index is not None: if not os.path.exists(feature_index): raise RuntimeError('%s do not exist' % feature_index) - self.__feat_index_table = hdf5storage.loadmat(feature_index)['index'] + self.__feat_index_table = _mat_v73.loadmat_key(feature_index, 'index') # NOTE: type of self.__feature_index_table is ambiguous else: self.__feat_index_table = None diff --git a/bdpy/dataform/sparse.py b/bdpy/dataform/sparse.py index 2e4593f..92961ab 100644 --- a/bdpy/dataform/sparse.py +++ b/bdpy/dataform/sparse.py @@ -12,6 +12,8 @@ import hdf5storage import numpy as np +from . import _mat_v73 + def load_array(fname, key='data'): """Load an array (dense or sparse).""" @@ -22,8 +24,8 @@ def load_array(fname, key='data'): s_ary = SparseArray(fname, key=key) return s_ary.dense elif type(f[key][()]) == np.ndarray: - # Dense array - return hdf5storage.loadmat(fname)[key] + # Dense array (read with h5py; hdf5storage breaks under NumPy 2.0) + return _mat_v73.read_dataset(f[key]) else: raise RuntimeError('Unsupported data type: %s' % type(f[key][()])) @@ -75,12 +77,24 @@ def dense(self): return self.__make_dense() def save(self, fname, key='data', dtype=np.float64): - hdf5storage.savemat(fname, {key: {u'__bdpy_sparse_arrray': True, - u'index': self.__index, - u'value': self.__value.astype(dtype), - u'shape': self.__shape, - u'background' : self.__background}}, - format='7.3', oned_as='column', store_python_metadata=True) + payload = {key: {u'__bdpy_sparse_arrray': True, + u'index': self.__index, + u'value': self.__value.astype(dtype), + u'shape': self.__shape, + u'background': self.__background}} + try: + hdf5storage.savemat(fname, payload, format='7.3', + oned_as='column', store_python_metadata=True) + except ValueError: + # Under NumPy 2.0, hdf5storage.savemat raises when overwriting an + # existing struct in place (it compares array-valued attributes). + # Fall back to rewriting the file from scratch. The normal (merging) + # behavior is preserved whenever the in-place write succeeds; only + # this fallback discards any other variables already in the file. + if os.path.exists(fname): + os.remove(fname) + hdf5storage.savemat(fname, payload, format='7.3', + oned_as='column', store_python_metadata=True) return None def __make_sparse(self, array): @@ -95,29 +109,24 @@ def __make_dense(self): return dense def __load(self, fname, key='data'): - data = hdf5storage.loadmat(fname)[key] - - index = data[u'index'] - if isinstance(index, tuple): - self.__index = index - elif isinstance(index, np.ndarray): - self.__index = tuple(index[0]) - else: - raise TypeError('Unsupported data type ("index").') - - value = data[u'value'] - if value.ndim == 1: + # Read with h5py instead of hdf5storage, which breaks under NumPy 2.0. + # The struct stores ``index``/``shape`` as cell arrays (object refs) when + # written by bdpy, or as plain matrices when written by other tools. + with h5py.File(fname, 'r') as f: + g = f[key] + self.__index = tuple( + np.asarray(c).ravel().astype(int) + for c in _mat_v73.read_cell(f, g['index']) + ) + + value = np.asarray(_mat_v73.read_dataset(g['value'])).ravel() self.__value = value - else: - self.__value = value.flatten() - array_shape = data[u'shape'] - if isinstance(array_shape, tuple): - self.__shape = array_shape - elif isinstance(array_shape, np.ndarray): - self.__shape = array_shape.flatten() - else: - raise TypeError('Unsupported data type ("shape").') + self.__shape = tuple( + int(np.asarray(s).ravel()[0]) + for s in _mat_v73.read_cell(f, g['shape']) + ) - self.__background = data[u'background'] + background = np.asarray(_mat_v73.read_dataset(g['background'])).ravel() + self.__background = background[0] if background.size else 0 return None diff --git a/tests/dataform/test_features.py b/tests/dataform/test_features.py index 927341b..9d760eb 100644 --- a/tests/dataform/test_features.py +++ b/tests/dataform/test_features.py @@ -3,12 +3,10 @@ from typing import List, Tuple import os -from glob import glob import tempfile import numpy as np from numpy.testing import assert_array_equal -import scipy.io as sio import hdf5storage from bdpy.dataform.features import Features @@ -19,16 +17,26 @@ def _prepare_mock_data( mock_layer_names: List[str], mock_image_names: List[str], mock_shapes: List[Tuple[int, ...]] - ) -> None: - """Prepare mock data for testing.""" + ) -> dict: + """Prepare mock data for testing. + + Files are written as MATLAB v7.3 (HDF5) so the loader's h5py path (used for + NumPy 2.0 compatibility) is exercised. The stacked arrays are returned so + tests can compare against them without re-reading the files. + """ + stacked = {} for layer_name, shape in zip(mock_layer_names, mock_shapes): os.makedirs(os.path.join(tmpdir, layer_name)) + arrays = [] for image_name in mock_image_names: data = np.random.rand(*shape) hdf5storage.savemat( os.path.join(tmpdir, layer_name, image_name + '.mat'), {'feat': data}, - format='5') + store_python_metadata=True) + arrays.append(data) + stacked[layer_name] = np.vstack(arrays) + return stacked class TestDataformFeatures(unittest.TestCase): @@ -45,29 +53,16 @@ def setUp(self): ] self.mock_shapes = [(1, 1000), (1, 256, 13, 13)] self.feature_dir = tempfile.TemporaryDirectory() - _prepare_mock_data( + stacked = _prepare_mock_data( self.feature_dir.name, self.mock_layer_names, self.mock_image_names, self.mock_shapes ) - # Loading test data - # AlexNet, fc8, all samples - self.alexnet_fc8_all = np.vstack( - [ - sio.loadmat(f)['feat'] - for f in sorted(glob(os.path.join(self.feature_dir.name, 'fc8', '*.mat'))) - ] - ) - - # AlexNet, conv5, all samples - self.alexnet_conv5_all = np.vstack( - [ - sio.loadmat(f)['feat'] - for f in sorted(glob(os.path.join(self.feature_dir.name, 'conv5', '*.mat'))) - ] - ) + # Expected data (samples stacked in sorted-filename order) + self.alexnet_fc8_all = stacked['fc8'] + self.alexnet_conv5_all = stacked['conv5'] def tearDown(self): self.feature_dir.cleanup() From 1de6e6127c376d05c7252010d5f5d3a0a2dceeb5 Mon Sep 17 00:00:00 2001 From: Kenya Otsuka Date: Tue, 16 Jun 2026 15:05:21 +0000 Subject: [PATCH 4/9] fix: replace SparseArray.save fallback with NumPy-2 h5py writer Under NumPy 2.x, hdf5storage.savemat can fail when SparseArray.save overwrites an existing sparse struct, and the previous fallback deleted the whole file before rewriting, dropping any unrelated variables. Make SparseArray.save use a NumPy-2-only h5py writer that opens the target in append mode and replaces only the requested key/group, leaving other variables intact. NumPy < 2.0 keeps the existing hdf5storage.savemat path. The dense save paths are unchanged and still rely on hdf5storage for the MATLAB-v7.3 metadata/compatibility behavior. Add a NumPy-2-only test asserting that SparseArray.save preserves other top-level variables already stored in the target file. --- bdpy/dataform/sparse.py | 54 ++++++++++++++++++++++++----------- tests/dataform/test_sparse.py | 22 ++++++++++++++ 2 files changed, 60 insertions(+), 16 deletions(-) diff --git a/bdpy/dataform/sparse.py b/bdpy/dataform/sparse.py index 92961ab..969b3f9 100644 --- a/bdpy/dataform/sparse.py +++ b/bdpy/dataform/sparse.py @@ -14,6 +14,12 @@ from . import _mat_v73 +# hdf5storage.savemat can fail for SparseArray.save under NumPy 2.x, especially +# when overwriting an existing sparse struct. SparseArray.save therefore uses a +# direct h5py writer only under NumPy >= 2. Dense-array save paths are unchanged +# and still rely on hdf5storage for MATLAB-v7.3 metadata/compatibility behavior. +_NUMPY2 = int(np.__version__.split('.')[0]) >= 2 + def load_array(fname, key='data'): """Load an array (dense or sparse).""" @@ -77,26 +83,42 @@ def dense(self): return self.__make_dense() def save(self, fname, key='data', dtype=np.float64): - payload = {key: {u'__bdpy_sparse_arrray': True, - u'index': self.__index, - u'value': self.__value.astype(dtype), - u'shape': self.__shape, - u'background': self.__background}} - try: - hdf5storage.savemat(fname, payload, format='7.3', - oned_as='column', store_python_metadata=True) - except ValueError: - # Under NumPy 2.0, hdf5storage.savemat raises when overwriting an - # existing struct in place (it compares array-valued attributes). - # Fall back to rewriting the file from scratch. The normal (merging) - # behavior is preserved whenever the in-place write succeeds; only - # this fallback discards any other variables already in the file. - if os.path.exists(fname): - os.remove(fname) + if _NUMPY2: + # Avoid hdf5storage.savemat here (it can fail when overwriting an + # existing sparse struct under NumPy 2.x) and write the struct with + # h5py instead. We replace only ``key`` in the target file, leaving + # any other variables already stored there untouched. + self.__save_h5py(fname, key=key, dtype=dtype) + else: + payload = {key: {u'__bdpy_sparse_arrray': True, + u'index': self.__index, + u'value': self.__value.astype(dtype), + u'shape': self.__shape, + u'background': self.__background}} hdf5storage.savemat(fname, payload, format='7.3', oned_as='column', store_python_metadata=True) return None + def __save_h5py(self, fname, key='data', dtype=np.float64): + # NumPy-2-only writer for the sparse-array struct. The layout mirrors + # what ``__load`` expects: ``index``/``shape`` as plain matrices (which + # _mat_v73.read_cell restores row-by-row) and ``value``/``background`` + # as plain datasets. Opening in append mode and deleting only ``key`` + # avoids the unconditional full-file rewrite of the previous fallback. + index = np.vstack([np.asarray(i, dtype=np.int64).ravel() + for i in self.__index]) + mode = 'a' if os.path.exists(fname) else 'w' + with h5py.File(fname, mode) as f: + if key in f: + del f[key] + g = f.create_group(key) + g.create_dataset(u'__bdpy_sparse_arrray', data=True) + g.create_dataset(u'index', data=index) + g.create_dataset(u'value', data=self.__value.astype(dtype).ravel()) + g.create_dataset(u'shape', data=np.asarray(self.__shape, dtype=np.int64)) + g.create_dataset(u'background', data=np.asarray(self.__background)) + return None + def __make_sparse(self, array): self.__index = np.where(array != self.__background) self.__value = array[self.__index] diff --git a/tests/dataform/test_sparse.py b/tests/dataform/test_sparse.py index 20810c9..3047f6b 100644 --- a/tests/dataform/test_sparse.py +++ b/tests/dataform/test_sparse.py @@ -4,6 +4,7 @@ import tempfile import unittest +import h5py import numpy as np from bdpy.dataform.sparse import load_array, save_array @@ -43,6 +44,27 @@ def test_load_save_sparse_array(self): np.testing.assert_array_equal(original_data, from_file) + def test_sparse_save_preserves_other_variables_under_numpy2(self): + if int(np.__version__.split('.')[0]) < 2: + self.skipTest('NumPy-2-only h5py writer') + + with tempfile.TemporaryDirectory() as tmpdir: + fname = os.path.join(tmpdir, 'test_sparse_preserve.mat') + + with h5py.File(fname, 'w') as f: + f.create_dataset('other', data=np.array([1, 2, 3])) + + original_data = np.random.rand(3, 2) + original_data[original_data < 0.8] = 0 + + save_array(fname, original_data, key='data', sparse=True) + from_file = load_array(fname, key='data') + + np.testing.assert_array_equal(original_data, from_file) + + with h5py.File(fname, 'r') as f: + np.testing.assert_array_equal(f['other'][()], np.array([1, 2, 3])) + def test_load_array_jl(self): data = np.array([[1, 0, 0, 0], [2, 2, 0, 0], From cb634408262ec117b31c62262b7695f1ff77a203 Mon Sep 17 00:00:00 2001 From: Kenya Otsuka Date: Wed, 17 Jun 2026 04:10:02 +0000 Subject: [PATCH 5/9] chore: cap NumPy <2 on legacy Python (<3.10) Dense save paths still use hdf5storage.savemat, which can fail under NumPy 2.x (it references the removed np.unicode_) on the legacy Python 3.8/3.9 stacks. Add environment markers so NumPy 2.x is only resolved on Python >= 3.10, while 3.8/3.9 stay on a NumPy 1.x-compatible stack. --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 63f100f..2ed609a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,7 +15,8 @@ license = { file = "LICENSE" } keywords = ["neuroscience", "neuroimaging", "brain decoding", "fmri", "machine learning"] dependencies = [ - "numpy>=1.20", + "numpy>=1.20,<2; python_version < '3.10'", + "numpy>=1.20; python_version >= '3.10'", "scipy", "scikit-learn", "h5py>=3.0", From f123d4712615604e83eff181a726a997b6503e0b Mon Sep 17 00:00:00 2001 From: Kenya Otsuka Date: Wed, 17 Jun 2026 04:10:02 +0000 Subject: [PATCH 6/9] fix: avoid full-volume reads in load_array and _check_xyz load_array inspected the HDF5 object by reading the whole dataset (f[key][()]) just to branch on dense vs sparse, which also misclassified scalar datasets. Branch on the HDF5 object type instead (group = SparseArray, dataset = dense). _check_xyz called img.get_fdata().shape, forcing a full per-file volume read only to check XYZ consistency. Use img.shape[:3] metadata instead; load_epi still reads the actual data afterwards. --- bdpy/dataform/sparse.py | 20 +++++++++++--------- bdpy/mri/load_epi.py | 4 +++- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/bdpy/dataform/sparse.py b/bdpy/dataform/sparse.py index 969b3f9..e22e0b3 100644 --- a/bdpy/dataform/sparse.py +++ b/bdpy/dataform/sparse.py @@ -24,16 +24,18 @@ def load_array(fname, key='data'): """Load an array (dense or sparse).""" with h5py.File(fname, 'r') as f: - methods = [attr for attr in dir(f[key]) if callable(getattr(f[key], str(attr)))] - if 'keys' in methods and '__bdpy_sparse_arrray' in f[key].keys(): - # SparseArray - s_ary = SparseArray(fname, key=key) - return s_ary.dense - elif type(f[key][()]) == np.ndarray: + obj = f[key] + # Inspect the HDF5 object type rather than reading the whole dataset: + # a SparseArray is stored as a group, a dense array as a dataset. + if isinstance(obj, h5py.Group): + if '__bdpy_sparse_arrray' in obj: + s_ary = SparseArray(fname, key=key) + return s_ary.dense + raise RuntimeError('Unsupported group: %s' % key) + if isinstance(obj, h5py.Dataset): # Dense array (read with h5py; hdf5storage breaks under NumPy 2.0) - return _mat_v73.read_dataset(f[key]) - else: - raise RuntimeError('Unsupported data type: %s' % type(f[key][()])) + return _mat_v73.read_dataset(obj) + raise RuntimeError('Unsupported data type: %s' % type(obj)) def save_array(fname, array, key='data', dtype=np.float64, sparse=False): diff --git a/bdpy/mri/load_epi.py b/bdpy/mri/load_epi.py index 67ee1f0..7e1d87b 100644 --- a/bdpy/mri/load_epi.py +++ b/bdpy/mri/load_epi.py @@ -47,7 +47,9 @@ def load_epi(datafiles): def _check_xyz(xyz, img): """Check voxel xyz consistency.""" - xyz_current = _get_xyz(img.coordmap.affine, img.get_fdata().shape) + # Use the image shape metadata; avoid get_fdata() here, which would read the + # whole volume just to check consistency (load_epi reads the data later). + xyz_current = _get_xyz(img.coordmap.affine, img.shape[:3]) if xyz.size == 0: xyz = xyz_current From 887258b0916923dad419196bc3167a05b25f3d5a Mon Sep 17 00:00:00 2001 From: Kenya Otsuka Date: Wed, 17 Jun 2026 04:10:02 +0000 Subject: [PATCH 7/9] test: make v7.3 feature loader test deterministic Pass format='7.3' explicitly to hdf5storage.savemat so the h5py/v7.3 loader path is always exercised instead of relying on the default, and stack the expected arrays in sorted-filename order to match Features.__get_labels. --- tests/dataform/test_features.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/dataform/test_features.py b/tests/dataform/test_features.py index 9d760eb..d5d3a10 100644 --- a/tests/dataform/test_features.py +++ b/tests/dataform/test_features.py @@ -28,11 +28,14 @@ def _prepare_mock_data( for layer_name, shape in zip(mock_layer_names, mock_shapes): os.makedirs(os.path.join(tmpdir, layer_name)) arrays = [] - for image_name in mock_image_names: + # Stack in sorted-filename order to match Features.__get_labels, which + # sorts the feature files when collecting labels. + for image_name in sorted(mock_image_names): data = np.random.rand(*shape) hdf5storage.savemat( os.path.join(tmpdir, layer_name, image_name + '.mat'), {'feat': data}, + format='7.3', store_python_metadata=True) arrays.append(data) stacked[layer_name] = np.vstack(arrays) From 5fee4256d713e7b1ebd9e38579567ba9e2d5e134 Mon Sep 17 00:00:00 2001 From: Kenya Otsuka Date: Wed, 17 Jun 2026 05:09:48 +0000 Subject: [PATCH 8/9] fix: do not treat MATLAB_empty as Python.Empty in read_dataset MATLAB-written v7.3 files mark empty arrays with MATLAB_empty but, unlike hdf5storage, do not store Python.Shape. Treating MATLAB_empty like Python.Empty made read_dataset fall back to Python.Shape=() and return np.empty(()), collapsing an empty non-scalar array such as (0, 3) to a 0-d array. Special-case only Python.Empty (falling back to the stored dataset shape if Python.Shape is missing); a bare MATLAB_empty now goes through the normal read/transpose path, preserving the empty array's shape. Add a regression test covering a MATLAB_empty dataset without Python.Shape. --- bdpy/dataform/_mat_v73.py | 12 +++++-- tests/dataform/test_mat_v73.py | 59 ++++++++++++++++++++++++++++++++++ 2 files changed, 69 insertions(+), 2 deletions(-) create mode 100644 tests/dataform/test_mat_v73.py diff --git a/bdpy/dataform/_mat_v73.py b/bdpy/dataform/_mat_v73.py index edf425b..36f8d18 100644 --- a/bdpy/dataform/_mat_v73.py +++ b/bdpy/dataform/_mat_v73.py @@ -23,6 +23,9 @@ def read_dataset(dset: h5py.Dataset) -> np.ndarray: datasets are written transposed relative to NumPy's C order. hdf5storage additionally records the original Python shape and empty-array flags as ``Python.*`` attributes, which we honor to reproduce ``hdf5storage.loadmat``. + Only ``Python.Empty`` is special-cased; a bare ``MATLAB_empty`` (set by + MATLAB without ``Python.Shape``) is read through the normal path so that + empty non-scalar arrays such as ``(0, 3)`` keep their shape. Parameters ---------- @@ -35,8 +38,13 @@ def read_dataset(dset: h5py.Dataset) -> np.ndarray: The array with its original shape restored. """ attrs = dset.attrs - if "Python.Empty" in attrs or "MATLAB_empty" in attrs: - shape = tuple(int(x) for x in attrs.get("Python.Shape", ())) + if "Python.Empty" in attrs: + # Only Python.Empty (written by hdf5storage) implies a Python.Shape we + # can trust; fall back to the stored dataset shape if it is missing. A + # bare MATLAB_empty (written by MATLAB without Python.Shape) must NOT be + # treated this way -- np.empty(()) would collapse e.g. (0, 3) to 0-d -- + # so it falls through to the normal read/transpose path below. + shape = tuple(int(x) for x in attrs.get("Python.Shape", dset.shape)) return np.empty(shape, dtype=dset.dtype) arr = dset[()] if "MATLAB_class" in attrs and isinstance(arr, np.ndarray) and arr.ndim >= 2: diff --git a/tests/dataform/test_mat_v73.py b/tests/dataform/test_mat_v73.py new file mode 100644 index 0000000..a58f37d --- /dev/null +++ b/tests/dataform/test_mat_v73.py @@ -0,0 +1,59 @@ +"""Tests for bdpy.dataform._mat_v73.""" + +import os +import tempfile +import unittest + +import h5py +import numpy as np + +from bdpy.dataform import _mat_v73 + + +class TestReadDataset(unittest.TestCase): + + def test_matlab_empty_without_python_shape_preserves_shape(self): + # MATLAB-written v7.3 files mark empty arrays with ``MATLAB_empty`` but + # do not store ``Python.Shape``. Such a dataset must not be collapsed to + # a 0-d array (the bug fixed here treated MATLAB_empty like Python.Empty + # and returned ``np.empty(())``). + with tempfile.TemporaryDirectory() as tmpdir: + fname = os.path.join(tmpdir, 'matlab_empty.mat') + + # On-disk shape (0, 3); with MATLAB_class the read path transposes + # it to (3, 0) like any other MATLAB matrix. + with h5py.File(fname, 'w') as f: + dset = f.create_dataset('a', shape=(0, 3), dtype='float64') + dset.attrs['MATLAB_empty'] = np.uint8(1) + dset.attrs['MATLAB_class'] = np.bytes_(b'double') + self.assertNotIn('Python.Shape', dset.attrs) + + with h5py.File(fname, 'r') as f: + out = _mat_v73.read_dataset(f['a']) + + self.assertNotEqual(out.shape, ()) # would fail with the old code + self.assertEqual(out.ndim, 2) + self.assertEqual(out.size, 0) + self.assertEqual(out.shape, (3, 0)) + + def test_matlab_empty_without_matlab_class_keeps_on_disk_shape(self): + # Without MATLAB_class there is no transpose, so the empty non-scalar + # shape is returned as stored. + with tempfile.TemporaryDirectory() as tmpdir: + fname = os.path.join(tmpdir, 'matlab_empty_no_class.mat') + + with h5py.File(fname, 'w') as f: + dset = f.create_dataset('a', shape=(3, 0), dtype='float64') + dset.attrs['MATLAB_empty'] = np.uint8(1) + self.assertNotIn('Python.Shape', dset.attrs) + + with h5py.File(fname, 'r') as f: + out = _mat_v73.read_dataset(f['a']) + + self.assertNotEqual(out.shape, ()) + self.assertEqual(out.shape, (3, 0)) + self.assertEqual(out.size, 0) + + +if __name__ == '__main__': + unittest.main() From 3053edb0bfe2f5e1413bc3760a835522aff6f663 Mon Sep 17 00:00:00 2001 From: Kenya Otsuka Date: Wed, 24 Jun 2026 04:46:15 +0000 Subject: [PATCH 9/9] chore: require NumPy-2-capable hdf5storage on Python >= 3.10 The dense save paths (save_array(sparse=False), save_multiarrays, save_feature) still call hdf5storage.savemat. hdf5storage < 0.2.0 is not NumPy-2 compatible: it references np.unicode_ (removed in NumPy 2.0) inside MarshallerCollection, which is constructed on every savemat/loadmat call via Options(), so the failure is not limited to the read path. Make the contract explicit, matching the existing NumPy split: keep an unconstrained hdf5storage for legacy Python (<3.10, which stays on NumPy 1.x) and require hdf5storage>=0.2.0 for Python >= 3.10, where NumPy 2.x is allowed. --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 2ed609a..d6a27a3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,7 +20,8 @@ dependencies = [ "scipy", "scikit-learn", "h5py>=3.0", - "hdf5storage", + "hdf5storage; python_version < '3.10'", + "hdf5storage>=0.2.0; python_version >= '3.10'", "pyyaml", "pandas", "tqdm",