Doxygen Source Code Documentation
3dbetafit2.c File Reference
#include "mrilib.h"
#include "betafit.c"
Go to the source code of this file.
Data Structures | |
struct | covmat |
struct | spanfit |
Defines | |
#define | IFREE(x) do{if((x)!=NULL)free(x);}while(0) |
#define | FREE_COVMAT(cc) |
#define | CM(i, j) cmat[(i)+(j)*ndim] |
#define | CH(i, j) cfac[(i)+(j)*ndim] |
#define | CCUT 3.5 |
#define | EPS 1.e-4 |
#define | DDD csum |
#define | OUT_THR 1 |
#define | OUT_BBB 2 |
#define | OUT_AAA 3 |
Functions | |
void | forward_solve_inplace (covmat *cv, float *vec) |
void | compute_choleski (covmat *cv) |
covmat * | robust_covar (int ndim, int nvec, float **vec) |
float | evaluate_span (int ndim, int nvec, int bot, int top, float *cvec, float **bvec) |
spanfit | find_best_span (int ndim, int nvec, int minspan, float *cvec, float **bvec) |
float | process_sample (float pcut, BFIT_data *bfd) |
int | main (int argc, char *argv[]) |
Variables | |
int | nran = 1000 |
float | abot = 0.5 |
float | atop = 4.0 |
float | bbot = 10.0 |
float | btop = 200.0 |
float | pbot = 50.0 |
float | ptop = 80.0 |
double | pthr = 1.e-4 |
int | sqr = 0 |
int | outmode = OUT_THR |
Define Documentation
|
Definition at line 95 of file 3dbetafit2.c. Referenced by robust_covar(). |
|
Definition at line 30 of file 3dbetafit2.c. Referenced by backward_solve_inplace(), compute_choleski(), and forward_solve_inplace(). |
|
Definition at line 29 of file 3dbetafit2.c. Referenced by compute_choleski(). |
|
|
|
Definition at line 96 of file 3dbetafit2.c. |
|
Value: Definition at line 24 of file 3dbetafit2.c. Referenced by evaluate_span(). |
|
Definition at line 22 of file 3dbetafit2.c. |
|
Definition at line 452 of file 3dbetafit2.c. Referenced by main(), and process_sample(). |
|
Definition at line 451 of file 3dbetafit2.c. Referenced by main(), and process_sample(). |
|
Definition at line 450 of file 3dbetafit2.c. Referenced by process_sample(). |
Function Documentation
|
Definition at line 66 of file 3dbetafit2.c. References covmat::cfac, CH, CM, covmat::cmat, free, malloc, and covmat::ndim.
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 } |
|
Definition at line 337 of file 3dbetafit2.c. References forward_solve_inplace(), free, FREE_COVMAT, malloc, robust_covar(), and top.
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 } |
|
Definition at line 404 of file 3dbetafit2.c. References spanfit::bot, evaluate_span(), spanfit::pval, spanfit::top, and top.
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 } |
|
Definition at line 34 of file 3dbetafit2.c. References covmat::cfac, CH, covmat::ndim, and vec.
|
|
compute the overall minimum and maximum voxel values for a dataset Definition at line 507 of file 3dbetafit2.c. References abot, argc, atop, bbot, BFIT_bootstrap_sample(), BFIT_free_data(), BFIT_prepare_dataset(), spanfit::bot, btop, DSET_BRICK_STATCODE, DSET_BRICK_TYPE, DSET_delete, DSET_load, DSET_LOADED, DSET_NVOX, evaluate_span(), find_best_span(), malloc, mmm, OUT_AAA, OUT_BBB, outmode, pbot, process_sample(), pthr, ptop, spanfit::pval, SQR, sqr, strtod(), THD_open_dataset(), spanfit::top, and xc.
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 } |
|
Definition at line 458 of file 3dbetafit2.c. References BFIT_result::a, abot, atop, BFIT_result::b, bbot, beta_init(), beta_p2t(), BFIT_compute(), BFIT_free_result(), btop, nran, OUT_AAA, OUT_BBB, OUT_THR, and pthr.
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 } |
|
Definition at line 98 of file 3dbetafit2.c. References CCUT, covmat::cfac, covmat::cmat, compute_choleski(), forward_solve_inplace(), free, malloc, covmat::mvec, covmat::ndim, and 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 } |
Variable Documentation
|
Definition at line 444 of file 3dbetafit2.c. Referenced by main(), and process_sample(). |
|
Definition at line 444 of file 3dbetafit2.c. Referenced by main(), and process_sample(). |
|
Definition at line 445 of file 3dbetafit2.c. Referenced by main(), and process_sample(). |
|
Definition at line 445 of file 3dbetafit2.c. Referenced by main(), and process_sample(). |
|
Definition at line 443 of file 3dbetafit2.c. Referenced by process_sample(). |
|
Definition at line 454 of file 3dbetafit2.c. Referenced by main(). |
|
Definition at line 446 of file 3dbetafit2.c. Referenced by main(). |
|
Definition at line 447 of file 3dbetafit2.c. Referenced by main(), and process_sample(). |
|
Definition at line 446 of file 3dbetafit2.c. Referenced by main(). |
|
Definition at line 448 of file 3dbetafit2.c. Referenced by main(). |