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_equalizer=0, deriv_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,25 @@ 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_equalizer: float
+            equalizer parameter of the Tikhonov regularization term
+        deriv_matrix: string
+            derivative matrix of the Tikhonov regularization term
+            If None (default), no regularization term is added to the model
+            If using this option, tikhonov_equalizer must be positive
         """
         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_equalizer < 0:
+            ERROR( 'Invalid value for Tikhonov equalizer parameter; value must be positive or zero' )
+        if tikhonov_equalizer > 0 and deriv_matrix == None:
+            ERROR( 'Tikhonov equalizer term given but derivative matrix was not selected; add "deriv_matrix" parameter in "build_operator()", valid options are \'L1\' (first  derivative with free boundary conditions), \'L2\' (second derivative with free boundary conditions), \'L1z\' (first  derivative with zero boundary conditions) and \'L2z\' (second derivative with zero boundary conditions)' )
+        if tikhonov_equalizer > 0 and deriv_matrix!='L1' and deriv_matrix!='L2' and deriv_matrix!='L1z' and deriv_matrix!='L2z':
+            ERROR( 'Invalid derivative matrix selection for regularization term; valid options are \'L1\' (first  derivative with free boundary conditions), \'L2\' (second derivative with free boundary conditions), \'L1z\' (first  derivative with zero boundary conditions) and \'L2z\' (second derivative with zero boundary conditions)' )
 
         if self.DICTIONARY['IC']['nF'] <= 0 :
             ERROR( 'No streamline found in the dictionary; check your data' )
@@ -710,7 +722,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_equalizer, deriv_matrix )
 
         LOG( '   [ %.1f seconds ]' % ( time.time() - tic ) )
 
@@ -724,7 +736,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_equalizer > 0:
+            if self.A.deriv_matrix == 'L1':
+                yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]-1, dtype=np.float64)
+            elif self.A.deriv_matrix == 'L2':
+                yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]-2, dtype=np.float64)
+            elif self.A.deriv_matrix == 'L1z':
+                yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]+1, dtype=np.float64)
+            elif self.A.deriv_matrix == 'L2z':
+                yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]  , dtype=np.float64)
+            else:
+                ERROR( 'Invalid derivative matrix selection for regularization term; valid options are \'L1\' (first  derivative with free boundary conditions), \'L2\' (second derivative with free boundary conditions), \'L1z\' (first  derivative with zero boundary conditions) and \'L2z\' (second derivative with zero boundary conditions)' )
+
+            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 ) :

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 COMMIT_L1(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void COMMIT_L2(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void COMMIT_L1z(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void COMMIT_L2z(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void COMMIT_L1t(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void COMMIT_L2t(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void COMMIT_L1zt(
+    int _nF, int _nIC, int _nV, int _nS, double _lambda,
+    double *_v_in, double *_v_out
+) nogil
+
+cdef extern void COMMIT_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_equalizer
+    cdef public char* deriv_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_equalizer=0, deriv_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_equalizer = tikhonov_equalizer # equalizer parameter of the Tikhonov regularization term
+        self.deriv_matrix       = deriv_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 deriv_matrix
+        if self.tikhonov_equalizer > 0.0:
+            if self.deriv_matrix == 0:
+                self.n1 = self.nV*self.nS + (self.nR-1)
+            elif self.deriv_matrix == 1:
+                self.n1 = self.nV*self.nS + (self.nR-2)
+            elif self.deriv_matrix == 2:
+                self.n1 = self.nV*self.nS + (self.nR+1)
+            else:
+                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_equalizer, self.deriv_matrix )
         C.adjoint = 1 - C.adjoint
         return C
 
@@ -188,4 +242,58 @@ cdef class LinearOperator :
                     self.ICthreadsT, self.ECthreadsT, self.ISOthreadsT
                 )
 
+        if self.tikhonov_equalizer > 0:
+            if not self.adjoint:
+                # DIRECT PRODUCT lambda*L*x
+                if self.deriv_matrix == 'L1':
+                    with nogil:
+                        COMMIT_L1(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_equalizer,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.deriv_matrix == 'L2':
+                    with nogil:
+                        COMMIT_L2(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_equalizer,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.deriv_matrix == 'L1z':
+                    with nogil:
+                        COMMIT_L1z(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_equalizer,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.deriv_matrix == 'L2z':
+                    with nogil:
+                        COMMIT_L2z(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_equalizer,
+                            &v_in[0], &v_out[0]
+                        )
+            else:
+                # INVERSE PRODUCT lambda*L'*y
+                if self.deriv_matrix == 'L1':
+                    with nogil:
+                        COMMIT_L1t(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_equalizer,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.deriv_matrix == 'L2':
+                    with nogil:
+                        COMMIT_L2t(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_equalizer,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.deriv_matrix == 'L1z':
+                    with nogil:
+                        COMMIT_L1zt(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_equalizer,
+                            &v_in[0], &v_out[0]
+                        )
+                elif self.deriv_matrix == 'L2z':
+                    with nogil:
+                        COMMIT_L2zt(
+                            self.nF, self.nR, self.nV, self.nS, self.tikhonov_equalizer,
+                            &v_in[0], &v_out[0]
+                        )
+
         return v_out

commit/operator/operator_noLUT.c

@@ -185,3 +185,62 @@ void COMMIT_At(
         pthread_join( threads[t], NULL );
     return;
 }
+
+void COMMIT_L(
+    int nF, int nIC, int nV, int nS, double regterm,
+    double *vIN, double *vOUT)
+{
+    /*for(int r = 0; r < nIC-1; r++){
+        for(int f = 0; f < nF; f++){
+            vOUT[nV*nS + r] += regterm*( -vIN[r*nF + f] + vIN[(r+1)*nF + f] );
+        }
+    }//*/
+}
+
+void COMMIT_Lt(
+    int nF, int nIC, int nV, int nS, double regterm,
+    double *vIN, double *vOUT)
+{
+    /*for(int f = 0; f < nF; f++){
+        vOUT[f] = -vIN[nV*nS];
+        vOUT[nF*(nIC-1) + f] = vIN[nV*nS + nIC-2];
+    }
+
+    for(int r = 0; r < nIC-2; r++){
+        for(int f = 0; f < nF; f++){
+            vOUT[nF*(r+1) + f] = vIN[nV*nS + r] + vIN[nV*nS + r+1];
+        }
+    }//*/
+}
+
+
+/*void COMMIT_L(
+    int nF, int nIC, int nV, int nS, double regterm,
+    double *vIN, double *vOUT)
+{
+    for(int f = 0; f < nF; f++){
+
+        vOUT[nV*nS] += regterm*( -2*vIN[f] + x[nF + f] );
+
+        for(int r = 1; r < nIC-1; r++){
+            vOUT[nV*nS + r] += regterm*( vIN[(r-1)*nF + f] -2*vIN[r*nF + f] + vIN[(r+1)*nF + f] );
+        }
+
+        vOUT[nV*nS + nIC - 1] += regterm*( vIN[(nIC-2)*nF + f] - 2*vIN[(nIC-1)*nF + f] );
+    }
+}
+
+void COMMIT_Lt(
+    int nF, int nIC, int nV, int nS, double regterm,
+    double *vIN, double *vOUT)
+{
+    for(int f = 0; f < nF; f++){
+        vOUT[f] += regterm*( -2*vIN[nV*nS] + vIN[nV*nS + 1] );
+
+        for (int r = 0; r < nIC; r++){
+            vOUT[r*nF + f] += regterm*( vIN[nV*nS + (r-1)] - 2*vIN[nV*nS + r] + vIN[nV*nS + (r+1)] );
+        }
+
+        vOUT[(nIC-1)*nF + f] += regterm*( vIN[nV*nS + (nIC-2)] - 2*vIN[nV*nS + (nIC-1)] );
+    }
+}//*/