view maths/pca/pca.c @ 483:fdaa63607c15

Untabify, indent, tidy

/*********************************/
/* Principal Components Analysis */
/*********************************/

/*********************************************************************/
/* Principal Components Analysis or the Karhunen-Loeve expansion is a
   classical method for dimensionality reduction or exploratory data
   analysis.  One reference among many is: F. Murtagh and A. Heck,
   Multivariate Data Analysis, Kluwer Academic, Dordrecht, 1987.

   Author:
   F. Murtagh
   Phone:        + 49 89 32006298 (work)
                 + 49 89 965307 (home)
   Earn/Bitnet:  fionn@dgaeso51,  fim@dgaipp1s,  murtagh@stsci
   Span:         esomc1::fionn
   Internet:     [e-mail]
   
   F. Murtagh, Munich, 6 June 1989                                   */   
/*********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "pca.h"

#define SIGN(a, b) ( (b) < 0 ? -fabs(a) : fabs(a) )

/**  Variance-covariance matrix: creation  *****************************/

/* Create m * m covariance matrix from given n * m data matrix. */
void covcol(double** data, int n, int m, double** symmat)
{
    double *mean;
    int i, j, j1, j2;

/* Allocate storage for mean vector */

    mean = (double*) malloc(m*sizeof(double));

/* Determine mean of column vectors of input data matrix */

    for (j = 0; j < m; j++)
    {
        mean[j] = 0.0;
        for (i = 0; i < n; i++)
        {
            mean[j] += data[i][j];
        }
        mean[j] /= (double)n;
    }

/*
  printf("\nMeans of column vectors:\n");
  for (j = 0; j < m; j++)  {
  printf("%12.1f",mean[j]);  }   printf("\n");
*/

/* Center the column vectors. */

    for (i = 0; i < n; i++)
    {
        for (j = 0; j < m; j++)
        {
            data[i][j] -= mean[j];
        }
    }

/* Calculate the m * m covariance matrix. */
    for (j1 = 0; j1 < m; j1++)
    {
        for (j2 = j1; j2 < m; j2++)
        {
            symmat[j1][j2] = 0.0;
            for (i = 0; i < n; i++)
            {
                symmat[j1][j2] += data[i][j1] * data[i][j2];
            }
            symmat[j2][j1] = symmat[j1][j2];
        }
    }

    free(mean);

    return;

}

/**  Error handler  **************************************************/

void erhand(char* err_msg)
{
    fprintf(stderr,"Run-time error:\n");
    fprintf(stderr,"%s\n", err_msg);
    fprintf(stderr,"Exiting to system.\n");
    exit(1);
}


/**  Reduce a real, symmetric matrix to a symmetric, tridiag. matrix. */

/* Householder reduction of matrix a to tridiagonal form.
   Algorithm: Martin et al., Num. Math. 11, 181-195, 1968.
   Ref: Smith et al., Matrix Eigensystem Routines -- EISPACK Guide
   Springer-Verlag, 1976, pp. 489-494.
   W H Press et al., Numerical Recipes in C, Cambridge U P,
   1988, pp. 373-374.  */
void tred2(double** a, int n, double* d, double* e)
{
    int l, k, j, i;
    double scale, hh, h, g, f;
        
    for (i = n-1; i >= 1; i--)
    {
        l = i - 1;
        h = scale = 0.0;
        if (l > 0)
        {
            for (k = 0; k <= l; k++)
                scale += fabs(a[i][k]);
            if (scale == 0.0)
                e[i] = a[i][l];
            else
            {
                for (k = 0; k <= l; k++)
                {
                    a[i][k] /= scale;
                    h += a[i][k] * a[i][k];
                }
                f = a[i][l];
                g = f>0 ? -sqrt(h) : sqrt(h);
                e[i] = scale * g;
                h -= f * g;
                a[i][l] = f - g;
                f = 0.0;
                for (j = 0; j <= l; j++)
                {
                    a[j][i] = a[i][j]/h;
                    g = 0.0;
                    for (k = 0; k <= j; k++)
                        g += a[j][k] * a[i][k];
                    for (k = j+1; k <= l; k++)
                        g += a[k][j] * a[i][k];
                    e[j] = g / h;
                    f += e[j] * a[i][j];
                }
                hh = f / (h + h);
                for (j = 0; j <= l; j++)
                {
                    f = a[i][j];
                    e[j] = g = e[j] - hh * f;
                    for (k = 0; k <= j; k++)
                        a[j][k] -= (f * e[k] + g * a[i][k]);
                }
            }
        }
        else
            e[i] = a[i][l];
        d[i] = h;
    }
    d[0] = 0.0;
    e[0] = 0.0;
    for (i = 0; i < n; i++)
    {
        l = i - 1;
        if (d[i])
        {
            for (j = 0; j <= l; j++)
            {
                g = 0.0;
                for (k = 0; k <= l; k++)
                    g += a[i][k] * a[k][j];
                for (k = 0; k <= l; k++)
                    a[k][j] -= g * a[k][i];
            }
        }
        d[i] = a[i][i];
        a[i][i] = 1.0;
        for (j = 0; j <= l; j++)
            a[j][i] = a[i][j] = 0.0;
    }
}

/**  Tridiagonal QL algorithm -- Implicit  **********************/

