00001 #include "mrilib.h"
00002 
00003 /* prototypes */
00004 
00005 void qsort_double( int n , double * a ) ;
00006 void gaussians( int n , double * g ) ;
00007 void normalize_doublevector( int n , double * far ) ;
00008 double rtest( int msam , int nvec , double * lamq , double * rho ) ;
00009 
00010 /*----------------------------------------------------------------------------
00011   Generate a bunch of N(0,1) deviates
00012 ------------------------------------------------------------------------------*/
00013 
00014 void gaussians( int n , double * g )
00015 {
00016    double u1,u2 ;
00017    int ii ;
00018 
00019    for( ii=0 ; ii < n ; ii+=2 ){
00020       do{ u1 = drand48() ; } while( u1==0.0 ) ;
00021       u1 = sqrt(-2.0*log(u1)) ;
00022       u2 = 6.2831853 * drand48() ;
00023       g[ii] = u1 * sin(u2) ;
00024       if( ii < n-1 ) g[ii+1] = u1 * cos(u2) ;
00025    }
00026    return ;
00027 }
00028 
00029 /*-----------------------------------------------------------------------------
00030    Remove the mean, and make it have L2 norm = 1
00031 -------------------------------------------------------------------------------*/
00032 
00033 void normalize_doublevector( int n , double * far )
00034 {
00035    int ii ;
00036    double ff,gg ;
00037 
00038    for( ff=0.0,ii=0 ; ii < n ; ii++ ) ff += far[ii] ;
00039    ff /= n ;
00040    for( gg=0.0,ii=0 ; ii < n ; ii++ ) gg += SQR( (far[ii]-ff) ) ;
00041    if( gg <= 0.0 ) return ;
00042    gg = 1.0 / sqrt(gg) ;
00043    for( ii=0 ; ii < n ; ii++ ) far[ii] = (far[ii]-ff)*gg ;
00044 
00045    return ;
00046 }
00047 
00048 /*-----------------------------------------------------------------------------
00049    El Supremo function!!!
00050 -------------------------------------------------------------------------------*/
00051 
00052 #define BBOT 10
00053 
00054 int main( int argc , char * argv[] )
00055 {
00056    int nopt , mdev , jsize , ii,jj , numj,jset[100] ;
00057    double * gset , * pset , * rr ;
00058    double pp , dm ;
00059 
00060    srand48((long)time(NULL)) ;  /* initialize randomness */
00061 
00062    /* read command line */
00063 
00064    if( argc < 3 || strcmp(argv[1],"-help") == 0 ){
00065       printf("Usage: rtest1 M J [j1 j2 j3 ...]\n"
00066              "   M = # of deviates per sample sets\n"
00067              "   J = # of sample sets to generate\n"
00068              "  j1 = first statistic to output, etc.\n"
00069             ) ;
00070       exit(0) ;
00071    }
00072 
00073    nopt = 1 ;
00074    mdev = (int) strtod(argv[nopt++],NULL) ;
00075    if( mdev < BBOT ){ fprintf(stderr,"M=%d is too small!\n",mdev); exit(1); }
00076    jsize = (int) strtod(argv[nopt++],NULL) ;
00077    if( jsize < BBOT ){ fprintf(stderr,"J=%d is too small!\n",jsize); exit(1); }
00078 
00079    for( numj=0 ; numj < 100 && nopt < argc ; ){
00080       ii = (int) strtod(argv[nopt++],NULL) ;
00081       if( ii < 0 || ii >= jsize )
00082          fprintf(stderr,"** j%d = %d is illegal!\n",nopt-2,ii) ;
00083       else
00084          jset[numj++] = ii ;
00085    }
00086 
00087    /* make space for sample sets and results */
00088 
00089    gset = (double *) malloc(sizeof(double)*mdev) ;  /* 1 sample set */
00090    pset = (double *) malloc(sizeof(double)*mdev) ;  /* cdf-inverses */
00091    rr   = (double *) malloc(sizeof(double)*jsize) ; /* results      */
00092 
00093    /* initialize pset to quantile points on the N(0,1) cdf */
00094 
00095    dm = 1.0 / (double) mdev ;
00096    for( ii=0 ; ii < mdev ; ii++ ){
00097       pp = dm * (ii+0.5) ; pset[ii] = qginv( 1.0 - pp ) ;
00098    }
00099    normalize_doublevector( mdev , pset ) ;  /* prepare for correlation */
00100 
00101    /* loop over realization of sample sets */
00102 
00103    for( jj=0 ; jj < jsize ; jj++ ){
00104       gaussians( mdev , gset ) ;              /* make a sample set       */
00105       qsort_double( mdev , gset ) ;           /* sort it                 */
00106       normalize_doublevector( mdev , gset ) ; /* prepare for correlation */
00107       for( pp=0.0,ii=0 ; ii < mdev ; ii++ )   /* compute 1.0-correlation */
00108          pp += dm - gset[ii]*pset[ii] ;
00109       rr[jj] = -pp ;                          /* save it (or minus it)   */
00110 
00111       if( jj > 0 && jj%10000 == 0 ) fprintf(stderr,"at #%d rr=%g\n",jj,pp) ;
00112    }
00113 
00114    free(gset) ; free(pset) ;                  /* toss some trash */
00115 
00116    /* sort results */
00117 
00118    qsort_double( jsize , rr ) ;
00119    for( jj=0 ; jj < jsize ; jj++ ) rr[jj] = -rr[jj] ;  /* back to 1-corr */
00120 
00121    /* print results */
00122 
00123    if( numj > 0 ){
00124       dm = 1.0 / jsize ;
00125       for( jj=0 ; jj < numj ; jj++ )
00126          printf(" %.4f %g\n", dm*(jset[jj]+1.0) , rr[jset[jj]] ) ;
00127    } else {
00128       for( jj=0 ; jj < jsize ; jj++ ) printf(" %g",rr[jj]) ;
00129       printf("\n");
00130    }
00131 
00132    exit(0) ;
00133 }
00134 
00135 double rtest( int msam , int nvec , double * lamq , double * rho )
00136 {
00137    int ii , jj ;
00138    double * zz , * gg ;
00139 
00140    zz = (double *) malloc(sizeof(double)*msam) ;
00141    gg = (double *) malloc(sizeof(double)*nvec) ;
00142 }
00143 
00144 /*------------------------------------------------------------------------------*/
00145 /*------------- insertion sort : sort an array of double in-place --------------*/
00146 
00147 void isort_double( int n , double * ar )
00148 {
00149    register int  j , p ;  /* array indices */
00150    register double temp ;  /* a[j] holding place */
00151    register double * a = ar ;
00152 
00153    if( n < 2 || ar == NULL ) return ;
00154 
00155    for( j=1 ; j < n ; j++ ){
00156 
00157      if( a[j] < a[j-1] ){  /* out of order */
00158         p    = j ;
00159         temp = a[j] ;
00160         do{
00161           a[p] = a[p-1] ;       /* at this point, a[p-1] > temp, so move it up */
00162           p-- ;
00163         } while( p > 0 && temp < a[p-1] ) ;
00164         a[p] = temp ;           /* finally, put temp in its place */
00165      }
00166    }
00167    return ;
00168 }
00169 
00170 #define QS_CUTOFF     20       /* cutoff to switch from qsort to isort */
00171 #define QS_STACK      1024     /* qsort stack size (way too big) */
00172 #define QS_SWAP(x,y)  (temp=(x), (x)=(y),(y)=temp )
00173 #define QS_SWAPI(x,y) (itemp=(x),(x)=(y),(y)=itemp)
00174 
00175 static int stack[QS_STACK] ;   /* stack for qsort "recursion" */
00176 
00177 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
00178 
00179 void qsrec_double( int n , double * ar , int cutoff )
00180 {
00181    register int i , j ;            /* scanning indices */
00182    register double temp , pivot ;  /* holding places */
00183    register double * a = ar ;
00184 
00185    int left , right , mst , itemp,nnew ;
00186 
00187    /* return if too short (insertion sort will clean up) */
00188 
00189    if( cutoff < 3 ) cutoff = 3 ;
00190    if( n < cutoff || ar == NULL ) return ;
00191 
00192    /* initialize stack to start with whole array */
00193 
00194    stack[0] = 0   ;
00195    stack[1] = n-1 ;
00196    mst      = 2   ;
00197 
00198    /* loop while the stack is nonempty */
00199 
00200    while( mst > 0 ){
00201       right = stack[--mst] ;  /* work on subarray from left -> right */
00202       left  = stack[--mst] ;
00203 
00204       i = ( left + right ) / 2 ;           /* middle of subarray */
00205 
00206       /*----- sort the left, middle, and right a[]'s -----*/
00207 
00208       if( a[left] > a[i]     ) QS_SWAP( a[left]  , a[i]     ) ;
00209       if( a[left] > a[right] ) QS_SWAP( a[left]  , a[right] ) ;
00210       if( a[i] > a[right]    ) QS_SWAP( a[right] , a[i]     ) ;
00211 
00212       pivot = a[i] ;                       /* a[i] is the median-of-3 pivot! */
00213       a[i]  = a[right] ;
00214 
00215       i = left ; j = right ;               /* initialize scanning */
00216 
00217       /*----- partition:  move elements bigger than pivot up and elements
00218                           smaller than pivot down, scanning in from ends -----*/
00219 
00220       do{
00221         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00222         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00223 
00224         if( j <= i ) break ;         /* if j meets i, quit */
00225 
00226         QS_SWAP( a[i] , a[j] ) ;
00227       } while( 1 ) ;
00228 
00229       /*----- at this point, the array is partitioned -----*/
00230 
00231       a[right] = a[i] ; a[i] = pivot ;  /* restore the pivot */
00232 
00233       /*----- signal the subarrays that need further work -----*/
00234 
00235       nnew = 0 ;
00236       if( (i-left)  > cutoff ){ stack[mst++] = left; stack[mst++] = i-1  ; nnew++; }
00237       if( (right-i) > cutoff ){ stack[mst++] = i+1 ; stack[mst++] = right; nnew++; }
00238 
00239       /* if just added two subarrays to stack, make sure shorter one comes first */
00240 
00241       if( nnew == 2 && stack[mst-3] - stack[mst-4] > stack[mst-1] - stack[mst-2] ){
00242          QS_SWAPI( stack[mst-4] , stack[mst-2] ) ;
00243          QS_SWAPI( stack[mst-3] , stack[mst-1] ) ;
00244       }
00245 
00246    }  /* end of while stack is non-empty */
00247    return ;
00248 }
00249 
00250 /* sort an array partially recursively, and partially insertion */
00251 
00252 void qsort_double( int n , double * a )
00253 {
00254    qsrec_double( n , a , QS_CUTOFF ) ;
00255    isort_double( n , a ) ;
00256    return ;
00257 }