1 #include "mrilib.h"
2
3 /* prototypes */
4
5 void qsort_double( int n , double * a ) ;
6 void gaussians( int n , double * g ) ;
7 void normalize_doublevector( int n , double * far ) ;
8 double rtest( int msam , int nvec , double * lamq , double * rho ) ;
9
10 /*----------------------------------------------------------------------------
11 Generate a bunch of N(0,1) deviates
12 ------------------------------------------------------------------------------*/
13
gaussians(int n,double * g)14 void gaussians( int n , double * g )
15 {
16 double u1,u2 ;
17 int ii ;
18
19 for( ii=0 ; ii < n ; ii+=2 ){
20 do{ u1 = drand48() ; } while( u1==0.0 ) ;
21 u1 = sqrt(-2.0*log(u1)) ;
22 u2 = 6.2831853 * drand48() ;
23 g[ii] = u1 * sin(u2) ;
24 if( ii < n-1 ) g[ii+1] = u1 * cos(u2) ;
25 }
26 return ;
27 }
28
29 /*-----------------------------------------------------------------------------
30 Remove the mean, and make it have L2 norm = 1
31 -------------------------------------------------------------------------------*/
32
normalize_doublevector(int n,double * far)33 void normalize_doublevector( int n , double * far )
34 {
35 int ii ;
36 double ff,gg ;
37
38 for( ff=0.0,ii=0 ; ii < n ; ii++ ) ff += far[ii] ;
39 ff /= n ;
40 for( gg=0.0,ii=0 ; ii < n ; ii++ ) gg += SQR( (far[ii]-ff) ) ;
41 if( gg <= 0.0 ) return ;
42 gg = 1.0 / sqrt(gg) ;
43 for( ii=0 ; ii < n ; ii++ ) far[ii] = (far[ii]-ff)*gg ;
44
45 return ;
46 }
47
48 /*-----------------------------------------------------------------------------
49 El Supremo function!!!
50 -------------------------------------------------------------------------------*/
51
52 #define BBOT 10
53
main(int argc,char * argv[])54 int main( int argc , char * argv[] )
55 {
56 int nopt , mdev , jsize , ii,jj , numj,jset[100] ;
57 double * gset , * pset , * rr ;
58 double pp , dm ;
59
60 srand48((long)time(NULL)) ; /* initialize randomness */
61
62 /* read command line */
63
64 if( argc < 3 || strcmp(argv[1],"-help") == 0 ){
65 printf("Usage: rtest1 M J [j1 j2 j3 ...]\n"
66 " M = # of deviates per sample sets\n"
67 " J = # of sample sets to generate\n"
68 " j1 = first statistic to output, etc.\n"
69 ) ;
70 exit(0) ;
71 }
72
73 nopt = 1 ;
74 mdev = (int) strtod(argv[nopt++],NULL) ;
75 if( mdev < BBOT ){ fprintf(stderr,"M=%d is too small!\n",mdev); exit(1); }
76 jsize = (int) strtod(argv[nopt++],NULL) ;
77 if( jsize < BBOT ){ fprintf(stderr,"J=%d is too small!\n",jsize); exit(1); }
78
79 for( numj=0 ; numj < 100 && nopt < argc ; ){
80 ii = (int) strtod(argv[nopt++],NULL) ;
81 if( ii < 0 || ii >= jsize )
82 fprintf(stderr,"** j%d = %d is illegal!\n",nopt-2,ii) ;
83 else
84 jset[numj++] = ii ;
85 }
86
87 /* make space for sample sets and results */
88
89 gset = (double *) malloc(sizeof(double)*mdev) ; /* 1 sample set */
90 pset = (double *) malloc(sizeof(double)*mdev) ; /* cdf-inverses */
91 rr = (double *) malloc(sizeof(double)*jsize) ; /* results */
92
93 /* initialize pset to quantile points on the N(0,1) cdf */
94
95 dm = 1.0 / (double) mdev ;
96 for( ii=0 ; ii < mdev ; ii++ ){
97 pp = dm * (ii+0.5) ; pset[ii] = qginv( 1.0 - pp ) ;
98 }
99 normalize_doublevector( mdev , pset ) ; /* prepare for correlation */
100
101 /* loop over realization of sample sets */
102
103 for( jj=0 ; jj < jsize ; jj++ ){
104 gaussians( mdev , gset ) ; /* make a sample set */
105 qsort_double( mdev , gset ) ; /* sort it */
106 normalize_doublevector( mdev , gset ) ; /* prepare for correlation */
107 for( pp=0.0,ii=0 ; ii < mdev ; ii++ ) /* compute 1.0-correlation */
108 pp += dm - gset[ii]*pset[ii] ;
109 rr[jj] = -pp ; /* save it (or minus it) */
110
111 if( jj > 0 && jj%10000 == 0 ) fprintf(stderr,"at #%d rr=%g\n",jj,pp) ;
112 }
113
114 free(gset) ; free(pset) ; /* toss some trash */
115
116 /* sort results */
117
118 qsort_double( jsize , rr ) ;
119 for( jj=0 ; jj < jsize ; jj++ ) rr[jj] = -rr[jj] ; /* back to 1-corr */
120
121 /* print results */
122
123 if( numj > 0 ){
124 dm = 1.0 / jsize ;
125 for( jj=0 ; jj < numj ; jj++ )
126 printf(" %.4f %g\n", dm*(jset[jj]+1.0) , rr[jset[jj]] ) ;
127 } else {
128 for( jj=0 ; jj < jsize ; jj++ ) printf(" %g",rr[jj]) ;
129 printf("\n");
130 }
131
132 exit(0) ;
133 }
134
rtest(int msam,int nvec,double * lamq,double * rho)135 double rtest( int msam , int nvec , double * lamq , double * rho )
136 {
137 int ii , jj ;
138 double * zz , * gg ;
139
140 zz = (double *) malloc(sizeof(double)*msam) ;
141 gg = (double *) malloc(sizeof(double)*nvec) ;
142 }
143
144 /*------------------------------------------------------------------------------*/
145 /*------------- insertion sort : sort an array of double in-place --------------*/
146
isort_double(int n,double * ar)147 void isort_double( int n , double * ar )
148 {
149 register int j , p ; /* array indices */
150 register double temp ; /* a[j] holding place */
151 register double * a = ar ;
152
153 if( n < 2 || ar == NULL ) return ;
154
155 for( j=1 ; j < n ; j++ ){
156
157 if( a[j] < a[j-1] ){ /* out of order */
158 p = j ;
159 temp = a[j] ;
160 do{
161 a[p] = a[p-1] ; /* at this point, a[p-1] > temp, so move it up */
162 p-- ;
163 } while( p > 0 && temp < a[p-1] ) ;
164 a[p] = temp ; /* finally, put temp in its place */
165 }
166 }
167 return ;
168 }
169
170 #undef QS_CUTOFF
171 #undef QS_SWAP
172 #undef QS_SWAPI
173 #define QS_CUTOFF 20 /* cutoff to switch from qsort to isort */
174 #define QS_SWAP(x,y) (temp=(x), (x)=(y),(y)=temp )
175 #define QS_SWAPI(x,y) (itemp=(x),(x)=(y),(y)=itemp)
176 #ifndef QS_STACK
177 # define QS_STACK 9999
178 #endif
179
180 static int stack[QS_STACK] ; /* stack for qsort "recursion" */
181
182 /*--------- qsrec : recursive part of quicksort (stack implementation) ----------*/
183
qsrec_double(int n,double * ar,int cutoff)184 void qsrec_double( int n , double * ar , int cutoff )
185 {
186 register int i , j ; /* scanning indices */
187 register double temp , pivot ; /* holding places */
188 register double * a = ar ;
189
190 int left , right , mst , itemp,nnew ;
191
192 /* return if too short (insertion sort will clean up) */
193
194 if( cutoff < 3 ) cutoff = 3 ;
195 if( n < cutoff || ar == NULL ) return ;
196
197 /* initialize stack to start with whole array */
198
199 stack[0] = 0 ;
200 stack[1] = n-1 ;
201 mst = 2 ;
202
203 /* loop while the stack is nonempty */
204
205 while( mst > 0 ){
206 right = stack[--mst] ; /* work on subarray from left -> right */
207 left = stack[--mst] ;
208
209 i = ( left + right ) / 2 ; /* middle of subarray */
210
211 /*----- sort the left, middle, and right a[]'s -----*/
212
213 if( a[left] > a[i] ) QS_SWAP( a[left] , a[i] ) ;
214 if( a[left] > a[right] ) QS_SWAP( a[left] , a[right] ) ;
215 if( a[i] > a[right] ) QS_SWAP( a[right] , a[i] ) ;
216
217 pivot = a[i] ; /* a[i] is the median-of-3 pivot! */
218 a[i] = a[right] ;
219
220 i = left ; j = right ; /* initialize scanning */
221
222 /*----- partition: move elements bigger than pivot up and elements
223 smaller than pivot down, scanning in from ends -----*/
224
225 do{
226 for( ; a[++i] < pivot ; ) ; /* scan i up, until a[i] >= pivot */
227 for( ; a[--j] > pivot ; ) ; /* scan j down, until a[j] <= pivot */
228
229 if( j <= i ) break ; /* if j meets i, quit */
230
231 QS_SWAP( a[i] , a[j] ) ;
232 } while( 1 ) ;
233
234 /*----- at this point, the array is partitioned -----*/
235
236 a[right] = a[i] ; a[i] = pivot ; /* restore the pivot */
237
238 /*----- signal the subarrays that need further work -----*/
239
240 nnew = 0 ;
241 if( (i-left) > cutoff ){ stack[mst++] = left; stack[mst++] = i-1 ; nnew++; }
242 if( (right-i) > cutoff ){ stack[mst++] = i+1 ; stack[mst++] = right; nnew++; }
243
244 /* if just added two subarrays to stack, make sure shorter one comes first */
245
246 if( nnew == 2 && stack[mst-3] - stack[mst-4] > stack[mst-1] - stack[mst-2] ){
247 QS_SWAPI( stack[mst-4] , stack[mst-2] ) ;
248 QS_SWAPI( stack[mst-3] , stack[mst-1] ) ;
249 }
250
251 } /* end of while stack is non-empty */
252 return ;
253 }
254
255 /* sort an array partially recursively, and partially insertion */
256
qsort_double(int n,double * a)257 void qsort_double( int n , double * a )
258 {
259 qsrec_double( n , a , QS_CUTOFF ) ;
260 isort_double( n , a ) ;
261 return ;
262 }
263