source: trunk/src/PCA.cc @ 13

#include "PCA.h"


using namespace std;
using namespace thep_cpp_tools;

PCA::PCA() : process_( false ), explained_calc_( false )
{ 
}

PCA::PCA( const thep_gsl_api::matrix& A ) : A_( A ), 
              process_( false ), 
              explained_calc_( false )
{
}

void PCA::process()
{
  process_ = true;
  // Row-center the data matrix
  thep_gsl_api::matrix A_center( A_.rows(), A_.cols() );
  this->row_center( A_center );
  

  // Single value decompose the data matrix
  auto_ptr<thep_cpp_tools::SVD> pSVD( new thep_cpp_tools::SVD( A_center ) );
  pSVD->process();
  thep_gsl_api::matrix U = pSVD->get_u(); 
  thep_gsl_api::matrix V = pSVD->get_v();

  // Read the eigenvectors and eigenvalues
  eigenvectors_ = U.transpose();
  eigenvalues_ = pSVD->get_s_vec();

  // T
  for( size_t i = 0; i < eigenvalues_.size(); ++i )
    {
     eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
    }

  eigenvalues_ /= A_center.rows();  
  
  // Sort the eigenvectors in order of eigenvalues
  // Simple (not efficient) algorithm that always 
  // make sure that the i:th element is in its correct 
  // position (N element --> Ordo( N*N ))
  for( size_t i = 0; i < eigenvalues_.size(); ++i )
    {
     for( size_t j = i + 1; j < eigenvalues_.size(); ++j )
       {
  if( eigenvalues_[ j ] > eigenvalues_[ i ] ) 
    {
     swap( eigenvalues_[ i ], eigenvalues_[ j ] ); 
     eigenvectors_.swap_rows( i, j );
    }
       }
    }
}


// This function will row-center the matrix A_,
// that is, A_ = A_ - M, where M is a matrix
// with the meanvalues of each row
void PCA::row_center( thep_gsl_api::matrix& A_center )
{
  meanvalues_ = thep_gsl_api::vector( A_.rows() );
  thep_gsl_api::vector A_row_sum = A_.row_sum();
  for( size_t i = 0; i < A_center.rows(); ++i ) 
    {
     meanvalues_[ i ] = A_row_sum.get( i ) / static_cast<double>( A_.cols() );
     for( size_t j = 0; j < A_center.cols(); ++j )
       {
  A_center.set(  i, j, A_.get( i , j ) - meanvalues_[ i ] );
       }
    }
}


thep_gsl_api::matrix PCA::projection( const thep_gsl_api::matrix& 
              samples ) const
{
  assert( process_ );
  const size_t Ncol = samples.cols();
  const size_t Nrow = samples.rows();
  thep_gsl_api::matrix projs( Ncol, Ncol );
  
  for( size_t j = 0; j < Ncol; ++j ) 
    {
     thep_gsl_api::matrix temp = samples.col( j );   
     thep_gsl_api::vector centered( Nrow );

     for( size_t i = 0; i < Nrow; ++i ) 
       {
  centered[ i ] = temp.get( i, 0 ) - meanvalues_[ i ];
       }

     thep_gsl_api::vector proj( eigenvectors_ * centered );
      for( size_t i = 0; i < samples.cols(); ++i ) 
       {
  projs.set( i, j,  proj[ i ] );
       }
    }  
  return projs;
}


void PCA::calculate_explained_intensity() 
{
  size_t DIM = eigenvalues_.size();
  explained_intensity_ = thep_gsl_api::vector( DIM );
  double sum = 0;
  for( size_t i = 0; i < DIM; ++i ) 
    {
     sum += eigenvalues_[ i ];
    }
  double exp_sum = 0;
  for( size_t i = 0; i < DIM; ++i ) 
    {
     exp_sum += eigenvalues_[ i ];
     explained_intensity_[ i ] = exp_sum / sum ;
    }  
}

double PCA::get_explained_intensity( const size_t& k )
{
  if( !explained_calc_ )
    {
     this->calculate_explained_intensity();
     explained_calc_ = true;
    }
  return explained_intensity_[ k ];
}


bool PCA::test()
{
  cout << "1. Print original matrix A = \n";//Testing row-centering of matrix\n";
  thep_gsl_api::matrix A( 3, 3 );
  for( size_t i = 0; i < 3; ++i )
    {
     for( size_t j = 0; j < 3; ++j )
       {
  A.set( i, j, sin( static_cast<double>( i + j + 2 + (i+1)*(j+1) ) ) );
       }
    }
  A_ = A;

// Print only matrix that is row-centered as PCA ( process() )
// will do that before calculation
  this->row_center( A );
  cout << A << endl;

  cout << "2. Testing PCA\n";
  this->process();

  for( size_t i = 0; i < eigenvalues_.size(); ++i )
    {
     cout << "lambda(" << i << ") = " << eigenvalues_[i] << " " 
    << "has eigenvector = " << endl;
     cout << eigenvectors_.row( i ) << endl;
    }
  
  // The eigenvalues and eigenvectors are
  // from matlab using eig(A'A)/3
  double answer_lambda[ 3 ] = { 1.8908, 0.0327, 0 };
  double answer_eigenvec[ 3 ][ 3 ] = 
    {
      { 0.4763, -0.6738, 0.5649  },
      { 0.8770, 0.3182,  -0.3599 },
      { 0.0627, 0.6669 , 0.7425  },
    };

  double ERROR_TOL = 0.0001;

  // error eigenvalues
  cout << "abs( lambda(i) - lambda_matlab(i) ) ) = \n( ";
  double max_err = 0;
  for( size_t i = 0; i < 3; ++i ) 
    {
     double err = fabs( eigenvalues_[ i ] - answer_lambda[ i ] );
     if( err > max_err ) max_err = err;
     cout << err << " ";
    }
  cout << ")\n";

  cerr << "max_error = " << max_err << ", to be compared with ERROR_TOL = " 
       << ERROR_TOL << endl; 
  if( max_err < ERROR_TOL ) 
    {
     cout << "ok, compared with precision in given matlab-references!\n";
    }
  else 
    { 
     cerr << "Not ok, error to big!\n"; 
     return false; 
    }

  // error eigenvectors
  max_err = 0;
  for( size_t i = 0; i < 3; ++i ) 
    {
     cout << "err_eigvec( " << i << " ) = \n";
     for( size_t j = 0; j < 3; ++j )
       {
  double err = fabs( fabs( eigenvectors_.get( i, j ) ) - fabs( answer_eigenvec[ i ][ j ] ) ); 
  if( err > max_err ) max_err = err;       
  cout << err << endl;
       }
     cout << ")\n\n";
    }

  cerr << "max_error = " << max_err << ", to be compared with ERROR_TOL = " 
       << ERROR_TOL << endl; 
  if( max_err < ERROR_TOL ) 
    {
     cout << "ok, compared with precision in given matlab-references!\n";
    }
  else 
    { 
     cerr << "Not ok, error to big!\n"; 
     return false; 
    }
    
  cout << "\nThe following projection is obtained:\n";
  cout << this->projection( A ) << endl;


  return true;
}