Add Tikhonov Regularization

.gitignore

@@ -15,4 +15,4 @@ __pycache__/
 *.cpp
 dist/
 
-trk2dictionary.c
+trk2dictionary.c

CHANGELOG.md

@@ -2,6 +2,16 @@
 # Change Log
 All notable changes to COMMIT will be documented in this file.
 
+## [1.5.0] - 2021-02-16
+
+### Added
+- core.pyx: Add to the function build_operator the parameter tikhonov_equalizer and deriv_matrix
+- operator.pyx: Extend Ax and A'y multiplications when Tikhonov regularization is enabled
+- operator_withLUC.c: Add C functions to extend Ax and A'y multiplications when Tikhonov regularization is enabled
+
+### Changed
+- core.pyx: The function get_y returns a larger vector y when Tikhonov regularization is enabled
+
 ## [1.4.5] - 2021-02-08
 
 ### Fixed

commit/core.pyx

@@ -639,7 +639,7 @@ cdef class Evaluation :
         LOG( '   [ %.1f seconds ]' % ( time.time() - tic ) )
 
 
-    def build_operator( self, build_dir=None ) :
+    def build_operator( self, build_dir=None, tikhonov_lambda=0, tikhonov_matrix=None ) :
         """Compile/build the operator for computing the matrix-vector multiplications by A and A'
         using the informations from self.DICTIONARY, self.KERNELS and self.THREADS.
         NB: needs to call this function to update pointers to data structures in case
@@ -652,13 +652,24 @@ cdef class Evaluation :
             If None (default), they will end up in the .pyxbld directory in the user’s home directory.
             If using this option, it is recommended to use a temporary directory, quit your python 
                 console between each build, and delete the content of the temporary directory.
+        tikhonov_lambda: float
+            Tikhonov lambda
+            If a positive value is given, tikhonov_matrix must not be None
+        tikhonov_matrix: string
+            Tikhonov matrix
         """
         if self.DICTIONARY is None :
             ERROR( 'Dictionary not loaded; call "load_dictionary()" first' )
         if self.KERNELS is None :
             ERROR( 'Response functions not generated; call "generate_kernels()" and "load_kernels()" first' )
         if self.THREADS is None :
             ERROR( 'Threads not set; call "set_threads()" first' )
+        if tikhonov_lambda < 0:
+            ERROR( 'Invalid lambda for Tikhonov regularization; value must be positive or zero' )
+        if tikhonov_lambda > 0 and tikhonov_matrix == None:
+            ERROR( 'Tikhonov lambda given but Tikhonov matrix was not selected; add "tikhonov_matrix" parameter in "build_operator()"' )
+        if tikhonov_lambda > 0 and tikhonov_matrix!='L1' and tikhonov_matrix!='L2' and tikhonov_matrix!='L1z' and tikhonov_matrix!='L2z':
+            ERROR( 'Invalid matrix selection for Tikhonov regularization term; check "tikhonov_matrix" parameter in "build_operator()"' )
 
         if self.DICTIONARY['IC']['nF'] <= 0 :
             ERROR( 'No streamline found in the dictionary; check your data' )
@@ -710,7 +721,7 @@ cdef class Evaluation :
             else :
                 reload( sys.modules['commit.operator.operator'] )
 
-        self.A = sys.modules['commit.operator.operator'].LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS )
+        self.A = sys.modules['commit.operator.operator'].LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS, tikhonov_lambda, tikhonov_matrix )
 
         LOG( '   [ %.1f seconds ]' % ( time.time() - tic ) )
 
@@ -724,7 +735,26 @@ cdef class Evaluation :
             ERROR( 'Dictionary not loaded; call "load_dictionary()" first' )
         if self.niiDWI is None :
             ERROR( 'Data not loaded; call "load_data()" first' )
