Marsyas
0.6.0-alpha
|
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 //}