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