pca.hpp

Go to the documentation of this file.

/*
 * Copyright(c):
 * Signal Image and Communications (SIC) Department 
 * http://www.sic.sp2mi.univ-poitiers.fr/ 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - XLIM  Institute UMR CNRS 7252  http://www.xlim.fr/ 
 * 
 * and 
 * 
 * D2 Fluid, Thermic and Combustion 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - PPRIME Institute -  UPR CNRS 3346  http://www.pprime.fr
 *  - ISAE-ENSMA http://www.ensma.fr
 *
 * Contributor(s):
 * The SLIP team,
 * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
 * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>, 
 * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>, 
 * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>, 
 * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>, 
 * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
 *
 * Description:
 * The Simple Library of Image Processing (SLIP) is a new image processing 
 * library. It is written in the C++ language following as much as possible 
 * the ISO/ANSI C++ standard. It is consequently compatible with any system  
 * satisfying the ANSI C++ complience. It works on different Unix , Linux ,  
 * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that  
 * was created by the Signal, Image and Communications (SIC) departement of  
 * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic  
 * and Combustion departement of the P', UPR 3346 CNRS Institute of the  
 * University of Poitiers.
 * 
 * The SLIP Library source code has been registered to the APP (French Agency 
 * for the Protection of Programs) by the University of Poitiers and CNRS, 
 * under  registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.

 * http://www.sic.sp2mi.univ-poitiers.fr/slip/
 *
 * This software is governed by the CeCILL-C license under French law and
 * abiding by the rules of distribution of free software.  You can  use, 
 * modify and/ or redistribute the software under the terms of the CeCILL-C 
 * license as circulated by CEA, CNRS and INRIA at the following URL 
 * http://www.cecill.info.
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 * 
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the 
 * same conditions as regards security. 
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-C license and that you accept its terms.
 */


#ifndef SLIP_PCA_HPP
#define SLIP_PCA_HPP

#include <cmath>
#include <map>
#include <complex>
#include <numeric>
#include <algorithm>
#include "macros.hpp"
#include "Matrix.hpp"
#include "Vector.hpp"
#include "linear_algebra.hpp"
#include "linear_algebra_eigen.hpp"
#include "statistics.hpp"

namespace slip
{


  template<class OutputMatrix> 
  inline
  void pca_uniform_metric(OutputMatrix& metric)
  {
    assert (metric.rows() >= metric.cols());
    assert (metric.rows() == metric.cols());
    slip::identity(metric);
  }

  template<class Matrix, class OutputMatrix> 
  inline
  void pca_adaptative_metric(const Matrix & data, 
                                                     OutputMatrix& metric)
  {
    assert (data.rows() >= data.cols());
    assert (metric.rows() == data.cols());
    assert (metric.rows() == metric.cols());
    typedef typename Matrix::value_type value_type; 
    std::size_t nbvector = data.rows();
    std::size_t size = data.cols();
    std::size_t ncmetric = metric.cols();
    value_type meanval = value_type(0.0);
    value_type stddev  = value_type(0.0);
    for (std::size_t j = 0; j < size; ++j)
      {
                meanval = value_type(0.0);
                for (std::size_t i = 0; i < nbvector; ++i)
                  {
                    meanval += data[i][j];
                  }
                stddev = value_type(0.0);
                for (std::size_t i = 0; i < nbvector; ++i)
                  {
                    stddev += (data[i][j] - meanval) * (data[i][j] - meanval);
                  }
                
                for (std::size_t i = 0; i < ncmetric; ++i)
                  {
                    metric[i][j] = value_type(0.0);
                  }
                metric[j][j] = value_type(1.0) / std::sqrt(stddev);
      }
  }

