Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/NormCut.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2010 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 "../common_source.h"
00020 #include "NormCut.h"
00021 #include <marsyas/NumericLib.h>
00022 
00023 #include <iostream>
00024 #include <cmath>
00025 #include <algorithm>
00026 
00027 //#define MTLB_DBG_LOG
00028 
00029 using std::ostringstream;
00030 using std::cout;
00031 using std::min;
00032 
00033 using std::endl;
00034 
00035 using namespace Marsyas;
00036 
00037 NormCut::NormCut(mrs_string name):MarSystem("NormCut", name)
00038 {
00039   numClusters_ = -1;
00040   paramOffset_ = 0.5;
00041   paramVerbose_ = 3;
00042   paramMaxiterations_ = 20;
00043   paramEigsErrorTolerance_ = 0.000001;
00044 
00045   addControls();
00046 }
00047 
00048 NormCut::NormCut(const NormCut& a) : MarSystem(a)
00049 {
00050   numClusters_ = a.numClusters_;
00051   paramOffset_ = a.paramOffset_;
00052   paramVerbose_ = a.paramVerbose_;
00053   paramMaxiterations_ = a.paramMaxiterations_;
00054   paramEigsErrorTolerance_ = a.paramEigsErrorTolerance_;
00055 
00056   ctrl_numClusters_ = getctrl("mrs_natural/numClusters");
00057 }
00058 
00059 NormCut::~NormCut()
00060 {
00061 
00062 }
00063 
00064 MarSystem*
00065 NormCut::clone() const
00066 {
00067   return new NormCut(*this);
00068 }
00069 
00070 void
00071 NormCut::addControls()
00072 {
00073   addctrl("mrs_natural/numClusters", 2, ctrl_numClusters_);
00074   setctrlState("mrs_natural/numClusters", true);
00075 
00076   addctrl("mrs_real/offset", 0.5);
00077   setctrlState("mrs_real/offset", true);
00078 
00079   addctrl("mrs_natural/verbose", 3);
00080   setctrlState("mrs_natural/verbose", true);
00081 
00082   addctrl("mrs_natural/maxIters", 20);
00083   setctrlState("mrs_natural/maxIters", true);
00084 
00085   addctrl("mrs_real/eigsErrorTol", 0.000001);
00086   setctrlState("mrs_real/eigsErrorTol", true);
00087 }
00088 
00089 void
00090 NormCut::myUpdate(MarControlPtr sender)
00091 {
00092   (void) sender;  //suppress warning of unused parameter(s)
00093   ctrl_onSamples_->setValue(ctrl_inSamples_, NOUPDATE);
00094   ctrl_onObservations_->setValue(1, NOUPDATE);
00095   ctrl_osrate_->setValue(ctrl_israte_, NOUPDATE);
00096   ctrl_onObsNames_->setValue("PeakLabels", NOUPDATE);
00097 
00098   //check if there was a change in the input similarity matrix (i.e. in the nr of elements)
00099   //or if the number of clusters to detect was changed
00100   // [AL]: [!] this means that the number of inObservations has to be set *before* the number of clusters,
00101   // [AL] otherwise nCutDiscrete_   and nCutEigVectors_ have the wrong dimension !!!
00102   if(numClusters_ != ctrl_numClusters_->to<mrs_natural>() || onSamples_ != inSamples_)
00103   {
00104     numClusters_ = ctrl_numClusters_->to<mrs_natural>();
00105 
00106     nCutDiscrete_.stretch(numClusters_ * ctrl_inObservations_->to<mrs_natural>());
00107     nCutEigVectors_.stretch(numClusters_ * ctrl_inObservations_->to<mrs_natural>());
00108     nCutEigValues_.stretch(numClusters_);
00109   }
00110 }
00111 
00112 void
00113 NormCut::myProcess(realvec& in, realvec& out)
00114 {
00115   mrs_natural t,o;
00116 
00117 #ifdef MARSYAS_MATLAB
00118 #ifdef MTLB_DBG_LOG
00119   MATLAB_PUT(in, "in");
00120   MATLAB_EVAL("figure(71),imagesc(in',[0 1]),colorbar");
00121 #endif
00122 #endif
00123   //check if there is any data at the input, otherwise do nothing
00124   if(in.getSize() == 0 || numClusters_ == 0)
00125   {
00126     //MRSWARN("NormCut::myProcess - empty input!");
00127     out.setval(-1.0);
00128     return;
00129   }
00130   if(in.getSize() == 1 || numClusters_ == 0)
00131   {
00132     //MRSWARN("NormCut::myProcess - empty input!");
00133     out.setval(0);
00134     return;
00135   }
00136 
00137   out.setval(0); //should it be -1.0 instead? [!] FIXME
00138 
00139   //ncut( n, W, nbcluster, NcutEigenvectors, NcutEigenvalues, params );
00140   //discretisation( n, nbcluster, NcutEigenvectors, NcutDiscrete, params );
00141 
00142   ncut(inObservations_, in, numClusters_, nCutEigVectors_, nCutDiscrete_ );
00143   discretisation(inObservations_, numClusters_, nCutEigVectors_, nCutDiscrete_ );
00144 
00145   for( o=0 ; o<inObservations_ ; o++ )
00146   {
00147     for( t=0 ; t<numClusters_ ; t++ )
00148     {
00149       // Initialize NcutDiscrete as we go
00150       if( nCutDiscrete_(o*(numClusters_)+t) == 1.0 )
00151         out(0,o) = t;
00152       //cout << "Peak " << o << " -> cluster = " << out(0,o) << endl;
00153     }
00154   }
00155 
00156   //cout << out << endl;
00157 
00158   //   //get a local copy of the current gain control value
00159   //   //(it will be used for this entire processing, even if it's
00160   //   //changed by someone else, e.g. by a different thread)
00161   //   mrs_real gainValue = ctrl_gain_->to<mrs_real>();
00162   //
00163   //   // It is important to loop over both observations
00164   //   // and channels so that for example a gain can be
00165   //   // applied to multi-channel signals
00166   //   for (o=0; o < inObservations_; o++)
00167   //   {
00168   //      for (t = 0; t < inSamples_; t++)
00169   //      {
00170   //         //apply gain to all channels
00171   //         out(o,t) = gainValue * in(o,t);
00172   //      }
00173   //   }
00174 }
00175 
00176 
00177 /*
00178  * Function :    ncutW (driver function)
00179  * Description : Clusters data items using the normalized cut algorithm given an n x n similarity matrix
00180  *
00181  * Arguments: mrs_natural *n                    -- size of data set (input)
00182  *            realvec &W                 -- n x n symmetric similarity matrix (input)
00183  *            mrs_natural *nbcluster            -- desired number of clusters (input)
00184  *            dataNcutPtr params                   -- algorithm parameters (input)
00185  *            realvec &NcutDiscrete      -- nbcluster x n clustering results (output)
00186  *            realvec &NcutEigenvectors  -- resulting eigenvectors (output)
00187  *            realvec &NcutEigenvalues   -- resulting eigenvalues (output)
00188  */
00189 //void
00190 //NormCut::ncutW(mrs_natural *n, realvec &W, mrs_natural *nbcluster, realvec &NcutDiscrete, realvec &NcutEigenvectors, realvec &NcutEigenvalues, dataNcutParamsPtr params)
00191 //{
00192 //   ncut( n, W, nbcluster, NcutEigenvectors, NcutEigenvalues, params );
00193 //  /* prmrs_natural(NcutEigenvectors, *n, *nbcluster);
00194 //   prmrs_natural(NcutEigenvalues, *nbcluster, 1);*/
00195 //   discretisation( n, nbcluster, NcutEigenvectors, NcutDiscrete, params );
00196 //  // prmrs_natural(NcutDiscrete, *n, *nbcluster);
00197 //}
00198 
00199 void
00200 NormCut::ncut(mrs_natural n, realvec &W, mrs_natural nbcluster, realvec &NcutEigenvectors, realvec &NcutEigenvalues)
00201 {
00202   //dataNcutParams param;
00203   realvec dinvsqrt(n);
00204   mrs_real ulp;
00205   realvec P(n*(n));
00206 
00207   mrs_real norm;
00208   mrs_real sqrtn = sqrt((mrs_real) n);
00209 
00210   nbcluster = min(nbcluster,n);
00211 
00212   // params required by dsyevr_ but are not initialized or used
00213   realvec evals(n);
00214   realvec mrs_naturalerm(n);
00215 
00216   mrs_natural i,j;
00217 
00218   //   if( params != NULL )
00219   //      param = *params;
00220   //   else
00221   //   {
00222   //      param.offset = 0.5;
00223   //      param.verbose = 3;
00224   //      param.maxiterations = 20;
00225   //      param.eigsErrorTolerance = 0.000001;
00226   //   }
00227   // SHOULD BE USED IN CONSTRUCTOR !!! JEN
00228 
00229   // Check for matrix symmetry
00230   for( i=0 ; i<n ; ++i )
00231     for( j=0 ; j<n ; j++ ) {
00232       if( W(i*(n)+j) > 1.0 )
00233       {
00234         //cout << W;
00235         //ERROR("W(i, j) values should be <= 1");
00236         MRSWARN("NormCut::ncut() - input values should be <= 1 : delta @(" << i << "," << j << ") = " << W(i*(n)+j)-1.0);
00237       }
00238       if( W(i*(n)+j) != W(j*(n)+i) )
00239       {
00240         //cout << W;
00241         //ERROR("W is not symmetric");
00242         MRSWARN("NormCut::ncut - input matrix is not symmetric!");
00243       }
00244       P(i*(n)+j)=0.;
00245     }
00246 
00247   //ulp = NumericLib::machp("Epsilon") * NumericLib::machp("Base"); // unit in last place
00248   ulp = std::numeric_limits<double>::epsilon() * std::numeric_limits<double>::radix;
00249 
00250   //sum each row
00251   for( i=0 ; i<n ; ++i )
00252   {
00253     // can sum the columns since it's symmetric... and its faster
00254     dinvsqrt(i) = 2.*paramOffset_;
00255     for( j=0 ; j<n ; j++ )
00256       dinvsqrt(i) += W(i*(n)+j);
00257     dinvsqrt(i) = 1./sqrt(dinvsqrt(i)+ulp);
00258     MRSASSERT(dinvsqrt(i) == dinvsqrt(i));
00259   }
00260 
00261   // P <- dinvsqrt*dinvsqrt'
00262   for( i=0 ; i<n ; ++i )
00263     for( j=i ; j<n ; j++ )
00264       P(i*(n)+j) = dinvsqrt(i)*dinvsqrt(j);
00265 
00266   // 1. Add offset to diagonal of W
00267   // 2. P <- P .* W (only lower triangle of P computed)
00268   for( j=0 ; j<n ; j++ ) {
00269     P(j*(n)+j) = P(j*(n)+j)*( W(j*(n)+j) + paramOffset_ );
00270     for( i=j+1 ; i<n ; ++i ) {
00271       P(j*(n)+i) = P(j*(n)+i)*W(j*(n)+i);
00272     }
00273   }
00274 
00275   NumericLib::tred2(P, n, evals, mrs_naturalerm );
00276   NumericLib::tqli( evals, mrs_naturalerm, n, P );
00277   /* evals now contains the eigenvalues,
00278      P now contains the associated eigenvectors. */
00279 
00280   // Must reverse eigenvectors and eigenvalues because
00281   // tqli returns them in ascending order, rather than
00282   // descending
00283   for( j=0 ; j<nbcluster ; j++ )
00284   {
00285     for( i=0 ; i<n ; ++i )
00286     {
00287       NcutEigenvectors(j*(n)+i) = P((n-j-1)*(n)+i);
00288       MRSASSERT(NcutEigenvectors(j*n + i) == NcutEigenvectors(j*n + i));
00289     }
00290     NcutEigenvalues(j) = evals(n-j-1);
00291   }
00292 
00293   for( j=0 ; j<nbcluster ; j++ )
00294     for( i=0 ; i<n ; ++i )
00295     {
00296       NcutEigenvectors(j*(n)+i) = NcutEigenvectors(j*(n)+i)*dinvsqrt(i);
00297       MRSASSERT(NcutEigenvectors(j*n + i) == NcutEigenvectors(j*n + i));
00298     }
00299 
00300 
00301   for( j=0 ; j<nbcluster ; j++ ) {
00302     norm=0.;
00303     for( i=0 ; i<n ; ++i )
00304       norm += NcutEigenvectors(j*(n) + i)*NcutEigenvectors(j*(n) + i);
00305     norm = sqrt(norm);
00306     norm = sqrtn / norm;
00307     if( NcutEigenvectors(j*(n)) >= 0. )
00308       for( i=0 ; i<n ; ++i )
00309       {
00310         NcutEigenvectors(j*(n) + i) *= -norm;
00311         MRSASSERT(NcutEigenvectors(j*n + i) == NcutEigenvectors(j*n + i));
00312       }
00313 
00314     else
00315       for( i=0 ; i<n ; ++i )
00316       {
00317         NcutEigenvectors(j*(n) + i) *= norm;
00318         MRSASSERT(NcutEigenvectors(j*n + i) == NcutEigenvectors(j*n + i));
00319       }
00320 
00321   }
00322 
00323 }
00324 
00325 void NormCut::discretisation(mrs_natural n,  mrs_natural nbcluster, realvec &NcutEigenvectors, realvec &NcutDiscrete)
00326 {
00327   realvec vm(n);
00328   realvec R(nbcluster*(nbcluster));
00329   realvec EVtimesR(n*(nbcluster));
00330   realvec c(n);
00331   realvec tmp(n);
00332   double minval;
00333   int mini  = 0; //AL: I don't really know what implications this has - but otherwise the variable might be used uninitialized...
00334   realvec EVDtimesEV(nbcluster*(nbcluster));
00335   double NcutValue;
00336   double lastObjectiveValue=0;
00337 
00338   // output params for dgesdd_
00339   realvec s(nbcluster+1); // array of singluar values in descending order
00340   realvec U(nbcluster*(nbcluster)); // U
00341   realvec V(nbcluster*(nbcluster)); // V
00342 
00343   //double ulp = NumericLib::machp("Epsilon") * NumericLib::machp("Base"); // unit in last place
00344   double ulp = std::numeric_limits<double>::epsilon() * std::numeric_limits<double>::radix;
00345 
00346   int nbIterDiscrete = 0;
00347   bool exitLoop = false;
00348 
00349   int randn;
00350   int i,j,k;
00351 
00352   //dataNcutParams param;
00353   //
00354   //   if( params != NULL )
00355   //      param = *params;
00356   //   else
00357   //   {
00358   //      param.offset = 0.5;
00359   //      param.verbose = 3;
00360   //      param.maxiterations = 20;
00361   //      param.eigsErrorTolerance = 0.000001;
00362   //   }
00363 
00364   for( i=0 ; i<n ; ++i ) {
00365     vm(i) = 0;
00366     for( j=0 ; j<nbcluster ; j++ )
00367     {
00368       vm(i) += NcutEigenvectors(j*(n)+i)*NcutEigenvectors(j*(n)+i);
00369       MRSASSERT(vm(i) == vm(i));
00370     }
00371     vm(i) = sqrt(vm(i));
00372     for( j=0 ; j<nbcluster ; j++ ) {
00373       if (vm(i) > 0)
00374         NcutEigenvectors(j*(n)+i) /= vm(i);
00375       else
00376         NcutEigenvectors(j*(n)+i)  = 0.;
00377       MRSASSERT(NcutEigenvectors(j*(n)+i) == NcutEigenvectors(j*(n)+i));
00378     }
00379     c(i) = 0;
00380   }
00381 
00382   // randn = (int)floor(0.5+(n-1)*(0.42863169099606)); // rand()/RAND_MAX for large decimal -- temporary while debugging
00383   randn = 0; // although I have absolutely no clue what this means, everything seems to work better with this being 0 (AL)
00384 
00385   for( i=0 ; i<nbcluster ; ++i ) {
00386     R(i) = NcutEigenvectors(i*(n)+randn);
00387     MRSASSERT(R(i) == R(i));
00388     for( j=0 ; j<nbcluster ; j++ )
00389       U(i*(nbcluster)+j) = 0.0;
00390   }
00391 
00392   for( j=1 ; j<nbcluster ; j++ )
00393   {
00394     minval = MAXREAL;
00395 
00396     for( i=0 ; i<n ; ++i )
00397     {
00398       tmp(i) = 0.;
00399       for( k=0 ; k<nbcluster ; k++ )
00400         tmp(i) += NcutEigenvectors(k*(n)+i)*R((j-1)*(nbcluster)+k);
00401     }
00402 
00403     for( i=0 ; i<n ; ++i ) {
00404       c(i) += fabs(tmp(i));
00405       if( c(i) < minval )
00406       {
00407         minval = c(i);
00408         mini = i;
00409       }
00410     }
00411 
00412     for( i=0 ; i<nbcluster ; ++i ) {
00413       R(j*(nbcluster)+i) = NcutEigenvectors(i*(n) + mini);
00414       MRSASSERT(R(j*(nbcluster)+i) == R(j*(nbcluster)+i));
00415     }
00416   }
00417 
00418   //     cout << endl;
00419   //  print(R, *nbcluster, *nbcluster);
00420   while( !exitLoop )
00421   {
00422     nbIterDiscrete += 1;
00423 
00424     for( i=0 ; i<n ; ++i ) {
00425       for( j=0 ; j<nbcluster ; j++ ) {
00426         EVtimesR(j*(n)+i) = 0.;
00427         for( k=0 ; k<nbcluster ; k++)
00428         {
00429           EVtimesR(j*(n)+i) += NcutEigenvectors(k*(n)+i)*R(j*(nbcluster)+k);
00430           MRSASSERT(EVtimesR(j*(n)+i) == EVtimesR(j*(n)+i));
00431         }
00432       }
00433     }
00434 
00435     discretisationEigenvectorData(n, nbcluster, EVtimesR,NcutDiscrete);
00436 
00437     for( i=0 ; i<nbcluster ; ++i ) {
00438       for( j=0 ; j<nbcluster ; j++ ) {
00439         EVDtimesEV(j*(nbcluster)+i) = 0.;
00440         for( k=0 ; k<n ; k++)
00441         {
00442           EVDtimesEV(j*(nbcluster)+i) += NcutDiscrete(k*(nbcluster)+i)*NcutEigenvectors(j*(n)+k);
00443           MRSASSERT(EVDtimesEV(j*(nbcluster)+i) == EVDtimesEV(j*(nbcluster)+i));
00444         }
00445       }
00446     }
00447     //  cout << endl;
00448     //  print(EVtimesR, *n, *nbcluster);
00449     // cout << endl;
00450     // print(EVDtimesEV, *nbcluster, *nbcluster);
00451     // Do SVD : store U,S,V
00452     NumericLib::svd( nbcluster, nbcluster, EVDtimesEV, U, V, s );
00453 
00454     // 2*( n-trace(S) )
00455     NcutValue = 0.;
00456     for( i=0 ; i<nbcluster ; ++i )
00457       NcutValue += s(i);
00458     NcutValue = 2*( n - NcutValue );
00459 
00460     if( fabs(NcutValue - lastObjectiveValue) < ulp || nbIterDiscrete > paramMaxiterations_ )
00461     {
00462       exitLoop = 1;
00463     }
00464     else
00465     {
00466       lastObjectiveValue = NcutValue;
00467       // R= V * U';
00468       for( i=0 ; i<nbcluster ; ++i ) {
00469         for( j=0 ; j<nbcluster ; j++ ) {
00470           R(j*(nbcluster)+i) = 0.;
00471           for( k=0 ; k<nbcluster ; k++)
00472             R(j*(nbcluster)+i) += V(k*(nbcluster)+i)*U(k*(nbcluster)+j);
00473         }
00474       }
00475 
00476     }
00477   }
00478 }
00479 
00480 void
00481 NormCut::discretisationEigenvectorData(mrs_natural n,  mrs_natural nbcluster, realvec &V, realvec &Vdiscrete)
00482 {
00483   mrs_real maxval;
00484   mrs_natural maxi = 0; //AL: I don't really know what implications this has - but otherwise the variable might be used uninitialized...;
00485 
00486   mrs_natural i,j;
00487 
00488   for( i=0 ; i<n ; ++i ) {
00489     maxval = -MAXREAL; // (Mathieu ??)
00490 
00491     // Find largest value in the ith row of the eigenvectors
00492     for( j=0 ; j<nbcluster ; j++ )
00493     {
00494       // Initialize NcutDiscrete as we go
00495       Vdiscrete(i*(nbcluster)+j) = 0.0;
00496 
00497       if( V(j*(n)+i) > maxval )
00498       {
00499         maxval = V(j*(n)+i);
00500         maxi = j;
00501       }
00502     }
00503     // Data i belongs to cluster maxi
00504     Vdiscrete(i*(nbcluster)+maxi) = 1.0;
00505     // Vdiscrete(i) = maxi;  // (Mathieu !!)
00506   }
00507 }
00508 
00509 
00510 /*
00511  * Display column-oriented matrices
00512  * Note : Matrices are stored in column-wise order for compatibility with CLAPACK (translated from FORTRAN)
00513  * Note2: n=1 --> display column vector         n=-1 --> display row vector
00514  */
00515 void
00516 NormCut::prmrs_natural( realvec &A, mrs_natural m , mrs_natural n )
00517 {
00518   mrs_natural i,j;
00519 
00520   if( n > 0 )
00521   {
00522     for(i=0; i < m; ++i)
00523     {
00524       for (j=0; j < n; j++)
00525         cout << A(j*m+i) << "\t";
00526       cout << endl;
00527     }
00528   }
00529   else if( n == -1 )
00530   {
00531     for( i=0 ; i<m ; ++i )
00532       cout << A(i) << "\t";
00533     cout << endl;
00534   }
00535 }
00536 
00537 /*
00538  * Display column-oriented matrices
00539  * Note : Matrices are stored in column-wise order for compatibility with CLAPACK (translated from FORTRAN)
00540  * Note2: n=1 --> display column vector         n=-1 --> display row vector
00541  */
00542 
00543 void
00544 NormCut::print( realvec &A, int m , int n )
00545 {
00546   int i,j;
00547 
00548   if( n > 0 )
00549   {
00550     for(i=0; i < m; ++i)
00551     {
00552       for (j=0; j < n; j++)
00553         cout << A(j*m+i) << "\t";
00554       cout << endl;
00555     }
00556   }
00557   else if( n == -1 )
00558   {
00559     for( i=0 ; i<m ; ++i )
00560       cout << A(i) << "\t";
00561     cout << endl;
00562   }
00563 }
00564 
00565 //
00566 //
00567 //void
00568 //Marsyas::discretisation(mrs_natural *n,  mrs_natural *nbcluster, realvec &NcutEigenvectors, realvec &NcutDiscrete, dataNcutParamsPtr params)
00569 //{
00570 //   realvec vm(*n);
00571 //   realvec R(*nbcluster*(*nbcluster));
00572 //   realvec EVtimesR(*n*(*nbcluster));
00573 //   realvec c(*n);
00574 //   realvec tmp(*n);
00575 //   mrs_real minval;
00576 //   mrs_natural mini;
00577 //   realvec EVDtimesEV(*nbcluster*(*nbcluster)); // (Mathieu ??)
00578 //   mrs_real NcutValue;
00579 //   mrs_real lastObjectiveValue=0;
00580 //
00581 //   // output params for dgesdd_
00582 //   realvec s(*nbcluster+1); // array of singluar values in descending order
00583 //   realvec U(*nbcluster*(*nbcluster)); // U
00584 //   realvec V(*nbcluster*(*nbcluster)); // V
00585 //
00586 //   mrs_real ulp = NumericLib::machp("Epsilon") * NumericLib::machp("Base"); // unit in last place
00587 //
00588 //   mrs_natural nbIterDiscrete = 0;
00589 //   bool exitLoop = false;
00590 //
00591 //   mrs_natural randn;
00592 //   mrs_natural i,j,k;
00593 //
00594 //   dataNcutParams param;
00595 //
00596 //   if( params != NULL )
00597 //      param = *params;
00598 //   else
00599 //   {
00600 //      param.offset = 0.5;
00601 //      param.verbose = 3;
00602 //      param.maxiterations = 20;
00603 //      param.eigsErrorTolerance = 0.000001;
00604 //   }
00605 //
00606 //   for( i=0 ; i<*n ; ++i ){
00607 //      vm(i) = 0;
00608 //      for( j=0 ; j<*nbcluster ; j++ )
00609 //         vm(i) += NcutEigenvectors(j*(*n)+i)*NcutEigenvectors(j*(*n)+i);
00610 //      vm(i) = sqrt(vm(i));
00611 //      for( j=0 ; j<*nbcluster ; j++ ){
00612 //         NcutEigenvectors(j*(*n)+i) /= vm(i);
00613 //      }
00614 //      c(i) = 0;
00615 //   }
00616 //
00617 //
00618 //   randn = (mrs_natural) floor(0.5+(*n-1)*(0.42863169099606)); // rand()/RAND_MAX for large decimal -- temporary while debugging
00620 //
00621 //
00622 //   for( i=0 ; i<*nbcluster ; ++i ){
00623 //      R(i) = NcutEigenvectors(i*(*n)+randn);
00624 //      for( j=0 ; j<*nbcluster ; j++ )
00625 //         U(i*(*nbcluster)+j) = 0.0;
00626 //   }
00627 //
00628 //   for( j=1 ; j<*nbcluster ; j++ )
00629 //   {
00630 //      minval = MAXREAL;
00631 //
00632 //      for( i=0 ; i<*n ; ++i )
00633 //      {
00634 //         tmp(i) = 0.;
00635 //         for( k=0 ; k<*nbcluster ; k++ )
00636 //            tmp(i) += NcutEigenvectors(k*(*n)+i)*R((j-1)*(*nbcluster)+k);
00637 //      }
00638 //
00639 //      for( i=0 ; i<*n ; ++i ){
00640 //         c(i) += fabs(tmp(i));
00641 //         if( c(i) < minval )
00642 //         {
00643 //            minval = c(i);
00644 //            mini = i;
00645 //         }
00646 //      }
00647 //
00648 //      for( i=0 ; i<*nbcluster ; ++i ){
00649 //         R(j*(*nbcluster)+i) = NcutEigenvectors(i*(*n) + mini);
00650 //      }
00651 //   }
00652 //   cout << endl;
00653 //   prmrs_natural(R, *nbcluster, *nbcluster);
00654 //
00655 //
00656 //   while( !exitLoop )
00657 //   {
00658 //      nbIterDiscrete += 1;
00659 //
00660 //      for( i=0 ; i<*n ; ++i ){
00661 //         for( j=0 ; j<*nbcluster ; j++ ){
00662 //            EVtimesR(j*(*n)+i) = 0.;
00663 //            for( k=0 ; k<*nbcluster ; k++)
00664 //               EVtimesR(j*(*n)+i) += NcutEigenvectors(k*(*n)+i)*R(j*(*nbcluster)+k);
00665 //         }
00666 //      }
00667 //          cout << endl;
00668 //      prmrs_natural(EVtimesR, *n, *nbcluster);
00669 //
00670 //      discretisationEigenvectorData(n, nbcluster, EVtimesR,NcutDiscrete);
00671 //
00672 //      for( i=0 ; i<*nbcluster ; ++i ){
00673 //         for( j=0 ; j<*nbcluster ; j++ ){ // (Mathieu n replaced by nbcluster)
00674 //            EVDtimesEV(j*(*nbcluster)+i) = 0.;
00675 //            for( k=0 ; k<*n ; k++)
00676 //               EVDtimesEV(j*(*nbcluster)+i) += NcutDiscrete(k*(*nbcluster)+i)*NcutEigenvectors(j*(*n)+k);
00677 //         }
00678 //      }
00679 //          cout << endl;
00680 //      prmrs_natural(EVDtimesEV, *nbcluster, *nbcluster);
00681 //      // Do SVD : store U,S,V
00682 //          NumericLib::svd( *nbcluster, *nbcluster, EVDtimesEV, U, V, s );
00683 //
00684 //      // 2*( n-trace(S) )
00685 //      NcutValue = 0.;
00686 //      for( i=0 ; i<*nbcluster ; ++i )
00687 //         NcutValue += s(i);
00688 //      NcutValue = 2*( *n - NcutValue );
00689 //
00690 //      if( fabs(NcutValue - lastObjectiveValue) < ulp || nbIterDiscrete > param.maxiterations )
00691 //      {
00692 //         exitLoop = 1;
00693 //      }
00694 //      else
00695 //      {
00696 //         lastObjectiveValue = NcutValue;
00697 //         // R= V * U';
00698 //         for( i=0 ; i<*nbcluster ; ++i ){
00699 //            for( j=0 ; j<*nbcluster ; j++ ){
00700 //               R(j*(*nbcluster)+i) = 0.;
00701 //               for( k=0 ; k<*nbcluster ; k++)
00702 //                  R(j*(*nbcluster)+i) += V(k*(*nbcluster)+i)*U(k*(*nbcluster)+j);
00703 //            }
00704 //         }
00705 //
00706 //      }
00707 //   }
00708 //}