diff --git a/bdpy/dataform/_mat_v73.py b/bdpy/dataform/_mat_v73.py new file mode 100644 index 00000000..36f8d18c --- /dev/null +++ b/bdpy/dataform/_mat_v73.py @@ -0,0 +1,127 @@ +"""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``. + 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 + ---------- + dset : h5py.Dataset + Dataset to read. + + Returns + ------- + numpy.ndarray + The array with its original shape restored. + """ + attrs = dset.attrs + 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: + 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 00e4d914..4098b157 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 46fee230..193537f0 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 2e4593fc..e22e0b3d 100644 --- a/bdpy/dataform/sparse.py +++ b/bdpy/dataform/sparse.py @@ -12,20 +12,30 @@ import hdf5storage import numpy as np +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).""" 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: - # Dense array - return hdf5storage.loadmat(fname)[key] - else: - raise RuntimeError('Unsupported data type: %s' % type(f[key][()])) + 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(obj) + raise RuntimeError('Unsupported data type: %s' % type(obj)) def save_array(fname, array, key='data', dtype=np.float64, sparse=False): @@ -75,12 +85,40 @@ 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) + 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): @@ -95,29 +133,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/bdpy/mri/fmriprep.py b/bdpy/mri/fmriprep.py index ff4e911a..6567a69d 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 d57e218a..2916cac1 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 ba1fd78b..7e1d87b8 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,9 @@ def load_epi(datafiles): def _check_xyz(xyz, img): """Check voxel xyz consistency.""" - xyz_current = _get_xyz(img.coordmap.affine, img.get_data().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 diff --git a/bdpy/mri/load_mri.py b/bdpy/mri/load_mri.py index c5601901..772d6b04 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/pyproject.toml b/pyproject.toml index 8f153111..d6a27a36 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,11 +15,13 @@ 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", - "hdf5storage", + "hdf5storage; python_version < '3.10'", + "hdf5storage>=0.2.0; python_version >= '3.10'", "pyyaml", "pandas", "tqdm", @@ -40,9 +42,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", diff --git a/tests/dataform/test_features.py b/tests/dataform/test_features.py index 927341be..d5d3a104 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,29 @@ 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)) - for image_name in mock_image_names: + arrays = [] + # 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='5') + format='7.3', + store_python_metadata=True) + arrays.append(data) + stacked[layer_name] = np.vstack(arrays) + return stacked class TestDataformFeatures(unittest.TestCase): @@ -45,29 +56,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() diff --git a/tests/dataform/test_mat_v73.py b/tests/dataform/test_mat_v73.py new file mode 100644 index 00000000..a58f37d6 --- /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() diff --git a/tests/dataform/test_sparse.py b/tests/dataform/test_sparse.py index 20810c9c..3047f6ba 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], diff --git a/tests/test_mri.py b/tests/test_mri.py index 3922f340..bee53ae7 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],