[RELEASE] v2.4.2

CHANGELOG.md

@@ -1,6 +1,19 @@
 # Change Log
 ### All notable changes to `COMMIT` will be documented in this file.
 
+## `v2.4.2`<br>_2025-10-06_
+
+### 🐛Fixed
+- Error when checking the number of streamlines in the weighted LASSO function
+- Error in verifying voxel correspondence between the lesion mask passed through the dictionary and the ISO weights map
+- Error when checking the number of streamlines in the grouping structure for group LASSO regularization (Fixes #154)
+- Added check on empty groups when setting group LASSO regularization (Fixes #153)
+- Tests fail using numpy >= 2.0.0
+- Typo in test script and reference results
+
+---
+---
+
 ## `v2.4.1`<br>_2025-09-18_
 
 ### 🛠️Changed

commit/core.pyx

@@ -796,10 +796,10 @@ cdef class Evaluation :
         if self.A is None :
             logger.error( 'Operator not built; call "build_operator()" first' )
 
-        if self.DICTIONARY['IC']['nF'] <= 0 :
+        if self.DICTIONARY['IC']['nSTR'] <= 0 :
             logger.error( 'No streamline found in the dictionary; check your data' )
 
-        if int( self.DICTIONARY['nV'] * self.KERNELS['iso'].shape[0] ) == 0 :
+        if int( self.DICTIONARY['ISO']['n'] * self.KERNELS['iso'].shape[0] ) == 0 :
             logger.error( 'Unable to set regularisation because no isotropic compartment found in the dictionary.' )
 
         # load image and check it
@@ -827,7 +827,8 @@ cdef class Evaluation :
 
         # compute array of weights from image
         array_weights = weights_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'] ].flatten().astype(np.float32)
-        if array_weights.size != self.DICTIONARY['nV']:
+        array_weights = array_weights[self.DICTIONARY['ISO']['vox']]
+        if array_weights.size != self.DICTIONARY['ISO']['n']:
             logger.error( 'Number of voxels in the weights image does not match the number of voxels in the dictionary' )
         if np.any(array_weights < 1):
             logger.error('All weights must be greater or equal to 1')
@@ -910,15 +911,16 @@ cdef class Evaluation :
             Available options:
                 'group_idx' - np.array(np.int32) :
                     group indices for the IC compartment (not implemented for EC and ISO).
-                    This field is necessary only if regularisers[0] is 'group_lasso' or 'sparse_group_lasso'.
+                    NB: Each element of the array corresponds to a group and contains the indices of the streamlines
+                    belonging to that group. All streamlines must belong to at least one group.
+                    NB: This field is necessary only if regularisers[0] is 'group_lasso' or 'sparse_group_lasso'.
                     Example:
-                        structureIC = np.array([[0,2,5],[1,3,4],[0,1,2,3,4,5],[6]], dtype=np.object_)
+                        group_idx = np.array([[0,2,5],[1,3,4],[0,1,2,3,4,5],[6]], dtype=np.object_)
                         that is equivalent to
                                     [0,1,2,3,4,5]        [6]
                                     /       \
                                 [0,2,5]       [1,3,4]
-                        which has two non-overlapping groups, one of which is the union
-                        of two other non-overlapping groups.
+                        which has two non-overlapping groups, one of which is the union of two other non-overlapping groups.
                 Group weights computation (if regularisers[0] is 'group_lasso' or 'sparse_group_lasso'):
                     w[k] = ( 1 * sqrt(|g[k]|) * group_weights_extra[k] ) / ( ||x_nnls[g[k]]||_2 )
                     where:
@@ -942,7 +944,7 @@ cdef class Evaluation :
                     This field can be specified only if regularisers[0] is 'group_lasso' or 'sparse_group_lasso'.
                     NB: this array must have the same size as the number of groups in the IC compartment and contain only non-negative values.
                 'coeff_weights' - np.array(np.float64) :
-                    weights associated to each individual element of the compartment (for the moment implemented only for IC compartment).
+                    weights associated to each individual element of the compartment (for the moment implemented only for IC and ISO compartments).
                     This field can be specified only if the chosen regulariser is 'lasso' or 'sparse_group_lasso'.
                     NB: this array must have the same size as the number of elements in the compartment and contain only non-negative values.
 
@@ -1114,10 +1116,23 @@ cdef class Evaluation :
         # check if group_idx contains all the indices of the input streamlines
         if regularisation['regIC'] == 'group_lasso' or regularisation['regIC'] == 'sparse_group_lasso':
             all_idx_in = np.sort(np.unique(np.concatenate(dictIC_params['group_idx'])))
-            all_idx = np.arange(self.DICTIONARY['IC']['nSTR'], dtype=np.int32)
+            all_idx = np.arange(self.DICTIONARY['TRK']['kept'].size, dtype=np.int32)
             if np.any(all_idx_in != all_idx):
                 logger.error('Group indices must contain all the indices of the input streamlines')
 
+        # check if in the provided group indices there is at least one empty group and remove it
+        if regularisation['regIC'] == 'group_lasso' or regularisation['regIC'] == 'sparse_group_lasso':
+            count_empty_groups = 0
+            list_non_empty_groups = []
+            for g in dictIC_params['group_idx']:
+                if g.size == 0:
+                    count_empty_groups += 1
+                else:
+                    list_non_empty_groups.append(g)
+            if count_empty_groups > 0:
+                logger.warning('At least one of the provided groups is empty; removing it from the group structure')
+            dictIC_params['group_idx'] = np.array(list_non_empty_groups, dtype=np.object_)
+
         # check if group_weights_extra is consistent with the number of groups
         if (regularisation['regIC'] == 'group_lasso' or regularisation['regIC'] == 'sparse_group_lasso') and 'group_weights_extra' in dictIC_params:
             if type(dictIC_params['group_weights_extra']) not in [list, np.ndarray]:

commit/trk2dictionary/trk2dictionary.pyx

@@ -410,7 +410,8 @@ cpdef run( filename_tractogram=None, path_out=None, filename_peaks=None, filenam
     cdef float* ptrToVOXMM
     if extension == ".tck":
         M = _get_affine( niiREF ).copy()
-        M[:3, :3] = M[:3, :3].dot( np.diag([1./Px,1./Py,1./Pz]) )
+        # float64 conversion added to comply with the new cast policy of numpy v2
+        M[:3, :3] = M[:3, :3].dot( np.diag([np.float64(1)/Px,np.float64(1)/Py,np.float64(1)/Pz]) )
         toVOXMM = np.ravel(np.linalg.inv(M)).astype('<f4')
         ptrToVOXMM = &toVOXMM[0]
 

pyproject.toml

@@ -1,18 +1,18 @@
 [build-system]
 requires = [
     "cython>=3.0.10",
-    "numpy>=1.24.4,<2.0.0",
+    "numpy>=1.24.4",
     "setuptools>=69.2.0"
 ]
 build-backend = "setuptools.build_meta"
 
 [project]
 name = "dmri-commit"
-version = "2.4.1"
+version = "2.4.2"
 dependencies = [
     "dmri-amico>=2.0.1",
     "dmri-dicelib>=1.1.0",
-    "numpy>=1.24.4,<2.0.0",
+    "numpy>=1.24.4",
     "scipy>=1.10.1"
 ]
 requires-python = ">=3.8"

tests/test_demo.py

@@ -19,23 +19,23 @@ def run_commit_StickZeppelinBall(local_path):
     trk2dictionary.run(
         filename_tractogram=os.path.join(local_path, 'demo01_fibers.tck'),
         filename_peaks=os.path.join(local_path, 'peaks.nii.gz'),
-        filename_mask=os.path.join(local_path, 'WM.nii.gz'),  # 'WM.nii.gz',
+        filename_mask=os.path.join(local_path, 'WM.nii.gz'), 
         fiber_shift=0.5,
         peaks_use_affine=True
     )
 
-    amico.util.fsl2scheme(os.path.join(local_path, 'bvals.txt'), os.path.join(
-        local_path, 'bvecs.txt'), os.path.join(local_path, 'DWI.scheme'))
+    amico.util.fsl2scheme(os.path.join(local_path, 'bvals.txt'), 
+                          os.path.join(local_path, 'bvecs.txt'), 
+                          os.path.join(local_path, 'DWI.scheme'))
 
     mit = commit.Evaluation(local_path, '.')
     mit.load_data(os.path.join(local_path, 'DWI.nii.gz'),
                   os.path.join(local_path, 'DWI.scheme'))
 
     mit.set_model('StickZeppelinBall')
-    d_par = 1.7E-3             # Parallel diffusivity [mm^2/s]
-    d_perps = [0.51E-3]        # Perpendicular diffusivity(s) [mm^2/s]
-    # d_perps = []        # Perpendicular diffusivity(s) [mm^2/s]
-    d_isos = [1.7E-3, 3.0E-3]  # Isotropic diffusivity(s) [mm^2/s]
+    d_par   = 1.7E-3            # Parallel diffusivity [mm^2/s]
+    d_perps = [0.51E-3]         # Perpendicular diffusivity(s) [mm^2/s]
+    d_isos  = [1.7E-3, 3.0E-3]  # Isotropic diffusivity(s) [mm^2/s]
     mit.model.set(d_par, d_perps, d_isos)
     mit.generate_kernels(regenerate=True)
     mit.load_kernels()
@@ -53,22 +53,23 @@ def run_commit_BallandStick(local_path):
     trk2dictionary.run(
         filename_tractogram=os.path.join(local_path, 'demo01_fibers.tck'),
         filename_peaks=os.path.join(local_path, 'peaks.nii.gz'),
-        filename_mask=os.path.join(local_path, 'WM.nii.gz'),  # 'WM.nii.gz',
+        filename_mask=os.path.join(local_path, 'WM.nii.gz'), 
         fiber_shift=0.5,
         peaks_use_affine=True
     )
 
-    amico.util.fsl2scheme(os.path.join(local_path, 'bvals.txt'), os.path.join(
-        local_path, 'bvecs.txt'), os.path.join(local_path, 'DWI.scheme'))
+    amico.util.fsl2scheme(os.path.join(local_path, 'bvals.txt'), 
+                          os.path.join(local_path, 'bvecs.txt'), 
+                          os.path.join(local_path, 'DWI.scheme'))
 
     mit = commit.Evaluation(local_path, '.')
     mit.load_data(os.path.join(local_path, 'DWI.nii.gz'),
                   os.path.join(local_path, 'DWI.scheme'))
 
     mit.set_model('StickZeppelinBall')
-    d_par = 1.7E-3             # Parallel diffusivity [mm^2/s]
-    d_perps = []        # Perpendicular diffusivity(s) [mm^2/s]
-    d_isos = [1.7E-3, 3.0E-3]  # Isotropic diffusivity(s) [mm^2/s]
+    d_par   = 1.7E-3            # Parallel diffusivity [mm^2/s]
+    d_perps = []                # Perpendicular diffusivity(s) [mm^2/s]
+    d_isos  = [1.7E-3, 3.0E-3]  # Isotropic diffusivity(s) [mm^2/s]
     mit.model.set(d_par, d_perps, d_isos)
     mit.generate_kernels(regenerate=True)
     mit.load_kernels()
@@ -86,22 +87,23 @@ def run_commit_VolumeFractions(local_path):
     trk2dictionary.run(
         filename_tractogram=os.path.join(local_path, 'demo01_fibers.tck'),
         filename_peaks=os.path.join(local_path, 'peaks.nii.gz'),
-        filename_mask=os.path.join(local_path, 'WM.nii.gz'),  # 'WM.nii.gz',
+        filename_mask=os.path.join(local_path, 'WM.nii.gz'),
         fiber_shift=0.5,
         peaks_use_affine=True
     )
 
-    amico.util.fsl2scheme(os.path.join(local_path, 'bvals.txt'), os.path.join(
-        local_path, 'bvecs.txt'), os.path.join(local_path, 'DWI.scheme'))
+    amico.util.fsl2scheme(os.path.join(local_path, 'bvals.txt'), 
+                          os.path.join(local_path, 'bvecs.txt'), 
+                          os.path.join(local_path, 'DWI.scheme'))
 
     mit = commit.Evaluation(local_path, '.')
     mit.load_data(os.path.join(local_path, 'DWI.nii.gz'),
                   os.path.join(local_path, 'DWI.scheme'))
 
     mit.set_model('StickZeppelinBall')
-    d_par = 1.7E-3             # Parallel diffusivity [mm^2/s]
-    d_perps = []        # Perpendicular diffusivity(s) [mm^2/s]
-    d_isos = [1.7E-3, 3.0E-3]  # Isotropic diffusivity(s) [mm^2/s]
+    d_par   = 1.7E-3            # Parallel diffusivity [mm^2/s]
+    d_perps = []                # Perpendicular diffusivity(s) [mm^2/s]
+    d_isos  = [1.7E-3, 3.0E-3]  # Isotropic diffusivity(s) [mm^2/s]
     mit.model.set(d_par, d_perps, d_isos)
     mit.generate_kernels(regenerate=True)
     mit.load_kernels()
@@ -131,7 +133,7 @@ def check_results(pickle_result, ref_pickle):
         assert abs(result_optimization["fit_details"]["abs_cost"] - ref_optimization["fit_details"]["abs_cost"]) < 1e-4
         assert abs(result_optimization["fit_details"]["rel_cost"] - ref_optimization["fit_details"]["rel_cost"]) < 1e-4
         assert abs(result_optimization["fit_details"]["abs_x"] - ref_optimization["fit_details"]["abs_x"]) < 1e-4
-        assert abs(result_optimization["fit_details"]["rel_x"] - ref_optimization["fit_details"]["rel _x"]) < 1e-4 #NOTE: There is a space in the key name
+        assert abs(result_optimization["fit_details"]["rel_x"] - ref_optimization["fit_details"]["rel_x"]) < 1e-4
         assert result_optimization["fit_details"]["iterations"] == ref_optimization["fit_details"]["iterations"]
 
     except AssertionError: