Fix occasional error in debias

CHANGELOG.md

@@ -1,6 +1,14 @@
 # Change Log
 ### All notable changes to `COMMIT` will be documented in this file.
 
+## `v2.4.3`<br>_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`<br>_2025-10-06_
 
 ### 🐛Fixed

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
 
@@ -724,9 +722,14 @@ cdef class Evaluation :
         logger.info( f'[ {format_time(time.time() - tic)} ]' )
 
 
-    def build_operator( self ) :
+    def build_operator( self, ic_debias_mask=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
+        ----------
+        ic_debias_mask : 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' )
@@ -743,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 ic_debias_mask is not None:
+            self.DICTIONARY['IC']['debias_mask'] = ic_debias_mask
+        else:
+            self.DICTIONARY['IC']['debias_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)} ]' )
 
@@ -760,9 +766,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):
@@ -1350,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 ) :
+    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_thr=0.0 ):
         """Fit the model to the data.
 
         Parameters
@@ -1371,6 +1374,15 @@ 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, a debiasing step will be performed after the main fitting
+            procedure. Highly suggested when using a regularisation. (default : False)
+        debias_thr : float
+            Threshold 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_thr in the main fitting procedure. (default : 0.0)
         """
         if self.niiDWI is None :
             logger.error( 'Data not loaded; call "load_data()" first' )
@@ -1380,9 +1392,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()
 
@@ -1473,39 +1487,31 @@ 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
-            logger.info( 'Running debias' )
+            logger.subinfo('')
+            logger.info( f'Running debias' )
+            logger.subinfo( f'Keeping columns in the IC compartment whose corresponding coefficient is greater than {debias_thr}', 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<0.000000000000001] = 0
-
-            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.set_regularisation()
-            self.set_verbose(temp_verb)
-
-            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
+            mask[xic<=debias_thr] = 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
-
-            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)
+            if np.sum(mask)==0:
+                self.set_verbose(temp_verb)
+                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:
+                # update the operator with the new mask for the IC compartment
+                self.build_operator(ic_debias_mask=mask)
+                # run fit on the debiased problem
+                self.set_regularisation()
+                self.set_verbose(temp_verb)
+                logger.subinfo('Recomputing coefficients', indent_lvl=1, indent_char='*', with_progress=True)
+                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')
@@ -1621,64 +1627,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) )

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 # debias mask for IC columns
         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']['debias_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],

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",