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

Last change on this file since 687 was 687, checked in by Jari Häkkinen, 16 years ago

Clean up of namespace comments.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.7 KB
Line 
1// $Id: PCA.cc 687 2006-10-16 23:51:10Z jari $
2
3/*
4  Copyright (C) The authors contributing to this file.
5
6  This file is part of the yat library, http://lev.thep.lu.se/trac/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 2 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with this program; if not, write to the Free Software
20  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
21  02111-1307, USA.
22*/
23
24#include "PCA.h"
25#include "SVD.h"
26
27#include <iostream>
28#include <cmath>
29#include <memory>
30
31namespace theplu {
32namespace yat {
33namespace utility {
34
35
36  void PCA::process(void)
37  {
38    process_ = true;
39    // Row-center the data matrix
40    utility::matrix A_center( A_.rows(), A_.columns() );
41    this->row_center( A_center );
42
43    // Single value decompose the data matrix
44    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
45    pSVD->decompose();
46    utility::matrix U = pSVD->U(); 
47    utility::matrix V = pSVD->V();
48
49    // Read the eigenvectors and eigenvalues
50    eigenvectors_ = U;
51    eigenvectors_ .transpose();
52    eigenvalues_ = pSVD->s();
53
54    // T
55    for( size_t i = 0; i < eigenvalues_.size(); ++i )
56      eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
57    eigenvalues_ *= 1.0/A_center.rows();
58
59    // Sort the eigenvectors in order of eigenvalues
60    // Simple (not efficient) algorithm that always
61    // make sure that the i:th element is in its correct
62    // position (N element --> Ordo( N*N ))
63    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
64      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
65        if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
66          std::swap( eigenvalues_[ i ], eigenvalues_[ j ] ); 
67          eigenvectors_.swap_rows(i,j);
68        }
69  }
70
71
72  void PCA::process_transposed_problem(void)
73  {
74    process_ = true;
75    // Row-center the data matrix
76    utility::matrix A_center( A_.rows(), A_.columns() );
77    this->row_center( A_center );
78
79    // Transform into SVD friendly dimension
80    A_.transpose();
81    A_center.transpose();
82
83    // Single value decompose the data matrix
84    std::auto_ptr<SVD> pSVD( new SVD( A_center ) );
85    pSVD->decompose();
86    utility::matrix U = pSVD->U(); 
87    utility::matrix V = pSVD->V();
88
89    // Read the eigenvectors and eigenvalues
90    eigenvectors_=V;
91    eigenvectors_.transpose();
92    eigenvalues_ = pSVD->s();
93
94    // Transform back when done with SVD!
95    // (used V insted of U now for eigenvectors)
96    A_.transpose();
97    A_center.transpose();
98
99    // T
100    for( size_t i = 0; i < eigenvalues_.size(); ++i )
101      eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
102    eigenvalues_ *= 1.0/A_center.rows();
103
104    // Sort the eigenvectors in order of eigenvalues
105    // Simple (not efficient) algorithm that always
106    // make sure that the i:th element is in its correct
107    // position (N element --> Ordo( N*N ))
108    for ( size_t i = 0; i < eigenvalues_.size(); ++i )
109      for ( size_t j = i + 1; j < eigenvalues_.size(); ++j )
110        if ( eigenvalues_[ j ] > eigenvalues_[ i ] ) {
111          std::swap( eigenvalues_[ i ], eigenvalues_[ j ] ); 
112          eigenvectors_.swap_rows(i,j);
113        }
114  }
115
116
117  // This function will row-center the matrix A_,
118  // that is, A_ = A_ - M, where M is a matrix
119  // with the meanvalues of each row
120  void PCA::row_center(utility::matrix& A_center)
121  {
122    meanvalues_ = utility::vector( A_.rows() );
123    utility::vector A_row_sum(A_.rows());
124    for (size_t i=0; i<A_row_sum.size(); ++i)
125      A_row_sum(i)=utility::vector(A_,i).sum();
126    for( size_t i = 0; i < A_center.rows(); ++i ) {
127      meanvalues_[i] = A_row_sum(i) / static_cast<double>(A_.columns());
128      for( size_t j = 0; j < A_center.columns(); ++j )
129        A_center(i,j) = A_(i,j) - meanvalues_(i);
130    }
131  }
132
133
134  utility::matrix PCA::projection(const utility::matrix& samples ) const
135  {
136    const size_t Ncol = samples.columns();
137    const size_t Nrow = samples.rows();
138    utility::matrix projs( Ncol, Ncol );
139
140    utility::vector temp(samples.rows());
141    for( size_t j = 0; j < Ncol; ++j ) {
142      for (size_t i=0; i<Ncol; ++i )
143        temp(i) = samples(i,j); 
144      utility::vector centered( Nrow );
145      for( size_t i = 0; i < Nrow; ++i ) 
146        centered(i) = temp(i) - meanvalues_(i);
147      utility::vector proj( eigenvectors_ * centered );
148      for( size_t i = 0; i < Ncol; ++i ) 
149        projs(i,j)=proj(i);
150    }
151    return projs;
152  }
153
154
155  utility::matrix
156  PCA::projection_transposed(const utility::matrix& samples) const
157  {
158    const size_t Ncol = samples.columns();
159    const size_t Nrow = samples.rows();
160    utility::matrix projs( Nrow, Ncol );
161
162    utility::vector temp(samples.rows());
163    for( size_t j = 0; j < Ncol; ++j ) {
164      for (size_t i=0; i<Ncol; ++i )
165        temp(i) = samples(i,j); 
166      utility::vector centered( Nrow );
167      for( size_t i = 0; i < Nrow; ++i ) 
168        centered(i)=temp(i)-meanvalues_(i);
169      utility::vector proj( eigenvectors_ * centered );
170      for( size_t i = 0; i < Nrow; ++i ) 
171        projs(i,j)=proj(i);
172    }
173    return projs;
174  }
175
176
177  void PCA::calculate_explained_intensity(void)
178  {
179    size_t DIM = eigenvalues_.size();
180    explained_intensity_ = utility::vector( DIM );
181    double sum = 0;
182    for( size_t i = 0; i < DIM; ++i ) 
183      sum += eigenvalues_[ i ];
184    double exp_sum = 0;
185    for( size_t i = 0; i < DIM; ++i ) {
186      exp_sum += eigenvalues_[ i ];
187      explained_intensity_[ i ] = exp_sum / sum ;
188    }
189  }
190
191
192  double PCA::get_explained_intensity( const size_t& k )
193  {
194    if( !explained_calc_ ) {
195      this->calculate_explained_intensity();
196      explained_calc_ = true;
197    }
198    return explained_intensity_[ k ];
199  }
200
201
202}}} // of namespace utility, yat, and theplu
Note: See TracBrowser for help on using the repository browser.