-        return self.niiDWI_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'], : ].flatten().astype(np.float64)
+
+        y = self.niiDWI_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'], : ].flatten().astype(np.float64)
+
+        # extend y for the tikhonov regularization term
+        if self.A.tikhonov_lambda > 0:
+            if self.A.tikhonov_matrix == 'L1':
+                yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]-1, dtype=np.float64)
+            elif self.A.tikhonov_matrix == 'L2':
+                yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]-2, dtype=np.float64)
+            elif self.A.tikhonov_matrix == 'L1z':
+                yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]+1, dtype=np.float64)
+            elif self.A.tikhonov_matrix == 'L2z':
+                yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]  , dtype=np.float64)
+            else:
+                ERROR( 'Invalid matrix selection for Tikhonov regularization term; check "tikhonov_matrix" parameter in "build_operator()"' )
+
+            yL[0:y.shape[0]] = y
+            return yL
+        else:
+            return y
 
 
     def fit( self, tol_fun=1e-3, tol_x=1e-6, max_iter=100, verbose=True, x0=None, regularisation=None ) :
@@ -852,6 +882,7 @@ cdef class Evaluation :
         nF = self.DICTIONARY['IC']['nF']
         nE = self.DICTIONARY['EC']['nE']
         nV = self.DICTIONARY['nV']
+        nS = self.KERNELS['iso'].shape[1]
         norm_fib = np.ones( nF )
         # x is the x of the original problem
         # self.x is the x preconditioned
@@ -884,7 +915,7 @@ cdef class Evaluation :
         niiMAP_hdr['descrip'] = 'Created with COMMIT %s'%self.get_config('version')
 
         y_mea = np.reshape( self.niiDWI_img[ self.DICTIONARY['MASK_ix'], self.DICTIONARY['MASK_iy'], self.DICTIONARY['MASK_iz'], : ].flatten().astype(np.float32), (nV,-1) )
-        y_est = np.reshape( self.A.dot(self.x), (nV,-1) ).astype(np.float32)
+        y_est = np.reshape( self.A.dot(self.x)[0:int(nV*nS)], (nV,-1) ).astype(np.float32)
 
         print( '\t\t- RMSE...  ', end='' )
         sys.stdout.flush()

commit/operator/operator.pyx

@@ -26,6 +26,45 @@ cdef extern void COMMIT_At(
     unsigned char *_ICthreadsT, unsigned int *_ECthreadsT, unsigned int *_ISOthreadsT
 ) nogil
 
