Doxygen Source Code Documentation
1dsvd.c File Reference
#include "mrilib.h"Go to the source code of this file.
Defines | |
| #define | A(i, j) amat[(i)+(j)*nx] |
| #define | U(i, j) umat[(i)+(j)*nx] |
| #define | V(i, j) vmat[(i)+(j)*nvec] |
| #define | X(i, j) amat[(i)+(j)*nvec] |
Functions | |
| int | dbc (const void *ap, const void *bp) |
| int | main (int argc, char *argv[]) |
Define Documentation
|
|
|
|
|
|
|
|
|
|
|
|
Function Documentation
|
||||||||||||
|
Definition at line 3 of file 1dsvd.c. References a. Referenced by main().
|
|
||||||||||||
|
convert three sub-briks to a raw dataset with consecutive triplets Definition at line 13 of file 1dsvd.c. References ADDTO_IMARR, argc, dbc(), DESTROY_IMARR, far, IMARR_COUNT, IMARR_SUBIM, INIT_IMARR, malloc, MAX, MRI_FLOAT_PTR, mri_read_1D(), MRI_IMAGE::nx, MRI_IMAGE::ny, and svd_double().
00014 {
00015 int iarg , ii,jj,kk,mm , nvec , do_one=0 , nx=0,ny , ff ;
00016 MRI_IMAGE *tim ;
00017 MRI_IMARR *tar ;
00018 double *amat , *sval , *umat , *vmat , smax,del,sum ;
00019 float *far ;
00020 int do_cond=0 ; /* 08 Nov 2004 */
00021 int do_sing=0 ;
00022 int do_1Drr=0 ; /* 05 Jan 2005 */
00023 int pall=1 ;
00024
00025 /* help? */
00026
00027 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00028 printf("Usage: 1dsvd [options] 1Dfile 1Dfile ...\n"
00029 "- Computes SVD of the matrix formed by the 1D file(s).\n"
00030 "- Output appears on stdout; to save it, use '>' redirection.\n"
00031 "\n"
00032 "Options:\n"
00033 " -one = Make 1st vector be all 1's.\n"
00034 " -cond = Only print condition number (ratio of extremes)\n"
00035 " -sing = Only print singular values\n"
00036 " -1Dright = Only output right eigenvectors, in a .1D format\n"
00037 " This can be useful for reducing the number of\n"
00038 " columns in a design matrix. The singular values\n"
00039 " are printed at the top of each vector column.\n"
00040 ) ;
00041 exit(0) ;
00042 }
00043
00044 /* options */
00045
00046 iarg = 1 ; nvec = 0 ;
00047 while( iarg < argc && argv[iarg][0] == '-' ){
00048
00049 if( strcmp(argv[iarg],"-1Dright") == 0 ){
00050 pall = 0 ; do_1Drr = 1 ; iarg++ ; continue ;
00051 }
00052
00053 if( strcmp(argv[iarg],"-one") == 0 ){
00054 do_one = 1 ; iarg++ ; continue ;
00055 }
00056
00057 if( strcmp(argv[iarg],"-cond") == 0 ){
00058 pall = 0 ; do_cond = 1 ; iarg++ ; continue ;
00059 }
00060
00061 if( strcmp(argv[iarg],"-sing") == 0 ){
00062 pall = 0 ; do_sing = 1 ; iarg++ ; continue ;
00063 }
00064
00065 fprintf(stderr,"** Unknown option: %s\n",argv[iarg]); exit(1);
00066 }
00067
00068 if( iarg == argc ){
00069 fprintf(stderr,"** No 1D files on command line!?\n"); exit(1);
00070 }
00071
00072 /* input 1D files */
00073
00074 ff = iarg ;
00075 INIT_IMARR(tar) ; if( do_one ) nvec = 1 ;
00076 for( ; iarg < argc ; iarg++ ){
00077 tim = mri_read_1D( argv[iarg] ) ;
00078 if( tim == NULL ){
00079 fprintf(stderr,"** Can't read 1D file %s\n",argv[iarg]); exit(1);
00080 }
00081 if( nx == 0 ){
00082 nx = tim->nx ;
00083 } else if( tim->nx != nx ){
00084 fprintf(stderr,"** 1D file %s doesn't match first file in length!\n",
00085 argv[iarg]); exit(1);
00086 }
00087 nvec += tim->ny ;
00088 ADDTO_IMARR(tar,tim) ;
00089 }
00090
00091 if( pall ){
00092 printf("\n") ;
00093 printf("++ 1dsvd input vectors:\n") ;
00094 jj = 0 ;
00095 if( do_one ){
00096 printf("00..00: all ones\n") ; jj = 1 ;
00097 }
00098 for( mm=0 ; mm < IMARR_COUNT(tar) ; mm++ ){
00099 tim = IMARR_SUBIM(tar,mm) ;
00100 printf("%02d..%02d: %s\n", jj,jj+tim->ny-1, argv[ff+mm] ) ;
00101 jj += tim->ny ;
00102 }
00103 }
00104
00105 /* create matrix from 1D files */
00106
00107 #define A(i,j) amat[(i)+(j)*nx] /* nx X nvec matrix */
00108 #define U(i,j) umat[(i)+(j)*nx] /* ditto */
00109 #define V(i,j) vmat[(i)+(j)*nvec] /* nvec X nvec matrix */
00110 #define X(i,j) amat[(i)+(j)*nvec] /* nvec X nx matrix */
00111
00112 amat = (double *)malloc( sizeof(double)*nx*nvec ) ;
00113 umat = (double *)malloc( sizeof(double)*nx*nvec ) ;
00114 vmat = (double *)malloc( sizeof(double)*nvec*nvec ) ;
00115 sval = (double *)malloc( sizeof(double)*nvec ) ;
00116
00117 kk = 0 ;
00118 if( do_one ){
00119 for( ii=0 ; ii < nx ; ii++ ) A(ii,kk) = 1.0l ;
00120 kk++ ;
00121 }
00122
00123 for( mm=0 ; mm < IMARR_COUNT(tar) ; mm++ ){
00124 tim = IMARR_SUBIM(tar,mm) ;
00125 far = MRI_FLOAT_PTR(tim) ;
00126 for( jj=0 ; jj < tim->ny ; jj++,kk++ ){
00127 for( ii=0 ; ii < nx ; ii++ ) A(ii,kk) = far[ii+jj*nx] ;
00128 }
00129 }
00130 DESTROY_IMARR(tar) ;
00131
00132 svd_double( nx , nvec , amat , sval , umat , vmat ) ;
00133
00134 if( do_cond ){
00135 double sbot,stop , cnum ;
00136 sbot = stop = MAX(sval[0],0.0) ;
00137 for( jj=1 ; jj < nvec ; jj++ ){
00138 if( sval[jj] < sbot ) sbot = sval[jj] ;
00139 if( sval[jj] > stop ) stop = sval[jj] ;
00140 }
00141 cnum = stop/sbot ;
00142 if( do_1Drr ) printf("# condition number = ") ;
00143 printf("%.7g\n",cnum) ;
00144 }
00145
00146 if( do_sing && !do_1Drr ){
00147 qsort( sval , nvec , sizeof(double) , dbc ) ;
00148 for( jj=0 ; jj < nvec ; jj++ ) printf(" %6g",sval[jj]) ;
00149 printf("\n") ;
00150 }
00151
00152 if( !pall && !do_1Drr ) exit(0) ;
00153
00154 if( !do_1Drr ){
00155 printf("\n"
00156 "++ Data vectors [A]:\n " ) ;
00157 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
00158 printf("\n") ;
00159 for( kk=0 ; kk < nx ; kk++ ){
00160 printf("%02d:",kk) ;
00161 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",A(kk,jj)) ;
00162 printf("\n") ;
00163 }
00164 }
00165
00166 if( !do_1Drr ) printf("\n++ Right Vectors [U]:\n " ) ;
00167
00168 if( do_1Drr ) printf("#") ;
00169 for( jj=0 ; jj < nvec ; jj++ ) printf(" %12.5g",sval[jj]) ;
00170 printf("\n") ;
00171 if( do_1Drr ) printf("#") ; else printf(" ") ;
00172 for( jj=0 ; jj < nvec ; jj++ ) printf(" ------------") ;
00173 printf("\n") ;
00174 for( kk=0 ; kk < nx ; kk++ ){
00175 if( !do_1Drr) printf("%02d:",kk) ;
00176 for( jj=0 ; jj < nvec ; jj++ ) printf(" %12.5g",U(kk,jj)) ;
00177 printf("\n") ;
00178 }
00179
00180 if( do_1Drr ) exit(0) ;
00181
00182 printf("\n"
00183 "++ Left Vectors [V]:\n " ) ;
00184 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",sval[jj]) ;
00185 printf("\n ") ;
00186 for( jj=0 ; jj < nvec ; jj++ ) printf(" ---------") ;
00187 printf("\n") ;
00188 for( kk=0 ; kk < nvec ; kk++ ){
00189 printf("%02d:",kk) ;
00190 for( jj=0 ; jj < nvec ; jj++ ) printf(" %9.5f",V(kk,jj)) ;
00191 printf("\n") ;
00192 }
00193
00194 smax = sval[0] ;
00195 for( ii=1 ; ii < nvec ; ii++ )
00196 if( sval[ii] > smax ) smax = sval[ii] ;
00197
00198 del = 1.e-12 * smax*smax ;
00199 for( ii=0 ; ii < nvec ; ii++ )
00200 sval[ii] = sval[ii] / ( sval[ii]*sval[ii] + del ) ;
00201
00202 /* create pseudo-inverse */
00203
00204 for( ii=0 ; ii < nvec ; ii++ ){
00205 for( jj=0 ; jj < nx ; jj++ ){
00206 sum = 0.0l ;
00207 for( kk=0 ; kk < nvec ; kk++ )
00208 sum += sval[kk] * V(ii,kk) * U(jj,kk) ;
00209 X(ii,jj) = sum ;
00210 }
00211 }
00212
00213 printf("\n"
00214 "++ Pseudo-inverse:\n " ) ;
00215 for( jj=0 ; jj < nx ; jj++ ) printf(" ---------") ;
00216 printf("\n") ;
00217 for( kk=0 ; kk < nvec ; kk++ ){
00218 printf("%02d:",kk) ;
00219 for( jj=0 ; jj < nx ; jj++ ) printf(" %9.5f",X(kk,jj)) ;
00220 printf("\n") ;
00221 }
00222
00223 /* done */
00224
00225 exit(0) ;
00226 }
|
