1 #include "mrilib.h"
2
3 /*------------------------------------------------------------
4 Set the one-sided tail probability at which we will cutoff
5 the unusuality test.
6 --------------------------------------------------------------*/
7
8 static float zstar = 0.0 ; /* the actual cutoff */
9 static float pstar = 0.0 ; /* tail probability */
10
set_unusuality_tail(float p)11 void set_unusuality_tail( float p )
12 {
13 if( p > 0.0 && p < 1.0 ){
14 zstar = qginv(p) ;
15 pstar = p ;
16 }
17 return ;
18 }
19
20 /*------------------------------------------------------------
21 Inputs: rr[0..nr-1] = array of correlation coefficients
22 (will not be altered)
23 Outputs: *up = unusuality index from positive tail of rr[]
24 *um = unusuality index from negative tail of rr[]
25 --------------------------------------------------------------*/
26
27 #undef NEAR1
28 #undef NEARM1
29 #undef NRBOT
30 #define NEAR1 0.999
31 #define NEARM1 -0.999
32 #define NRBOT 999
33
unusuality(int nr,float * rr,float * up,float * um)34 void unusuality( int nr , float * rr , float * up , float * um )
35 {
36 int ii,jj , nzero , mzero ;
37 float zmid,zsig,zmed, uval, fac, zrat, ff,fp, ss,ds,pp,ee , sigma ;
38 float rmid , rcut , zsigt ;
39
40 static int nzz=0 ; /* workspace array and its size */
41 static float * zz=NULL ;
42 float * aa ;
43
44 if( up == NULL || um == NULL ) return ; /* bad inputs */
45 if( nr < NRBOT || rr == NULL ){ *up=*um=0.0; return; } /* bad inputs */
46
47 /*-- make workspace, if needed --*/
48
49 if( nzz < 2*nr ){
50 if( zz != NULL ) free(zz) ;
51 nzz = 2*nr ;
52 zz = (float *) malloc(sizeof(float)*nzz) ;
53 }
54 aa = zz + nr ; /* second half */
55
56 /*-- set cutoff tail, if needed --*/
57
58 if( zstar <= 0.0 ){
59 char * cp = getenv("UTAIL") ;
60 float pp = 0.01 ; /* default */
61 if( cp != NULL ){
62 float xx = strtod( cp , NULL ) ;
63 if( xx > 0.0 && xx < 1.0 ) pp = xx ;
64 }
65 set_unusuality_tail( pp ) ;
66 }
67
68 /*-- copy data into workspace, eliding values near 1 --*/
69
70 for( ii=jj=0 ; ii < nr ; ii++ )
71 if( rr[ii] <= NEAR1 && rr[ii] >= NEARM1 ) zz[jj++] = rr[ii] ;
72
73 nr = jj ;
74 if( nr < NRBOT ){ *up=*um=0.0; return; } /* shouldn't happen */
75
76 /*-- find median of atanh(zz) --*/
77
78 rmid = qmed_float( nr , zz ) ; /* median of correlations */
79 zmid = atanh(rmid) ; /* median of atanh(zz) */
80
81 /*-- find MAD of atanh(zz) = median{fabs(atanh(zz)-zmid)} --*/
82 /* [be tricky -> use the subtraction formula for tanh] */
83 /* [tanh(atanh(zz)-zmid) = (zz-rmid)/(1-zz*rmid), and] */
84 /* [since tanh/atanh are monotonic increasing, atanh] */
85 /* [of median{fabs(tanh(atanh(zz)-zmid))} is the same] */
86 /* [as median{fabs(atanh(zz)-zmid)}. ] */
87
88 for( ii=0 ; ii < nr ; ii++ )
89 aa[ii] = fabs( (zz[ii]-rmid)/(1.0-zz[ii]*rmid) ) ;
90
91 zmed = qmed_float( nr , aa ) ; /* median of aa */
92 zmed = atanh(zmed) ; /* MAD of atanh(zz) */
93
94 zsig = 1.4826 * zmed ; /* estimate standard deviation of zz */
95 /* 1/1.4826 = sqrt(2)*erfinv(0.5) */
96
97 if( zsig <= 0.0 ){ *up=*um=0.0; return; } /* shouldn't happen */
98
99 #undef SQRT_2PI
100 #define SQRT_2PI 2.5066283 /* sqrt(2*pi) */
101
102 #undef PHI
103 #define PHI(s) (1.0-0.5*normal_t2p(ss)) /* N(0,1) cdf */
104
105 /****************************************/
106 /*** Compute positive tail unusuality ***/
107
108 fac = 1.0 / zsig ;
109
110 /* Find values >= zstar, compute sum of squares */
111
112 rcut = tanh( zsig * zstar + zmid ) ; /* (atanh(zz)-zmid)/zsig >= zstar */
113 zsigt = 0.0 ;
114 for( mzero=ii=0 ; ii < nr ; ii++ ){
115 if( zz[ii] >= rcut ){
116 ff = fac * ( atanh(zz[ii]) - zmid ) ; /* normalized Fisher */
117 zsigt += ff*ff ; /* sum of squares */
118 mzero++ ; /* how many we get */
119 }
120 }
121 nzero = nr - mzero ;
122
123 /* if we don't have decent data, output is 0 */
124
125 if( nzero < 2 || mzero < MAX(1.0,pstar*nr) ){ /* too weird for words */
126 *up = 0.0 ;
127 } else { /* have decent data here */
128 zsigt = zsigt / mzero ; /* sigma-tilde squared */
129
130 /* set up to compute f(s) */
131
132 zrat = zstar*zstar / zsigt ;
133 fac = ( zrat * nzero ) / ( SQRT_2PI * mzero ) ;
134 ss = zstar ; /* initial guess for s = zstar/sigma */
135
136 /* Newton's method [almost] */
137
138 for( ii=0 ; ii < 5 ; ii++ ){
139 pp = PHI(ss) ; /* Phi(ss) \approx 1 */
140 ee = exp(-0.5*ss*ss) ;
141 ff = ss*ss - (fac/pp) * ss * ee - zrat ; /* f(s) */
142 fp = 2.0*ss + (fac/pp) * ee * (ss*ss-1.0) ; /* f'(s) */
143 ds = ff / fp ; /* Newton step */
144 ss -= ds ; /* update */
145 }
146
147 sigma = zstar / ss ; /* estimate of sigma */
148 /* from upper tail data */
149
150 if( sigma <= 1.0 ){ /* the boring case */
151 *up = 0.0 ;
152 } else {
153
154 /* compute the log-likelihood difference */
155
156 uval = nzero * log( PHI(ss)/(1.0-pstar) )
157 - mzero * ( log(sigma) + 0.5 * zsigt * (1.0/(sigma*sigma)-1.0) ) ;
158
159 *up = uval ;
160 }
161 }
162
163 /****************************************/
164 /*** Compute negative tail unusuality ***/
165
166 fac = 1.0 / zsig ;
167
168 /* Find values <= -zstar, compute sum of squares */
169
170 rcut = tanh( zmid - zsig * zstar ) ; /* (atanh(zz)-zmid)/zsig <= -zstar */
171 zsigt = 0.0 ;
172 for( mzero=ii=0 ; ii < nr ; ii++ ){
173 if( zz[ii] <= rcut ){
174 ff = fac * ( atanh(zz[ii]) - zmid ) ; /* normalized Fisher */
175 zsigt += ff*ff ; /* sum of squares */
176 mzero++ ; /* how many we get */
177 }
178 }
179 nzero = nr - mzero ;
180
181 /* if we don't have decent data, output is 0 */
182
183 if( nzero < 2 || mzero < MAX(1.0,pstar*nr) ){ /* too weird for words */
184 *um = 0.0 ;
185 } else { /* have decent data here */
186 zsigt = zsigt / mzero ; /* sigma-tilde squared */
187
188 /* set up to compute f(s) */
189
190 zrat = zstar*zstar / zsigt ;
191 fac = ( zrat * nzero ) / ( SQRT_2PI * mzero ) ;
192 ss = zstar ; /* initial guess for s = zstar/sigma */
193
194 /* Newton's method [almost] */
195
196 for( ii=0 ; ii < 5 ; ii++ ){
197 pp = PHI(ss) ; /* Phi(ss) \approx 1 */
198 ee = exp(-0.5*ss*ss) ;
199 ff = ss*ss - (fac/pp) * ss * ee - zrat ; /* f(s) */
200 fp = 2.0*ss + (fac/pp) * ee * (ss*ss-1.0) ; /* f'(s) */
201 ds = ff / fp ; /* Newton step */
202 ss -= ds ; /* update */
203 }
204
205 sigma = zstar / ss ; /* estimate of sigma */
206 /* from upper tail data */
207
208 if( sigma <= 1.0 ){ /* the boring case */
209 *um = 0.0 ;
210 } else {
211
212 /* compute the log-likelihood difference */
213
214 uval = nzero * log( PHI(ss)/(1.0-pstar) )
215 - mzero * ( log(sigma) + 0.5 * zsigt * (1.0/(sigma*sigma)-1.0) ) ;
216
217 *um = uval ;
218 }
219 }
220
221 /*-- done! --*/
222
223 return ;
224 }
225
226 /***************************************************************************/
227
228 static int nt = 0 ; /* length of vectors [bitvec and float] */
229 static int nfv = 0 ; /* number of float vectors */
230 static int nlev = 2 ; /* default number of levels in bitvecs */
231
232 typedef struct {
233 byte * bv ; /* bitvector [nt] */
234 float * dp ; /* dot products [nfv] */
235
236 float up , um , ubest ; /* unusualities */
237 int nlev ; /* # levels */
238 } bitvec ;
239
240 #define bv_free(b) \
241 do{ if((b)!=NULL){free((b)->bv);free((b)->dp);free((b));} }while(0)
242
243 typedef struct {
244 int num , nall ;
245 bitvec ** bvarr ;
246 } bvarr ;
247
248 static bvarr * bvar = NULL ;
249 static float ** fvar = NULL ; /* nfv of these */
250
251 #define BITVEC_IN_BVARR(name,nn) ((name)->bvarr[(nn)])
252 #define BVARR_SUB BITVEC_IN_BVARR
253 #define BVARR_COUNT(name) ((name)->num)
254
255 #define INC_BVARR 32
256
257 #define INIT_BVARR(name) \
258 do{ int iq ; (name) = (bvarr *) malloc(sizeof(bvarr)) ; \
259 (name)->num = 0 ; (name)->nall = INC_BVARR ; \
260 (name)->bvarr = (bitvec **)malloc(sizeof(bitvec *)*INC_BVARR) ; \
261 for( iq=(name)->num ; iq < (name)->nall ; iq++ ) (name)->bvarr[iq] = NULL ; \
262 break ; } while(0)
263
264 #define ADDTO_BVARR(name,imm) \
265 do{ int nn , iq ; \
266 if( (name)->num == (name)->nall ){ \
267 nn = (name)->nall = 1.1*(name)->nall + INC_BVARR ; \
268 (name)->bvarr = realloc( (name)->bvarr,sizeof(bitvec *)*nn ); \
269 for( iq=(name)->num ; iq < (name)->nall ; iq++ ) (name)->bvarr[iq] = NULL ; } \
270 nn = (name)->num ; ((name)->num)++ ; \
271 (name)->bvarr[nn] = (imm) ; break ; } while(0)
272
273 #define FREE_BVARR(name) \
274 do{ if( (name) != NULL ){ \
275 free((name)->bvarr); free((name)); (name) = NULL; } break; } while(0)
276
277 #define DESTROY_BVARR(name) \
278 do{ int nn ; \
279 if( (name) != NULL ){ \
280 for( nn=0 ; nn < (name)->num ; nn++ ) bv_free((name)->bvarr[nn]) ; \
281 free((name)->bvarr); free((name)); (name) = NULL; } break; } while(0)
282
283 #define TRUNCATE_BVARR(name,qq) \
284 do{ int nn ; \
285 if( (name) != NULL && qq < (name)->num ){ \
286 for( nn=qq ; nn < (name)->num ; nn++ ) bv_free((name)->bvarr[nn]); \
287 (name)->num = qq ; \
288 } } while(0)
289
290 /***************************************************************************/
291
equal_bitvector_piece(bitvec * b,bitvec * c,int aa,int bb)292 int equal_bitvector_piece( bitvec * b , bitvec * c , int aa , int bb )
293 {
294 int ii ; byte * bv=b->bv , * cv=c->bv ;
295
296 if( aa < 0 ) aa = 0 ;
297 if( bb >= nt ) bb = nt-1 ;
298 for( ii=aa ; ii <= bb ; ii++ ) if( bv[ii] != cv[ii] ) return 0 ;
299 return 1;
300 }
301
equal_bitvector(bitvec * b,bitvec * c)302 int equal_bitvector( bitvec * b , bitvec * c )
303 {
304 return equal_bitvector_piece( b , c , 0 , nt-1 ) ;
305 }
306
307 /*--------------------------------------------------------------------------*/
308
randomize_bitvector_piece(bitvec * b,int aa,int bb)309 void randomize_bitvector_piece( bitvec * b , int aa , int bb )
310 {
311 int ii ; byte * bv=b->bv ;
312
313 if( aa < 0 ) aa = 0 ;
314 if( bb >= nt ) bb = nt-1 ;
315 for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = (byte)( b->nlev*drand48() ) ;
316 return ;
317 }
318
randomize_bitvector(bitvec * b)319 void randomize_bitvector( bitvec * b )
320 {
321 randomize_bitvector_piece( b , 0 , nt-1 ) ;
322 return ;
323 }
324
325 /*--------------------------------------------------------------------------*/
326
zero_bitvector_piece(bitvec * b,int aa,int bb)327 void zero_bitvector_piece( bitvec * b , int aa , int bb )
328 {
329 int ii ; byte * bv=b->bv ;
330
331 if( aa < 0 ) aa = 0 ;
332 if( bb >= nt ) bb = nt-1 ;
333 for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = 0 ;
334 return ;
335 }
336
zero_bitvector(bitvec * b)337 void zero_bitvector( bitvec * b )
338 {
339 zero_bitvector_piece( b , 0 , nt-1 ) ;
340 return ;
341 }
342
343 /*--------------------------------------------------------------------------*/
344
fix_bitvector_piece(bitvec * b,int aa,int bb,int val)345 void fix_bitvector_piece( bitvec * b , int aa , int bb , int val )
346 {
347 int ii ; byte * bv=b->bv , vv=(byte)val ;
348
349 if( aa < 0 ) aa = 0 ;
350 if( bb >= nt ) bb = nt-1 ;
351 for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = vv ;
352 return ;
353 }
354
fix_bitvector(bitvec * b,int val)355 void fix_bitvector( bitvec * b , int val )
356 {
357 fix_bitvector_piece( b , 0 , nt-1 , val ) ;
358 return ;
359 }
360
361 /*--------------------------------------------------------------------------*/
362
invert_bitvector_piece(bitvec * b,int aa,int bb)363 void invert_bitvector_piece( bitvec * b , int aa , int bb )
364 {
365 int ii,nl=b->nlev-1 ; byte * bv=b->bv ;
366
367 if( aa < 0 ) aa = 0 ;
368 if( bb >= nt ) bb = nt-1 ;
369 for( ii=aa ; ii <= bb ; ii++ ) bv[ii] = nl - bv[ii] ;
370 return ;
371 }
372
invert_bitvector(bitvec * b)373 void invert_bitvector( bitvec * b )
374 {
375 invert_bitvector_piece( b , 0 , nt-1 ) ;
376 return ;
377 }
378
379 /*--------------------------------------------------------------------------*/
380
new_bitvector(void)381 bitvec * new_bitvector(void)
382 {
383 bitvec * b ;
384 b = (bitvec *) malloc(sizeof(bitvec) ) ;
385 b->bv = (byte *) malloc(sizeof(byte) *nt ) ;
386 b->dp = (float *) malloc(sizeof(float)*nfv) ;
387
388 b->nlev = nlev ;
389 return b ;
390 }
391
copy_bitvector(bitvec * b)392 bitvec * copy_bitvector( bitvec * b )
393 {
394 int ii ;
395 bitvec * c = new_bitvector() ;
396 memcpy( c->bv , b->bv , sizeof(byte)*nt ) ;
397 #if 0
398 memcpy( c->dp , b->dp , sizeof(float)*nfv ) ;
399 c->up = b->up ; c->um = b->um ; c->ubest = b->ubest ;
400 #endif
401 c->nlev = b->nlev ;
402 return c ;
403 }
404
405 /*--------------------------------------------------------------------------*/
406
count_bitvector(bitvec * b)407 int count_bitvector( bitvec * b )
408 {
409 int ii,ss ;
410 byte * bv = b->bv ;
411 for( ii=ss=0 ; ii < nt ; ii++ ) if( bv[ii] ) ss++ ;
412 return ss ;
413 }
414
415 /*--------------------------------------------------------------------------*/
416
417 #define LINEAR_DETREND
418
normalize_floatvector(float * far)419 void normalize_floatvector( float * far )
420 {
421 int ii ;
422 float ff,gg ;
423
424 #ifdef LINEAR_DETREND
425 THD_linear_detrend( nt , far , NULL , NULL ) ; /* remove trend */
426 for( ff=0.0,ii=0 ; ii < nt ; ii++ ) ff += far[ii]*far[ii] ; /* and normalize */
427 if( ff <= 0.0 ) return ;
428 ff = 1.0 / sqrt(ff) ;
429 for( ii=0 ; ii < nt ; ii++ ) far[ii] *= ff ;
430 #else
431 for( ff=0.0,ii=0 ; ii < nt ; ii++ ) ff += far[ii] ;
432 ff /= nt ;
433 for( gg=0.0,ii=0 ; ii < nt ; ii++ ) gg += SQR( (far[ii]-ff) ) ;
434 if( gg <= 0.0 ) return ;
435 gg = 1.0 / sqrt(gg) ;
436 for( ii=0 ; ii < nt ; ii++ ) far[ii] = (far[ii]-ff)*gg ;
437 #endif
438 return ;
439 }
440
441 /*--------------------------------------------------------------------------*/
442
corr_floatbit(float * far,bitvec * b)443 float corr_floatbit( float * far , bitvec * b )
444 {
445 int ii , ns ;
446 float ss , bb,bq ;
447 byte * bar=b->bv ;
448
449 if( b->nlev == 2 ){ /* binary case */
450 for( ss=0.0,ns=ii=0 ; ii < nt ; ii++ )
451 if( bar[ii] ){ ns++ ; ss += far[ii] ; }
452 if( ns == 0 || ns == nt ) return 0.0 ;
453 ss *= sqrt( ((float) nt) / (float)(ns*(nt-ns)) ) ;
454
455 } else { /* multilevel case */
456 for( ss=bb=bq=0.0,ii=0 ; ii < nt ; ii++ ){
457 ss += bar[ii]*far[ii] ;
458 bb += bar[ii] ; bq += bar[ii]*bar[ii] ;
459 }
460 bq -= bb*bb/nt ; if( bq <= 0.0 ) return 0.0 ;
461 ss /= sqrt(bq) ;
462 }
463
464 return ss ;
465 }
466
467 /*--------------------------------------------------------------------------*/
468
evaluate_bitvec(bitvec * bim)469 void evaluate_bitvec( bitvec * bim )
470 {
471 int jj ;
472
473 for( jj=0 ; jj < nfv ; jj++ )
474 bim->dp[jj] = corr_floatbit( fvar[jj] , bim ) ;
475
476 unusuality( nfv , bim->dp , &(bim->up) , &(bim->um) ) ;
477
478 if( bim->up < bim->um ){
479 float tt ;
480 invert_bitvector( bim ) ;
481 tt = bim->um ; bim->um = bim->up ; bim->up = tt ;
482 }
483
484 bim->ubest = bim->up ;
485 return ;
486 }
487
488 /*--------------------------------------------------------------------------*/
489
490 #define DERR(s) fprintf(stderr,"** %s\n",(s))
491
init_floatvector_array(char * dname,char * mname)492 void init_floatvector_array( char * dname , char * mname )
493 {
494 THD_3dim_dataset * dset ;
495 byte * mask = NULL ;
496 int ii,jj , nvox ;
497 MRI_IMAGE * im ;
498
499 dset = THD_open_dataset( dname ) ;
500 if( dset == NULL ){ DERR("can't open dataset"); return; }
501 nt = DSET_NVALS(dset) ;
502 if( nt < 20 ){ DSET_delete(dset); DERR("dataset too short"); return; }
503 DSET_load(dset) ;
504 if( !DSET_LOADED(dset) ){ DSET_delete(dset); DERR("can't load dataset"); return; }
505 nvox = DSET_NVOX(dset) ;
506
507 if( mname != NULL ){
508 THD_3dim_dataset * mset ;
509 mset = THD_open_dataset( mname ) ;
510 if( mset == NULL ){ DSET_delete(dset); DERR("can't open mask"); return; }
511 if( DSET_NVOX(mset) != nvox ){
512 DSET_delete(mset); DSET_delete(dset); DERR("mask size mismatch"); return;
513 }
514 mask = THD_makemask( mset , 0 , 1.0,0.0 ) ;
515 DSET_delete(mset) ;
516 if( mask == NULL ){ DSET_delete(dset); DERR("mask is empty"); return; }
517 nfv = THD_countmask( nvox , mask ) ;
518 if( nfv < nt ){ DSET_delete(dset); DERR("mask is too small"); return; }
519 } else {
520 nfv = nvox ;
521 }
522
523 fvar = (float **) malloc(sizeof(float *)*nfv) ;
524 for( jj=ii=0 ; ii < nvox ; ii++ ){
525 if( mask != NULL && mask[ii] == 0 ) continue ; /* skip */
526
527 im = THD_extract_series( ii , dset , 0 ) ;
528 fvar[jj] = MRI_FLOAT_PTR(im) ;
529 normalize_floatvector( fvar[jj] ) ;
530 mri_fix_data_pointer(NULL,im) ; mri_free(im) ; jj++ ;
531 }
532
533 if( mask != NULL ) free(mask) ;
534 DSET_delete(dset) ; return ;
535 }
536
537 /*--------------------------------------------------------------------------*/
538
init_bitvector_array(int nbv)539 void init_bitvector_array( int nbv )
540 {
541 bitvec * bim ;
542 int ii , jj ;
543 byte * bar ;
544
545 INIT_BVARR(bvar) ;
546
547 for( ii=0 ; ii < nbv ; ii++ ){
548 bim = new_bitvector() ;
549 randomize_bitvector( bim ) ;
550 evaluate_bitvec( bim ) ;
551 ADDTO_BVARR(bvar,bim) ;
552 }
553
554 return ;
555 }
556
557 /*--------------------------------------------------------------------------*/
558
559 #define IRAN(j) (lrand48() % (j))
560
561 #define PROMO_MAX 4
562
563 static int promo_ok = 0 ;
564
evolve_bitvector_array(void)565 void evolve_bitvector_array(void)
566 {
567 static int nrvec=-1 ;
568 static float * rvec=NULL ;
569
570 int ii , nbv,nbv1 , aa,bb,vv , qbv , kup ;
571 bitvec * bim , * qim ;
572 float aup,aum ;
573
574 /* promote a few to a higher plane of being */
575
576 nbv1=nbv = BVARR_COUNT(bvar) ;
577
578 if( promo_ok ){
579 for( aa=ii=0 ; ii < nbv ; ii++ )
580 if( BVARR_SUB(bvar,ii)->nlev > nlev ) aa++ ;
581
582 if( aa < nbv/4 ){
583 for( kup=ii=0 ; ii < nbv && kup < PROMO_MAX ; ii++ ){
584
585 bim = BVARR_SUB(bvar,ii) ;
586 if( bim->nlev > nlev ) continue ; /* skip */
587
588 qim = copy_bitvector(bim) ;
589 qim->nlev *= 2 ;
590 for( aa=0 ; aa < nt ; aa++ )
591 if( qim->bv[aa] < nlev-1 ) qim->bv[aa] *= 2 ;
592 else qim->bv[aa] = 2*nlev-1 ;
593 evaluate_bitvec( qim ) ;
594 ADDTO_BVARR(bvar,qim) ;
595 kup++ ;
596 }
597 fprintf(stderr,"%d PROMO up\n",kup) ;
598 } else if( aa > 3*nbv/4 ){
599 for( kup=ii=0 ; ii < nbv && kup < PROMO_MAX ; ii++ ){
600
601 bim = BVARR_SUB(bvar,ii) ;
602 if( bim->nlev == nlev ) continue ; /* skip */
603
604 qim = copy_bitvector(bim) ;
605 qim->nlev /= 2 ;
606 for( aa=0 ; aa < nt ; aa++ ) qim->bv[aa] /= 2 ;
607 evaluate_bitvec( qim ) ;
608 ADDTO_BVARR(bvar,qim) ;
609 kup++ ;
610 }
611 fprintf(stderr,"%d PROMO down\n",kup) ;
612 }
613 }
614 /* create mutants */
615
616 qim = copy_bitvector(BVARR_SUB(bvar,0)) ; /* add copy of first one */
617 evaluate_bitvec( qim ) ;
618 ADDTO_BVARR(bvar,qim) ;
619
620 nbv = BVARR_COUNT(bvar) ;
621
622 for( ii=0 ; ii < nbv ; ii++ ){
623
624 bim = BVARR_SUB(bvar,ii) ;
625
626 aa = IRAN(nt) ; bb = aa + IRAN(5) ; if( bb >= nt ) bb = nt-1 ;
627
628 qim = copy_bitvector(bim) ;
629 zero_bitvector_piece( qim , aa , bb ) ;
630 if( equal_bitvector_piece(bim,qim,aa,bb) ){
631 bv_free(qim) ;
632 } else {
633 evaluate_bitvec( qim ) ;
634 ADDTO_BVARR(bvar,qim) ;
635 }
636
637 vv = (bim->nlev == 2) ? 1 : 1+IRAN(bim->nlev-1) ;
638 qim = copy_bitvector(bim) ;
639 fix_bitvector_piece( qim , aa , bb , vv ) ;
640 if( equal_bitvector_piece(bim,qim,aa,bb) ){
641 bv_free(qim) ;
642 } else {
643 evaluate_bitvec( qim ) ;
644 ADDTO_BVARR(bvar,qim) ;
645 }
646
647 qim = copy_bitvector(bim) ;
648 randomize_bitvector_piece( qim , aa , bb ) ;
649 if( equal_bitvector_piece(bim,qim,aa,bb) ){
650 bv_free(qim) ;
651 } else {
652 evaluate_bitvec( qim ) ;
653 ADDTO_BVARR(bvar,qim) ;
654 }
655
656 qim = copy_bitvector(bim) ;
657 invert_bitvector_piece( qim , aa , bb ) ;
658 if( equal_bitvector_piece(bim,qim,aa,bb) ){
659 bv_free(qim) ;
660 } else {
661 evaluate_bitvec( qim ) ;
662 ADDTO_BVARR(bvar,qim) ;
663 }
664 }
665
666 /* sort everybody */
667
668 qbv = BVARR_COUNT(bvar) ;
669 if( nrvec < qbv ){
670 if( rvec != NULL ) free(rvec) ;
671 rvec = (float *) malloc(sizeof(float)*qbv) ;
672 nrvec = qbv ;
673 }
674 for( ii=0 ; ii < qbv ; ii++ )
675 rvec[ii] = - BVARR_SUB(bvar,ii)->ubest ;
676
677 qsort_floatstuff( qbv , rvec , (void **) bvar->bvarr ) ;
678
679 TRUNCATE_BVARR( bvar , nbv1 ) ;
680 return ;
681 }
682
683 /*--------------------------------------------------------------------------*/
684
main(int argc,char * argv[])685 int main( int argc , char * argv[] )
686 {
687 int ii , nv , nite=0 , neq=0 , nopt , nbv=64 ;
688 float fold , fnew ;
689 char * mname=NULL , * dname ;
690 bitvec * bim ;
691
692 if( argc < 2 ){printf("Usage: unu [-mask mset] [-lev n] [-nbv n] dset > bname.1D\n");exit(0);}
693
694 nopt = 1 ;
695 while( nopt < argc && argv[nopt][0] == '-' ){
696
697 if( strcmp(argv[nopt],"-mask") == 0 ){
698 mname = argv[++nopt] ;
699 nopt++ ; continue ;
700 }
701
702 if( strcmp(argv[nopt],"-lev") == 0 ){
703 nlev = (int)strtod(argv[++nopt],NULL) ;
704 if( nlev < 2 || nlev > 8 ){ DERR("bad -nlev"); exit(1); }
705 nopt++ ; continue ;
706 }
707
708 if( strcmp(argv[nopt],"-nbv") == 0 ){
709 nbv = (int)strtod(argv[++nopt],NULL) ;
710 if( nbv < 8 || nbv > 999 ){ DERR("bad -nbv"); exit(1); }
711 nopt++ ; continue ;
712 }
713
714 fprintf(stderr,"** Illegal option %s\n",argv[nopt]); exit(1);
715 }
716 if( nopt >= argc ){fprintf(stderr,"** No dataset!?\n");exit(1);}
717
718 dname = argv[nopt] ;
719
720 init_floatvector_array( dname , mname ) ;
721 if( fvar == NULL ){
722 fprintf(stderr,"** Couldn't init floatvector!?\n") ; exit(1) ;
723 } else {
724 fprintf(stderr,"nt=%d nfv=%d\n",nt,nfv) ;
725 }
726
727 srand48((long)time(NULL)) ;
728
729 init_bitvector_array( nbv ) ;
730 fold = BVARR_SUB(bvar,0)->ubest ;
731 fprintf(stderr,"fold = %7.4f\n",fold) ;
732
733 while(1){
734 evolve_bitvector_array() ;
735 nv = BVARR_COUNT(bvar) ;
736 nite++ ;
737 #if 1
738 fprintf(stderr,"---nite=%d\n",nite) ;
739 for( nopt=ii=0 ; ii < nv ; ii++ ){
740 fprintf(stderr," %7.4f[%d]",BVARR_SUB(bvar,ii)->ubest,BVARR_SUB(bvar,ii)->nlev) ;
741 if( BVARR_SUB(bvar,ii)->nlev > nlev ) nopt++ ;
742 }
743 fprintf(stderr," *%d\n",nopt) ;
744 #endif
745
746 fnew = fabs(BVARR_SUB(bvar,0)->ubest) ;
747 if( fnew <= fold ){
748 neq++ ;
749 if( neq == 8 && promo_ok ) break ;
750 if( neq == 8 && !promo_ok ){ promo_ok = 1 ; neq = 0 ; }
751 } else {
752 neq = 0 ;
753 fold = fnew ;
754 fprintf(stderr,"%d: %7.4f\n",nite,fold) ;
755 }
756 }
757
758 bim = BVARR_SUB(bvar,0) ;
759 for( ii=0 ; ii < nt ; ii++ )
760 printf(" %d\n",(int)bim->bv[ii]) ;
761
762 exit(0) ;
763 }
764