00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008 #include "betafit.c"
00009
00010 #undef SINGLET
00011
00012 typedef struct {
00013 int bot , top ;
00014 float pval ;
00015 } spanfit ;
00016
00017 typedef struct {
00018 int ndim ;
00019 float * cmat , * cfac , * mvec ;
00020 } covmat ;
00021
00022 #define IFREE(x) do{if((x)!=NULL)free(x);}while(0)
00023
00024 #define FREE_COVMAT(cc) \
00025 do{ if(cc != NULL){ \
00026 IFREE(cc->cmat); IFREE(cc->cfac); \
00027 IFREE(cc->mvec); free(cc); } } while(0)
00028
00029 #define CM(i,j) cmat[(i)+(j)*ndim]
00030 #define CH(i,j) cfac[(i)+(j)*ndim]
00031
00032
00033
00034 void forward_solve_inplace( covmat * cv , float * vec )
00035 {
00036 register int ndim=cv->ndim , ii,jj ;
00037 register float * cfac=cv->cfac , sum ;
00038
00039 for( ii=0 ; ii < ndim ; ii++ ){
00040 sum = vec[ii] ;
00041 for( jj=0 ; jj < ii ; jj++ ) sum -= CH(ii,jj) * vec[jj] ;
00042 vec[ii] = sum / CH(ii,ii) ;
00043 }
00044 return ;
00045 }
00046
00047
00048
00049 void backward_solve_inplace( covmat * cv , float * vec )
00050 {
00051 register int ndim=cv->ndim , ii,jj ;
00052 register float * cfac=cv->cfac , sum ;
00053
00054 for( ii=ndim-1 ; ii >= 0 ; ii-- ){
00055 sum = vec[ii] ;
00056 for( jj=ii+1 ; jj < ndim ; jj++ ) sum -= CH(jj,ii) * vec[jj] ;
00057 vec[ii] = sum / CH(ii,ii) ;
00058 }
00059 return ;
00060 }
00061
00062
00063
00064 void compute_choleski( covmat * cv )
00065 {
00066 register int ndim=cv->ndim , ii,jj,kk ;
00067 register float * cmat=cv->cmat , * cfac , sum ;
00068
00069 if( ndim < 1 || cmat == NULL ) return ;
00070
00071 if( cv->cfac == NULL )
00072 cv->cfac = (float *) malloc(sizeof(float)*ndim*ndim) ;
00073
00074 cfac = cv->cfac ;
00075
00076 if( ndim == 1 ){
00077 if( cmat[0] <= 0.0 ){ free(cv->cfac); cv->cfac = NULL; return; }
00078 cfac[0] = sqrt(cmat[0]) ;
00079 return ;
00080 }
00081
00082 for( ii=0 ; ii < ndim ; ii++ ){
00083 for( jj=0 ; jj < ii ; jj++ ){
00084 sum = CM(ii,jj) ;
00085 for( kk=0 ; kk < jj ; kk++ ) sum -= CH(ii,kk) * CH(jj,kk) ;
00086 CH(ii,jj) = sum / CH(jj,jj) ;
00087 }
00088 sum = CM(ii,ii) ;
00089 for( kk=0 ; kk < ii ; kk++ ) sum -= CH(ii,kk) * CH(ii,kk) ;
00090 if( sum <= 0.0 ){ free(cv->cfac); cv->cfac = NULL; return; }
00091 CH(ii,ii) = sqrt(sum) ;
00092 for( jj=ii+1 ; jj < ndim ; jj++ ) CH(ii,jj) = 0.0 ;
00093 }
00094 return ;
00095 }
00096
00097
00098
00099 #define CCUT 3.5
00100 #define EPS 1.e-4
00101
00102 covmat * robust_covar( int ndim , int nvec , float ** vec )
00103 {
00104 covmat * cv ;
00105 float *nmat, *cmat , fnvec,fndim,cnorm,csum , *tv , *vv , *mv , *wv ;
00106 int ii , jj , kk , nite ;
00107 float bcut , cwt ;
00108
00109 #if 0
00110 fprintf(stderr,"Enter robust_covar: ndim=%d nvec=%d\n",ndim,nvec) ;
00111 #endif
00112
00113 if( ndim < 2 || nvec < ndim || vec == NULL ) return NULL ;
00114
00115 cv = (covmat *) malloc(sizeof(covmat)) ;
00116 cv->ndim = ndim ;
00117 cv->cmat = NULL ;
00118 cv->cfac = NULL ;
00119 cv->mvec = NULL ;
00120
00121 nmat = (float *) malloc(sizeof(float)*ndim*ndim) ;
00122 tv = (float *) malloc(sizeof(float)*ndim) ;
00123 mv = (float *) malloc(sizeof(float)*ndim) ;
00124 wv = (float *) malloc(sizeof(float)*nvec) ;
00125
00126 fnvec = 1.0/nvec ; fndim = 1.0/ndim ;
00127 bcut = 1.0 + CCUT*sqrt(fndim) ;
00128
00129
00130
00131 for( jj=0 ; jj < ndim ; jj++ ) mv[jj] = 0.0 ;
00132
00133 for( kk=0 ; kk < nvec ; kk++ ){
00134 vv = vec[kk] ;
00135 for( jj=0 ; jj < ndim ; jj++ ) mv[jj] += vv[jj] ;
00136 }
00137 for( jj=0 ; jj < ndim ; jj++ ) mv[jj] *= fnvec ;
00138
00139 for( jj=0 ; jj < ndim ; jj++ )
00140 for( ii=0 ; ii < ndim ; ii++ ) nmat[ii+jj*ndim] = 0.0 ;
00141
00142 for( kk=0 ; kk < nvec ; kk++ ){
00143 vv = vec[kk] ;
00144 for( jj=0 ; jj < ndim ; jj++ ){
00145 for( ii=0 ; ii <= jj ; ii++ )
00146 nmat[ii+jj*ndim] += (vv[ii]-mv[ii])*(vv[jj]-mv[jj]) ;
00147 }
00148 }
00149 for( jj=0 ; jj < ndim ; jj++ ){
00150 for( ii=0 ; ii < jj ; ii++ )
00151 nmat[jj+ii*ndim] = (nmat[ii+jj*ndim] *= fnvec) ;
00152 nmat[jj+jj*ndim] *= fnvec ;
00153 }
00154
00155
00156
00157 nite = 0 ;
00158
00159 while(1){
00160
00161 nite++ ;
00162
00163 #if 0
00164 fprintf(stderr,"\niteration %2d:\n",nite) ;
00165 #endif
00166
00167 cmat = cv->cmat = nmat ;
00168 cv->mvec = mv ;
00169 compute_choleski(cv) ;
00170
00171 if( cv->cfac == NULL ){
00172 free(cv->cmat); free(cv->mvec); free(cv); free(tv); free(wv);
00173 return NULL ;
00174 }
00175
00176 nmat = (float *) malloc(sizeof(float)*ndim*ndim) ;
00177 mv = (float *) malloc(sizeof(float)*ndim) ;
00178
00179 for( jj=0 ; jj < ndim ; jj++ ){
00180 mv[jj] = 0.0 ;
00181 for( ii=0 ; ii < ndim ; ii++ ) nmat[ii+jj*ndim] = 0.0 ;
00182 }
00183
00184
00185
00186 csum = 0.0 ;
00187 for( kk=0 ; kk < nvec ; kk++ ){
00188 vv = vec[kk] ;
00189
00190
00191
00192
00193 for( jj=0 ; jj < ndim ; jj++ ) tv[jj] = vv[jj] - cv->mvec[jj] ;
00194 forward_solve_inplace(cv,tv) ;
00195
00196
00197
00198 cnorm = 0.0 ; for( ii=0 ; ii < ndim ; ii++ ) cnorm += tv[ii]*tv[ii] ;
00199 cnorm = cnorm*fndim ;
00200 cnorm = (cnorm <= bcut) ? 1.0 : bcut/cnorm ;
00201 wv[kk] = cnorm ; csum += cnorm ;
00202
00203
00204
00205 for( jj=0 ; jj < ndim ; jj++ ) mv[jj] += cnorm*vv[jj] ;
00206 }
00207 csum = 1.0 / csum ; cwt = nvec*csum ;
00208 for( jj=0 ; jj < ndim ; jj++ ) mv[jj] *= csum ;
00209
00210
00211
00212 for( kk=0 ; kk < nvec ; kk++ ){
00213 vv = vec[kk] ; cnorm = wv[kk] ;
00214 for( jj=0 ; jj < ndim ; jj++ ){
00215 for( ii=0 ; ii <= jj ; ii++ )
00216 nmat[ii+jj*ndim] +=
00217 cnorm*(vv[ii]-cv->mvec[ii])*(vv[jj]-cv->mvec[jj]) ;
00218 }
00219 }
00220 #define DDD csum
00221 for( jj=0 ; jj < ndim ; jj++ ){
00222 for( ii=0 ; ii < jj ; ii++ )
00223 nmat[jj+ii*ndim] = (nmat[ii+jj*ndim] *= DDD) ;
00224 nmat[jj+jj*ndim] *= DDD ;
00225 }
00226
00227
00228
00229 cnorm = csum = 0.0 ;
00230 for( jj=0 ; jj < ndim ; jj++ ){
00231 for( ii=0 ; ii <= jj ; ii++ ){
00232 cnorm += fabs( nmat[ii+jj*ndim] - cmat[ii+jj*ndim] ) ;
00233 csum += fabs( nmat[ii+jj*ndim] ) ;
00234 }
00235 }
00236
00237 #if 0
00238 fprintf(stderr," |dif|=%12.4g |mat|=%12.4g cwt=%12.4g\n",cnorm,csum,cwt) ;
00239 fprintf(stderr," matrix:\n") ;
00240 for( ii=0 ; ii < ndim ; ii++ ){
00241 fprintf(stderr," Row%2d: %12.4g ",ii,mv[ii]) ;
00242 for( jj=0 ; jj < ndim ; jj++ )
00243 fprintf(stderr," %12.4g",
00244 (jj<=ii) ? nmat[ii+jj*ndim] :
00245 nmat[ii+jj*ndim]/sqrt(nmat[ii+ii*ndim]*nmat[jj+jj*ndim]) );
00246 fprintf(stderr,"\n") ;
00247 }
00248 #endif
00249
00250 free(cv->cmat) ; free(cv->mvec) ;
00251 if( cnorm <= EPS*csum || nite > 3*ndim ){
00252 cv->cmat = nmat; cv->mvec = mv; break;
00253 }
00254 }
00255
00256 free(wv) ; free(tv) ; compute_choleski(cv) ; return cv ;
00257 }
00258
00259
00260
00261 float evaluate_span( int ndim, int nvec,
00262 int bot , int top ,
00263 int nbasis , float ** basis ,
00264 float * cvec , float ** bvec )
00265 {
00266 int ii,kk,jj , npt=top-bot+1 , nbd ;
00267 float ** svec , **ee,*xx,*cc,s,t,xd,tinv , bd,dx ;
00268 covmat * cv , * ce ;
00269
00270 if( nbasis >= npt ) return 0.0 ;
00271
00272
00273
00274 svec = (float **) malloc(sizeof(float *)*nvec) ;
00275 for( kk=0 ; kk < nvec ; kk++ ) svec[kk] = bvec[kk] + bot ;
00276
00277
00278
00279 cv = robust_covar( npt , nvec , svec ) ;
00280 free(svec) ;
00281 if( cv == NULL ) return 0.0 ;
00282
00283
00284
00285 #if 0
00286 fprintf(stderr,"Normalizing cvec, basis\n") ;
00287 #endif
00288
00289 ee = (float **) malloc(sizeof(float *)*nbasis) ;
00290 for( kk=0 ; kk < nbasis ; kk++ )
00291 ee[kk] = (float *) malloc(sizeof(float)*npt) ;
00292 xx = (float *) malloc(sizeof(float)*npt) ;
00293 for( ii=0 ; ii < npt ; ii++ ){
00294 xx[ii] = cvec[ii+bot] ;
00295 for( kk=0 ; kk < nbasis ; kk++ ) ee[kk][ii] = basis[kk][ii+bot] ;
00296 }
00297 forward_solve_inplace( cv , xx ) ;
00298 for( kk=0 ; kk < nbasis ; kk++ )
00299 forward_solve_inplace( cv , ee[kk] ) ;
00300
00301 #if 0
00302 fprintf(stderr,"Computing [ee]**T [ee]\n") ;
00303 #endif
00304
00305
00306
00307
00308 ce = (covmat *) malloc(sizeof(covmat)) ;
00309 ce->ndim = nbasis ;
00310 ce->cmat = (float *) malloc(sizeof(float)*nbasis*nbasis) ;
00311 ce->cfac = NULL ;
00312 ce->mvec = NULL ;
00313 for( kk=0 ; kk < nbasis ; kk++ ){
00314 for( jj=0 ; jj <= kk ; jj++ ){
00315 s = 0.0 ;
00316 for( ii=0 ; ii < npt ; ii++ ) s += ee[kk][ii] * ee[jj][ii] ;
00317 ce->cmat[jj+kk*nbasis] = s ;
00318 if( jj < kk ) ce->cmat[kk+jj*nbasis] = s ;
00319 }
00320 }
00321
00322
00323
00324 compute_choleski(ce) ;
00325 if( ce->cfac == NULL )fprintf(stderr,"choleski failed!\n") ;
00326 cc = (float *) malloc(sizeof(float)*nbasis) ;
00327
00328 for( kk=0 ; kk < nbasis ; kk++ ){
00329 s = 0.0 ;
00330 for( ii=0 ; ii < npt ; ii++ ) s += ee[kk][ii] * xx[ii] ;
00331 cc[kk] = s ;
00332 }
00333 forward_solve_inplace( ce , cc ) ;
00334 backward_solve_inplace( ce , cc ) ;
00335 for( ii=0 ; ii < npt ; ii++ ){
00336 s = xx[ii] ;
00337 for( kk=0 ; kk < nbasis ; kk++ ) s -= ee[kk][ii] * cc[kk] ;
00338 xx[ii] = s ;
00339 }
00340
00341
00342
00343
00344
00345 nbd = 0 ;
00346 for( kk=0 ; kk < nvec ; kk++ ){
00347 memcpy( ee[0] , bvec[kk]+bot , sizeof(float)*npt ) ;
00348 forward_solve_inplace( cv , ee[0] ) ;
00349 bd = 0.0 ;
00350 for( ii=0 ; ii < npt ; ii++ ) bd += xx[ii] * ee[0][ii] ;
00351 if( bd <= 0.0 ) nbd++ ;
00352 #if 0
00353 fprintf(stderr," %12.4g",bd);
00354 #endif
00355 }
00356 s = (float)nbd / (float)nvec ;
00357 #if 0
00358 fprintf(stderr," => nbd=%d\n",nbd) ;
00359 #endif
00360
00361 for( kk=0 ; kk < nbasis ; kk++ ) free(ee[kk]) ;
00362 free(ee) ; free(cc) ; free(xx); FREE_COVMAT(cv) ; FREE_COVMAT(ce) ;
00363 return s ;
00364 }
00365
00366
00367
00368 spanfit find_best_span( int ndim , int nvec , int minspan ,
00369 int nbasis , float ** basis ,
00370 float * cvec , float ** bvec )
00371 {
00372 spanfit result = {0,0,0.0} ;
00373 int ii,kk , bot,top , bot_best,top_best ;
00374 float val , val_best ;
00375
00376 if( minspan < 3 || ndim < minspan || nvec < 100 ) return result ;
00377 if( cvec == NULL || bvec == NULL ) return result ;
00378 if( nbasis < 1 || basis == NULL ) return result ;
00379
00380 val_best = -1.0 ;
00381 for( bot=0 ; bot < ndim+1-minspan ; bot++ ){
00382
00383 for( top=0 ; top < bot+minspan-1 ; top++ ) printf(" 0") ;
00384
00385 for( top=bot+minspan-1 ; top < ndim ; top++ ){
00386 val = evaluate_span( ndim,nvec , bot,top , nbasis,basis , cvec,bvec ) ;
00387
00388 printf(" %g",val) ;
00389
00390 if( val > val_best ){
00391 val_best = val ; bot_best = bot ; top_best = top ;
00392 }
00393 #if 1
00394 if( val >= 0.10 ) fprintf(stderr,"bot=%2d top=%2d: %.4f\n",bot,top,val) ;
00395 #endif
00396 }
00397
00398 printf("\n") ;
00399 }
00400
00401 result.bot = bot_best; result.top = top_best; result.pval = val_best;
00402 return result ;
00403 }
00404
00405
00406
00407 static int nran=1000 ;
00408 static float abot= 0.5 , atop= 4.0 ;
00409 static float bbot=10.0 , btop=200.0 ;
00410 static float pbot=50.0 , ptop= 80.0 ;
00411 static double pthr=1.e-4 ;
00412 static int sqr=0 ;
00413
00414 #define OUT_THR 1
00415 #define OUT_BBB 2
00416 #define OUT_AAA 3
00417 #define OUT_AB 4
00418
00419 static int outmode = OUT_THR ;
00420 static int nbasis = 1 ;
00421
00422 #define FVECTOR_DIM 2
00423 typedef struct {
00424 float v[FVECTOR_DIM] ;
00425 } fvector ;
00426
00427
00428
00429 fvector process_sample( float pcut , BFIT_data * bfd )
00430 {
00431 BFIT_result * bfr ;
00432 double xth ;
00433 fvector result ;
00434
00435 static double aold,bold ;
00436 static BFIT_data * bfdold=NULL ;
00437
00438 #if 1
00439 if( bfd == bfdold ){
00440 beta_init( aold , bold ) ;
00441 nran = 400 ;
00442 abot = aold*0.5 ; atop = aold*2.0 ; if( abot <= 0.1 ) abot = 0.101 ;
00443 bbot = bold*0.5 ; btop = bold*2.0 ; if( bbot <= 9.9 ) bbot = 9.999 ;
00444 } else {
00445 beta_init( 0.0 , 0.0 ) ;
00446 bfdold = bfd ;
00447 nran = 1000 ;
00448 abot = 0.5 ; atop = 4.0 ;
00449 bbot = 10.0 ; btop =200.0 ;
00450 }
00451 #endif
00452
00453 bfr = BFIT_compute( bfd , pcut , abot,atop , bbot,btop , nran,0 ) ;
00454 if( bfr == NULL ){
00455 fprintf(stderr,"*** Can't compute betafit at pcut=%f\n",pcut) ;
00456 exit(1) ;
00457 }
00458
00459 aold = bfr->a ; bold = bfr->b ;
00460
00461 switch( outmode ){
00462
00463 default:
00464 case OUT_THR:
00465 xth = beta_p2t( pthr , bfr->a,bfr->b ) ;
00466 if( sqr ) xth = sqrt(xth) ;
00467 result.v[0] = xth ;
00468 break ;
00469
00470 case OUT_BBB: result.v[0] = bold ; break ;
00471 case OUT_AAA: result.v[0] = aold ; break ;
00472
00473 case OUT_AB: result.v[0] = aold ; result.v[1] = bold ; break ;
00474 }
00475
00476 BFIT_free_result(bfr) ;
00477 return result ;
00478 }
00479
00480
00481
00482 int main( int argc , char * argv[] )
00483 {
00484 BFIT_data * bfd , * nfd ;
00485 float * bf_tvec , ** boot_tvec ;
00486 int ndim , nvec ;
00487
00488 int nvals,ival , nvox , nbin , miv ;
00489 float pcut , eps,eps1 ;
00490 float *bval , *cval ;
00491 double aa,bb,xc,xth ;
00492 fvector bfit ;
00493 float ** basis ;
00494
00495 int mcount,mgood , ii , jj , kk , ibot,itop , qq ;
00496
00497 int narg=1 ;
00498 int nboot=0 ;
00499 double aboot,bboot,tboot , pthr=1.e-4 ;
00500 float asig , bsig , tsig , abcor ;
00501
00502 THD_3dim_dataset * input_dset , * mask_dset=NULL ;
00503 float mask_bot=666.0 , mask_top=-666.0 ;
00504 byte * mmm=NULL ;
00505
00506 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00507 fprintf(stderr,"Usage: 3dbetafit3 [options] dataset\n"
00508 "Fits a beta distribution to the values in a brick.\n"
00509 "\n"
00510 "Options:\n"
00511 " -arange abot atop = Sets the search range for parameter\n"
00512 " 'a' to abot..atop.\n"
00513 " [default is 0.5 .. 4.0]\n"
00514 "\n"
00515 " -brange bbot btop = Sets the search range for parameter\n"
00516 " 'b' to bbot..btop\n"
00517 " [default is 10 .. 200]\n"
00518 "\n"
00519 " -prange pbot ptop = Will evaluate for percent cutoffs\n"
00520 " from pbot to ptop (steps of 1%%)\n"
00521 " [default is 50 .. 80]\n"
00522 "\n"
00523 " -bootstrap N = Does N bootstrap evaluations\n"
00524 "\n"
00525 " -mask mset = A mask dataset to indicate which\n"
00526 " voxels are to be used\n"
00527 " -mrange b t = Use only mask values in range from\n"
00528 " 'b' to 't' (inclusive)\n"
00529 "\n"
00530 " -sqr = Flag to square the data from the dataset\n"
00531 " -pthr p = Sets p-value of cutoff for threshold evaluation\n"
00532 " [default = 1.e-4]\n"
00533 " -bout = Use 'b' for the output, instead of thr\n"
00534 " -aout = Use 'a' for the output, instead of thr\n"
00535 " -about = Use 'a' AND 'b' for the ouput.\n"
00536 ) ;
00537 exit(0) ;
00538 }
00539
00540
00541
00542 while( narg < argc && argv[narg][0] == '-' ){
00543
00544 if( strcmp(argv[narg],"-aout") == 0 ){
00545 outmode = OUT_AAA ; nbasis = 1 ; narg++ ; continue ;
00546 }
00547
00548 if( strcmp(argv[narg],"-bout") == 0 ){
00549 outmode = OUT_BBB ; nbasis = 1 ; narg++ ; continue ;
00550 }
00551
00552 if( strcmp(argv[narg],"-about") == 0 ){
00553 outmode = OUT_AB ; nbasis = 2 ; narg++ ; continue ;
00554 }
00555
00556 if( strcmp(argv[narg],"-pthr") == 0 ){
00557 pthr = strtod(argv[++narg],NULL) ;
00558 if( pthr <= 0.0 || pthr >= 1.0 ){
00559 fprintf(stderr,"*** Illegal value after -pthr!\n");exit(1);
00560 }
00561 narg++ ; continue;
00562 }
00563
00564 if( strcmp(argv[narg],"-sqr") == 0 ){
00565 sqr = 1 ; narg++ ; continue;
00566 }
00567
00568 if( strcmp(argv[narg],"-arange") == 0 ){
00569 abot = strtod(argv[++narg],NULL) ;
00570 atop = strtod(argv[++narg],NULL) ;
00571 if( abot < 0.1 || abot > atop ){
00572 fprintf(stderr,"*** Illegal value after -arange!\n");exit(1);
00573 }
00574 narg++ ; continue;
00575 }
00576
00577 if( strcmp(argv[narg],"-brange") == 0 ){
00578 bbot = strtod(argv[++narg],NULL) ;
00579 btop = strtod(argv[++narg],NULL) ;
00580 if( bbot < 0.1 || bbot > btop ){
00581 fprintf(stderr,"*** Illegal value after -brange!\n");exit(1);
00582 }
00583 narg++ ; continue;
00584 }
00585
00586 if( strcmp(argv[narg],"-prange") == 0 ){
00587 pbot = (int) strtod(argv[++narg],NULL) ;
00588 ptop = (int) strtod(argv[++narg],NULL) ;
00589 if( pbot < 30.0 || pbot > ptop || ptop > 99.0 ){
00590 fprintf(stderr,"*** Illegal value after -prange!\n");exit(1);
00591 }
00592 narg++ ; continue;
00593 }
00594
00595 if( strcmp(argv[narg],"-bootstrap") == 0 ){
00596 nboot = (int) strtod(argv[++narg],NULL) ;
00597 if( nboot < 100 ){
00598 fprintf(stderr,"*** Illegal value after -bootstrap!\n");exit(1);
00599 }
00600 narg++ ; continue;
00601 }
00602
00603 if( strncmp(argv[narg],"-mask",5) == 0 ){
00604 if( mask_dset != NULL ){
00605 fprintf(stderr,"*** Cannot have two -mask options!\n") ; exit(1) ;
00606 }
00607 if( narg+1 >= argc ){
00608 fprintf(stderr,"*** -mask option requires a following argument!\n");
00609 exit(1) ;
00610 }
00611 mask_dset = THD_open_dataset( argv[++narg] ) ;
00612 if( mask_dset == NULL ){
00613 fprintf(stderr,"*** Cannot open mask dataset!\n") ; exit(1) ;
00614 }
00615 if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex ){
00616 fprintf(stderr,"*** Cannot deal with complex-valued mask dataset!\n");
00617 exit(1) ;
00618 }
00619 narg++ ; continue ;
00620 }
00621
00622 if( strncmp(argv[narg],"-mrange",5) == 0 ){
00623 if( narg+2 >= argc ){
00624 fprintf(stderr,"*** -mrange option requires 2 following arguments!\n") ;
00625 exit(1) ;
00626 }
00627 mask_bot = strtod( argv[++narg] , NULL ) ;
00628 mask_top = strtod( argv[++narg] , NULL ) ;
00629 if( mask_top < mask_top ){
00630 fprintf(stderr,"*** -mrange inputs are illegal!\n") ; exit(1) ;
00631 }
00632 narg++ ; continue ;
00633 }
00634
00635 fprintf(stderr,"*** Illegal option: %s\n",argv[narg]) ; exit(1) ;
00636 }
00637
00638 if( nboot < 100 ){
00639 fprintf(stderr,"*** Must use -bootstrap 'option'!\n"); exit(1);
00640 }
00641
00642 if( narg >= argc ){
00643 fprintf(stderr,"*** No dataset argument on command line!?\n"); exit(1);
00644 }
00645
00646 input_dset = THD_open_dataset( argv[narg] ) ;
00647 if( input_dset == NULL ){
00648 fprintf(stderr,"*** Can't open dataset %s\n",argv[narg]); exit(1);
00649 }
00650
00651 nvox = DSET_NVOX(input_dset) ;
00652
00653
00654
00655 DSET_load(input_dset) ;
00656 if( !DSET_LOADED(input_dset) ){
00657 fprintf(stderr,"*** Couldn't load dataset brick!\n");exit(1);
00658 }
00659
00660 if( DSET_BRICK_STATCODE(input_dset,0) == FUNC_COR_TYPE ) sqr = 1 ;
00661
00662 bfd = BFIT_prepare_dataset( input_dset , 0 , sqr ,
00663 mask_dset , 0 , mask_bot , mask_top ) ;
00664
00665 if( bfd == NULL ){
00666 fprintf(stderr,"*** Couldn't prepare data from input dataset!\n");
00667 exit(1) ;
00668 }
00669
00670 DSET_delete(mask_dset) ; DSET_delete(input_dset) ;
00671
00672
00673
00674 fprintf(stderr,"Computing bootstrap") ;
00675
00676 ndim = (ptop - pbot + 1.0)*nbasis ;
00677
00678 basis = (float **) malloc(sizeof(float *)*nbasis) ;
00679 for( qq=0 ; qq < nbasis ; qq++ )
00680 basis[qq] = (float *) malloc(sizeof(float)*ndim) ;
00681 switch(nbasis){
00682 case 1:
00683 for( jj=0 ; jj < ndim ; jj++ ) basis[0][jj] = 1.0 ;
00684 break ;
00685
00686 case 2:
00687 for( jj=0 ; jj < ndim ; jj++ ){
00688 basis[0][jj] = ((jj%2) == 0) ;
00689 basis[1][jj] = ((jj%2) == 1) ;
00690 }
00691 break ;
00692 }
00693
00694 bf_tvec = (float *) malloc(sizeof(float)*ndim) ;
00695 for( pcut=pbot ; pcut <= ptop ; pcut += 1.0 ){
00696 bfit = process_sample( pcut , bfd ) ;
00697 qq = (int)(pcut-pbot) ;
00698 switch(nbasis){
00699 case 1: bf_tvec[qq] = bfit.v[0] ; break ;
00700 case 2: bf_tvec[2*qq] = bfit.v[0] ; bf_tvec[2*qq+1] = bfit.v[1] ; break ;
00701 }
00702 }
00703
00704 nvec = nboot ;
00705 boot_tvec = (float **) malloc(sizeof(float *)*nvec) ;
00706 for( jj=0 ; jj < nboot ; jj++ ){
00707 boot_tvec[jj] = (float *) malloc(sizeof(float)*ndim) ;
00708 nfd = BFIT_bootstrap_sample( bfd ) ;
00709 for( pcut=pbot ; pcut <= ptop ; pcut += 1.0 ){
00710 bfit = process_sample( pcut , nfd ) ;
00711 qq = (int)(pcut-pbot) ;
00712 switch(nbasis){
00713 case 1: boot_tvec[jj][qq] = bfit.v[0] ; break ;
00714 case 2: boot_tvec[jj][2*qq] = bfit.v[0] ; boot_tvec[jj][2*qq+1] = bfit.v[1] ; break ;
00715 }
00716 }
00717 BFIT_free_data(nfd) ;
00718 if( jj%10 == 0 ) fprintf(stderr,".") ;
00719 }
00720 fprintf(stderr,"\n") ;
00721 BFIT_free_data(bfd) ;
00722
00723 #ifdef SINGLET
00724 while(1){
00725 fprintf(stderr,"Enter ibot itop [ndim=%d]: ",ndim) ;
00726 ibot = itop = 0 ;
00727 fscanf(stdin,"%d%d",&ibot,&itop) ;
00728 if( itop < 0 || itop-ibot+1 < 3 || itop >= ndim ) continue ;
00729
00730 eps = evaluate_span( ndim,nvec , ibot,itop , nbasis,basis , bf_tvec,boot_tvec ) ;
00731 fprintf(stderr,"Evaluate = %f\n\n",eps) ;
00732 }
00733 #else
00734 { spanfit sf = find_best_span( ndim,nvec , 10*nbasis , nbasis,basis , bf_tvec,boot_tvec ) ;
00735 float tbar = 0.0 ;
00736 for( ii=sf.bot ; ii <= sf.top ; ii++ ) tbar += bf_tvec[ii] ;
00737 tbar /= (sf.top-sf.bot+1.0) ;
00738 fprintf(stderr,"\nBEST bot=%2d top=%2d: %.4f %12.4g\n",sf.bot,sf.top,sf.pval,tbar) ;
00739 }
00740 #endif
00741
00742 #if 0
00743 { float xx,ss ;
00744 for( pcut=pbot ; pcut <= ptop ; pcut += 1.0 ){
00745 kk = (int)(pcut-pbot) ;
00746 xx = bf_tvec[kk] ;
00747 ss = 0.0 ;
00748 for( jj=0 ; jj < nboot ; jj++ ) ss += SQR((boot_tvec[jj][kk]-xx)) ;
00749 ss = sqrt(ss/nboot) ;
00750 printf("%.1f %12.4g %12.4g\n",pcut,xx,ss) ;
00751 }
00752 }
00753 #endif
00754
00755 exit(0) ;
00756 }