Marsyas
0.6.0-alpha
|
00001 /* 00002 ** Copyright (C) 1998-2006 George Tzanetakis <gtzan@cs.uvic.ca> 00003 ** 00004 ** This program is free software; you can redistribute it and/or modify 00005 ** it under the terms of the GNU General Public License as published by 00006 ** the Free Software Foundation; either version 2 of the License, or 00007 ** (at your option) any later version. 00008 ** 00009 ** This program is distributed in the hope that it will be useful, 00010 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00012 ** GNU General Public License for more details. 00013 ** 00014 ** You should have received a copy of the GNU General Public License 00015 ** along with this program; if not, write to the Free Software 00016 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 00017 */ 00018 00019 #include "PCA.h" 00020 #include "../common_source.h" 00021 00022 #include "config.h" 00023 00024 using std::ostringstream; 00025 using namespace Marsyas; 00026 00027 #define SIGN(a, b) ( (b) < 0 ? -fabs(a) : fabs(a) ) 00028 00029 00030 PCA::PCA(mrs_string name):MarSystem("PCA",name) 00031 { 00032 //type_ = "PCA"; 00033 //name_ = name; 00034 00035 addControls(); 00036 } 00037 00038 00039 PCA::PCA(const PCA& a): MarSystem(a) 00040 { 00041 evals_ = NULL; 00042 interm_ = NULL; 00043 } 00044 00045 00046 PCA::~PCA() 00047 { 00048 delete [] evals_; 00049 delete [] interm_; 00050 } 00051 00052 MarSystem* 00053 PCA::clone() const 00054 { 00055 return new PCA(*this); 00056 } 00057 00058 void 00059 PCA::addControls() 00060 { 00061 npcs_.create(3,3); 00062 00063 addctrl("mrs_natural/npc",4); 00064 setctrlState("mrs_natural/npc",true); 00065 addctrl("mrs_realvec/pcs",npcs_); 00066 dims_ = 0; 00067 evals_ = NULL; 00068 interm_ = NULL; 00069 } 00070 00071 void 00072 PCA::myUpdate(MarControlPtr sender) 00073 { 00074 (void) sender; //suppress warning of unused parameter(s) 00075 MRSDIAG("PCA.cpp - PCA:myUpdate"); 00076 00077 setctrl("mrs_natural/onSamples", getctrl("mrs_natural/inSamples")); 00078 setctrl("mrs_natural/onObservations", getctrl("mrs_natural/npc")); 00079 setctrl("mrs_real/osrate", getctrl("mrs_real/israte")); 00080 00081 inObservations_ = getctrl("mrs_natural/inObservations")->to<mrs_natural>(); 00082 onObservations_ = getctrl("mrs_natural/onObservations")->to<mrs_natural>(); 00083 00084 npc_ = getctrl("mrs_natural/npc")->to<mrs_natural>(); 00085 00086 if( npcs_.getRows() != inObservations_-1 || npcs_.getCols() != npc_ ) 00087 npcs_.create(inObservations_-1,npc_); 00088 00089 if( npc_ != onObservations_-1 ) { 00090 00091 updControl("mrs_natural/onObservations", npc_ +1); 00092 onObservations_ = npc_+1; 00093 } 00094 00095 if( dims_ != inObservations_-1 ) 00096 { 00097 dims_ = inObservations_-1; 00098 corr_matrix_.create(dims_,dims_); 00099 temp_matrix_.create(dims_, inSamples_); 00100 evals_ = new mrs_real[dims_]; 00101 interm_ = new mrs_real[dims_]; 00102 } 00103 00104 ostringstream oss; 00105 for (int i=1 ; i <= npc_ ; ++i ) 00106 { 00107 oss << "PC_" << i << ","; 00108 } 00109 setctrl("mrs_string/onObsNames", oss.str()); 00110 } 00111 00112 void 00113 PCA::myProcess(realvec& in, realvec& out) 00114 { 00115 mrs_natural t,o; 00116 mrs_natural o1,o2; 00117 00118 00119 for (o=0; o < inObservations_-1; o++) 00120 for (t = 0; t < inSamples_; t++) 00121 { 00122 temp_matrix_(o, t) = in(o,t); 00123 } 00124 00125 00126 00127 //in.meanSample(means_);//original code 00128 //in.stdSample(means_,stds_); //original code 00129 //means_.create(in.getRows());//not needed anymore if using new realvec::operator=() ;-) 00130 //stds_.create(in.getRows());//not needed anymore if using new realvec::operator=() ;-) 00131 temp_matrix_.meanObs(means_); 00132 temp_matrix_.stdObs(stds_); 00133 00134 00135 00136 00137 00138 00139 // Adjust data : ( X - means(X) ) / ( sqrt(n) * stds(X) ) 00140 for (o=0; o < inObservations_-1; o++) 00141 for (t = 0; t < inSamples_; t++) 00142 temp_matrix_(o,t) = ( temp_matrix_(o,t) - means_(o) ) / ( sqrt((mrs_real)inSamples_) * stds_(o) ) ; 00143 00144 00145 // Calculate the correlation matrix 00146 for ( o1 = 0 ; o1 < inObservations_-1; o1++ ) 00147 { 00148 corr_matrix_(o1,o1) = 1.0; 00149 for ( o2 = o1; o2 < inObservations_-1 ; o2++ ) 00150 { 00151 corr_matrix_(o1,o2) = 0.0; 00152 for( t=0 ; t < inSamples_ ; t++ ) 00153 corr_matrix_(o1,o2) += temp_matrix_(o1,t) * temp_matrix_(o2,t); 00154 corr_matrix_(o2,o1) = corr_matrix_(o1,o2); 00155 } 00156 } 00157 corr_matrix_(inObservations_-2, inObservations_-2) = 1.0; 00158 00159 00160 00161 // Triangular decomposition 00162 tred2(corr_matrix_, inObservations_-1, evals_, interm_); 00163 00164 00165 // Reduction of symmetric tridiagonal matrix 00166 tqli( evals_, interm_, inObservations_-1, corr_matrix_); 00167 00168 00169 00170 mrs_real percent_eig = 0.0; 00171 mrs_real sum_eig = 0.0; 00172 for (int m = inObservations_-2; m >= 0; m--) 00173 sum_eig += evals_[m]; 00174 00175 00176 for (int m = inObservations_-2; m >= 0; m--) 00177 { 00178 percent_eig += evals_[m]; 00179 std::cout << evals_[m] / sum_eig << "\t"; 00180 std::cout << percent_eig / sum_eig << std::endl; 00181 } 00182 00183 /* evals now contains the eigenvalues, 00184 corr_matrix_ now contains the associated eigenvectors. */ 00185 00186 /* Project row data onto the top "npc_" principal components. */ 00187 for( t=0 ; t<inSamples_ ; t++ ) 00188 { 00189 for( o=0 ; o<inObservations_ -1; o++ ) 00190 { 00191 interm_[o] = in(o,t); 00192 } 00193 for( o=0 ; o<onObservations_ -1; o++ ) 00194 { 00195 out(o,t) = 0.0; 00196 for(o2=0 ; o2 < inObservations_ -1; o2++) 00197 { 00198 out(o,t) += interm_[o2] * corr_matrix_(o2,inObservations_-o-2); 00199 npcs_(o2,o) = corr_matrix_(o2,inObservations_-o-2); 00200 } 00201 } 00202 } 00203 00204 // copy the labels 00205 for( t=0 ; t<inSamples_ ; t++ ) 00206 out(onObservations_-1, t) = in(inObservations_-1, t); 00207 setctrl("mrs_realvec/pcs",npcs_); 00208 00209 } 00210 00211 /* Reduce a real, symmetric matrix to a symmetric, tridiag. matrix. */ 00212 void 00213 PCA::tred2(realvec &a, mrs_natural m, mrs_real *d, mrs_real *e) 00214 /* Householder reductiom of matrix a to tridiagomal form. 00215 Algorithm: Martim et al., Num. Math. 11, 181-195, 1968. 00216 Ref: Smith et al., Matrix Eigemsystem Routimes -- EISPACK Guide 00217 Sprimger-Verlag, 1976, pp. 489-494. 00218 W H Press et al., Numerical Recipes im C, Cambridge U P, 00219 1988, pp. 373-374. 00220 00221 Source code adapted from F. Murtagh, Munich, 6 June 1989 00222 http://astro.u-strasbg.fr/~fmurtagh/mda-sw/pca.c 00223 */ 00224 { 00225 mrs_natural l, k, j, i; 00226 mrs_real scale, hh, h, g, f; 00227 00228 for (i = m-1; i > 0; i--) 00229 { 00230 l = i - 1; 00231 h = scale = 0.0; 00232 if (l > 0) 00233 { 00234 for (k = 0; k <= l; k++) 00235 scale += fabs(a(i,k)); 00236 if (scale == 0.0) 00237 e[i] = a(i,l); 00238 else 00239 { 00240 for (k = 0; k <= l; k++) 00241 { 00242 a(i,k) /= scale; 00243 h += a(i,k) * a(i,k); 00244 } 00245 f = a(i,l); 00246 g = f>0 ? -sqrt(h) : sqrt(h); 00247 e[i] = scale * g; 00248 00249 h -= f * g; 00250 a(i,l) = f - g; 00251 f = 0.0; 00252 for (j = 0; j <= l; j++) 00253 { 00254 a(j,i) = a(i,j)/h; 00255 g = 0.0; 00256 for (k = 0; k <= j; k++) 00257 g += a(j,k) * a(i,k); 00258 for (k = j+1; k <= l; k++) 00259 g += a(k,j) * a(i,k); 00260 e[j] = g / h; 00261 f += e[j] * a(i,j); 00262 } 00263 hh = f / (h + h); 00264 for (j = 0; j <= l; j++) 00265 { 00266 f = a(i,j); 00267 e[j] = g = e[j] - hh * f; 00268 for (k = 0; k <= j; k++) 00269 a(j,k) -= (f * e[k] + g * a(i,k)); 00270 } 00271 } 00272 } 00273 else 00274 e[i] = a(i,l); 00275 d[i] = h; 00276 } 00277 d[0] = 0.0; 00278 e[0] = 0.0; 00279 for (i = 0; i < m; ++i) 00280 { 00281 l = i - 1; 00282 if (d[i]) 00283 { 00284 for (j = 0; j <= l; j++) 00285 { 00286 g = 0.0; 00287 for (k = 0; k <= l; k++) 00288 g += a(i,k) * a(k,j); 00289 for (k = 0; k <= l; k++) 00290 a(k,j) -= g * a(k,i); 00291 } 00292 } 00293 d[i] = a(i,i); 00294 a(i,i) = 1.0; 00295 00296 for (j = 0; j <= l; j++) 00297 a(j,i) = a(i,j) = 0.0; 00298 } 00299 } 00300 00301 /* Tridiagonal QL algorithm -- Implicit */ 00302 void 00303 PCA::tqli(mrs_real d[], mrs_real e[], mrs_natural m, realvec &z) 00304 /* 00305 Source code adapted from F. Murtagh, Munich, 6 June 1989 00306 http://astro.u-strasbg.fr/~fmurtagh/mda-sw/pca.c 00307 */ 00308 { 00309 mrs_natural n, l, i, k; 00310 mrs_real s, r, p, g, f, dd, c, b; 00311 00312 for (i = 1; i < m; ++i) 00313 e[i-1] = e[i]; 00314 e[m-1] = 0.0; 00315 for (l = 0; l < m; l++) 00316 { 00317 #ifdef MARSYAS_ASSERTS 00318 mrs_natural iter = 0; 00319 #endif 00320 do 00321 { 00322 for (n = l; n < m-1; n++) 00323 { 00324 dd = fabs(d[n]) + fabs(d[n+1]); 00325 if (fabs(e[n]) + dd == dd) break; 00326 } 00327 if (n != l) 00328 { 00329 MRSASSERT( iter++ != 30 ); // No convergence 00330 00331 g = (d[l+1] - d[l]) / (2.0 * e[l]); 00332 r = sqrt((g * g) + 1.0); 00333 g = d[n] - d[l] + e[l] / (g + SIGN(r, g)); 00334 s = c = 1.0; 00335 p = 0.0; 00336 for (i = n-1; i >= l; i--) 00337 { 00338 f = s * e[i]; 00339 b = c * e[i]; 00340 if (fabs(f) >= fabs(g)) 00341 { 00342 c = g / f; 00343 r = sqrt((c * c) + 1.0); 00344 e[i+1] = f * r; 00345 c *= (s = 1.0/r); 00346 } 00347 else 00348 { 00349 s = f / g; 00350 r = sqrt((s * s) + 1.0); 00351 e[i+1] = g * r; 00352 s *= (c = 1.0/r); 00353 } 00354 g = d[i+1] - p; 00355 r = (d[i] - g) * s + 2.0 * c * b; 00356 p = s * r; 00357 d[i+1] = g + p; 00358 g = c * r - b; 00359 for (k = 0; k < m; k++) 00360 { 00361 f = z(k,i+1); 00362 z(k,i+1) = s * z(k,i) + c * f; 00363 z(k,i) = c * z(k,i) - s * f; 00364 } 00365 } 00366 d[l] = d[l] - p; 00367 e[l] = g; 00368 e[n] = 0.0; 00369 } 00370 } while (n != l); 00371 } 00372 } 00373 00374 00375 00376 00377 00378 00379 00380