1 #include "mrilib.h"
2
usage_1ddot(int detail)3 void usage_1ddot(int detail) {
4 printf("Usage: 1ddot [options] 1Dfile 1Dfile ...\n"
5 "* Prints out correlation matrix of the 1D files and\n"
6 " their inverse correlation matrix.\n"
7 "* Output appears on stdout.\n"
8 "* Program 1dCorrelate does something similar-ish.\n"
9 "\n"
10 "Options:\n"
11 " -one = Make 1st vector be all 1's.\n"
12 " -dem = Remove mean from all vectors (conflicts with '-one')\n"
13 " -cov = Compute with covariance matrix instead of correlation\n"
14 " -inn = Computed with inner product matrix instead\n"
15 " -rank = Compute Spearman rank correlation instead\n"
16 " (also implies '-terse')\n"
17 " -terse= Output only the correlation or covariance matrix\n"
18 " and without any of the garnish. \n"
19 " -okzero= Do not quit if a vector is all zeros.\n"
20 " The correlation matrix will have 0 where NaNs ought to go.\n"
21 " Expect rubbish in the inverse matrices if all zero \n"
22 " vectors exist.\n"
23 ) ;
24 PRINT_COMPILE_DATE ;
25 return;
26 }
27
28 /*----------------------------------------------------------------------------*/
29
main(int argc,char * argv[])30 int main( int argc , char *argv[] )
31 {
32 int iarg , ii,jj,kk,mm , nvec , do_one=0 , nx=0,ny , ff, doterse = 0 ;
33 MRI_IMAGE *tim ;
34 MRI_IMARR *tar ;
35 double sum , *eval , *amat , **tvec , *bmat , *svec ;
36 float *far ;
37 int demean=0 , docov=0 ;
38 char *matname ;
39 int okzero = 0;
40
41 mainENTRY("1ddot main"); machdep();
42 /* options */
43
44 iarg = 1 ; nvec = 0 ;
45 while( iarg < argc && argv[iarg][0] == '-' ){
46
47 if (strcmp(argv[iarg],"-help") == 0 || strcmp(argv[iarg],"-h") == 0){
48 usage_1ddot(strlen(argv[iarg])>3 ? 2:1);
49 exit(0);
50 }
51
52 if( strcmp(argv[iarg],"-one") == 0 ){
53 demean = 0 ; do_one = 1 ; iarg++ ; continue ;
54 }
55
56 if( strcmp(argv[iarg],"-okzero") == 0 ){
57 okzero = 1 ; iarg++ ; continue ;
58 }
59
60 if( strncmp(argv[iarg],"-dem",4) == 0 ){
61 demean = 1 ; do_one = 0 ; iarg++ ; continue ;
62 }
63
64 if( strncmp(argv[iarg],"-cov",4) == 0 ){
65 docov = 1 ; iarg++ ; continue ;
66 }
67
68 if( strncmp(argv[iarg],"-inn",4) == 0 ){
69 docov = 2 ; iarg++ ; continue ;
70 }
71
72 if( strcasecmp(argv[iarg],"-rank") == 0 ||
73 strcasecmp(argv[iarg],"-spearman") == 0 ){
74 do_one = 0; docov = 3; demean = 0; doterse = 1; iarg++; continue;
75 }
76
77 if( strncmp(argv[iarg],"-terse",4) == 0 ){
78 doterse = 1 ; iarg++ ; continue ;
79 }
80
81 fprintf(stderr,"** Unknown option: %s\n",argv[iarg]);
82 suggest_best_prog_option(argv[0], argv[iarg]);
83 exit(1);
84 }
85
86 if( argc < 2 ){ usage_1ddot(1); exit(0) ; }
87
88 if( iarg == argc ) ERROR_exit("No 1D files on command line!?") ;
89
90 /* input 1D files */
91
92 ff = iarg ;
93 INIT_IMARR(tar) ; if( do_one ) nvec = 1 ;
94 for( ; iarg < argc ; iarg++ ){
95 tim = mri_read_1D( argv[iarg] ) ;
96 if( tim == NULL ){
97 fprintf(stderr,"** Can't read 1D file %s\n",argv[iarg]); exit(1);
98 }
99 if( nx == 0 ){
100 nx = tim->nx ;
101 } else if( tim->nx != nx ){
102 fprintf(stderr,"** 1D file %s doesn't match first file in length!\n",
103 argv[iarg]); exit(1);
104 }
105 nvec += tim->ny ;
106 ADDTO_IMARR(tar,tim) ;
107 }
108
109 if (!doterse) {
110 printf("\n") ;
111 printf("++ 1ddot input vectors:\n") ;
112 }
113 jj = 0 ;
114 if( do_one ){
115 if (!doterse) printf("00..00: all ones\n") ;
116 jj = 1 ;
117 }
118 for( mm=0 ; mm < IMARR_COUNT(tar) ; mm++ ){
119 tim = IMARR_SUBIM(tar,mm) ;
120 if (!doterse) printf("%02d..%02d: %s\n", jj,jj+tim->ny-1, argv[ff+mm] ) ;
121 jj += tim->ny ;
122 }
123
124 /* create vectors from 1D files */
125
126 tvec = (double **) malloc( sizeof(double *)*nvec ) ;
127 svec = (double * ) malloc( sizeof(double )*nvec ) ;
128 for( jj=0 ; jj < nvec ; jj++ )
129 tvec[jj] = (double *) malloc( sizeof(double)*nx ) ;
130
131 kk = 0 ;
132 if( do_one ){
133 svec[0] = 1.0 / sqrt((double)nx) ;
134 for( ii=0 ; ii < nx ; ii++ ) tvec[0][ii] = 1.0 ;
135 kk = 1 ;
136 }
137
138 for( mm=0 ; mm < IMARR_COUNT(tar) ; mm++ ){
139 tim = IMARR_SUBIM(tar,mm) ;
140 far = MRI_FLOAT_PTR(tim) ;
141 for( jj=0 ; jj < tim->ny ; jj++,kk++ ){
142 for( ii=0 ; ii < nx ; ii++ ) tvec[kk][ii] = far[ii+jj*nx] ;
143 if( demean ){
144 sum = 0.0 ;
145 for( ii=0 ; ii < nx ; ii++ ) sum += tvec[kk][ii] ;
146 sum /= nx ;
147 for( ii=0 ; ii < nx ; ii++ ) tvec[kk][ii] -= sum ;
148 }
149 sum = 0.0 ;
150 for( ii=0 ; ii < nx ; ii++ ) sum += tvec[kk][ii] * tvec[kk][ii] ;
151 if( sum == 0.0 ) {
152 if (okzero) svec[kk] = 0.0;
153 else ERROR_exit("Input column %02d is all zero!",kk) ;
154 } else {
155 svec[kk] = 1.0 / sqrt(sum) ;
156 }
157 }
158 }
159 DESTROY_IMARR(tar) ;
160
161 /* normalize vectors? (for ordinary correlation) */
162
163 if( !docov ){
164 for( kk=0 ; kk < nvec ; kk++ ){
165 sum = svec[kk] ;
166 for( ii=0 ; ii < nx ; ii++ ) tvec[kk][ii] *= sum ;
167 }
168 }
169
170 switch(docov){
171 default:
172 case 3: matname = "Spearman" ; break ;
173 case 2: matname = "InnerProduct" ; break ;
174 case 1: matname = "Covariance" ; break ;
175 case 0: matname = "Correlation" ; break ;
176 }
177
178 /* create matrix from dot product of vectors */
179
180 amat = (double *) calloc( sizeof(double) , nvec*nvec ) ;
181
182 if( docov != 3 ){
183 for( kk=0 ; kk < nvec ; kk++ ){
184 for( jj=0 ; jj <= kk ; jj++ ){
185 sum = 0.0 ;
186 for( ii=0 ; ii < nx ; ii++ ) sum += tvec[jj][ii] * tvec[kk][ii] ;
187 amat[jj+nvec*kk] = sum ;
188 if( jj < kk ) amat[kk+nvec*jj] = sum ;
189 }
190 }
191 } else { /* Spearman */
192 for( kk=0 ; kk < nvec ; kk++ ){
193 for( jj=0 ; jj <= kk ; jj++ ){
194 amat[jj+nvec*kk] = THD_spearman_corr_dble( nx , tvec[jj] , tvec[kk] ) ;
195 if( jj < kk ) amat[kk+nvec*jj] = amat[jj+nvec*kk] ;
196 }
197 }
198 }
199
200 /* normalize */
201 if (docov==1) {
202 for( kk=0 ; kk < nvec ; kk++ ){
203 for( jj=0 ; jj <= kk ; jj++ ){
204 sum = amat[jj+nvec*kk] / (double) (nx-1);
205 amat[jj+nvec*kk] = sum;
206 if( jj < kk ) amat[kk+nvec*jj] = sum ;
207 }
208 }
209 }
210
211 /* print matrix out */
212
213 if (!doterse) {
214 printf("\n"
215 "++ %s Matrix:\n ",matname) ;
216 for( jj=0 ; jj < nvec ; jj++ ) printf(" %02d ",jj) ;
217 printf("\n ") ;
218 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
219 printf("\n") ;
220 }
221 for( kk=0 ; kk < nvec ; kk++ ){
222 if (!doterse) printf("%02d:",kk) ;
223 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",amat[jj+kk*nvec]) ;
224 printf("\n") ;
225 }
226
227 if (doterse) exit(0) ; /* au revoir */
228
229 /* compute eigendecomposition */
230
231 eval = (double *) malloc( sizeof(double)*nvec ) ;
232 symeig_double( nvec , amat , eval ) ;
233
234 printf("\n"
235 "++ Eigensolution of %s Matrix:\n " , matname ) ;
236 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",eval[jj]) ;
237 printf("\n ") ;
238 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
239 printf("\n") ;
240 for( kk=0 ; kk < nvec ; kk++ ){
241 printf("%02d:",kk) ;
242 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",amat[kk+jj*nvec]) ;
243 printf("\n") ;
244 }
245
246 /* compute matrix inverse */
247 if ( eval[0]/eval[nvec-1] < 1.0e-10) {
248 printf("\n"
249 "-- WARNING: Matrix is near singular,\n"
250 " rubbish likely for inverses ahead.\n");
251 }
252 for( jj=0 ; jj < nvec ; jj++ ) eval[jj] = 1.0 / eval[jj] ;
253
254 bmat = (double *) calloc( sizeof(double) , nvec*nvec ) ;
255
256 for( ii=0 ; ii < nvec ; ii++ ){
257 for( jj=0 ; jj < nvec ; jj++ ){
258 sum = 0.0 ;
259 for( kk=0 ; kk < nvec ; kk++ )
260 sum += amat[ii+kk*nvec] * amat[jj+kk*nvec] * eval[kk] ;
261 bmat[ii+jj*nvec] = sum ;
262 }
263 }
264
265 printf("\n") ;
266 printf("++ %s Matrix Inverse:\n " , matname ) ;
267 for( jj=0 ; jj < nvec ; jj++ ) printf(" %02d ",jj) ;
268 printf("\n ") ;
269 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
270 printf("\n") ;
271 for( kk=0 ; kk < nvec ; kk++ ){
272 printf("%02d:",kk) ;
273 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",bmat[jj+kk*nvec]) ;
274 printf("\n") ;
275 }
276
277 /* square roots of diagonals of the above */
278
279 printf("\n") ;
280 printf("++ %s sqrt(diagonal)\n ",matname) ;
281 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",sqrt(bmat[jj+jj*nvec])) ;
282 printf("\n") ;
283
284 /* normalize matrix inverse */
285
286 for( ii=0 ; ii < nvec ; ii++ ){
287 for( jj=0 ; jj < nvec ; jj++ ){
288 sum = bmat[ii+ii*nvec] * bmat[jj+jj*nvec] ;
289 if( sum > 0.0 )
290 amat[ii+jj*nvec] = bmat[ii+jj*nvec] / sqrt(sum) ;
291 else
292 amat[ii+jj*nvec] = 0.0 ;
293 }
294 }
295
296 printf("\n") ;
297 printf("++ %s Matrix Inverse Normalized:\n " , matname ) ;
298 for( jj=0 ; jj < nvec ; jj++ ) printf(" %02d ",jj) ;
299 printf("\n ") ;
300 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
301 printf("\n") ;
302 for( kk=0 ; kk < nvec ; kk++ ){
303 printf("%02d:",kk) ;
304 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",amat[jj+kk*nvec]) ;
305 printf("\n") ;
306 }
307
308 /* done */
309
310 exit(0) ;
311 }
312