1 /***** This code is part of the AFNI software package, which is *****
2 ***** partly in the public domain and partly covered by the GPL. *****
3 ***** See https://afni.nimh.nih.gov/afni for more information. *****/
4
5 #include "mrilib.h"
6
7 static floatvec * THD_deconvolve_autopen( int npt , float *far ,
8 int minlag , int maxlag , float *kern,
9 int nbase , float *base[],
10 int meth , float *ccon , int dcon ,
11 int pencode, float *penfac_used );
12
13 /* 'AO' changes (from Aomp.h) are for compatibility with OpenMP */
14
15 AO_DEFINE_SCALAR(int,voxid) ;
THD_fitter_voxid(int i)16 void THD_fitter_voxid( int i ){ AO_VALUE(voxid) = i; }
17
18 /*------------------------------------------------------------------*/
19 /* Least squares fitting without constraints. (cf mri_matrix.c) */
20 /*------------------------------------------------------------------*/
21
new_lsqfit(int npt,float * far,int nref,float * ref[])22 static float * new_lsqfit( int npt , float *far ,
23 int nref , float *ref[] )
24 {
25 int jj ;
26 MRI_IMAGE *rmat,*pmat,*smat ;
27 float *rar;
28
29 ENTRY("new_lsqfit") ;
30
31 /* compute pseudo-inverse of matrix into pmat */
32
33 rmat = mri_new(npt,nref,MRI_float ); rar = MRI_FLOAT_PTR(rmat);
34 for( jj=0 ; jj < nref ; jj++ )
35 memcpy( rar+jj*npt , ref[jj] , sizeof(float)*npt ) ;
36 pmat = mri_matrix_psinv(rmat,NULL,0.0f) ;
37 mri_free(rmat) ;
38 if( pmat == NULL ) RETURN(NULL) ; /* should not happen */
39
40 /* create vector of data and multiply by pseudo-inverse */
41
42 rmat = mri_new_vol_empty( npt , 1 , 1 , MRI_float ) ;
43 mri_fix_data_pointer( far , rmat ) ;
44 smat = mri_matrix_mult( pmat , rmat ) ;
45 mri_free(pmat); mri_clear_data_pointer(rmat); mri_free(rmat);
46 if( smat == NULL ) RETURN(NULL) ; /* should not happen */
47
48 /* get pointer to results array and return it */
49
50 rar = MRI_FLOAT_PTR(smat);
51 mri_clear_data_pointer(smat); mri_free(smat);
52 RETURN(rar) ;
53 }
54
55 /*------------------------------------------------------------------*/
56
57 #undef GOOD_METH
58 #define GOOD_METH(m) ( (m)==1 || (m)==2 || (m)==-2 || (m)==-1 )
59
60 /**--- 05 Mar 2008: stuff to get the fitted model back ---**/
61
62 /* this variable is set once-and-for-all so doesn't need to be thread-ized */
63 static int do_fitv = 0 ; /* if 1, will compute fitts into fitv */
THD_fitter_do_fitts(int qq)64 void THD_fitter_do_fitts(int qq){ do_fitv = qq; }
65
66 AO_DEFINE_SCALAR(floatvec*,gfitv) ;
67 #if 0
68 static floatvec *gfitv= NULL ;
69 #endif
70
71 /* get the fitted time series vector (immediately after getting params) */
THD_retrieve_fitts(void)72 floatvec * THD_retrieve_fitts(void){ return AO_VALUE(gfitv); }
73
74 /* these variables are set once-and-for-all so don't need to be thread-ized */
75 static int use_rcmat = 0 ;
76
77 static float vthresh = 0.0f ; /* 18 May 2010: vector threshold */
THD_fitter_set_vthresh(float vvv)78 void THD_fitter_set_vthresh( float vvv ){
79 vthresh = (vvv > 0.0f && vvv < 0.1f) ? vvv : 0.0f ;
80 }
81
82 /*------------------------------------------------------------------*/
83 /* Fit the npt-long vector far[] to the nref vectors in ref[].
84 - meth = 1 ==> L1 fit
85 - meth = 2 ==> L2 fit
86 - meth = -2 ==> L2 fit with LASSO [11 Mar 2011]
87 - meth = -1 ==> L2 fit with SQRT(LASSO) [not yet implemented]
88 - ccon != NULL ==> ccon[i] is sign constraint on coef #i
89 ccon[i] = 0 == no constraint
90 > 0 == coef #i must be >= 0
91 < 0 == coef #i must be <= 0
92 - Output is vector of coefficients (nref of them).
93 - Output is NULL if some error transpired.
94 *//*----------------------------------------------------------------*/
95
THD_fitter(int npt,float * far,int nref,float * ref[],int meth,float * ccon)96 floatvec * THD_fitter( int npt , float *far ,
97 int nref, float *ref[],
98 int meth, float *ccon )
99 {
100 int ii,jj , nbad ;
101 float *qfit=NULL, val , qmax ;
102 floatvec *fv=NULL ;
103 float *qmag=NULL ; /* 18 May 2010 */
104
105 ENTRY("THD_fitter") ;
106
107 KILL_floatvec(AO_VALUE(gfitv)) ; /* 05 Mar 2008 */
108
109 /* check inputs for stupid users */
110
111 if( npt <= 1 || far == NULL ||
112 nref <= 0 || ref == NULL ||
113 (nref >= npt-1 && meth > 0) ) RETURN(NULL) ;
114
115 for( jj=0 ; jj < nref ; jj++ ) if( ref[jj] == NULL ) RETURN(NULL) ;
116
117 /*--- 08 Apr 2008: check if some columns are way small;
118 if so, excise them and solve smaller problem ---*/
119
120 qmag = (float *)malloc(sizeof(float)*nref) ;
121 for( qmax=0.0f,jj=0 ; jj < nref ; jj++ ){
122 for( val=0.0f,ii=0 ; ii < npt ; ii++ ) val += fabsf(ref[jj][ii]) ;
123 qmag[jj] = val ; if( val > qmax ) qmax = val ;
124 }
125 if( qmax == 0.0f ){ free(qmag); RETURN(NULL); } /* all zero?! */
126
127 qmax *= vthresh ; /* vthresh used to be fixed at 0.000333 [18 May 2010] */
128 for( nbad=jj=0 ; jj < nref ; jj++ ) if( qmag[jj] <= qmax ) nbad++ ;
129
130 if( nbad > 0 ){ /*-- must excise the tiny columns before solving --*/
131 int ngood = nref-nbad , nc ;
132 float **qref = (float **)calloc(sizeof(float *),ngood) ;
133 floatvec *qv ; float *qcon=NULL ;
134 if( ccon != NULL ) qcon = (float *)calloc(sizeof(float),ngood) ;
135 for( nc=ii=jj=0 ; jj < nref ; jj++ ){
136 if( qmag[jj] > qmax ){
137 if( ccon != NULL && ccon[jj] != 0.0f ){ qcon[ii] = ccon[jj] ; nc++ ; }
138 qref[ii++] = ref[jj] ;
139 }
140 }
141 if( nc == 0 ){ free(qcon) ; qcon = NULL ; }
142 qv = THD_fitter( npt , far , ngood , qref , meth , qcon ) ; /* recurse! */
143 if( qcon != NULL ) free((void *)qcon) ;
144 free((void *)qref) ;
145 if( qv == NULL ){ free(qmag); RETURN(NULL); } /* bad solve? */
146 MAKE_floatvec(fv,nref) ; /* will be all zero */
147 for( ii=jj=0 ; jj < nref ; jj++ ){ /* load results into output */
148 if( qmag[jj] > qmax ) fv->ar[jj] = qv->ar[ii++] ;
149 }
150 KILL_floatvec(qv) ; free(qmag) ; RETURN(fv) ;
151 }
152
153 /*------- actually solve now -------*/
154
155 switch( meth ){
156
157 default: free(qmag) ; RETURN(NULL) ; /* stupid user */
158
159 /*-- least squares with LASSO-ish things [experimental] --*/
160
161 case -1:{ /* not implemented yet */
162 floatvec *lfit ;
163 lfit = THD_sqrtlasso_L2fit( npt , far , nref , ref , NULL , ccon ) ;
164 if( lfit == NULL ){ free(qmag) ; RETURN(NULL) ; }
165 qfit = lfit->ar ; lfit->ar = NULL ; KILL_floatvec(lfit) ;
166 }
167 break ;
168
169 case -2:{
170 floatvec *lfit ;
171 lfit = THD_lasso_L2fit( npt , far , nref , ref , NULL , ccon ) ;
172 if( lfit == NULL ){ free(qmag) ; RETURN(NULL) ; }
173 qfit = lfit->ar ; lfit->ar = NULL ; KILL_floatvec(lfit) ;
174 }
175 break ;
176
177 /*-- straightforward least squares --*/
178
179 case 2:
180 if( ccon == NULL ){ /* unconstrained */
181 if( use_rcmat && nref > 9 ){
182 int first=1 ;
183 qfit = rcmat_lsqfit( npt, far, nref, ref ) ; /* 30 Dec 2008 */
184 if( qfit == NULL && first ){
185 WARNING_message("sparse matrix least squares fails: trying pseudo-inverse") ;
186 first = 0 ;
187 }
188 }
189 if( qfit == NULL )
190 qfit = new_lsqfit ( npt, far, nref, ref ) ;
191 } else { /* constrained */
192 qfit = (float *)malloc(sizeof(float)*nref); /* output array */
193 memcpy(qfit,ccon,sizeof(float)*nref) ;
194 val = cl2_solve( npt, nref, far, ref, qfit, 1 ) ; /* cf cl2.c */
195 if( val < 0.0f ){ free(qfit); qfit = NULL; } /* bad */
196 }
197 break ;
198
199 /*-- L1 fitting --*/
200
201 case 1:
202 qfit = (float *)malloc(sizeof(float)*nref) ; /* output array */
203 if( ccon != NULL ) memcpy(qfit,ccon,sizeof(float)*nref) ;
204 val = cl1_solve( npt,nref, far,ref, qfit, (ccon!=NULL) ); /* cf cl1.c */
205 if( val < 0.0f ){ free(qfit); qfit = NULL; } /* bad */
206 break ;
207 }
208
209 if( qfit == NULL ){ free(qmag); RETURN(NULL); } /* didn't get output array */
210
211 MAKE_floatvec(fv,nref) ; /* copy output array */
212 memcpy( fv->ar, qfit, sizeof(float)*nref ) ; /* into floatvec and */
213 free(qfit) ; /* free the trashola */
214 if( do_fitv ) /* compute fitts? */
215 AO_VALUE(gfitv) = THD_fitter_fitts( npt,fv,nref,ref,NULL ); /* 05 Mar 2008 */
216
217 free(qmag) ; RETURN(fv) ; /* return to caller */
218 }
219
220 /*-------------------------------------------------------------------------*/
221 /* Get the fitted time series (length npt), given the nref parameters
222 in fv (presumably from THD_fitter itself) and the reference functions
223 in ref[0..nref-1]. If far is not NULL, then its contents are
224 subtracted from the fit -- thus, giving the residuals.
225 *//*-----------------------------------------------------------------------*/
226
THD_fitter_fitts(int npt,floatvec * fv,int nref,float * ref[],float * far)227 floatvec * THD_fitter_fitts( int npt , floatvec *fv ,
228 int nref, float *ref[] , float *far )
229 {
230 int ii , jj ;
231 float sum , *qar , pval ;
232 floatvec *qv ;
233
234 ENTRY("THD_fitter_fitts") ;
235
236 if( npt < 1 || fv == NULL || fv->nar < nref ||
237 nref < 1 || ref == NULL ) RETURN(NULL) ;
238
239 MAKE_floatvec(qv,npt) ; qar = qv->ar ;
240 for( jj=0 ; jj < nref ; jj++ ){
241 pval = fv->ar[jj] ;
242 for( ii=0 ; ii < npt ; ii++ ) qar[ii] += ref[jj][ii] * pval ;
243 }
244 if( far != NULL )
245 for( ii=0 ; ii < npt ; ii++ ) qar[ii] -= far[ii] ;
246
247 RETURN(qv) ;
248 }
249
250 /*-------------------------------------------------------------------------*/
251
252 #if 0
253 # define ERREX(s) do { ERROR_message(s); RETURN(NULL); } while(0)
254 #else
255 # define ERREX(s) RETURN(NULL)
256 #endif
257
258 /* these fitted time series vectors are for deconvolution */
259
260 AO_DEFINE_SCALAR(int,nggfitvv) ;
261 AO_DEFINE_SCALAR(floatvec**,ggfitvv) ;
262 #if 0
263 static int nggfitvv = 0 ;
264 static floatvec **ggfitvv = NULL ;
265 #endif
266
267 /*------------------------------------------------------------------------*/
268 /* Fit a deconvolution model, using THD_fitter() as the workhorse.
269 - npt = length of time series
270 - far = time series data [0..npt-1]
271 - minlag = minimum lag into the past (can be negative) -- usually 0
272 - maxlag = maximum lag into the past
273 - kern = kernel to deconvolve = kern[0..maxlag-minlag]
274 - nbase = number of baseline parameters (can be 0)
275 - base = baseline reference functions (can be NULL if nbase==0)
276 - meth = 1 or 2 for L1 or L2 regresstion
277 - ccon = constraints on signs of baseline parameters
278 - dcon = constraint on sign of deconvolved time series s(t)
279 - pencode= penalty function for s(t) -- cf. 3dTfitter -help
280 - npfac = number of penalty factors to use
281 - pfac = array of penalty factors (if NULL, program makes them up)
282
283 Return value is an array of npfac floatvec-s, each of which has
284 npt+nbase values -- the first npt of which are s(t) and the last
285 nbase of which are the baseline parameters. A NULL return indicates
286 some bad input. An individual NULL floatvec indicates that particular
287 value of pfac caused the regression to fail for some hideous reason.
288
289 The purpose of evaluating multiple pfac-s is to implement the L-curve
290 methodology for choosing a pfac -- cf. THD_deconvolve_autopen()
291 *//*----------------------------------------------------------------------*/
292
THD_deconvolve_multipen(int npt,float * far,int minlag,int maxlag,float * kern,int nbase,float * base[],int meth,float * ccon,int dcon,int pencode,int npfac,float * pfac,float * pres,float * psiz)293 floatvec ** THD_deconvolve_multipen( int npt , float *far ,
294 int minlag , int maxlag , float *kern,
295 int nbase , float *base[],
296 int meth , float *ccon , int dcon ,
297 int pencode, int npfac , float *pfac,
298 float *pres, float *psiz )
299 {
300 int ii , jj , kk ;
301 float val,kernsum, *zar , *zcon=NULL ;
302 floatvec **fvv ;
303 int nref,nlag,npen,nplu ; float **ref ;
304 int p0,p1,p2,p3 , np0,np1,np2,np3 , rp0,rp1,rp2,rp3 , ipf ;
305 float penfac,fmed,fsig , *qfac , p1scl,p2scl,p3scl,p1fac,p2fac,p3fac ;
306
307 ENTRY("THD_deconvolve_multipen") ;
308
309 /* delete any leftover junk */
310
311 if( AO_VALUE(nggfitvv) > 0 ){
312 for( ii=0 ; ii < AO_VALUE(nggfitvv) ; ii++ ) KILL_floatvec(AO_VALUE(ggfitvv)[ii]) ;
313 free((void *)AO_VALUE(ggfitvv)) ; AO_VALUE(nggfitvv) = 0 ; AO_VALUE(ggfitvv) = NULL ;
314 }
315
316 /* check inputs for stupid users */
317
318 if( npt <= 3 || far == NULL || npfac < 1 ) ERREX("e1") ;
319 nlag = maxlag-minlag+1 ;
320 if( nlag <= 1 || kern == NULL || !GOOD_METH(meth) ) ERREX("e2") ;
321 if( minlag <= -npt+1 || maxlag >= npt-1 ) ERREX("e3") ;
322 if( nbase > 0 ){
323 if( base == NULL ) ERREX("e4") ;
324 for( jj=0 ; jj < nbase ; jj++ ) if( base[jj] == NULL ) ERREX("e5") ;
325 } else if( nbase < 0 ){ /* user = dumb as a brick */
326 nbase = 0 ;
327 }
328
329 for( kernsum=0.0f,jj=0 ; jj < nlag ; jj++ ){
330 val = kern[jj] ; kernsum += val*val ;
331 }
332 if( kernsum <= 0.0f ) ERREX("e6") ;
333 kernsum = sqrtf(kernsum) ;
334
335 /* count penalty equations */
336
337 p0 = (pencode & 1) ; /* which penalty functions are enabled */
338 p1 = (pencode & 2) ;
339 p2 = (pencode & 4) ;
340 p3 = (pencode & 8) ;
341 if( p0==0 && p1==0 && p2==0 && p3==0 ) p0 = 1 ; /* some penalty */
342 np0 = (p0) ? npt : 0 ; /* number of equations for each case */
343 np1 = (p1) ? npt-1 : 0 ;
344 np2 = (p2) ? npt-2 : 0 ;
345 np3 = (p3) ? npt-3 : 0 ;
346 rp0 = npt ; /* row offset for p0 functions */
347 rp1 = rp0 + np0 ; /* row offset for p1 functions */
348 rp2 = rp1 + np1 ; /* row offset for p2 functions */
349 rp3 = rp2 + np2 ; /* row offset for p3 functions */
350 npen = np0+np1+np2+np3 ; /* total number of penalty equations */
351
352 /* set scale level for negative pfac values */
353
354 #if 0
355 qmedmad_float ( npt , far , NULL , &fmed ) ; fmed *= 1.777f ;
356 meansigma_float( npt , far , &val , &fsig ) ;
357 if( fmed < fsig ) fmed = fsig ;
358 if( fmed == 0.0f ) fmed = fabsf(val) ; /* data is constant? */
359 if( fmed == 0.0f ) ERREX("e7") ; /* data is all zero? */
360 fmed = fmed * 1.777 / kernsum ;
361 #else
362 fmed = 0.333f * kernsum ;
363 #endif
364
365 #if 1
366 if( AFNI_yesenv("AFNI_TFITTER_VERBOSE") )
367 ININFO_message("default penalty scale factor=%g",fmed) ;
368 #endif
369
370 /* number of equations and number of references */
371
372 nplu = npt + npen ; /* number of equations */
373 nref = npt + nbase ; /* number of references */
374 if( nref > nplu ) ERREX("e8") ; /* only if user is really stupid */
375
376 /** make new reference vectors (that are nplu long) **/
377
378 ref = (float **)malloc(sizeof(float *)*nref) ;
379
380 /* deconvolution equations (rows and columns #0..npt-1) */
381
382 STATUS("setup deconvolution") ;
383
384 for( jj=0 ; jj < npt ; jj++ ){
385 ref[jj] = (float *)calloc(sizeof(float),nplu) ;
386 for( ii=0 ; ii < npt ; ii++ ){
387 kk = ii-jj ; if( kk < minlag ) continue; if( kk > maxlag ) break ;
388 ref[jj][ii] = kern[kk-minlag] ;
389 }
390 }
391
392 STATUS("setup baseline") ;
393
394 /* baseline equations, if any (columns #npt..nref-1, rows #0..npt-1) */
395
396 for( jj=0 ; jj < nbase ; jj++ ){
397 ref[jj+npt] = (float *)calloc(sizeof(float),nplu) ;
398 memcpy( ref[jj+npt] , base[jj] , sizeof(float)*npt ) ;
399 }
400
401 STATUS("copy data") ;
402
403 /* copy data (rows #0..npt-1) */
404
405 zar = (float *)calloc(sizeof(float),nplu) ;
406 memcpy(zar,far,sizeof(float)*npt) ; /* the rest will be zero */
407
408 /* copy constraints? */
409
410 STATUS("copy constraints?") ;
411
412 if( ccon != NULL || dcon != 0 ){
413 zcon = (float *)calloc(sizeof(float),nref) ;
414 if( dcon != 0 ) /* sign constraints on deconvolution output */
415 for( ii=0 ; ii < npt ; ii++ ) zcon[ii] = (float)dcon ;
416 if( nbase > 0 && ccon != NULL ) /* sign constraints on baseline params */
417 memcpy(zcon+npt,ccon,sizeof(float)*nbase) ;
418 }
419
420 /* make up some penalty factors if not supplied by user */
421
422 if( pfac == NULL ){
423 STATUS("make up penalty factors") ;
424 qfac = (float *)malloc(sizeof(float)*npfac) ;
425 if( npfac == 1 ){
426 qfac[0] = -1.0f ;
427 } else {
428 float pb,pt,dp ;
429 pt = sqrtf((float)npfac); pb = 1.0f/pt;
430 dp = powf(pt/pb,1.0f/(npfac-1.0f));
431 qfac[0] = -pb ;
432 for( ii=1 ; ii < npfac ; ii++ ) qfac[ii] = qfac[ii-1] * dp ;
433 }
434 } else {
435 qfac = pfac ;
436 }
437
438 /** loop over different penalty factors **/
439
440 fvv = (floatvec **)calloc(sizeof(floatvec *),npfac) ;
441
442 /* allocate space for multiple fit vectors, for multiple penalty factors */
443
444 if( do_fitv ){
445 AO_VALUE(ggfitvv) = (floatvec **)calloc(sizeof(floatvec *),npfac) ;
446 AO_VALUE(nggfitvv) = npfac ;
447 }
448
449 p1scl = AFNI_numenv("AFNI_TFITTER_P1SCALE"); if( p1scl <= 0.0f ) p1scl=1.0f ;
450 p2scl = AFNI_numenv("AFNI_TFITTER_P2SCALE"); if( p2scl <= 0.0f ) p2scl=1.0f ;
451 p3scl = AFNI_numenv("AFNI_TFITTER_P3SCALE"); if( p3scl <= 0.0f ) p3scl=1.0f ;
452
453 STATUS("loop over ipf") ;
454
455 for( ipf=0 ; ipf < npfac ; ipf++ ){
456
457 /* penalty eqations for deconv (columns #0..npt-1, rows #npt..nplu-1) */
458
459 penfac = qfac[ipf] ;
460 if( penfac == 0.0f ) penfac = -0.999f ;
461 if( penfac < 0.0f ) penfac = -penfac * fmed ;
462 p1fac = p1scl * penfac ;
463 p2fac = p2scl * penfac * 2.0f ;
464 p3fac = p3scl * penfac * 3.0f ;
465
466 if( p0 ){
467 STATUS(" setup p0") ;
468 for( jj=0 ; jj < npt ; jj++ ) ref[jj][rp0+jj] = penfac ;
469 }
470 if( p1 ){
471 STATUS(" setup p1") ;
472 for( jj=0 ; jj < npt-1 ; jj++ ) ref[jj][rp1+jj] = -p1fac ;
473 for( jj=1 ; jj < npt ; jj++ ) ref[jj][rp1+jj-1] = p1fac ;
474 }
475 if( p2 ){
476 STATUS(" setup p2") ;
477 val = 0.5f*p2fac ;
478 for( jj=0 ; jj < npt-2 ; jj++ ) ref[jj][rp2+jj] = -val ;
479 for( jj=1 ; jj < npt-1 ; jj++ ) ref[jj][rp2+jj-1] = p2fac ;
480 for( jj=2 ; jj < npt ; jj++ ) ref[jj][rp2+jj-2] = -val ;
481 }
482 if( p3 ){
483 STATUS(" setup p3") ;
484 val = 0.3333333f*p3fac ;
485 for( jj=0 ; jj < npt-3 ; jj++ ) ref[jj][rp3+jj] = -val ;
486 for( jj=1 ; jj < npt-2 ; jj++ ) ref[jj][rp3+jj-1] = p3fac ;
487 for( jj=2 ; jj < npt-1 ; jj++ ) ref[jj][rp3+jj-2] = -p3fac ;
488 for( jj=3 ; jj < npt ; jj++ ) ref[jj][rp3+jj-3] = val ;
489 }
490
491 /***** actually fit the parameters (deconvolution + baseline) *****/
492
493 fvv[ipf] = THD_fitter( nplu , zar , nref , ref , meth , zcon ) ;
494
495 /* fvv = fitted params; gfitv = fitted time series */
496
497 if( do_fitv && fvv[ipf] != NULL && AO_VALUE(gfitv) != NULL ){
498 STATUS("copy fitvec") ;
499 COPY_floatvec( AO_VALUE(ggfitvv)[ipf] , AO_VALUE(gfitv) ) ;
500 STATUS(" end copy fitvec") ;
501 }
502
503 if( pres != NULL && psiz != NULL && fvv[ipf] != NULL ){
504 floatvec *fitv ; float *rar , *sar , rsum=0.0f,ssum=0.0f ;
505 STATUS("computing pres + psiz") ;
506 fitv = THD_fitter_fitts( nplu,fvv[ipf] , nref,ref , zar ) ;
507 rar = fitv->ar ; sar = fvv[ipf]->ar ;
508 switch(meth){
509 case 1: case -1:
510 for( ii=0; ii < npt; ii++ ){
511 rsum += fabsf(rar[ii]); ssum += fabsf(sar[ii]);
512 }
513 if( p1 || p2 || p3 )
514 for( ii=1 ; ii < npt ; ii++ ) ssum += fabsf(sar[ii]-sar[ii-1]) ;
515 break ;
516 case 2: case -2:
517 for( ii=0; ii < npt; ii++ ){
518 rsum += rar[ii]*rar[ii] ; ssum += sar[ii]*sar[ii] ;
519 }
520 if( p1 || p2 || p3 )
521 for( ii=1 ; ii < npt ; ii++ ) ssum += SQR(sar[ii]-sar[ii-1]) ;
522 rsum = sqrtf(rsum) ; ssum = sqrtf(ssum) ;
523 break ;
524 }
525 KILL_floatvec(fitv) ; pres[ipf] = rsum ; psiz[ipf] = ssum ;
526 #if 1
527 if( AFNI_yesenv("AFNI_TFITTER_VERBOSE") )
528 ININFO_message("qfac=%g penfac=%g resid=%g norm=%g",
529 qfac[ipf],penfac,rsum,ssum) ;
530 #endif
531 }
532
533 STATUS("at end of ipf loop") ;
534 }
535
536 /* free the enslaved memory and return to the user */
537
538 STATUS("free-ing") ;
539
540 if( qfac != pfac ) free((void *)qfac) ;
541 if( zcon != NULL ) free((void *)zcon) ;
542 free((void *)zar) ;
543 for( jj=0 ; jj < nref ; jj++ ) free((void *)ref[jj]) ;
544 free((void *)ref) ;
545
546 STATUS("done") ;
547 RETURN(fvv) ;
548 }
549
550 /*-------------------------------------------------------------------------*/
551 /* The real work is outsourced to the function above, or the one below. */
552
THD_deconvolve(int npt,float * far,int minlag,int maxlag,float * kern,int nbase,float * base[],int meth,float * ccon,int dcon,int pencode,float penfac)553 floatvec * THD_deconvolve( int npt , float *far ,
554 int minlag , int maxlag , float *kern,
555 int nbase , float *base[],
556 int meth , float *ccon , int dcon ,
557 int pencode, float penfac )
558 {
559 floatvec **fvv , *fv=NULL ; float pfac=penfac ;
560
561 ENTRY("THD_deconvolve") ;
562
563 use_rcmat = !AFNI_noenv("AFNI_FITTER_RCMAT") ; /* 30 Dec 2008 */
564
565 if( pfac == -666.0f || pfac == 0.0f ){
566 fv = THD_deconvolve_autopen( npt , far , minlag , maxlag , kern ,
567 nbase , base , meth , ccon , dcon ,
568 pencode , NULL ) ;
569 } else {
570 fvv = THD_deconvolve_multipen( npt , far , minlag , maxlag , kern ,
571 nbase , base , meth , ccon , dcon ,
572 pencode , 1 , &pfac , NULL,NULL ) ;
573
574 if( fvv != NULL ){ fv = fvv[0]; free((void *)fvv); }
575 }
576
577 use_rcmat = 0 ;
578
579 RETURN(fv) ;
580 }
581
582 /*-------------------------------------------------------------------------*/
583 /* L curving. */
584
585 #undef NPFAC
586 #define NPFAC 11
587
fillerup(float bot,float top,int nval,float * val)588 static void fillerup( float bot, float top, int nval, float *val )
589 {
590 int ii ; double fac ;
591
592 if( nval < 1 || val == NULL ) return ;
593 if( nval == 1 ){ val[0] = sqrtf(bot*top); return; }
594 if( nval == 2 ){ val[0] = bot; val[1] = top; return; }
595 fac = pow( fabs(top/bot) , 1.0/(nval-1.0) ) ;
596 #if 0
597 ININFO_message("fillerup: bot=%g top=%g nval=%d fac=%g",bot,top,nval,fac) ;
598 #endif
599 val[0] = bot ;
600 for( ii=1 ; ii < nval-1 ; ii++ ) val[ii] = bot * pow(fac,(double)ii) ;
601 val[nval-1] = top ; return ;
602 }
603
604 /*-------------------------------------------------------------------------*/
605
ellish(float a,float b,float c,float d,float e,float f)606 static float ellish( float a , float b , float c ,
607 float d , float e , float f )
608 {
609 float dd , ee , ss=0.0f ;
610 dd = SQR(b-a)+SQR(e-d) ;
611 ee = SQR(c-b)+SQR(f-e) ;
612 if( dd > 0.0f && ee > 0.0f )
613 ss = fabsf((b-a)*(f-e)-(c-b)*(e-d)) / sqrtf(dd*ee) ;
614 return ss ;
615 }
616
617 /*-------------------------------------------------------------------------*/
618 /* Like THD_deconvolve, but choose the penalty factor automatically.
619 Will be significantly slower as it searches through penfac space.
620 *//*-----------------------------------------------------------------------*/
621
THD_deconvolve_autopen(int npt,float * far,int minlag,int maxlag,float * kern,int nbase,float * base[],int meth,float * ccon,int dcon,int pencode,float * penfac_used)622 static floatvec * THD_deconvolve_autopen( int npt , float *far ,
623 int minlag , int maxlag , float *kern,
624 int nbase , float *base[],
625 int meth , float *ccon , int dcon ,
626 int pencode, float *penfac_used )
627 {
628 floatvec **fvv=NULL , *fv=NULL , *gv=NULL ;
629 float pfac[NPFAC] , pres[NPFAC] , psiz[NPFAC] ;
630 float pbot,ptop , ppk , val , rpk,spk ;
631 int ii , ipk ;
632
633 ENTRY("THD_deconvolve_autopen") ;
634
635 /*--- solve many problems, using a crude mesh in pfac ---*/
636
637 fillerup( -0.01f , -10.0f , NPFAC , pfac ) ;
638 memset( (void *)pres , 0 , sizeof(float)*NPFAC ) ;
639 memset( (void *)psiz , 0 , sizeof(float)*NPFAC ) ;
640
641 fvv = THD_deconvolve_multipen( npt , far , minlag , maxlag , kern ,
642 nbase , base , meth , ccon , dcon ,
643 pencode , NPFAC , pfac , pres,psiz ) ;
644
645 if( fvv == NULL ){
646 ERROR_message(
647 "THD_deconvolve_autopen failed to solve initial problems: voxel ID=%d",
648 AO_VALUE(voxid) ) ;
649 RETURN(NULL) ;
650 }
651
652 /* find the best combination of residual and solution size */
653
654 ipk = -1 ; ppk = 0.0f ;
655 #if 0
656 for( ii=0 ; ii < NPFAC ; ii++ ){
657 val = pres[ii]*psiz[ii] ;
658 if( val > ppk ){ ipk = ii; ppk = val; }
659 }
660 #else
661 rpk = spk = 0.0f ;
662 for( ii=0 ; ii < NPFAC ; ii++ ){
663 rpk = MAX(rpk,pres[ii]) ; spk = MAX(rpk,psiz[ii]) ;
664 }
665 if( rpk > 0.0f ) rpk = 1.0f / rpk ;
666 if( spk > 0.0f ) spk = 1.0f / spk ;
667 for( ii=0 ; ii < NPFAC ; ii++ ){
668 pres[ii] *= rpk ; psiz[ii] *= spk ;
669 }
670 for( ii=1 ; ii < NPFAC-1 ; ii++ ){
671 val = ellish( pres[ii-1],pres[ii],pres[ii+1] ,
672 psiz[ii-1],psiz[ii],psiz[ii+1] ) ;
673 #if 1
674 if( AFNI_yesenv("AFNI_TFITTER_VERBOSE") ) ININFO_message("ellish[%d] = %g",ii,val) ;
675 #endif
676 if( val > ppk && psiz[ii] > 1.e-5 ){ ipk = ii; ppk = val; }
677 }
678 if( ipk > 0 && psiz[ipk+1] > 1.e-5 ) ipk++ ;
679 #endif
680
681 if( ipk < 0 || ppk == 0.0f ){ /* all fits failed?! */
682 for( ii=0 ; ii < NPFAC ; ii++ ) KILL_floatvec(fvv[ii]) ;
683 free((void *)fvv) ;
684 ERROR_message(
685 "THD_deconvolve_autopen fails to find initial good fit: voxel ID=%d",
686 AO_VALUE(voxid));
687 RETURN(NULL) ;
688 }
689
690 fv = fvv[ipk] ; if( do_fitv ){ COPY_floatvec(gv,AO_VALUE(ggfitvv)[ipk]); }
691 for( ii=0 ; ii < NPFAC ; ii++ ) if( ii != ipk ) KILL_floatvec(fvv[ii]) ;
692 free((void *)fvv) ;
693 if( penfac_used != NULL ) *penfac_used = pfac[ipk] ;
694 #if 1
695 if( AFNI_yesenv("AFNI_TFITTER_VERBOSE") )
696 ININFO_message("optimal penfac_used#%d = %g",ipk,pfac[ipk]) ;
697 #endif
698
699 /*--- refinement step in pfac: scan around best one found above ---*/
700
701 pbot = (ipk > 0) ? 0.8f*pfac[ipk-1] : 0.5f*pfac[0] ;
702 ptop = (ipk < NPFAC-1) ? 1.0f*pfac[ipk+1] : 2.0f*pfac[NPFAC-1] ;
703 fillerup( pbot , ptop , NPFAC , pfac ) ;
704 memset( (void *)pres , 0 , sizeof(float)*NPFAC ) ;
705 memset( (void *)psiz , 0 , sizeof(float)*NPFAC ) ;
706
707 fvv = THD_deconvolve_multipen( npt , far , minlag , maxlag , kern ,
708 nbase , base , meth , ccon , dcon ,
709 pencode , NPFAC , pfac , pres,psiz ) ;
710
711 if( fvv == NULL ){ /* failed ==> return what we found earlier */
712 ERROR_message("THD_deconvolve_autopen semi-failed: voxel ID=%d",AO_VALUE(voxid)) ;
713 if( do_fitv ){ KILL_floatvec(AO_VALUE(gfitv)); AO_VALUE(gfitv) = gv; }
714 RETURN(fv) ;
715 }
716
717 ipk = -1 ; ppk = 0.0f ;
718 #if 0
719 for( ii=0 ; ii < NPFAC ; ii++ ){
720 val = pres[ii]*psiz[ii] ;
721 if( val > ppk ){ ipk = ii; ppk = val; }
722 }
723 #else
724 rpk = spk = 0.0f ;
725 for( ii=0 ; ii < NPFAC ; ii++ ){
726 rpk = MAX(rpk,pres[ii]) ; spk = MAX(rpk,psiz[ii]) ;
727 }
728 if( rpk > 0.0f ) rpk = 1.0f / rpk ;
729 if( spk > 0.0f ) spk = 1.0f / spk ;
730 for( ii=0 ; ii < NPFAC ; ii++ ){
731 pres[ii] *= rpk ; psiz[ii] *= spk ;
732 }
733 for( ii=1 ; ii < NPFAC-1 ; ii++ ){
734 val = ellish( pres[ii-1],pres[ii],pres[ii+1] ,
735 psiz[ii-1],psiz[ii],psiz[ii+1] ) ;
736 #if 1
737 if( AFNI_yesenv("AFNI_TFITTER_VERBOSE") ) ININFO_message("ellish[%d] = %g",ii,val) ;
738 #endif
739 if( val > ppk && psiz[ii] > 1.e-5 ){ ipk = ii; ppk = val; }
740 }
741 if( ipk > 0 && psiz[ipk+1] > 1.e-5 ) ipk++ ;
742 #endif
743
744 if( ipk < 0 || ppk == 0.0f ){ /* all failed? use old result */
745 for( ii=0 ; ii < NPFAC ; ii++ ) KILL_floatvec(fvv[ii]) ;
746 free((void *)fvv) ;
747 ERROR_message("THD_deconvolve_autopen semi-fails: voxel ID=%d",AO_VALUE(voxid)) ;
748 if( do_fitv ){ KILL_floatvec(AO_VALUE(gfitv)); AO_VALUE(gfitv) = gv; }
749 RETURN(fv) ;
750 }
751
752 KILL_floatvec(fv) ; fv = fvv[ipk] ;
753 if( do_fitv ){ KILL_floatvec(gv); COPY_floatvec(gv,AO_VALUE(ggfitvv)[ipk]); }
754 for( ii=0 ; ii < NPFAC ; ii++ ) if( ii != ipk ) KILL_floatvec(fvv[ii]) ;
755 free((void *)fvv) ;
756 if( penfac_used != NULL ) *penfac_used = pfac[ipk] ;
757 #if 1
758 if( AFNI_yesenv("AFNI_TFITTER_VERBOSE") )
759 ININFO_message("Optimal penfac_used#%d = %g",ipk,pfac[ipk]) ;
760 #endif
761
762 if( do_fitv ){ KILL_floatvec(AO_VALUE(gfitv)); AO_VALUE(gfitv) = gv; }
763 RETURN(fv) ;
764 }
765