+cdef extern void Tikhonov_L1(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void Tikhonov_L2(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void Tikhonov_L1z(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void Tikhonov_L2z(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void Tikhonov_L1t(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void Tikhonov_L2t(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void Tikhonov_L1zt(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void Tikhonov_L2zt(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
 
 
 cdef class LinearOperator :
@@ -35,6 +74,8 @@ cdef class LinearOperator :
     """
     cdef int nS, nF, nR, nE, nT, nV, nI, n, ndirs
     cdef public int adjoint, n1, n2
+    cdef public float tikhonov_lambda
+    cdef public tikhonov_matrix
 
     cdef DICTIONARY
     cdef KERNELS
@@ -61,20 +102,22 @@ cdef class LinearOperator :
     cdef unsigned int*   ISOthreadsT
 
 
-    def __init__( self, DICTIONARY, KERNELS, THREADS ) :
+    def __init__( self, DICTIONARY, KERNELS, THREADS, tikhonov_lambda=0, tikhonov_matrix=None ) :
         """Set the pointers to the data structures used by the C code."""
         self.DICTIONARY = DICTIONARY
         self.KERNELS    = KERNELS
         self.THREADS    = THREADS
 
-        self.nF         = DICTIONARY['IC']['nF']    # number of FIBERS
-        self.nR         = KERNELS['wmr'].shape[0]   # number of FIBER RADII
-        self.nE         = DICTIONARY['EC']['nE']    # number of EC segments
-        self.nT         = KERNELS['wmh'].shape[0]   # number of EC TORTUOSITY values
-        self.nV         = DICTIONARY['nV']          # number of VOXELS
-        self.nI         = KERNELS['iso'].shape[0]   # number of ISO contributions
-        self.n          = DICTIONARY['IC']['n']     # numbner of IC segments
-        self.ndirs      = KERNELS['wmr'].shape[1]   # number of directions
+        self.nF         = DICTIONARY['IC']['nF']     # number of FIBERS
+        self.nR         = KERNELS['wmr'].shape[0]    # number of FIBER RADII
+        self.nE         = DICTIONARY['EC']['nE']     # number of EC segments
+        self.nT         = KERNELS['wmh'].shape[0]    # number of EC TORTUOSITY values
+        self.nV         = DICTIONARY['nV']           # number of VOXELS
+        self.nI         = KERNELS['iso'].shape[0]    # number of ISO contributions
+        self.n          = DICTIONARY['IC']['n']      # numbner of IC segments
+        self.ndirs      = KERNELS['wmr'].shape[1]    # number of directions
+        self.tikhonov_lambda = tikhonov_lambda # equalizer parameter of the Tikhonov regularization term
+        self.tikhonov_matrix       = tikhonov_matrix       # derivative matrix of the Tikhonov regularization term
 
         if KERNELS['wmr'].size > 0 :
             self.nS = KERNELS['wmr'].shape[2]       # number of SAMPLES
@@ -85,7 +128,18 @@ cdef class LinearOperator :
 
         self.adjoint    = 0                         # direct of inverse product
 
-        self.n1 = self.nV*self.nS
+        # set shape of the operator according to tikhonov_matrix
+        if self.tikhonov_lambda > 0:
+            if self.tikhonov_matrix == 'L1':
+                self.n1 = self.nV*self.nS + (self.nR-1)
+            elif self.tikhonov_matrix == 'L2':
+                self.n1 = self.nV*self.nS + (self.nR-2)
+            elif self.tikhonov_matrix == 'L1z':
+                self.n1 = self.nV*self.nS + (self.nR+1)
+            elif self.tikhonov_matrix == 'L2z':
+                self.n1 = self.nV*self.nS + (self.nR)
+        else:
+            self.n1 = self.nV*self.nS
         self.n2 = self.nR*self.nF + self.nT*self.nE + self.nI*self.nV
 
         # get C pointers to arrays in DICTIONARY
@@ -131,7 +185,7 @@ cdef class LinearOperator :
     @property
     def T( self ) :
         """Transpose of the explicit matrix."""
-        C = LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS )
+        C = LinearOperator( self.DICTIONARY, self.KERNELS, self.THREADS, self.tikhonov_lambda, self.tikhonov_matrix )
         C.adjoint = 1 - C.adjoint
         return C
 
@@ -188,4 +242,58 @@ cdef class LinearOperator :
                     self.ICthreadsT, self.ECthreadsT, self.ISOthreadsT
                 )
 
+        if self.tikhonov_lambda > 0:
+            if not self.adjoint:
+                # DIRECT PRODUCT lambda*L*x
+                if self.tikhonov_matrix == 'L1':
+                    with nogil:
+                        Tikhonov_L1(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.tikhonov_matrix == 'L2':
+                    with nogil:
+                        Tikhonov_L2(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.tikhonov_matrix == 'L1z':
+                    with nogil:
+                        Tikhonov_L1z(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.tikhonov_matrix == 'L2z':
+                    with nogil:
+                        Tikhonov_L2z(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda,
+                            &v_in[0], &v_out[0]
+                        )
+            else:
+                # INVERSE PRODUCT lambda*L'*y
+                if self.tikhonov_matrix == 'L1':
+                    with nogil:
+                        Tikhonov_L1t(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.tikhonov_matrix == 'L2':
+                    with nogil:
+                        Tikhonov_L2t(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.tikhonov_matrix == 'L1z':
+                    with nogil:
+                        Tikhonov_L1zt(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.tikhonov_matrix == 'L2z':
+                    with nogil:
+                        Tikhonov_L2zt(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_lambda,
+                            &v_in[0], &v_out[0]
+                        )
+
         return v_out

commit/operator/operator_noLUT.c

@@ -185,3 +185,5 @@ void COMMIT_At(
         pthread_join( threads[t], NULL );
     return;
 }
+
+//TODO: Add tikhonov regularization when no LUT is required