PLearn 0.1
|
00001 // -*- C++ -*- 00002 00003 // ManifoldParzen2.cc 00004 // 00005 // Copyright (C) 2003 Pascal Vincent, Julien Keable 00006 // Copyright (C) 2005 University of Montreal 00007 // 00008 // Redistribution and use in source and binary forms, with or without 00009 // modification, are permitted provided that the following conditions are met: 00010 // 00011 // 1. Redistributions of source code must retain the above copyright 00012 // notice, this list of conditions and the following disclaimer. 00013 // 00014 // 2. Redistributions in binary form must reproduce the above copyright 00015 // notice, this list of conditions and the following disclaimer in the 00016 // documentation and/or other materials provided with the distribution. 00017 // 00018 // 3. The name of the authors may not be used to endorse or promote 00019 // products derived from this software without specific prior written 00020 // permission. 00021 // 00022 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR 00023 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 00024 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN 00025 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 00027 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00028 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00029 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00030 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00031 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00032 // 00033 // This file is part of the PLearn library. For more information on the PLearn 00034 // library, go to the PLearn Web site at www.plearn.org 00035 00036 00038 #include "ManifoldParzen2.h" 00039 00040 #include <plearn/math/plapack.h> 00041 #include <plearn/base/general.h> 00042 #include <plearn/math/TMat.h> 00043 #include <plearn/math/TMat_maths.h> 00044 #include <plearn/math/BottomNI.h> 00045 #include <plearn/vmat/VMat_computeNearestNeighbors.h> 00046 #include <plearn/vmat/MemoryVMatrix.h> 00047 00048 namespace PLearn { 00049 00050 PLEARN_IMPLEMENT_OBJECT(ManifoldParzen2, 00051 "Manifold Parzen density estimate ", 00052 "Parzen Window algorithm, where the covariance matrices of the gaussians are\n" 00053 "computed using a knn kernel estimate. Also, only the ncomponents principal\n" 00054 "eigen vector and eigen values of the covariance matrices are stored.\n" 00055 ); 00056 00058 // ManifoldParzen2 // 00060 ManifoldParzen2::ManifoldParzen2() 00061 : nneighbors(4), 00062 ncomponents(1), 00063 use_last_eigenval(true), 00064 scale_factor(1), 00065 learn_mu(false) 00066 { 00067 type = "general"; 00068 nstages = 1; 00069 } 00070 00071 ManifoldParzen2::ManifoldParzen2(int the_nneighbors, int the_ncomponents, bool use_last_eigenvalue, real the_scale_factor) 00072 : nneighbors(the_nneighbors),ncomponents(the_ncomponents),use_last_eigenval(true),scale_factor(the_scale_factor) 00073 { 00074 type = "general"; 00075 nstages = 1; 00076 } 00077 00078 void ManifoldParzen2::build() 00079 { 00080 inherited::build(); 00081 build_(); 00082 } 00083 00085 // declareOptions // 00087 void ManifoldParzen2::declareOptions(OptionList& ol) 00088 { 00089 declareOption(ol,"nneighbors", &ManifoldParzen2::nneighbors, OptionBase::buildoption, 00090 "Number of neighbors for covariance matrix estimation."); 00091 00092 declareOption(ol,"ncomponents", &ManifoldParzen2::ncomponents, OptionBase::buildoption, 00093 "Number of components to store from the PCA."); 00094 00095 declareOption(ol,"use_last_eigenval", &ManifoldParzen2::use_last_eigenval, OptionBase::buildoption, 00096 "Indication that the last eigenvalue should be used for the remaining directions' variance."); 00097 00098 declareOption(ol,"learn_mu", &ManifoldParzen2::learn_mu, OptionBase::buildoption, 00099 "Indication that the difference vector between the training points and the gaussians should be learned.\n" 00100 "By default, the gaussians are centered at the training points."); 00101 00102 declareOption(ol,"global_lambda0", &ManifoldParzen2::global_lambda0, OptionBase::buildoption, 00103 "If use_last_eigenvalue is false, used value for the minimum variance in all directions"); 00104 00105 declareOption(ol,"scale_factor", &ManifoldParzen2::scale_factor, OptionBase::buildoption, 00106 "Scale factor"); 00107 00108 // Now call the parent class' declareOptions 00109 inherited::declareOptions(ol); 00110 00111 // Redeclare some parent's options. 00112 00113 redeclareOption(ol,"L", &ManifoldParzen2::L, OptionBase::learntoption, 00114 "Automatically set (to train_set->length())."); 00115 00116 redeclareOption(ol,"n_eigen", &ManifoldParzen2::n_eigen, OptionBase::learntoption, 00117 "Automatically set during training."); 00118 00119 redeclareOption(ol,"type", &ManifoldParzen2::type, OptionBase::nosave, 00120 "Automatically set (to 'general')."); 00121 00122 redeclareOption(ol,"alpha_min", &ManifoldParzen2::alpha_min, OptionBase::nosave, 00123 "Not used."); 00124 00125 redeclareOption(ol,"kmeans_iterations", &ManifoldParzen2::kmeans_iterations, OptionBase::nosave, 00126 "Not used."); 00127 00128 } 00129 00130 void ManifoldParzen2::build_() 00131 {} 00132 00133 void ManifoldParzen2::makeDeepCopyFromShallowCopy(CopiesMap& copies) 00134 { 00135 inherited::makeDeepCopyFromShallowCopy(copies); 00136 00137 // ### Call deepCopyField on all "pointer-like" fields 00138 // ### that you wish to be deepCopied rather than 00139 // ### shallow-copied. 00140 // ### ex: 00141 // deepCopyField(trainvec, copies); 00142 00143 // ### Remove this line when you have fully implemented this method. 00144 //PLERROR("ManifoldParzen2::makeDeepCopyFromShallowCopy not fully (correctly) implemented yet!"); 00145 } 00146 00147 // This is an efficient version of the most basic nearest neighbor search, using a Mat and euclidean distance 00148 void computeNearestNeighbors(Mat dataset, Vec x, Mat& neighbors, int ignore_row=-1) 00149 { 00150 int K = neighbors.length(); // how many neighbors do we want? 00151 BottomNI<real> neighbs(K); 00152 for(int i=0; i<dataset.length(); i++) 00153 if(i!=ignore_row) 00154 neighbs.update(powdistance(dataset(i),x), i); 00155 neighbs.sort(); 00156 TVec< pair<real,int> > indices = neighbs.getBottomN(); 00157 int nonzero=0; 00158 for(int k=0; k<K; k++) 00159 { 00160 if(indices[k].first>0) 00161 nonzero++; 00162 neighbors(k) << dataset(indices[k].second); 00163 } 00164 if(nonzero==0) 00165 PLERROR("All neighbors had 0 distance. Use more neighbors. (There were %i other patterns with same values)",neighbs.nZeros()); 00166 } 00167 00168 // computes the first components.length() principal components of the rows of datset. 00169 // Result will be put in the rows of components (which is expected to have the appropriate size). 00170 // Computed components do not have a norm of 1, but rather a norm corresponding to the eigenvalues 00171 // (indicating how important the component is...) 00172 void computePrincipalComponents(Mat dataset, Vec& eig_values, Mat& eig_vectors, real global_lambda0) 00173 { 00174 #ifdef BOUNDCHECK 00175 if(eig_vectors.width()!=dataset.width()) 00176 PLERROR("In computePrincipalComponents eig_vectors and dataset must have same width"); 00177 if(eig_values.length() != eig_vectors.length()) 00178 PLERROR("In computePrincipalComponents eig_values vec and eig_vectors mat must have same length"); 00179 #endif 00180 00181 static Mat covar; 00182 int ncomp = eig_values.length(); // number of components we want 00183 covar.resize(dataset.width(), dataset.width()); 00184 transposeProduct(covar, dataset,dataset); 00185 eigenVecOfSymmMat(covar, ncomp, eig_values, eig_vectors); 00186 for (int i=0;i<eig_values.length();i++) 00187 { 00188 if (eig_values[i]<0) 00189 eig_values[i] = 0; 00190 eig_values[i] = eig_values[i]/dataset.length() + global_lambda0; 00191 } 00192 } 00193 00194 void computeLocalPrincipalComponents(Mat& dataset, int which_pattern, Mat& delta_neighbors, Vec& eig_values, Mat& eig_vectors, Vec& mean, bool learn_mu=false,real global_lambda0=0) 00195 { 00196 Vec center_ = dataset(which_pattern); 00197 if (center_.hasMissing()) 00198 PLERROR("dataset row %d has missing values!", which_pattern); 00199 computeNearestNeighbors(dataset, center_, delta_neighbors, which_pattern); 00200 if(learn_mu) 00201 { 00202 mean.resize(delta_neighbors.width()); // Hugo: the mean should be the current point... 00203 columnMean(delta_neighbors, mean); 00204 } 00205 delta_neighbors -= mean; 00206 computePrincipalComponents(delta_neighbors, eig_values, eig_vectors,global_lambda0); 00207 } 00208 00209 void ManifoldParzen2::train() 00210 { 00211 Mat trainset(train_set); 00212 int l = train_set.length(); 00213 int w = train_set.width(); 00214 eigenval_copy.resize(ncomponents+1); 00215 row.resize(w); 00216 00217 L = l; 00218 // D = ncomponents; // This appears to cause a bug, D is expected by GaussMix to be the total input dimensionality, which may be different from ncomponents (Pascal V.) 00219 n_eigen = ncomponents; 00220 GaussMix::build(); // TODO Still needed? 00221 resizeDataBeforeTraining(); 00222 // setMixtureTypeGeneral(l, ncomponents, w); // TODO Remove this line when it works. 00223 00224 // storage for neighbors 00225 Mat delta_neighbors(nneighbors, w); 00226 Vec eigvals(ncomponents+1); 00227 Mat components_eigenvecs(ncomponents+1,w); 00228 for(int i=0; i<l; i++) 00229 { 00230 if(i%100==0) 00231 cerr << "[SEQUENTIAL TRAIN: processing pattern #" << i << "/" << l << "]\n"; 00232 00233 // center is sample 00234 Vec center_ = center(i); 00235 00236 if(!learn_mu) 00237 center_ << trainset(i); 00238 00239 if(use_last_eigenval) 00240 computeLocalPrincipalComponents(trainset, i, delta_neighbors, 00241 eigvals, components_eigenvecs, center_, learn_mu); 00242 else 00243 computeLocalPrincipalComponents(trainset, i, delta_neighbors, 00244 eigvals, components_eigenvecs, center_, learn_mu, global_lambda0); 00245 00246 eigvals *= scale_factor; 00247 00248 // cout<<delta_neighbors<<endl; 00249 00250 /* What is this d all about??? (-- Pascal V.) 00251 real d=0; 00252 for(int k=0;k<delta_neighbors.length();k++) 00253 d+=dist(delta_neighbors(k),Vec(D,0.0),2); 00254 d/=delta_neighbors.length(); 00255 */ 00256 00257 // find out eigenvalue (a.k.a lambda0) that will be used for all inputsize_-K directions 00258 real lambda0; 00259 if(use_last_eigenval) 00260 { 00261 // take last (smallest) eigenvalue as a variance in the non-principal directions 00262 // (but if it is 0 because of linear dependencies in the data, take the 00263 // last, i.e. smallest, non-zero eigenvalue). 00264 int last=min(ncomponents,eigvals.length()-1); 00265 lambda0 = eigvals[last]; 00266 while (fast_exact_is_equal(lambda0, 0) && last>0) 00267 lambda0 = eigvals[--last]; 00268 // the sigma-square for all remaining dimensions 00269 if (fast_exact_is_equal(lambda0, 0)) 00270 PLERROR("All (%i) principal components have zero variance!?",eigvals.length()); 00271 } 00272 else 00273 lambda0 = global_lambda0; 00274 00275 alpha[i] = 1.0 / l; 00276 n_eigen = eigvals.length() - 1; 00277 //GaussMix::build(); 00278 //resizeStuffBeforeTraining(); 00279 eigenvalues(i) << eigvals; 00280 eigenvalues(i, n_eigen_computed - 1) = lambda0; 00281 eigenvectors[i] << components_eigenvecs; 00282 // setGaussianGeneral(i, 1.0/l, center, eigvals.subVec(0,eigvals.length()-1), components_eigenvecs.subMatRows(0,eigvals.length()-1), lambda0); 00283 } 00284 stage = 1; 00285 // precomputeStuff(); TODO Put back? 00286 build(); 00287 } 00288 00289 int ManifoldParzen2::find_nearest_neighbor(VMat data, Vec x) const 00290 { 00291 int ret = -1; 00292 real distance = MISSING_VALUE; 00293 real temp; 00294 for(int i=0; i<data->length(); i++) 00295 { 00296 data->getRow(i,row); 00297 temp = dist(row,x,2); 00298 if(is_missing(distance) || temp<distance) 00299 { 00300 distance = temp; 00301 ret = i; 00302 } 00303 } 00304 return ret; 00305 } 00306 00307 real ManifoldParzen2::evaluate(Vec x1,Vec x2,real scale) 00308 { 00309 real ret; 00310 int i = find_nearest_neighbor(train_set,x2); 00311 00312 if(fast_exact_is_equal(scale, 1)) 00313 ret = computeLogLikelihood(x1,i); 00314 else 00315 { 00316 eigenval_copy << eigenvalues(i); 00317 eigenvalues(i) *= scale; 00318 00319 // Maybe sigma_min should be adjusted! 00320 00321 ret = computeLogLikelihood(x1,i); 00322 00323 eigenvalues(i) << eigenval_copy; 00324 } 00325 return exp(ret); 00326 00327 /* 00328 row = x1 - mu(i); 00329 ret = 0.5 * scale * pownorm(row)/eigenvalues(i,n_eigen_computed - 1); 00330 for (int k = 0; k < n_eigen_computed - 1; k++) { 00331 ret += 0.5 * scale * (1.0 / eigenvalues(i,k) - 1.0/eigenvalues(i,n_eigen_computed-1)) * square(dot(eigenvectors[i](k), row)); 00332 } 00333 return ret; 00334 */ 00335 } 00336 00337 real ManifoldParzen2::evaluate_i_j(int i,int j,real scale) 00338 { 00339 real ret; 00340 00341 if(fast_exact_is_equal(scale, 1)) 00342 ret = computeLogLikelihood(center(i),j); 00343 else 00344 { 00345 eigenval_copy << eigenvalues(j); 00346 eigenvalues(j) *= scale; 00347 00348 // Maybe sigma_min should be adjusted! 00349 00350 ret = computeLogLikelihood(center(i),j); 00351 00352 eigenvalues(j) << eigenval_copy; 00353 } 00354 return exp(ret); 00355 00356 /* 00357 row = mu(i) - mu(j); 00358 ret = scale * pownorm(row)/eigenvalues(j,n_eigen_computed - 1); 00359 for (int k = 0; k < n_eigen_computed - 1; k++) { 00360 ret += scale * (1.0 / eigenvalues(j,k) - 1.0/eigenvalues(j,n_eigen_computed-1)) * square(dot(eigenvectors[j](k), row)); 00361 } 00362 return ret; 00363 */ 00364 } 00365 00366 00367 00369 // computeOutput // 00371 void ManifoldParzen2::computeOutput(const Vec& input, Vec& output) const 00372 { 00373 switch(outputs_def[0]) 00374 { 00375 case 'r': 00376 { 00377 int i, last_i=-1; 00378 int nstep = 100000; 00379 real step = 0.001; 00380 int save_every = 100; 00381 string fsave = ""; 00382 string musave = ""; 00383 VMat temp; 00384 t_row.resize(input.length()); 00385 row.resize(input.length()); 00386 t_row << input; 00387 mu_temp.resize(center.length(),center.width()); 00388 temp_eigv.resize(input.length()); 00389 mu_temp << center; 00390 for(int s=0; s<nstep;s++) 00391 { 00392 i = find_nearest_neighbor(new MemoryVMatrix(mu_temp),t_row); 00393 00394 if(s % save_every == 0) 00395 { 00396 fsave = "mp_walk_" + tostring(s) + ".amat"; 00397 temp = new MemoryVMatrix(t_row.toMat(1,t_row.length())); 00398 temp->saveAMAT(fsave,false,true); 00399 00400 musave = "mp_mu_" + tostring(s) + ".amat"; 00401 temp = new MemoryVMatrix(mu_temp(i).toMat(1,mu_temp(i).length())); 00402 temp->saveAMAT(musave,false,true); 00403 00404 } 00405 temp_eigv << eigenvectors[i](0); 00406 real sign = (last_i == -1 || dot(eigenvectors[i](0),eigenvectors[last_i](0)) >= 0 ? 1 : -1); 00407 t_row += step*sign*temp_eigv ; 00408 last_i = i; 00409 } 00410 output << t_row; 00411 break; 00412 } 00413 default: 00414 inherited::computeOutput(input,output); 00415 } 00416 } 00417 00419 // outputsize // 00421 int ManifoldParzen2::outputsize() const 00422 { 00423 switch(outputs_def[0]) 00424 { 00425 case 'r': 00426 return eigenvectors[0].width(); 00427 default: 00428 return inherited::outputsize(); 00429 } 00430 } 00431 00432 } // end of namespace PLearn 00433 00434 00435 /* 00436 Local Variables: 00437 mode:c++ 00438 c-basic-offset:4 00439 c-file-style:"stroustrup" 00440 c-file-offsets:((innamespace . 0)(inline-open . 0)) 00441 indent-tabs-mode:nil 00442 fill-column:79 00443 End: 00444 */ 00445 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :