/*****************************************************************************
   Major portions of this software are copyrighted by the Medical College
   of Wisconsin, 1994-2000, and are released under the Gnu General Public
   License, Version 2.  See the file README.Copyright for details.
******************************************************************************/
   
#include "mrilib.h"
#include "betafit.c"

/*-----------------------------------------------------------------------*/

int main( int argc , char * argv[] )
{
   BFIT_data   * bfd , * nfd ;
   BFIT_result * bfr , * nfr ;

   int nvals,ival , nvox , nbin , miv , sqr=0 ;
   float pcut , eps,eps1 ;
   float *bval , *cval ;
   double aa,bb,xc,xth ;
   double chq,ccc,cdf ;
   int    ihqbot,ihqtop ;

   int mcount,mgood , ii , jj , ibot,itop ;

   int narg=1 , nran=1000 ;
   float abot= 0.5 , atop=  4.0 ;
   float bbot=10.0 , btop=200.0 ;
   float pbot=50.0 , ptop= 80.0 ;

   int nboot=0 ;
   double  aboot,bboot,tboot , pthr=1.e-4 ;
   float   asig , bsig , tsig , abcor ;

   THD_3dim_dataset * input_dset , * mask_dset=NULL ;
   float mask_bot=666.0 , mask_top=-666.0 ;
   byte * mmm=NULL ;

   if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
      printf("Usage: 3dbetafit [options] dataset\n"
             "Fits a beta distribution to the values in a brick.\n"
             "\n"
             "Options:\n"
             "  -arange abot atop = Sets the search range for parameter\n"
             "                        'a' to abot..atop.\n"
             "                        [default is 0.5 .. 4.0]\n"
             "\n"
             "  -brange bbot btop = Sets the search range for parameter\n"
             "                        'b' to bbot..btop\n"
             "                        [default is 10 .. 200]\n"
             "\n"
             "  -prange pbot ptop = Will evaluate for percent cutoffs\n"
             "                        from pbot to ptop (steps of 1%%)\n"
             "                        [default is 50 .. 80]\n"
             "\n"
             "  -bootstrap N      = Does N bootstrap evaluations to\n"
             "                        compute variance of 'a' and 'b'\n"
             "                        estimates [default is no bootstrap]\n"
             "\n"
             "  -mask mset  = A mask dataset to indicate which\n"
             "                 voxels are to be used\n"
             "  -mrange b t = Use only mask values in range from\n"
             "                 'b' to 't' (inclusive)\n"
             "\n"
             "  -sqr    = Flag to square the data from the dataset\n"
             "  -pthr p = Sets p-value of cutoff for threshold evaluation\n"
             "              [default = 1.e-4]\n"
         ) ;
         exit(0) ;
   }

   /* scan command-line args */

   while( narg < argc && argv[narg][0] == '-' ){

      if( strcmp(argv[narg],"-pthr") == 0 ){
         pthr = strtod(argv[++narg],NULL) ;
         if( pthr <= 0.0 || pthr >= 1.0 ){
            fprintf(stderr,"*** Illegal value after -pthr!\n");exit(1);
         }
         narg++ ; continue;
      }

      if( strcmp(argv[narg],"-sqr") == 0 ){
         sqr = 1 ; narg++ ; continue;
      }

      if( strcmp(argv[narg],"-arange") == 0 ){
         abot = strtod(argv[++narg],NULL) ;
         atop = strtod(argv[++narg],NULL) ;
         if( abot < 0.1 || abot > atop ){
            fprintf(stderr,"*** Illegal value after -arange!\n");exit(1);
         }
         narg++ ; continue;
      }

      if( strcmp(argv[narg],"-brange") == 0 ){
         bbot = strtod(argv[++narg],NULL) ;
         btop = strtod(argv[++narg],NULL) ;
         if( bbot < 0.1 || bbot > btop ){
            fprintf(stderr,"*** Illegal value after -brange!\n");exit(1);
         }
         narg++ ; continue;
      }

      if( strcmp(argv[narg],"-prange") == 0 ){
         pbot = strtod(argv[++narg],NULL) ;
         ptop = strtod(argv[++narg],NULL) ;
         if( pbot < 30.0 || pbot > ptop ){
            fprintf(stderr,"*** Illegal value after -prange!\n");exit(1);
         }
         narg++ ; continue;
      }

      if( strcmp(argv[narg],"-bootstrap") == 0 ){
         nboot = strtod(argv[++narg],NULL) ;
         if( nboot < 10 ){
            fprintf(stderr,"*** Illegal value after -bootstrap!\n");exit(1);
         }
         narg++ ; continue;
      }

      if( strncmp(argv[narg],"-mask",5) == 0 ){
         if( mask_dset != NULL ){
            fprintf(stderr,"*** Cannot have two -mask options!\n") ; exit(1) ;
         }
         if( narg+1 >= argc ){
            fprintf(stderr,"*** -mask option requires a following argument!\n");
            exit(1) ;
         }
         mask_dset = THD_open_dataset( argv[++narg] ) ;
         if( mask_dset == NULL ){
            fprintf(stderr,"*** Cannot open mask dataset!\n") ; exit(1) ;
         }
         if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex ){
            fprintf(stderr,"*** Cannot deal with complex-valued mask dataset!\n");
            exit(1) ;
         }
         narg++ ; continue ;
      }

      if( strncmp(argv[narg],"-mrange",5) == 0 ){
         if( narg+2 >= argc ){
            fprintf(stderr,"*** -mrange option requires 2 following arguments!\n")
;
             exit(1) ;
         }
         mask_bot = strtod( argv[++narg] , NULL ) ;
         mask_top = strtod( argv[++narg] , NULL ) ;
         if( mask_top < mask_top ){
            fprintf(stderr,"*** -mrange inputs are illegal!\n") ; exit(1) ;
         }
         narg++ ; continue ;
      }

      fprintf(stderr,"*** Illegal option: %s\n",argv[narg]) ; exit(1) ;
   }

   if( narg >= argc ){
      fprintf(stderr,"*** No dataset argument on command line!?\n");exit(1);
   }

   input_dset = THD_open_dataset( argv[narg] ) ;
   if( input_dset == NULL ){
      fprintf(stderr,"*** Can't open dataset %s\n",argv[narg]); exit(1);
   }

   nvox = DSET_NVOX(input_dset) ;

   /* load data from dataset */

   DSET_load(input_dset) ; CHECK_LOAD_ERROR(input_dset) ;

   if( DSET_BRICK_STATCODE(input_dset,0) == FUNC_COR_TYPE ) sqr = 1 ;

   bfd = BFIT_prepare_dataset( input_dset , 0 , sqr ,
                               mask_dset , 0 , mask_bot , mask_top ) ;

   if( bfd == NULL ){
      fprintf(stderr,"*** Couldn't prepare data from input dataset!\n");
      exit(1) ;
   }

   DSET_delete(mask_dset) ; DSET_delete(input_dset) ;

   for( pcut=pbot ; pcut <= ptop ; pcut += 1.0 ){
      bfr = BFIT_compute( bfd ,
                          pcut , abot,atop , bbot,btop , nran,200 ) ;
      if( bfr == NULL ){
         fprintf(stderr,"*** Can't compute betafit at pcut=%f\n",pcut) ;
         exit(1) ;
      }

      itop  = bfr->itop ;
      mgood = bfr->mgood ;

      ibot   = bfd->ibot ;
      bval   = bfd->bval ;
      cval   = bfd->cval ;
      mcount = bfd->mcount ;

      xc   = bfr->xcut ;
      aa   = bfr->a ;
      bb   = bfr->b ;
      eps  = bfr->eps ;

      xth  = beta_p2t( pthr , aa,bb ) ;
      if( sqr ) xth = sqrt(xth) ;
      if( sqr ) xc  = sqrt(xc ) ;

      if( nboot > 0 ){
         asig = bsig = tsig = abcor = 0.0 ;
         for( ii=0 ; ii < nboot ; ii++ ){
            nfd = BFIT_bootstrap_sample( bfd ) ;
            nfr = BFIT_compute( nfd ,
                                pcut , abot,atop , bbot,btop , nran,200 ) ;
            aboot = nfr->a ;
            bboot = nfr->b ;
            tboot = beta_p2t( pthr , aboot,bboot ) ;
            if( sqr ) tboot = sqrt(tboot) ;
            BFIT_free_result(nfr) ;
            BFIT_free_data  (nfd) ;
            asig += SQR( aboot-aa ) ;
            bsig += SQR( bboot-bb ) ;
            tsig += SQR( tboot-xth) ;
            abcor+= (aboot-aa)*(bboot-bb) ;
         }
         asig = sqrt(asig/nboot) ;
         bsig = sqrt(bsig/nboot) ;
         tsig = sqrt(tsig/nboot) ;
         abcor = abcor/(nboot*asig*bsig) ;
      }

      if( nboot <= 0 ){
         printf("%3.0f%%: a=%.3f b=%.3f eps=%.3f cutoff=%.3f qfit=%8.2e thr=%.3f\n",
                pcut , aa,bb,eps , xc , bfr->q_chisq , xth ) ;
      } else {
         printf("%3.0f%%: a=%.3f[%.3f] b=%.3f[%.3f] [%.3f] eps=%.3f cutoff=%.3f qfit=%8.2e thr=%.3f[%.3f]\n",
                pcut , aa,asig,bb,bsig,abcor,eps , xc , bfr->q_chisq , xth,tsig ) ;
      }
      fflush(stdout) ;

      BFIT_free_result(bfr) ;
   }
   exit(0) ;
}