do{ if(cc != NULL){                            \
        IFREE(cc->cmat); IFREE(cc->cfac);        \
        IFREE(cc->mvec); free(cc); } } while(0)

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 }

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 }

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    /* make pointers to subvectors of bvec */
00273 
00274    svec = (float **) malloc(sizeof(float *)*nvec) ;
00275    for( kk=0 ; kk < nvec ; kk++ ) svec[kk] = bvec[kk] + bot ;
00276 
00277    /* estimate covariance of subvectors */
00278 
00279    cv = robust_covar( npt , nvec , svec ) ;
00280    free(svec) ;
00281    if( cv == NULL ) return 0.0 ;                  /* shouldn't happen */
00282 
00283    /* compute normalized cvec and basis into xx, ee */
00284 
00285 #if 0
00286 fprintf(stderr,"Normalizing cvec, basis\n") ;
00287 #endif
00288 
00289    ee = (float **) malloc(sizeof(float *)*nbasis) ;    /* make space */
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++ ){                      /* copy in */
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 ) ;                  /* normalize */
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    /*             T     */
00306    /* compute [ee] [ee] */
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 ;                  /* won't be used */
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    /* project [xx] onto space orthogonal to [ee] */
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    /* normalize each bvec, dot into residual vector, count negatives */
00342    /* (don't have to project bvec onto space orthog)                */
00343    /* (to [ee], since [xx] is already in that space)               */
00344 
00345    nbd = 0 ;
00346    for( kk=0 ; kk < nvec ; kk++ ){
00347       memcpy( ee[0] , bvec[kk]+bot , sizeof(float)*npt ) ; /* bvec */
00348       forward_solve_inplace( cv , ee[0] ) ;                /* normalized */
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 }

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 }

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 }

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    /* scan command-line args */
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    /* load data from dataset */
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 }

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: /* use the threshold as the output parameter */
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 }

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) ;  /* matrix     */
00122    tv   = (float *) malloc(sizeof(float)*ndim) ;       /* temp vector */
00123    mv   = (float *) malloc(sizeof(float)*ndim) ;       /* mean vector  */
00124    wv   = (float *) malloc(sizeof(float)*nvec) ;       /* weight vector */
00125 
00126    fnvec = 1.0/nvec ; fndim = 1.0/ndim ;
00127    bcut  = 1.0 + CCUT*sqrt(fndim) ;
00128 
00129    /* compute initial mean & covariance matrix with all weights = 1 */
00130 
00131    for( jj=0 ; jj < ndim ; jj++ ) mv[jj] = 0.0 ;
00132 
00133    for( kk=0 ; kk < nvec ; kk++ ){   /* mean vector sum */
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 ;  /* scale mean vector */
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++ ){   /* covariance matrix sum */
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++ ){   /* scale covariance matrix */
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    /* now iterate until convergence, or something */
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 ;  /* put old matrix into cv */
00168       cv->mvec = mv ;           /* and old mean vector   */
00169       compute_choleski(cv) ;    /* decompose matrix     */
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) ; /* new matrix */
00177       mv   = (float *) malloc(sizeof(float)*ndim) ;      /* new mean vector */
00178 
00179       for( jj=0 ; jj < ndim ; jj++ ){  /* initialize new things to zero */
00180          mv[jj] = 0.0 ;
00181          for( ii=0 ; ii < ndim ; ii++ ) nmat[ii+jj*ndim] = 0.0 ;
00182       }
00183 
00184       /* update mean */
00185 
00186       csum = 0.0 ;
00187       for( kk=0 ; kk < nvec ; kk++ ){
00188          vv = vec[kk] ;
00189 
00190          /*                    -1/2          */
00191          /* compute tv = [cmat]    (vv-mvec) */
00192 
00193          for( jj=0 ; jj < ndim ; jj++ ) tv[jj] = vv[jj] - cv->mvec[jj] ;
00194          forward_solve_inplace(cv,tv) ;
00195 
00196          /* compute norm of tv, then weighting factor for this vector */
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          /* add vv into accumulating mean, with weight cnorm */
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 ;  /* scale new mean */
00209 
00210       /* update covariance */
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       /* check for convergence - L1 norm */
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;  /* exit loop */
00253       }
00254    }
00255 
00256    free(wv) ; free(tv) ; compute_choleski(cv) ; return cv ;
00257 }