From c83437d3f4a209de089e40f2dcca9e8f482e854a Mon Sep 17 00:00:00 2001 From: Ilaria Gabusi Date: Tue, 12 May 2026 14:50:31 +0200 Subject: [PATCH 1/7] Add warning if all IC columns removed in debias --- commit/core.pyx | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/commit/core.pyx b/commit/core.pyx index e8595c5..0e5bb0d 100644 --- a/commit/core.pyx +++ b/commit/core.pyx @@ -1473,6 +1473,7 @@ cdef class Evaluation : logger.subinfo(f'Stopped after {opt_details["iterations"]} iterations', indent_lvl=1, indent_char='*', with_progress=True) logger.subinfo(f'Stopping condition: \"{opt_details["stopping_criterion"]}\"', indent_lvl=1, indent_char='*') + # DEBIAS if (self.regularisation_params['regIC']!=None or self.regularisation_params['regEC']!= None or self.regularisation_params['regISO']!= None) and debias: from commit.operator import operator temp_verb = self.verbose @@ -1484,28 +1485,31 @@ cdef class Evaluation : mask = np.ones(offset, dtype=np.uint32) mask[xic<0.000000000000001] = 0 - self.DICTIONARY["IC"]["eval"] = mask + if np.sum(mask)==0: + logger.warning('All coefficients of the IC compartment are below the debias threshold. The debias step will not be performed. Note: consider softening the regularisation by decreasing the lambda value(s).') + else: + self.DICTIONARY["IC"]["eval"] = mask - self.A = operator.LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS, nolut=True if hasattr(self.model, 'nolut') else False ) + self.A = operator.LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS, nolut=True if hasattr(self.model, 'nolut') else False ) - self.set_regularisation() - self.set_verbose(temp_verb) + self.set_regularisation() + self.set_verbose(temp_verb) - logger.subinfo('Recomputing coefficients', indent_lvl=1, indent_char='*', with_progress=True) + logger.subinfo('Recomputing coefficients', indent_lvl=1, indent_char='*', with_progress=True) - x_debias = self.x.copy() - x_debias[:offset] *= mask - x_debias[offset:] = 0 + x_debias = self.x.copy() + x_debias[:offset] *= mask + x_debias[offset:] = 0 - y_mask = np.asarray(self.A.dot(x_debias)) - # binarize y_debias - y_mask[y_mask<0] = 0 - y_mask[y_mask>0] = 1 + y_mask = np.asarray(self.A.dot(x_debias)) + # binarize + y_mask[y_mask<0] = 0 + y_mask[y_mask>0] = 1 - self.debias_mask = y_mask + self.debias_mask = y_mask - with ProgressBar(disable=self.verbose!=3, hide_on_exit=True, subinfo=True) as pbar: - self.x, opt_details = commit.solvers.solve(self.get_y(), self.A, self.A.T, tol_fun=tol_fun, tol_x=tol_x, max_iter=max_iter, verbose=self.verbose, x0=x0, regularisation=self.regularisation_params, confidence_array=confidence_array) + with ProgressBar(disable=self.verbose!=3, hide_on_exit=True, subinfo=True) as pbar: + self.x, opt_details = commit.solvers.solve(self.get_y(), self.A, self.A.T, tol_fun=tol_fun, tol_x=tol_x, max_iter=max_iter, verbose=self.verbose, x0=x0, regularisation=self.regularisation_params, confidence_array=confidence_array) elif (self.regularisation_params['regIC']!=None or self.regularisation_params['regEC']!= None or self.regularisation_params['regISO']!= None) and not debias: logger.warning('Fitting with regularisation but without debiasing. The coefficients will be biased, use "debias=True" to debias the coefficients') From 0f959a5c2917baf8a29f0a0a59db3e96229c6c45 Mon Sep 17 00:00:00 2001 From: Ilaria Gabusi Date: Tue, 12 May 2026 17:29:24 +0200 Subject: [PATCH 2/7] Add threshold parameter for debias in fit --- commit/core.pyx | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/commit/core.pyx b/commit/core.pyx index 0e5bb0d..6ab78e3 100644 --- a/commit/core.pyx +++ b/commit/core.pyx @@ -1350,7 +1350,7 @@ cdef class Evaluation : logger.info( f'[ {format_time(time.time() - tr)} ]' ) - def fit( self, tol_fun=1e-3, tol_x=1e-6, max_iter=100, x0=None, confidence_map_filename=None, confidence_map_rescale=False, debias=False ) : + def fit( self, tol_fun=1e-3, tol_x=1e-6, max_iter=100, x0=None, confidence_map_filename=None, confidence_map_rescale=False, debias=False, thr_debias=1e-15 ): """Fit the model to the data. Parameters @@ -1371,6 +1371,10 @@ cdef class Evaluation : If true, the values of the confidence map will be rescaled to the range [0.0,1.0]. Only the voxels considered in the mask will be affected. (default : False) + debias : boolean + If true, the debias step will be performed (default : False) + thr_debias : float + Threshold for the debiasing step (default : 1e-15) """ if self.niiDWI is None : logger.error( 'Data not loaded; call "load_data()" first' ) @@ -1477,13 +1481,13 @@ cdef class Evaluation : if (self.regularisation_params['regIC']!=None or self.regularisation_params['regEC']!= None or self.regularisation_params['regISO']!= None) and debias: from commit.operator import operator temp_verb = self.verbose - logger.info( 'Running debias' ) + logger.info( f'Running debias (with threshold={thr_debias:.2e})' ) self.set_verbose(0) offset = self.DICTIONARY['IC']['nSTR'] * self.KERNELS['wmr'].shape[0] xic = self.x[:offset] mask = np.ones(offset, dtype=np.uint32) - mask[xic<0.000000000000001] = 0 + mask[xic Date: Tue, 12 May 2026 18:49:00 +0200 Subject: [PATCH 3/7] Fix computation of RMSE and NRMSE after debias and remove debias_mask from Evaluation class --- commit/core.pyx | 94 +++++++++++-------------------------------------- 1 file changed, 20 insertions(+), 74 deletions(-) diff --git a/commit/core.pyx b/commit/core.pyx index 6ab78e3..212c9f6 100644 --- a/commit/core.pyx +++ b/commit/core.pyx @@ -76,7 +76,6 @@ cdef class Evaluation : cdef public CONFIG cdef public temp_data cdef public confidence_map_img - cdef public debias_mask cdef public verbose def __init__( self, study_path='.', subject='.', dictionary_path='COMMIT'): @@ -102,7 +101,6 @@ cdef class Evaluation : self.regularisation_params = None # set by "set_regularisation" method self.x = None # set by "fit" method self.confidence_map_img = None # set by "fit" method - self.debias_mask = None # set by "fit" method self.x_nnls = None # set by "fit" method (coefficients of IC compartment estimated without regularization) self.verbose = 3 @@ -760,9 +758,6 @@ cdef class Evaluation : y = self.niiDWI_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'], : ].flatten().astype(np.float64) - if self.debias_mask is not None : - y *= self.debias_mask - return y def set_wLasso_ISO(self, img_weights_filename, lambda_perc_iso): @@ -1501,17 +1496,6 @@ cdef class Evaluation : logger.subinfo('Recomputing coefficients', indent_lvl=1, indent_char='*', with_progress=True) - x_debias = self.x.copy() - x_debias[:offset] *= mask - x_debias[offset:] = 0 - - y_mask = np.asarray(self.A.dot(x_debias)) - # binarize - y_mask[y_mask<0] = 0 - y_mask[y_mask>0] = 1 - - self.debias_mask = y_mask - with ProgressBar(disable=self.verbose!=3, hide_on_exit=True, subinfo=True) as pbar: self.x, opt_details = commit.solvers.solve(self.get_y(), self.A, self.A.T, tol_fun=tol_fun, tol_x=tol_x, max_iter=max_iter, verbose=self.verbose, x0=x0, regularisation=self.regularisation_params, confidence_array=confidence_array) @@ -1629,64 +1613,26 @@ cdef class Evaluation : niiMAP_hdr['descrip'] = f'Created with COMMIT {self.get_config("version")}' niiMAP_hdr['db_name'] = '' - if self.debias_mask is not None: - nVOX = int(np.sum(self.debias_mask)/self.niiDWI_img.shape[3]) - ind_mask = np.where(self.debias_mask>0)[0] - - y_mea = np.reshape( self.get_y()[ind_mask], (nVOX,-1) ) - y_est = np.reshape( np.asarray(self.A.dot(self.x))[ind_mask], (nVOX,-1) ) - tmp = np.sqrt( np.mean((y_mea-y_est)**2,axis=1) ) - - logger.subinfo(f'RMSE: {tmp.mean():.3f} +/- {tmp.std():.3f}', indent_lvl=2, indent_char='-') - - tmp = np.sum(y_mea**2,axis=1) - idx = np.where( tmp < 1E-12 ) - tmp[ idx ] = 1 - tmp = np.sqrt( np.sum((y_mea-y_est)**2,axis=1) / tmp ) - tmp[ idx ] = 0 - logger.subinfo(f'NRMSE: {tmp.mean():.3f} +/- {tmp.std():.3f}', indent_lvl=2, indent_char='-') - - y_mea = np.reshape( self.get_y(), (self.DICTIONARY['IC']['nVOX'],-1) ) - y_est = np.reshape( self.A.dot(self.x), (self.DICTIONARY['IC']['nVOX'],-1) ).astype(np.float32) - tmp = np.sqrt( np.mean((y_mea-y_est)**2,axis=1) ) - - niiMAP_img[self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz']] = tmp - niiMAP_hdr['cal_min'] = 0 - niiMAP_hdr['cal_max'] = tmp.max() - nibabel.save( niiMAP, pjoin(RESULTS_path,'fit_RMSE.nii.gz') ) - - tmp = np.sum(y_mea**2,axis=1) - idx = np.where( tmp < 1E-12 ) - tmp[ idx ] = 1 - tmp = np.sqrt( np.sum((y_mea-y_est)**2,axis=1) / tmp ) - tmp[ idx ] = 0 - - niiMAP_img[self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz']] = tmp - niiMAP_hdr['cal_min'] = 0 - niiMAP_hdr['cal_max'] = 1 - nibabel.save( niiMAP, pjoin(RESULTS_path,'fit_NRMSE.nii.gz') ) - - else: - nVOX = self.DICTIONARY['IC']['nVOX'] - y_mea = np.reshape( self.niiDWI_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'], : ].flatten().astype(np.float32), (nVOX,-1) ) - y_est = np.reshape( self.A.dot(self.x), (nVOX,-1) ).astype(np.float32) - tmp = np.sqrt( np.mean((y_mea-y_est)**2,axis=1) ) - logger.subinfo(f'RMSE: {tmp.mean():.3f} +/- {tmp.std():.3f}', indent_lvl=2, indent_char='-') - niiMAP_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'] ] = tmp - niiMAP_hdr['cal_min'] = 0 - niiMAP_hdr['cal_max'] = tmp.max() - nibabel.save( niiMAP, pjoin(RESULTS_path,'fit_RMSE.nii.gz') ) - - tmp = np.sum(y_mea**2,axis=1) - idx = np.where( tmp < 1E-12 ) - tmp[ idx ] = 1 - tmp = np.sqrt( np.sum((y_mea-y_est)**2,axis=1) / tmp ) - tmp[ idx ] = 0 - logger.subinfo(f'NRMSE: {tmp.mean():.3f} +/- {tmp.std():.3f}', indent_lvl=2, indent_char='-') - niiMAP_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'] ] = tmp - niiMAP_hdr['cal_min'] = 0 - niiMAP_hdr['cal_max'] = 1 - nibabel.save( niiMAP, pjoin(RESULTS_path,'fit_NRMSE.nii.gz') ) + nVOX = self.DICTIONARY['IC']['nVOX'] + y_mea = np.reshape( self.niiDWI_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'], : ].flatten().astype(np.float32), (nVOX,-1) ) + y_est = np.reshape( self.A.dot(self.x), (nVOX,-1) ).astype(np.float32) + tmp = np.sqrt( np.mean((y_mea-y_est)**2,axis=1) ) + logger.subinfo(f'RMSE: {tmp.mean():.3f} +/- {tmp.std():.3f}', indent_lvl=2, indent_char='-') + niiMAP_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'] ] = tmp + niiMAP_hdr['cal_min'] = 0 + niiMAP_hdr['cal_max'] = tmp.max() + nibabel.save( niiMAP, pjoin(RESULTS_path,'fit_RMSE.nii.gz') ) + + tmp = np.sum(y_mea**2,axis=1) + idx = np.where( tmp < 1E-12 ) + tmp[ idx ] = 1 + tmp = np.sqrt( np.sum((y_mea-y_est)**2,axis=1) / tmp ) + tmp[ idx ] = 0 + logger.subinfo(f'NRMSE: {tmp.mean():.3f} +/- {tmp.std():.3f}', indent_lvl=2, indent_char='-') + niiMAP_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'] ] = tmp + niiMAP_hdr['cal_min'] = 0 + niiMAP_hdr['cal_max'] = 1 + nibabel.save( niiMAP, pjoin(RESULTS_path,'fit_NRMSE.nii.gz') ) if self.confidence_map_img is not None: confidence_array = np.reshape( self.confidence_map_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'], : ].flatten().astype(np.float32), (nVOX,-1) ) From 786f57e6256b396ad72c762529b3c8dfc810bbbe Mon Sep 17 00:00:00 2001 From: Ilaria Gabusi Date: Tue, 12 May 2026 18:53:03 +0200 Subject: [PATCH 4/7] Update changelog and version --- CHANGELOG.md | 8 ++++++++ pyproject.toml | 2 +- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e137e4e..970d0eb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,14 @@ # Change Log ### All notable changes to `COMMIT` will be documented in this file. +## `v2.4.3`
_2026-05-12_ + +### 🐛Fixed +- Occasional error in saving results (RMSE and NRMSE maps) after debias when using BS or SZB models (Fixes #159) + +--- +--- + ## `v2.4.2`
_2025-10-06_ ### 🐛Fixed diff --git a/pyproject.toml b/pyproject.toml index 4d6f5cb..e0c11a6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,7 +8,7 @@ build-backend = "setuptools.build_meta" [project] name = "dmri-commit" -version = "2.4.2" +version = "2.4.3" dependencies = [ "dmri-amico>=2.0.1", "dmri-dicelib>=1.1.0", From 4b043e104f89562ebe8d6e985f5fea2c576bfc94 Mon Sep 17 00:00:00 2001 From: Ilaria Gabusi Date: Wed, 13 May 2026 12:15:44 +0200 Subject: [PATCH 5/7] Call build operator at the beginning of fit() --- commit/core.pyx | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/commit/core.pyx b/commit/core.pyx index 212c9f6..38385bc 100644 --- a/commit/core.pyx +++ b/commit/core.pyx @@ -722,9 +722,14 @@ cdef class Evaluation : logger.info( f'[ {format_time(time.time() - tic)} ]' ) - def build_operator( self ) : + def build_operator( self, mask_ic=None ) : """Build the operator for computing the matrix-vector multiplications by A and A' using the informations from self.DICTIONARY, self.KERNELS and self.THREADS. + + Parameters + ---------- + mask_ic : np.array + Binary mask to restrict the evaluation on a subset of columns of the IC compartment. """ if self.DICTIONARY is None : logger.error( 'Dictionary not loaded; call "load_dictionary()" first' ) @@ -741,7 +746,10 @@ cdef class Evaluation : tic = time.time() logger.subinfo('') logger.info( 'Building linear operator A' ) - self.DICTIONARY["IC"]["eval"] = np.ones( int(self.DICTIONARY['IC']['nSTR'] * self.KERNELS['wmr'].shape[0]), dtype=np.uint32) + if mask_ic is not None: + self.DICTIONARY["IC"]["eval"] = mask_ic + else: + self.DICTIONARY["IC"]["eval"] = np.ones( int(self.DICTIONARY['IC']['nSTR'] * self.KERNELS['wmr'].shape[0]), dtype=np.uint32) self.A = operator.LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS, True if hasattr(self.model, 'nolut') else False ) logger.info( f'[ {format_time(time.time() - tic)} ]' ) @@ -1379,9 +1387,11 @@ cdef class Evaluation : logger.error( 'Response functions not generated; call "generate_kernels()" and "load_kernels()" first' ) if self.THREADS is None : logger.error( 'Threads not set; call "set_threads()" first' ) - if self.A is None : - logger.error( 'Operator not built; call "build_operator()" first' ) + # Build operator + self.build_operator() + + # Set default regularisation parameters if not set by the user if self.regularisation_params is None: self.set_regularisation() @@ -1474,7 +1484,6 @@ cdef class Evaluation : # DEBIAS if (self.regularisation_params['regIC']!=None or self.regularisation_params['regEC']!= None or self.regularisation_params['regISO']!= None) and debias: - from commit.operator import operator temp_verb = self.verbose logger.info( f'Running debias (with threshold={thr_debias:.2e})' ) self.set_verbose(0) @@ -1485,17 +1494,15 @@ cdef class Evaluation : mask[xic Date: Wed, 13 May 2026 12:43:23 +0200 Subject: [PATCH 6/7] Rename variables regarding IC 'eval' to 'mask' --- commit/core.pyx | 8 ++++--- commit/operator/operator.pyx | 12 +++++----- setup_operator.py | 46 ++++++++++++++++++------------------ 3 files changed, 34 insertions(+), 32 deletions(-) diff --git a/commit/core.pyx b/commit/core.pyx index 38385bc..e17149d 100644 --- a/commit/core.pyx +++ b/commit/core.pyx @@ -747,9 +747,9 @@ cdef class Evaluation : logger.subinfo('') logger.info( 'Building linear operator A' ) if mask_ic is not None: - self.DICTIONARY["IC"]["eval"] = mask_ic + self.DICTIONARY["IC"]["mask"] = mask_ic else: - self.DICTIONARY["IC"]["eval"] = np.ones( int(self.DICTIONARY['IC']['nSTR'] * self.KERNELS['wmr'].shape[0]), dtype=np.uint32) + self.DICTIONARY["IC"]["mask"] = np.ones( int(self.DICTIONARY['IC']['nSTR'] * self.KERNELS['wmr'].shape[0]), dtype=np.uint32) self.A = operator.LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS, True if hasattr(self.model, 'nolut') else False ) logger.info( f'[ {format_time(time.time() - tic)} ]' ) @@ -1485,7 +1485,9 @@ cdef class Evaluation : # DEBIAS if (self.regularisation_params['regIC']!=None or self.regularisation_params['regEC']!= None or self.regularisation_params['regISO']!= None) and debias: temp_verb = self.verbose - logger.info( f'Running debias (with threshold={thr_debias:.2e})' ) + logger.subinfo('') + logger.info( f'Running debias' ) + logger.subinfo( f'Creating mask for IC compartment with threshold={thr_debias:.2e}', indent_lvl=1, indent_char='*' ) self.set_verbose(0) offset = self.DICTIONARY['IC']['nSTR'] * self.KERNELS['wmr'].shape[0] diff --git a/commit/operator/operator.pyx b/commit/operator/operator.pyx index 8bb5ebc..67fcfec 100755 --- a/commit/operator/operator.pyx +++ b/commit/operator/operator.pyx @@ -21,7 +21,7 @@ cdef class LinearOperator : unsigned int [::1] ICv unsigned int [::1] ECv unsigned int [::1] ISOv - unsigned int [::1] ICeval + unsigned int [::1] ICm float [::1] ICl float [:, :, ::1] LUT_IC float [:, :, ::1] LUT_EC @@ -74,7 +74,7 @@ cdef class LinearOperator : self.ICl = DICTIONARY['IC']['len'] self.ICv = DICTIONARY['IC']['vox'] self.ICo = DICTIONARY['IC']['dir'] - self.ICeval = DICTIONARY['IC']['eval'] + self.ICm = DICTIONARY['IC']['mask'] self.ECv = DICTIONARY['EC']['vox'] self.ECo = DICTIONARY['EC']['dir'] self.ISOv = DICTIONARY['ISO']['vox'] @@ -140,7 +140,7 @@ cdef class LinearOperator : COMMIT_A_nolut( &v_in[0], &v_out[0], self.ICnSTR, - &self.ICf[0], &self.ICeval[0], &self.ICv[0], &self.ICl[0], + &self.ICf[0], &self.ICm[0], &self.ICv[0], &self.ICl[0], &self.ISOv[0], &self.ICthreads[0], &self.ISOthreads[0], self.ISOnRF, self.nThreads @@ -150,7 +150,7 @@ cdef class LinearOperator : &v_in[0], &v_out[0], self.nSAMPLES, self.ndirs, self.ICnSTR, self.ECn, self.ISOn, - &self.ICf[0], &self.ICeval[0], &self.ICv[0], &self.ICo[0], &self.ICl[0], + &self.ICf[0], &self.ICm[0], &self.ICv[0], &self.ICo[0], &self.ICl[0], &self.ECv[0], &self.ECo[0], &self.ISOv[0], &self.LUT_IC[0,0,0], &self.LUT_EC[0,0,0], &self.LUT_ISO[0,0], @@ -163,7 +163,7 @@ cdef class LinearOperator : COMMIT_At_nolut( &v_in[0], &v_out[0], self.ICnSTR, self.ICn, - &self.ICf[0], &self.ICeval[0], &self.ICv[0], &self.ICl[0], + &self.ICf[0], &self.ICm[0], &self.ICv[0], &self.ICl[0], &self.ISOv[0], &self.ICthreadsT[0], &self.ISOthreadsT[0], self.ISOnRF, self.nThreads @@ -173,7 +173,7 @@ cdef class LinearOperator : &v_in[0], &v_out[0], self.nSAMPLES, self.ndirs, self.ICnSTR, self.ICn, self.ECn, self.ISOn, - &self.ICf[0], &self.ICeval[0], &self.ICv[0], &self.ICo[0], &self.ICl[0], + &self.ICf[0], &self.ICm[0], &self.ICv[0], &self.ICo[0], &self.ICl[0], &self.ECv[0], &self.ECo[0], &self.ISOv[0], &self.LUT_IC[0,0,0], &self.LUT_EC[0,0,0], &self.LUT_ISO[0,0], diff --git a/setup_operator.py b/setup_operator.py index 5ed9506..62db56e 100644 --- a/setup_operator.py +++ b/setup_operator.py @@ -35,7 +35,7 @@ def add_global_variables() -> str: uint32_t *ICthreads, *ECthreads, *ISOthreads; uint8_t *ICthreadsT; uint32_t *ECthreadsT, *ISOthreadsT; -uint32_t *ICf, *ICeval, *ICv, *ECv, *ISOv; +uint32_t *ICf, *ICm, *ICv, *ECv, *ISOv; uint16_t *ICo, *ECo; float *ICl; float *wmrSFP0, *wmrSFP1, *wmrSFP2, *wmrSFP3, *wmrSFP4, *wmrSFP5, *wmrSFP6, *wmrSFP7, *wmrSFP8, *wmrSFP9, *wmrSFP10, *wmrSFP11, *wmrSFP12, *wmrSFP13, *wmrSFP14, *wmrSFP15, *wmrSFP16, *wmrSFP17, *wmrSFP18, *wmrSFP19; @@ -61,8 +61,8 @@ def add_intracellular_compartments() -> str: {''' s1: str = '''\ xPtr0 = x + (*t_f); - eval0 = ICeval + *t_f; - x0 = *xPtr0 * (double)(*eval0);''' + m0 = ICm + *t_f; + x0 = *xPtr0 * (double)(*m0);''' s2: str = 'if (x0 != 0' s3: str = 'SFP0ptr = wmrSFP0 + offset;' s4: str = 'x0 * (*SFP0ptr++)' @@ -72,8 +72,8 @@ def add_intracellular_compartments() -> str: s1 += f'''\ xPtr{i} = xPtr{i - 1} + nF; - eval{i} = eval{i - 1} + nF; - x{i} = *xPtr{i} * (double)(*eval{i});''' + m{i} = m{i - 1} + nF; + x{i} = *xPtr{i} * (double)(*m{i});''' s2 += f' || x{i} != 0' s3 += f'''\ @@ -225,7 +225,7 @@ def add_isotropic_compartments() -> str: { int id = (long)ptr; int offset; - uint32_t *eval0, *eval1, *eval2, *eval3, *eval4, *eval5, *eval6, *eval7, *eval8, *eval9, *eval10, *eval11, *eval12, *eval13, *eval14, *eval15, *eval16, *eval17, *eval18, *eval19; + uint32_t *m0, *m1, *m2, *m3, *m4, *m5, *m6, *m7, *m8, *m9, *m10, *m11, *m12, *m13, *m14, *m15, *m16, *m17, *m18, *m19; double x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, w; double *xPtr0, *xPtr1, *xPtr2, *xPtr3, *xPtr4, *xPtr5, *xPtr6, *xPtr7, *xPtr8, *xPtr9, *xPtr10, *xPtr11, *xPtr12, *xPtr13, *xPtr14, *xPtr15, *xPtr16, *xPtr17, *xPtr18, *xPtr19; double *YPtr, *YPtrEnd; @@ -315,7 +315,7 @@ def add_isotropic_compartments() -> str: double *_vIN, double *_vOUT, int _nS, int _ndirs, int _nF, int _nE, int _nV, - uint32_t *_ICf, uint32_t *_ICeval, uint32_t *_ICv, uint16_t *_ICo, float *_ICl, + uint32_t *_ICf, uint32_t *_ICm, uint32_t *_ICv, uint16_t *_ICo, float *_ICl, uint32_t *_ECv, uint16_t *_ECo, uint32_t *_ISOv, float *_wmrSFP, float *_wmhSFP, float *_isoSFP, @@ -333,7 +333,7 @@ def add_isotropic_compartments() -> str: Y = _vOUT; ICf = _ICf; - ICeval = _ICeval; + ICm = _ICm; ICv = _ICv; ICo = _ICo; ICl = _ICl; @@ -383,10 +383,10 @@ def add_intracellular_compartments() -> str: s1: str = 'SFP0ptr = wmrSFP0 + offset;' s2: str = '''\ x0 = (*SFP0ptr++) * YTmp; - eval0 = ICeval + *t_f; + m0 = ICm + *t_f; ''' s3: str = 'x0 += (*SFP0ptr++) * YTmp;' - s4: str = 'x[*t_f] += w * x0 * (double)(*eval0);' + s4: str = 'x[*t_f] += w * x0 * (double)(*m0);' s5: str = '' for i in range(0, 20): @@ -399,13 +399,13 @@ def add_intracellular_compartments() -> str: s2 += f'''\ x{i} = (*SFP{i}ptr++) * YTmp; - eval{i} = eval{i - 1} + nF;''' + m{i} = m{i - 1} + nF;''' s3 += f'''\ x{i} += (*SFP{i}ptr++) * YTmp;''' s4 += f'''\ - x[*t_f+{s5}nF] += w * x{i} * (double)(*eval{i});''' + x[*t_f+{s5}nF] += w * x{i} * (double)(*m{i});''' s += f'''\ case {i + 1}: @@ -570,7 +570,7 @@ def add_isotropic_compartments() -> str: { int id = (long)ptr; int offset; - uint32_t *eval0, *eval1, *eval2, *eval3, *eval4, *eval5, *eval6, *eval7, *eval8, *eval9, *eval10, *eval11, *eval12, *eval13, *eval14, *eval15, *eval16, *eval17, *eval18, *eval19; + uint32_t *m0, *m1, *m2, *m3, *m4, *m5, *m6, *m7, *m8, *m9, *m10, *m11, *m12, *m13, *m14, *m15, *m16, *m17, *m18, *m19; double x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, w, YTmp; double *xPtr0, *xPtr1, *xPtr2, *xPtr3, *xPtr4, *xPtr5, *xPtr6, *xPtr7, *xPtr8, *xPtr9, *xPtr10, *xPtr11, *xPtr12, *xPtr13, *xPtr14, *xPtr15, *xPtr16, *xPtr17, *xPtr18, *xPtr19; double *YPtr, *YPtrEnd; @@ -660,7 +660,7 @@ def add_isotropic_compartments() -> str: void COMMIT_At( double *_vIN, double *_vOUT, int _nS, int _ndirs, int _nF, int _n, int _nE, int _nV, - uint32_t *_ICf, uint32_t *_ICeval, uint32_t *_ICv, uint16_t *_ICo, float *_ICl, + uint32_t *_ICf, uint32_t *_ICm, uint32_t *_ICv, uint16_t *_ICo, float *_ICl, uint32_t *_ECv, uint16_t *_ECo, uint32_t *_ISOv, float *_wmrSFP, float *_wmhSFP, float *_isoSFP, @@ -679,7 +679,7 @@ def add_isotropic_compartments() -> str: Y = _vIN; ICf = _ICf; - ICeval = _ICeval; + ICm = _ICm; ICv = _ICv; ICo = _ICo; ICl = _ICl; @@ -722,7 +722,7 @@ def add_intracellular_compartments() -> str: while( t_v != t_vEnd ) { - x0 = x[*t_f] * (double)(ICeval[*t_f]); + x0 = x[*t_f] * (double)(ICm[*t_f]); if ( x0 != 0 ) Y[*t_v] += (double)(*t_l) * x0; t_f++; @@ -757,7 +757,7 @@ def add_isotropic_compartments() -> str: void* COMMIT_A__block_nolut( void *ptr ) { int id = (long)ptr; - uint32_t *eval0; + uint32_t *m0; double x0; double *xPtr; uint32_t *t_v, *t_vEnd, *t_f; @@ -779,7 +779,7 @@ def add_commit_a_nolut() -> str: void COMMIT_A_nolut( double *_vIN, double *_vOUT, int _nF, - uint32_t *_ICf, uint32_t *_ICeval, uint32_t *_ICv, float *_ICl, + uint32_t *_ICf, uint32_t *_ICm, uint32_t *_ICv, float *_ICl, uint32_t *_ISOv, uint32_t* _ICthreads, uint32_t* _ISOthreads, uint32_t _nISO, uint32_t _nThreads @@ -791,7 +791,7 @@ def add_commit_a_nolut() -> str: Y = _vOUT; ICf = _ICf; - ICeval = _ICeval; + ICm = _ICm; ICv = _ICv; ICl = _ICl; ISOv = _ISOv; @@ -828,7 +828,7 @@ def add_intracellular_compartments() -> str: { // in this case, I need to walk throug because the segments are ordered in "voxel order" if ( *t_t == id ) - x[*t_f] += (double)(*t_l) * Y[*t_v] * (double)(ICeval[*t_f]); + x[*t_f] += (double)(*t_l) * Y[*t_v] * (double)(ICm[*t_f]); t_t++; t_f++; t_v++; @@ -858,7 +858,7 @@ def add_isotropic_compartments() -> str: { int id = (long)ptr; double *xPtr; - uint32_t *eval0; + uint32_t *m0; uint32_t *t_v, *t_vEnd, *t_f; float *t_l; uint8_t *t_t;\n\n''' @@ -880,7 +880,7 @@ def add_commit_at_nolut() -> str: void COMMIT_At_nolut( double *_vIN, double *_vOUT, int _nF, int _n, - uint32_t *_ICf, uint32_t *_ICeval, uint32_t *_ICv, float *_ICl, + uint32_t *_ICf, uint32_t *_ICm, uint32_t *_ICv, float *_ICl, uint32_t *_ISOv, uint8_t* _ICthreadsT, uint32_t* _ISOthreadsT, uint32_t _nISO, uint32_t _nThreads @@ -893,7 +893,7 @@ def add_commit_at_nolut() -> str: Y = _vIN; ICf = _ICf; - ICeval = _ICeval; + ICm = _ICm; ICv = _ICv; ICl = _ICl; ISOv = _ISOv; From 4724edcad71fca18757417e37ad0a541f7d65b8e Mon Sep 17 00:00:00 2001 From: Ilaria Gabusi Date: Wed, 13 May 2026 14:50:34 +0200 Subject: [PATCH 7/7] Update debias documentation and condition for creating the mask ('debias_cond') --- commit/core.pyx | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/commit/core.pyx b/commit/core.pyx index e17149d..7b36001 100644 --- a/commit/core.pyx +++ b/commit/core.pyx @@ -1353,7 +1353,7 @@ cdef class Evaluation : logger.info( f'[ {format_time(time.time() - tr)} ]' ) - def fit( self, tol_fun=1e-3, tol_x=1e-6, max_iter=100, x0=None, confidence_map_filename=None, confidence_map_rescale=False, debias=False, thr_debias=1e-15 ): + def fit( self, tol_fun=1e-3, tol_x=1e-6, max_iter=100, x0=None, confidence_map_filename=None, confidence_map_rescale=False, debias=False, debias_cond=0.0 ): """Fit the model to the data. Parameters @@ -1375,9 +1375,14 @@ cdef class Evaluation : range [0.0,1.0]. Only the voxels considered in the mask will be affected. (default : False) debias : boolean - If true, the debias step will be performed (default : False) - thr_debias : float - Threshold for the debiasing step (default : 1e-15) + If true, a debiasing step will be performed after the main fitting + procedure. Highly suggested when using a regularisation. (default : False) + debias_cond : float + Condition used to select the coefficients to be debiased. + A second fit (without regularisation) will be performed on the reduced + problem defined by as a subset of the original linear operator. This is obtained + by selecting only the columns corresponding to estimated coefficients greater + than debias_cond in the main fitting procedure. (default : 0.0) """ if self.niiDWI is None : logger.error( 'Data not loaded; call "load_data()" first' ) @@ -1487,17 +1492,17 @@ cdef class Evaluation : temp_verb = self.verbose logger.subinfo('') logger.info( f'Running debias' ) - logger.subinfo( f'Creating mask for IC compartment with threshold={thr_debias:.2e}', indent_lvl=1, indent_char='*' ) + logger.subinfo( f'Creating mask for IC compartment with the condition: <={debias_cond:.2e}', indent_lvl=1, indent_char='*' ) self.set_verbose(0) offset = self.DICTIONARY['IC']['nSTR'] * self.KERNELS['wmr'].shape[0] xic = self.x[:offset] mask = np.ones(offset, dtype=np.uint32) - mask[xic