source: trunk/src/PCA.cc @ 14

Last change on this file since 14 was 14, checked in by daniel, 19 years ago

One correction in the projection-method.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.8 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  //cout << temp.get( i, 0 ) << " " << meanvalues_[ i ] << " " << centered[ i ] << endl;
96       }
97
98     thep_gsl_api::vector proj( eigenvectors_ * centered );
99      for( size_t i = 0; i < Ncol; ++i ) 
100       {
101  projs.set( i, j,  proj[ i ] );
102       }
103    } 
104  return projs;
105}
106
107
108void PCA::calculate_explained_intensity() 
109{
110  size_t DIM = eigenvalues_.size();
111  explained_intensity_ = thep_gsl_api::vector( DIM );
112  double sum = 0;
113  for( size_t i = 0; i < DIM; ++i ) 
114    {
115     sum += eigenvalues_[ i ];
116    }
117  double exp_sum = 0;
118  for( size_t i = 0; i < DIM; ++i ) 
119    {
120     exp_sum += eigenvalues_[ i ];
121     explained_intensity_[ i ] = exp_sum / sum ;
122    } 
123}
124
125double PCA::get_explained_intensity( const size_t& k )
126{
127  if( !explained_calc_ )
128    {
129     this->calculate_explained_intensity();
130     explained_calc_ = true;
131    }
132  return explained_intensity_[ k ];
133}
134
135
136bool PCA::test()
137{
138  cout << "1. Print original matrix A = \n";//Testing row-centering of matrix\n";
139  thep_gsl_api::matrix A( 4, 3 );
140  for( size_t i = 0; i < 4; ++i )
141    {
142     for( size_t j = 0; j < 3; ++j )
143       {
144  A.set( i, j, sin( static_cast<double>( i + j + 2 + (i+1)*(j+1) ) ) );
145       }
146    }
147  A_ = A;
148
149// Print only matrix that is row-centered as PCA ( process() )
150// will do that before calculation
151  this->row_center( A );
152  cout << A << endl;
153
154  cout << "2. Testing PCA\n";
155  this->process();
156
157  for( size_t i = 0; i < eigenvalues_.size(); ++i )
158    {
159     cout << "lambda(" << i << ") = " << eigenvalues_[i] << " " 
160    << "has eigenvector = " << endl;
161     cout << eigenvectors_.row( i ) << endl;
162    }
163 
164  // The eigenvalues and eigenvectors are
165  // from matlab using eig(A'A)/3
166//   double answer_lambda[ 3 ] = { 1.8908, 0.0327, 0 };
167//   double answer_eigenvec[ 3 ][ 3 ] =
168//     {
169//       { 0.4763, -0.6738, 0.5649  },
170//       { 0.8770, 0.3182,  -0.3599 },
171//       { 0.0627, 0.6669 , 0.7425  },
172//     };
173
174//   double ERROR_TOL = 0.0001;
175
176//   // error eigenvalues
177//   cout << "abs( lambda(i) - lambda_matlab(i) ) ) = \n( ";
178//   double max_err = 0;
179//   for( size_t i = 0; i < 3; ++i )
180//     {
181//      double err = fabs( eigenvalues_[ i ] - answer_lambda[ i ] );
182//      if( err > max_err ) max_err = err;
183//      cout << err << " ";
184//     }
185//   cout << ")\n";
186
187//   cerr << "max_error = " << max_err << ", to be compared with ERROR_TOL = "
188//        << ERROR_TOL << endl;
189//   if( max_err < ERROR_TOL )
190//     {
191//      cout << "ok, compared with precision in given matlab-references!\n";
192//     }
193//   else
194//     {
195//      cerr << "Not ok, error to big!\n";
196//      return false;
197//     }
198
199//   // error eigenvectors
200//   max_err = 0;
201//   for( size_t i = 0; i < 3; ++i )
202//     {
203//      cout << "err_eigvec( " << i << " ) = \n";
204//      for( size_t j = 0; j < 3; ++j )
205//        {
206//  double err = fabs( fabs( eigenvectors_.get( i, j ) ) - fabs( answer_eigenvec[ i ][ j ] ) );
207//  if( err > max_err ) max_err = err;       
208//  cout << err << endl;
209//        }
210//      cout << ")\n\n";
211//     }
212
213//   cerr << "max_error = " << max_err << ", to be compared with ERROR_TOL = "
214//        << ERROR_TOL << endl;
215//   if( max_err < ERROR_TOL )
216//     {
217//      cout << "ok, compared with precision in given matlab-references!\n";
218//     }
219//   else
220//     {
221//      cerr << "Not ok, error to big!\n";
222//      return false;
223//     }
224   
225  cout << "\nThe following projection is obtained:\n";
226  cout << this->projection( A_ ) << endl;
227 
228
229  return true;
230}
Note: See TracBrowser for help on using the repository browser.