From 0abdf927cecd07682b46f5fdb12d7a7ac341d7d0 Mon Sep 17 00:00:00 2001 From: fullbat Date: Mon, 24 Jun 2024 11:10:08 +0200 Subject: [PATCH] Refactor code and fix print statements, force numpy version to 1.26.4 --- .github/workflows/run_tests.yml | 2 ++ commit/core.pyx | 17 ----------------- commit/operator/operator_c.c | 14 +++++++------- commit/solvers.py | 3 --- setup.py | 8 ++++---- setup_operator.py | 8 ++------ 6 files changed, 15 insertions(+), 37 deletions(-) diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml index 16a2dbf..3607ace 100644 --- a/.github/workflows/run_tests.yml +++ b/.github/workflows/run_tests.yml @@ -47,6 +47,8 @@ jobs: python-version: ${{ matrix.config.py }} - name: Install dmri-commit run: pip install . --no-cache-dir + - name: Force numpy version + run: pip install numpy==1.26.4 - name: Run test id: run_test working-directory: tests diff --git a/commit/core.pyx b/commit/core.pyx index 994fc57..64eb34e 100755 --- a/commit/core.pyx +++ b/commit/core.pyx @@ -1368,28 +1368,16 @@ cdef class Evaluation : logger.subinfo('Recomputing coefficients', indent_lvl=1, indent_char='*', with_progress=True) x_debias = self.x.copy() - - logger.debug( f'positive values of x before debias: {np.sum(x_debias>0)}' ) - logger.debug( f'positive values of mask: {np.sum(mask>0)}' ) - x_debias[:offset1] *= mask x_debias[offset1:] = 0 - logger.debug( f"positive values of masked x: {np.sum(x_debias[:offset1]>0)}" ) - - logger.debug( f'Shape of y: {self.get_y().size} Number of non zero values in y before: {np.sum(self.get_y()>0)}' ) - y_mask = np.asarray(self.A.dot(x_debias)) - print(f"number of non zero values in y_mask before bin: {np.sum(y_mask>0)}") # binarize y_debias y_mask[y_mask<0] = 0 y_mask[y_mask>0] = 1 - print(f"number of non zero values in y_mask after bin: {np.sum(y_mask>0)}") self.debias_mask = y_mask - logger.debug( f'Shape of y: {self.get_y().size} Number of non zero values in y after: {np.sum(self.get_y()>0)}' ) - # print the first 10 non zero values of y_debias 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) @@ -1506,18 +1494,13 @@ cdef class Evaluation : ind_mask = np.where(self.debias_mask>0)[0] vox_mask = np.reshape( self.debias_mask[ind_mask], (nV,-1) ) - print(f'number of voxels in debias mask: {nV}') - print(f'ind_mask shape: {ind_mask.shape}') y_mea = np.reshape( self.get_y()[ind_mask], (nV,-1) ) y_est_ = np.asarray(self.A.dot(self.x)) y_est = np.reshape( y_est_[ind_mask], (nV,-1) ) - print(f"y_mea shape: {y_mea.shape}, y_est shape: {y_est.shape}") tmp = np.sqrt( np.mean((y_mea-y_est)**2,axis=1) ) - print(f'tmp shape: {tmp.shape}') - print(f'number of non-zero elements in tmp: {np.count_nonzero(tmp)}') logger.subinfo(f'RMSE: {tmp.mean():.3f} +/- {tmp.std():.3f}', indent_lvl=2, indent_char='-') tmp = np.sum(y_mea**2,axis=1) diff --git a/commit/operator/operator_c.c b/commit/operator/operator_c.c index 1694963..2e5dae1 100644 --- a/commit/operator/operator_c.c +++ b/commit/operator/operator_c.c @@ -10,7 +10,7 @@ double *x, *Y; uint32_t *ICthreads, *ECthreads, *ISOthreads; uint8_t *ICthreadsT; uint32_t *ECthreadsT, *ISOthreadsT; -uint32_t *ICf, *ICv, *ECv, *ISOv, *ICeval; +uint32_t *ICf, *ICeval, *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; @@ -49,8 +49,8 @@ void* COMMIT_A__block( void *ptr ) while (t_v != t_vEnd) { xPtr0 = x + (*t_f); - // eval0 = ICeval + *t_f; - x0 = *xPtr0 * (double)(ICeval[*t_f]); + eval0 = ICeval + *t_f; + x0 = *xPtr0 * (double)(*eval0); if (x0 != 0) { YPtr = Y + nS * (*t_v); @@ -3907,14 +3907,14 @@ void* COMMIT_At__block( void *ptr ) YTmp = *YPtr; SFP0ptr = wmrSFP0 + offset; x0 = (*SFP0ptr++) * YTmp; - // eval0 = ICeval + *t_f; + eval0 = ICeval + *t_f; while (++YPtr != YPtrEnd) { YTmp = *YPtr; x0 += (*SFP0ptr++) * YTmp; } - x[*t_f] += w * x0 * (double)(ICeval[*t_f]); + x[*t_f] += w * x0 * (double)(*eval0); } t_f++; t_v++; @@ -8060,7 +8060,7 @@ void* COMMIT_At__block( void *ptr ) void COMMIT_At( int _nF, int _n, int _nE, int _nV, int _nS, int _ndirs, double *_vIN, double *_vOUT, - uint32_t *_ICf, unsigned int *_ICeval, uint32_t *_ICv, uint16_t *_ICo, float *_ICl, + uint32_t *_ICf, uint32_t *_ICeval, uint32_t *_ICv, uint16_t *_ICo, float *_ICl, uint32_t *_ECv, uint16_t *_ECo, uint32_t *_ISOv, float *_wmrSFP, float *_wmhSFP, float *_isoSFP, @@ -9032,4 +9032,4 @@ void COMMIT_At_nolut( for(t=0; t<_nThreads ; t++) pthread_join( threads[t], NULL ); return; -} \ No newline at end of file +} diff --git a/commit/solvers.py b/commit/solvers.py index 2257db2..86fce26 100755 --- a/commit/solvers.py +++ b/commit/solvers.py @@ -273,9 +273,6 @@ def fista(y, A, At, omega, prox, sqrt_W=None, tol_fun=1e-4, tol_x=1e-6, max_iter else: res = A.dot(xhat) - y grad = np.asarray(At.dot(res)) - - print("grad", grad[grad!=0][:50]) - print("res", res[res!=0][:50]) prox( xhat, 1.0 ) reg_term = omega( xhat ) diff --git a/setup.py b/setup.py index 9b98d0f..1c4d09b 100644 --- a/setup.py +++ b/setup.py @@ -53,10 +53,10 @@ def run(self): build_ext.finalize_options(self) build_ext.run(self) -# # generate the operator_c.c file -# sys.path.insert(0, os.path.dirname(__file__)) -# from setup_operator import write_operator_c_file -# write_operator_c_file() +# generate the operator_c.c file +sys.path.insert(0, os.path.dirname(__file__)) +from setup_operator import write_operator_c_file +write_operator_c_file() # create the 'build' directory if not os.path.exists('build'): diff --git a/setup_operator.py b/setup_operator.py index 4697c2e..97b8ecc 100644 --- a/setup_operator.py +++ b/setup_operator.py @@ -701,8 +701,7 @@ def add_intracellular_compartments() -> str: while( t_v != t_vEnd ) { - eval0 = ICeval[*t_f]; - x0 = x[*t_f]* (double)(*eval0); + x0 = x[*t_f] * (double)(ICeval[*t_f]); if ( x0 != 0 ) Y[*t_v] += (double)(*t_l) * x0; t_f++; @@ -737,7 +736,6 @@ def add_isotropic_compartments() -> str: void* COMMIT_A__block_nolut( void *ptr ) { int id = (long)ptr; - uint32_t *eval0; double x0; double *xPtr; uint32_t *t_v, *t_vEnd, *t_f; @@ -806,8 +804,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 ) - eval0 = ICeval + *t_f; - x[*t_f] += (double)(*t_l) * Y[*t_v]* (double)(*eval0); + x[*t_f] += (double)(*t_l) * Y[*t_v] * (double)(ICeval[*t_f]); t_t++; t_f++; t_v++; @@ -837,7 +834,6 @@ def add_isotropic_compartments() -> str: { int id = (long)ptr; double *xPtr; - uint32_t *eval0; uint32_t *t_v, *t_vEnd, *t_f; float *t_l; uint8_t *t_t;\n\n'''