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)

do{ int nn ;                                                                 \
       if( (name) != NULL ){                                                    \
          for( nn=0 ; nn < (name)->num ; nn++ ) bv_free((name)->bvarr[nn]) ;    \
          free((name)->bvarr); free((name)); (name) = NULL; } break; } while(0)

do{ if( (name) != NULL ){                                                    \
          free((name)->bvarr); free((name)); (name) = NULL; } break; } while(0)

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)

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)

00393 {
00394    int ii ;
00395    bitvec * c = new_bitvector() ;
00396    memcpy( c->bv , b->bv , sizeof(byte)*nt   ) ;
00397 #if 0
00398    memcpy( c->dp , b->dp , sizeof(float)*nfv ) ;
00399    c->up = b->up ; c->um = b->um ; c->ubest = b->ubest ;
00400 #endif
00401    c->nlev = b->nlev ;
00402    return c ;
00403 }

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 }

00408 {
00409    int ii,ss ;
00410    byte * bv = b->bv ;
00411    for( ii=ss=0 ; ii < nt ; ii++ ) if( bv[ii] ) ss++ ;
00412    return ss ;
00413 }

00303 {
00304    return equal_bitvector_piece( b , c , 0 , nt-1 ) ;
00305 }

00293 {
00294    int ii ; byte * bv=b->bv , * cv=c->bv ;
00295 
00296    if( aa <  0  ) aa = 0 ;
00297    if( bb >= nt ) bb = nt-1 ;
00298    for( ii=aa ; ii <= bb ; ii++ ) if( bv[ii] != cv[ii] ) return 0 ;
00299    return 1;
00300 }

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 }

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 }

00356 {
00357    fix_bitvector_piece( b , 0 , nt-1 , val ) ;
00358    return ;
00359 }

00346 {
00347    int ii ; byte * bv=b->bv , vv=(byte)val ;
00348 
00349    if( aa <  0  ) aa = 0 ;
00350    if( bb >= nt ) bb = nt-1 ;
00351    for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = vv ;
00352    return ;
00353 }

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 }

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 }

00374 {
00375    invert_bitvector_piece( b , 0 , nt-1 ) ;
00376    return ;
00377 }

00364 {
00365    int ii,nl=b->nlev-1 ; byte * bv=b->bv ;
00366 
00367    if( aa <  0  ) aa = 0 ;
00368    if( bb >= nt ) bb = nt-1 ;
00369    for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = nl - bv[ii] ;
00370    return ;
00371 }

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 }

00382 {
00383    bitvec * b ;
00384    b     = (bitvec *) malloc(sizeof(bitvec)   ) ;
00385    b->bv = (byte *)   malloc(sizeof(byte) *nt ) ;
00386    b->dp = (float *)  malloc(sizeof(float)*nfv) ;
00387 
00388    b->nlev = nlev ;
00389    return b ;
00390 }

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 }

00320 {
00321    randomize_bitvector_piece( b , 0 , nt-1 ) ;
00322    return ;
00323 }

00310 {
00311    int ii ; byte * bv=b->bv ;
00312 
00313    if( aa <  0  ) aa = 0 ;
00314    if( bb >= nt ) bb = nt-1 ;
00315    for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = (byte)( b->nlev*drand48() ) ;
00316    return ;
00317 }

00012 {
00013    if( p > 0.0 && p < 1.0 ){
00014       zstar = qginv(p) ;
00015       pstar = p ;
00016    }
00017    return ;
00018 }

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 }

00338 {
00339    zero_bitvector_piece( b , 0 , nt-1 ) ;
00340    return ;
00341 }

00328 {
00329    int ii ; byte * bv=b->bv ;
00330 
00331    if( aa <  0  ) aa = 0 ;
00332    if( bb >= nt ) bb = nt-1 ;
00333    for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = 0 ;
00334    return ;
00335 }