00001 #include "mrilib.h"
00002
00003
00004
00005
00006
00007
00008 static float zstar = 0.0 ;
00009 static float pstar = 0.0 ;
00010
00011 void set_unusuality_tail( float p )
00012 {
00013 if( p > 0.0 && p < 1.0 ){
00014 zstar = qginv(p) ;
00015 pstar = p ;
00016 }
00017 return ;
00018 }
00019
00020
00021
00022
00023
00024
00025
00026
00027 #undef NEAR1
00028 #undef NEARM1
00029 #undef NRBOT
00030 #define NEAR1 0.999
00031 #define NEARM1 -0.999
00032 #define NRBOT 999
00033
00034 void unusuality( int nr , float * rr , float * up , float * um )
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 ;
00041 static float * zz=NULL ;
00042 float * aa ;
00043
00044 if( up == NULL || um == NULL ) return ;
00045 if( nr < NRBOT || rr == NULL ){ *up=*um=0.0; return; }
00046
00047
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 ;
00055
00056
00057
00058 if( zstar <= 0.0 ){
00059 char * cp = getenv("UTAIL") ;
00060 float pp = 0.01 ;
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
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; }
00075
00076
00077
00078 rmid = qmed_float( nr , zz ) ;
00079 zmid = atanh(rmid) ;
00080
00081
00082
00083
00084
00085
00086
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 ) ;
00092 zmed = atanh(zmed) ;
00093
00094 zsig = 1.4826 * zmed ;
00095
00096
00097 if( zsig <= 0.0 ){ *up=*um=0.0; return; }
00098
00099 #undef SQRT_2PI
00100 #define SQRT_2PI 2.5066283
00101
00102 #undef PHI
00103 #define PHI(s) (1.0-0.5*normal_t2p(ss))
00104
00105
00106
00107
00108 fac = 1.0 / zsig ;
00109
00110
00111
00112 rcut = tanh( zsig * zstar + zmid ) ;
00113 zsigt = 0.0 ;
00114 for( mzero=ii=0 ; ii < nr ; ii++ ){
00115 if( zz[ii] >= rcut ){
00116 ff = fac * ( atanh(zz[ii]) - zmid ) ;
00117 zsigt += ff*ff ;
00118 mzero++ ;
00119 }
00120 }
00121 nzero = nr - mzero ;
00122
00123
00124
00125 if( nzero < 2 || mzero < MAX(1.0,pstar*nr) ){
00126 *up = 0.0 ;
00127 } else {
00128 zsigt = zsigt / mzero ;
00129
00130
00131
00132 zrat = zstar*zstar / zsigt ;
00133 fac = ( zrat * nzero ) / ( SQRT_2PI * mzero ) ;
00134 ss = zstar ;
00135
00136
00137
00138 for( ii=0 ; ii < 5 ; ii++ ){
00139 pp = PHI(ss) ;
00140 ee = exp(-0.5*ss*ss) ;
00141 ff = ss*ss - (fac/pp) * ss * ee - zrat ;
00142 fp = 2.0*ss + (fac/pp) * ee * (ss*ss-1.0) ;
00143 ds = ff / fp ;
00144 ss -= ds ;
00145 }
00146
00147 sigma = zstar / ss ;
00148
00149
00150 if( sigma <= 1.0 ){
00151 *up = 0.0 ;
00152 } else {
00153
00154
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
00165
00166 fac = 1.0 / zsig ;
00167
00168
00169
00170 rcut = tanh( 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 ) ;
00175 zsigt += ff*ff ;
00176 mzero++ ;
00177 }
00178 }
00179 nzero = nr - mzero ;
00180
00181
00182
00183 if( nzero < 2 || mzero < MAX(1.0,pstar*nr) ){
00184 *um = 0.0 ;
00185 } else {
00186 zsigt = zsigt / mzero ;
00187
00188
00189
00190 zrat = zstar*zstar / zsigt ;
00191 fac = ( zrat * nzero ) / ( SQRT_2PI * mzero ) ;
00192 ss = zstar ;
00193
00194
00195
00196 for( ii=0 ; ii < 5 ; ii++ ){
00197 pp = PHI(ss) ;
00198 ee = exp(-0.5*ss*ss) ;
00199 ff = ss*ss - (fac/pp) * ss * ee - zrat ;
00200 fp = 2.0*ss + (fac/pp) * ee * (ss*ss-1.0) ;
00201 ds = ff / fp ;
00202 ss -= ds ;
00203 }
00204
00205 sigma = zstar / ss ;
00206
00207
00208 if( sigma <= 1.0 ){
00209 *um = 0.0 ;
00210 } else {
00211
00212
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
00222
00223 return ;
00224 }
00225
00226
00227
00228 static int nt = 0 ;
00229 static int nfv = 0 ;
00230 static int nlev = 2 ;
00231
00232 typedef struct {
00233 byte * bv ;
00234 float * dp ;
00235
00236 float up , um , ubest ;
00237 int nlev ;
00238 } bitvec ;
00239
00240 #define bv_free(b) \
00241 do{ if((b)!=NULL){free((b)->bv);free((b)->dp);free((b));} }while(0)
00242
00243 typedef struct {
00244 int num , nall ;
00245 bitvec ** bvarr ;
00246 } bvarr ;
00247
00248 static bvarr * bvar = NULL ;
00249 static float ** fvar = NULL ;
00250
00251 #define BITVEC_IN_BVARR(name,nn) ((name)->bvarr[(nn)])
00252 #define BVARR_SUB BITVEC_IN_BVARR
00253 #define BVARR_COUNT(name) ((name)->num)
00254
00255 #define INC_BVARR 32
00256
00257 #define INIT_BVARR(name) \
00258 do{ int iq ; (name) = (bvarr *) malloc(sizeof(bvarr)) ; \
00259 (name)->num = 0 ; (name)->nall = INC_BVARR ; \
00260 (name)->bvarr = (bitvec **)malloc(sizeof(bitvec *)*INC_BVARR) ; \
00261 for( iq=(name)->num ; iq < (name)->nall ; iq++ ) (name)->bvarr[iq] = NULL ; \
00262 break ; } while(0)
00263
00264 #define ADDTO_BVARR(name,imm) \
00265 do{ int nn , iq ; \
00266 if( (name)->num == (name)->nall ){ \
00267 nn = (name)->nall = 1.1*(name)->nall + INC_BVARR ; \
00268 (name)->bvarr = realloc( (name)->bvarr,sizeof(bitvec *)*nn ); \
00269 for( iq=(name)->num ; iq < (name)->nall ; iq++ ) (name)->bvarr[iq] = NULL ; } \
00270 nn = (name)->num ; ((name)->num)++ ; \
00271 (name)->bvarr[nn] = (imm) ; break ; } while(0)
00272
00273 #define FREE_BVARR(name) \
00274 do{ if( (name) != NULL ){ \
00275 free((name)->bvarr); free((name)); (name) = NULL; } break; } while(0)
00276
00277 #define DESTROY_BVARR(name) \
00278 do{ int nn ; \
00279 if( (name) != NULL ){ \
00280 for( nn=0 ; nn < (name)->num ; nn++ ) bv_free((name)->bvarr[nn]) ; \
00281 free((name)->bvarr); free((name)); (name) = NULL; } break; } while(0)
00282
00283 #define TRUNCATE_BVARR(name,qq) \
00284 do{ int nn ; \
00285 if( (name) != NULL && qq < (name)->num ){ \
00286 for( nn=qq ; nn < (name)->num ; nn++ ) bv_free((name)->bvarr[nn]); \
00287 (name)->num = qq ; \
00288 } } while(0)
00289
00290
00291
00292 int equal_bitvector_piece( bitvec * b , bitvec * c , int aa , int bb )
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 }
00301
00302 int equal_bitvector( bitvec * b , bitvec * c )
00303 {
00304 return equal_bitvector_piece( b , c , 0 , nt-1 ) ;
00305 }
00306
00307
00308
00309 void randomize_bitvector_piece( bitvec * b , int aa , int bb )
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 }
00318
00319 void randomize_bitvector( bitvec * b )
00320 {
00321 randomize_bitvector_piece( b , 0 , nt-1 ) ;
00322 return ;
00323 }
00324
00325
00326
00327 void zero_bitvector_piece( bitvec * b , int aa , int bb )
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 }
00336
00337 void zero_bitvector( bitvec * b )
00338 {
00339 zero_bitvector_piece( b , 0 , nt-1 ) ;
00340 return ;
00341 }
00342
00343
00344
00345 void fix_bitvector_piece( bitvec * b , int aa , int bb , int val )
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 }
00354
00355 void fix_bitvector( bitvec * b , int val )
00356 {
00357 fix_bitvector_piece( b , 0 , nt-1 , val ) ;
00358 return ;
00359 }
00360
00361
00362
00363 void invert_bitvector_piece( bitvec * b , int aa , int bb )
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 }
00372
00373 void invert_bitvector( bitvec * b )
00374 {
00375 invert_bitvector_piece( b , 0 , nt-1 ) ;
00376 return ;
00377 }
00378
00379
00380
00381 bitvec * new_bitvector(void)
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 }
00391
00392 bitvec * copy_bitvector( bitvec * b )
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 }
00404
00405
00406
00407 int count_bitvector( bitvec * b )
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 }
00414
00415
00416
00417 #define LINEAR_DETREND
00418
00419 void normalize_floatvector( float * far )
00420 {
00421 int ii ;
00422 float ff,gg ;
00423
00424 #ifdef LINEAR_DETREND
00425 THD_linear_detrend( nt , far , NULL , NULL ) ;
00426 for( ff=0.0,ii=0 ; ii < nt ; ii++ ) ff += far[ii]*far[ii] ;
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 }
00440
00441
00442
00443 float corr_floatbit( float * far , bitvec * b )
00444 {
00445 int ii , ns ;
00446 float ss , bb,bq ;
00447 byte * bar=b->bv ;
00448
00449 if( b->nlev == 2 ){
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 {
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 }
00466
00467
00468
00469 void evaluate_bitvec( bitvec * bim )
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 }
00487
00488
00489
00490 #define DERR(s) fprintf(stderr,"** %s\n",(s))
00491
00492 void init_floatvector_array( char * dname , char * mname )
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 ;
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 }
00536
00537
00538
00539 void init_bitvector_array( int nbv )
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 }
00556
00557
00558
00559 #define IRAN(j) (lrand48() % (j))
00560
00561 #define PROMO_MAX 4
00562
00563 static int promo_ok = 0 ;
00564
00565 void evolve_bitvector_array(void)
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
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 ;
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 ;
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
00615
00616 qim = copy_bitvector(BVARR_SUB(bvar,0)) ;
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
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 }
00682
00683
00684
00685 int main( int argc , char * argv[] )
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 }