#include "mrilib.h"
/*------------------------------------------------------------
Set the one-sided tail probability at which we will cutoff
the unusuality test.
--------------------------------------------------------------*/
static float zstar = 0.0 ; /* the actual cutoff */
static float pstar = 0.0 ; /* tail probability */
void set_unusuality_tail( float p )
{
if( p > 0.0 && p < 1.0 ){
zstar = qginv(p) ;
pstar = p ;
}
return ;
}
/*------------------------------------------------------------
Inputs: rr[0..nr-1] = array of correlation coefficients
(will not be altered)
Outputs: *up = unusuality index from positive tail of rr[]
*um = unusuality index from negative tail of rr[]
--------------------------------------------------------------*/
#undef NEAR1
#undef NEARM1
#undef NRBOT
#define NEAR1 0.999
#define NEARM1 -0.999
#define NRBOT 999
void unusuality( int nr , float * rr , float * up , float * um )
{
int ii,jj , nzero , mzero ;
float zmid,zsig,zmed, uval, fac, zrat, ff,fp, ss,ds,pp,ee , sigma ;
float rmid , rcut , zsigt ;
static int nzz=0 ; /* workspace array and its size */
static float * zz=NULL ;
float * aa ;
if( up == NULL || um == NULL ) return ; /* bad inputs */
if( nr < NRBOT || rr == NULL ){ *up=*um=0.0; return; } /* bad inputs */
/*-- make workspace, if needed --*/
if( nzz < 2*nr ){
if( zz != NULL ) free(zz) ;
nzz = 2*nr ;
zz = (float *) malloc(sizeof(float)*nzz) ;
}
aa = zz + nr ; /* second half */
/*-- set cutoff tail, if needed --*/
if( zstar <= 0.0 ){
char * cp = getenv("UTAIL") ;
float pp = 0.01 ; /* default */
if( cp != NULL ){
float xx = strtod( cp , NULL ) ;
if( xx > 0.0 && xx < 1.0 ) pp = xx ;
}
set_unusuality_tail( pp ) ;
}
/*-- copy data into workspace, eliding values near 1 --*/
for( ii=jj=0 ; ii < nr ; ii++ )
if( rr[ii] <= NEAR1 && rr[ii] >= NEARM1 ) zz[jj++] = rr[ii] ;
nr = jj ;
if( nr < NRBOT ){ *up=*um=0.0; return; } /* shouldn't happen */
/*-- find median of atanh(zz) --*/
rmid = qmed_float( nr , zz ) ; /* median of correlations */
zmid = atanh(rmid) ; /* median of atanh(zz) */
/*-- find MAD of atanh(zz) = median{fabs(atanh(zz)-zmid)} --*/
/* [be tricky -> use the subtraction formula for tanh] */
/* [tanh(atanh(zz)-zmid) = (zz-rmid)/(1-zz*rmid), and] */
/* [since tanh/atanh are monotonic increasing, atanh] */
/* [of median{fabs(tanh(atanh(zz)-zmid))} is the same] */
/* [as median{fabs(atanh(zz)-zmid)}. ] */
for( ii=0 ; ii < nr ; ii++ )
aa[ii] = fabs( (zz[ii]-rmid)/(1.0-zz[ii]*rmid) ) ;
zmed = qmed_float( nr , aa ) ; /* median of aa */
zmed = atanh(zmed) ; /* MAD of atanh(zz) */
zsig = 1.4826 * zmed ; /* estimate standard deviation of zz */
/* 1/1.4826 = sqrt(2)*erfinv(0.5) */
if( zsig <= 0.0 ){ *up=*um=0.0; return; } /* shouldn't happen */
#undef SQRT_2PI
#define SQRT_2PI 2.5066283 /* sqrt(2*pi) */
#undef PHI
#define PHI(s) (1.0-0.5*normal_t2p(ss)) /* N(0,1) cdf */
/****************************************/
/*** Compute positive tail unusuality ***/
fac = 1.0 / zsig ;
/* Find values >= zstar, compute sum of squares */
rcut = tanh( zsig * zstar + zmid ) ; /* (atanh(zz)-zmid)/zsig >= zstar */
zsigt = 0.0 ;
for( mzero=ii=0 ; ii < nr ; ii++ ){
if( zz[ii] >= rcut ){
ff = fac * ( atanh(zz[ii]) - zmid ) ; /* normalized Fisher */
zsigt += ff*ff ; /* sum of squares */
mzero++ ; /* how many we get */
}
}
nzero = nr - mzero ;
/* if we don't have decent data, output is 0 */
if( nzero < 2 || mzero < MAX(1.0,pstar*nr) ){ /* too weird for words */
*up = 0.0 ;
} else { /* have decent data here */
zsigt = zsigt / mzero ; /* sigma-tilde squared */
/* set up to compute f(s) */
zrat = zstar*zstar / zsigt ;
fac = ( zrat * nzero ) / ( SQRT_2PI * mzero ) ;
ss = zstar ; /* initial guess for s = zstar/sigma */
/* Newton's method [almost] */
for( ii=0 ; ii < 5 ; ii++ ){
pp = PHI(ss) ; /* Phi(ss) \approx 1 */
ee = exp(-0.5*ss*ss) ;
ff = ss*ss - (fac/pp) * ss * ee - zrat ; /* f(s) */
fp = 2.0*ss + (fac/pp) * ee * (ss*ss-1.0) ; /* f'(s) */
ds = ff / fp ; /* Newton step */
ss -= ds ; /* update */
}
sigma = zstar / ss ; /* estimate of sigma */
/* from upper tail data */
if( sigma <= 1.0 ){ /* the boring case */
*up = 0.0 ;
} else {
/* compute the log-likelihood difference */
uval = nzero * log( PHI(ss)/(1.0-pstar) )
- mzero * ( log(sigma) + 0.5 * zsigt * (1.0/(sigma*sigma)-1.0) ) ;
*up = uval ;
}
}
/****************************************/
/*** Compute negative tail unusuality ***/
fac = 1.0 / zsig ;
/* Find values <= -zstar, compute sum of squares */
rcut = tanh( zmid - zsig * zstar ) ; /* (atanh(zz)-zmid)/zsig <= -zstar */
zsigt = 0.0 ;
for( mzero=ii=0 ; ii < nr ; ii++ ){
if( zz[ii] <= rcut ){
ff = fac * ( atanh(zz[ii]) - zmid ) ; /* normalized Fisher */
zsigt += ff*ff ; /* sum of squares */
mzero++ ; /* how many we get */
}
}
nzero = nr - mzero ;
/* if we don't have decent data, output is 0 */
if( nzero < 2 || mzero < MAX(1.0,pstar*nr) ){ /* too weird for words */
*um = 0.0 ;
} else { /* have decent data here */
zsigt = zsigt / mzero ; /* sigma-tilde squared */
/* set up to compute f(s) */
zrat = zstar*zstar / zsigt ;
fac = ( zrat * nzero ) / ( SQRT_2PI * mzero ) ;
ss = zstar ; /* initial guess for s = zstar/sigma */
/* Newton's method [almost] */
for( ii=0 ; ii < 5 ; ii++ ){
pp = PHI(ss) ; /* Phi(ss) \approx 1 */
ee = exp(-0.5*ss*ss) ;
ff = ss*ss - (fac/pp) * ss * ee - zrat ; /* f(s) */
fp = 2.0*ss + (fac/pp) * ee * (ss*ss-1.0) ; /* f'(s) */
ds = ff / fp ; /* Newton step */
ss -= ds ; /* update */
}
sigma = zstar / ss ; /* estimate of sigma */
/* from upper tail data */
if( sigma <= 1.0 ){ /* the boring case */
*um = 0.0 ;
} else {
/* compute the log-likelihood difference */
uval = nzero * log( PHI(ss)/(1.0-pstar) )
- mzero * ( log(sigma) + 0.5 * zsigt * (1.0/(sigma*sigma)-1.0) ) ;
*um = uval ;
}
}
/*-- done! --*/
return ;
}
/***************************************************************************/
static int nt = 0 ; /* length of vectors [bitvec and float] */
static int nfv = 0 ; /* number of float vectors */
static int nlev = 2 ; /* default number of levels in bitvecs */
typedef struct {
byte * bv ; /* bitvector [nt] */
float * dp ; /* dot products [nfv] */
float up , um , ubest ; /* unusualities */
int nlev ; /* # levels */
} bitvec ;
#define bv_free(b) \
do{ if((b)!=NULL){free((b)->bv);free((b)->dp);free((b));} }while(0)
typedef struct {
int num , nall ;
bitvec ** bvarr ;
} bvarr ;
static bvarr * bvar = NULL ;
static float ** fvar = NULL ; /* nfv of these */
#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) \
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)
#define ADDTO_BVARR(name,imm) \
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)
#define FREE_BVARR(name) \
do{ if( (name) != NULL ){ \
free((name)->bvarr); free((name)); (name) = NULL; } break; } while(0)
#define DESTROY_BVARR(name) \
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)
#define TRUNCATE_BVARR(name,qq) \
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)
/***************************************************************************/
int equal_bitvector_piece( bitvec * b , bitvec * c , int aa , int bb )
{
int ii ; byte * bv=b->bv , * cv=c->bv ;
if( aa < 0 ) aa = 0 ;
if( bb >= nt ) bb = nt-1 ;
for( ii=aa ; ii <= bb ; ii++ ) if( bv[ii] != cv[ii] ) return 0 ;
return 1;
}
int equal_bitvector( bitvec * b , bitvec * c )
{
return equal_bitvector_piece( b , c , 0 , nt-1 ) ;
}
/*--------------------------------------------------------------------------*/
void randomize_bitvector_piece( bitvec * b , int aa , int bb )
{
int ii ; byte * bv=b->bv ;
if( aa < 0 ) aa = 0 ;
if( bb >= nt ) bb = nt-1 ;
for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = (byte)( b->nlev*drand48() ) ;
return ;
}
void randomize_bitvector( bitvec * b )
{
randomize_bitvector_piece( b , 0 , nt-1 ) ;
return ;
}
/*--------------------------------------------------------------------------*/
void zero_bitvector_piece( bitvec * b , int aa , int bb )
{
int ii ; byte * bv=b->bv ;
if( aa < 0 ) aa = 0 ;
if( bb >= nt ) bb = nt-1 ;
for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = 0 ;
return ;
}
void zero_bitvector( bitvec * b )
{
zero_bitvector_piece( b , 0 , nt-1 ) ;
return ;
}
/*--------------------------------------------------------------------------*/
void fix_bitvector_piece( bitvec * b , int aa , int bb , int val )
{
int ii ; byte * bv=b->bv , vv=(byte)val ;
if( aa < 0 ) aa = 0 ;
if( bb >= nt ) bb = nt-1 ;
for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = vv ;
return ;
}
void fix_bitvector( bitvec * b , int val )
{
fix_bitvector_piece( b , 0 , nt-1 , val ) ;
return ;
}
/*--------------------------------------------------------------------------*/
void invert_bitvector_piece( bitvec * b , int aa , int bb )
{
int ii,nl=b->nlev-1 ; byte * bv=b->bv ;
if( aa < 0 ) aa = 0 ;
if( bb >= nt ) bb = nt-1 ;
for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = nl - bv[ii] ;
return ;
}
void invert_bitvector( bitvec * b )
{
invert_bitvector_piece( b , 0 , nt-1 ) ;
return ;
}
/*--------------------------------------------------------------------------*/
bitvec * new_bitvector(void)
{
bitvec * b ;
b = (bitvec *) malloc(sizeof(bitvec) ) ;
b->bv = (byte *) malloc(sizeof(byte) *nt ) ;
b->dp = (float *) malloc(sizeof(float)*nfv) ;
b->nlev = nlev ;
return b ;
}
bitvec * copy_bitvector( bitvec * b )
{
int ii ;
bitvec * c = new_bitvector() ;
memcpy( c->bv , b->bv , sizeof(byte)*nt ) ;
#if 0
memcpy( c->dp , b->dp , sizeof(float)*nfv ) ;
c->up = b->up ; c->um = b->um ; c->ubest = b->ubest ;
#endif
c->nlev = b->nlev ;
return c ;
}
/*--------------------------------------------------------------------------*/
int count_bitvector( bitvec * b )
{
int ii,ss ;
byte * bv = b->bv ;
for( ii=ss=0 ; ii < nt ; ii++ ) if( bv[ii] ) ss++ ;
return ss ;
}
/*--------------------------------------------------------------------------*/
#define LINEAR_DETREND
void normalize_floatvector( float * far )
{
int ii ;
float ff,gg ;
#ifdef LINEAR_DETREND
THD_linear_detrend( nt , far , NULL , NULL ) ; /* remove trend */
for( ff=0.0,ii=0 ; ii < nt ; ii++ ) ff += far[ii]*far[ii] ; /* and normalize */
if( ff <= 0.0 ) return ;
ff = 1.0 / sqrt(ff) ;
for( ii=0 ; ii < nt ; ii++ ) far[ii] *= ff ;
#else
for( ff=0.0,ii=0 ; ii < nt ; ii++ ) ff += far[ii] ;
ff /= nt ;
for( gg=0.0,ii=0 ; ii < nt ; ii++ ) gg += SQR( (far[ii]-ff) ) ;
if( gg <= 0.0 ) return ;
gg = 1.0 / sqrt(gg) ;
for( ii=0 ; ii < nt ; ii++ ) far[ii] = (far[ii]-ff)*gg ;
#endif
return ;
}
/*--------------------------------------------------------------------------*/
float corr_floatbit( float * far , bitvec * b )
{
int ii , ns ;
float ss , bb,bq ;
byte * bar=b->bv ;
if( b->nlev == 2 ){ /* binary case */
for( ss=0.0,ns=ii=0 ; ii < nt ; ii++ )
if( bar[ii] ){ ns++ ; ss += far[ii] ; }
if( ns == 0 || ns == nt ) return 0.0 ;
ss *= sqrt( ((float) nt) / (float)(ns*(nt-ns)) ) ;
} else { /* multilevel case */
for( ss=bb=bq=0.0,ii=0 ; ii < nt ; ii++ ){
ss += bar[ii]*far[ii] ;
bb += bar[ii] ; bq += bar[ii]*bar[ii] ;
}
bq -= bb*bb/nt ; if( bq <= 0.0 ) return 0.0 ;
ss /= sqrt(bq) ;
}
return ss ;
}
/*--------------------------------------------------------------------------*/
void evaluate_bitvec( bitvec * bim )
{
int jj ;
for( jj=0 ; jj < nfv ; jj++ )
bim->dp[jj] = corr_floatbit( fvar[jj] , bim ) ;
unusuality( nfv , bim->dp , &(bim->up) , &(bim->um) ) ;
if( bim->up < bim->um ){
float tt ;
invert_bitvector( bim ) ;
tt = bim->um ; bim->um = bim->up ; bim->up = tt ;
}
bim->ubest = bim->up ;
return ;
}
/*--------------------------------------------------------------------------*/
#define DERR(s) fprintf(stderr,"** %s\n",(s))
void init_floatvector_array( char * dname , char * mname )
{
THD_3dim_dataset * dset ;
byte * mask = NULL ;
int ii,jj , nvox ;
MRI_IMAGE * im ;
dset = THD_open_dataset( dname ) ;
if( dset == NULL ){ DERR("can't open dataset"); return; }
nt = DSET_NVALS(dset) ;
if( nt < 20 ){ DSET_delete(dset); DERR("dataset too short"); return; }
DSET_load(dset) ;
if( !DSET_LOADED(dset) ){ DSET_delete(dset); DERR("can't load dataset"); return; }
nvox = DSET_NVOX(dset) ;
if( mname != NULL ){
THD_3dim_dataset * mset ;
mset = THD_open_dataset( mname ) ;
if( mset == NULL ){ DSET_delete(dset); DERR("can't open mask"); return; }
if( DSET_NVOX(mset) != nvox ){
DSET_delete(mset); DSET_delete(dset); DERR("mask size mismatch"); return;
}
mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
DSET_delete(mset) ;
if( mask == NULL ){ DSET_delete(dset); DERR("mask is empty"); return; }
nfv = THD_countmask( nvox , mask ) ;
if( nfv < nt ){ DSET_delete(dset); DERR("mask is too small"); return; }
} else {
nfv = nvox ;
}
fvar = (float **) malloc(sizeof(float *)*nfv) ;
for( jj=ii=0 ; ii < nvox ; ii++ ){
if( mask != NULL && mask[ii] == 0 ) continue ; /* skip */
im = THD_extract_series( ii , dset , 0 ) ;
fvar[jj] = MRI_FLOAT_PTR(im) ;
normalize_floatvector( fvar[jj] ) ;
mri_fix_data_pointer(NULL,im) ; mri_free(im) ; jj++ ;
}
if( mask != NULL ) free(mask) ;
DSET_delete(dset) ; return ;
}
/*--------------------------------------------------------------------------*/
void init_bitvector_array( int nbv )
{
bitvec * bim ;
int ii , jj ;
byte * bar ;
INIT_BVARR(bvar) ;
for( ii=0 ; ii < nbv ; ii++ ){
bim = new_bitvector() ;
randomize_bitvector( bim ) ;
evaluate_bitvec( bim ) ;
ADDTO_BVARR(bvar,bim) ;
}
return ;
}
/*--------------------------------------------------------------------------*/
#define IRAN(j) (lrand48() % (j))
#define PROMO_MAX 4
static int promo_ok = 0 ;
void evolve_bitvector_array(void)
{
static int nrvec=-1 ;
static float * rvec=NULL ;
int ii , nbv,nbv1 , aa,bb,vv , qbv , kup ;
bitvec * bim , * qim ;
float aup,aum ;
/* promote a few to a higher plane of being */
nbv1=nbv = BVARR_COUNT(bvar) ;
if( promo_ok ){
for( aa=ii=0 ; ii < nbv ; ii++ )
if( BVARR_SUB(bvar,ii)->nlev > nlev ) aa++ ;
if( aa < nbv/4 ){
for( kup=ii=0 ; ii < nbv && kup < PROMO_MAX ; ii++ ){
bim = BVARR_SUB(bvar,ii) ;
if( bim->nlev > nlev ) continue ; /* skip */
qim = copy_bitvector(bim) ;
qim->nlev *= 2 ;
for( aa=0 ; aa < nt ; aa++ )
if( qim->bv[aa] < nlev-1 ) qim->bv[aa] *= 2 ;
else qim->bv[aa] = 2*nlev-1 ;
evaluate_bitvec( qim ) ;
ADDTO_BVARR(bvar,qim) ;
kup++ ;
}
fprintf(stderr,"%d PROMO up\n",kup) ;
} else if( aa > 3*nbv/4 ){
for( kup=ii=0 ; ii < nbv && kup < PROMO_MAX ; ii++ ){
bim = BVARR_SUB(bvar,ii) ;
if( bim->nlev == nlev ) continue ; /* skip */
qim = copy_bitvector(bim) ;
qim->nlev /= 2 ;
for( aa=0 ; aa < nt ; aa++ ) qim->bv[aa] /= 2 ;
evaluate_bitvec( qim ) ;
ADDTO_BVARR(bvar,qim) ;
kup++ ;
}
fprintf(stderr,"%d PROMO down\n",kup) ;
}
}
/* create mutants */
qim = copy_bitvector(BVARR_SUB(bvar,0)) ; /* add copy of first one */
evaluate_bitvec( qim ) ;
ADDTO_BVARR(bvar,qim) ;
nbv = BVARR_COUNT(bvar) ;
for( ii=0 ; ii < nbv ; ii++ ){
bim = BVARR_SUB(bvar,ii) ;
aa = IRAN(nt) ; bb = aa + IRAN(5) ; if( bb >= nt ) bb = nt-1 ;
qim = copy_bitvector(bim) ;
zero_bitvector_piece( qim , aa , bb ) ;
if( equal_bitvector_piece(bim,qim,aa,bb) ){
bv_free(qim) ;
} else {
evaluate_bitvec( qim ) ;
ADDTO_BVARR(bvar,qim) ;
}
vv = (bim->nlev == 2) ? 1 : 1+IRAN(bim->nlev-1) ;
qim = copy_bitvector(bim) ;
fix_bitvector_piece( qim , aa , bb , vv ) ;
if( equal_bitvector_piece(bim,qim,aa,bb) ){
bv_free(qim) ;
} else {
evaluate_bitvec( qim ) ;
ADDTO_BVARR(bvar,qim) ;
}
qim = copy_bitvector(bim) ;
randomize_bitvector_piece( qim , aa , bb ) ;
if( equal_bitvector_piece(bim,qim,aa,bb) ){
bv_free(qim) ;
} else {
evaluate_bitvec( qim ) ;
ADDTO_BVARR(bvar,qim) ;
}
qim = copy_bitvector(bim) ;
invert_bitvector_piece( qim , aa , bb ) ;
if( equal_bitvector_piece(bim,qim,aa,bb) ){
bv_free(qim) ;
} else {
evaluate_bitvec( qim ) ;
ADDTO_BVARR(bvar,qim) ;
}
}
/* sort everybody */
qbv = BVARR_COUNT(bvar) ;
if( nrvec < qbv ){
if( rvec != NULL ) free(rvec) ;
rvec = (float *) malloc(sizeof(float)*qbv) ;
nrvec = qbv ;
}
for( ii=0 ; ii < qbv ; ii++ )
rvec[ii] = - BVARR_SUB(bvar,ii)->ubest ;
qsort_floatstuff( qbv , rvec , (void **) bvar->bvarr ) ;
TRUNCATE_BVARR( bvar , nbv1 ) ;
return ;
}
/*--------------------------------------------------------------------------*/
int main( int argc , char * argv[] )
{
int ii , nv , nite=0 , neq=0 , nopt , nbv=64 ;
float fold , fnew ;
char * mname=NULL , * dname ;
bitvec * bim ;
if( argc < 2 ){printf("Usage: unu [-mask mset] [-lev n] [-nbv n] dset > bname.1D\n");exit(0);}
nopt = 1 ;
while( nopt < argc && argv[nopt][0] == '-' ){
if( strcmp(argv[nopt],"-mask") == 0 ){
mname = argv[++nopt] ;
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-lev") == 0 ){
nlev = (int)strtod(argv[++nopt],NULL) ;
if( nlev < 2 || nlev > 8 ){ DERR("bad -nlev"); exit(1); }
nopt++ ; continue ;
}
if( strcmp(argv[nopt],"-nbv") == 0 ){
nbv = (int)strtod(argv[++nopt],NULL) ;
if( nbv < 8 || nbv > 999 ){ DERR("bad -nbv"); exit(1); }
nopt++ ; continue ;
}
fprintf(stderr,"** Illegal option %s\n",argv[nopt]); exit(1);
}
if( nopt >= argc ){fprintf(stderr,"** No dataset!?\n");exit(1);}
dname = argv[nopt] ;
init_floatvector_array( dname , mname ) ;
if( fvar == NULL ){
fprintf(stderr,"** Couldn't init floatvector!?\n") ; exit(1) ;
} else {
fprintf(stderr,"nt=%d nfv=%d\n",nt,nfv) ;
}
srand48((long)time(NULL)) ;
init_bitvector_array( nbv ) ;
fold = BVARR_SUB(bvar,0)->ubest ;
fprintf(stderr,"fold = %7.4f\n",fold) ;
while(1){
evolve_bitvector_array() ;
nv = BVARR_COUNT(bvar) ;
nite++ ;
#if 1
fprintf(stderr,"---nite=%d\n",nite) ;
for( nopt=ii=0 ; ii < nv ; ii++ ){
fprintf(stderr," %7.4f[%d]",BVARR_SUB(bvar,ii)->ubest,BVARR_SUB(bvar,ii)->nlev) ;
if( BVARR_SUB(bvar,ii)->nlev > nlev ) nopt++ ;
}
fprintf(stderr," *%d\n",nopt) ;
#endif
fnew = fabs(BVARR_SUB(bvar,0)->ubest) ;
if( fnew <= fold ){
neq++ ;
if( neq == 8 && promo_ok ) break ;
if( neq == 8 && !promo_ok ){ promo_ok = 1 ; neq = 0 ; }
} else {
neq = 0 ;
fold = fnew ;
fprintf(stderr,"%d: %7.4f\n",nite,fold) ;
}
}
bim = BVARR_SUB(bvar,0) ;
for( ii=0 ; ii < nt ; ii++ )
printf(" %d\n",(int)bim->bv[ii]) ;
exit(0) ;
}
syntax highlighted by Code2HTML, v. 0.9.1