pca_tool.h

Go to the documentation of this file.

//
// Principal component analysis tool
//


#ifndef PCA_TOOL_H
#define PCA_TOOL_H

#include <iostream>
#include <cmath>
#include "object.h"

namespace mimas {
  class pca_tool: public object
{
  private:
    matrix_tool scores, loadings, variances;
    bool pca_performed;
  
  public:
    pca_tool(void); 
    void pca(matrix_tool& X); 
    void project(matrix_tool& Y, int n); 
};

pca_tool::pca_tool(void)
{
  pca_performed=false;
}

void pca_tool::pca(matrix_tool& X)
{
  matrix_tool u, d, v, centeredX;
  
   MMERROR( X.initialised(), exception, ,
            "pca_tool::pca - input matrix not initialised" );

  centeredX=X;

  // remove column means from centeredX
  for (int j=0; j<centeredX.getCols(); j++)
  {
    double mean=0.0;

    for (int i=0; i<centeredX.getRows(); i++)
      mean+=centeredX(i,j);

    mean/=centeredX.getRows();

    for (int i=0; i<centeredX.getRows(); i++)
      centeredX(i,j)=centeredX(i,j)-mean;

  }
  
  centeredX.svd(u,d,loadings);
  
  scores=u*d;

  variances.redim(d.getCols(),1);
  for(int i=0; i<d.getCols(); i++)
    variances(i,0)=d(i,i);

  pca_performed=true;
  
}

void pca_tool::project(matrix_tool& Y, int n)
{
  if (!pca_performed) return;

   MMERROR( n>=0 && n<loadings.getCols(), exception, ,
            "pca_tool::reduce - number of components out of range" );
  
  Y.redim(scores.getRows(),loadings.getRows());

  matrix_tool t,p;

  t.redim(scores.getRows(),n);

  for (int i=0; i<scores.getRows(); i++)
    for (int j=0; j<n; j++)
      t(i,j)=scores(i,j);
  
  p.redim(loadings.getRows(),n);

  for (int i=0; i<loadings.getRows(); i++)
    for (int j=0; j<n; j++)
      p(i,j)=loadings(i,j);


  p.save("p.dat");
  t.save("t.dat");
  
  p.transpose();

  Y=t*p;
  
}
}

#endif

---