#include "mrilib.h" /*==============================================================================*/ /*========== The following functions were moved from afni_fimfunc.c -===========*/ /*===========================================================================*/ /*---------------------------------------------------------------------------*/ /*! Rank-order a float array, with ties getting the average rank. The output overwrites the input. [27 Jun 2010: modified to create/destroy workspace on each call] *//*-------------------------------------------------------------------------*/ void rank_order_float( int n , float *a ) { register int ii , ns , n1 , ib ; int *b ; /* workspaces */ float *c ; float cs ; /*- handle special cases -*/ if( a == NULL || n < 1 ) return ; /* meaningless input */ if( n == 1 ){ a[0] = 0.0f ; return ; } /* only one point!? */ /*- make workspaces -*/ b = (int *) malloc(sizeof(int )*n) ; c = (float *) malloc(sizeof(float)*n) ; /*- initialize position b[] and rank [c] -*/ for( ii=0 ; ii < n ; ii++ ) c[ii] = b[ii] = ii ; /*- sort input a, carrying b along; when sorted, a[i] is in a[b[i]] -*/ qsort_floatint( n , a , b ) ; /* see cs_sort_fi.c */ /* compute tied ranks into c[] */ n1 = n-1 ; for( ii=0 ; ii < n1 ; ii++ ){ if( a[ii] == a[ii+1] ){ /*- handle ties -*/ cs = 2*ii+1 ; /* sum of tied ranks */ ns = 2 ; /* number of tied ranks */ ib = ii ; ii++ ; /* scan ii up until no ties */ while( ii < n1 && a[ii] == a[ii+1] ){ ii++ ; ns++ ; cs += ii ; } /* rescan to fill tied ranks with average rank */ for( cs/=ns ; ib <= ii ; ib++ ) c[ib] = cs ; } } /*- store the rank c[] back into input a[] (if no ties, c[i] == i) -*/ for( ii=0 ; ii < n ; ii++ ) a[b[ii]] = c[ii] ; free(c) ; free(b) ; return ; } /*---------------------------------------------------------------------------*/ /*! Rank orders a[], subtracts the mean rank, and returns the sum-of-squares. -----------------------------------------------------------------------------*/ float spearman_rank_prepare( int n , float *a ) { register int ii ; register float rb , rs ; rank_order_float( n , a ) ; rb = 0.5f*(n-1) ; rs=0.0f ; for( ii=0 ; ii < n ; ii++ ){ a[ii] -= rb ; /* remove mean rank */ rs += a[ii]*a[ii] ; /* sum squares */ } return rs ; } /*---------------------------------------------------------------------------*/ /*! Prepare for quadrant correlation with a[]. -----------------------------------------------------------------------------*/ float quadrant_corr_prepare( int n , float *a ) { register int ii ; register float rb , rs ; rank_order_float( n , a ) ; rb = 0.5f*(n-1) ; rs=0.0f ; for( ii=0 ; ii < n ; ii++ ){ a[ii] = (a[ii] > rb) ? 1.0 : (a[ii] < rb) ? -1.0f : 0.0f ; rs += a[ii]*a[ii] ; } return rs ; } /*-----------------------------------------------------------------------------*/ /*! To Spearman (rank-order) correlate x[] with r[], first do rv = spearman_rank_prepare(n,r) ; then corr = spearman_rank_corr(n,x,rv,r) ; Note that these 2 routines are destructive (r and x are replaced by ranks). -------------------------------------------------------------------------------*/ float spearman_rank_corr( int n , float *x , float rv , float *r ) { register int ii ; register float ss ; float xv ; xv = spearman_rank_prepare( n , x ) ; if( xv <= 0.0f ) return 0.0f ; for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii] * r[ii] ; return ( ss/sqrtf(rv*xv) ) ; } /*------------------------------------------------------------------------------*/ /*! To do quadrant correlation of x[] with r[], first do rv = quadrant_corr_prepare(n,r) ; then corr = quadrant_corr(n,x,rv,r) ; Note that these 2 routines are destructive (r and x are modified). -------------------------------------------------------------------------------*/ float quadrant_corr( int n , float *x , float rv , float *r ) { register int ii ; register float ss ; float xv ; xv = quadrant_corr_prepare( n , x ) ; if( xv <= 0.0f ) return 0.0f ; for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii] * r[ii] ; return ( ss/sqrtf(rv*xv) ) ; } /*------------------------------------------------------------------------------*/ static int num_quantile = 9 ; void THD_quantile_corr_setup( int nq ) { if( nq > 1 && nq < 100 ) num_quantile = nq ; } float quantile_prepare( int n , float *a ) { int ii ; float rb , rs , jf ; jf = 0.001f + 1.00001f * (n-0.5f) / (float)num_quantile ; if( jf <= 2.0f ) return spearman_rank_prepare(n,a) ; jf = 1.0f / jf ; rank_order_float(n,a) ; /* convert to ranks */ for( rb=0.0f,ii=0 ; ii < n ; ii++ ){ a[ii] = (int)( (a[ii]+0.333f)*jf ) ; rb += a[ii] ; } rb /= n ; for( rs=0.0f,ii=0 ; ii < n ; ii++ ){ a[ii] -= rb ; rs += a[ii]*a[ii] ; } return rs ; } float THD_quantile_corr( int n , float *x , float *y ) { float xv,yv,ss ; int ii ; if( n < 2 ) return 0.0f ; xv = quantile_prepare(n,x) ; if( xv <= 0.0f ) return 0.0f ; yv = quantile_prepare(n,y) ; if( yv <= 0.0f ) return 0.0f ; for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii] * y[ii] ; return ( ss/sqrtf(yv*xv) ) ; } float quantile_corr( int n , float *x , float rv , float *r ) { register int ii ; register float ss ; float xv ; xv = quantile_prepare( n , x ) ; if( xv <= 0.0f ) return 0.0f ; for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii] * r[ii] ; return ( ss/sqrtf(rv*xv) ) ; } /*---------------------------------------------------------------------------*/ #if 1 #undef ttt_bot #undef ttt_top #define ttt_bot 0.3333333f #define ttt_top 0.6666667f void tictactoe_set_thresh( float bb , float tt ){ return; } #else static float ttt_bot = 0.3333333f ; static float ttt_top = 0.6666667f ; void tictactoe_set_thresh( float bb , float tt ) { if( bb >= 0.0f && bb < tt && tt <= 1.0f ){ ttt_bot = bb; ttt_top = tt; } else { ttt_bot = 0.3333333f; ttt_top = 0.6666667f; } } #endif /*---------------------------------------------------------------------------*/ /*! Prepare for tictactoe correlation with a[]. -----------------------------------------------------------------------------*/ float tictactoe_corr_prepare( int n , float *a ) { register int ii ; register float rb , rs , rt ; rank_order_float( n , a ) ; rb = ttt_bot * (n-1) ; rt = ttt_top * (n-1) ; rs = 0.0f ; for( ii=0 ; ii < n ; ii++ ){ if( a[ii] > rt ){ a[ii] = 1.0f ; rs += 1.0f ; } else if( a[ii] < rb ){ a[ii] = -1.0f ; rs += 1.0f ; } else { a[ii] = 0.0f ; } } return rs ; } /*------------------------------------------------------------------------------*/ /*! To do tictactoe correlation of x[] with r[], first do rv = tictactoe_corr_prepare(n,r) ; then corr = tictactoe_corr(n,x,rv,r) ; Note that these 2 routines are destructive (r and x are modified). -------------------------------------------------------------------------------*/ float tictactoe_corr( int n , float *x , float rv , float *r ) { register int ii ; register float ss ; float xv ; xv = tictactoe_corr_prepare( n , x ) ; if( xv <= 0.0f ) return 0.0f ; for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii] * r[ii] ; return ( ss/sqrtf(rv*xv) ) ; } /*------------------------------------------------------------------------------*/ float THD_tictactoe_corr( int n , float *x , float *y ) { float yv , xv ; register float ss ; register int ii ; if( n < 3 ) return 0.0f ; xv = tictactoe_corr_prepare(n,x) ; if( xv <= 0.0f ) return 0.0f ; yv = tictactoe_corr_prepare(n,y) ; if( yv <= 0.0f ) return 0.0f ; for( ii=0,ss=0.0f ; ii < n ; ii++ ) ss += x[ii]*y[ii] ; return( ss / sqrtf(xv*yv) ) ; } /*============================================================================= Compute correlations, destructively (i.e., mangling the input arrays) ===============================================================================*/ /*--------------------------------------------------------------------------*/ /*! Spearman rank-order correlation of x[] and y[] (x and y are modified). */ float THD_spearman_corr( int n , float *x , float *y ) { float xv ; if( n < 2 ) return 0.0f ; xv = spearman_rank_prepare(n,x) ; if( xv <= 0.0f ) return 0.0f ; return spearman_rank_corr( n,y,xv,x ) ; } /*--------------------------------------------------------------------------*/ /*! Spearman rank-order correlation of x[] and y[] (nondestructive). */ float THD_spearman_corr_nd( int n , float *x , float *y ) { float *qx, *qy , cv=0.0f ; qx = (float *)malloc(sizeof(float)*n); memcpy(qx,x,sizeof(float)*n); qy = (float *)malloc(sizeof(float)*n); memcpy(qy,y,sizeof(float)*n); cv = THD_spearman_corr(n,qx,qy) ; free((void *)qy); free((void *)qx); return cv ; } double THD_spearman_corr_dble( int n , double *x , double *y ) { float *qx, *qy , cv=0.0f ; int ii ; qx = (float *)malloc(sizeof(float)*n); qy = (float *)malloc(sizeof(float)*n); for( ii=0 ; ii < n ; ii++ ){ qx[ii] = x[ii] ; qy[ii] = y[ii] ; } cv = THD_spearman_corr(n,qx,qy) ; free((void *)qy); free((void *)qx); return (double)cv ; } /*--------------------------------------------------------------------------*/ /*! Kendall Tau_b (x and y are modified) */ float THD_ktaub_corr( int n , float *x , float *y ) { if( n < 2 ) return 0.0f ; qsort_floatfloat( n , x , y ) ; /* preliminary sorting of x, carrying y */ return kendallNlogN( x , y , n ) ; /* the actual work */ } /*--------------------------------------------------------------*/ /*! Quadrant correlation of x[] and y[] (x and y are modified). */ float THD_quadrant_corr( int n , float *x , float *y ) { float xv ; if( n < 2 ) return 0.0f ; xv = quadrant_corr_prepare(n,x) ; if( xv <= 0.0f ) return 0.0f ; return quadrant_corr( n,y,xv,x ) ; } /*--------------------------------------------------------------------------*/ /*! Quadrant rank-order correlation of x[] and y[] (nondestructive). */ #if 1 float THD_quadrant_corr_nd( int n , float *x , float *y ) { float *z ; register float xm,ym,qc ; register int ii ; z = (float *)malloc(sizeof(float)*n) ; memcpy( z , x , sizeof(float)*n ) ; xm = qmed_float( n , z ) ; z = (float *)malloc(sizeof(float)*n) ; memcpy( z , y , sizeof(float)*n ) ; ym = qmed_float( n , z ) ; free(z) ; qc = 0.0f ; for( ii=0 ; ii < n ; ii++ ) qc += (x[ii] > xm) * (y[ii] > ym) ; qc = (4.0f*qc) / n - 1.0f ; if( qc < -1.0f ) qc = -1.0f; else if( qc > 1.0f ) qc = 1.0f; qc = sinf(1.570796f*qc) ; /* adjust to normal model */ return qc ; } #else /*--------------------------------------------------------------------------*/ float THD_quadrant_corr_nd( int n , float *x , float *y ) { float *qx, *qy , cv=0.0f ; qx = (float *)malloc(sizeof(float)*n); memcpy(qx,x,sizeof(float)*n); qy = (float *)malloc(sizeof(float)*n); memcpy(qy,y,sizeof(float)*n); cv = THD_quadrant_corr(n,qx,qy) ; free((void *)qy); free((void *)qx); return cv ; } #endif /*--------------------------------------------------------------------------*/ /*! Pearson correlation of x[] and y[] (x and y are NOT modified. */ float THD_pearson_corr( int n, float *x , float *y ) { float xv=0.0f , yv=0.0f , xy=0.0f , vv,ww ; float xm=0.0f , ym=0.0f ; register int ii ; if( n < 2 ) return 0.0f ; for( ii=0 ; ii < n ; ii++ ){ xm += x[ii] ; ym += y[ii] ; } xm /= n ; ym /= n ; for( ii=0 ; ii < n ; ii++ ){ vv = x[ii]-xm; ww = y[ii]-ym; xv += vv*vv; yv += ww*ww; xy += vv*ww; } if( xv <= 0.0f || yv <= 0.0f ) return 0.0f ; return xy/sqrtf(xv*yv) ; } /*--------------------------------------------------------------------------*/ /*! Double Pearson correlation of x[] and y[] (x and y are NOT modified) */ double THD_pearson_corrd( int n, double *x , double *y ) { double xv=0.0 , yv=0.0 , xy=0.0 , vv,ww ; double xm=0.0 , ym=0.0 ; register int ii ; if( n < 2 ) return 0.0 ; for( ii=0 ; ii < n ; ii++ ){ xm += x[ii] ; ym += y[ii] ; } xm /= n ; ym /= n ; for( ii=0 ; ii < n ; ii++ ){ vv = x[ii]-xm; ww = y[ii]-ym; xv += vv*vv; yv += ww*ww; xy += vv*ww; } if( xv <= 0.0 || yv <= 0.0 ) return 0.0 ; return xy/sqrt(xv*yv) ; } /*-----------------------------------------------------------------------*/ /* Partial Correlation of A & B, adjusting for C. */ float THD_pearson_partial_corr( int n, float *a, float *b, float *c ) { float AB, AC, BC; //calculate correlation AB AB = THD_pearson_corr(n, a, b); //calculate correlation AC AC = THD_pearson_corr(n, a, c); //calculate correlation BC BC = THD_pearson_corr(n, b, c); return (AB - (AC*BC)) / sqrt( (1-(AC*AC)) * (1-(BC*BC)) ); } /*-----------------------------------------------------------------------*/ /* Double Partial Correlation of A & B, adjusting for C. */ double THD_pearson_partial_corrd( int n, double *a, double *b, double *c ) { double AB, AC, BC; //calculate correlation AB AB = THD_pearson_corrd(n, a, b); //calculate correlation AC AC = THD_pearson_corrd(n, a, c); //calculate correlation BC BC = THD_pearson_corrd(n, b, c); return (AB - (AC*BC)) / sqrt( (1-(AC*AC)) * (1-(BC*BC)) ); } /*--------------------------------------------------------------------------*/ /*! Covariance of x[] and y[] (x and y are NOT modified). ZSS May 18 2012*/ float THD_covariance( int n, float *x , float *y ) { float xy=0.0f , vv,ww ; float xm=0.0f , ym=0.0f ; register int ii ; if( n < 2 ) return 0.0f ; for( ii=0 ; ii < n ; ii++ ){ xm += x[ii] ; ym += y[ii] ; } xm /= n ; ym /= n ; for( ii=0 ; ii < n ; ii++ ){ xy += (x[ii]-xm)*(y[ii]-ym); } return xy/(float)(n-1) ; } /*-------------------------------------------------------------------------*/ /*! Returns a float_triple with (a,b,r) where y = a*x + b is the L2 regression result and r = Pearson correlation coeff. For bootstrapping, ix[i] is the i-th index in x[] and y[] to use. For non-bootstrapping, pass in ix==NULL. *//*-----------------------------------------------------------------------*/ #undef IX #define IX(i) ( ((ix) == NULL) ? (i) : ix[(i)] ) float_triple THD_pearson_indexed( int nix, int *ix, float *x, float *y ) { float xbar=0,ybar=0, xq=0,yq=0,xyq=0, a=0,b=0,r=0; int jj,ii; float_triple abr; for( jj=0 ; jj < nix ; jj++ ){ ii = IX(jj); xbar += x[ii]; ybar += y[ii]; } xbar /= nix ; ybar /= nix ; for( jj=0 ; jj < nix ; jj++ ){ ii = IX(jj) ; xq += (x[ii]-xbar)*(x[ii]-xbar) ; yq += (y[ii]-ybar)*(y[ii]-ybar) ; xyq += (x[ii]-xbar)*(y[ii]-ybar) ; } if( xq > 0.0f && yq > 0.0f ){ r = xyq/sqrtf(xq*yq); a = xyq/xq; b = (xq*ybar-xbar*xyq)/xq; } abr.a = a ; abr.b = b ; abr.c = r ; return abr ; } #undef IX /*---------------------------------------------------------------------------*/ /* Correlates and also returns 2.5%..97.5% confidence interval, via bootstrap. rrr = correlation coefficient, 2.5% level, 97.5% level [in that order] aaa = regression 'a' coefficient, in y=ax+b [as .a .b .c ] bbb = regression 'b' coefficient [components ] *//*-------------------------------------------------------------------------*/ #undef NBOOT #undef NB5 #define NBOOT 960 #define NB5 24 /* must be 2.5% of the above */ void THD_pearson_corr_boot( int n , float *x , float *y , float_triple *rrr , float_triple *aaa , float_triple *bbb ) { int ii,kk ; float aa,bb,rr ; float_triple abr ; int *ix ; float ax[NBOOT] , bx[NBOOT] , rx[NBOOT] ; if( n < 5 || x == NULL || y == NULL ) return ; if( rrr == NULL && aaa == NULL && bbb == NULL ) return ; /* compute standard results */ abr = THD_pearson_indexed( n , NULL , x, y ) ; aa = abr.a ; bb = abr.b ; rr = abr.c ; /* non-bootstrapped answers */ ix = (int *)malloc(sizeof(int)*n) ; /* compute bootstrap results */ for( kk=0 ; kk < NBOOT ; kk++ ){ for( ii=0 ; ii < n ; ii++ ) ix[ii] = lrand48() % n ; abr = THD_pearson_indexed( n , ix , x, y ) ; ax[kk] = abr.a ; bx[kk] = abr.b ; rx[kk] = abr.c ; } free(ix) ; /* sort, then find 2.5% and 97.5% points, save into output structs */ /* [bootstrap confidence intervals are now bias-corrected!] */ if( rrr != NULL ){ float_triple qqq = THD_bootstrap_confinv( rr , 0.05f , NBOOT , rx ) ; rrr->a = rr ; /* would use qqq.b for bias-corrected estimate */ rrr->b = qqq.a ; /* lower edge of confidence interval */ rrr->c = qqq.c ; /* upper edge */ } if( aaa != NULL ){ float_triple qqq = THD_bootstrap_confinv( aa , 0.05f , NBOOT , ax ) ; aaa->a = aa ; aaa->b = qqq.a ; aaa->c = qqq.c ; } if( bbb != NULL ){ float_triple qqq = THD_bootstrap_confinv( bb , 0.05f , NBOOT , bx ) ; bbb->a = bb ; bbb->b = qqq.a ; bbb->c = qqq.c ; } return ; } /*---------------------------------------------------------------------------*/ float THD_spearman_indexed( int nix, int *ix, float *x, float *y ) { float *xt , *yt ; float r ; xt = (float *)malloc(sizeof(float)*nix) ; yt = (float *)malloc(sizeof(float)*nix) ; if( ix == NULL ){ memcpy(xt,x,sizeof(float)*nix) ; memcpy(yt,y,sizeof(float)*nix) ; } else { int ii , jj ; for( jj=0 ; jj < nix ; jj++ ){ ii = ix[jj] ; xt[jj] = x[ii] ; yt[jj] = y[ii] ; } } r = THD_spearman_corr( nix , xt , yt ) ; free(yt) ; free(xt) ; return r ; } /*------------------------------------------------------------------------------*/ float_pair THD_l1_fit_to_line( int n , float *x , float *y ) { float_pair ab={0.0f,0.0f} ; float *A[2] , coef[2] , val ; int ii ; if( n < 3 || x == NULL || y == NULL ) return ab ; A[0] = x ; A[1] = (float *)malloc(sizeof(float)*n) ; for( ii=0 ; ii < n ; ii++ ) A[1][ii] = 1.0f ; val = cl1_solve( n , 2 , y , A , coef , 0 ) ; free(A[1]) ; if( val >= 0.0f ){ ab.a = coef[0] ; ab.b = coef[1] ; } return ab ; } /*------------------------------------------------------------------------------*/ void THD_spearman_corr_boot( int n , float *x , float *y , float_triple *rrr ) { int ii,kk ; float rr ; float_triple qqq ; int *ix ; float rx[NBOOT] ; ENTRY("THD_spearman_corr_boot") ; if( n < 5 || x == NULL || y == NULL ) EXRETURN ; if( rrr == NULL ) EXRETURN ; /* compute standard result */ rr = THD_spearman_indexed( n , NULL , x , y ) ; /* non-bootstrapped answer */ ix = (int *)malloc(sizeof(int)*n) ; /* compute bootstrap results */ for( kk=0 ; kk < NBOOT ; kk++ ){ for( ii=0 ; ii < n ; ii++ ) ix[ii] = lrand48() % n ; rx[kk] = THD_spearman_indexed( n , ix , x , y ) ; } free(ix) ; /* sort, then find 2.5% and 97.5% points, save into output structs */ /* [bootstrap confidence intervals are now bias-corrected!] */ qqq = THD_bootstrap_confinv( rr , 0.05f , NBOOT , rx ) ; rrr->a = rr ; /* would use qqq.b for bias-corrected estimate */ rrr->b = qqq.a ; /* lower edge of confidence interval */ rrr->c = qqq.c ; /* upper edge */ EXRETURN ; } /*----------------------------------------------------------------*/ float THD_pearson_corr_wt( int n, float *x , float *y , float *wt ) { float xv=0.0f , yv=0.0f , xy=0.0f , vv,ww ; float xm=0.0f , ym=0.0f , ws=0.0f ; int ii ; if( wt == NULL ) return THD_pearson_corr(n,x,y) ; for( ii=0 ; ii < n ; ii++ ){ xm += wt[ii]*x[ii]; ym += wt[ii]*y[ii]; ws += wt[ii]; } xm /= ws ; ym /= ws ; for( ii=0 ; ii < n ; ii++ ){ vv = x[ii]-xm ; ww = y[ii]-ym ; xv += wt[ii]*vv*vv ; yv += wt[ii]*ww*ww ; xy += wt[ii]*vv*ww ; } #if 0 ININFO_message("THD_pearson_corr_wt: n=%d ws=%g xm=%g ym=%g xv=%g yv=%g xy=%g", n,ws,xm,ym,xv,yv,xy ) ; #endif if( xv <= 0.0f || yv <= 0.0f ) return 0.0f ; return xy/sqrtf(xv*yv) ; } /*----------------------------------------------------------------*/ // [PT: Nov 9, 2022] weighted Pearson corr--double time! double THD_pearson_corrd_wt( int n, double *x , double *y , float *wt ) { double xv=0.0 , yv=0.0 , xy=0.0 , vv,ww ; double xm=0.0 , ym=0.0 , ws=0.0 ; int ii ; if( wt == NULL ) return THD_pearson_corrd(n,x,y) ; for( ii=0 ; ii < n ; ii++ ){ xm += wt[ii]*x[ii]; ym += wt[ii]*y[ii]; ws += wt[ii]; } xm /= ws ; ym /= ws ; for( ii=0 ; ii < n ; ii++ ){ vv = x[ii]-xm ; ww = y[ii]-ym ; xv += wt[ii]*vv*vv ; yv += wt[ii]*ww*ww ; xy += wt[ii]*vv*ww ; } #if 0 ININFO_message("THD_pearson_corrd_wt: n=%d ws=%g xm=%g ym=%g xv=%g yv=%g xy=%g", n,ws,xm,ym,xv,yv,xy ) ; #endif if( xv <= 0.0f || yv <= 0.0f ) return 0.0f ; return xy/sqrt(xv*yv) ; } /*----------------------------------------------------------------*/ /* compute distance between two arrays ZSS May 04 2012 */ float THD_distance( int n, float *x , float *y, int abs) { float dp=0.0, a1, a2; int ii, n1 = n-1; for( ii=0 ; ii < n1 ; ii+=2 ) { a1 = x[ii]-y[ii]; a2 = x[ii+1]-y[ii+1]; if (!abs) dp += (a1*a1+a2*a2); else dp += (ABS(a1)+ABS(a2)); } if( ii == n1 ) { a1 = x[ii]-y[ii]; if (!abs) dp += a1*a1; else dp += ABS(a1); } if (!abs) dp = sqrt(dp) ; return (dp); } /*--------------------------------------------------------------------------*/ /*! Compute the rank-order correlation between 2 images [08 Mar 2006]. ----------------------------------------------------------------------------*/ float mri_spearman_corr( MRI_IMAGE *im , MRI_IMAGE *jm ) { float *far , *gar , cc ; MRI_IMAGE *fim , *gim ; if( im == NULL || jm == NULL || im->nvox != jm->nvox ) return 0.0f ; fim = mri_to_float(im) ; far = mri_data_pointer(fim) ; gim = mri_to_float(jm) ; gar = mri_data_pointer(gim) ; cc = THD_spearman_corr( fim->nvox , far , gar ) ; mri_free(gim) ; mri_free(fim) ; return cc ; } /*----------------------------------------------------------------*/ /*! eta^2 (Cohen, NeuroImage 2008) 25 Jun 2010 [rickr] * * eta^2 = 1 - SUM[ (a_i - m_i)^2 + (b_i - m_i)^2 ] * ------------------------------------ * SUM[ (a_i - M )^2 + (b_i - M )^2 ] * * where o a_i and b_i are the vector elements * o m_i = (a_i + b_i)/2 * o M = mean across both vectors -----------------------------------------------------------------*/ float THD_eta_squared( int n, float *x , float *y ) { double num=0.0f , denom = 0.0f ; float gm=0.0f , lm, vv, ww; int ii ; for( ii=0 ; ii < n ; ii++ ){ gm += x[ii] + y[ii] ; } gm /= (2*n) ; for( ii=0 ; ii < n ; ii++ ){ lm = 0.5 * ( x[ii] + y[ii] ) ; vv = (x[ii]-lm); ww = (y[ii]-lm); num += ( vv*vv + ww*ww ); vv = (x[ii]-gm); ww = (y[ii]-gm); denom += ( vv*vv + ww*ww ); } if( num < 0.0f || denom <= 0.0f || num >= denom ) return 0.0f ; return 1.0 - num/denom ; } /*----------------------------------------------------------------*/ /*! eta^2 (Cohen, NeuroImage 2008) 16 Jun 2011 [rickr] * * Same as THD_eta_squared, but allow a mask and return a double. * * eta^2 = 1 - SUM[ (a_i - m_i)^2 + (b_i - m_i)^2 ] * ------------------------------------ * SUM[ (a_i - M )^2 + (b_i - M )^2 ] * * where o a_i and b_i are the vector elements * o m_i = (a_i + b_i)/2 * o M = mean across both vectors -----------------------------------------------------------------*/ double THD_eta_squared_masked( int n, float *x , float *y, byte *mask ) { double num=0.0f , denom = 0.0f ; float gm=0.0f , lm, vv, ww; int ii, nm ; for( ii=0, nm=0 ; ii < n ; ii++ ) if( !mask || mask[ii] ) { gm += x[ii] + y[ii] ; nm++ ; } if( nm == 0 ) return 0.0 ; /* bail on empty mask */ gm /= (2*nm) ; for( ii=0 ; ii < n ; ii++ ){ if( mask && !mask[ii] ) continue; /* skip any not in mask */ lm = 0.5f * ( x[ii] + y[ii] ) ; vv = (x[ii]-lm); ww = (y[ii]-lm); num += ( vv*vv + ww*ww ); vv = (x[ii]-gm); ww = (y[ii]-gm); denom += ( vv*vv + ww*ww ); } if( num < 0.0 || denom <= 0.0 || num >= denom ) return 0.0 ; return 1.0 - num/denom ; } /*----------------------------------------------------------------*/ /*! dice coefficient (allow for mask) 27 Oct 2015 [rickr] * (float input) * * dice = 2 * size(intersection) / ( size(A) + size(B) ) * * Simply count voxels for the 3 cases. -----------------------------------------------------------------*/ float THD_dice_coef_f_masked( int n, float *x , float *y, byte *mask ) { int nA=0, nB=0, nAB=0; int ii; for( ii=0 ; ii < n ; ii++ ) { if( mask && !mask[ii] ) continue; if( x[ii] ) { nA++; if( y[ii] ) { nB++; nAB++; } } else if ( y[ii] ) nB++; } if( (nA+nB) > 0 ) return 2.0*nAB/(nA+nB); return 0.0; } /****************************************************************************/ /*** Histogram-based measurements of dependence between two float arrays. ***/ /****************************************************************************/ /*--------------------------------------------------------------------------*/ /* Extensive changes from the olde method: * values below bot and above top are not used * histogram can have bot and top bins unequal, and in-between ones equal (if use_xyclip != 0) * histogram can have all unequal bin sizes (if nxybin > 2) ----------------------------------------------------------------------------*/ #define LINHIST /* do linear spread in histogram */ #undef WW #define WW(i) ((w==NULL) ? 1.0f : w[i]) /* weight function for i'th datum */ static float *xc=NULL , *yc=NULL , *xyc=NULL , nww=0.0f ; static int nbin=0 , nbp=0 , nbm=0 ; static float *xbin=NULL , *ybin=NULL ; static int nxybin=0 ; static int use_xyclip=0 ; static float xclip_bot, xclip_top ; static float yclip_bot, yclip_top ; #undef XYC #define XYC(p,q) xyc[(p)+(q)*nbp] #ifndef WAY_BIG # define WAY_BIG 1.e+10 #endif #undef GOODVAL #define GOODVAL(x) ((x) < WAY_BIG) /* x is not preposterous */ #undef RANGVAL #define RANGVAL(x,b,t) ((x) >= (b) && (x) <= (t)) /* x between b and t */ #undef FREEIF #define FREEIF(x) \ do{ if((x)!=NULL){ free((void *)(x)); (x)=NULL; } } while(0) /*--------------------------------------------------------------------------*/ /*! Clear the internal 2D histogram (but not the bin settings). ----------------------------------------------------------------------------*/ void clear_2Dhist(void) { FREEIF(xc); FREEIF(yc); FREEIF(xyc); nbin = nbp = nbm = 0; nww = 0.0f; return; } /*--------------------------------------------------------------------------*/ static double hpow = 0.33333333333 ; void set_2Dhist_hpower( double hh ) { hpow = (hh > 0.0 && hh < 1.0) ? hh : 0.33333333333 ; clear_2Dhist() ; } /*--------------------------------------------------------------------------*/ static int nhbin = 0 ; void set_2Dhist_hbin( int nn ){ nhbin = nn; clear_2Dhist(); } int get_2Dhist_hbin( void ){ return nhbin ; } /*--------------------------------------------------------------------------*/ /*! Retrieve the 2D histogram built previously in build_2Dhist(). - Return value is the number of bins in each direction (may be 0). - *xyhist is points to internal 2D array (may be NULL). ----------------------------------------------------------------------------*/ int retrieve_2Dhist( float **xyhist ) { if( xyhist == NULL ) return 0 ; *xyhist = xyc ; return nbp ; } /*--------------------------------------------------------------------------*/ int retrieve_2Dhist1( float **xh , float **yh ) { if( xh == NULL || yh == NULL ) return 0 ; *xh = xc ; *yh = yc ; return nbp ; } /*--------------------------------------------------------------------------*/ /*! Set the bins to be used in 2D histogram-izing in build_2Dhist(). - nb = number of bins; if nb < 3, equal sized bins will be used and the internal bin arrays will be cleared - xb = [0..nb] (nb+1 elements) = boundary pts for x-axis bins - yb = [0..nb] (nb+1 elements) = boundary pts for y-axis bins - Both xb[] and yb[] must be strictly increasing arrays! - Also see get_2Dhist_xybin() ----------------------------------------------------------------------------*/ void set_2Dhist_xybin( int nb , float *xb , float *yb ) { int ii ; FREEIF(xbin) ; FREEIF(ybin) ; nxybin = 0 ; if( nb > 2 && xb != NULL && yb != NULL ){ for( ii=1 ; ii <= nb ; ii++ ) /* check that bin sizes are OK */ if( xb[ii-1] >= xb[ii] || yb[ii-1] > yb[ii] ) break ; if( ii > nb ){ nxybin = nb ; xbin = (float *)malloc(sizeof(float)*(nb+1)) ; ybin = (float *)malloc(sizeof(float)*(nb+1)) ; memcpy( xbin , xb , sizeof(float)*(nb+1) ) ; memcpy( ybin , yb , sizeof(float)*(nb+1) ) ; } else { WARNING_message("set_2Dhist_xybin: illegal inputs!") ; } } return ; } /*--------------------------------------------------------------------------*/ void set_2Dhist_xybin_eqwide( int nb, float xbot,float xtop,float ybot,float ytop ) { int ii ; float dx , dy ; FREEIF(xbin) ; FREEIF(ybin) ; nxybin = 0 ; if( nb > 2 && xbot < xtop && ybot < ytop ){ nxybin = nb ; xbin = (float *)malloc(sizeof(float)*(nb+1)) ; ybin = (float *)malloc(sizeof(float)*(nb+1)) ; dx = (xtop-xbot) / nb ; dy = (ytop-ybot) / nb ; for( ii=0 ; ii < nb ; ii++ ){ xbin[ii] = xbot + ii*dx ; ybin[ii] = ybot + ii*dy ; } xbin[nb] = xtop ; ybin[nb] = ytop ; #if 0 INFO_message("set_2Dhist_xybin_eqwide: %d %f..%f %f..%f", nb,xbot,xtop,ybot,ytop ) ; #endif } return ; } /*--------------------------------------------------------------------------*/ static float_pair clipate( int nval , float *xar ) { MRI_IMAGE *qim; float cbot,ctop, mmm , *qar; float_pair rr; int ii,nq; ENTRY("clipate") ; qim = mri_new_vol( nval,1,1 , MRI_float ) ; qar = MRI_FLOAT_PTR(qim) ; for( ii=nq=0 ; ii < nval ; ii++ ) if( GOODVAL(xar[ii]) ) qar[nq++] = xar[ii]; qim->nx = qim->nvox = nq ; if( nq < 666 ){ rr.a = 1.0f; rr.b = 0.0f; mri_free(qim); RETURN(rr); } mmm = mri_min( qim ) ; if( mmm >= 0.0f ){ /* for positive images */ cbot = THD_cliplevel( qim , 0.345f ) ; ctop = mri_quantile ( qim , 0.966f ) ; if( ctop > 4.321f*cbot ) ctop = 4.321f*cbot ; } else { /* for images including negative values: no go */ cbot = 1.0f; ctop = 0.0f; } mri_free(qim) ; rr.a = cbot ; rr.b = ctop ; RETURN(rr) ; } /*--------------------------------------------------------------------------*/ void set_2Dhist_xyclip( int nval , float *xval , float *yval ) { float_pair xcc , ycc ; ENTRY("set_2Dhist_xyclip") ; use_xyclip = 0 ; if( nval < 666 || xval == NULL || yval == NULL ) EXRETURN ; xcc = clipate( nval , xval ) ; ycc = clipate( nval , yval ) ; if( xcc.a >= xcc.b || ycc.a >= ycc.b ) EXRETURN ; use_xyclip = 1 ; nxybin = 0 ; xclip_bot = xcc.a ; xclip_top = xcc.b ; yclip_bot = ycc.a ; yclip_top = ycc.b ; EXRETURN ; } /*--------------------------------------------------------------------------*/ int get_2Dhist_xyclip( float *xbc , float *xtc , float *ybc , float *ytc ) { *xbc = xclip_bot ; *xtc = xclip_top ; *ybc = yclip_bot ; *ytc = yclip_top ; return use_xyclip ; } /*--------------------------------------------------------------------------*/ static int eqhighate( int nb , int nval , float *xar , float *xb ) { float *xd , xbot,xtop , frac,fi , xtest ; int ii, pp, pbot,ptop , bsiz , pstart,pend , nxd ; xd = (float *)malloc(sizeof(float)*nval) ; for( ii=nxd=0 ; ii < nval ; ii++ ) if( GOODVAL(xar[ii]) ) xd[nxd++] = xar[ii]; if( nxd < 7*nb ){ free(xd); return 0; } /* bad */ qsort_float( nxd , xd ) ; /** scan for plateaus = runs of constant values longer than nominal bin width, at the bottom and top **/ xbot = xb[0] = xd[0] ; xtop = xb[nb] = xd[nxd-1] ; if( xtop <= xbot ){ free(xd); return 0; } bsiz = (nxd/nb) ; xtest = xbot + (xtop-xbot)/(100.0f*nb) ; for( pp=1 ; pp < nxd && xd[pp] < xtest ; pp++ ) ; /*nada*/ if( pp == nxd ){ free(xd); return 0; } /* data is constant? */ pbot = (pp > bsiz) ? pp : 0 ; xtest = xtop - (xtop-xbot)/(100.0f*nb) ; for( pp=nxd-2 ; pp > 0 && xd[pp] > xtest ; pp-- ) ; /*nada*/ if( pp <= pbot ){ free(xd); return 0; } /* something screwy */ ptop = (nxd-1-pp > bsiz) ? pp : nxd-1 ; if( pbot > 0 ){ xb[1] = 0.999999f*xd[pbot] + 0.000001f*xbot ; pstart = 2 ; } else { pstart = 1 ; } if( ptop < nxd-1 ){ xb[nb-1] = 0.999999f*xd[ptop] + 0.000001f*xtop ; pend = nb-2 ; } else { pend = nb-1 ; } frac = (ptop-pbot)/(pend-pstart+2.0f) ; for( pp=pstart ; pp <= pend ; pp++ ){ fi = pbot + frac*(pp-pstart+1.0f) ; ii = (int)fi ; fi = fi - ii ; xb[pp] = (1.0f-fi) * xd[ii] + fi * xd[ii+1] ; } free(xd) ; return nb ; } /*--------------------------------------------------------------------------*/ void set_2Dhist_xybin_eqhigh( int nb, int nval, float *xval, float *yval ) { int ii,jj ; FREEIF(xbin) ; FREEIF(ybin) ; nxybin = 0 ; if( nb < 3 || nval < 9*nb || xval == NULL || yval == NULL ) return ; nxybin = nb ; xbin = (float *)malloc(sizeof(float)*(nb+1)) ; ybin = (float *)malloc(sizeof(float)*(nb+1)) ; ii = eqhighate( nb , nval , xval , xbin ) ; jj = eqhighate( nb , nval , yval , ybin ) ; if( ii == 0 || jj == 0 ){ FREEIF(xbin) ; FREEIF(ybin) ; nxybin = 0 ; /* bad things happened */ } return ; } /*--------------------------------------------------------------------------*/ /*! Sets pointers to the bins in 2D histogram-ization, and return value is the number of bins (if 0, equal size bins are in use). Do NOT screw with the contents of the pointers: these are internal arrays! Also see set_2Dhist_xybin(). ----------------------------------------------------------------------------*/ int get_2Dhist_xybin( float **xb , float **yb ) { if( xb != NULL ) *xb = xbin ; if( yb != NULL ) *yb = ybin ; return nxybin ; } /*--------------------------------------------------------------------------*/ /*! Load 2D histogram of x[0..n-1] and y[0..n-1], each point optionally weighted by w[0..n-1] (weights are all 1 if w==NULL). Used in the histogram-based measures of dependence between x[] and y[i]. If something is bad on input, nbin is set to 0. Otherwise, these global variables are set: - nbin = # of bins, nbp = nbin+1 - nww = sum of the weights used - xc = marginal histogram of x[], for xc[0..nbin] (nbp points in) - yc = marginal histogram of y[], for yc[0..nbin] (each direction) - xyc = joint histogram of (x[],y[]), for XYC(0..nbin,0..nbin) - The histograms are normalized (by 1/nww) to have sum==1 - Histogram can be retrieved by retrieve_2Dhist() and can be erased by clear_2Dhist() - Default number of equal-spaced bins in each direction is n^(1/3) - the exponent can be changed with set_2Dhist_hpower() - you can set the number of bins with set_2Dhist_hbin() - you can set unequal bins with set_2Dhist_xybin() - x[] values outside the range xbot..xtop (inclusive) or outside the unequal bins set in set_2Dhist_xybin() (if applicable) will not be used in the histogram; mutatis mutandum for y[] ----------------------------------------------------------------------------*/ static void stuff_2Dhist( int n, float xbot,float xtop,float *x, float ybot,float ytop,float *y, float *w, int addon ) { register int ii,jj,kk ; float xb,xi , yb,yi , xx,yy , x1,y1 , ww ; byte *good ; int ngood , xyclip ; static float xxbot=0.0f,xxtop=0.0f , yybot=0.0f,yytop=0.0f ; ENTRY("stuff_2Dhist") ; if( n <= 0 || x == NULL || y == NULL ) EXRETURN ; /* get the min..max range for x and y data? */ good = (byte *)malloc(sizeof(byte)*n) ; /* 28 Feb 2007 */ for( ii=0 ; ii < n ; ii++ ) good[ii] = GOODVAL(x[ii]) && GOODVAL(y[ii]) ; addon = addon && (xxbot < xxtop) && (yybot < yytop) && (nww > 0.0f) ; if( addon ){ xbot = xxbot ; ybot = yybot ; xtop = xxtop ; ytop = yytop ; } else { clear_2Dhist() ; if( nxybin > 2 ){ /* unequal bins */ xbot = xbin[0] ; xtop = xbin[nxybin] ; } else if( xbot >= xtop ){ /* equal bins, and must find range */ xbot = WAY_BIG ; xtop = -WAY_BIG ; for( ii=0 ; ii < n ; ii++ ) if( good[ii] ){ if( x[ii] > xtop ) xtop = x[ii] ; else if( x[ii] < xbot ) xbot = x[ii] ; } } if( xbot >= xtop ){ clear_2Dhist(); free(good); EXRETURN; } if( nxybin > 2 ){ ybot = ybin[0] ; ytop = ybin[nxybin] ; } else if( ybot >= ytop ){ ybot = WAY_BIG ; ytop = -WAY_BIG ; for( ii=0 ; ii < n ; ii++ ) if( good[ii] ){ if( y[ii] > ytop ) ytop = y[ii] ; else if( y[ii] < ybot ) ybot = y[ii] ; } } if( ybot >= ytop ){ free(good); EXRETURN; } xxbot = xbot ; xxtop = xtop ; yybot = ybot ; yytop = ytop ; } /*-- count number of good values left in range (in both x and y) --*/ memset(good,0,n) ; for( ngood=ii=0 ; ii < n ; ii++ ){ if( RANGVAL(x[ii],xbot,xtop) && RANGVAL(y[ii],ybot,ytop) && WW(ii) > 0.0f ){ good[ii] = 1 ; ngood++ ; } } if( ngood == 0 ){ free(good) ; EXRETURN ; } /*--- allocate space for 2D and 1D histograms ---*/ if( !addon ){ if( nxybin > 2 ){ /* array size is fixed by bins */ nbin = nxybin ; } else { /* compute new histo array size */ nbin = (nhbin > 2) ? nhbin : (int)pow((double)n,hpow) ; if( nbin > 255 ) nbin = 255; else if( nbin < 3 ) nbin = 3; } nbp = nbin+1 ; nbm = nbin-1 ; FREEIF(xc) ; FREEIF(yc) ; FREEIF(xyc) ; xc = (float *)calloc(sizeof(float),nbp) ; yc = (float *)calloc(sizeof(float),nbp) ; xyc = (float *)calloc(sizeof(float),nbp*nbp) ; nww = 0.0f ; } /*--------------- make the 2D and 1D histograms ---------------*/ xyclip = (nxybin <= 2) && use_xyclip && (xbot < xclip_bot) && (xclip_bot < xclip_top) && (xclip_top < xtop) && (ybot < yclip_bot) && (yclip_bot < yclip_top) && (yclip_top < ytop) ; if( nxybin <= 0 && !xyclip ){ /*------------ equal size bins ------------*/ #if 0 if(PRINT_TRACING){ char str[256]; sprintf(str,"equal size bins: xbot=%g xtop=%g ybot=%g ytop=%g nbin=%d nval=%d ngood=%d", xbot,xtop,ybot,ytop,nbin,n,ngood); STATUS(str); } #endif xb = xbot ; xi = nbm/(xtop-xbot) ; yb = ybot ; yi = nbm/(ytop-ybot) ; for( ii=0 ; ii < n ; ii++ ){ if( !good[ii] ) continue ; xx = (x[ii]-xb)*xi ; jj = (int)xx ; xx = xx - jj ; x1 = 1.0f-xx ; yy = (y[ii]-yb)*yi ; kk = (int)yy ; yy = yy - kk ; y1 = 1.0f-yy ; ww = WW(ii) ; nww += ww ; #ifdef LINHIST xc[jj] += x1*ww ; xc[jj+1] += xx*ww ; yc[kk] += (y1*ww); yc[kk+1] += (yy*ww); XYC(jj ,kk ) += x1*(y1*ww) ; XYC(jj+1,kk ) += xx*(y1*ww) ; XYC(jj ,kk+1) += x1*(yy*ww) ; XYC(jj+1,kk+1) += xx*(yy*ww) ; #else xc[jj] += ww ; yc[kk] += ww ; XYC(jj,kk) += ww ; #endif } } else if( xyclip ){ /*------------ mostly equal bins ----------------*/ float xbc=xclip_bot , xtc=xclip_top , ybc=yclip_bot , ytc=yclip_top ; #if 0 if(PRINT_TRACING){ char str[256]; sprintf(str,"mostly equal bins: xbc=%g xtc=%g ybc=%g ytc=%g nbin=%d", xbc,xtc,ybc,ytc,nbin) ; STATUS(str); } #endif xi = (nbin-2.000001f)/(xtc-xbc) ; yi = (nbin-2.000001f)/(ytc-ybc) ; for( ii=0 ; ii < n ; ii++ ){ if( !good[ii] ) continue ; xx = x[ii] ; if( xx < xbc ){ jj = 0 ; xx = 0.0f ; } else if( xx > xtc ){ jj = nbm ; xx = 1.0f ; } else { xx = 1.0f+(xx-xbc)*xi; jj = (int)xx; xx = xx - jj; } yy = y[ii] ; if( yy < ybc ){ kk = 0 ; yy = 0.0f ; } else if( yy > ytc ){ kk = nbm ; yy = 1.0f ; } else { yy = 1.0f+(yy-ybc)*yi; kk = (int)yy; yy = yy - kk; } x1 = 1.0f-xx ; y1 = 1.0f-yy ; ww = WW(ii) ; nww += ww ; #ifdef LINHIST xc[jj] += x1*ww ; xc[jj+1] += xx*ww ; yc[kk] += (y1*ww); yc[kk+1] += (yy*ww); XYC(jj ,kk ) += x1*(y1*ww) ; XYC(jj+1,kk ) += xx*(y1*ww) ; XYC(jj ,kk+1) += x1*(yy*ww) ; XYC(jj+1,kk+1) += xx*(yy*ww) ; #else xc[jj] += ww ; yc[kk] += ww ; XYC(jj,kk) += ww ; #endif } } else { /*-------------------- unequal size bins --------------------*/ float *xdup, *ydup, *xv, *yv, *wv, frac,val,xt,yt ; int *xn, *yn, *xin, *yin , ib,nb , maxsort ; maxsort = 32*nxybin ; if( maxsort > 2048 ) maxsort = 2048 ; xdup = (float *)malloc(sizeof(float)*maxsort) ; ydup = (float *)malloc(sizeof(float)*maxsort) ; xv = (float *)malloc(sizeof(float)*maxsort) ; yv = (float *)malloc(sizeof(float)*maxsort) ; wv = (float *)malloc(sizeof(float)*maxsort) ; xn = (int *) malloc(sizeof(float)*maxsort) ; yn = (int *) malloc(sizeof(float)*maxsort) ; xin = (int *) malloc(sizeof(float)*maxsort) ; yin = (int *) malloc(sizeof(float)*maxsort) ; /* not yang */ /* outer loop: process points starting at index = ib */ for( ib=0 ; ib < n ; ){ /* extract up to maxsort good points from the data arrays */ for( nb=0,ii=ib ; ii < n && nb < maxsort ; ii++ ){ if( good[ii] ){ wv[nb]=WW(ii); xdup[nb]=x[ii]; ydup[nb]=y[ii]; xn[nb]=yn[nb]=nb; nb++; /* separate for clarity */ } } #if 0 INFO_message("extracted %d data points up to index=%d",nb,ii) ; #endif ib = ii ; /* where to start extracting next outer loop */ if( nb == 0 ) break ; /* didn't find any good data ==> done! */ /* sort the extracted x data values in xdup[], and keep track of where they came from in xn[] */ qsort_floatint(nb,xdup,xn) ; /* scan through sorted xdup[] values, which will go into bin #kk which is xb..xt; store the bin index and the fractional location within the bin; the reason for doing the sorting is that finding the bin index kk will be efficient (don't have to search from start as we would for unordered values) */ xb = xbin[0] ; xt = xbin[1] ; kk=0 ; frac = 1.0f/(xt-xb) ; for( ii=0 ; ii < nb ; ii++ ){ val = xdup[ii] ; if( val > xt ){ /* not in the xb..xt bin, so scan up until it is */ for( kk++ ; kk+1 < nxybin && val > xbin[kk+1] ; kk++ ) ; /*nada*/ xb = xbin[kk] ; xt = xbin[kk+1] ; frac = 1.0f/(xt-xb) ; } jj = xn[ii] ; /* index in unsorted array */ xv[jj] = frac*(val-xb) ; /* fractional position in bin */ xin[jj] = kk ; /* bin index */ } /* repeat the above for the y arrays */ qsort_floatint(nb,ydup,yn) ; yb = ybin[0] ; yt = ybin[1] ; kk=0 ; frac = 1.0f/(yt-yb) ; for( ii=0 ; ii < nb ; ii++ ){ val = ydup[ii] ; if( val > yt ){ /* not in the yb..yt bin, so scan up until it is */ for( kk++ ; kk+1 < nxybin && val > ybin[kk+1] ; kk++ ) ; /*nada*/ yb = ybin[kk] ; yt = ybin[kk+1] ; frac = 1.0f/(yt-yb) ; } jj = yn[ii] ; /* index in unsorted array */ yv[jj] = frac*(val-yb) ; /* fractional position in bin */ yin[jj] = kk ; /* bin index */ } /* now distribute the values into the 1D and 2D histograms */ for( ii=0 ; ii < nb ; ii++ ){ ww = wv[ii] ; nww += ww ; jj = xin[ii]; kk = yin[ii] ; xx = xv[ii] ; x1 = 1.0f-xx ; yy = yv[ii] ; y1 = 1.0f-yy ; #ifdef LINHIST xc[jj] += x1*ww ; xc[jj+1] += xx*ww ; yc[kk] += (y1*ww); yc[kk+1] += (yy*ww); XYC(jj ,kk ) += x1*(y1*ww) ; XYC(jj+1,kk ) += xx*(y1*ww) ; XYC(jj ,kk+1) += x1*(yy*ww) ; XYC(jj+1,kk+1) += xx*(yy*ww) ; #else xc[jj] += ww ; yc[kk] += ww ; XYC(jj,kk) += ww ; #endif } } /* end of outer loop (ib) over blocks */ free(yin); free(xin); free(yn); free(xn); free(wv); free(yv); free(xv); free(ydup); free(xdup); } /*----- end of test on equal or unequal size bins -----*/ free(good); EXRETURN; } /*............................................................................*/ void build_2Dhist( int n , float xbot,float xtop,float *x , float ybot,float ytop,float *y , float *w ) { if( n > 9 && x != NULL && y != NULL ) stuff_2Dhist( n , xbot,xtop,x , ybot,ytop,y , w , 0 ) ; return ; } /*............................................................................*/ void addto_2Dhist( int n , float xbot,float xtop,float *x , float ybot,float ytop,float *y , float *w ) { if( n > 0 && x != NULL && y != NULL ) stuff_2Dhist( n , xbot,xtop,x , ybot,ytop,y , w , 1 ) ; return ; } /*............................................................................*/ /* scale histogram to have sum==1 */ void normalize_2Dhist(void) { if( nww > 0.0f && xyc != NULL && xc != NULL && yc != NULL ){ register float ni ; register int nbq , ii ; ni = 1.0f / nww ; for( ii=0 ; ii < nbp ; ii++ ){ xc[ii] *= ni; yc[ii] *= ni; } nbq = nbp*nbp ; for( ii=0 ; ii < nbq ; ii++ ){ xyc[ii] *= ni; } } return ; } /*--------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------*/ MRI_IMAGE *build_byteized_vectors( int n , float xbot,float xtop,float *x , float ybot,float ytop,float *y ) { register int ii,jj,kk , gg ; float xb,xi , yb,yi , xx,yy , x1,y1 ; byte *good ; int ngood , xyclip , nbm ; MRI_IMAGE *bim ; byte *xbar,*ybar ; ENTRY("build_byteized_vectors") ; /* bad inputs? */ if( n <= 1 || x == NULL || y == NULL ) RETURN(NULL) ; /* get the min..max range for x data? */ STATUS("compute good[]") ; good = (byte *)malloc(sizeof(byte)*n) ; for( ii=0 ; ii < n ; ii++ ) good[ii] = GOODVAL(x[ii]) && GOODVAL(y[ii]) ; if( nxybin > 2 ){ /* unequal bins */ xbot = xbin[0] ; xtop = xbin[nxybin] ; } else if( xbot >= xtop ){ /* equal bins, and must find range */ xbot = WAY_BIG ; xtop = -WAY_BIG ; for( ii=0 ; ii < n ; ii++ ) if( good[ii] ){ if( x[ii] > xtop ) xtop = x[ii] ; else if( x[ii] < xbot ) xbot = x[ii] ; } } if( xbot >= xtop ){ free(good); RETURN(NULL); } /* get the min..max range for y data? */ if( nxybin > 0 ){ ybot = ybin[0] ; ytop = ybin[nxybin] ; } else if( ybot >= ytop ){ ybot = WAY_BIG ; ytop = -WAY_BIG ; for( ii=0 ; ii < n ; ii++ ) if( good[ii] ){ if( y[ii] > ytop ) ytop = y[ii] ; else if( y[ii] < ybot ) ybot = y[ii] ; } } if( ybot >= ytop ){ free(good); RETURN(NULL); } /*--- figure out quantization levels ---*/ if( nxybin > 2 ){ /* array size is fixed by bins */ nbin = nxybin ; } else { /* compute new histo array size */ nbin = (nhbin > 2) ? nhbin : (int)pow((double)n,hpow) ; } if( nbin > 255 ) nbin = 255; else if( nbin < 16 ) nbin = 16 ; nbm = nbin-1 ; /*-- count number of good values left in range (in both x and y) --*/ memset(good,0,n) ; for( ngood=ii=0 ; ii < n ; ii++ ){ if( RANGVAL(x[ii],xbot,xtop) && RANGVAL(y[ii],ybot,ytop) ){ good[ii] = 1 ; ngood++ ; } } if( ngood < 3*nbin ){ free(good); RETURN(NULL); } /*--------------- make the output data -----------------*/ bim = mri_new( ngood , 2 , MRI_byte ) ; /* zero filled */ xbar = MRI_BYTE_PTR(bim) ; ybar = xbar + ngood ; xyclip = (nxybin <= 0) && use_xyclip && (xbot < xclip_bot) && (xclip_bot < xclip_top) && (xclip_top < xtop) && (ybot < yclip_bot) && (yclip_bot < yclip_top) && (yclip_top < ytop) ; if( nxybin <= 0 && !xyclip ){ /*------------ equal size bins ------------*/ xb = xbot ; xi = (nbin-0.01f)/(xtop-xbot) ; yb = ybot ; yi = (nbin-0.01f)/(ytop-ybot) ; for( gg=ii=0 ; ii < n ; ii++ ){ if( !good[ii] ) continue ; xbar[gg] = (byte)((x[ii]-xb)*xi) ; ybar[gg] = (byte)((y[ii]-yb)*yi) ; gg++ ; } } else if( xyclip ){ /*------------ mostly equal bins ----------------*/ float xbc=xclip_bot , xtc=xclip_top , ybc=yclip_bot , ytc=yclip_top ; xi = (nbin-2.01f)/(xtc-xbc) ; yi = (nbin-2.01f)/(ytc-ybc) ; for( gg=ii=0 ; ii < n ; ii++ ){ if( !good[ii] ) continue ; xx = x[ii] ; if( xx < xbc ){ xbar[gg] = 0 ; } else if( xx > xtc ){ xbar[gg] = nbm ; } else { xbar[gg] = (byte)(1.0f+(xx-xbc)*xi) ; } yy = y[ii] ; if( yy < ybc ){ ybar[gg] = 0 ; } else if( yy > ytc ){ ybar[gg] = nbm ; } else { ybar[gg] = (byte)(1.0f+(yy-ybc)*yi) ; } gg++ ; } } else { /*-------------------- unequal size bins --------------------*/ float *xdup, *ydup, val,xt,yt ; int *xn, *yn, *xin, *yin , ib,nb , maxsort ; maxsort = 16*nxybin ; if( maxsort > 512 ) maxsort = 512 ; xdup = (float *)malloc(sizeof(float)*maxsort) ; ydup = (float *)malloc(sizeof(float)*maxsort) ; xn = (int *) malloc(sizeof(float)*maxsort) ; yn = (int *) malloc(sizeof(float)*maxsort) ; xin = (int *) malloc(sizeof(float)*maxsort) ; yin = (int *) malloc(sizeof(float)*maxsort) ; /* not yang */ /* outer loop: process points starting at index = ib */ for( gg=ib=0 ; ib < n ; ){ /* extract up to maxsort good points from the data arrays */ for( nb=0,ii=ib ; ii < n && nb < maxsort ; ii++ ){ if( good[ii] ){ xdup[nb]=x[ii]; ydup[nb]=y[ii]; xn[nb]=yn[nb]=nb; nb++; } } ib = ii ; /* where to start extracting next outer loop */ if( nb == 0 ) break ; /* didn't find any good data ==> done done! */ /* sort the extracted x data values in xdup[], and keep track of where they came from in xn[] */ qsort_floatint(nb,xdup,xn) ; /* scan through sorted xdup[] values, which will go into bin #kk which is xb..xt; store the bin index; the reason for doing the sorting is that finding the bin index kk will be efficient (don't have to search from start as we would for unordered values) */ xb = xbin[0] ; xt = xbin[1] ; kk=0 ; for( ii=0 ; ii < nb ; ii++ ){ val = xdup[ii] ; if( val > xt ){ /* not in the xb..xt bin, so scan up until it is */ for( kk++ ; kk+1 < nxybin && val > xbin[kk+1] ; kk++ ) ; /*nada*/ xb = xbin[kk] ; xt = xbin[kk+1] ; } jj = xn[ii] ; /* index in unsorted array */ xin[jj] = kk ; /* bin index */ } /* repeat the above for the y arrays */ qsort_floatint(nb,ydup,yn) ; yb = ybin[0] ; yt = ybin[1] ; kk=0 ; for( ii=0 ; ii < nb ; ii++ ){ val = ydup[ii] ; if( val > yt ){ /* not in the yb..yt bin, so scan up until it is */ for( kk++ ; kk+1 < nxybin && val > ybin[kk+1] ; kk++ ) ; /*nada*/ yb = ybin[kk] ; yt = ybin[kk+1] ; } jj = yn[ii] ; /* index in unsorted array */ yin[jj] = kk ; /* bin index */ } /* now distribute the values into the 1D and 2D histograms */ for( ii=0 ; ii < nb ; ii++ ){ xbar[gg] = (byte)xin[ii] ; ybar[gg] = (byte)yin[ii] ; gg++ ; } } /* end of outer loop (ib) over blocks */ free(yin); free(xin); free(yn); free(xn); free(ydup); free(xdup); } /*----- end of test on equal or unequal size bins -----*/ free(good); RETURN(bim); } /*--------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------*/ static int ignore_zz = 0 ; void THD_correlate_ignore_zerozero( int z ){ ignore_zz = z ; } /*-------------------- Shannon entropy function --------------------*/ #define SHANENT(z) ( ( (z) <= 0.0f ) ? 0.0f : (z)*logf(z) ) /*--------------------------------------------------------------------------*/ /*! Compute the mutual info between two vectors, sort of. [16 Aug 2006] ----------------------------------------------------------------------------*/ float THD_mutual_info_scl( int n , float xbot,float xtop,float *x , float ybot,float ytop,float *y , float *w ) { register int ii,jj ; register float val ; /*-- build 2D histogram --*/ build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ; if( nbin <= 0 || nww <= 0 ) return 0.0f ; /* something bad happened! */ normalize_2Dhist() ; /*-- compute MI from histogram --*/ val = 0.0f ; for( ii=0 ; ii < nbp ; ii++ ){ for( jj=0 ; jj < nbp ; jj++ ){ #if 0 if( ii==0 && jj==0 && ignore_zz ) continue ; #endif val += SHANENT( XYC(ii,jj) ) ; }} return (1.4427f*val) ; /* units are bits, just for fun */ } /*--------------------------------------------------------------------------*/ /*! Compute MI of x[] and y[i], with autoscaling. */ float THD_mutual_info( int n , float *x , float *y ) { return THD_mutual_info_scl( n, 1.0f,-1.0f, x, 1.0f,-1.0f, y, NULL ) ; } /*--------------------------------------------------------------------------*/ /*! Compute the normalized mutual info between two vectors, sort of. Actually, returns H(x,y) / [ H(x)+H(y) ], which should be small if x and y are redundant and should be large if they are independent. ----------------------------------------------------------------------------*/ float THD_norm_mutinf_scl( int n , float xbot,float xtop,float *x , float ybot,float ytop,float *y , float *w ) { register int ii,jj ; register float numer , denom ; /*-- build 2D histogram --*/ build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ; if( nbin <= 0 || nww <= 0 ) return 0.0f ; /* something bad happened! */ normalize_2Dhist() ; /*-- compute NMI from histogram --*/ denom = numer = 0.0f ; for( ii=0 ; ii < nbp ; ii++ ){ denom += SHANENT( xc[ii] ) + SHANENT( yc[ii] ) ; for( jj=0 ; jj < nbp ; jj++ ){ #if 0 if( ii==0 && jj==0 && ignore_zz ) continue ; #endif numer += SHANENT( XYC(ii,jj) ) ; } } if( denom != 0.0f ) denom = numer / denom ; return denom ; } /*--------------------------------------------------------------------------*/ /*! Compute NMI of x[] & y[i], with autoscaling: cf. THD_norm_mutinf_scl(). */ float THD_norm_mutinf( int n , float *x , float *y ) { return THD_norm_mutinf_scl( n, 1.0f,-1.0f, x, 1.0f,-1.0f, y, NULL ) ; } /*--------------------------------------------------------------------------*/ /*! Compute the joint entropy between two vectors, sort of. ----------------------------------------------------------------------------*/ float THD_jointentrop_scl( int n , float xbot,float xtop,float *x , float ybot,float ytop,float *y , float *w ) { register int ii,jj ; register float val ; /*-- build 2D histogram --*/ build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ; if( nbin <= 0 || nww <= 0 ) return 0.0f ; /* something bad happened! */ normalize_2Dhist() ; /*-- compute MI from histogram --*/ val = 0.0f ; for( ii=0 ; ii < nbp ; ii++ ){ for( jj=0 ; jj < nbp ; jj++ ){ val -= SHANENT( XYC(ii,jj) ) ; }} return (1.4427f*val) ; /* units are bits, just for fun */ } /*--------------------------------------------------------------------------*/ /*! Compute joint entropy of x[] and y[i], with autoscaling. */ float THD_jointentrop( int n , float *x , float *y ) { return THD_jointentrop_scl( n, 1.0f,-1.0f, x, 1.0f,-1.0f, y, NULL ) ; } /*--------------------------------------------------------------------------*/ /* Decide if THD_corr_ratio_scl() computes symmetric or unsymmetric. */ static int cr_mode = 1 ; /* 0=unsym 1=sym mult 2=sym add */ void THD_corr_ratio_mode( int mm ){ cr_mode = mm ; } /*--------------------------------------------------------------------------*/ /*! Compute the correlation ratio between two vectors, sort of. [23 Aug 2006] ----------------------------------------------------------------------------*/ float THD_corr_ratio_scl( int n , float xbot,float xtop,float *x , float ybot,float ytop,float *y , float *w ) { register int ii,jj ; register float vv,mm ; float val , cyvar , uyvar , yrat,xrat ; /*-- build 2D histogram --*/ build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ; if( nbin <= 0 ) return 0.0f ; /* something bad happened! */ normalize_2Dhist() ; /*-- compute CR(y|x) from histogram --*/ cyvar = 0.0f ; for( ii=0 ; ii < nbp ; ii++ ){ if( xc[ii] > 0.0f ){ vv = mm = 0.0f ; /* mm=E(y|x) vv=E(y^2|x) */ for( jj=1 ; jj < nbp ; jj++ ){ mm += (jj * XYC(ii,jj)) ; vv += jj * (jj * XYC(ii,jj)) ; } cyvar += (vv - mm*mm/xc[ii] ) ; /* Var(y|x) */ } } vv = mm = uyvar = 0.0f ; for( jj=1 ; jj < nbp ; jj++ ){ /* mm=E(y) vv=E(y^2) */ mm += (jj * yc[jj]) ; vv += jj * (jj * yc[jj]) ; } uyvar = vv - mm*mm ; /* Var(y) */ yrat = (uyvar > 0.0f) ? cyvar/uyvar /* Var(y|x) / Var(y) */ : 1.0f ; if( cr_mode == 0 ) return (1.0f-yrat) ; /** unsymmetric **/ /** compute CR(x|y) also, for symmetrization **/ cyvar = 0.0f ; for( jj=0 ; jj < nbp ; jj++ ){ if( yc[jj] > 0.0f ){ vv = mm = 0.0f ; /* mm=E(x|y) vv=E(x^2|y) */ for( ii=1 ; ii < nbp ; ii++ ){ mm += (ii * XYC(ii,jj)) ; vv += ii * (ii * XYC(ii,jj)) ; } cyvar += (vv - mm*mm/yc[jj] ) ; /* Var(x|y) */ } } vv = mm = uyvar = 0.0f ; for( ii=1 ; ii < nbp ; ii++ ){ /* mm=E(x) vv=E(x^2) */ mm += (ii * xc[ii]) ; vv += ii * (ii * xc[ii]) ; } uyvar = vv - mm*mm ; /* Var(x) */ xrat = (uyvar > 0.0f) ? cyvar/uyvar /* Var(x|y) / Var(x) */ : 1.0f ; if( cr_mode == 2 ) return (1.0f - 0.5f*(xrat+yrat)) ; /** additive **/ return (1.0f - xrat*yrat) ; /** multiplicative **/ } /*--------------------------------------------------------------------------*/ /*! Compute CR of x[] and y[], using autoscaling. */ float THD_corr_ratio( int n , float *x , float *y ) { return THD_corr_ratio_scl( n , 1.0f,-1.0f , x, 1.0f,-1.0f , y , NULL ) ; } /*--------------------------------------------------------------------------*/ /*! Compute the Hellinger metric between two vectors, sort of. ----------------------------------------------------------------------------*/ float THD_hellinger_scl( int n , float xbot,float xtop,float *x , float ybot,float ytop,float *y , float *w ) { register int ii,jj ; register float val , pq ; /*-- build 2D histogram --*/ build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ; if( nbin <= 0 || nww <= 0 ) return 0.0f ; /* something bad happened! */ normalize_2Dhist() ; /*-- compute metric from histogram --*/ val = 0.0f ; for( ii=0 ; ii < nbp ; ii++ ){ for( jj=0 ; jj < nbp ; jj++ ){ pq = XYC(ii,jj) ; if( pq > 0.0f ) val += sqrtf( pq * xc[ii] * yc[jj] ) ; }} return (1.0f-val) ; } /*--------------------------------------------------------------------------*/ /*! Compute Hellinger metric between x[] and y[i], with autoscaling. */ float THD_hellinger( int n , float *x , float *y ) { return THD_hellinger_scl( n, 1.0f,-1.0f, x, 1.0f,-1.0f, y, NULL ) ; } /*--------------------------------------------------------------------------*/ /*! Compute the Hellinger metric, mutual info, normalized MI, and symmetrized correlation ratio, and return all 4 (in that order) using the 1D and 2D histograms from build_2Dhist(). The first 3 values all measure the closeness of the joint histogram to the product of the marginals: - Hellinger is smaller when the joint is closer to the marginals' product - MI is also smaller when the joint is closer to the marginal's product - NMI is larger when the joint is closer to the marginal's product Correlation ratio (symmetrized by addition == CRA) is larger when the two variables are nonlinearly correlated. As measures of association (generalized correlation): more closely associated variables correspond to larger Hellinger and MI and CRA, and to smaller NMI. *//*------------------------------------------------------------------------*/ float_quad THD_helmicra_scl( int n , float xbot,float xtop,float *x , float ybot,float ytop,float *y , float *w ) { register int ii,jj ; register float hel , pq , vv,uu ; float val , cyvar , uyvar , yrat,xrat ; float_quad hmc = {0.0f,0.0f,0.0f,0.f} ; /*-- build 2D histogram --*/ build_2Dhist( n,xbot,xtop,x,ybot,ytop,y,w ) ; if( nbin <= 0 || nww <= 0 ) return hmc ; /* something bad happened! */ normalize_2Dhist() ; /*-- compute Hel, MI, NMI from histogram --*/ hel = vv = uu = 0.0f ; for( ii=0 ; ii < nbp ; ii++ ){ vv += SHANENT( xc[ii] ) + SHANENT( yc[ii] ) ; for( jj=0 ; jj < nbp ; jj++ ){ pq = XYC(ii,jj) ; hel += sqrtf( pq * xc[ii] * yc[jj] ) ; uu += SHANENT( pq ) ; } } hmc.a = 1.0f - hel ; /* Hellinger */ hmc.b = uu - vv ; /* MI */ hmc.c = (vv != 0.0f) ? uu/vv : 0.0f ; /* NMI */ /*-- compute CR(y|x) from histogram --*/ cyvar = 0.0f ; for( ii=0 ; ii < nbp ; ii++ ){ if( xc[ii] > 0.0f ){ vv = uu = 0.0f ; /* uu=E(y|x) vv=E(y^2|x) */ for( jj=1 ; jj < nbp ; jj++ ){ uu += (jj * XYC(ii,jj)) ; vv += jj * (jj * XYC(ii,jj)) ; } cyvar += (vv - uu*uu/xc[ii] ) ; /* Var(y|x) */ } } vv = uu = uyvar = 0.0f ; for( jj=1 ; jj < nbp ; jj++ ){ /* uu=E(y) vv=E(y^2) */ uu += (jj * yc[jj]) ; vv += jj * (jj * yc[jj]) ; } uyvar = vv - uu*uu ; /* Var(y) */ yrat = (uyvar > 0.0f) ? cyvar/uyvar /* Var(y|x) / Var(y) */ : 1.0f ; /** compute CR(x|y) also, for symmetrization **/ cyvar = 0.0f ; for( jj=0 ; jj < nbp ; jj++ ){ if( yc[jj] > 0.0f ){ vv = uu = 0.0f ; /* uu=E(x|y) vv=E(x^2|y) */ for( ii=1 ; ii < nbp ; ii++ ){ uu += (ii * XYC(ii,jj)) ; vv += ii * (ii * XYC(ii,jj)) ; } cyvar += (vv - uu*uu/yc[jj] ) ; /* Var(x|y) */ } } vv = uu = uyvar = 0.0f ; for( ii=1 ; ii < nbp ; ii++ ){ /* uu=E(x) vv=E(x^2) */ uu += (ii * xc[ii]) ; vv += ii * (ii * xc[ii]) ; } uyvar = vv - uu*uu ; /* Var(x) */ xrat = (uyvar > 0.0f) ? cyvar/uyvar /* Var(x|y) / Var(x) */ : 1.0f ; hmc.d = 1.0f - 0.5f*(xrat+yrat) ; /** additive symmetrization **/ return hmc ; } /*--------------------------------------------------------------------------*/ /*! see THD_helmicra_scl(). */ float_quad THD_helmicra( int n , float *x , float *y ) { return THD_helmicra_scl( n, 1.0f,-1.0f, x, 1.0f,-1.0f, y, NULL ) ; } /*---------------------------------------------------------------------------*/ /* Return the MI and NMI of the 2 arrays y0[] and y1[] as 'correlated' with the binary 0-1 variable. That is, does knowing a value is in group 0 or 1 help predict it? *//*-------------------------------------------------------------------------*/ float_pair THD_binary_mutinfo( int n0 , float *y0 , int n1 , float *y1 ) { float_pair mmii = {-666.6f,666.6f} ; float *y0h , *y1h ; int nbin=64 , ii,kk ; float ymin,ymax , dy,dc, hsum,hjoint ; if( n0 < nbin || y0 == NULL || n1 < nbin || y1 == NULL ) return mmii ; /* data range */ ymin = ymax = y0[0] ; for( ii=1 ; ii < n0 ; ii++ ){ if( ymin > y0[ii] ) ymin = y0[ii] ; else if( ymax < y0[ii] ) ymax = y0[ii] ; } for( ii=0 ; ii < n1 ; ii++ ){ if( ymin > y1[ii] ) ymin = y1[ii] ; else if( ymax < y1[ii] ) ymax = y1[ii] ; } if( ymin >= ymax ) return mmii ; /* histograms */ y0h = (float *)calloc(sizeof(float),nbin) ; y1h = (float *)calloc(sizeof(float),nbin) ; dy = nbin / (ymax-ymin) ; dc = 1.0f / (float)(n0+n1) ; for( ii=0 ; ii < n0 ; ii++ ){ kk = ( y0[ii] - ymin ) * dy ; if( kk < 0 ) kk = 0 ; else if( kk >= nbin ) kk = nbin-1 ; y0h[kk] += dc ; } for( ii=0 ; ii < n1 ; ii++ ){ kk = (int)( ( y1[ii] - ymin ) * dy ) ; if( kk < 0 ) kk = 0 ; else if( kk >= nbin ) kk = nbin-1 ; y1h[kk] += dc ; } /* calculate mutual info vs. 0/1 */ /* entropy of the 0/1 'variable' */ hsum = - ( SHANENT(n0*dc) + SHANENT(n1*dc) ) ; hjoint = 0.0f ; for( kk=0 ; kk < nbin ; kk++ ){ hjoint -= ( SHANENT( y0h[kk] ) + SHANENT( y1h[kk] ) ) ; hsum -= SHANENT( y0h[kk] + y0h[kk] ) ; } free(y0h) ; free(y1h) ; mmii.a = hsum - hjoint ; if( hjoint > 0.0f ) mmii.b = hsum / hjoint ; return mmii ; } /*---------------------------------------------------------------------------*/ /*! Rank orders a collection of float arrays; ties get the average rank. The output overwrites the input. *//*-------------------------------------------------------------------------*/ void rank_order_float_arrays( int nar , int *nn , float **aa ) { int ii,jj,kk , jbase , ntot,nmax,n1,ns,ib ; int *b ; /* workspaces */ float *a , *c , cs ; /*- handle special cases -*/ if( nar < 1 || nn == NULL || aa == NULL ) return ; /* junk inputs */ if( nar == 1 ){ rank_order_float( nn[0] , aa[0] ) ; return ; /* just one input */ } ntot = nmax = 0 ; for( jj=0 ; jj < nar ; jj++ ){ ntot += nn[jj] ; if( nn[jj] > nmax ) nmax = nn[jj] ; } if( ntot < nar ) return ; /* bad inputs */ /*- make workspaces -*/ a = (float *)malloc(sizeof(float)*ntot) ; b = (int *)malloc(sizeof(int )*ntot) ; c = (float *)malloc(sizeof(float)*ntot) ; for( kk=jj=0 ; jj < nar ; jj++ ){ jbase = jj*nmax ; for( ii=0 ; ii < nn[jj] ; ii++,kk++ ){ a[kk] = aa[jj][ii] ; /* data */ b[kk] = ii + jbase ; /* where it came from */ c[kk] = (float)kk ; /* default rank of kk-th value after sorting */ } } /*- sort input, carrying b along -*/ qsort_floatint( ntot , a , b ) ; /* c[kk] now contains the global rank of the kk-th element, but c[] perhaps needs to modified to allow for ties ==> average rank */ n1 = ntot-1 ; for( ii=0 ; ii < n1 ; ii++ ){ if( a[ii] == a[ii+1] ){ /* handle ties */ cs = (ii)+(ii+1) ; ns = 2 ; ib = ii ; ii++ ; while( ii < n1 && a[ii] == a[ii+1] ){ ii++ ; ns++ ; cs += ii ; } for( cs/=ns ; ib <= ii ; ib++ ) c[ib] = cs ; } } /* put c[] results back into the original arrays */ for( kk=0 ; kk < ntot ; kk++ ){ jj = b[kk] / nmax ; /* which array did it come from? */ ii = b[kk] % nmax ; /* which element in that array? */ aa[jj][ii] = c[kk] ; /* replace data with rank */ } free(c) ; free(b) ; free(a) ; return ; } /*---------------------------------------------------------------------------*/ void rank_order_2floats( int n1 , float *a1 , int n2 , float *a2 ) { int nn[2] ; float *aa[2] ; if( n1 <= 0 || n2 <= 0 || a1 == NULL || a2 == NULL ) return ; nn[0] = n1 ; nn[1] = n2 ; aa[0] = a1 ; aa[1] = a2 ; rank_order_float_arrays( 2 , nn , aa ) ; return ; } /*----------------------------------------------------------------------------*/ /* Given a set of bootstrap replicates, find the bias-corrected (BC) estimates of the 1-alpha confidence interval and the central value. * estim = actual estimate from the data * alpha = 0.05 is a typical value * nboot = number of replicates in eboot[] -- at least MAX(100,10/alpha) * eboot = array of replicates; will be sorted on output The return value (rval) is a triple: * rval.a = lower edge of confidence interval * rval.b = middle value = bias-corrected estimate * rval.c = upper edge of confidence interval If all values are returned as 0, the inputs were bad bad bad. *//*--------------------------------------------------------------------------*/ #undef PHI #define PHI(x) (1.0-qg(x)) /* CDF of N(0,1) */ #undef PHINV #define PHINV(x) qginv(1.0-(x)) /* inverse CDF of N(0,1) */ #undef ZLIM #define ZLIM 0.5f float_triple THD_bootstrap_confinv( float estim , float alpha , int nboot , float *eboot ) { float_triple rval = {0.0f,0.0f,0.0f} ; int ii ; float z0 , zal , pp ; ENTRY("THD_bootstrap_confinv") ; if( nboot < 100 || eboot == NULL ) RETURN( rval ) ; /* bad user */ if( alpha <= 0.001f || alpha >= 0.9f ) alpha = 0.05f ; /* stupid user */ alpha *= 0.5f ; /* 1-sided tail */ if( (int)(alpha*nboot) < 5 ) alpha = 5.0f / nboot ; /* upward adjust */ zal = PHINV(alpha) ; /* should be negative */ qsort_float( nboot , eboot ) ; /* increasing order */ for( ii=0 ; ii < nboot && eboot[ii] < estim ; ii++ ) ; /*nada*/ if( ii <= 1 || ii >= nboot-1 ) RETURN( rval ) ; /* crummy data */ z0 = PHINV( (ii+0.5f) / nboot ) ; /* ii = #values < estim */ if( z0 < -ZLIM ) z0 = -ZLIM ; else if( z0 > ZLIM ) z0 = ZLIM ; /* limits */ pp = PHI( 2.0*z0 + zal ) * nboot ; /* lower edge */ ii = (int)pp ; pp = pp - ii ; if( ii >= nboot-1 ) ii = nboot-2 ; rval.a = (1.0f-pp)*eboot[ii] + pp*eboot[ii+1] ; pp = PHI( 2.0*z0 - zal ) * nboot ; /* upper edge */ ii = (int)pp ; pp = pp - ii ; if( ii >= nboot-1 ) ii = nboot-2 ; rval.c = (1.0f-pp)*eboot[ii] + pp*eboot[ii+1] ; pp = PHI( 2.0*z0 ) * nboot ; /* bias-corrected central estimate */ ii = (int)pp ; pp = pp - ii ; if( ii >= nboot-1 ) ii = nboot-2 ; rval.b = (1.0f-pp)*eboot[ii] + pp*eboot[ii+1] ; RETURN( rval ) ; } /*----------------------------------------------------------------------------*/ float THD_bootstrap_biascorr( float estim , int nboot , float *eboot ) { int ii ; float z0 , pp ; if( nboot < 50 || eboot == NULL ) return estim ; /* bad user */ qsort_float( nboot , eboot ) ; /* increasing order */ for( ii=0 ; ii < nboot && eboot[ii] < estim ; ii++ ) ; /*nada*/ if( ii <= 1 || ii >= nboot-1 ) return estim ; /* crummy data */ z0 = PHINV( (ii+0.5f) / nboot ) ; /* ii = #values < estim */ if( z0 < -ZLIM ) z0 = -ZLIM ; else if( z0 > ZLIM ) z0 = ZLIM ; /* limits */ pp = PHI( 2.0*z0 ) * nboot ; ii = (int)pp ; pp = pp - ii ; if( ii >= nboot-1 ) ii = nboot-2 ; pp = (1.0f-pp)*eboot[ii] + pp*eboot[ii+1] ; return pp ; } /*----------------------------------------------------------------------------*/ #undef DEMEAN #define DEMEAN(n,v) do{ register int i ; register float s ; \ for( s=i=0 ; i < n ; i++ ) s += v[i] ; \ s /= n ; for( i=0 ; i < n ; i++ ) v[i] -= s ; \ } while(0) #undef XPT #define XPT(q) ( (xtyp<=0) ? xx+(q)*nlen : xpt[q] ) #undef YPT #define YPT(q) ( (xtyp<=0) ? yy+(q)*nlen : ypt[q] ) /*----------------------------------------------------------------------------*/ /* Return bootstrap BC estimate for correlation of a bunch of vectors. *//*--------------------------------------------------------------------------*/ float THD_bootstrap_vectcorr( int nlen, int nboot, int use_pv, int xtyp, int xnum, void *xp , int ynum , void *yp ) { float rval, rxy, **xar, **yar, *rboot, *xbar, *ybar, *uvec, *vvec ; int kb,ii,jj , dox,doy ; unsigned short xran[3] ; float *xx=NULL, **xpt=NULL, *yy=NULL, **ypt=NULL ; void *pvw=NULL ; ENTRY("THD_bootstrap_vectcorr") ; if( nlen < 3 || xnum < 1 || ynum < 1 || xp == NULL || yp == NULL ) RETURN(0.0f) ; if( xtyp <= 0 ){ xx = (float *) xp ; yy = (float *) yp ; } else { xpt = (float **)xp ; ypt = (float **)yp ; } /* trivial case */ if( xnum == 1 && ynum == 1 ){ rval = THD_pearson_corr(nlen,XPT(0),YPT(0)) ; RETURN(rval) ; } /* compute mean vectors */ #pragma omp critical (MALLOC) { xbar = (float *)malloc(sizeof(float)*nlen) ; ybar = (float *)malloc(sizeof(float)*nlen) ; } (void)mean_vector( nlen , xnum , xtyp , xp , xbar ) ; (void)mean_vector( nlen , ynum , xtyp , yp , ybar ) ; dox = (xnum > 5) ; doy = (ynum > 5) ; /* which to bootstrap */ if( !dox && !doy ){ rval = THD_pearson_corr(nlen,xbar,ybar) ; #pragma omp critical (MALLOC) { free(ybar) ; free(xbar) ; } RETURN(rval) ; } if( nboot < 50 ) nboot = 50 ; #pragma omp critical (MALLOC) { xar = (float **)malloc(sizeof(float *)*xnum ) ; yar = (float **)malloc(sizeof(float *)*ynum ) ; rboot = (float * )malloc(sizeof(float) *nboot) ; uvec = (float * )malloc(sizeof(float) *nlen ) ; vvec = (float * )malloc(sizeof(float) *nlen ) ; } /* correlation with no resampling at all */ if( use_pv ){ pvw = pv_get_workspace( nlen , MAX(xnum,ynum) ) ; xran[0] = (unsigned short)(xnum+ynum+nlen) ; xran[1] = 42731 ; xran[2] = 23172 ; (void)principal_vector( nlen, xnum, xtyp, xp, uvec, xbar, pvw, xran) ; (void)principal_vector( nlen, ynum, xtyp, yp, vvec, ybar, pvw, xran) ; } else { for( ii=0 ; ii < nlen ; ii++ ){ uvec[ii] = xbar[ii]; vvec[ii] = ybar[ii]; } } rxy = THD_pearson_corr( nlen , uvec , vvec ) ; /* bootstrap correlations [selecting subsets from each set of vectors] */ for( kb=0 ; kb < nboot ; kb++ ){ if( dox ){ for( ii=0 ; ii < xnum ; ii++ ){ /* resample xx vectors */ jj = nrand48(xran) % xnum ; xar[ii] = XPT(jj) ; } if( use_pv ) (void)principal_vector( nlen, xnum, 1, xar, uvec, xbar, pvw, xran) ; else mean_vector( nlen , xnum , 1 , xar , uvec ) ; } if( doy ){ for( ii=0 ; ii < ynum ; ii++ ){ /* resample yy vectors */ jj = nrand48(xran) % ynum ; yar[ii] = YPT(jj) ; } if( use_pv ) (void)principal_vector( nlen, ynum, 1, yar, vvec, ybar, pvw, xran) ; else mean_vector( nlen , ynum , 1 , yar , vvec ) ; } rboot[kb] = THD_pearson_corr( nlen , uvec , vvec ) ; } /* final answer */ rval = THD_bootstrap_biascorr( rxy , nboot , rboot ) ; /* cleanup the trash */ #pragma omp critical (MALLOC) { free(vvec); free(uvec); free(rboot); free(yar); free(xar); free(ybar); free(xbar); if(use_pv)free(pvw); } RETURN(rval) ; }