Support NumPy 2.x and NiBabel 5.x

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)

bdpy/dataform/datastore.py

@@ -1,7 +1,7 @@
 """DataStore class
 
 This file is a part of BdPy.
 """
 
 
 from __future__ import print_function
@@ -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']
 
 
@@ -52,8 +53,8 @@
     - Add default file name pattern (`pattern`).
     """
 
     def __init__(self,
                  dpath=None, file_type=None,
                  pattern=None, extractor=None):
         self.__key_sep = '/'
 
@@ -126,7 +127,7 @@
 
     def __get_key_num(self, pat):
         """Return no. of keys."""
         n_keys = len(re.findall('\(.*?\)', pat))
         return n_keys
 
     def __parse_datafiles(self, dpath):
@@ -250,7 +251,8 @@
         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

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

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

bdpy/mri/fmriprep.py

@@ -95,13 +95,13 @@
         # File name patterns
         # FIXME
         if self.__fmriprep_version == '1.2':
             file_pattern = {'volume_native'   : '.*_space-T1w_desc-preproc_bold\.nii\.gz$',
                             'volume_standard' : '.*_space-MNI152NLin2009cAsym_desc-preproc_bold\.nii\.gz$',
                             'surf_native_left'  : '.*_space-fsnative_hemi-L\.func\.gii$',
                             'surf_native_right' : '.*_space-fsnative_hemi-R\.func\.gii$',
                             'surf_standard_left'  : '.*_space-fsaverage_hemi-L\.func\.gii$',
                             'surf_standard_right' : '.*_space-fsaverage_hemi-R\.func\.gii$',
                             'surf_standard_41k_left'  : '.*_space-fsaverage6_hemi-L\.func\.gii$',
                             'surf_standard_41k_right' : '.*_space-fsaverage6_hemi-R\.func\.gii$',
                             'surf_standard_10k_left'  : '.*_space-fsaverage5_hemi-L\.func\.gii$',
                             'surf_standard_10k_right' : '.*_space-fsaverage5_hemi-R\.func\.gii$',
@@ -408,7 +408,7 @@
         """
         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 @@
         # 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 @@
     """
     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

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

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

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