# source:trunk/src/PCA.cc@13

Last change on this file since 13 was 13, checked in by daniel, 20 years ago

* empty log message *

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Author Date Id Revision`
File size: 5.6 KB
Line
1#include "PCA.h"
2
3
4using namespace std;
5using namespace thep_cpp_tools;
6
7PCA::PCA() : process_( false ), explained_calc_( false )
8{
9}
10
11PCA::PCA( const thep_gsl_api::matrix& A ) : A_( A ),
12              process_( false ),
13              explained_calc_( false )
14{
15}
16
17void PCA::process()
18{
19  process_ = true;
20  // Row-center the data matrix
21  thep_gsl_api::matrix A_center( A_.rows(), A_.cols() );
22  this->row_center( A_center );
23
24
25  // Single value decompose the data matrix
26  auto_ptr<thep_cpp_tools::SVD> pSVD( new thep_cpp_tools::SVD( A_center ) );
27  pSVD->process();
28  thep_gsl_api::matrix U = pSVD->get_u();
29  thep_gsl_api::matrix V = pSVD->get_v();
30
31  // Read the eigenvectors and eigenvalues
32  eigenvectors_ = U.transpose();
33  eigenvalues_ = pSVD->get_s_vec();
34
35  // T
36  for( size_t i = 0; i < eigenvalues_.size(); ++i )
37    {
38     eigenvalues_[ i ] = eigenvalues_[ i ]*eigenvalues_[ i ];
39    }
40
41  eigenvalues_ /= A_center.rows();
42
43  // Sort the eigenvectors in order of eigenvalues
44  // Simple (not efficient) algorithm that always
45  // make sure that the i:th element is in its correct
46  // position (N element --> Ordo( N*N ))
47  for( size_t i = 0; i < eigenvalues_.size(); ++i )
48    {
49     for( size_t j = i + 1; j < eigenvalues_.size(); ++j )
50       {
51  if( eigenvalues_[ j ] > eigenvalues_[ i ] )
52    {
53     swap( eigenvalues_[ i ], eigenvalues_[ j ] );
54     eigenvectors_.swap_rows( i, j );
55    }
56       }
57    }
58}
59
60
61// This function will row-center the matrix A_,
62// that is, A_ = A_ - M, where M is a matrix
63// with the meanvalues of each row
64void PCA::row_center( thep_gsl_api::matrix& A_center )
65{
66  meanvalues_ = thep_gsl_api::vector( A_.rows() );
67  thep_gsl_api::vector A_row_sum = A_.row_sum();
68  for( size_t i = 0; i < A_center.rows(); ++i )
69    {
70     meanvalues_[ i ] = A_row_sum.get( i ) / static_cast<double>( A_.cols() );
71     for( size_t j = 0; j < A_center.cols(); ++j )
72       {
73  A_center.set(  i, j, A_.get( i , j ) - meanvalues_[ i ] );
74       }
75    }
76}
77
78
79thep_gsl_api::matrix PCA::projection( const thep_gsl_api::matrix&
80              samples ) const
81{
82  assert( process_ );
83  const size_t Ncol = samples.cols();
84  const size_t Nrow = samples.rows();
85  thep_gsl_api::matrix projs( Ncol, Ncol );
86
87  for( size_t j = 0; j < Ncol; ++j )
88    {
89     thep_gsl_api::matrix temp = samples.col( j );
90     thep_gsl_api::vector centered( Nrow );
91
92     for( size_t i = 0; i < Nrow; ++i )
93       {
94  centered[ i ] = temp.get( i, 0 ) - meanvalues_[ i ];
95       }
96
97     thep_gsl_api::vector proj( eigenvectors_ * centered );
98      for( size_t i = 0; i < samples.cols(); ++i )
99       {
100  projs.set( i, j,  proj[ i ] );
101       }
102    }
103  return projs;
104}
105
106
107void PCA::calculate_explained_intensity()
108{
109  size_t DIM = eigenvalues_.size();
110  explained_intensity_ = thep_gsl_api::vector( DIM );
111  double sum = 0;
112  for( size_t i = 0; i < DIM; ++i )
113    {
114     sum += eigenvalues_[ i ];
115    }
116  double exp_sum = 0;
117  for( size_t i = 0; i < DIM; ++i )
118    {
119     exp_sum += eigenvalues_[ i ];
120     explained_intensity_[ i ] = exp_sum / sum ;
121    }
122}
123
124double PCA::get_explained_intensity( const size_t& k )
125{
126  if( !explained_calc_ )
127    {
128     this->calculate_explained_intensity();
129     explained_calc_ = true;
130    }
131  return explained_intensity_[ k ];
132}
133
134
135bool PCA::test()
136{
137  cout << "1. Print original matrix A = \n";//Testing row-centering of matrix\n";
138  thep_gsl_api::matrix A( 3, 3 );
139  for( size_t i = 0; i < 3; ++i )
140    {
141     for( size_t j = 0; j < 3; ++j )
142       {
143  A.set( i, j, sin( static_cast<double>( i + j + 2 + (i+1)*(j+1) ) ) );
144       }
145    }
146  A_ = A;
147
148// Print only matrix that is row-centered as PCA ( process() )
149// will do that before calculation
150  this->row_center( A );
151  cout << A << endl;
152
153  cout << "2. Testing PCA\n";
154  this->process();
155
156  for( size_t i = 0; i < eigenvalues_.size(); ++i )
157    {
158     cout << "lambda(" << i << ") = " << eigenvalues_[i] << " "
159    << "has eigenvector = " << endl;
160     cout << eigenvectors_.row( i ) << endl;
161    }
162
163  // The eigenvalues and eigenvectors are
164  // from matlab using eig(A'A)/3
165  double answer_lambda[ 3 ] = { 1.8908, 0.0327, 0 };
166  double answer_eigenvec[ 3 ][ 3 ] =
167    {
168      { 0.4763, -0.6738, 0.5649  },
169      { 0.8770, 0.3182,  -0.3599 },
170      { 0.0627, 0.6669 , 0.7425  },
171    };
172
173  double ERROR_TOL = 0.0001;
174
175  // error eigenvalues
176  cout << "abs( lambda(i) - lambda_matlab(i) ) ) = \n( ";
177  double max_err = 0;
178  for( size_t i = 0; i < 3; ++i )
179    {
180     double err = fabs( eigenvalues_[ i ] - answer_lambda[ i ] );
181     if( err > max_err ) max_err = err;
182     cout << err << " ";
183    }
184  cout << ")\n";
185
186  cerr << "max_error = " << max_err << ", to be compared with ERROR_TOL = "
187       << ERROR_TOL << endl;
188  if( max_err < ERROR_TOL )
189    {
190     cout << "ok, compared with precision in given matlab-references!\n";
191    }
192  else
193    {
194     cerr << "Not ok, error to big!\n";
195     return false;
196    }
197
198  // error eigenvectors
199  max_err = 0;
200  for( size_t i = 0; i < 3; ++i )
201    {
202     cout << "err_eigvec( " << i << " ) = \n";
203     for( size_t j = 0; j < 3; ++j )
204       {
205  double err = fabs( fabs( eigenvectors_.get( i, j ) ) - fabs( answer_eigenvec[ i ][ j ] ) );
206  if( err > max_err ) max_err = err;
207  cout << err << endl;
208       }
209     cout << ")\n\n";
210    }
211
212  cerr << "max_error = " << max_err << ", to be compared with ERROR_TOL = "
213       << ERROR_TOL << endl;
214  if( max_err < ERROR_TOL )
215    {
216     cout << "ok, compared with precision in given matlab-references!\n";
217    }
218  else
219    {
220     cerr << "Not ok, error to big!\n";
221     return false;
222    }
223
224  cout << "\nThe following projection is obtained:\n";
225  cout << this->projection( A ) << endl;
226
227
228  return true;
229}
Note: See TracBrowser for help on using the repository browser.