1 #include "mrilib.h"
2
3 /*==============================================================================*/
4 /*========== The following functions were moved from afni_fimfunc.c -===========*/
5 /*===========================================================================*/
6
7 /*---------------------------------------------------------------------------*/
8 /*! Rank-order a float array, with ties getting the average rank.
9 The output overwrites the input.
10 [27 Jun 2010: modified to create/destroy workspace on each call]
11 *//*-------------------------------------------------------------------------*/
12
rank_order_float(int n,float * a)13 void rank_order_float( int n , float *a )
14 {
15 register int ii , ns , n1 , ib ;
16 int *b ; /* workspaces */
17 float *c ;
18 float cs ;
19
20 /*- handle special cases -*/
21
22 if( a == NULL || n < 1 ) return ; /* meaningless input */
23 if( n == 1 ){ a[0] = 0.0f ; return ; } /* only one point!? */
24
25 /*- make workspaces -*/
26
27 b = (int *) malloc(sizeof(int )*n) ;
28 c = (float *) malloc(sizeof(float)*n) ;
29
30 for( ii=0 ; ii < n ; ii++ ) c[ii] = b[ii] = ii ;
31
32 /*- sort input, carrying b along -*/
33
34 qsort_floatint( n , a , b ) ; /* see cs_sort_fi.c */
35
36 /* compute ranks into c[] */
37
38 n1 = n-1 ;
39 for( ii=0 ; ii < n1 ; ii++ ){
40 if( a[ii] == a[ii+1] ){ /* handle ties */
41 cs = 2*ii+1 ; ns = 2 ; ib = ii ; ii++ ;
42 while( ii < n1 && a[ii] == a[ii+1] ){ ii++ ; ns++ ; cs += ii ; }
43 for( cs/=ns ; ib <= ii ; ib++ ) c[ib] = cs ;
44 }
45 }
46
47 for( ii=0 ; ii < n ; ii++ ) a[b[ii]] = c[ii] ;
48
49 free(c) ; free(b) ; return ;
50 }
51
52 /*---------------------------------------------------------------------------*/
53 /*! Rank orders a[], subtracts the mean rank, and returns the sum-of-squares.
54 -----------------------------------------------------------------------------*/
55
spearman_rank_prepare(int n,float * a)56 float spearman_rank_prepare( int n , float *a )
57 {
58 register int ii ;
59 register float rb , rs ;
60
61 rank_order_float( n , a ) ;
62
63 rb = 0.5f*(n-1) ; rs=0.0f ;
64 for( ii=0 ; ii < n ; ii++ ){
65 a[ii] -= rb ; /* remove mean rank */
66 rs += a[ii]*a[ii] ; /* sum squares */
67 }
68
69 return rs ;
70 }
71
72 /*---------------------------------------------------------------------------*/
73 /*! Prepare for quadrant correlation with a[].
74 -----------------------------------------------------------------------------*/
75
quadrant_corr_prepare(int n,float * a)76 float quadrant_corr_prepare( int n , float *a )
77 {
78 register int ii ;
79 register float rb , rs ;
80
81 rank_order_float( n , a ) ;
82
83 rb = 0.5f*(n-1) ; rs=0.0f ;
84 for( ii=0 ; ii < n ; ii++ ){
85 a[ii] = (a[ii] > rb) ? 1.0
86 : (a[ii] < rb) ? -1.0f : 0.0f ;
87 rs += a[ii]*a[ii] ;
88 }
89
90 return rs ;
91 }
92
93 /*-----------------------------------------------------------------------------*/
94 /*! To Spearman (rank-order) correlate x[] with r[], first do
95 rv = spearman_rank_prepare(n,r) ;
96 then
97 corr = spearman_rank_corr(n,x,rv,r) ;
98 Note that these 2 routines are destructive (r and x are replaced by ranks).
99 -------------------------------------------------------------------------------*/
100
spearman_rank_corr(int n,float * x,float rv,float * r)101 float spearman_rank_corr( int n , float *x , float rv , float *r )
102 {
103 register int ii ;
104 register float ss ; float xv ;
105
106 xv = spearman_rank_prepare( n , x ) ; if( xv <= 0.0f ) return 0.0f ;
107
108 for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii] * r[ii] ;
109
110 return ( ss/sqrtf(rv*xv) ) ;
111 }
112
113 /*------------------------------------------------------------------------------*/
114 /*! To do quadrant correlation of x[] with r[], first do
115 rv = quadrant_corr_prepare(n,r) ;
116 then
117 corr = quadrant_corr(n,x,rv,r) ;
118 Note that these 2 routines are destructive (r and x are modified).
119 -------------------------------------------------------------------------------*/
120
quadrant_corr(int n,float * x,float rv,float * r)121 float quadrant_corr( int n , float *x , float rv , float *r )
122 {
123 register int ii ;
124 register float ss ; float xv ;
125
126 xv = quadrant_corr_prepare( n , x ) ; if( xv <= 0.0f ) return 0.0f ;
127
128 for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii] * r[ii] ;
129
130 return ( ss/sqrtf(rv*xv) ) ;
131 }
132
133 /*------------------------------------------------------------------------------*/
134
135 static int num_quantile = 9 ;
THD_quantile_corr_setup(int nq)136 void THD_quantile_corr_setup( int nq )
137 {
138 if( nq > 1 && nq < 100 ) num_quantile = nq ;
139 }
140
quantile_prepare(int n,float * a)141 float quantile_prepare( int n , float *a )
142 {
143 int ii ;
144 float rb , rs , jf ;
145
146 jf = 0.001f + 1.00001f * (n-0.5f) / (float)num_quantile ;
147 if( jf <= 2.0f ) return spearman_rank_prepare(n,a) ;
148 jf = 1.0f / jf ;
149
150 rank_order_float(n,a) ; /* convert to ranks */
151
152 for( rb=0.0f,ii=0 ; ii < n ; ii++ ){
153 a[ii] = (int)( (a[ii]+0.333f)*jf ) ; rb += a[ii] ;
154 }
155 rb /= n ;
156 for( rs=0.0f,ii=0 ; ii < n ; ii++ ){
157 a[ii] -= rb ; rs += a[ii]*a[ii] ;
158 }
159
160 return rs ;
161 }
162
THD_quantile_corr(int n,float * x,float * y)163 float THD_quantile_corr( int n , float *x , float *y )
164 {
165 float xv,yv,ss ; int ii ;
166
167 if( n < 2 ) return 0.0f ;
168
169 xv = quantile_prepare(n,x) ; if( xv <= 0.0f ) return 0.0f ;
170 yv = quantile_prepare(n,y) ; if( yv <= 0.0f ) return 0.0f ;
171
172 for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii] * y[ii] ;
173
174 return ( ss/sqrtf(yv*xv) ) ;
175 }
176
quantile_corr(int n,float * x,float rv,float * r)177 float quantile_corr( int n , float *x , float rv , float *r )
178 {
179 register int ii ;
180 register float ss ; float xv ;
181
182 xv = quantile_prepare( n , x ) ; if( xv <= 0.0f ) return 0.0f ;
183
184 for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii] * r[ii] ;
185
186 return ( ss/sqrtf(rv*xv) ) ;
187 }
188
189 /*---------------------------------------------------------------------------*/
190
191 #if 1
192 #undef ttt_bot
193 #undef ttt_top
194 #define ttt_bot 0.3333333f
195 #define ttt_top 0.6666667f
196
tictactoe_set_thresh(float bb,float tt)197 void tictactoe_set_thresh( float bb , float tt ){ return; }
198
199 #else
200
201 static float ttt_bot = 0.3333333f ;
202 static float ttt_top = 0.6666667f ;
203
tictactoe_set_thresh(float bb,float tt)204 void tictactoe_set_thresh( float bb , float tt )
205 {
206 if( bb >= 0.0f && bb < tt && tt <= 1.0f ){ ttt_bot = bb; ttt_top = tt; }
207 else { ttt_bot = 0.3333333f; ttt_top = 0.6666667f; }
208 }
209
210 #endif
211
212 /*---------------------------------------------------------------------------*/
213 /*! Prepare for tictactoe correlation with a[].
214 -----------------------------------------------------------------------------*/
215
tictactoe_corr_prepare(int n,float * a)216 float tictactoe_corr_prepare( int n , float *a )
217 {
218 register int ii ;
219 register float rb , rs , rt ;
220
221 rank_order_float( n , a ) ;
222
223 rb = ttt_bot * (n-1) ;
224 rt = ttt_top * (n-1) ;
225 rs = 0.0f ;
226 for( ii=0 ; ii < n ; ii++ ){
227 if( a[ii] > rt ){ a[ii] = 1.0f ; rs += 1.0f ; }
228 else if( a[ii] < rb ){ a[ii] = -1.0f ; rs += 1.0f ; }
229 else { a[ii] = 0.0f ; }
230 }
231
232 return rs ;
233 }
234
235 /*------------------------------------------------------------------------------*/
236 /*! To do tictactoe correlation of x[] with r[], first do
237 rv = tictactoe_corr_prepare(n,r) ;
238 then
239 corr = tictactoe_corr(n,x,rv,r) ;
240 Note that these 2 routines are destructive (r and x are modified).
241 -------------------------------------------------------------------------------*/
242
tictactoe_corr(int n,float * x,float rv,float * r)243 float tictactoe_corr( int n , float *x , float rv , float *r )
244 {
245 register int ii ;
246 register float ss ; float xv ;
247
248 xv = tictactoe_corr_prepare( n , x ) ; if( xv <= 0.0f ) return 0.0f ;
249
250 for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii] * r[ii] ;
251
252 return ( ss/sqrtf(rv*xv) ) ;
253 }
254
255 /*------------------------------------------------------------------------------*/
256
THD_tictactoe_corr(int n,float * x,float * y)257 float THD_tictactoe_corr( int n , float *x , float *y )
258 {
259 float yv , xv ; register float ss ; register int ii ;
260 if( n < 3 ) return 0.0f ;
261 xv = tictactoe_corr_prepare(n,x) ; if( xv <= 0.0f ) return 0.0f ;
262 yv = tictactoe_corr_prepare(n,y) ; if( yv <= 0.0f ) return 0.0f ;
263 for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii]*y[ii] ;
264 return( ss / sqrtf(xv*yv) ) ;
265 }
266
267 /*=============================================================================
268 Compute correlations, destructively (i.e., mangling the input arrays)
269 ===============================================================================*/
270
271 /*--------------------------------------------------------------------------*/
272 /*! Spearman rank-order correlation of x[] and y[] (x and y are modified). */
273
THD_spearman_corr(int n,float * x,float * y)274 float THD_spearman_corr( int n , float *x , float *y )
275 {
276 float xv ;
277 if( n < 2 ) return 0.0f ;
278 xv = spearman_rank_prepare(n,x) ;
279 if( xv <= 0.0f ) return 0.0f ;
280 return spearman_rank_corr( n,y,xv,x ) ;
281 }
282
283 /*--------------------------------------------------------------------------*/
284 /*! Spearman rank-order correlation of x[] and y[] (nondestructive). */
285
THD_spearman_corr_nd(int n,float * x,float * y)286 float THD_spearman_corr_nd( int n , float *x , float *y )
287 {
288 float *qx, *qy , cv=0.0f ;
289 qx = (float *)malloc(sizeof(float)*n); memcpy(qx,x,sizeof(float)*n);
290 qy = (float *)malloc(sizeof(float)*n); memcpy(qy,y,sizeof(float)*n);
291 cv = THD_spearman_corr(n,qx,qy) ;
292 free((void *)qy); free((void *)qx);
293 return cv ;
294 }
295
THD_spearman_corr_dble(int n,double * x,double * y)296 double THD_spearman_corr_dble( int n , double *x , double *y )
297 {
298 float *qx, *qy , cv=0.0f ; int ii ;
299 qx = (float *)malloc(sizeof(float)*n);
300 qy = (float *)malloc(sizeof(float)*n);
301 for( ii=0 ; ii < n ; ii++ ){ qx[ii] = x[ii] ; qy[ii] = y[ii] ; }
302 cv = THD_spearman_corr(n,qx,qy) ;
303 free((void *)qy); free((void *)qx);
304 return (double)cv ;
305 }
306
307 /*--------------------------------------------------------------------------*/
308 /*! Kendall Tau_b (x and y are modified) */
309
THD_ktaub_corr(int n,float * x,float * y)310 float THD_ktaub_corr( int n , float *x , float *y )
311 {
312 if( n < 2 ) return 0.0f ;
313 qsort_floatfloat( n , x , y ) ; /* preliminary sorting of x, carrying y */
314 return kendallNlogN( x , y , n ) ; /* the actual work */
315 }
316
317 /*--------------------------------------------------------------*/
318 /*! Quadrant correlation of x[] and y[] (x and y are modified). */
319
THD_quadrant_corr(int n,float * x,float * y)320 float THD_quadrant_corr( int n , float *x , float *y )
321 {
322 float xv ;
323 if( n < 2 ) return 0.0f ;
324 xv = quadrant_corr_prepare(n,x) ;
325 if( xv <= 0.0f ) return 0.0f ;
326 return quadrant_corr( n,y,xv,x ) ;
327 }
328
329 /*--------------------------------------------------------------------------*/
330 /*! Quadrant rank-order correlation of x[] and y[] (nondestructive). */
331
332 #if 1
THD_quadrant_corr_nd(int n,float * x,float * y)333 float THD_quadrant_corr_nd( int n , float *x , float *y )
334 {
335 float *z ; register float xm,ym,qc ; register int ii ;
336
337 z = (float *)malloc(sizeof(float)*n) ;
338 memcpy( z , x , sizeof(float)*n ) ;
339 xm = qmed_float( n , z ) ;
340 z = (float *)malloc(sizeof(float)*n) ;
341 memcpy( z , y , sizeof(float)*n ) ;
342 ym = qmed_float( n , z ) ;
343 free(z) ;
344
345 qc = 0.0f ;
346 for( ii=0 ; ii < n ; ii++ )
347 qc += (x[ii] > xm) * (y[ii] > ym) ;
348 qc = (4.0f*qc) / n - 1.0f ;
349 if( qc < -1.0f ) qc = -1.0f; else if( qc > 1.0f ) qc = 1.0f;
350 qc = sinf(1.570796f*qc) ; /* adjust to normal model */
351 return qc ;
352 }
353 #else
354 /*--------------------------------------------------------------------------*/
THD_quadrant_corr_nd(int n,float * x,float * y)355 float THD_quadrant_corr_nd( int n , float *x , float *y )
356 {
357 float *qx, *qy , cv=0.0f ;
358 qx = (float *)malloc(sizeof(float)*n); memcpy(qx,x,sizeof(float)*n);
359 qy = (float *)malloc(sizeof(float)*n); memcpy(qy,y,sizeof(float)*n);
360 cv = THD_quadrant_corr(n,qx,qy) ;
361 free((void *)qy); free((void *)qx);
362 return cv ;
363 }
364 #endif
365
366 /*--------------------------------------------------------------------------*/
367 /*! Pearson correlation of x[] and y[] (x and y are NOT modified. */
368
THD_pearson_corr(int n,float * x,float * y)369 float THD_pearson_corr( int n, float *x , float *y )
370 {
371 float xv=0.0f , yv=0.0f , xy=0.0f , vv,ww ;
372 float xm=0.0f , ym=0.0f ;
373 register int ii ;
374
375 if( n < 2 ) return 0.0f ;
376 for( ii=0 ; ii < n ; ii++ ){ xm += x[ii] ; ym += y[ii] ; }
377 xm /= n ; ym /= n ;
378 for( ii=0 ; ii < n ; ii++ ){
379 vv = x[ii]-xm; ww = y[ii]-ym; xv += vv*vv; yv += ww*ww; xy += vv*ww;
380 }
381
382 if( xv <= 0.0f || yv <= 0.0f ) return 0.0f ;
383 return xy/sqrtf(xv*yv) ;
384 }
385
386 /*--------------------------------------------------------------------------*/
387 /*! Double Pearson correlation of x[] and y[] (x and y are NOT modified) */
388
THD_pearson_corrd(int n,double * x,double * y)389 double THD_pearson_corrd( int n, double *x , double *y )
390 {
391 double xv=0.0 , yv=0.0 , xy=0.0 , vv,ww ;
392 double xm=0.0 , ym=0.0 ;
393 register int ii ;
394
395 if( n < 2 ) return 0.0 ;
396 for( ii=0 ; ii < n ; ii++ ){ xm += x[ii] ; ym += y[ii] ; }
397 xm /= n ; ym /= n ;
398 for( ii=0 ; ii < n ; ii++ ){
399 vv = x[ii]-xm; ww = y[ii]-ym; xv += vv*vv; yv += ww*ww; xy += vv*ww;
400 }
401
402 if( xv <= 0.0 || yv <= 0.0 ) return 0.0 ;
403 return xy/sqrt(xv*yv) ;
404 }
405
406 /*-----------------------------------------------------------------------*/
407 /* Partial Correlation of A & B, adjusting for C. */
408
THD_pearson_partial_corr(int n,float * a,float * b,float * c)409 float THD_pearson_partial_corr( int n, float *a, float *b, float *c )
410 {
411 float AB, AC, BC;
412 //calculate correlation AB
413 AB = THD_pearson_corr(n, a, b);
414
415 //calculate correlation AC
416 AC = THD_pearson_corr(n, a, c);
417
418 //calculate correlation BC
419 BC = THD_pearson_corr(n, b, c);
420
421 return (AB - (AC*BC)) / sqrt( (1-(AC*AC)) * (1-(BC*BC)) );
422 }
423
424 /*-----------------------------------------------------------------------*/
425 /* Double Partial Correlation of A & B, adjusting for C. */
426
THD_pearson_partial_corrd(int n,double * a,double * b,double * c)427 double THD_pearson_partial_corrd( int n, double *a, double *b, double *c )
428 {
429 double AB, AC, BC;
430 //calculate correlation AB
431 AB = THD_pearson_corrd(n, a, b);
432
433 //calculate correlation AC
434 AC = THD_pearson_corrd(n, a, c);
435
436 //calculate correlation BC
437 BC = THD_pearson_corrd(n, b, c);
438
439 return (AB - (AC*BC)) / sqrt( (1-(AC*AC)) * (1-(BC*BC)) );
440 }
441
442 /*--------------------------------------------------------------------------*/
443 /*! Covariance of x[] and y[] (x and y are NOT modified). ZSS May 18 2012*/
444
THD_covariance(int n,float * x,float * y)445 float THD_covariance( int n, float *x , float *y )
446 {
447 float xy=0.0f , vv,ww ;
448 float xm=0.0f , ym=0.0f ;
449 register int ii ;
450
451 if( n < 2 ) return 0.0f ;
452 for( ii=0 ; ii < n ; ii++ ){ xm += x[ii] ; ym += y[ii] ; }
453 xm /= n ; ym /= n ;
454 for( ii=0 ; ii < n ; ii++ ){
455 xy += (x[ii]-xm)*(y[ii]-ym);
456 }
457
458 return xy/(float)(n-1) ;
459 }
460
461
462 /*-------------------------------------------------------------------------*/
463 /*! Returns a float_triple with (a,b,r) where
464 y = a*x + b
465 is the L2 regression result and r = Pearson correlation coeff.
466 For bootstrapping, ix[i] is the i-th index in x[] and y[] to use.
467 For non-bootstrapping, pass in ix==NULL.
468 *//*-----------------------------------------------------------------------*/
469
470 #undef IX
471 #define IX(i) ( ((ix) == NULL) ? (i) : ix[(i)] )
472
THD_pearson_indexed(int nix,int * ix,float * x,float * y)473 float_triple THD_pearson_indexed( int nix, int *ix, float *x, float *y )
474 {
475 float xbar=0,ybar=0, xq=0,yq=0,xyq=0, a=0,b=0,r=0;
476 int jj,ii; float_triple abr;
477
478 for( jj=0 ; jj < nix ; jj++ ){
479 ii = IX(jj); xbar += x[ii]; ybar += y[ii];
480 }
481 xbar /= nix ; ybar /= nix ;
482 for( jj=0 ; jj < nix ; jj++ ){
483 ii = IX(jj) ;
484 xq += (x[ii]-xbar)*(x[ii]-xbar) ;
485 yq += (y[ii]-ybar)*(y[ii]-ybar) ;
486 xyq += (x[ii]-xbar)*(y[ii]-ybar) ;
487 }
488 if( xq > 0.0f && yq > 0.0f ){
489 r = xyq/sqrtf(xq*yq); a = xyq/xq; b = (xq*ybar-xbar*xyq)/xq;
490 }
491 abr.a = a ; abr.b = b ; abr.c = r ; return abr ;
492 }
493
494 #undef IX
495
496 /*---------------------------------------------------------------------------*/
497 /* Correlates and also returns 2.5%..97.5% confidence interval, via bootstrap.
498 rrr = correlation coefficient, 2.5% level, 97.5% level [in that order]
499 aaa = regression 'a' coefficient, in y=ax+b [as .a .b .c ]
500 bbb = regression 'b' coefficient [components ]
501 *//*-------------------------------------------------------------------------*/
502
503 #undef NBOOT
504 #undef NB5
505
506 #define NBOOT 960
507 #define NB5 24 /* must be 2.5% of the above */
508
THD_pearson_corr_boot(int n,float * x,float * y,float_triple * rrr,float_triple * aaa,float_triple * bbb)509 void THD_pearson_corr_boot( int n , float *x , float *y ,
510 float_triple *rrr ,
511 float_triple *aaa ,
512 float_triple *bbb )
513 {
514 int ii,kk ; float aa,bb,rr ; float_triple abr ;
515 int *ix ; float ax[NBOOT] , bx[NBOOT] , rx[NBOOT] ;
516
517 if( n < 5 || x == NULL || y == NULL ) return ;
518 if( rrr == NULL && aaa == NULL && bbb == NULL ) return ;
519
520 /* compute standard results */
521
522 abr = THD_pearson_indexed( n , NULL , x, y ) ;
523 aa = abr.a ; bb = abr.b ; rr = abr.c ; /* non-bootstrapped answers */
524 ix = (int *)malloc(sizeof(int)*n) ;
525
526 /* compute bootstrap results */
527
528 for( kk=0 ; kk < NBOOT ; kk++ ){
529 for( ii=0 ; ii < n ; ii++ ) ix[ii] = lrand48() % n ;
530 abr = THD_pearson_indexed( n , ix , x, y ) ;
531 ax[kk] = abr.a ; bx[kk] = abr.b ; rx[kk] = abr.c ;
532 }
533 free(ix) ;
534
535 /* sort, then find 2.5% and 97.5% points, save into output structs */
536 /* [bootstrap confidence intervals are now bias-corrected!] */
537
538 if( rrr != NULL ){
539 float_triple qqq = THD_bootstrap_confinv( rr , 0.05f , NBOOT , rx ) ;
540 rrr->a = rr ; /* would use qqq.b for bias-corrected estimate */
541 rrr->b = qqq.a ; /* lower edge of confidence interval */
542 rrr->c = qqq.c ; /* upper edge */
543 }
544
545 if( aaa != NULL ){
546 float_triple qqq = THD_bootstrap_confinv( aa , 0.05f , NBOOT , ax ) ;
547 aaa->a = aa ;
548 aaa->b = qqq.a ;
549 aaa->c = qqq.c ;
550 }
551
552 if( bbb != NULL ){
553 float_triple qqq = THD_bootstrap_confinv( bb , 0.05f , NBOOT , bx ) ;
554 bbb->a = bb ;
555 bbb->b = qqq.a ;
556 bbb->c = qqq.c ;
557 }
558
559 return ;
560 }
561
562 /*---------------------------------------------------------------------------*/
563
THD_spearman_indexed(int nix,int * ix,float * x,float * y)564 float THD_spearman_indexed( int nix, int *ix, float *x, float *y )
565 {
566 float *xt , *yt ; float r ;
567
568 xt = (float *)malloc(sizeof(float)*nix) ;
569 yt = (float *)malloc(sizeof(float)*nix) ;
570 if( ix == NULL ){
571 memcpy(xt,x,sizeof(float)*nix) ;
572 memcpy(yt,y,sizeof(float)*nix) ;
573 } else {
574 int ii , jj ;
575 for( jj=0 ; jj < nix ; jj++ ){
576 ii = ix[jj] ; xt[jj] = x[ii] ; yt[jj] = y[ii] ;
577 }
578 }
579
580 r = THD_spearman_corr( nix , xt , yt ) ;
581 free(yt) ; free(xt) ; return r ;
582 }
583
584 /*------------------------------------------------------------------------------*/
585
THD_l1_fit_to_line(int n,float * x,float * y)586 float_pair THD_l1_fit_to_line( int n , float *x , float *y )
587 {
588 float_pair ab={0.0f,0.0f} ;
589 float *A[2] , coef[2] , val ;
590 int ii ;
591
592 if( n < 3 || x == NULL || y == NULL ) return ab ;
593
594 A[0] = x ;
595 A[1] = (float *)malloc(sizeof(float)*n) ;
596 for( ii=0 ; ii < n ; ii++ ) A[1][ii] = 1.0f ;
597
598 val = cl1_solve( n , 2 , y , A , coef , 0 ) ;
599 free(A[1]) ;
600 if( val >= 0.0f ){
601 ab.a = coef[0] ; ab.b = coef[1] ;
602 }
603 return ab ;
604 }
605
606 /*------------------------------------------------------------------------------*/
607
THD_spearman_corr_boot(int n,float * x,float * y,float_triple * rrr)608 void THD_spearman_corr_boot( int n , float *x , float *y , float_triple *rrr )
609 {
610 int ii,kk ; float rr ; float_triple qqq ;
611 int *ix ; float rx[NBOOT] ;
612
613 ENTRY("THD_spearman_corr_boot") ;
614
615 if( n < 5 || x == NULL || y == NULL ) EXRETURN ;
616 if( rrr == NULL ) EXRETURN ;
617
618 /* compute standard result */
619
620 rr = THD_spearman_indexed( n , NULL , x , y ) ; /* non-bootstrapped answer */
621 ix = (int *)malloc(sizeof(int)*n) ;
622
623 /* compute bootstrap results */
624
625 for( kk=0 ; kk < NBOOT ; kk++ ){
626 for( ii=0 ; ii < n ; ii++ ) ix[ii] = lrand48() % n ;
627 rx[kk] = THD_spearman_indexed( n , ix , x , y ) ;
628 }
629 free(ix) ;
630
631 /* sort, then find 2.5% and 97.5% points, save into output structs */
632 /* [bootstrap confidence intervals are now bias-corrected!] */
633
634 qqq = THD_bootstrap_confinv( rr , 0.05f , NBOOT , rx ) ;
635 rrr->a = rr ; /* would use qqq.b for bias-corrected estimate */
636 rrr->b = qqq.a ; /* lower edge of confidence interval */
637 rrr->c = qqq.c ; /* upper edge */
638
639 EXRETURN ;
640 }
641
642 /*----------------------------------------------------------------*/
643
THD_pearson_corr_wt(int n,float * x,float * y,float * wt)644 float THD_pearson_corr_wt( int n, float *x , float *y , float *wt )
645 {
646 float xv=0.0f , yv=0.0f , xy=0.0f , vv,ww ;
647 float xm=0.0f , ym=0.0f , ws=0.0f ;
648 int ii ;
649
650 if( wt == NULL ) return THD_pearson_corr(n,x,y) ;
651
652 for( ii=0 ; ii < n ; ii++ ){
653 xm += wt[ii]*x[ii]; ym += wt[ii]*y[ii]; ws += wt[ii];
654 }
655 xm /= ws ; ym /= ws ;
656 for( ii=0 ; ii < n ; ii++ ){
657 vv = x[ii]-xm ; ww = y[ii]-ym ;
658 xv += wt[ii]*vv*vv ; yv += wt[ii]*ww*ww ; xy += wt[ii]*vv*ww ;
659 }
660
661 #if 0
662 ININFO_message("THD_pearson_corr_wt: n=%d ws=%g xm=%g ym=%g xv=%g yv=%g xy=%g",
663 n,ws,xm,ym,xv,yv,xy ) ;
664 #endif
665
666 if( xv <= 0.0f || yv <= 0.0f ) return 0.0f ;
667 return xy/sqrtf(xv*yv) ;
668 }
669
670 /*----------------------------------------------------------------*/
671 /* compute distance between two arrays ZSS May 04 2012 */
THD_distance(int n,float * x,float * y,int abs)672 float THD_distance( int n, float *x , float *y, int abs)
673 {
674 float dp=0.0, a1, a2;
675 int ii, n1 = n-1;
676
677 for( ii=0 ; ii < n1 ; ii+=2 ) {
678 a1 = x[ii]-y[ii]; a2 = x[ii+1]-y[ii+1];
679 if (!abs) dp += (a1*a1+a2*a2);
680 else dp += (ABS(a1)+ABS(a2));
681 }
682 if( ii == n1 ) {
683 a1 = x[ii]-y[ii];
684 if (!abs) dp += a1*a1;
685 else dp += ABS(a1);
686 }
687
688 if (!abs) dp = sqrt(dp) ;
689
690 return (dp);
691 }
692
693 /*--------------------------------------------------------------------------*/
694 /*! Compute the rank-order correlation between 2 images [08 Mar 2006].
695 ----------------------------------------------------------------------------*/
696
mri_spearman_corr(MRI_IMAGE * im,MRI_IMAGE * jm)697 float mri_spearman_corr( MRI_IMAGE *im , MRI_IMAGE *jm )
698 {
699 float *far , *gar , cc ;
700 MRI_IMAGE *fim , *gim ;
701
702 if( im == NULL || jm == NULL || im->nvox != jm->nvox ) return 0.0f ;
703
704 fim = mri_to_float(im) ; far = mri_data_pointer(fim) ;
705 gim = mri_to_float(jm) ; gar = mri_data_pointer(gim) ;
706 cc = THD_spearman_corr( fim->nvox , far , gar ) ;
707 mri_free(gim) ; mri_free(fim) ; return cc ;
708 }
709
710 /*----------------------------------------------------------------*/
711 /*! eta^2 (Cohen, NeuroImage 2008) 25 Jun 2010 [rickr]
712 *
713 * eta^2 = 1 - SUM[ (a_i - m_i)^2 + (b_i - m_i)^2 ]
714 * ------------------------------------
715 * SUM[ (a_i - M )^2 + (b_i - M )^2 ]
716 *
717 * where o a_i and b_i are the vector elements
718 * o m_i = (a_i + b_i)/2
719 * o M = mean across both vectors
720 -----------------------------------------------------------------*/
THD_eta_squared(int n,float * x,float * y)721 float THD_eta_squared( int n, float *x , float *y )
722 {
723 double num=0.0f , denom = 0.0f ;
724 float gm=0.0f , lm, vv, ww;
725 int ii ;
726
727 for( ii=0 ; ii < n ; ii++ ){ gm += x[ii] + y[ii] ; }
728 gm /= (2*n) ;
729
730 for( ii=0 ; ii < n ; ii++ ){
731 lm = 0.5 * ( x[ii] + y[ii] ) ;
732 vv = (x[ii]-lm); ww = (y[ii]-lm);
733 num += ( vv*vv + ww*ww );
734 vv = (x[ii]-gm); ww = (y[ii]-gm);
735 denom += ( vv*vv + ww*ww );
736 }
737
738 if( num < 0.0f || denom <= 0.0f || num >= denom ) return 0.0f ;
739 return 1.0 - num/denom ;
740 }
741
742 /*----------------------------------------------------------------*/
743 /*! eta^2 (Cohen, NeuroImage 2008) 16 Jun 2011 [rickr]
744 *
745 * Same as THD_eta_squared, but allow a mask and return a double.
746 *
747 * eta^2 = 1 - SUM[ (a_i - m_i)^2 + (b_i - m_i)^2 ]
748 * ------------------------------------
749 * SUM[ (a_i - M )^2 + (b_i - M )^2 ]
750 *
751 * where o a_i and b_i are the vector elements
752 * o m_i = (a_i + b_i)/2
753 * o M = mean across both vectors
754 -----------------------------------------------------------------*/
THD_eta_squared_masked(int n,float * x,float * y,byte * mask)755 double THD_eta_squared_masked( int n, float *x , float *y, byte *mask )
756 {
757 double num=0.0f , denom = 0.0f ;
758 float gm=0.0f , lm, vv, ww;
759 int ii, nm ;
760
761 for( ii=0, nm=0 ; ii < n ; ii++ )
762 if( !mask || mask[ii] ) { gm += x[ii] + y[ii] ; nm++ ; }
763
764 if( nm == 0 ) return 0.0 ; /* bail on empty mask */
765
766 gm /= (2*nm) ;
767
768 for( ii=0 ; ii < n ; ii++ ){
769 if( mask && !mask[ii] ) continue; /* skip any not in mask */
770
771 lm = 0.5f * ( x[ii] + y[ii] ) ;
772 vv = (x[ii]-lm); ww = (y[ii]-lm); num += ( vv*vv + ww*ww );
773 vv = (x[ii]-gm); ww = (y[ii]-gm); denom += ( vv*vv + ww*ww );
774 }
775
776 if( num < 0.0 || denom <= 0.0 || num >= denom ) return 0.0 ;
777
778 return 1.0 - num/denom ;
779 }
780
781 /*----------------------------------------------------------------*/
782 /*! dice coefficient (allow for mask) 27 Oct 2015 [rickr]
783 * (float input)
784 *
785 * dice = 2 * size(interection) / ( size(A) + size(B) )
786 *
787 * Simply count voxels for the 3 cases.
788 -----------------------------------------------------------------*/
THD_dice_coef_f_masked(int n,float * x,float * y,byte * mask)789 float THD_dice_coef_f_masked( int n, float *x , float *y, byte *mask )
790 {
791 int nA=0, nB=0, nAB=0;
792 int ii;
793
794 for( ii=0 ; ii < n ; ii++ ) {
795 if( mask && !mask[ii] ) continue;
796 if( x[ii] ) {
797 nA++;
798 if( y[ii] ) { nB++; nAB++; }
799 } else
800 if ( y[ii] ) nB++;
801 }
802
803 if( (nA+nB) > 0 ) return 2.0*nAB/(nA+nB);
804 return 0.0;
805 }
806
807 /****************************************************************************/
808 /*** Histogram-based measurements of dependence between two float arrays. ***/
809 /****************************************************************************/
810 /*--------------------------------------------------------------------------*/
811 /* Extensive changes from the olde method:
812 * values below bot and above top are not used
813 * histogram can have bot and top bins unequal, and in-between ones
814 equal (if use_xyclip != 0)
815 * histogram can have all unequal bin sizes (if nxybin > 2)
816 ----------------------------------------------------------------------------*/
817 #define LINHIST /* do linear spread in histogram */
818 #undef WW
819 #define WW(i) ((w==NULL) ? 1.0f : w[i]) /* weight function for i'th datum */
820
821 static float *xc=NULL , *yc=NULL , *xyc=NULL , nww=0.0f ;
822 static int nbin=0 , nbp=0 , nbm=0 ;
823
824 static float *xbin=NULL , *ybin=NULL ;
825 static int nxybin=0 ;
826
827 static int use_xyclip=0 ;
828 static float xclip_bot, xclip_top ;
829 static float yclip_bot, yclip_top ;
830
831 #undef XYC
832 #define XYC(p,q) xyc[(p)+(q)*nbp]
833
834 #ifndef WAY_BIG
835 # define WAY_BIG 1.e+10
836 #endif
837
838 #undef GOODVAL
839 #define GOODVAL(x) ((x) < WAY_BIG) /* x is not preposterous */
840
841 #undef RANGVAL
842 #define RANGVAL(x,b,t) ((x) >= (b) && (x) <= (t)) /* x between b and t */
843
844 #undef FREEIF
845 #define FREEIF(x) \
846 do{ if((x)!=NULL){ free((void *)(x)); (x)=NULL; } } while(0)
847
848 /*--------------------------------------------------------------------------*/
849 /*! Clear the internal 2D histogram (but not the bin settings).
850 ----------------------------------------------------------------------------*/
851
clear_2Dhist(void)852 void clear_2Dhist(void)
853 {
854 FREEIF(xc); FREEIF(yc); FREEIF(xyc); nbin = nbp = nbm = 0; nww = 0.0f; return;
855 }
856
857 /*--------------------------------------------------------------------------*/
858
859 static double hpow = 0.33333333333 ;
set_2Dhist_hpower(double hh)860 void set_2Dhist_hpower( double hh )
861 {
862 hpow = (hh > 0.0 && hh < 1.0) ? hh
863 : 0.33333333333 ;
864 clear_2Dhist() ;
865 }
866
867 /*--------------------------------------------------------------------------*/
868
869 static int nhbin = 0 ;
set_2Dhist_hbin(int nn)870 void set_2Dhist_hbin( int nn ){ nhbin = nn; clear_2Dhist(); }
get_2Dhist_hbin(void)871 int get_2Dhist_hbin( void ){ return nhbin ; }
872
873 /*--------------------------------------------------------------------------*/
874 /*! Retrieve the 2D histogram built previously in build_2Dhist().
875 - Return value is the number of bins in each direction (may be 0).
876 - *xyhist is points to internal 2D array (may be NULL).
877 ----------------------------------------------------------------------------*/
878
retrieve_2Dhist(float ** xyhist)879 int retrieve_2Dhist( float **xyhist )
880 {
881 if( xyhist == NULL ) return 0 ;
882 *xyhist = xyc ; return nbp ;
883 }
884
885 /*--------------------------------------------------------------------------*/
886
retrieve_2Dhist1(float ** xh,float ** yh)887 int retrieve_2Dhist1( float **xh , float **yh )
888 {
889 if( xh == NULL || yh == NULL ) return 0 ;
890 *xh = xc ; *yh = yc ; return nbp ;
891 }
892
893 /*--------------------------------------------------------------------------*/
894 /*! Set the bins to be used in 2D histogram-izing in build_2Dhist().
895 - nb = number of bins; if nb < 3, equal sized bins will be used
896 and the internal bin arrays will be cleared
897 - xb = [0..nb] (nb+1 elements) = boundary pts for x-axis bins
898 - yb = [0..nb] (nb+1 elements) = boundary pts for y-axis bins
899 - Both xb[] and yb[] must be strictly increasing arrays!
900 - Also see get_2Dhist_xybin()
901 ----------------------------------------------------------------------------*/
902
set_2Dhist_xybin(int nb,float * xb,float * yb)903 void set_2Dhist_xybin( int nb , float *xb , float *yb )
904 {
905 int ii ;
906 FREEIF(xbin) ; FREEIF(ybin) ; nxybin = 0 ;
907 if( nb > 2 && xb != NULL && yb != NULL ){
908 for( ii=1 ; ii <= nb ; ii++ ) /* check that bin sizes are OK */
909 if( xb[ii-1] >= xb[ii] || yb[ii-1] > yb[ii] ) break ;
910 if( ii > nb ){
911 nxybin = nb ;
912 xbin = (float *)malloc(sizeof(float)*(nb+1)) ;
913 ybin = (float *)malloc(sizeof(float)*(nb+1)) ;
914 memcpy( xbin , xb , sizeof(float)*(nb+1) ) ;
915 memcpy( ybin , yb , sizeof(float)*(nb+1) ) ;
916 } else {
917 WARNING_message("set_2Dhist_xybin: illegal inputs!") ;
918 }
919 }
920 return ;
921 }
922
923 /*--------------------------------------------------------------------------*/
924
set_2Dhist_xybin_eqwide(int nb,float xbot,float xtop,float ybot,float ytop)925 void set_2Dhist_xybin_eqwide( int nb, float xbot,float xtop,float ybot,float ytop )
926 {
927 int ii ; float dx , dy ;
928
929 FREEIF(xbin) ; FREEIF(ybin) ; nxybin = 0 ;
930 if( nb > 2 && xbot < xtop && ybot < ytop ){
931 nxybin = nb ;
932 xbin = (float *)malloc(sizeof(float)*(nb+1)) ;
933 ybin = (float *)malloc(sizeof(float)*(nb+1)) ;
934 dx = (xtop-xbot) / nb ; dy = (ytop-ybot) / nb ;
935 for( ii=0 ; ii < nb ; ii++ ){
936 xbin[ii] = xbot + ii*dx ;
937 ybin[ii] = ybot + ii*dy ;
938 }
939 xbin[nb] = xtop ; ybin[nb] = ytop ;
940 #if 0
941 INFO_message("set_2Dhist_xybin_eqwide: %d %f..%f %f..%f",
942 nb,xbot,xtop,ybot,ytop ) ;
943 #endif
944 }
945 return ;
946 }
947
948 /*--------------------------------------------------------------------------*/
949
clipate(int nval,float * xar)950 static float_pair clipate( int nval , float *xar )
951 {
952 MRI_IMAGE *qim; float cbot,ctop, mmm , *qar; float_pair rr; int ii,nq;
953
954 ENTRY("clipate") ;
955
956 qim = mri_new_vol( nval,1,1 , MRI_float ) ; qar = MRI_FLOAT_PTR(qim) ;
957 for( ii=nq=0 ; ii < nval ; ii++ ) if( GOODVAL(xar[ii]) ) qar[nq++] = xar[ii];
958 qim->nx = qim->nvox = nq ;
959 if( nq < 666 ){ rr.a = 1.0f; rr.b = 0.0f; mri_free(qim); RETURN(rr); }
960 mmm = mri_min( qim ) ;
961 if( mmm >= 0.0f ){ /* for positive images */
962 cbot = THD_cliplevel( qim , 0.345f ) ;
963 ctop = mri_quantile ( qim , 0.966f ) ;
964 if( ctop > 4.321f*cbot ) ctop = 4.321f*cbot ;
965 } else { /* for images including negative values: no go */
966 cbot = 1.0f; ctop = 0.0f;
967 }
968 mri_free(qim) ;
969 rr.a = cbot ; rr.b = ctop ; RETURN(rr) ;
970 }
971
972 /*--------------------------------------------------------------------------*/
973
set_2Dhist_xyclip(int nval,float * xval,float * yval)974 void set_2Dhist_xyclip( int nval , float *xval , float *yval )
975 {
976 float_pair xcc , ycc ;
977
978 ENTRY("set_2Dhist_xyclip") ;
979
980 use_xyclip = 0 ;
981 if( nval < 666 || xval == NULL || yval == NULL ) EXRETURN ;
982
983 xcc = clipate( nval , xval ) ;
984 ycc = clipate( nval , yval ) ;
985
986 if( xcc.a >= xcc.b || ycc.a >= ycc.b ) EXRETURN ;
987
988 use_xyclip = 1 ; nxybin = 0 ;
989 xclip_bot = xcc.a ; xclip_top = xcc.b ;
990 yclip_bot = ycc.a ; yclip_top = ycc.b ;
991 EXRETURN ;
992 }
993
994 /*--------------------------------------------------------------------------*/
995
get_2Dhist_xyclip(float * xbc,float * xtc,float * ybc,float * ytc)996 int get_2Dhist_xyclip( float *xbc , float *xtc , float *ybc , float *ytc )
997 {
998 *xbc = xclip_bot ; *xtc = xclip_top ;
999 *ybc = yclip_bot ; *ytc = yclip_top ; return use_xyclip ;
1000 }
1001
1002 /*--------------------------------------------------------------------------*/
1003
eqhighate(int nb,int nval,float * xar,float * xb)1004 static int eqhighate( int nb , int nval , float *xar , float *xb )
1005 {
1006 float *xd , xbot,xtop , frac,fi , xtest ;
1007 int ii, pp, pbot,ptop , bsiz , pstart,pend , nxd ;
1008
1009 xd = (float *)malloc(sizeof(float)*nval) ;
1010 for( ii=nxd=0 ; ii < nval ; ii++ ) if( GOODVAL(xar[ii]) ) xd[nxd++] = xar[ii];
1011 if( nxd < 7*nb ){ free(xd); return 0; } /* bad */
1012 qsort_float( nxd , xd ) ;
1013
1014 /** scan for plateaus = runs of constant values longer than
1015 nominal bin width, at the bottom and top **/
1016
1017 xbot = xb[0] = xd[0] ;
1018 xtop = xb[nb] = xd[nxd-1] ; if( xtop <= xbot ){ free(xd); return 0; }
1019 bsiz = (nxd/nb) ;
1020
1021 xtest = xbot + (xtop-xbot)/(100.0f*nb) ;
1022 for( pp=1 ; pp < nxd && xd[pp] < xtest ; pp++ ) ; /*nada*/
1023 if( pp == nxd ){ free(xd); return 0; } /* data is constant? */
1024 pbot = (pp > bsiz) ? pp : 0 ;
1025
1026 xtest = xtop - (xtop-xbot)/(100.0f*nb) ;
1027 for( pp=nxd-2 ; pp > 0 && xd[pp] > xtest ; pp-- ) ; /*nada*/
1028 if( pp <= pbot ){ free(xd); return 0; } /* something screwy */
1029 ptop = (nxd-1-pp > bsiz) ? pp : nxd-1 ;
1030
1031 if( pbot > 0 ){
1032 xb[1] = 0.999999f*xd[pbot] + 0.000001f*xbot ; pstart = 2 ;
1033 } else {
1034 pstart = 1 ;
1035 }
1036 if( ptop < nxd-1 ){
1037 xb[nb-1] = 0.999999f*xd[ptop] + 0.000001f*xtop ; pend = nb-2 ;
1038 } else {
1039 pend = nb-1 ;
1040 }
1041
1042 frac = (ptop-pbot)/(pend-pstart+2.0f) ;
1043 for( pp=pstart ; pp <= pend ; pp++ ){
1044 fi = pbot + frac*(pp-pstart+1.0f) ; ii = (int)fi ; fi = fi - ii ;
1045 xb[pp] = (1.0f-fi) * xd[ii] + fi * xd[ii+1] ;
1046 }
1047
1048 free(xd) ; return nb ;
1049 }
1050
1051 /*--------------------------------------------------------------------------*/
1052
set_2Dhist_xybin_eqhigh(int nb,int nval,float * xval,float * yval)1053 void set_2Dhist_xybin_eqhigh( int nb, int nval, float *xval, float *yval )
1054 {
1055 int ii,jj ;
1056
1057 FREEIF(xbin) ; FREEIF(ybin) ; nxybin = 0 ;
1058 if( nb < 3 || nval < 9*nb || xval == NULL || yval == NULL ) return ;
1059
1060 nxybin = nb ;
1061 xbin = (float *)malloc(sizeof(float)*(nb+1)) ;
1062 ybin = (float *)malloc(sizeof(float)*(nb+1)) ;
1063
1064 ii = eqhighate( nb , nval , xval , xbin ) ;
1065 jj = eqhighate( nb , nval , yval , ybin ) ;
1066
1067 if( ii == 0 || jj == 0 ){
1068 FREEIF(xbin) ; FREEIF(ybin) ; nxybin = 0 ; /* bad things happened */
1069 }
1070 return ;
1071 }
1072
1073 /*--------------------------------------------------------------------------*/
1074 /*! Sets pointers to the bins in 2D histogram-ization, and return value
1075 is the number of bins (if 0, equal size bins are in use).
1076 Do NOT screw with the contents of the pointers: these are internal arrays!
1077 Also see set_2Dhist_xybin().
1078 ----------------------------------------------------------------------------*/
1079
get_2Dhist_xybin(float ** xb,float ** yb)1080 int get_2Dhist_xybin( float **xb , float **yb )
1081 {
1082 if( xb != NULL ) *xb = xbin ;
1083 if( yb != NULL ) *yb = ybin ;
1084 return nxybin ;
1085 }
1086
1087 /*--------------------------------------------------------------------------*/
1088 /*! Load 2D histogram of x[0..n-1] and y[0..n-1], each point optionally
1089 weighted by w[0..n-1] (weights are all 1 if w==NULL).
1090 Used in the histogram-based measures of dependence between x[] and y[i].
1091 If something is bad on input, nbin is set to 0. Otherwise, these global
1092 variables are set:
1093 - nbin = # of bins, nbp = nbin+1
1094 - nww = sum of the weights used
1095 - xc = marginal histogram of x[], for xc[0..nbin] (nbp points in)
1096 - yc = marginal histogram of y[], for yc[0..nbin] (each direction)
1097 - xyc = joint histogram of (x[],y[]), for XYC(0..nbin,0..nbin)
1098 - The histograms are normalized (by 1/nww) to have sum==1
1099 - Histogram can be retrieved by retrieve_2Dhist() and can be
1100 erased by clear_2Dhist()
1101 - Default number of equal-spaced bins in each direction is n^(1/3)
1102 - the exponent can be changed with set_2Dhist_hpower()
1103 - you can set the number of bins with set_2Dhist_hbin()
1104 - you can set unequal bins with set_2Dhist_xybin()
1105 - x[] values outside the range xbot..xtop (inclusive) or outside
1106 the unequal bins set in set_2Dhist_xybin() (if applicable) will
1107 not be used in the histogram; mutatis mutandum for y[]
1108 ----------------------------------------------------------------------------*/
1109
stuff_2Dhist(int n,float xbot,float xtop,float * x,float ybot,float ytop,float * y,float * w,int addon)1110 static void stuff_2Dhist( int n, float xbot,float xtop,float *x,
1111 float ybot,float ytop,float *y, float *w, int addon )
1112 {
1113 register int ii,jj,kk ;
1114 float xb,xi , yb,yi , xx,yy , x1,y1 , ww ;
1115 byte *good ; int ngood , xyclip ;
1116 static float xxbot=0.0f,xxtop=0.0f , yybot=0.0f,yytop=0.0f ;
1117
1118 ENTRY("stuff_2Dhist") ;
1119
1120 if( n <= 0 || x == NULL || y == NULL ) EXRETURN ;
1121
1122 /* get the min..max range for x and y data? */
1123
1124 good = (byte *)malloc(sizeof(byte)*n) ; /* 28 Feb 2007 */
1125 for( ii=0 ; ii < n ; ii++ )
1126 good[ii] = GOODVAL(x[ii]) && GOODVAL(y[ii]) ;
1127
1128 addon = addon && (xxbot < xxtop) && (yybot < yytop) && (nww > 0.0f) ;
1129
1130 if( addon ){
1131
1132 xbot = xxbot ; ybot = yybot ; xtop = xxtop ; ytop = yytop ;
1133
1134 } else {
1135
1136 clear_2Dhist() ;
1137
1138 if( nxybin > 2 ){ /* unequal bins */
1139 xbot = xbin[0] ; xtop = xbin[nxybin] ;
1140 } else if( xbot >= xtop ){ /* equal bins, and must find range */
1141 xbot = WAY_BIG ; xtop = -WAY_BIG ;
1142 for( ii=0 ; ii < n ; ii++ )
1143 if( good[ii] ){
1144 if( x[ii] > xtop ) xtop = x[ii] ;
1145 else if( x[ii] < xbot ) xbot = x[ii] ;
1146 }
1147 }
1148 if( xbot >= xtop ){ clear_2Dhist(); free(good); EXRETURN; }
1149
1150 if( nxybin > 2 ){
1151 ybot = ybin[0] ; ytop = ybin[nxybin] ;
1152 } else if( ybot >= ytop ){
1153 ybot = WAY_BIG ; ytop = -WAY_BIG ;
1154 for( ii=0 ; ii < n ; ii++ )
1155 if( good[ii] ){
1156 if( y[ii] > ytop ) ytop = y[ii] ;
1157 else if( y[ii] < ybot ) ybot = y[ii] ;
1158 }
1159 }
1160 if( ybot >= ytop ){ free(good); EXRETURN; }
1161
1162 xxbot = xbot ; xxtop = xtop ; yybot = ybot ; yytop = ytop ;
1163
1164 }
1165
1166 /*-- count number of good values left in range (in both x and y) --*/
1167
1168 memset(good,0,n) ;
1169 for( ngood=ii=0 ; ii < n ; ii++ ){
1170 if( RANGVAL(x[ii],xbot,xtop) && RANGVAL(y[ii],ybot,ytop) && WW(ii) > 0.0f ){
1171 good[ii] = 1 ; ngood++ ;
1172 }
1173 }
1174 if( ngood == 0 ){ free(good) ; EXRETURN ; }
1175
1176 /*--- allocate space for 2D and 1D histograms ---*/
1177
1178 if( !addon ){
1179 if( nxybin > 2 ){ /* array size is fixed by bins */
1180 nbin = nxybin ;
1181 } else { /* compute new histo array size */
1182 nbin = (nhbin > 2) ? nhbin : (int)pow((double)n,hpow) ;
1183 if( nbin > 255 ) nbin = 255; else if( nbin < 3 ) nbin = 3;
1184 }
1185 nbp = nbin+1 ; nbm = nbin-1 ;
1186
1187 FREEIF(xc) ; FREEIF(yc) ; FREEIF(xyc) ;
1188
1189 xc = (float *)calloc(sizeof(float),nbp) ;
1190 yc = (float *)calloc(sizeof(float),nbp) ;
1191 xyc = (float *)calloc(sizeof(float),nbp*nbp) ;
1192 nww = 0.0f ;
1193 }
1194
1195 /*--------------- make the 2D and 1D histograms ---------------*/
1196
1197 xyclip = (nxybin <= 2) && use_xyclip &&
1198 (xbot < xclip_bot) && (xclip_bot < xclip_top) && (xclip_top < xtop) &&
1199 (ybot < yclip_bot) && (yclip_bot < yclip_top) && (yclip_top < ytop) ;
1200
1201 if( nxybin <= 0 && !xyclip ){ /*------------ equal size bins ------------*/
1202
1203 #if 0
1204 if(PRINT_TRACING){
1205 char str[256];
1206 sprintf(str,"equal size bins: xbot=%g xtop=%g ybot=%g ytop=%g nbin=%d nval=%d ngood=%d",
1207 xbot,xtop,ybot,ytop,nbin,n,ngood);
1208 STATUS(str);
1209 }
1210 #endif
1211
1212 xb = xbot ; xi = nbm/(xtop-xbot) ;
1213 yb = ybot ; yi = nbm/(ytop-ybot) ;
1214 for( ii=0 ; ii < n ; ii++ ){
1215 if( !good[ii] ) continue ;
1216 xx = (x[ii]-xb)*xi ;
1217 jj = (int)xx ; xx = xx - jj ; x1 = 1.0f-xx ;
1218 yy = (y[ii]-yb)*yi ;
1219 kk = (int)yy ; yy = yy - kk ; y1 = 1.0f-yy ;
1220 ww = WW(ii) ; nww += ww ;
1221
1222 #ifdef LINHIST
1223 xc[jj] += x1*ww ; xc[jj+1] += xx*ww ;
1224 yc[kk] += (y1*ww); yc[kk+1] += (yy*ww);
1225
1226 XYC(jj ,kk ) += x1*(y1*ww) ;
1227 XYC(jj+1,kk ) += xx*(y1*ww) ;
1228 XYC(jj ,kk+1) += x1*(yy*ww) ;
1229 XYC(jj+1,kk+1) += xx*(yy*ww) ;
1230 #else
1231 xc[jj] += ww ; yc[kk] += ww ; XYC(jj,kk) += ww ;
1232 #endif
1233 }
1234
1235 } else if( xyclip ){ /*------------ mostly equal bins ----------------*/
1236
1237 float xbc=xclip_bot , xtc=xclip_top , ybc=yclip_bot , ytc=yclip_top ;
1238
1239 #if 0
1240 if(PRINT_TRACING){
1241 char str[256];
1242 sprintf(str,"mostly equal bins: xbc=%g xtc=%g ybc=%g ytc=%g nbin=%d",
1243 xbc,xtc,ybc,ytc,nbin) ; STATUS(str);
1244 }
1245 #endif
1246
1247 xi = (nbin-2.000001f)/(xtc-xbc) ;
1248 yi = (nbin-2.000001f)/(ytc-ybc) ;
1249 for( ii=0 ; ii < n ; ii++ ){
1250 if( !good[ii] ) continue ;
1251 xx = x[ii] ;
1252 if( xx < xbc ){ jj = 0 ; xx = 0.0f ; }
1253 else if( xx > xtc ){ jj = nbm ; xx = 1.0f ; }
1254 else { xx = 1.0f+(xx-xbc)*xi; jj = (int)xx; xx = xx - jj; }
1255 yy = y[ii] ;
1256 if( yy < ybc ){ kk = 0 ; yy = 0.0f ; }
1257 else if( yy > ytc ){ kk = nbm ; yy = 1.0f ; }
1258 else { yy = 1.0f+(yy-ybc)*yi; kk = (int)yy; yy = yy - kk; }
1259
1260 x1 = 1.0f-xx ; y1 = 1.0f-yy ; ww = WW(ii) ; nww += ww ;
1261
1262 #ifdef LINHIST
1263 xc[jj] += x1*ww ; xc[jj+1] += xx*ww ;
1264 yc[kk] += (y1*ww); yc[kk+1] += (yy*ww);
1265
1266 XYC(jj ,kk ) += x1*(y1*ww) ;
1267 XYC(jj+1,kk ) += xx*(y1*ww) ;
1268 XYC(jj ,kk+1) += x1*(yy*ww) ;
1269 XYC(jj+1,kk+1) += xx*(yy*ww) ;
1270 #else
1271 xc[jj] += ww ; yc[kk] += ww ; XYC(jj,kk) += ww ;
1272 #endif
1273 }
1274
1275 } else { /*-------------------- unequal size bins --------------------*/
1276
1277 float *xdup, *ydup, *xv, *yv, *wv, frac,val,xt,yt ;
1278 int *xn, *yn, *xin, *yin , ib,nb , maxsort ;
1279
1280 maxsort = 32*nxybin ; if( maxsort > 2048 ) maxsort = 2048 ;
1281 xdup = (float *)malloc(sizeof(float)*maxsort) ;
1282 ydup = (float *)malloc(sizeof(float)*maxsort) ;
1283 xv = (float *)malloc(sizeof(float)*maxsort) ;
1284 yv = (float *)malloc(sizeof(float)*maxsort) ;
1285 wv = (float *)malloc(sizeof(float)*maxsort) ;
1286 xn = (int *) malloc(sizeof(float)*maxsort) ;
1287 yn = (int *) malloc(sizeof(float)*maxsort) ;
1288 xin = (int *) malloc(sizeof(float)*maxsort) ;
1289 yin = (int *) malloc(sizeof(float)*maxsort) ; /* not yang */
1290
1291 /* outer loop: process points starting at index = ib */
1292
1293 for( ib=0 ; ib < n ; ){
1294
1295 /* extract up to maxsort good points from the data arrays */
1296
1297 for( nb=0,ii=ib ; ii < n && nb < maxsort ; ii++ ){
1298 if( good[ii] ){
1299 wv[nb]=WW(ii); xdup[nb]=x[ii]; ydup[nb]=y[ii]; xn[nb]=yn[nb]=nb;
1300 nb++; /* separate for clarity */
1301 }
1302 }
1303 #if 0
1304 INFO_message("extracted %d data points up to index=%d",nb,ii) ;
1305 #endif
1306 ib = ii ; /* where to start extracting next outer loop */
1307 if( nb == 0 ) break ; /* didn't find any good data ==> done! */
1308
1309 /* sort the extracted x data values in xdup[],
1310 and keep track of where they came from in xn[] */
1311
1312 qsort_floatint(nb,xdup,xn) ;
1313
1314 /* scan through sorted xdup[] values,
1315 which will go into bin #kk which is xb..xt;
1316 store the bin index and the fractional location within the bin;
1317 the reason for doing the sorting is that finding the bin index
1318 kk will be efficient
1319 (don't have to search from start as we would for unordered values) */
1320
1321 xb = xbin[0] ; xt = xbin[1] ; kk=0 ; frac = 1.0f/(xt-xb) ;
1322 for( ii=0 ; ii < nb ; ii++ ){
1323 val = xdup[ii] ;
1324 if( val > xt ){ /* not in the xb..xt bin, so scan up until it is */
1325 for( kk++ ; kk+1 < nxybin && val > xbin[kk+1] ; kk++ ) ; /*nada*/
1326 xb = xbin[kk] ; xt = xbin[kk+1] ; frac = 1.0f/(xt-xb) ;
1327 }
1328 jj = xn[ii] ; /* index in unsorted array */
1329 xv[jj] = frac*(val-xb) ; /* fractional position in bin */
1330 xin[jj] = kk ; /* bin index */
1331 }
1332
1333 /* repeat the above for the y arrays */
1334
1335 qsort_floatint(nb,ydup,yn) ;
1336 yb = ybin[0] ; yt = ybin[1] ; kk=0 ; frac = 1.0f/(yt-yb) ;
1337 for( ii=0 ; ii < nb ; ii++ ){
1338 val = ydup[ii] ;
1339 if( val > yt ){ /* not in the yb..yt bin, so scan up until it is */
1340 for( kk++ ; kk+1 < nxybin && val > ybin[kk+1] ; kk++ ) ; /*nada*/
1341 yb = ybin[kk] ; yt = ybin[kk+1] ; frac = 1.0f/(yt-yb) ;
1342 }
1343 jj = yn[ii] ; /* index in unsorted array */
1344 yv[jj] = frac*(val-yb) ; /* fractional position in bin */
1345 yin[jj] = kk ; /* bin index */
1346 }
1347
1348 /* now distribute the values into the 1D and 2D histograms */
1349
1350 for( ii=0 ; ii < nb ; ii++ ){
1351 ww = wv[ii] ; nww += ww ;
1352 jj = xin[ii]; kk = yin[ii] ;
1353 xx = xv[ii] ; x1 = 1.0f-xx ;
1354 yy = yv[ii] ; y1 = 1.0f-yy ;
1355
1356 #ifdef LINHIST
1357 xc[jj] += x1*ww ; xc[jj+1] += xx*ww ;
1358 yc[kk] += (y1*ww); yc[kk+1] += (yy*ww);
1359
1360 XYC(jj ,kk ) += x1*(y1*ww) ;
1361 XYC(jj+1,kk ) += xx*(y1*ww) ;
1362 XYC(jj ,kk+1) += x1*(yy*ww) ;
1363 XYC(jj+1,kk+1) += xx*(yy*ww) ;
1364 #else
1365 xc[jj] += ww ; yc[kk] += ww ; XYC(jj,kk) += ww ;
1366 #endif
1367 }
1368
1369 } /* end of outer loop (ib) over blocks */
1370
1371 free(yin); free(xin); free(yn); free(xn);
1372 free(wv); free(yv); free(xv); free(ydup); free(xdup);
1373
1374 } /*----- end of test on equal or unequal size bins -----*/
1375
1376 free(good); EXRETURN;
1377 }
1378
1379 /*............................................................................*/
1380
build_2Dhist(int n,float xbot,float xtop,float * x,float ybot,float ytop,float * y,float * w)1381 void build_2Dhist( int n , float xbot,float xtop,float *x ,
1382 float ybot,float ytop,float *y , float *w )
1383 {
1384 if( n > 9 && x != NULL && y != NULL )
1385 stuff_2Dhist( n , xbot,xtop,x , ybot,ytop,y , w , 0 ) ;
1386 return ;
1387 }
1388
1389 /*............................................................................*/
1390
addto_2Dhist(int n,float xbot,float xtop,float * x,float ybot,float ytop,float * y,float * w)1391 void addto_2Dhist( int n , float xbot,float xtop,float *x ,
1392 float ybot,float ytop,float *y , float *w )
1393 {
1394 if( n > 0 && x != NULL && y != NULL )
1395 stuff_2Dhist( n , xbot,xtop,x , ybot,ytop,y , w , 1 ) ;
1396 return ;
1397 }
1398
1399 /*............................................................................*/
1400 /* scale histogram to have sum==1 */
1401
normalize_2Dhist(void)1402 void normalize_2Dhist(void)
1403 {
1404 if( nww > 0.0f && xyc != NULL && xc != NULL && yc != NULL ){
1405 register float ni ; register int nbq , ii ;
1406 ni = 1.0f / nww ;
1407 for( ii=0 ; ii < nbp ; ii++ ){ xc[ii] *= ni; yc[ii] *= ni; }
1408 nbq = nbp*nbp ;
1409 for( ii=0 ; ii < nbq ; ii++ ){ xyc[ii] *= ni; }
1410 }
1411 return ;
1412 }
1413
1414 /*--------------------------------------------------------------------------*/
1415 /*--------------------------------------------------------------------------*/
1416
build_byteized_vectors(int n,float xbot,float xtop,float * x,float ybot,float ytop,float * y)1417 MRI_IMAGE *build_byteized_vectors( int n ,
1418 float xbot,float xtop,float *x ,
1419 float ybot,float ytop,float *y )
1420 {
1421 register int ii,jj,kk , gg ;
1422 float xb,xi , yb,yi , xx,yy , x1,y1 ;
1423 byte *good ; int ngood , xyclip , nbm ;
1424 MRI_IMAGE *bim ; byte *xbar,*ybar ;
1425
1426 ENTRY("build_byteized_vectors") ;
1427
1428 /* bad inputs? */
1429
1430 if( n <= 1 || x == NULL || y == NULL ) RETURN(NULL) ;
1431
1432 /* get the min..max range for x data? */
1433
1434 STATUS("compute good[]") ;
1435 good = (byte *)malloc(sizeof(byte)*n) ;
1436 for( ii=0 ; ii < n ; ii++ )
1437 good[ii] = GOODVAL(x[ii]) && GOODVAL(y[ii]) ;
1438
1439 if( nxybin > 2 ){ /* unequal bins */
1440 xbot = xbin[0] ; xtop = xbin[nxybin] ;
1441 } else if( xbot >= xtop ){ /* equal bins, and must find range */
1442 xbot = WAY_BIG ; xtop = -WAY_BIG ;
1443 for( ii=0 ; ii < n ; ii++ )
1444 if( good[ii] ){
1445 if( x[ii] > xtop ) xtop = x[ii] ;
1446 else if( x[ii] < xbot ) xbot = x[ii] ;
1447 }
1448 }
1449 if( xbot >= xtop ){ free(good); RETURN(NULL); }
1450
1451 /* get the min..max range for y data? */
1452
1453 if( nxybin > 0 ){
1454 ybot = ybin[0] ; ytop = ybin[nxybin] ;
1455 } else if( ybot >= ytop ){
1456 ybot = WAY_BIG ; ytop = -WAY_BIG ;
1457 for( ii=0 ; ii < n ; ii++ )
1458 if( good[ii] ){
1459 if( y[ii] > ytop ) ytop = y[ii] ;
1460 else if( y[ii] < ybot ) ybot = y[ii] ;
1461 }
1462 }
1463 if( ybot >= ytop ){ free(good); RETURN(NULL); }
1464
1465 /*--- figure out quantization levels ---*/
1466
1467 if( nxybin > 2 ){ /* array size is fixed by bins */
1468 nbin = nxybin ;
1469 } else { /* compute new histo array size */
1470 nbin = (nhbin > 2) ? nhbin : (int)pow((double)n,hpow) ;
1471 }
1472 if( nbin > 255 ) nbin = 255; else if( nbin < 16 ) nbin = 16 ;
1473 nbm = nbin-1 ;
1474
1475 /*-- count number of good values left in range (in both x and y) --*/
1476
1477 memset(good,0,n) ;
1478 for( ngood=ii=0 ; ii < n ; ii++ ){
1479 if( RANGVAL(x[ii],xbot,xtop) && RANGVAL(y[ii],ybot,ytop) ){
1480 good[ii] = 1 ; ngood++ ;
1481 }
1482 }
1483 if( ngood < 3*nbin ){ free(good); RETURN(NULL); }
1484
1485 /*--------------- make the output data -----------------*/
1486
1487 bim = mri_new( ngood , 2 , MRI_byte ) ; /* zero filled */
1488 xbar = MRI_BYTE_PTR(bim) ;
1489 ybar = xbar + ngood ;
1490
1491 xyclip = (nxybin <= 0) && use_xyclip &&
1492 (xbot < xclip_bot) && (xclip_bot < xclip_top) && (xclip_top < xtop) &&
1493 (ybot < yclip_bot) && (yclip_bot < yclip_top) && (yclip_top < ytop) ;
1494
1495 if( nxybin <= 0 && !xyclip ){ /*------------ equal size bins ------------*/
1496
1497 xb = xbot ; xi = (nbin-0.01f)/(xtop-xbot) ;
1498 yb = ybot ; yi = (nbin-0.01f)/(ytop-ybot) ;
1499 for( gg=ii=0 ; ii < n ; ii++ ){
1500 if( !good[ii] ) continue ;
1501 xbar[gg] = (byte)((x[ii]-xb)*xi) ;
1502 ybar[gg] = (byte)((y[ii]-yb)*yi) ; gg++ ;
1503 }
1504
1505 } else if( xyclip ){ /*------------ mostly equal bins ----------------*/
1506
1507 float xbc=xclip_bot , xtc=xclip_top , ybc=yclip_bot , ytc=yclip_top ;
1508
1509 xi = (nbin-2.01f)/(xtc-xbc) ;
1510 yi = (nbin-2.01f)/(ytc-ybc) ;
1511 for( gg=ii=0 ; ii < n ; ii++ ){
1512 if( !good[ii] ) continue ;
1513 xx = x[ii] ;
1514 if( xx < xbc ){ xbar[gg] = 0 ; }
1515 else if( xx > xtc ){ xbar[gg] = nbm ; }
1516 else { xbar[gg] = (byte)(1.0f+(xx-xbc)*xi) ; }
1517 yy = y[ii] ;
1518 if( yy < ybc ){ ybar[gg] = 0 ; }
1519 else if( yy > ytc ){ ybar[gg] = nbm ; }
1520 else { ybar[gg] = (byte)(1.0f+(yy-ybc)*yi) ; }
1521 gg++ ;
1522 }
1523
1524 } else { /*-------------------- unequal size bins --------------------*/
1525
1526 float *xdup, *ydup, val,xt,yt ;
1527 int *xn, *yn, *xin, *yin , ib,nb , maxsort ;
1528
1529 maxsort = 16*nxybin ; if( maxsort > 512 ) maxsort = 512 ;
1530 xdup = (float *)malloc(sizeof(float)*maxsort) ;
1531 ydup = (float *)malloc(sizeof(float)*maxsort) ;
1532 xn = (int *) malloc(sizeof(float)*maxsort) ;
1533 yn = (int *) malloc(sizeof(float)*maxsort) ;
1534 xin = (int *) malloc(sizeof(float)*maxsort) ;
1535 yin = (int *) malloc(sizeof(float)*maxsort) ; /* not yang */
1536
1537 /* outer loop: process points starting at index = ib */
1538
1539 for( gg=ib=0 ; ib < n ; ){
1540
1541 /* extract up to maxsort good points from the data arrays */
1542
1543 for( nb=0,ii=ib ; ii < n && nb < maxsort ; ii++ ){
1544 if( good[ii] ){
1545 xdup[nb]=x[ii]; ydup[nb]=y[ii]; xn[nb]=yn[nb]=nb;
1546 nb++;
1547 }
1548 }
1549 ib = ii ; /* where to start extracting next outer loop */
1550 if( nb == 0 ) break ; /* didn't find any good data ==> done done! */
1551
1552 /* sort the extracted x data values in xdup[],
1553 and keep track of where they came from in xn[] */
1554
1555 qsort_floatint(nb,xdup,xn) ;
1556
1557 /* scan through sorted xdup[] values,
1558 which will go into bin #kk which is xb..xt;
1559 store the bin index;
1560 the reason for doing the sorting is that finding the bin index
1561 kk will be efficient
1562 (don't have to search from start as we would for unordered values) */
1563
1564 xb = xbin[0] ; xt = xbin[1] ; kk=0 ;
1565 for( ii=0 ; ii < nb ; ii++ ){
1566 val = xdup[ii] ;
1567 if( val > xt ){ /* not in the xb..xt bin, so scan up until it is */
1568 for( kk++ ; kk+1 < nxybin && val > xbin[kk+1] ; kk++ ) ; /*nada*/
1569 xb = xbin[kk] ; xt = xbin[kk+1] ;
1570 }
1571 jj = xn[ii] ; /* index in unsorted array */
1572 xin[jj] = kk ; /* bin index */
1573 }
1574
1575 /* repeat the above for the y arrays */
1576
1577 qsort_floatint(nb,ydup,yn) ;
1578 yb = ybin[0] ; yt = ybin[1] ; kk=0 ;
1579 for( ii=0 ; ii < nb ; ii++ ){
1580 val = ydup[ii] ;
1581 if( val > yt ){ /* not in the yb..yt bin, so scan up until it is */
1582 for( kk++ ; kk+1 < nxybin && val > ybin[kk+1] ; kk++ ) ; /*nada*/
1583 yb = ybin[kk] ; yt = ybin[kk+1] ;
1584 }
1585 jj = yn[ii] ; /* index in unsorted array */
1586 yin[jj] = kk ; /* bin index */
1587 }
1588
1589 /* now distribute the values into the 1D and 2D histograms */
1590
1591 for( ii=0 ; ii < nb ; ii++ ){
1592 xbar[gg] = (byte)xin[ii] ; ybar[gg] = (byte)yin[ii] ; gg++ ;
1593 }
1594
1595 } /* end of outer loop (ib) over blocks */
1596
1597 free(yin); free(xin); free(yn); free(xn); free(ydup); free(xdup);
1598
1599 } /*----- end of test on equal or unequal size bins -----*/
1600
1601 free(good); RETURN(bim);
1602 }
1603
1604 /*--------------------------------------------------------------------------*/
1605 /*--------------------------------------------------------------------------*/
1606
1607 static int ignore_zz = 0 ;
THD_correlate_ignore_zerozero(int z)1608 void THD_correlate_ignore_zerozero( int z ){ ignore_zz = z ; }
1609
1610 /*-------------------- Shannon entropy function --------------------*/
1611
1612 #define SHANENT(z) ( ( (z) <= 0.0f ) ? 0.0f : (z)*logf(z) )
1613
1614 /*--------------------------------------------------------------------------*/
1615 /*! Compute the mutual info between two vectors, sort of. [16 Aug 2006]
1616 ----------------------------------------------------------------------------*/
1617
THD_mutual_info_scl(int n,float xbot,float xtop,float * x,float ybot,float ytop,float * y,float * w)1618 float THD_mutual_info_scl( int n , float xbot,float xtop,float *x ,
1619 float ybot,float ytop,float *y , float *w )
1620 {
1621 register int ii,jj ;
1622 register float val ;
1623
1624 /*-- build 2D histogram --*/
1625
1626 build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ;
1627 if( nbin <= 0 || nww <= 0 ) return 0.0f ; /* something bad happened! */
1628 normalize_2Dhist() ;
1629
1630 /*-- compute MI from histogram --*/
1631
1632 val = 0.0f ;
1633 for( ii=0 ; ii < nbp ; ii++ ){
1634 for( jj=0 ; jj < nbp ; jj++ ){
1635 #if 0
1636 if( ii==0 && jj==0 && ignore_zz ) continue ;
1637 #endif
1638 val += SHANENT( XYC(ii,jj) ) ;
1639 }}
1640 return (1.4427f*val) ; /* units are bits, just for fun */
1641 }
1642
1643 /*--------------------------------------------------------------------------*/
1644 /*! Compute MI of x[] and y[i], with autoscaling. */
1645
THD_mutual_info(int n,float * x,float * y)1646 float THD_mutual_info( int n , float *x , float *y )
1647 {
1648 return THD_mutual_info_scl( n, 1.0f,-1.0f, x, 1.0f,-1.0f, y, NULL ) ;
1649 }
1650
1651 /*--------------------------------------------------------------------------*/
1652 /*! Compute the normalized mutual info between two vectors, sort of.
1653 Actually, returns H(x,y) / [ H(x)+H(y) ], which should be small if
1654 x and y are redundant and should be large if they are independent.
1655 ----------------------------------------------------------------------------*/
1656
THD_norm_mutinf_scl(int n,float xbot,float xtop,float * x,float ybot,float ytop,float * y,float * w)1657 float THD_norm_mutinf_scl( int n , float xbot,float xtop,float *x ,
1658 float ybot,float ytop,float *y , float *w )
1659 {
1660 register int ii,jj ;
1661 register float numer , denom ;
1662
1663 /*-- build 2D histogram --*/
1664
1665 build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ;
1666 if( nbin <= 0 || nww <= 0 ) return 0.0f ; /* something bad happened! */
1667 normalize_2Dhist() ;
1668
1669 /*-- compute NMI from histogram --*/
1670
1671 denom = numer = 0.0f ;
1672 for( ii=0 ; ii < nbp ; ii++ ){
1673 denom += SHANENT( xc[ii] ) + SHANENT( yc[ii] ) ;
1674 for( jj=0 ; jj < nbp ; jj++ ){
1675 #if 0
1676 if( ii==0 && jj==0 && ignore_zz ) continue ;
1677 #endif
1678 numer += SHANENT( XYC(ii,jj) ) ;
1679 }
1680 }
1681 if( denom != 0.0f ) denom = numer / denom ;
1682 return denom ;
1683 }
1684
1685 /*--------------------------------------------------------------------------*/
1686 /*! Compute NMI of x[] & y[i], with autoscaling: cf. THD_norm_mutinf_scl(). */
1687
THD_norm_mutinf(int n,float * x,float * y)1688 float THD_norm_mutinf( int n , float *x , float *y )
1689 {
1690 return THD_norm_mutinf_scl( n, 1.0f,-1.0f, x, 1.0f,-1.0f, y, NULL ) ;
1691 }
1692
1693 /*--------------------------------------------------------------------------*/
1694 /*! Compute the joint entropy between two vectors, sort of.
1695 ----------------------------------------------------------------------------*/
1696
THD_jointentrop_scl(int n,float xbot,float xtop,float * x,float ybot,float ytop,float * y,float * w)1697 float THD_jointentrop_scl( int n , float xbot,float xtop,float *x ,
1698 float ybot,float ytop,float *y , float *w )
1699 {
1700 register int ii,jj ;
1701 register float val ;
1702
1703 /*-- build 2D histogram --*/
1704
1705 build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ;
1706 if( nbin <= 0 || nww <= 0 ) return 0.0f ; /* something bad happened! */
1707 normalize_2Dhist() ;
1708
1709 /*-- compute MI from histogram --*/
1710
1711 val = 0.0f ;
1712 for( ii=0 ; ii < nbp ; ii++ ){
1713 for( jj=0 ; jj < nbp ; jj++ ){
1714 val -= SHANENT( XYC(ii,jj) ) ;
1715 }}
1716 return (1.4427f*val) ; /* units are bits, just for fun */
1717 }
1718
1719 /*--------------------------------------------------------------------------*/
1720 /*! Compute joint entropy of x[] and y[i], with autoscaling. */
1721
THD_jointentrop(int n,float * x,float * y)1722 float THD_jointentrop( int n , float *x , float *y )
1723 {
1724 return THD_jointentrop_scl( n, 1.0f,-1.0f, x, 1.0f,-1.0f, y, NULL ) ;
1725 }
1726
1727
1728 /*--------------------------------------------------------------------------*/
1729 /* Decide if THD_corr_ratio_scl() computes symmetric or unsymmetric. */
1730
1731 static int cr_mode = 1 ; /* 0=unsym 1=sym mult 2=sym add */
1732
THD_corr_ratio_mode(int mm)1733 void THD_corr_ratio_mode( int mm ){ cr_mode = mm ; }
1734
1735 /*--------------------------------------------------------------------------*/
1736 /*! Compute the correlation ratio between two vectors, sort of. [23 Aug 2006]
1737 ----------------------------------------------------------------------------*/
1738
THD_corr_ratio_scl(int n,float xbot,float xtop,float * x,float ybot,float ytop,float * y,float * w)1739 float THD_corr_ratio_scl( int n , float xbot,float xtop,float *x ,
1740 float ybot,float ytop,float *y , float *w )
1741 {
1742 register int ii,jj ;
1743 register float vv,mm ;
1744 float val , cyvar , uyvar , yrat,xrat ;
1745
1746 /*-- build 2D histogram --*/
1747
1748 build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ;
1749 if( nbin <= 0 ) return 0.0f ; /* something bad happened! */
1750 normalize_2Dhist() ;
1751
1752 /*-- compute CR(y|x) from histogram --*/
1753
1754 cyvar = 0.0f ;
1755 for( ii=0 ; ii < nbp ; ii++ ){
1756 if( xc[ii] > 0.0f ){
1757 vv = mm = 0.0f ; /* mm=E(y|x) vv=E(y^2|x) */
1758 for( jj=1 ; jj < nbp ; jj++ ){
1759 mm += (jj * XYC(ii,jj)) ; vv += jj * (jj * XYC(ii,jj)) ;
1760 }
1761 cyvar += (vv - mm*mm/xc[ii] ) ; /* Var(y|x) */
1762 }
1763 }
1764 vv = mm = uyvar = 0.0f ;
1765 for( jj=1 ; jj < nbp ; jj++ ){ /* mm=E(y) vv=E(y^2) */
1766 mm += (jj * yc[jj]) ; vv += jj * (jj * yc[jj]) ;
1767 }
1768 uyvar = vv - mm*mm ; /* Var(y) */
1769 yrat = (uyvar > 0.0f) ? cyvar/uyvar /* Var(y|x) / Var(y) */
1770 : 1.0f ;
1771
1772 if( cr_mode == 0 ) return (1.0f-yrat) ; /** unsymmetric **/
1773
1774 /** compute CR(x|y) also, for symmetrization **/
1775
1776 cyvar = 0.0f ;
1777 for( jj=0 ; jj < nbp ; jj++ ){
1778 if( yc[jj] > 0.0f ){
1779 vv = mm = 0.0f ; /* mm=E(x|y) vv=E(x^2|y) */
1780 for( ii=1 ; ii < nbp ; ii++ ){
1781 mm += (ii * XYC(ii,jj)) ; vv += ii * (ii * XYC(ii,jj)) ;
1782 }
1783 cyvar += (vv - mm*mm/yc[jj] ) ; /* Var(x|y) */
1784 }
1785 }
1786 vv = mm = uyvar = 0.0f ;
1787 for( ii=1 ; ii < nbp ; ii++ ){ /* mm=E(x) vv=E(x^2) */
1788 mm += (ii * xc[ii]) ; vv += ii * (ii * xc[ii]) ;
1789 }
1790 uyvar = vv - mm*mm ; /* Var(x) */
1791 xrat = (uyvar > 0.0f) ? cyvar/uyvar /* Var(x|y) / Var(x) */
1792 : 1.0f ;
1793
1794 if( cr_mode == 2 ) return (1.0f - 0.5f*(xrat+yrat)) ; /** additive **/
1795
1796 return (1.0f - xrat*yrat) ; /** multiplicative **/
1797 }
1798
1799 /*--------------------------------------------------------------------------*/
1800 /*! Compute CR of x[] and y[], using autoscaling. */
1801
THD_corr_ratio(int n,float * x,float * y)1802 float THD_corr_ratio( int n , float *x , float *y )
1803 {
1804 return THD_corr_ratio_scl( n , 1.0f,-1.0f , x, 1.0f,-1.0f , y , NULL ) ;
1805 }
1806
1807 /*--------------------------------------------------------------------------*/
1808 /*! Compute the Hellinger metric between two vectors, sort of.
1809 ----------------------------------------------------------------------------*/
1810
THD_hellinger_scl(int n,float xbot,float xtop,float * x,float ybot,float ytop,float * y,float * w)1811 float THD_hellinger_scl( int n , float xbot,float xtop,float *x ,
1812 float ybot,float ytop,float *y , float *w )
1813 {
1814 register int ii,jj ;
1815 register float val , pq ;
1816
1817 /*-- build 2D histogram --*/
1818
1819 build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ;
1820 if( nbin <= 0 || nww <= 0 ) return 0.0f ; /* something bad happened! */
1821 normalize_2Dhist() ;
1822
1823 /*-- compute metric from histogram --*/
1824
1825 val = 0.0f ;
1826 for( ii=0 ; ii < nbp ; ii++ ){
1827 for( jj=0 ; jj < nbp ; jj++ ){
1828 pq = XYC(ii,jj) ;
1829 if( pq > 0.0f ) val += sqrtf( pq * xc[ii] * yc[jj] ) ;
1830 }}
1831 return (1.0f-val) ;
1832 }
1833
1834 /*--------------------------------------------------------------------------*/
1835 /*! Compute Hellinger metric between x[] and y[i], with autoscaling. */
1836
THD_hellinger(int n,float * x,float * y)1837 float THD_hellinger( int n , float *x , float *y )
1838 {
1839 return THD_hellinger_scl( n, 1.0f,-1.0f, x, 1.0f,-1.0f, y, NULL ) ;
1840 }
1841
1842 /*--------------------------------------------------------------------------*/
1843 /*! Compute the Hellinger metric, mutual info, normalized MI, and
1844 symmetrized correlation ratio, and return all 4 (in that order)
1845 using the 1D and 2D histograms from build_2Dhist().
1846
1847 The first 3 values all measure the closeness of the joint histogram to
1848 the product of the marginals:
1849 - Hellinger is smaller when the joint is closer to the marginals' product
1850 - MI is also smaller when the joint is closer to the marginal's product
1851 - NMI is larger when the joint is closer to the marginal's product
1852 Correlation ratio (symmetrized by addition == CRA) is larger when
1853 the two variables are nonlinearly correlated.
1854
1855 As measures of association (generalized correlation): more closely
1856 associated variables correspond to larger Hellinger and MI and CRA,
1857 and to smaller NMI.
1858 *//*------------------------------------------------------------------------*/
1859
THD_helmicra_scl(int n,float xbot,float xtop,float * x,float ybot,float ytop,float * y,float * w)1860 float_quad THD_helmicra_scl( int n , float xbot,float xtop,float *x ,
1861 float ybot,float ytop,float *y , float *w )
1862 {
1863 register int ii,jj ;
1864 register float hel , pq , vv,uu ;
1865 float val , cyvar , uyvar , yrat,xrat ;
1866 float_quad hmc = {0.0f,0.0f,0.0f,0.f} ;
1867
1868 /*-- build 2D histogram --*/
1869
1870 build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ;
1871 if( nbin <= 0 || nww <= 0 ) return hmc ; /* something bad happened! */
1872 normalize_2Dhist() ;
1873
1874 /*-- compute Hel, MI, NMI from histogram --*/
1875
1876 hel = vv = uu = 0.0f ;
1877 for( ii=0 ; ii < nbp ; ii++ ){
1878 vv += SHANENT( xc[ii] ) + SHANENT( yc[ii] ) ;
1879 for( jj=0 ; jj < nbp ; jj++ ){
1880 pq = XYC(ii,jj) ;
1881 hel += sqrtf( pq * xc[ii] * yc[jj] ) ;
1882 uu += SHANENT( pq ) ;
1883 }
1884 }
1885 hmc.a = 1.0f - hel ; /* Hellinger */
1886 hmc.b = uu - vv ; /* MI */
1887 hmc.c = (vv != 0.0f) ? uu/vv : 0.0f ; /* NMI */
1888
1889 /*-- compute CR(y|x) from histogram --*/
1890
1891 cyvar = 0.0f ;
1892 for( ii=0 ; ii < nbp ; ii++ ){
1893 if( xc[ii] > 0.0f ){
1894 vv = uu = 0.0f ; /* uu=E(y|x) vv=E(y^2|x) */
1895 for( jj=1 ; jj < nbp ; jj++ ){
1896 uu += (jj * XYC(ii,jj)) ; vv += jj * (jj * XYC(ii,jj)) ;
1897 }
1898 cyvar += (vv - uu*uu/xc[ii] ) ; /* Var(y|x) */
1899 }
1900 }
1901 vv = uu = uyvar = 0.0f ;
1902 for( jj=1 ; jj < nbp ; jj++ ){ /* uu=E(y) vv=E(y^2) */
1903 uu += (jj * yc[jj]) ; vv += jj * (jj * yc[jj]) ;
1904 }
1905 uyvar = vv - uu*uu ; /* Var(y) */
1906 yrat = (uyvar > 0.0f) ? cyvar/uyvar /* Var(y|x) / Var(y) */
1907 : 1.0f ;
1908
1909 /** compute CR(x|y) also, for symmetrization **/
1910
1911 cyvar = 0.0f ;
1912 for( jj=0 ; jj < nbp ; jj++ ){
1913 if( yc[jj] > 0.0f ){
1914 vv = uu = 0.0f ; /* uu=E(x|y) vv=E(x^2|y) */
1915 for( ii=1 ; ii < nbp ; ii++ ){
1916 uu += (ii * XYC(ii,jj)) ; vv += ii * (ii * XYC(ii,jj)) ;
1917 }
1918 cyvar += (vv - uu*uu/yc[jj] ) ; /* Var(x|y) */
1919 }
1920 }
1921 vv = uu = uyvar = 0.0f ;
1922 for( ii=1 ; ii < nbp ; ii++ ){ /* uu=E(x) vv=E(x^2) */
1923 uu += (ii * xc[ii]) ; vv += ii * (ii * xc[ii]) ;
1924 }
1925 uyvar = vv - uu*uu ; /* Var(x) */
1926 xrat = (uyvar > 0.0f) ? cyvar/uyvar /* Var(x|y) / Var(x) */
1927 : 1.0f ;
1928
1929 hmc.d = 1.0f - 0.5f*(xrat+yrat) ; /** additive symmetrization **/
1930 return hmc ;
1931 }
1932
1933 /*--------------------------------------------------------------------------*/
1934 /*! see THD_helmicra_scl(). */
1935
THD_helmicra(int n,float * x,float * y)1936 float_quad THD_helmicra( int n , float *x , float *y )
1937 {
1938 return THD_helmicra_scl( n, 1.0f,-1.0f, x, 1.0f,-1.0f, y, NULL ) ;
1939 }
1940
1941 /*---------------------------------------------------------------------------*/
1942 /* Return the MI and NMI of the 2 arrays y0[] and y1[] as 'correlated'
1943 with the binary 0-1 variable. That is, does knowing a value is in
1944 group 0 or 1 help predict it?
1945 *//*-------------------------------------------------------------------------*/
1946
THD_binary_mutinfo(int n0,float * y0,int n1,float * y1)1947 float_pair THD_binary_mutinfo( int n0 , float *y0 , int n1 , float *y1 )
1948 {
1949 float_pair mmii = {-666.6f,666.6f} ;
1950 float *y0h , *y1h ; int nbin=64 , ii,kk ;
1951 float ymin,ymax , dy,dc, hsum,hjoint ;
1952
1953 if( n0 < nbin || y0 == NULL || n1 < nbin || y1 == NULL ) return mmii ;
1954
1955 /* data range */
1956
1957 ymin = ymax = y0[0] ;
1958 for( ii=1 ; ii < n0 ; ii++ ){
1959 if( ymin > y0[ii] ) ymin = y0[ii] ;
1960 else if( ymax < y0[ii] ) ymax = y0[ii] ;
1961 }
1962 for( ii=0 ; ii < n1 ; ii++ ){
1963 if( ymin > y1[ii] ) ymin = y1[ii] ;
1964 else if( ymax < y1[ii] ) ymax = y1[ii] ;
1965 }
1966 if( ymin >= ymax ) return mmii ;
1967
1968 /* histograms */
1969
1970 y0h = (float *)calloc(sizeof(float),nbin) ;
1971 y1h = (float *)calloc(sizeof(float),nbin) ;
1972
1973 dy = nbin / (ymax-ymin) ;
1974 dc = 1.0f / (float)(n0+n1) ;
1975 for( ii=0 ; ii < n0 ; ii++ ){
1976 kk = ( y0[ii] - ymin ) * dy ;
1977 if( kk < 0 ) kk = 0 ; else if( kk >= nbin ) kk = nbin-1 ;
1978 y0h[kk] += dc ;
1979 }
1980 for( ii=0 ; ii < n1 ; ii++ ){
1981 kk = (int)( ( y1[ii] - ymin ) * dy ) ;
1982 if( kk < 0 ) kk = 0 ; else if( kk >= nbin ) kk = nbin-1 ;
1983 y1h[kk] += dc ;
1984 }
1985
1986 /* calculate mutual info vs. 0/1 */
1987
1988 /* entropy of the 0/1 'variable' */
1989 hsum = - ( SHANENT(n0*dc) + SHANENT(n1*dc) ) ;
1990
1991 hjoint = 0.0f ;
1992 for( kk=0 ; kk < nbin ; kk++ ){
1993 hjoint -= ( SHANENT( y0h[kk] ) + SHANENT( y1h[kk] ) ) ;
1994 hsum -= SHANENT( y0h[kk] + y0h[kk] ) ;
1995 }
1996
1997 free(y0h) ; free(y1h) ;
1998
1999 mmii.a = hsum - hjoint ;
2000 if( hjoint > 0.0f ) mmii.b = hsum / hjoint ;
2001 return mmii ;
2002 }
2003
2004 /*---------------------------------------------------------------------------*/
2005 /*! Rank orders a collection of float arrays; ties get the average rank.
2006 The output overwrites the input.
2007 *//*-------------------------------------------------------------------------*/
2008
rank_order_float_arrays(int nar,int * nn,float ** aa)2009 void rank_order_float_arrays( int nar , int *nn , float **aa )
2010 {
2011 int ii,jj,kk , jbase , ntot,nmax,n1,ns,ib ;
2012 int *b ; /* workspaces */
2013 float *a , *c , cs ;
2014
2015 /*- handle special cases -*/
2016
2017 if( nar < 1 || nn == NULL || aa == NULL ) return ; /* junk inputs */
2018 if( nar == 1 ){
2019 rank_order_float( nn[0] , aa[0] ) ; return ; /* just one input */
2020 }
2021
2022 ntot = nmax = 0 ;
2023 for( jj=0 ; jj < nar ; jj++ ){
2024 ntot += nn[jj] ; if( nn[jj] > nmax ) nmax = nn[jj] ;
2025 }
2026 if( ntot < nar ) return ; /* bad inputs */
2027
2028 /*- make workspaces -*/
2029
2030 a = (float *)malloc(sizeof(float)*ntot) ;
2031 b = (int *)malloc(sizeof(int )*ntot) ;
2032 c = (float *)malloc(sizeof(float)*ntot) ;
2033
2034 for( kk=jj=0 ; jj < nar ; jj++ ){
2035 jbase = jj*nmax ;
2036 for( ii=0 ; ii < nn[jj] ; ii++,kk++ ){
2037 a[kk] = aa[jj][ii] ; /* data */
2038 b[kk] = ii + jbase ; /* where it came from */
2039 c[kk] = (float)kk ; /* default rank of kk-th value after sorting */
2040 }
2041 }
2042
2043 /*- sort input, carrying b along -*/
2044
2045 qsort_floatint( ntot , a , b ) ;
2046
2047 /* c[kk] now contains the global rank of the kk-th element,
2048 but c[] perhaps needs to modified to allow for ties ==> average rank */
2049
2050 n1 = ntot-1 ;
2051 for( ii=0 ; ii < n1 ; ii++ ){
2052 if( a[ii] == a[ii+1] ){ /* handle ties */
2053 cs = (ii)+(ii+1) ; ns = 2 ; ib = ii ; ii++ ;
2054 while( ii < n1 && a[ii] == a[ii+1] ){ ii++ ; ns++ ; cs += ii ; }
2055 for( cs/=ns ; ib <= ii ; ib++ ) c[ib] = cs ;
2056 }
2057 }
2058
2059 /* put c[] results back into the original arrays */
2060
2061 for( kk=0 ; kk < ntot ; kk++ ){
2062 jj = b[kk] / nmax ; /* which array did it come from? */
2063 ii = b[kk] % nmax ; /* which element in that arrary? */
2064 aa[jj][ii] = c[kk] ; /* replace data with rank */
2065 }
2066
2067 free(c) ; free(b) ; free(a) ; return ;
2068 }
2069
2070 /*---------------------------------------------------------------------------*/
2071
rank_order_2floats(int n1,float * a1,int n2,float * a2)2072 void rank_order_2floats( int n1 , float *a1 , int n2 , float *a2 )
2073 {
2074 int nn[2] ; float *aa[2] ;
2075
2076 if( n1 <= 0 || n2 <= 0 || a1 == NULL || a2 == NULL ) return ;
2077
2078 nn[0] = n1 ; nn[1] = n2 ;
2079 aa[0] = a1 ; aa[1] = a2 ;
2080 rank_order_float_arrays( 2 , nn , aa ) ;
2081 return ;
2082 }
2083
2084 /*----------------------------------------------------------------------------*/
2085 /* Given a set of bootstrap replicates, find the bias-corrected (BC)
2086 estimates of the 1-alpha confidence interval and the central value.
2087 * estim = actual estimate from the data
2088 * alpha = 0.05 is a typical value
2089 * nboot = number of replicates in eboot[] -- at least MAX(100,10/alpha)
2090 * eboot = array of replicates; will be sorted on output
2091 The return value (rval) is a triple:
2092 * rval.a = lower edge of confidence interval
2093 * rval.b = middle value = bias-corrected estimate
2094 * rval.c = upper edge of confidence interval
2095 If all values are returned as 0, the inputs were bad bad bad.
2096 *//*--------------------------------------------------------------------------*/
2097
2098
2099 #undef PHI
2100 #define PHI(x) (1.0-qg(x)) /* CDF of N(0,1) */
2101
2102 #undef PHINV
2103 #define PHINV(x) qginv(1.0-(x)) /* inverse CDF of N(0,1) */
2104
2105 #undef ZLIM
2106 #define ZLIM 0.5f
2107
THD_bootstrap_confinv(float estim,float alpha,int nboot,float * eboot)2108 float_triple THD_bootstrap_confinv( float estim , float alpha ,
2109 int nboot , float *eboot )
2110 {
2111 float_triple rval = {0.0f,0.0f,0.0f} ;
2112 int ii ; float z0 , zal , pp ;
2113
2114 ENTRY("THD_bootstrap_confinv") ;
2115
2116 if( nboot < 100 || eboot == NULL ) RETURN( rval ) ; /* bad user */
2117
2118 if( alpha <= 0.001f || alpha >= 0.9f ) alpha = 0.05f ; /* stupid user */
2119 alpha *= 0.5f ; /* 1-sided tail */
2120 if( (int)(alpha*nboot) < 5 ) alpha = 5.0f / nboot ; /* upward adjust */
2121 zal = PHINV(alpha) ; /* should be negative */
2122
2123 qsort_float( nboot , eboot ) ; /* increasing order */
2124
2125 for( ii=0 ; ii < nboot && eboot[ii] < estim ; ii++ ) ; /*nada*/
2126 if( ii <= 1 || ii >= nboot-1 ) RETURN( rval ) ; /* crummy data */
2127 z0 = PHINV( (ii+0.5f) / nboot ) ; /* ii = #values < estim */
2128 if( z0 < -ZLIM ) z0 = -ZLIM ; else if( z0 > ZLIM ) z0 = ZLIM ; /* limits */
2129
2130 pp = PHI( 2.0*z0 + zal ) * nboot ; /* lower edge */
2131 ii = (int)pp ; pp = pp - ii ; if( ii >= nboot-1 ) ii = nboot-2 ;
2132 rval.a = (1.0f-pp)*eboot[ii] + pp*eboot[ii+1] ;
2133
2134 pp = PHI( 2.0*z0 - zal ) * nboot ; /* upper edge */
2135 ii = (int)pp ; pp = pp - ii ; if( ii >= nboot-1 ) ii = nboot-2 ;
2136 rval.c = (1.0f-pp)*eboot[ii] + pp*eboot[ii+1] ;
2137
2138 pp = PHI( 2.0*z0 ) * nboot ; /* bias-corrected central estimate */
2139 ii = (int)pp ; pp = pp - ii ; if( ii >= nboot-1 ) ii = nboot-2 ;
2140 rval.b = (1.0f-pp)*eboot[ii] + pp*eboot[ii+1] ;
2141
2142 RETURN( rval ) ;
2143 }
2144
2145 /*----------------------------------------------------------------------------*/
2146
THD_bootstrap_biascorr(float estim,int nboot,float * eboot)2147 float THD_bootstrap_biascorr( float estim , int nboot , float *eboot )
2148 {
2149 int ii ; float z0 , pp ;
2150
2151 if( nboot < 50 || eboot == NULL ) return estim ; /* bad user */
2152
2153 qsort_float( nboot , eboot ) ; /* increasing order */
2154
2155 for( ii=0 ; ii < nboot && eboot[ii] < estim ; ii++ ) ; /*nada*/
2156 if( ii <= 1 || ii >= nboot-1 ) return estim ; /* crummy data */
2157 z0 = PHINV( (ii+0.5f) / nboot ) ; /* ii = #values < estim */
2158 if( z0 < -ZLIM ) z0 = -ZLIM ; else if( z0 > ZLIM ) z0 = ZLIM ; /* limits */
2159
2160 pp = PHI( 2.0*z0 ) * nboot ;
2161 ii = (int)pp ; pp = pp - ii ; if( ii >= nboot-1 ) ii = nboot-2 ;
2162 pp = (1.0f-pp)*eboot[ii] + pp*eboot[ii+1] ;
2163
2164 return pp ;
2165 }
2166
2167 /*----------------------------------------------------------------------------*/
2168 #undef DEMEAN
2169 #define DEMEAN(n,v) do{ register int i ; register float s ; \
2170 for( s=i=0 ; i < n ; i++ ) s += v[i] ; \
2171 s /= n ; for( i=0 ; i < n ; i++ ) v[i] -= s ; \
2172 } while(0)
2173
2174 #undef XPT
2175 #define XPT(q) ( (xtyp<=0) ? xx+(q)*nlen : xpt[q] )
2176
2177 #undef YPT
2178 #define YPT(q) ( (xtyp<=0) ? yy+(q)*nlen : ypt[q] )
2179
2180 /*----------------------------------------------------------------------------*/
2181 /* Return bootstrap BC estimate for correlation of a bunch of vectors.
2182 *//*--------------------------------------------------------------------------*/
2183
THD_bootstrap_vectcorr(int nlen,int nboot,int use_pv,int xtyp,int xnum,void * xp,int ynum,void * yp)2184 float THD_bootstrap_vectcorr( int nlen, int nboot, int use_pv, int xtyp,
2185 int xnum, void *xp , int ynum , void *yp )
2186 {
2187 float rval, rxy, **xar, **yar, *rboot, *xbar, *ybar, *uvec, *vvec ;
2188 int kb,ii,jj , dox,doy ;
2189 unsigned short xran[3] ;
2190 float *xx=NULL, **xpt=NULL, *yy=NULL, **ypt=NULL ; void *pvw=NULL ;
2191
2192 ENTRY("THD_bootstrap_vectcorr") ;
2193
2194 if( nlen < 3 || xnum < 1 || ynum < 1 || xp == NULL || yp == NULL )
2195 RETURN(0.0f) ;
2196
2197 if( xtyp <= 0 ){ xx = (float *) xp ; yy = (float *) yp ; }
2198 else { xpt = (float **)xp ; ypt = (float **)yp ; }
2199
2200 /* trivial case */
2201
2202 if( xnum == 1 && ynum == 1 ){
2203 rval = THD_pearson_corr(nlen,XPT(0),YPT(0)) ; RETURN(rval) ;
2204 }
2205
2206 /* compute mean vectors */
2207
2208 #pragma omp critical (MALLOC)
2209 { xbar = (float *)malloc(sizeof(float)*nlen) ;
2210 ybar = (float *)malloc(sizeof(float)*nlen) ; }
2211
2212 (void)mean_vector( nlen , xnum , xtyp , xp , xbar ) ;
2213 (void)mean_vector( nlen , ynum , xtyp , yp , ybar ) ;
2214
2215 dox = (xnum > 5) ; doy = (ynum > 5) ; /* which to bootstrap */
2216
2217 if( !dox && !doy ){
2218 rval = THD_pearson_corr(nlen,xbar,ybar) ;
2219 #pragma omp critical (MALLOC)
2220 { free(ybar) ; free(xbar) ; }
2221 RETURN(rval) ;
2222 }
2223
2224 if( nboot < 50 ) nboot = 50 ;
2225
2226 #pragma omp critical (MALLOC)
2227 { xar = (float **)malloc(sizeof(float *)*xnum ) ;
2228 yar = (float **)malloc(sizeof(float *)*ynum ) ;
2229 rboot = (float * )malloc(sizeof(float) *nboot) ;
2230 uvec = (float * )malloc(sizeof(float) *nlen ) ;
2231 vvec = (float * )malloc(sizeof(float) *nlen ) ;
2232 }
2233
2234 /* correlation with no resampling at all */
2235
2236 if( use_pv ){
2237 pvw = pv_get_workspace( nlen , MAX(xnum,ynum) ) ;
2238 xran[0] = (unsigned short)(xnum+ynum+nlen) ;
2239 xran[1] = 42731 ; xran[2] = 23172 ;
2240 (void)principal_vector( nlen, xnum, xtyp, xp, uvec, xbar, pvw, xran) ;
2241 (void)principal_vector( nlen, ynum, xtyp, yp, vvec, ybar, pvw, xran) ;
2242 } else {
2243 for( ii=0 ; ii < nlen ; ii++ ){ uvec[ii] = xbar[ii]; vvec[ii] = ybar[ii]; }
2244 }
2245 rxy = THD_pearson_corr( nlen , uvec , vvec ) ;
2246
2247 /* bootstrap correlations [selecting subsets from each set of vectors] */
2248
2249 for( kb=0 ; kb < nboot ; kb++ ){
2250 if( dox ){
2251 for( ii=0 ; ii < xnum ; ii++ ){ /* resample xx vectors */
2252 jj = nrand48(xran) % xnum ; xar[ii] = XPT(jj) ;
2253 }
2254 if( use_pv )
2255 (void)principal_vector( nlen, xnum, 1, xar, uvec, xbar, pvw, xran) ;
2256 else
2257 mean_vector( nlen , xnum , 1 , xar , uvec ) ;
2258 }
2259 if( doy ){
2260 for( ii=0 ; ii < ynum ; ii++ ){ /* resample yy vectors */
2261 jj = nrand48(xran) % ynum ; yar[ii] = YPT(jj) ;
2262 }
2263 if( use_pv )
2264 (void)principal_vector( nlen, ynum, 1, yar, vvec, ybar, pvw, xran) ;
2265 else
2266 mean_vector( nlen , ynum , 1 , yar , vvec ) ;
2267 }
2268 rboot[kb] = THD_pearson_corr( nlen , uvec , vvec ) ;
2269 }
2270
2271 /* final answer */
2272
2273 rval = THD_bootstrap_biascorr( rxy , nboot , rboot ) ;
2274
2275 /* cleanup the trash */
2276
2277 #pragma omp critical (MALLOC)
2278 { free(vvec); free(uvec); free(rboot);
2279 free(yar); free(xar); free(ybar); free(xbar); if(use_pv)free(pvw); }
2280
2281 RETURN(rval) ;
2282 }
2283