source: trunk/yat/utility/PCA.h @ 703

#ifndef _theplu_yat_utility_pca_ 
#define _theplu_yat_utility_pca_ 

// $Id: PCA.h 703 2006-12-18 00:47:44Z jari $

/*
  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 "matrix.h"
#include "vector.h"

namespace theplu {
namespace yat {
namespace utility {

  /**
     Class performing PCA using SVD. This class assumes that 
     the columns corresponds to the dimenension of the problem.
     That means if data has dimension NxM (M=columns) the number
     of principal-axes will equal M-1. When projecting data into 
     this space, all Nx1 vectors will have dimension Mx1. Hence
     the projection will have dimension MxM where each column is 
     a point in the new space. Also, it assumes that M>N. The opposite
     problem is added in the functions: process_transposed_problem and
     projection_transposed()...
  */
  class PCA
  {
  public:
    /**
       Constructor taking the data-matrix as input. No row-centering
       should have been performed and no products. 
     */
    explicit PCA(const utility::matrix&);
 
    /**
       Will perform PCA according to the following scheme: \n
       1: Rowcenter A  \n
       2: SVD(A)  --> USV' \n
       3: Calculate eigenvalues according to \n 
          \f$ \lambda_{ii} = s_{ii}/N_{rows} \f$ \n
       4: Sort eigenvectors (from matrix V) according to descending eigenvalues\n
    */
    void process(void);

    /**
       If M<N use this method instead. Using the same format as before 
       where rows in the matrix corresponds to the dimensional coordinate. 
       The only difference is in the SVD step where the matrix V is used 
       after running the transposed matrix. For projections, see 
       projection_transposed() method.
     */
    void process_transposed_problem(void);

    /**
       @return Eigenvector \a i.
    */
    inline utility::vector get_eigenvector(size_t i) const
      { return utility::vector(eigenvectors_,i); }

    /**
       Returns eigenvalues to covariance matrix 
       \f$ C = \frac{1}{N^2}A^TA \f$
    */
    inline double get_eigenvalue(size_t i) const { return eigenvalues_[i]; }

    /**
       Returns the explained intensity of component \a K \n
       \f$I = \frac{ \sum^{K}_{i=1} \lambda_i }{ \sum^{N}_{j=1} \lambda_j }\f$ \n
       where \f$N\f$ is the dimension 
    */
    double get_explained_intensity( const size_t& k );

    /**
       This function will project data onto the new coordinate-system
       where the axes are the calculated eigenvectors. This means that
       PCA must have been run before this function can be used!
       Output is presented as coordinates in the N-dimensional room
       spanned by the eigenvectors.
    */
    utility::matrix projection( const utility::matrix& ) const;

    /**
       Same as projection() but works when used
       process_transposed_problem().
    */
    utility::matrix projection_transposed( const utility::matrix& ) const;


  private:
    utility::matrix A_; 
    utility::matrix eigenvectors_;
    utility::vector eigenvalues_;
    utility::vector explained_intensity_;
    utility::vector meanvalues_;
    bool process_, explained_calc_;
    
    /**
       Private function that will row-center the matrix A,
       that is, A = A - M, where M is a matrix
       with the meanvalues of each row
    */
    void row_center( utility::matrix& A_center );

    /**
       Private function that will calculate the explained 
       intensity
    */
    void calculate_explained_intensity();
  }; // class PCA  


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

#endif