diff --git a/CHANGELOG.md b/CHANGELOG.md index fb8ce0c..e137e4e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,19 @@ # Change Log ### All notable changes to `COMMIT` will be documented in this file. +## `v2.4.2`
_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`
_2025-09-18_ ### 🛠️Changed diff --git a/commit/core.pyx b/commit/core.pyx index 20f12d6..e8595c5 100644 --- a/commit/core.pyx +++ b/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]: diff --git a/commit/trk2dictionary/trk2dictionary.pyx b/commit/trk2dictionary/trk2dictionary.pyx index 4b60d23..072968d 100644 --- a/commit/trk2dictionary/trk2dictionary.pyx +++ b/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('=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" diff --git a/tests/demo_data/ref_results/ref_results_BallandStick.pickle b/tests/demo_data/ref_results/ref_results_BallandStick.pickle index c0b7bde..9ca004a 100644 Binary files a/tests/demo_data/ref_results/ref_results_BallandStick.pickle and b/tests/demo_data/ref_results/ref_results_BallandStick.pickle differ diff --git a/tests/demo_data/ref_results/ref_results_StickZeppelinBall.pickle b/tests/demo_data/ref_results/ref_results_StickZeppelinBall.pickle index 7435916..7cceab8 100644 Binary files a/tests/demo_data/ref_results/ref_results_StickZeppelinBall.pickle and b/tests/demo_data/ref_results/ref_results_StickZeppelinBall.pickle differ diff --git a/tests/test_demo.py b/tests/test_demo.py index ef2ab28..e6118d6 100644 --- a/tests/test_demo.py +++ b/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: