00001 /*********************************/
00002 /* Principal Components Analysis */
00003 /*********************************/
00004 
00005 /*** Modified 1993, 1994 by RW Cox.
00006      Now reads input matrix from file, instead of data vectors.
00007 ***/
00008 
00009 /*********************************************************************/
00010 /* Principal Components Analysis or the Karhunen-Loeve expansion is a
00011    classical method for dimensionality reduction or exploratory data
00012    analysis.  One reference among many is: F. Murtagh and A. Heck,
00013    Multivariate Data Analysis, Kluwer Academic, Dordrecht, 1987.
00014 
00015    Author: F. Murtagh Phone:        + 49 89 32006298 (work) + 49 89
00016                  965307 (home) Earn/Bitnet:  fionn@dgaeso51,
00017    fim@dgaipp1s,  murtagh@stsci Span:         esomc1::fionn Internet:
00018    murtagh@scivax.stsci.edu
00019 
00020    F. Murtagh, Munich, 6 June 1989                                   */
00021 /*********************************************************************/
00022 
00023 #include <stdio.h>
00024 #include <string.h>
00025 #include <math.h>
00026 
00027 #define SIGN(a, b) ( (b) < 0 ? -fabs(a) : fabs(a) )
00028 
00029 /**************************************************************************/
00030 
00031 void vector_write( char *fff , int nvec , float *vvv )
00032 {
00033    FILE *fffile ;
00034    float scale , va ;
00035    int ii , vi ;
00036 
00037    fffile = fopen( fff , "w" ) ;
00038    if( fffile == NULL ){
00039       fprintf( stderr , "cannot open output file %s\n" , fff ) ;
00040       exit(1) ;
00041    }
00042    scale = 0.0 ;
00043    for( ii=1 ; ii <= nvec ; ii++ ){
00044       va = fabs(vvv[ii]) ;
00045       if( va > scale ) scale = va ;
00046    }
00047    if( scale > 0.0 ) scale = 10000.0 / scale ;
00048    for( ii=1 ; ii <= nvec ; ii++ ){
00049       vi = scale * vvv[ii] + 0.499 ;
00050       fprintf( fffile , "%6d\n" , vi ) ;
00051    }
00052    fclose( fffile ) ;
00053    return ;
00054 }
00055 
00056 /**************************************************************************/
00057 
00058 main(argc, argv)
00059 int argc;
00060 char *argv[];
00061 
00062 {
00063 FILE *stream;
00064 int  n, m,  i, j, k, k2;
00065 float **data, **matrix(), **symmat, **symmat2, *vector(), *evals, *interm;
00066 void free_matrix(), free_vector(), corcol(), covcol(), scpcol();
00067 void tred2(), tqli();
00068 float in_value;
00069 char option, *strncpy();
00070 
00071 int ii, jbot ;
00072 float *vsum ;
00073 double dtemp ;   /* entries on disk are doubles */
00074 
00075 int cvtype ;
00076 char prefix[64] = "pcmat" , fname[128] ;
00077 
00078 /*********************************************************************
00079    Get from command line:
00080    covmat_filename option
00081 *********************************************************************/
00082 
00083    if (argc < 3)
00084       {
00085          fprintf( stderr ,
00086        "Usage: %s covfile option [prefix] \n  where option=R,V,S\n",argv[0]);
00087          exit(1) ;
00088       }
00089 
00090    if( argc >= 4 ){
00091       strcpy( prefix , argv[3] ) ;
00092    }
00093 
00094    stream = fopen( argv[1] , "r" ) ;
00095    if( stream == NULL ){fprintf(stderr,"can't open covfile\n");exit(1);}
00096 
00097 #define COVERR(n) if(ii<(n)){fprintf(stderr,"read error\n");exit(1);}
00098 
00099    ii = fread( &n      , sizeof(int) , 1 , stream ) ; COVERR(1) ;
00100    ii = fread( &m      , sizeof(int) , 1 , stream ) ; COVERR(1) ;
00101    ii = fread( &cvtype , sizeof(int) , 1 , stream ) ; COVERR(1) ;
00102 
00103    option = argv[2][0] ;   /* Analysis option */
00104 
00105    printf("No. of rows: %d, no. of columns: %d.\n",n,m);
00106    printf("Input file: %s.\n",argv[1]);
00107 
00108    vsum = vector(m) ;
00109 
00110    for( i=1 ; i <= m ; i++ ){
00111       ii = fread( &dtemp , sizeof(double) , 1 , stream ) ; COVERR(1) ;
00112       vsum[i] = dtemp / n ;
00113    }
00114 
00115    strcpy( fname , prefix ) ;  strcat( fname , ".ave" ) ;
00116    vector_write( fname , m , vsum ) ;
00117 
00118    if( cvtype <= 1 ) exit(0) ;
00119 
00120    /********************************************************************/
00121 
00122    symmat = matrix(m, m);  /* Allocation of correlation (etc.) matrix */
00123    for( i=1 ; i <= m ; i++ ){
00124       for( j=1 ; j <= m ; j++ ){
00125          ii = fread( &dtemp , sizeof(double) , 1 , stream ) ; COVERR(1) ;
00126          symmat[i][j] = dtemp / n ;
00127       }
00128    }
00129    fclose( stream ) ;
00130 
00131    /* Look at analysis option; branch in accordance with this. */
00132 
00133    switch(option){
00134       case 'R':
00135       case 'r':
00136           printf("Analysis of correlations chosen.\n");
00137           corcol(vsum, n, m, symmat);
00138         break;
00139       case 'V':
00140       case 'v':
00141           printf("Analysis of variances-covariances chosen.\n");
00142           covcol(vsum, n, m, symmat);
00143         break;
00144       case 'S':
00145       case 's':
00146           printf("Analysis of sums-of-squares-cross-products");
00147           printf(" matrix chosen.\n");
00148           scpcol(vsum, n, m, symmat);
00149         break;
00150       default:
00151           printf("Option: %c\n",option);
00152           printf("For option, please type R, V, or S\n");
00153           printf("(upper or lower case).\n");
00154           exit(1);
00155         break;
00156    }
00157 
00158 /*********************************************************************
00159     Eigen-reduction
00160 **********************************************************************/
00161 
00162     /* Allocate storage for dummy and new vectors. */
00163 
00164     evals = vector(m);     /* Storage alloc. for vector of eigenvalues */
00165     interm = vector(m);    /* Storage alloc. for 'intermediate' vector */
00166 
00167     tred2(symmat, m, evals, interm);  /* Triangular decomposition */
00168     tqli(evals, interm, m, symmat);   /* Reduction of sym. trid. matrix */
00169 
00170     /* evals now contains the eigenvalues,
00171        columns of symmat now contain the associated eigenvectors. */
00172 
00173      printf("\nEigenvalues:\n");
00174      jbot = (m<10) ? 1 : m-9 ;
00175      for (j = m; j >= jbot; j--) {
00176        printf("%18.5f\n", evals[j]);
00177      }
00178 
00179      for( i=1 ; i <= 3 ; i++ ){
00180         for( j=1 ; j <= m ; j++ ){
00181            interm[j] = symmat[j][m-i+1] ;
00182         }
00183         sprintf( fname , "%s.pc%d" , prefix , i ) ;
00184         vector_write( fname , m , interm ) ;
00185      }
00186 
00187    exit(0) ;
00188 }
00189 
00190 /**  Correlation matrix: creation  ***********************************/
00191 
00192 void corcol(vsum, n, m, symmat)
00193 float *vsum, **symmat;
00194 int n, m;
00195 {
00196 float x, *mean, *stddev, *vector();
00197 int i, j, j1, j2;
00198 void covcol() ;
00199 
00200 covcol( vsum, n , m , symmat ) ;  /* form covariance matrix first */
00201 
00202 /* Allocate storage for mean and std. dev. vectors */
00203 
00204    stddev = vector(m);
00205 
00206    for (j = 1; j <= m; j++) stddev[j] = sqrt( symmat[j][j] ) ;
00207 
00208    for( i=1 ; i <= m ; i++){
00209       for( j=1 ; j <= m ; j++ ){
00210          symmat[i][j] /= (stddev[i]*stddev[j]) ;
00211       }
00212    }
00213    return;
00214 }
00215 
00216 /**  Covariance matrix: creation  *****************************/
00217 
00218 void covcol(vsum, n, m, symmat)
00219 float *vsum, **symmat;
00220 int n, m;
00221 {
00222 float *mean, *vector();
00223 int i, j, j1, j2;
00224 
00225    for( i=1 ; i <= m ; i++ ){
00226       for( j=1 ; j <= m ; j++ ){
00227          symmat[i][j] -= vsum[i] * vsum[j] ;
00228       }
00229    }
00230    return;
00231 }
00232 
00233 /**  Sums-of-squares-and-cross-products matrix: creation  **************/
00234 
00235 void scpcol(data, n, m, symmat)
00236 float **data, **symmat;
00237 int n, m;
00238 /* Create m * m sums-of-cross-products matrix from n * m data matrix. */
00239 {
00240 int i, j1, j2;
00241 
00242 return;
00243 
00244 }
00245 
00246 /**  Error handler  **************************************************/
00247 
00248 void erhand(err_msg)
00249 char err_msg[];
00250 /* Error handler */
00251 {
00252     fprintf(stderr,"Run-time error:\n");
00253     fprintf(stderr,"%s\n", err_msg);
00254     fprintf(stderr,"Exiting to system.\n");
00255     exit(1);
00256 }
00257 
00258 /**  Allocation of vector storage  ***********************************/
00259 
00260 float *vector(n)
00261 int n;
00262 /* Allocates a float vector with range [1..n]. */
00263 {
00264 
00265     float *v;
00266 
00267     v = (float *) malloc ((unsigned) n*sizeof(float));
00268     if (!v) erhand("Allocation failure in vector().");
00269     return v-1;
00270 
00271 }
00272 
00273 /**  Allocation of float matrix storage  *****************************/
00274 
00275 float **matrix(n,m)
00276 int n, m;
00277 /* Allocate a float matrix with range [1..n][1..m]. */
00278 {
00279     int i;
00280     float **mat;
00281 
00282     /* Allocate pointers to rows. */
00283     mat = (float **) malloc((unsigned) (n)*sizeof(float*));
00284     if (!mat) erhand("Allocation failure 1 in matrix().");
00285     mat -= 1;
00286 
00287     /* Allocate rows and set pointers to them. */
00288     for (i = 1; i <= n; i++)
00289         {
00290         mat[i] = (float *) malloc((unsigned) (m)*sizeof(float));
00291         if (!mat[i]) erhand("Allocation failure 2 in matrix().");
00292         mat[i] -= 1;
00293         }
00294 
00295      /* Return pointer to array of pointers to rows. */
00296      return mat;
00297 
00298 }
00299 
00300 /**  Deallocate vector storage  *********************************/
00301 
00302 void free_vector(v,n)
00303 float *v;
00304 int n;
00305 /* Free a float vector allocated by vector(). */
00306 {
00307    free((char*) (v+1));
00308 }
00309 
00310 /**  Deallocate float matrix storage  ***************************/
00311 
00312 void free_matrix(mat,n,m)
00313 float **mat;
00314 int n, m;
00315 /* Free a float matrix allocated by matrix(). */
00316 {
00317    int i;
00318 
00319    for (i = n; i >= 1; i--)
00320        {
00321        free ((char*) (mat[i]+1));
00322        }
00323    free ((char*) (mat+1));
00324 }
00325 
00326 /**  Reduce a real, symmetric matrix to a symmetric, tridiag. matrix. */
00327 
00328 void tred2(a, n, d, e)
00329 float **a, *d, *e;
00330 /* float **a, d[], e[]; */
00331 int n;
00332 /* Householder reduction of matrix a to tridiagonal form.
00333    Algorithm: Martin et al., Num. Math. 11, 181-195, 1968.
00334    Ref: Smith et al., Matrix Eigensystem Routines -- EISPACK Guide
00335         Springer-Verlag, 1976, pp. 489-494.
00336         W H Press et al., Numerical Recipes in C, Cambridge U P,
00337         1988, pp. 373-374.  */
00338 {
00339 int l, k, j, i;
00340 float scale, hh, h, g, f;
00341 
00342 for (i = n; i >= 2; i--)
00343     {
00344     l = i - 1;
00345     h = scale = 0.0;
00346     if (l > 1)
00347        {
00348        for (k = 1; k <= l; k++)
00349            scale += fabs(a[i][k]);
00350        if (scale == 0.0)
00351           e[i] = a[i][l];
00352        else
00353           {
00354           for (k = 1; k <= l; k++)
00355               {
00356               a[i][k] /= scale;
00357               h += a[i][k] * a[i][k];
00358               }
00359           f = a[i][l];
00360           g = f>0 ? -sqrt(h) : sqrt(h);
00361           e[i] = scale * g;
00362           h -= f * g;
00363           a[i][l] = f - g;
00364           f = 0.0;
00365           for (j = 1; j <= l; j++)
00366               {
00367               a[j][i] = a[i][j]/h;
00368               g = 0.0;
00369               for (k = 1; k <= j; k++)
00370                   g += a[j][k] * a[i][k];
00371               for (k = j+1; k <= l; k++)
00372                   g += a[k][j] * a[i][k];
00373               e[j] = g / h;
00374               f += e[j] * a[i][j];
00375               }
00376           hh = f / (h + h);
00377           for (j = 1; j <= l; j++)
00378               {
00379               f = a[i][j];
00380               e[j] = g = e[j] - hh * f;
00381               for (k = 1; k <= j; k++)
00382                   a[j][k] -= (f * e[k] + g * a[i][k]);
00383               }
00384          }
00385     }
00386     else
00387         e[i] = a[i][l];
00388     d[i] = h;
00389     }
00390 d[1] = 0.0;
00391 e[1] = 0.0;
00392 for (i = 1; i <= n; i++)
00393     {
00394     l = i - 1;
00395     if (d[i])
00396        {
00397        for (j = 1; j <= l; j++)
00398            {
00399            g = 0.0;
00400            for (k = 1; k <= l; k++)
00401                g += a[i][k] * a[k][j];
00402            for (k = 1; k <= l; k++)
00403                a[k][j] -= g * a[k][i];
00404            }
00405        }
00406        d[i] = a[i][i];
00407        a[i][i] = 1.0;
00408        for (j = 1; j <= l; j++)
00409            a[j][i] = a[i][j] = 0.0;
00410     }
00411 }
00412 
00413 /**  Tridiagonal QL algorithm -- Implicit  **********************/
00414 
00415 void tqli(d, e, n, z)
00416 float d[], e[], **z;
00417 int n;
00418 {
00419 int m, l, iter, i, k;
00420 float s, r, p, g, f, dd, c, b;
00421 void erhand();
00422 
00423 for (i = 2; i <= n; i++)
00424     e[i-1] = e[i];
00425 e[n] = 0.0;
00426 for (l = 1; l <= n; l++)
00427     {
00428     iter = 0;
00429     do
00430       {
00431       for (m = l; m <= n-1; m++)
00432           {
00433           dd = fabs(d[m]) + fabs(d[m+1]);
00434           if (fabs(e[m]) + dd == dd) break;
00435           }
00436           if (m != l)
00437              {
00438              if (iter++ == 30) erhand("No convergence in TLQI.");
00439              g = (d[l+1] - d[l]) / (2.0 * e[l]);
00440              r = sqrt((g * g) + 1.0);
00441              g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
00442              s = c = 1.0;
00443              p = 0.0;
00444              for (i = m-1; i >= l; i--)
00445                  {
00446                  f = s * e[i];
00447                  b = c * e[i];
00448                  if (fabs(f) >= fabs(g))
00449                     {
00450                     c = g / f;
00451                     r = sqrt((c * c) + 1.0);
00452                     e[i+1] = f * r;
00453                     c *= (s = 1.0/r);
00454                     }
00455                  else
00456                     {
00457                     s = f / g;
00458                     r = sqrt((s * s) + 1.0);
00459                     e[i+1] = g * r;
00460                     s *= (c = 1.0/r);
00461                     }
00462                  g = d[i+1] - p;
00463                  r = (d[i] - g) * s + 2.0 * c * b;
00464                  p = s * r;
00465                  d[i+1] = g + p;
00466                  g = c * r - b;
00467                  for (k = 1; k <= n; k++)
00468                      {
00469                      f = z[k][i+1];
00470                      z[k][i+1] = s * z[k][i] + c * f;
00471                      z[k][i] = c * z[k][i] - s * f;
00472                      }
00473                  }
00474                  d[l] = d[l] - p;
00475                  e[l] = g;
00476                  e[m] = 0.0;
00477              }
00478           }  while (m != l);
00479       }
00480  }
00481