source: trunk/yat/utility/PCA.cc

Last change on this file was 3938, checked in by Peter, 15 months ago

let configure fail if compiler is not a C++11 compiler. YAT_HAVE_RVALUE and friends are no longer defined in 'config.h' and guards araound declarations and implementations are removed. refs #949

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.2 KB
Line 
1// $Id: PCA.cc 3938 2020-07-16 13:16:56Z peter $
2
3/*
4  Copyright (C) 2003 Daniel Dalevi
5  Copyright (C) 2004 Jari Häkkinen
6  Copyright (C) 2005 Jari Häkkinen, Peter Johansson
7  Copyright (C) 2006 Jari Häkkinen
8  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
9  Copyright (C) 2010, 2012, 2018, 2019 Peter Johansson
10
11  This file is part of the yat library, http://dev.thep.lu.se/yat
12
13  The yat library is free software; you can redistribute it and/or
14  modify it under the terms of the GNU General Public License as
15  published by the Free Software Foundation; either version 3 of the
16  License, or (at your option) any later version.
17
18  The yat library is distributed in the hope that it will be useful,
19  but WITHOUT ANY WARRANTY; without even the implied warranty of
20  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21  General Public License for more details.
22
23  You should have received a copy of the GNU General Public License
24  along with yat. If not, see <http://www.gnu.org/licenses/>.
25*/
26
27#include <config.h>
28
29#include "PCA.h"
30#include "SVD.h"
31#include "utility.h"
32#include "Vector.h"
33
34#include <boost/smart_ptr/scoped_ptr.hpp>
35
36#include <iostream>
37#include <cmath>
38#include <memory>
39#include <utility>
40
41namespace theplu {
42namespace yat {
43namespace utility {
44
45
46  PCA::PCA(const Matrix& A)
47    : A_(A)
48  {
49    process();
50  }
51
52
53  PCA::PCA(Matrix&& A)
54    : A_(std::move(A))
55  {
56    process();
57  }
58
59
60  const Vector& PCA::eigenvalues(void) const
61  {
62    return eigenvalues_;
63  }
64
65
66  const Matrix& PCA::eigenvectors(void) const
67  {
68    return eigenvectors_;
69  }
70
71
72  void PCA::process(void)
73  {
74    // Row-center the data matrix
75    Matrix A_center( A_.rows(), A_.columns() );
76    this->row_center( A_center );
77
78    // Single value decompose the data matrix
79    boost::scoped_ptr<SVD> pSVD( new SVD( A_center ) );
80    pSVD->decompose();
81    Matrix U(pSVD->U());
82    Matrix V(pSVD->V());
83
84    // Read the eigenvectors and eigenvalues
85    eigenvectors_=U;
86
87    eigenvectors_ .transpose();
88    eigenvalues_ = pSVD->s();
89
90    // T
91    for( size_t i = 0; i < eigenvalues_.size(); ++i )
92      eigenvalues_(i) = eigenvalues_(i)*eigenvalues_(i);
93    eigenvalues_ *= 1.0/A_center.rows();
94  }
95
96
97  Matrix PCA::projection(const Matrix& samples ) const
98  {
99    const size_t Ncol = samples.columns();
100    const size_t Nrow = samples.rows();
101    Matrix projs( Ncol, Ncol );
102
103    Vector temp(samples.rows());
104    for( size_t j = 0; j < Ncol; ++j ) {
105      for (size_t i=0; i<Ncol; ++i )
106        temp(i) = samples(i,j);
107      Vector centered( Nrow );
108      for( size_t i = 0; i < Nrow; ++i )
109        centered(i) = temp(i) - meanvalues_(i);
110      Vector proj( eigenvectors_ * centered );
111      for( size_t i = 0; i < Ncol; ++i )
112        projs(i,j)=proj(i);
113    }
114    return projs;
115  }
116
117
118  // This function will row-center the matrix A_,
119  // that is, A_ = A_ - M, where M is a matrix
120  // with the meanvalues of each row
121  void PCA::row_center(Matrix& A_center)
122  {
123    meanvalues_ = Vector(A_.rows());
124    Vector A_row_sum(A_.rows());
125    for (size_t i=0; i<A_row_sum.size(); ++i){
126      A_row_sum(i) = sum(A_.row_const_view(i));
127    }
128    for( size_t i = 0; i < A_center.rows(); ++i ) {
129      meanvalues_(i) = A_row_sum(i) / static_cast<double>(A_.columns());
130      for( size_t j = 0; j < A_center.columns(); ++j )
131        A_center(i,j) = A_(i,j) - meanvalues_(i);
132    }
133  }
134
135}}} // of namespace utility, yat, and theplu
Note: See TracBrowser for help on using the repository browser.