PLearn 0.1
ManifoldParzen2.cc
Go to the documentation of this file.
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 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines