Doxygen Source Code Documentation
unu.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Data Structures | |
struct | bitvec |
struct | bvarr |
Defines | |
#define | NEAR1 0.999 |
#define | NEARM1 -0.999 |
#define | NRBOT 999 |
#define | SQRT_2PI 2.5066283 |
#define | PHI(s) (1.0-0.5*normal_t2p(ss)) |
#define | bv_free(b) do{ if((b)!=NULL){free((b)->bv);free((b)->dp);free((b));} }while(0) |
#define | BITVEC_IN_BVARR(name, nn) ((name)->bvarr[(nn)]) |
#define | BVARR_SUB BITVEC_IN_BVARR |
#define | BVARR_COUNT(name) ((name)->num) |
#define | INC_BVARR 32 |
#define | INIT_BVARR(name) |
#define | ADDTO_BVARR(name, imm) |
#define | FREE_BVARR(name) |
#define | DESTROY_BVARR(name) |
#define | TRUNCATE_BVARR(name, qq) |
#define | LINEAR_DETREND |
#define | DERR(s) fprintf(stderr,"** %s\n",(s)) |
#define | IRAN(j) (lrand48() % (j)) |
#define | PROMO_MAX 4 |
Functions | |
void | set_unusuality_tail (float p) |
void | unusuality (int nr, float *rr, float *up, float *um) |
int | equal_bitvector_piece (bitvec *b, bitvec *c, int aa, int bb) |
int | equal_bitvector (bitvec *b, bitvec *c) |
void | randomize_bitvector_piece (bitvec *b, int aa, int bb) |
void | randomize_bitvector (bitvec *b) |
void | zero_bitvector_piece (bitvec *b, int aa, int bb) |
void | zero_bitvector (bitvec *b) |
void | fix_bitvector_piece (bitvec *b, int aa, int bb, int val) |
void | fix_bitvector (bitvec *b, int val) |
void | invert_bitvector_piece (bitvec *b, int aa, int bb) |
void | invert_bitvector (bitvec *b) |
bitvec * | new_bitvector (void) |
bitvec * | copy_bitvector (bitvec *b) |
int | count_bitvector (bitvec *b) |
void | normalize_floatvector (float *far) |
float | corr_floatbit (float *far, bitvec *b) |
void | evaluate_bitvec (bitvec *bim) |
void | init_floatvector_array (char *dname, char *mname) |
void | init_bitvector_array (int nbv) |
void | evolve_bitvector_array (void) |
int | main (int argc, char *argv[]) |
Variables | |
float | zstar = 0.0 |
float | pstar = 0.0 |
int | nt = 0 |
int | nfv = 0 |
int | nlev = 2 |
bvarr * | bvar = NULL |
float ** | fvar = NULL |
int | promo_ok = 0 |
Define Documentation
|
Value: do{ int nn , iq ; \ if( (name)->num == (name)->nall ){ \ nn = (name)->nall = 1.1*(name)->nall + INC_BVARR ; \ (name)->bvarr = realloc( (name)->bvarr,sizeof(bitvec *)*nn ); \ for( iq=(name)->num ; iq < (name)->nall ; iq++ ) (name)->bvarr[iq] = NULL ; } \ nn = (name)->num ; ((name)->num)++ ; \ (name)->bvarr[nn] = (imm) ; break ; } while(0) Definition at line 264 of file unu.c. Referenced by evolve_bitvector_array(), and init_bitvector_array(). |
|
|
|
Definition at line 240 of file unu.c. Referenced by evolve_bitvector_array(). |
|
Definition at line 253 of file unu.c. Referenced by evolve_bitvector_array(), and main(). |
|
Definition at line 252 of file unu.c. Referenced by evolve_bitvector_array(), and main(). |
|
Definition at line 490 of file unu.c. Referenced by init_floatvector_array(), and main(). |
|
Value: |
|
Value: |
|
|
|
Value: do{ int iq ; (name) = (bvarr *) malloc(sizeof(bvarr)) ; \ (name)->num = 0 ; (name)->nall = INC_BVARR ; \ (name)->bvarr = (bitvec **)malloc(sizeof(bitvec *)*INC_BVARR) ; \ for( iq=(name)->num ; iq < (name)->nall ; iq++ ) (name)->bvarr[iq] = NULL ; \ break ; } while(0) Definition at line 257 of file unu.c. Referenced by init_bitvector_array(). |
|
|
|
|
|
Definition at line 30 of file unu.c. Referenced by unusuality(). |
|
Definition at line 31 of file unu.c. Referenced by unusuality(). |
|
Definition at line 32 of file unu.c. Referenced by unusuality(). |
|
|
|
Definition at line 561 of file unu.c. Referenced by evolve_bitvector_array(). |
|
|
|
Value: do{ int nn ; \ if( (name) != NULL && qq < (name)->num ){ \ for( nn=qq ; nn < (name)->num ; nn++ ) bv_free((name)->bvarr[nn]); \ (name)->num = qq ; \ } } while(0) Definition at line 283 of file unu.c. Referenced by evolve_bitvector_array(). |
Function Documentation
|
Definition at line 392 of file unu.c. References bitvec::bv, c, bitvec::dp, new_bitvector(), nfv, bitvec::nlev, nt, bitvec::ubest, bitvec::um, and bitvec::up. Referenced by evolve_bitvector_array().
|
|
Definition at line 443 of file unu.c. References bitvec::bv, far, bitvec::nlev, and nt. Referenced by evaluate_bitvec(), evolve_bitvector_array(), and init_bitvector_array().
00444 { 00445 int ii , ns ; 00446 float ss , bb,bq ; 00447 byte * bar=b->bv ; 00448 00449 if( b->nlev == 2 ){ /* binary case */ 00450 for( ss=0.0,ns=ii=0 ; ii < nt ; ii++ ) 00451 if( bar[ii] ){ ns++ ; ss += far[ii] ; } 00452 if( ns == 0 || ns == nt ) return 0.0 ; 00453 ss *= sqrt( ((float) nt) / (float)(ns*(nt-ns)) ) ; 00454 00455 } else { /* multilevel case */ 00456 for( ss=bb=bq=0.0,ii=0 ; ii < nt ; ii++ ){ 00457 ss += bar[ii]*far[ii] ; 00458 bb += bar[ii] ; bq += bar[ii]*bar[ii] ; 00459 } 00460 bq -= bb*bb/nt ; if( bq <= 0.0 ) return 0.0 ; 00461 ss /= sqrt(bq) ; 00462 } 00463 00464 return ss ; 00465 } |
|
Definition at line 407 of file unu.c. References bitvec::bv, and nt.
|
|
Definition at line 302 of file unu.c. References c, equal_bitvector_piece(), and nt.
00303 { 00304 return equal_bitvector_piece( b , c , 0 , nt-1 ) ; 00305 } |
|
Definition at line 292 of file unu.c. References bitvec::bv, c, and nt. Referenced by equal_bitvector(), and evolve_bitvector_array().
|
|
Definition at line 469 of file unu.c. References corr_floatbit(), bitvec::dp, fvar, invert_bitvector(), nfv, tt, bitvec::ubest, bitvec::um, unusuality(), and bitvec::up. Referenced by evolve_bitvector_array(), and init_bitvector_array().
00470 { 00471 int jj ; 00472 00473 for( jj=0 ; jj < nfv ; jj++ ) 00474 bim->dp[jj] = corr_floatbit( fvar[jj] , bim ) ; 00475 00476 unusuality( nfv , bim->dp , &(bim->up) , &(bim->um) ) ; 00477 00478 if( bim->up < bim->um ){ 00479 float tt ; 00480 invert_bitvector( bim ) ; 00481 tt = bim->um ; bim->um = bim->up ; bim->up = tt ; 00482 } 00483 00484 bim->ubest = bim->up ; 00485 return ; 00486 } |
|
Definition at line 565 of file unu.c. References ADDTO_BVARR, bitvec::bv, bv_free, bvarr::bvarr, BVARR_COUNT, BVARR_SUB, copy_bitvector(), equal_bitvector_piece(), evaluate_bitvec(), fix_bitvector_piece(), free, invert_bitvector_piece(), IRAN, malloc, bitvec::nlev, nlev, nt, PROMO_MAX, qsort_floatstuff(), randomize_bitvector_piece(), TRUNCATE_BVARR, and zero_bitvector_piece(). Referenced by main().
00566 { 00567 static int nrvec=-1 ; 00568 static float * rvec=NULL ; 00569 00570 int ii , nbv,nbv1 , aa,bb,vv , qbv , kup ; 00571 bitvec * bim , * qim ; 00572 float aup,aum ; 00573 00574 /* promote a few to a higher plane of being */ 00575 00576 nbv1=nbv = BVARR_COUNT(bvar) ; 00577 00578 if( promo_ok ){ 00579 for( aa=ii=0 ; ii < nbv ; ii++ ) 00580 if( BVARR_SUB(bvar,ii)->nlev > nlev ) aa++ ; 00581 00582 if( aa < nbv/4 ){ 00583 for( kup=ii=0 ; ii < nbv && kup < PROMO_MAX ; ii++ ){ 00584 00585 bim = BVARR_SUB(bvar,ii) ; 00586 if( bim->nlev > nlev ) continue ; /* skip */ 00587 00588 qim = copy_bitvector(bim) ; 00589 qim->nlev *= 2 ; 00590 for( aa=0 ; aa < nt ; aa++ ) 00591 if( qim->bv[aa] < nlev-1 ) qim->bv[aa] *= 2 ; 00592 else qim->bv[aa] = 2*nlev-1 ; 00593 evaluate_bitvec( qim ) ; 00594 ADDTO_BVARR(bvar,qim) ; 00595 kup++ ; 00596 } 00597 fprintf(stderr,"%d PROMO up\n",kup) ; 00598 } else if( aa > 3*nbv/4 ){ 00599 for( kup=ii=0 ; ii < nbv && kup < PROMO_MAX ; ii++ ){ 00600 00601 bim = BVARR_SUB(bvar,ii) ; 00602 if( bim->nlev == nlev ) continue ; /* skip */ 00603 00604 qim = copy_bitvector(bim) ; 00605 qim->nlev /= 2 ; 00606 for( aa=0 ; aa < nt ; aa++ ) qim->bv[aa] /= 2 ; 00607 evaluate_bitvec( qim ) ; 00608 ADDTO_BVARR(bvar,qim) ; 00609 kup++ ; 00610 } 00611 fprintf(stderr,"%d PROMO down\n",kup) ; 00612 } 00613 } 00614 /* create mutants */ 00615 00616 qim = copy_bitvector(BVARR_SUB(bvar,0)) ; /* add copy of first one */ 00617 evaluate_bitvec( qim ) ; 00618 ADDTO_BVARR(bvar,qim) ; 00619 00620 nbv = BVARR_COUNT(bvar) ; 00621 00622 for( ii=0 ; ii < nbv ; ii++ ){ 00623 00624 bim = BVARR_SUB(bvar,ii) ; 00625 00626 aa = IRAN(nt) ; bb = aa + IRAN(5) ; if( bb >= nt ) bb = nt-1 ; 00627 00628 qim = copy_bitvector(bim) ; 00629 zero_bitvector_piece( qim , aa , bb ) ; 00630 if( equal_bitvector_piece(bim,qim,aa,bb) ){ 00631 bv_free(qim) ; 00632 } else { 00633 evaluate_bitvec( qim ) ; 00634 ADDTO_BVARR(bvar,qim) ; 00635 } 00636 00637 vv = (bim->nlev == 2) ? 1 : 1+IRAN(bim->nlev-1) ; 00638 qim = copy_bitvector(bim) ; 00639 fix_bitvector_piece( qim , aa , bb , vv ) ; 00640 if( equal_bitvector_piece(bim,qim,aa,bb) ){ 00641 bv_free(qim) ; 00642 } else { 00643 evaluate_bitvec( qim ) ; 00644 ADDTO_BVARR(bvar,qim) ; 00645 } 00646 00647 qim = copy_bitvector(bim) ; 00648 randomize_bitvector_piece( qim , aa , bb ) ; 00649 if( equal_bitvector_piece(bim,qim,aa,bb) ){ 00650 bv_free(qim) ; 00651 } else { 00652 evaluate_bitvec( qim ) ; 00653 ADDTO_BVARR(bvar,qim) ; 00654 } 00655 00656 qim = copy_bitvector(bim) ; 00657 invert_bitvector_piece( qim , aa , bb ) ; 00658 if( equal_bitvector_piece(bim,qim,aa,bb) ){ 00659 bv_free(qim) ; 00660 } else { 00661 evaluate_bitvec( qim ) ; 00662 ADDTO_BVARR(bvar,qim) ; 00663 } 00664 } 00665 00666 /* sort everybody */ 00667 00668 qbv = BVARR_COUNT(bvar) ; 00669 if( nrvec < qbv ){ 00670 if( rvec != NULL ) free(rvec) ; 00671 rvec = (float *) malloc(sizeof(float)*qbv) ; 00672 nrvec = qbv ; 00673 } 00674 for( ii=0 ; ii < qbv ; ii++ ) 00675 rvec[ii] = - BVARR_SUB(bvar,ii)->ubest ; 00676 00677 qsort_floatstuff( qbv , rvec , (void **) bvar->bvarr ) ; 00678 00679 TRUNCATE_BVARR( bvar , nbv1 ) ; 00680 return ; 00681 } |
|
Definition at line 355 of file unu.c. References fix_bitvector_piece(), and nt.
00356 { 00357 fix_bitvector_piece( b , 0 , nt-1 , val ) ; 00358 return ; 00359 } |
|
Definition at line 345 of file unu.c. References bitvec::bv, and nt. Referenced by evolve_bitvector_array(), and fix_bitvector().
|
|
Definition at line 539 of file unu.c. References ADDTO_BVARR, evaluate_bitvec(), INIT_BVARR, new_bitvector(), and randomize_bitvector(). Referenced by main().
00540 { 00541 bitvec * bim ; 00542 int ii , jj ; 00543 byte * bar ; 00544 00545 INIT_BVARR(bvar) ; 00546 00547 for( ii=0 ; ii < nbv ; ii++ ){ 00548 bim = new_bitvector() ; 00549 randomize_bitvector( bim ) ; 00550 evaluate_bitvec( bim ) ; 00551 ADDTO_BVARR(bvar,bim) ; 00552 } 00553 00554 return ; 00555 } |
|
Definition at line 492 of file unu.c. References DERR, DSET_delete, DSET_load, DSET_LOADED, DSET_NVALS, DSET_NVOX, free, fvar, malloc, mri_fix_data_pointer(), MRI_FLOAT_PTR, mri_free(), nfv, normalize_floatvector(), nt, THD_countmask(), THD_extract_series(), THD_makemask(), and THD_open_dataset(). Referenced by main().
00493 { 00494 THD_3dim_dataset * dset ; 00495 byte * mask = NULL ; 00496 int ii,jj , nvox ; 00497 MRI_IMAGE * im ; 00498 00499 dset = THD_open_dataset( dname ) ; 00500 if( dset == NULL ){ DERR("can't open dataset"); return; } 00501 nt = DSET_NVALS(dset) ; 00502 if( nt < 20 ){ DSET_delete(dset); DERR("dataset too short"); return; } 00503 DSET_load(dset) ; 00504 if( !DSET_LOADED(dset) ){ DSET_delete(dset); DERR("can't load dataset"); return; } 00505 nvox = DSET_NVOX(dset) ; 00506 00507 if( mname != NULL ){ 00508 THD_3dim_dataset * mset ; 00509 mset = THD_open_dataset( mname ) ; 00510 if( mset == NULL ){ DSET_delete(dset); DERR("can't open mask"); return; } 00511 if( DSET_NVOX(mset) != nvox ){ 00512 DSET_delete(mset); DSET_delete(dset); DERR("mask size mismatch"); return; 00513 } 00514 mask = THD_makemask( mset , 0 , 1.0,0.0 ) ; 00515 DSET_delete(mset) ; 00516 if( mask == NULL ){ DSET_delete(dset); DERR("mask is empty"); return; } 00517 nfv = THD_countmask( nvox , mask ) ; 00518 if( nfv < nt ){ DSET_delete(dset); DERR("mask is too small"); return; } 00519 } else { 00520 nfv = nvox ; 00521 } 00522 00523 fvar = (float **) malloc(sizeof(float *)*nfv) ; 00524 for( jj=ii=0 ; ii < nvox ; ii++ ){ 00525 if( mask != NULL && mask[ii] == 0 ) continue ; /* skip */ 00526 00527 im = THD_extract_series( ii , dset , 0 ) ; 00528 fvar[jj] = MRI_FLOAT_PTR(im) ; 00529 normalize_floatvector( fvar[jj] ) ; 00530 mri_fix_data_pointer(NULL,im) ; mri_free(im) ; jj++ ; 00531 } 00532 00533 if( mask != NULL ) free(mask) ; 00534 DSET_delete(dset) ; return ; 00535 } |
|
Definition at line 373 of file unu.c. References invert_bitvector_piece(), and nt. Referenced by evaluate_bitvec(), evolve_bitvector_array(), and init_bitvector_array().
00374 { 00375 invert_bitvector_piece( b , 0 , nt-1 ) ; 00376 return ; 00377 } |
|
Definition at line 363 of file unu.c. References bitvec::bv, bitvec::nlev, and nt. Referenced by evolve_bitvector_array(), and invert_bitvector().
|
|
Definition at line 685 of file unu.c. References argc, bitvec::bv, BVARR_COUNT, BVARR_SUB, DERR, evolve_bitvector_array(), fold(), fvar, init_bitvector_array(), init_floatvector_array(), nfv, nlev, nt, promo_ok, and strtod().
00686 { 00687 int ii , nv , nite=0 , neq=0 , nopt , nbv=64 ; 00688 float fold , fnew ; 00689 char * mname=NULL , * dname ; 00690 bitvec * bim ; 00691 00692 if( argc < 2 ){printf("Usage: unu [-mask mset] [-lev n] [-nbv n] dset > bname.1D\n");exit(0);} 00693 00694 nopt = 1 ; 00695 while( nopt < argc && argv[nopt][0] == '-' ){ 00696 00697 if( strcmp(argv[nopt],"-mask") == 0 ){ 00698 mname = argv[++nopt] ; 00699 nopt++ ; continue ; 00700 } 00701 00702 if( strcmp(argv[nopt],"-lev") == 0 ){ 00703 nlev = (int)strtod(argv[++nopt],NULL) ; 00704 if( nlev < 2 || nlev > 8 ){ DERR("bad -nlev"); exit(1); } 00705 nopt++ ; continue ; 00706 } 00707 00708 if( strcmp(argv[nopt],"-nbv") == 0 ){ 00709 nbv = (int)strtod(argv[++nopt],NULL) ; 00710 if( nbv < 8 || nbv > 999 ){ DERR("bad -nbv"); exit(1); } 00711 nopt++ ; continue ; 00712 } 00713 00714 fprintf(stderr,"** Illegal option %s\n",argv[nopt]); exit(1); 00715 } 00716 if( nopt >= argc ){fprintf(stderr,"** No dataset!?\n");exit(1);} 00717 00718 dname = argv[nopt] ; 00719 00720 init_floatvector_array( dname , mname ) ; 00721 if( fvar == NULL ){ 00722 fprintf(stderr,"** Couldn't init floatvector!?\n") ; exit(1) ; 00723 } else { 00724 fprintf(stderr,"nt=%d nfv=%d\n",nt,nfv) ; 00725 } 00726 00727 srand48((long)time(NULL)) ; 00728 00729 init_bitvector_array( nbv ) ; 00730 fold = BVARR_SUB(bvar,0)->ubest ; 00731 fprintf(stderr,"fold = %7.4f\n",fold) ; 00732 00733 while(1){ 00734 evolve_bitvector_array() ; 00735 nv = BVARR_COUNT(bvar) ; 00736 nite++ ; 00737 #if 1 00738 fprintf(stderr,"---nite=%d\n",nite) ; 00739 for( nopt=ii=0 ; ii < nv ; ii++ ){ 00740 fprintf(stderr," %7.4f[%d]",BVARR_SUB(bvar,ii)->ubest,BVARR_SUB(bvar,ii)->nlev) ; 00741 if( BVARR_SUB(bvar,ii)->nlev > nlev ) nopt++ ; 00742 } 00743 fprintf(stderr," *%d\n",nopt) ; 00744 #endif 00745 00746 fnew = fabs(BVARR_SUB(bvar,0)->ubest) ; 00747 if( fnew <= fold ){ 00748 neq++ ; 00749 if( neq == 8 && promo_ok ) break ; 00750 if( neq == 8 && !promo_ok ){ promo_ok = 1 ; neq = 0 ; } 00751 } else { 00752 neq = 0 ; 00753 fold = fnew ; 00754 fprintf(stderr,"%d: %7.4f\n",nite,fold) ; 00755 } 00756 } 00757 00758 bim = BVARR_SUB(bvar,0) ; 00759 for( ii=0 ; ii < nt ; ii++ ) 00760 printf(" %d\n",(int)bim->bv[ii]) ; 00761 00762 exit(0) ; 00763 } |
|
Definition at line 381 of file unu.c. References bitvec::bv, bitvec::dp, malloc, nfv, nlev, bitvec::nlev, and nt. Referenced by copy_bitvector(), and init_bitvector_array().
|
|
Definition at line 419 of file unu.c. References far, nt, SQR, and THD_linear_detrend(). Referenced by init_bitvector_array(), and init_floatvector_array().
00420 { 00421 int ii ; 00422 float ff,gg ; 00423 00424 #ifdef LINEAR_DETREND 00425 THD_linear_detrend( nt , far , NULL , NULL ) ; /* remove trend */ 00426 for( ff=0.0,ii=0 ; ii < nt ; ii++ ) ff += far[ii]*far[ii] ; /* and normalize */ 00427 if( ff <= 0.0 ) return ; 00428 ff = 1.0 / sqrt(ff) ; 00429 for( ii=0 ; ii < nt ; ii++ ) far[ii] *= ff ; 00430 #else 00431 for( ff=0.0,ii=0 ; ii < nt ; ii++ ) ff += far[ii] ; 00432 ff /= nt ; 00433 for( gg=0.0,ii=0 ; ii < nt ; ii++ ) gg += SQR( (far[ii]-ff) ) ; 00434 if( gg <= 0.0 ) return ; 00435 gg = 1.0 / sqrt(gg) ; 00436 for( ii=0 ; ii < nt ; ii++ ) far[ii] = (far[ii]-ff)*gg ; 00437 #endif 00438 return ; 00439 } |
|
Definition at line 319 of file unu.c. References nt, and randomize_bitvector_piece(). Referenced by init_bitvector_array().
00320 { 00321 randomize_bitvector_piece( b , 0 , nt-1 ) ; 00322 return ; 00323 } |
|
Definition at line 309 of file unu.c. References bitvec::bv, bitvec::nlev, and nt. Referenced by evolve_bitvector_array(), and randomize_bitvector().
|
|
Definition at line 11 of file unu.c. References p, pstar, qginv(), and zstar.
|
|
Definition at line 34 of file unu.c. References free, getenv(), malloc, MAX, NEAR1, NEARM1, nr, NRBOT, pstar, qmed_float(), set_unusuality_tail(), strtod(), and zstar.
00035 { 00036 int ii,jj , nzero , mzero ; 00037 float zmid,zsig,zmed, uval, fac, zrat, ff,fp, ss,ds,pp,ee , sigma ; 00038 float rmid , rcut , zsigt ; 00039 00040 static int nzz=0 ; /* workspace array and its size */ 00041 static float * zz=NULL ; 00042 float * aa ; 00043 00044 if( up == NULL || um == NULL ) return ; /* bad inputs */ 00045 if( nr < NRBOT || rr == NULL ){ *up=*um=0.0; return; } /* bad inputs */ 00046 00047 /*-- make workspace, if needed --*/ 00048 00049 if( nzz < 2*nr ){ 00050 if( zz != NULL ) free(zz) ; 00051 nzz = 2*nr ; 00052 zz = (float *) malloc(sizeof(float)*nzz) ; 00053 } 00054 aa = zz + nr ; /* second half */ 00055 00056 /*-- set cutoff tail, if needed --*/ 00057 00058 if( zstar <= 0.0 ){ 00059 char * cp = getenv("UTAIL") ; 00060 float pp = 0.01 ; /* default */ 00061 if( cp != NULL ){ 00062 float xx = strtod( cp , NULL ) ; 00063 if( xx > 0.0 && xx < 1.0 ) pp = xx ; 00064 } 00065 set_unusuality_tail( pp ) ; 00066 } 00067 00068 /*-- copy data into workspace, eliding values near 1 --*/ 00069 00070 for( ii=jj=0 ; ii < nr ; ii++ ) 00071 if( rr[ii] <= NEAR1 && rr[ii] >= NEARM1 ) zz[jj++] = rr[ii] ; 00072 00073 nr = jj ; 00074 if( nr < NRBOT ){ *up=*um=0.0; return; } /* shouldn't happen */ 00075 00076 /*-- find median of atanh(zz) --*/ 00077 00078 rmid = qmed_float( nr , zz ) ; /* median of correlations */ 00079 zmid = atanh(rmid) ; /* median of atanh(zz) */ 00080 00081 /*-- find MAD of atanh(zz) = median{fabs(atanh(zz)-zmid)} --*/ 00082 /* [be tricky -> use the subtraction formula for tanh] */ 00083 /* [tanh(atanh(zz)-zmid) = (zz-rmid)/(1-zz*rmid), and] */ 00084 /* [since tanh/atanh are monotonic increasing, atanh] */ 00085 /* [of median{fabs(tanh(atanh(zz)-zmid))} is the same] */ 00086 /* [as median{fabs(atanh(zz)-zmid)}. ] */ 00087 00088 for( ii=0 ; ii < nr ; ii++ ) 00089 aa[ii] = fabs( (zz[ii]-rmid)/(1.0-zz[ii]*rmid) ) ; 00090 00091 zmed = qmed_float( nr , aa ) ; /* median of aa */ 00092 zmed = atanh(zmed) ; /* MAD of atanh(zz) */ 00093 00094 zsig = 1.4826 * zmed ; /* estimate standard deviation of zz */ 00095 /* 1/1.4826 = sqrt(2)*erfinv(0.5) */ 00096 00097 if( zsig <= 0.0 ){ *up=*um=0.0; return; } /* shouldn't happen */ 00098 00099 #undef SQRT_2PI 00100 #define SQRT_2PI 2.5066283 /* sqrt(2*pi) */ 00101 00102 #undef PHI 00103 #define PHI(s) (1.0-0.5*normal_t2p(ss)) /* N(0,1) cdf */ 00104 00105 /****************************************/ 00106 /*** Compute positive tail unusuality ***/ 00107 00108 fac = 1.0 / zsig ; 00109 00110 /* Find values >= zstar, compute sum of squares */ 00111 00112 rcut = tanh( zsig * zstar + zmid ) ; /* (atanh(zz)-zmid)/zsig >= zstar */ 00113 zsigt = 0.0 ; 00114 for( mzero=ii=0 ; ii < nr ; ii++ ){ 00115 if( zz[ii] >= rcut ){ 00116 ff = fac * ( atanh(zz[ii]) - zmid ) ; /* normalized Fisher */ 00117 zsigt += ff*ff ; /* sum of squares */ 00118 mzero++ ; /* how many we get */ 00119 } 00120 } 00121 nzero = nr - mzero ; 00122 00123 /* if we don't have decent data, output is 0 */ 00124 00125 if( nzero < 2 || mzero < MAX(1.0,pstar*nr) ){ /* too weird for words */ 00126 *up = 0.0 ; 00127 } else { /* have decent data here */ 00128 zsigt = zsigt / mzero ; /* sigma-tilde squared */ 00129 00130 /* set up to compute f(s) */ 00131 00132 zrat = zstar*zstar / zsigt ; 00133 fac = ( zrat * nzero ) / ( SQRT_2PI * mzero ) ; 00134 ss = zstar ; /* initial guess for s = zstar/sigma */ 00135 00136 /* Newton's method [almost] */ 00137 00138 for( ii=0 ; ii < 5 ; ii++ ){ 00139 pp = PHI(ss) ; /* Phi(ss) \approx 1 */ 00140 ee = exp(-0.5*ss*ss) ; 00141 ff = ss*ss - (fac/pp) * ss * ee - zrat ; /* f(s) */ 00142 fp = 2.0*ss + (fac/pp) * ee * (ss*ss-1.0) ; /* f'(s) */ 00143 ds = ff / fp ; /* Newton step */ 00144 ss -= ds ; /* update */ 00145 } 00146 00147 sigma = zstar / ss ; /* estimate of sigma */ 00148 /* from upper tail data */ 00149 00150 if( sigma <= 1.0 ){ /* the boring case */ 00151 *up = 0.0 ; 00152 } else { 00153 00154 /* compute the log-likelihood difference */ 00155 00156 uval = nzero * log( PHI(ss)/(1.0-pstar) ) 00157 - mzero * ( log(sigma) + 0.5 * zsigt * (1.0/(sigma*sigma)-1.0) ) ; 00158 00159 *up = uval ; 00160 } 00161 } 00162 00163 /****************************************/ 00164 /*** Compute negative tail unusuality ***/ 00165 00166 fac = 1.0 / zsig ; 00167 00168 /* Find values <= -zstar, compute sum of squares */ 00169 00170 rcut = tanh( zmid - zsig * zstar ) ; /* (atanh(zz)-zmid)/zsig <= -zstar */ 00171 zsigt = 0.0 ; 00172 for( mzero=ii=0 ; ii < nr ; ii++ ){ 00173 if( zz[ii] <= rcut ){ 00174 ff = fac * ( atanh(zz[ii]) - zmid ) ; /* normalized Fisher */ 00175 zsigt += ff*ff ; /* sum of squares */ 00176 mzero++ ; /* how many we get */ 00177 } 00178 } 00179 nzero = nr - mzero ; 00180 00181 /* if we don't have decent data, output is 0 */ 00182 00183 if( nzero < 2 || mzero < MAX(1.0,pstar*nr) ){ /* too weird for words */ 00184 *um = 0.0 ; 00185 } else { /* have decent data here */ 00186 zsigt = zsigt / mzero ; /* sigma-tilde squared */ 00187 00188 /* set up to compute f(s) */ 00189 00190 zrat = zstar*zstar / zsigt ; 00191 fac = ( zrat * nzero ) / ( SQRT_2PI * mzero ) ; 00192 ss = zstar ; /* initial guess for s = zstar/sigma */ 00193 00194 /* Newton's method [almost] */ 00195 00196 for( ii=0 ; ii < 5 ; ii++ ){ 00197 pp = PHI(ss) ; /* Phi(ss) \approx 1 */ 00198 ee = exp(-0.5*ss*ss) ; 00199 ff = ss*ss - (fac/pp) * ss * ee - zrat ; /* f(s) */ 00200 fp = 2.0*ss + (fac/pp) * ee * (ss*ss-1.0) ; /* f'(s) */ 00201 ds = ff / fp ; /* Newton step */ 00202 ss -= ds ; /* update */ 00203 } 00204 00205 sigma = zstar / ss ; /* estimate of sigma */ 00206 /* from upper tail data */ 00207 00208 if( sigma <= 1.0 ){ /* the boring case */ 00209 *um = 0.0 ; 00210 } else { 00211 00212 /* compute the log-likelihood difference */ 00213 00214 uval = nzero * log( PHI(ss)/(1.0-pstar) ) 00215 - mzero * ( log(sigma) + 0.5 * zsigt * (1.0/(sigma*sigma)-1.0) ) ; 00216 00217 *um = uval ; 00218 } 00219 } 00220 00221 /*-- done! --*/ 00222 00223 return ; 00224 } |
|
Definition at line 337 of file unu.c. References nt, and zero_bitvector_piece().
00338 { 00339 zero_bitvector_piece( b , 0 , nt-1 ) ; 00340 return ; 00341 } |
|
Definition at line 327 of file unu.c. References bitvec::bv, and nt. Referenced by evolve_bitvector_array(), and zero_bitvector().
|
Variable Documentation
|
|
|
Definition at line 249 of file unu.c. Referenced by evaluate_bitvec(), init_floatvector_array(), and main(). |
|
Definition at line 229 of file unu.c. Referenced by copy_bitvector(), evaluate_bitvec(), init_floatvector_array(), main(), and new_bitvector(). |
|
Definition at line 230 of file unu.c. Referenced by evolve_bitvector_array(), main(), and new_bitvector(). |
|
|
Definition at line 563 of file unu.c. Referenced by main(). |
|
Definition at line 9 of file unu.c. Referenced by set_unusuality_tail(), and unusuality(). |
|
Definition at line 8 of file unu.c. Referenced by set_unusuality_tail(), and unusuality(). |