1 #include "mrilib.h"
2
3 /*-----------------------------------------------------------------*/
4
main(int argc,char * argv[])5 int main( int argc , char *argv[] )
6 {
7 int iarg , ii,jj,kk,mm , nvec , do_one=0 , nx=0,ny , ff ;
8 MRI_IMAGE *tim ;
9 MRI_IMARR *tar ;
10 double *amat , *sval , *umat , *vmat , smax,del,sum ;
11 float *far ;
12 int do_cond=0 ; /* 08 Nov 2004 */
13 int do_sing=0 ;
14 int do_1Dll=0 ; /* 05 Jan 2005 */
15 int do_vmean=0 , do_vnorm=0 ; /* 25 Feb 2008 */
16 int pall=1 , nev=0 ;
17 double nevper=0.0 , svsum ;
18
19 /*---------- help? ----------*/
20
21 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ){
22 printf(
23 "Usage: 1dsvd [options] 1Dfile 1Dfile ...\n"
24 "- Computes SVD of the matrix formed by the 1D file(s).\n"
25 "- Output appears on stdout; to save it, use '>' redirection.\n"
26 "\n"
27 "OPTIONS:\n"
28 " -one = Make 1st vector be all 1's.\n"
29 " -vmean = Remove mean from each vector (can't be used with -one).\n"
30 " -vnorm = Make L2-norm of each vector = 1 before SVD.\n"
31 " * The above 2 options mirror those in 3dpc.\n"
32 " -cond = Only print condition number (ratio of extremes)\n"
33 " -sing = Only print singular values\n"
34 " * To compare the singular values from 1dsvd with those from\n"
35 " 3dDeconvolve you must use the -vnorm option with 1dsvd.\n"
36 " For example, try\n"
37 " 3dDeconvolve -nodata 200 1 -polort 5 -num_stimts 1 \\\n"
38 " -stim_times 1 '1D: 30 130' 'BLOCK(50,1)' -singvals\n"
39 " 1dsvd -sing -vnorm nodata.xmat.1D\n"
40 " -sort = Sort singular values (descending) [the default]\n"
41 " -nosort = Don't bother to sort the singular values\n"
42 " -asort = Sort singular values (ascending)\n"
43 " -1Dleft = Only output left eigenvectors, in a .1D format\n"
44 " This might be useful for reducing the number of\n"
45 " columns in a design matrix. The singular values\n"
46 " are printed at the top of each vector column,\n"
47 " as a '#...' comment line.\n"
48 " -nev n = If -1Dleft is used, '-nev' specifies to output only\n"
49 " the first 'n' eigenvectors, rather than all of them.\n"
50 " * If you are a tricky person, such as Souheil, you can\n"
51 " put a '%%' after the value, and then you are saying\n"
52 " keep eigenvectors until at least n%% of the sum of\n"
53 " singular values is accounted for. In this usage,\n"
54 " 'n' must be a number less than 100; for example, to\n"
55 " reduce a matrix down to a smaller set of columns that\n"
56 " capture most of its column space, try something like\n"
57 " 1dsvd -1Dleft -nev 99%% Xorig.1D > X99.1D\n"
58 "EXAMPLE:\n"
59 " 1dsvd -vmean -vnorm -1Dleft fred.1D'[1..6]' | 1dplot -stdin\n"
60 "NOTES:\n"
61 "* Call the input n X m matrix [A] (n rows, m columns). The SVD\n"
62 " is the factorization [A] = [U] [S] [V]' ('=transpose), where\n"
63 " - [U] is an n x m matrix (whose columns are the 'Left vectors')\n"
64 " - [S] is a diagonal m x m matrix (the 'singular values')\n"
65 " - [V] is an m x m matrix (whose columns are the 'Right vectors')\n"
66 "* The default output of the program is\n"
67 " - An echo of the input [A]\n"
68 " - The [U] matrix, each column headed by its singular value\n"
69 " - The [V] matrix, each column headed by its singular value\n"
70 " (please note that [V] is output, not [V]')\n"
71 " - The pseudo-inverse of [A]\n"
72 "* This program was written simply for some testing purposes,\n"
73 " but is distributed with AFNI because it might be useful-ish.\n"
74 "* Recall that you can transpose a .1D file on input by putting\n"
75 " an escaped ' character after the filename. For example,\n"
76 " 1dsvd fred.1D\\'\n"
77 " You can use this feature to get around the fact that there\n"
78 " is no '-1Dright' option. If you understand.\n"
79 "* For more information on the SVD, you can start at\n"
80 " http://en.wikipedia.org/wiki/Singular_value_decomposition\n"
81 "* Author: Zhark the Algebraical (Linear).\n"
82 ) ;
83 PRINT_COMPILE_DATE ; exit(0) ;
84 }
85
86 /*---------- options ----------*/
87
88 iarg = 1 ; nvec = 0 ; set_svd_sort(-1) ;
89 while( iarg < argc && argv[iarg][0] == '-' ){
90
91 if( strcasecmp(argv[iarg],"-nev") == 0 ){ /* 04 Mar 2009 */
92 char *ppp ; double val ;
93 val = strtod(argv[++iarg],&ppp) ; iarg++ ;
94 if( *ppp == '%' ){
95 nevper = val ;
96 if( nevper < 1.0 ) nevper = 1.0 ;
97 if( nevper > 100.0 ) nevper = 100.0 ;
98 nev = -1 ;
99 } else {
100 nev = (int)val ; if( nev <= 0 ) nev = 1 ;
101 }
102 continue ;
103 }
104
105 if( strcasecmp(argv[iarg],"-sort") == 0 ){
106 #if 0
107 INFO_message("1dsvd: '-sort' option is now the default");
108 #endif
109 set_svd_sort(-1) ; iarg++ ; continue ;
110 }
111 if( strcasecmp(argv[iarg],"-nosort") == 0 ){
112 set_svd_sort(0) ; iarg++ ; continue ;
113 }
114 if( strcasecmp(argv[iarg],"-asort") == 0 ){
115 set_svd_sort(1) ; iarg++ ; continue ;
116 }
117
118 if( strcasecmp(argv[iarg],"-1Dright") == 0 ){
119 ERROR_exit("-1Dright has been replaced by -1Dleft") ;
120 }
121
122 if( strcasecmp(argv[iarg],"-1Dleft") == 0 ){
123 pall = 0 ; do_1Dll = 1 ; iarg++ ; continue ;
124 }
125
126 if( strcasecmp(argv[iarg],"-one") == 0 ){
127 do_one = 1 ; iarg++ ; continue ;
128 }
129
130 if( strcasecmp(argv[iarg],"-vmean") == 0 ){
131 do_vmean = 1 ; iarg++ ; continue ;
132 }
133 if( strcasecmp(argv[iarg],"-vnorm") == 0 ){
134 do_vnorm = 1 ; iarg++ ; continue ;
135 }
136
137 if( strcasecmp(argv[iarg],"-cond") == 0 ){
138 pall = 0 ; do_cond = 1 ; iarg++ ; continue ;
139 }
140
141 if( strcasecmp(argv[iarg],"-sing") == 0 ){
142 pall = 0 ; do_sing = 1 ; iarg++ ; continue ;
143 }
144
145 ERROR_exit("Unknown option: %s",argv[iarg]) ;
146 }
147
148 if( iarg == argc ) ERROR_exit("No 1D files on command line!?\n") ;
149
150 if( do_vmean && do_one )
151 ERROR_exit("Can't use -vmean and -one in the same run!") ;
152
153 /* input 1D files */
154
155 ff = iarg ;
156 INIT_IMARR(tar) ; if( do_one ) nvec = 1 ;
157 for( ; iarg < argc ; iarg++ ){
158 tim = mri_read_1D( argv[iarg] ) ;
159 if( tim == NULL ) ERROR_exit("Can't read 1D file %s\n",argv[iarg]) ;
160 if( nx == 0 )
161 nx = tim->nx ;
162 else if( tim->nx != nx )
163 ERROR_exit("1D file %s doesn't match first file in length!",argv[iarg]);
164 nvec += tim->ny ;
165 ADDTO_IMARR(tar,tim) ;
166 }
167
168 if( !do_1Dll || nev == 0 || nev > nvec ) nev = nvec ;
169
170 if( pall ){
171 printf("\n") ;
172 printf("++ 1dsvd input vectors:\n") ;
173 jj = 0 ;
174 if( do_one ){
175 printf("00..00: all ones\n") ; jj = 1 ;
176 }
177 for( mm=0 ; mm < IMARR_COUNT(tar) ; mm++ ){
178 tim = IMARR_SUBIM(tar,mm) ;
179 printf("%02d..%02d: %s\n", jj,jj+tim->ny-1, argv[ff+mm] ) ;
180 jj += tim->ny ;
181 }
182 }
183
184 /* create matrix from 1D files */
185
186 #define A(i,j) amat[(i)+(j)*nx] /* nx X nvec matrix */
187 #define U(i,j) umat[(i)+(j)*nx] /* ditto */
188 #define V(i,j) vmat[(i)+(j)*nvec] /* nvec X nvec matrix */
189 #define X(i,j) amat[(i)+(j)*nvec] /* nvec X nx matrix */
190
191 amat = (double *)malloc( sizeof(double)*nx*nvec ) ;
192 umat = (double *)malloc( sizeof(double)*nx*nvec ) ;
193 vmat = (double *)malloc( sizeof(double)*nvec*nvec ) ;
194 sval = (double *)malloc( sizeof(double)*nvec ) ;
195
196 kk = 0 ;
197 if( do_one ){
198 for( ii=0 ; ii < nx ; ii++ ) A(ii,kk) = 1.0 ;
199 kk++ ;
200 }
201
202 for( mm=0 ; mm < IMARR_COUNT(tar) ; mm++ ){
203 tim = IMARR_SUBIM(tar,mm) ;
204 far = MRI_FLOAT_PTR(tim) ;
205 for( jj=0 ; jj < tim->ny ; jj++,kk++ ){
206 for( ii=0 ; ii < nx ; ii++ ) A(ii,kk) = far[ii+jj*nx] ;
207 }
208 }
209 DESTROY_IMARR(tar) ; /* done with input data images */
210
211 if( do_vmean ){ /* 25 Feb 2008 */
212 double sum ;
213 for( jj=0 ; jj < nvec ; jj++ ){
214 sum = 0.0 ;
215 for( ii=0 ; ii < nx ; ii++ ) sum += A(ii,jj) ;
216 sum /= nx ;
217 for( ii=0 ; ii < nx ; ii++ ) A(ii,jj) -= sum ;
218 }
219 }
220 if( do_vnorm ){ /* 25 Feb 2008 */
221 double sum ;
222 for( jj=0 ; jj < nvec ; jj++ ){
223 sum = 0.0 ;
224 for( ii=0 ; ii < nx ; ii++ ) sum += A(ii,jj)*A(ii,jj) ;
225 if( sum > 0.0 ){
226 sum = 1.0 / sqrt(sum) ;
227 for( ii=0 ; ii < nx ; ii++ ) A(ii,jj) *= sum ;
228 }
229 }
230 }
231
232 /**----- the core of the program -----**/
233
234 svd_double( nx , nvec , amat , sval , umat , vmat ) ;
235
236 svsum = 0.0 ;
237 for( jj=0 ; jj < nvec ; jj++ )
238 if( sval[jj] > 0.0 ) svsum += sval[jj] ;
239
240 if( nevper > 0.0 && svsum > 0.0 ){
241 double ss=0.0 , sst=nevper*svsum/100.0 ;
242 for( jj=0 ; jj < nvec && ss < sst ; jj++ ) ss += sval[jj] ;
243 nev = jj ;
244 }
245
246 /*-- various outputs now go to stdout --*/
247
248 if( do_cond ){
249 double sbot,stop , cnum ;
250 sbot = stop = MAX(sval[0],0.0) ;
251 for( jj=1 ; jj < nvec ; jj++ ){
252 if( sval[jj] < sbot ) sbot = sval[jj] ;
253 if( sval[jj] > stop ) stop = sval[jj] ;
254 }
255 cnum = stop/sbot ;
256 if( do_1Dll ) printf("# condition number = ") ;
257 printf("%.7g\n",cnum) ;
258 }
259
260 if( do_sing && !do_1Dll ){
261 for( jj=0 ; jj < nvec ; jj++ ) printf(" %6g",sval[jj]) ;
262 printf("\n") ;
263 }
264
265 if( !pall && !do_1Dll ) exit(0) ;
266
267 if( !do_1Dll ){
268 printf("\n"
269 "++ Data vectors [A]:\n " ) ;
270 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
271 printf("\n") ;
272 for( kk=0 ; kk < nx ; kk++ ){
273 printf("%02d:",kk) ;
274 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",A(kk,jj)) ;
275 printf("\n") ;
276 }
277 }
278
279 if( !do_1Dll ) printf("\n++ Left Vectors [U]:\n " ) ;
280 else printf("# Left Vectors\n#") ;
281 for( jj=0 ; jj < nev ; jj++ ) printf(" %12.6g",sval[jj]) ;
282 printf("\n") ;
283 if( do_1Dll ) printf("#") ; else printf(" ") ;
284 for( jj=0 ; jj < nev ; jj++ ) printf(" ------------") ;
285 printf("\n") ;
286 for( kk=0 ; kk < nx ; kk++ ){
287 if( !do_1Dll) printf("%02d:",kk) ;
288 for( jj=0 ; jj < nev ; jj++ ) printf(" %12.6g",U(kk,jj)) ;
289 printf("\n") ;
290 }
291
292 if( do_1Dll ) exit(0) ;
293
294 printf("\n"
295 "++ Right Vectors [V]:\n " ) ;
296 for( jj=0 ; jj < nvec ; jj++ ) printf(" %12.6g",sval[jj]) ;
297 printf("\n ") ;
298 for( jj=0 ; jj < nvec ; jj++ ) printf(" ------------") ;
299 printf("\n") ;
300 for( kk=0 ; kk < nvec ; kk++ ){
301 printf("%02d:",kk) ;
302 for( jj=0 ; jj < nvec ; jj++ ) printf(" %12.6g",V(kk,jj)) ;
303 printf("\n") ;
304 }
305
306 smax = sval[0] ;
307 for( ii=1 ; ii < nvec ; ii++ )
308 if( sval[ii] > smax ) smax = sval[ii] ;
309
310 del = 1.e-12 * smax*smax ;
311 for( ii=0 ; ii < nvec ; ii++ )
312 sval[ii] = sval[ii] / ( sval[ii]*sval[ii] + del ) ;
313
314 /* create pseudo-inverse */
315
316 for( ii=0 ; ii < nvec ; ii++ ){
317 for( jj=0 ; jj < nx ; jj++ ){
318 sum = 0.0 ;
319 for( kk=0 ; kk < nvec ; kk++ )
320 sum += sval[kk] * V(ii,kk) * U(jj,kk) ;
321 X(ii,jj) = sum ;
322 }
323 }
324
325 printf("\n"
326 "++ Pseudo-inverse:\n " ) ;
327 for( jj=0 ; jj < nx ; jj++ ) printf(" ------------") ;
328 printf("\n") ;
329 for( kk=0 ; kk < nvec ; kk++ ){
330 printf("%02d:",kk) ;
331 for( jj=0 ; jj < nx ; jj++ ) printf(" %12.6g",X(kk,jj)) ;
332 printf("\n") ;
333 }
334
335 /* done */
336
337 exit(0) ;
338 }
339