source: trunk/yat/utility/PCA.cc @ 735

// $Id: PCA.cc 735 2007-01-06 17:11:28Z peter $

/*
  Copyright (C) The authors contributing to this file.

  This file is part of the yat library, http://lev.thep.lu.se/trac/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "PCA.h"
#include "SVD.h"
#include "utility.h"

#include <iostream>
#include <cmath>
#include <memory>

namespace theplu {
namespace yat {
namespace utility {


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


  utility::vector PCA::get_eigenvector(size_t i) const
  {
    return utility::vector(eigenvectors_,i);
  }


  double PCA::get_eigenvalue(size_t i) const
  {
    return eigenvalues_[i];
  }


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

    // Single value decompose the data matrix
    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
    pSVD->decompose();
    utility::matrix U = pSVD->U(); 
    utility::matrix V = pSVD->V();

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

    // T
    for( size_t i = 0; i < eigenvalues_.size(); ++i )
      eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
    eigenvalues_ *= 1.0/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 ] ) {
          std::swap( eigenvalues_[ i ], eigenvalues_[ j ] ); 
          eigenvectors_.swap_rows(i,j);
        }
  }


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

    // Transform into SVD friendly dimension
    A_.transpose();
    A_center.transpose();

    // Single value decompose the data matrix
    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
    pSVD->decompose();
    utility::matrix U = pSVD->U(); 
    utility::matrix V = pSVD->V();

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

    // Transform back when done with SVD!
    // (used V insted of U now for eigenvectors)
    A_.transpose();
    A_center.transpose();

    // T
    for( size_t i = 0; i < eigenvalues_.size(); ++i )
      eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
    eigenvalues_ *= 1.0/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 ] ) {
          std::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(utility::matrix& A_center)
  {
    meanvalues_ = utility::vector( A_.rows() );
    utility::vector A_row_sum(A_.rows());
    for (size_t i=0; i<A_row_sum.size(); ++i)
      A_row_sum(i)=utility::vector(A_,i).sum();
    for( size_t i = 0; i < A_center.rows(); ++i ) {
      meanvalues_[i] = A_row_sum(i) / static_cast<double>(A_.columns());
      for( size_t j = 0; j < A_center.columns(); ++j )
        A_center(i,j) = A_(i,j) - meanvalues_(i);
    }
  }


  utility::matrix PCA::projection(const utility::matrix& samples ) const
  {
    const size_t Ncol = samples.columns();
    const size_t Nrow = samples.rows();
    utility::matrix projs( Ncol, Ncol );

    utility::vector temp(samples.rows());
    for( size_t j = 0; j < Ncol; ++j ) {
      for (size_t i=0; i<Ncol; ++i )
        temp(i) = samples(i,j); 
      utility::vector centered( Nrow );
      for( size_t i = 0; i < Nrow; ++i ) 
        centered(i) = temp(i) - meanvalues_(i);
      utility::vector proj( eigenvectors_ * centered );
      for( size_t i = 0; i < Ncol; ++i ) 
        projs(i,j)=proj(i);
    }
    return projs;
  }


  utility::matrix
  PCA::projection_transposed(const utility::matrix& samples) const
  {
    const size_t Ncol = samples.columns();
    const size_t Nrow = samples.rows();
    utility::matrix projs( Nrow, Ncol );

    utility::vector temp(samples.rows());
    for( size_t j = 0; j < Ncol; ++j ) {
      for (size_t i=0; i<Ncol; ++i )
        temp(i) = samples(i,j); 
      utility::vector centered( Nrow );
      for( size_t i = 0; i < Nrow; ++i ) 
        centered(i)=temp(i)-meanvalues_(i);
      utility::vector proj( eigenvectors_ * centered );
      for( size_t i = 0; i < Nrow; ++i ) 
        projs(i,j)=proj(i);
    }
    return projs;
  }


  void PCA::calculate_explained_intensity(void)
  {
    size_t DIM = eigenvalues_.size();
    explained_intensity_ = utility::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 ];
  }


}}} // of namespace utility, yat, and theplu