Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/PCA.cpp
Go to the documentation of this file.
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