void tqli(double* d, double* e, int n, double** z)
{
    int m, l, iter, i, k;
    double s, r, p, g, f, dd, c, b;
        
    for (i = 1; i < n; i++)
        e[i-1] = e[i];
    e[n-1] = 0.0;
    for (l = 0; l < n; l++)
    {
        iter = 0;
        do
        {
            for (m = l; m < n-1; m++)
            {
                dd = fabs(d[m]) + fabs(d[m+1]);
                if (fabs(e[m]) + dd == dd) break;
            }
            if (m != l)
            {
                if (iter++ == 30) erhand("No convergence in TLQI.");
                g = (d[l+1] - d[l]) / (2.0 * e[l]);
                r = sqrt((g * g) + 1.0);
                g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
                s = c = 1.0;
                p = 0.0;
                for (i = m-1; i >= l; i--)
                {
                    f = s * e[i];
                    b = c * e[i];
                    if (fabs(f) >= fabs(g))
                    {
                        c = g / f;
                        r = sqrt((c * c) + 1.0);
                        e[i+1] = f * r;
                        c *= (s = 1.0/r);
                    }
                    else
                    {
                        s = f / g;
                        r = sqrt((s * s) + 1.0);
                        e[i+1] = g * r;
                        s *= (c = 1.0/r);
                    }
                    g = d[i+1] - p;
                    r = (d[i] - g) * s + 2.0 * c * b;
                    p = s * r;
                    d[i+1] = g + p;
                    g = c * r - b;
                    for (k = 0; k < n; k++)
                    {
                        f = z[k][i+1];
                        z[k][i+1] = s * z[k][i] + c * f;
                        z[k][i] = c * z[k][i] - s * f;
                    }
                }
                d[l] = d[l] - p;
                e[l] = g;
                e[m] = 0.0;
            }
        }  while (m != l);
    }
}

/* In place projection onto basis vectors */
void pca_project(double** data, int n, int m, int ncomponents)
{
    int  i, j, k, k2;
    double  **symmat, /* **symmat2, */ *evals, *interm;
        
    //TODO: assert ncomponents < m
        
    symmat = (double**) malloc(m*sizeof(double*));
    for (i = 0; i < m; i++)
        symmat[i] = (double*) malloc(m*sizeof(double));
                
    covcol(data, n, m, symmat);
        
    /*********************************************************************
                Eigen-reduction
    **********************************************************************/
        
    /* Allocate storage for dummy and new vectors. */
    evals = (double*) malloc(m*sizeof(double));     /* Storage alloc. for vector of eigenvalues */
    interm = (double*) malloc(m*sizeof(double));    /* Storage alloc. for 'intermediate' vector */
    //MALLOC_ARRAY(symmat2,m,m,double);    
    //for (i = 0; i < m; i++) {
    //      for (j = 0; j < m; j++) {
    //              symmat2[i][j] = symmat[i][j]; /* Needed below for col. projections */
    //      }
    //}
    tred2(symmat, m, evals, interm);  /* Triangular decomposition */
    tqli(evals, interm, m, symmat);   /* Reduction of sym. trid. matrix */
/* evals now contains the eigenvalues,
   columns of symmat now contain the associated eigenvectors. */   

/*
  printf("\nEigenvalues:\n");
  for (j = m-1; j >= 0; j--) {
  printf("%18.5f\n", evals[j]); }
  printf("\n(Eigenvalues should be strictly positive; limited\n");
  printf("precision machine arithmetic may affect this.\n");
  printf("Eigenvalues are often expressed as cumulative\n");
  printf("percentages, representing the 'percentage variance\n");
  printf("explained' by the associated axis or principal component.)\n");
        
  printf("\nEigenvectors:\n");
  printf("(First three; their definition in terms of original vbes.)\n");
  for (j = 0; j < m; j++) {
  for (i = 1; i <= 3; i++)  {
  printf("%12.4f", symmat[j][m-i]);  }
  printf("\n");  }
*/

/* Form projections of row-points on prin. components. */
/* Store in 'data', overwriting original data. */
    for (i = 0; i < n; i++) {
        for (j = 0; j < m; j++) {
            interm[j] = data[i][j]; }   /* data[i][j] will be overwritten */
        for (k = 0; k < ncomponents; k++) {
            data[i][k] = 0.0;
            for (k2 = 0; k2 < m; k2++) {
                data[i][k] += interm[k2] * symmat[k2][m-k-1]; }
        }
    }

/*      
        printf("\nProjections of row-points on first 3 prin. comps.:\n");
        for (i = 0; i < n; i++) {
        for (j = 0; j < 3; j++)  {
        printf("%12.4f", data[i][j]);  }
        printf("\n");  }
*/

/* Form projections of col.-points on first three prin. components. */
/* Store in 'symmat2', overwriting what was stored in this. */
//for (j = 0; j < m; j++) {
//       for (k = 0; k < m; k++) {
//               interm[k] = symmat2[j][k]; }  /*symmat2[j][k] will be overwritten*/
//  for (i = 0; i < 3; i++) {
//      symmat2[j][i] = 0.0;
//              for (k2 = 0; k2 < m; k2++) {
//                      symmat2[j][i] += interm[k2] * symmat[k2][m-i-1]; }
//              if (evals[m-i-1] > 0.0005)   /* Guard against zero eigenvalue */
//                      symmat2[j][i] /= sqrt(evals[m-i-1]);   /* Rescale */
//              else
//                      symmat2[j][i] = 0.0;    /* Standard kludge */
//    }
// }

/*
  printf("\nProjections of column-points on first 3 prin. comps.:\n");
  for (j = 0; j < m; j++) {
  for (k = 0; k < 3; k++)  {
  printf("%12.4f", symmat2[j][k]);  }
  printf("\n");  }
*/


    for (i = 0; i < m; i++)
        free(symmat[i]);
    free(symmat);
//FREE_ARRAY(symmat2,m);
    free(evals);
    free(interm);

}