Skip to content

Commit

Permalink
testing Tikhonov implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
@fullbat
fullbat committed Oct 17, 2024
1 parent 21b4467 commit 38a8f83
Show file tree
Hide file tree
Showing 3 changed files with 29 additions and 51 deletions.
15 changes: 2 additions & 13 deletions commit/core.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -757,15 +757,8 @@ cdef class Evaluation :

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:
yL = np.zeros(y.shape[0] + self.KERNELS['wmr'].shape[0]-1, dtype=np.float64)
print( f'yl shape: {yL.shape}' )
print( f'y shape: {y.shape}' )
print( f'wmr shape: {self.KERNELS["wmr"].shape}' )
yL[:y.shape[0]] = y
return yL
else:
return y

return y


def set_regularisation(self, regularisers=(None, None, None), lambdas=(None, None, None), is_nonnegative=(True, True, True), params=(None, None, None)):
Expand Down Expand Up @@ -1204,12 +1197,8 @@ cdef class Evaluation :
# if dictISO_params is not None and 'coeff_weights' in dictISO_params:
# regularisation['lambdaISO_max'] = compute_lambda_max_lasso(regularisation['startISO'], regularisation['sizeISO'], dictISO_params['coeff_weights'])
# else:
print('before compute_lambda_max_lasso')
regularisation['lambdaISO_max'] = compute_lambda_max_lasso(regularisation['startISO'], regularisation['sizeISO'], np.ones(regularisation['sizeISO'], dtype=np.float64))
print('here AAAAAAAAAA')
print('lambdaISO_max', regularisation['lambdaISO_max'])
regularisation['lambdaISO'] = regularisation['lambdaISO_perc'] * regularisation['lambdaISO_max']
print('here')

# print
if regularisation['regISO'] is not None:
Expand Down
53 changes: 24 additions & 29 deletions commit/operator/operator.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,8 +139,6 @@ def precompute_neighbours(mask, mask_ix, mask_iy, mask_iz):
# The neighbour indices array
neighbours = neighbor_local_indices.astype(np.uint32)

print(f"neighbours: {neighbours[:10]}")
print(f"indptr: {indptr[:10]}")

return neighbours, indptr

Expand Down Expand Up @@ -173,16 +171,19 @@ cdef class LinearOperator :
cdef float* LUT_EC
cdef float* LUT_ISO

cdef unsigned int* ICthreads
cdef unsigned int* ECthreads
cdef unsigned int* ISOthreads
cdef unsigned int* ICthreads
cdef unsigned int* ECthreads
cdef unsigned int* ISOthreads

cdef unsigned char* ICthreadsT
cdef unsigned int* ECthreadsT
cdef unsigned int* ISOthreadsT
cdef unsigned char* ICthreadsT
cdef unsigned int* ECthreadsT
cdef unsigned int* ISOthreadsT

cdef unsigned int* neighbours
cdef unsigned int* indptr
cdef unsigned int* neighbours
cdef unsigned int* indptr

cdef public neigh
cdef public neigh_nptr

def __init__( self, DICTIONARY, KERNELS, THREADS, tikhonov_lambda=0, nolut=False ) :
"""Set the pointers to the data structures used by the C code."""
Expand Down Expand Up @@ -210,12 +211,7 @@ cdef class LinearOperator :

self.adjoint = 0 # direct of inverse product

# set shape of the operator according to tikhonov_matrix
if self.tikhonov_lambda > 0:
self.n1 = self.nV*self.nS + (self.nR-1)
else:
self.n1 = self.nV*self.nS

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
Expand Down Expand Up @@ -258,19 +254,18 @@ cdef class LinearOperator :
cdef unsigned int [::1] ISOthreadsT = THREADS['ISOt']
self.ISOthreadsT = &ISOthreadsT[0]

# precompute neighbours
print(f"nV: {DICTIONARY['nV']}")
n, np = precompute_neighbours(DICTIONARY['MASK'], DICTIONARY['MASK_ix'], DICTIONARY['MASK_iy'], DICTIONARY['MASK_iz'])
cdef unsigned int [::1] neighbours = n
self.neighbours = &neighbours[0]
cdef unsigned int [::1] indptr = np
self.indptr = &indptr[0]

# print neigbours and indptr as memoryview
for i in range(10):
print(f"neighbours: {self.neighbours[i]}")
for i in range(10):
print(f"indptr: {self.indptr[i]}")
cdef unsigned int [::1] neighbours
cdef unsigned int [::1] indptr

if tikhonov_lambda > 0:
# precompute neighbours
self.neigh, self.neigh_nptr = precompute_neighbours(DICTIONARY['MASK'], DICTIONARY['MASK_ix'], DICTIONARY['MASK_iy'], DICTIONARY['MASK_iz'])

neighbours = self.neigh
self.neighbours = &neighbours[0]

indptr = self.neigh_nptr
self.indptr = &indptr[0]


@property
Expand Down
12 changes: 3 additions & 9 deletions commit/operator/operator_c.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8604,7 +8604,6 @@ void* Tikhonov_block(void *ptr)
// Isotropic compartments
if (nISO > 0)
{
printf("Tikhonov_block\n");
t_v = ISOv + ISOthreads[id];
t_vEnd = ISOv + ISOthreads[id + 1];
xPtr = x + nF + ISOthreads[id]; // Pointer to x values for the thread
Expand All @@ -8614,16 +8613,13 @@ void* Tikhonov_block(void *ptr)
voxel_index = *t_v++;
neigh_s = neighbour_ptr[voxel_index];
neigh_e = neighbour_ptr[voxel_index + 1];
printf("voxel_index: %d\n", voxel_index);
printf("neigh_s: %d\n", neigh_s);
printf("neigh_e: %d\n", neigh_e);
for (int j = neigh_s; j < neigh_e; j++)
{
// printf("j: %d\n", j);
// printf("neighbours[j]: %d\n", neighbours[j]);
neighbor_index = neighbours[j];
// printf("x[voxel_index]: %f\n", x[voxel_index]);
diff = x[voxel_index] - x[neighbor_index];
diff = xPtr[voxel_index] - xPtr[neighbor_index];

// Update the gradient with respect to x_i
Y[voxel_index] += 2.0 * lambda * diff;
Expand Down Expand Up @@ -8684,7 +8680,6 @@ void* Tikhonov_t_block(void *ptr)
// Isotropic compartments
if (nISO > 0)
{
printf("Tikhonov_t_block\n");
t_v = ISOv + ISOthreadsT[id];
t_vEnd = ISOv + ISOthreadsT[id + 1];
xPtr = x + nF + ISOthreadsT[id]; // Pointer to x values for the thread
Expand All @@ -8695,11 +8690,10 @@ void* Tikhonov_t_block(void *ptr)
voxel_index = t_v[i];
neigh_s = neighbour_ptr[voxel_index];
neigh_e = neighbour_ptr[voxel_index + 1];
// printf("voxel_index: %d\n", voxel_index);
// printf("neigh_s: %d\n", neigh_s);
// printf("neigh_e: %d\n", neigh_e);

for (int j = neigh_s; j < neigh_e; j++)
{

(*xPtr) += 2.0 * lambda * (Y[voxel_index] - Y[neighbours[j]]);
}
xPtr++;
Expand Down

0 comments on commit 38a8f83

Please sign in to comment.