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