  template<class Matrix, class OutputMatrix> 
  inline
  void pca_center_data(const Matrix & data, 
                                       OutputMatrix& centerdata, 
                                       const bool normalize_data = false)
  {
                assert (data.rows() >= data.cols());
                assert (data.rows() == centerdata.rows());
                assert (data.cols() == centerdata.cols());
                typedef typename Matrix::value_type  value_type;
                std::size_t nbvector = data.rows();
                std::size_t size = data.cols();
                // Copy datas into matrix structure
                std::copy(data.begin(),data.end(),centerdata.begin());

                // compute means and center datas
                value_type meanval = value_type(0.0);
                value_type stddev = value_type(0.0);
                for (std::size_t j = 0; j < size; ++j)
                  {
                    meanval = value_type(0.0);
                    for (std::size_t i = 0; i < nbvector; ++i)
                      {
                                meanval += data[i][j];
                      }
                    for (std::size_t i = 0; i < nbvector; ++i)
                      {
                                centerdata[i][j] -= meanval;
                      }
                    if (normalize_data)
                      {

                                stddev = value_type(0.0);
                                for (std::size_t i = 0; i < nbvector; ++i)
                                  {
                                    stddev += (data[i][j] - meanval) * (data[i][j] - meanval);
                                  }
                                
                                for (std::size_t i = 0; i < nbvector; ++i)
                                  {
                                    centerdata[i][j] /= stddev;
                                  }
                      }
                  }
  }

  template<class Matrix, class Matrix2, class OutputMatrix,class OutputVector> 
  inline
  void pca_computation(const Matrix & data, 
                                       OutputMatrix &acp, 
                                       const Matrix2& metric, 
                                       OutputVector & eigenval)
  {
                assert( acp.rows() == data.rows());
                assert( acp.cols() == data.cols());
                assert( data.rows() >= data.cols());
                assert( eigenval.size() == data.cols());
                assert( metric.rows() == data.cols());
                assert( metric.rows() == metric.cols());
                typedef typename Matrix::value_type value_type;
                std::size_t nbvector = data.rows();
                std::size_t size = data.cols();
                // Compute variance / covariance matric
        slip::Matrix<value_type> varcovar(size,size);
        slip::Matrix<value_type> T(nbvector,size);
                slip::Matrix<value_type> Ttranspose(size,nbvector);
                slip::Matrix<value_type> temp(size,nbvector);
        slip::Matrix<value_type> P(nbvector,nbvector,value_type(0.0));
                for (std::size_t i = 0; i < nbvector; ++i)
                  {
                    P[i][i] = value_type(1.0) / ((value_type) nbvector);
                  }

                slip::matrix_matrix_multiplies(data,metric,T);
                slip::transpose(T,Ttranspose);
                slip::matrix_matrix_multiplies(Ttranspose,P,temp);
                slip::matrix_matrix_multiplies(temp,T,varcovar);

                //computes the eigen values and eigen vectors of varcovar
                slip::Matrix<value_type> eigenvectors(size,size);
                slip::hermitian_eigen(varcovar,eigenval,eigenvectors);

                // Compute ACP
                slip::matrix_matrix_multiplies(T,eigenvectors,acp);
  }


  template<class Matrix, class OutputMatrix> 
  inline
  void pca_simple_data_reduction_axis(const Matrix & acp, 
                                                                      OutputMatrix &result, 
                                                                      const std::size_t & nbaxis)
  {
                assert(acp.rows() == result.rows());
                assert(nbaxis <= acp.cols());
                assert(acp.rows() >= acp.cols());
                std::size_t nbvector = acp.rows();
                for (std::size_t i = 0; i < nbvector; ++i)
                  {
                    for (std::size_t j = 0; j < nbaxis; ++j)
                      {
                                result[i][j] = acp[i][j];
                      }
                  }
  }

  template<class Matrix, class Matrix2, class OutputMatrix> 
  inline
  void pca_data_reduction(const Matrix & data, 
                                                  const Matrix2& metric, 
                                                  OutputMatrix &result, 
                                                  const std::size_t & nbaxis,
                                                  const bool normalize_data = false)
  {
    typedef typename Matrix::value_type real;
    std::size_t nbvector = data.rows();
    std::size_t size = data.cols();
    
    slip::Matrix<real> acp(nbvector,size);
    slip::Matrix<real> centerdata(nbvector,size);
    slip::Vector<real> eigenval(size);

    slip::pca_center_data(data,centerdata,normalize_data);
                
    slip::pca_computation(data, 
                                                  acp, 
                                                  metric,
                                                  eigenval);
    
    slip::pca_simple_data_reduction_axis(acp, 
                                                                                 result,nbaxis);

  }

  template<class Matrix, class Matrix2, class OutputMatrix> 
  inline
  void pca_data_reduction(const Matrix & data, 
                                                  const Matrix2& metric, 
                                                  OutputMatrix &result, 
                                                  const bool normalize_data = false)
  {
    slip::pca_data_reduction(data, 
                                                     metric, 
                                                     result, result.cols(),normalize_data);
  }




}


#endif //SLIP_PCA_HPP