PLearn 0.1
NonLocalManifoldParzen.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // NonLocalManifoldParzen.cc
00004 //
00005 // Copyright (C) 2004 Yoshua Bengio & Hugo Larochelle
00006 //
00007 // Redistribution and use in source and binary forms, with or without
00008 // modification, are permitted provided that the following conditions are met:
00009 //
00010 //  1. Redistributions of source code must retain the above copyright
00011 //     notice, this list of conditions and the following disclaimer.
00012 //
00013 //  2. Redistributions in binary form must reproduce the above copyright
00014 //     notice, this list of conditions and the following disclaimer in the
00015 //     documentation and/or other materials provided with the distribution.
00016 //
00017 //  3. The name of the authors may not be used to endorse or promote
00018 //     products derived from this software without specific prior written
00019 //     permission.
00020 //
00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031 //
00032 // This file is part of the PLearn library. For more information on the PLearn
00033 // library, go to the PLearn Web site at www.plearn.org
00034 
00035 /* *******************************************************
00036  * $Id: NonLocalManifoldParzen.cc 8172 2007-10-10 22:41:33Z larocheh $
00037  ******************************************************* */
00038 
00039 // Authors: Yoshua Bengio & Martin Monperrus
00040 
00044 #include "NonLocalManifoldParzen.h"
00045 #include <plearn/display/DisplayUtils.h>
00046 #include <plearn/math/plapack.h>
00047 #include <plearn/var/AffineTransformVariable.h>
00048 #include <plearn/var/AffineTransformWeightPenalty.h>
00049 #include <plearn/var/ColumnSumVariable.h>
00050 #include <plearn/var/NllGeneralGaussianVariable.h>
00051 #include <plearn/var/NoBpropVariable.h>
00052 #include <plearn/var/ReshapeVariable.h>
00053 #include <plearn/var/SourceVariable.h>
00054 #include <plearn/var/SquareVariable.h>
00055 #include <plearn/var/SumOfVariable.h>
00056 #include <plearn/var/TanhVariable.h>
00057 #include <plearn/var/ThresholdBpropVariable.h>
00058 #include <plearn/var/Var_operators.h>
00059 #include <plearn/vmat/AppendNeighborsVMatrix.h>
00060 #include <plearn/vmat/ConcatColumnsVMatrix.h>
00061 
00062 
00063 namespace PLearn {
00064 using namespace std;
00065 
00066 
00067 NonLocalManifoldParzen::NonLocalManifoldParzen()
00068     :  
00069     reference_set(0), 
00070     ncomponents(1), 
00071     nneighbors(5), 
00072     nneighbors_density(-1), 
00073     store_prediction(false),
00074     learn_mu(false),
00075     sigma_init(0.1), 
00076     sigma_min(0.00001), 
00077     mu_nneighbors(2), 
00078     sigma_threshold_factor(-1), 
00079     svd_threshold(1e-8), 
00080     nhidden(10), 
00081     weight_decay(0),
00082     penalty_type("L2_square"),
00083     batch_size(1)
00084 {
00085 }
00086 
00087 PLEARN_IMPLEMENT_OBJECT(NonLocalManifoldParzen, 
00088                         "Non-Local version of Manifold Parzen Windows",
00089                         "Manifold Parzen Windows density model, where the\n"
00090                         "parameters of the kernel for each training point\n"
00091                         "are predicted by a neural network.\n"
00092     );
00093 
00094 
00095 void NonLocalManifoldParzen::declareOptions(OptionList& ol)
00096 {
00097 
00098     declareOption(ol, "parameters", &NonLocalManifoldParzen::parameters, 
00099                   OptionBase::learntoption,
00100                   "Parameters of the tangent_predictor function.\n"
00101         );
00102 
00103     declareOption(ol, "reference_set", &NonLocalManifoldParzen::reference_set, 
00104                   OptionBase::learntoption,
00105                   "Reference points for density computation.\n"
00106         );
00107 
00108     declareOption(ol, "ncomponents", &NonLocalManifoldParzen::ncomponents, 
00109                   OptionBase::buildoption,
00110                   "Number of \"principal components\" to predict\n"
00111                   "for kernel parameters prediction.\n"
00112         );
00113 
00114     declareOption(ol, "nneighbors", &NonLocalManifoldParzen::nneighbors, 
00115                   OptionBase::buildoption,
00116                   "Number of nearest neighbors to consider in training procedure.\n"
00117         );
00118 
00119     declareOption(ol, "nneighbors_density", 
00120                   &NonLocalManifoldParzen::nneighbors_density, 
00121                   OptionBase::buildoption,
00122                   "Number of nearest neighbors to consider for\n"
00123                   "p(x) density estimation.\n"
00124         );
00125 
00126     declareOption(ol, "store_prediction", 
00127                   &NonLocalManifoldParzen::store_prediction, 
00128                   OptionBase::buildoption,
00129                   "Indication that the predicted parameters should be stored.\n"
00130                   "This may make testing faster. Note that the predictions are\n"
00131                   "stored after the last training stage\n"
00132         );
00133 
00134 
00135     declareOption(ol, "paramsvalues", 
00136                   &NonLocalManifoldParzen::paramsvalues, 
00137                   OptionBase::learntoption,
00138                   "The learned parameter vector.\n"
00139         );
00140 
00141     // ** Gaussian kernel options
00142 
00143     declareOption(ol, "learn_mu", &NonLocalManifoldParzen::learn_mu, 
00144                   OptionBase::buildoption,
00145                   "Indication that the deviation from the training point\n"
00146                   "in a Gaussian kernel (called mu) should be learned.\n"
00147         );
00148 
00149     declareOption(ol, "sigma_init", &NonLocalManifoldParzen::sigma_init, 
00150                   OptionBase::buildoption,
00151                   "Initial minimum value for sigma noise.\n"
00152         );
00153 
00154     declareOption(ol, "sigma_min", &NonLocalManifoldParzen::sigma_min, 
00155                   OptionBase::buildoption,
00156                   "The minimum value for sigma noise.\n"
00157         );
00158 
00159     declareOption(ol, "mu_nneighbors", &NonLocalManifoldParzen::mu_nneighbors, 
00160                   OptionBase::buildoption,
00161                   "Number of nearest neighbors to learn the mus \n"
00162                   "(if < 0, mu_nneighbors = nneighbors).\n"
00163         );
00164 
00165     declareOption(ol, "sigma_threshold_factor", 
00166                   &NonLocalManifoldParzen::sigma_threshold_factor, 
00167                   OptionBase::buildoption,
00168                   "Threshold factor of the gradient on the sigma noise\n"
00169                   "parameter of the Gaussian kernel. If < 0, then\n"
00170                   "no threshold is used."
00171         );
00172 
00173     declareOption(ol, "svd_threshold", 
00174                   &NonLocalManifoldParzen::svd_threshold, OptionBase::buildoption,
00175                   "Threshold to accept singular values of F in solving for\n"
00176                   "linear combination weights on tangent subspace.\n"
00177         );
00178 
00179     // ** Neural network predictor **
00180 
00181     declareOption(ol, "nhidden", 
00182                   &NonLocalManifoldParzen::nhidden, OptionBase::buildoption,
00183                   "Number of hidden units of the neural network.\n"
00184         );
00185 
00186     declareOption(ol, "weight_decay", &NonLocalManifoldParzen::weight_decay, 
00187                   OptionBase::buildoption,
00188                   "Global weight decay for all layers.\n");
00189 
00190     declareOption(ol, "penalty_type", &NonLocalManifoldParzen::penalty_type,
00191                   OptionBase::buildoption,
00192                   "Penalty to use on the weights (for weight and bias decay).\n"
00193                   "Can be any of:\n"
00194                   "  - \"L1\": L1 norm,\n"
00195                   "  - \"L2_square\" (default): square of the L2 norm.\n");
00196 
00197     declareOption(ol, "optimizer", &NonLocalManifoldParzen::optimizer, 
00198                   OptionBase::buildoption,
00199                   "Optimizer that optimizes the cost function.\n"
00200         );
00201 
00202     declareOption(ol, "batch_size", 
00203                   &NonLocalManifoldParzen::batch_size, OptionBase::buildoption,
00204                   "How many samples to use to estimate the average gradient\n"
00205                   "before updating the weights. If <= 0, is equivalent to\n"
00206                   "specifying training_set->length() \n");
00207 
00208 
00209     // ** Stored outputs of neural network
00210 
00211     declareOption(ol, "mus", 
00212                   &NonLocalManifoldParzen::mus, OptionBase::learntoption,
00213                   "The stored mu vectors for the reference set.\n"
00214         );
00215 
00216     declareOption(ol, "sns", &NonLocalManifoldParzen::sns, 
00217                   OptionBase::learntoption,
00218                   "The stored sigma noise values for the reference set.\n"
00219         );
00220 
00221     declareOption(ol, "sms", &NonLocalManifoldParzen::sms, 
00222                   OptionBase::learntoption,
00223                   "The stored sigma manifold values for the reference set.\n"
00224         );
00225 
00226     declareOption(ol, "Fs", &NonLocalManifoldParzen::Fs, OptionBase::learntoption,
00227                   "The storaged \"principal components\" (F) values for\n"
00228                   "the reference set.\n"
00229         );
00230 
00231 
00232     // Now call the parent class' declareOptions
00233     inherited::declareOptions(ol);
00234 }
00235 
00236 void NonLocalManifoldParzen::build_()
00237 {
00238 
00239     if (inputsize_>0)
00240     {
00241         if (nhidden <= 0) 
00242             PLERROR("NonLocalManifoldParzen::Number of hidden units "
00243                     "should be positive, now %d\n",nhidden);
00244 
00245         Var log_n_examples(1,1,"log(n_examples)");
00246         if(train_set)
00247         {
00248             L = train_set->length();
00249             reference_set = train_set; 
00250         }
00251 
00252         log_L= pl_log((real) L);
00253         parameters.resize(0);
00254         
00255         // Neural network prediction of principal components
00256 
00257         x = Var(inputsize_);
00258         x->setName("x");
00259 
00260         W = Var(nhidden+1,inputsize_,"W");
00261         parameters.append(W);
00262 
00263         Var a; // outputs of hidden layer
00264         a = affine_transform(x,W);
00265         a->setName("a");
00266 
00267         V = Var(ncomponents*(inputsize_+1),nhidden,"V");
00268         parameters.append(V);
00269 
00270         // TODO: instead, make NllGeneralGaussianVariable use vector... (DONE)
00271         //components = reshape(affine_transform(V,a),ncomponents,n);
00272         components = affine_transform(V,a);
00273         components->setName("components");
00274 
00275         // Gaussian kernel parameters prediction
00276 
00277         muV = Var(inputsize_+1,nhidden,"muV");
00278         snV = Var(2,nhidden,"snV");
00279     
00280         parameters.append(muV);
00281         parameters.append(snV);
00282 
00283         if(learn_mu)
00284             mu = affine_transform(muV,a);
00285         else
00286         {
00287             mu = new SourceVariable(inputsize_,1);
00288             mu->value.clear();
00289         }
00290         mu->setName("mu");
00291 
00292         min_sig = new SourceVariable(1,1);
00293         min_sig->value[0] = sigma_min;
00294         min_sig->setName("min_sig");
00295         init_sig = Var(1,1);
00296         init_sig->setName("init_sig");
00297         parameters.append(init_sig);
00298 
00299         sn = square(affine_transform(snV,a)) + min_sig + square(init_sig);
00300         sn->setName("sn");
00301         
00302         if(sigma_threshold_factor > 0)
00303             sn = threshold_bprop(sn,sigma_threshold_factor);
00304 
00305         predictor = Func(x, parameters , components & mu & sn );
00306     
00307         Var target_index = Var(1,1);
00308         target_index->setName("target_index");
00309         Var neighbor_indexes = Var(nneighbors,1);
00310         neighbor_indexes->setName("neighbor_indexes");
00311 
00312         tangent_targets = Var(nneighbors,inputsize_);
00313         if(mu_nneighbors < 0 ) mu_nneighbors = nneighbors;
00314 
00315         Var nll;
00316         nll = nll_general_gaussian(components, mu, sn, tangent_targets, 
00317                                    log_L, learn_mu, mu_nneighbors); 
00318 
00319         Var knn = new SourceVariable(1,1);
00320         knn->setName("knn");
00321         knn->value[0] = nneighbors;
00322         sum_nll = new ColumnSumVariable(nll) / knn;
00323 
00324         // Weight decay penalty
00325         if(weight_decay > 0 )
00326         {
00327             sum_nll += affine_transform_weight_penalty(
00328                 W,weight_decay,0,penalty_type) + 
00329                 affine_transform_weight_penalty(
00330                 V,weight_decay,0,penalty_type) + 
00331                 affine_transform_weight_penalty(
00332                 muV,weight_decay,0,penalty_type) + 
00333                 affine_transform_weight_penalty(
00334                 snV,weight_decay,0,penalty_type);
00335         }
00336 
00337         cost_of_one_example = Func(x & tangent_targets & target_index & 
00338                                    neighbor_indexes, parameters, sum_nll);
00339 
00340         if(nneighbors_density >= L || nneighbors_density < 0) 
00341             nneighbors_density = L;
00342 
00343         // Output storage variables
00344         t_row.resize(inputsize_);
00345         Ut_svd.resize(inputsize_,inputsize_);
00346         V_svd.resize(ncomponents,ncomponents);
00347         F.resize(components->length(),components->width());
00348         z.resize(inputsize_);
00349         x_minus_neighbor.resize(inputsize_);
00350         neighbor_row.resize(inputsize_);
00351 
00352         // log_density and Kernel methods variables
00353         U_temp.resize(ncomponents,inputsize_);
00354         mu_temp.resize(inputsize_);
00355         sm_temp.resize(ncomponents);
00356         sn_temp.resize(1);
00357         diff.resize(inputsize_);
00358 
00359         mus.resize(L, inputsize_);
00360         sns.resize(L);
00361         sms.resize(L,ncomponents);
00362         Fs.resize(L);
00363         for(int i=0; i<L; i++)
00364         {
00365             Fs[i].resize(ncomponents,inputsize_);
00366         }
00367 
00368         if(paramsvalues.length() == parameters.nelems())
00369             parameters << paramsvalues;
00370         else
00371         {
00372             paramsvalues.resize(parameters.nelems());
00373             initializeParams();
00374             if(optimizer)
00375                 optimizer->reset();
00376         }
00377         parameters.makeSharedValue(paramsvalues);
00378     }
00379 
00380 }
00381 
00382 void NonLocalManifoldParzen::knn(const VMat& vm, const Vec& x, const int& k, TVec<int>& neighbors, bool sortk) const
00383 {
00384     int n = vm->length();
00385     distances.resize(n,2);
00386     distances.column(1) << Vec(0, n-1, 1);
00387     dk.setDataForKernelMatrix(vm);
00388     t_dist.resize(n);
00389     dk.evaluate_all_i_x(x, t_dist);
00390     distances.column(0) << t_dist;
00391     partialSortRows(distances, k, sortk);
00392     neighbors.resize(k);
00393 
00394     for (int i=0; i < k  && i<n; i++)
00395     {
00396         neighbors[i] = int(distances(i,1));
00397     }
00398 
00399 }
00400 
00401 // ### Nothing to add here, simply calls build_
00402 void NonLocalManifoldParzen::build()
00403 {
00404     inherited::build();
00405     build_();
00406 }
00407 
00408 #ifdef __INTEL_COMPILER
00409 #pragma warning(disable:1419)  // Get rid of compiler warning.
00410 #endif
00411 extern void varDeepCopyField(Var& field, CopiesMap& copies);
00412 #ifdef __INTEL_COMPILER
00413 #pragma warning(default:1419)
00414 #endif
00415 
00416 void NonLocalManifoldParzen::makeDeepCopyFromShallowCopy(CopiesMap& copies)
00417 {  
00418     inherited::makeDeepCopyFromShallowCopy(copies);
00419 
00420     // Protected
00421 
00422     deepCopyField(cost_of_one_example, copies);
00423     varDeepCopyField(x, copies);
00424     varDeepCopyField(W, copies);
00425     varDeepCopyField(V, copies);
00426     varDeepCopyField(muV, copies);
00427     varDeepCopyField(snV, copies);
00428     varDeepCopyField(tangent_targets, copies);
00429     varDeepCopyField(components, copies);
00430     varDeepCopyField(mu, copies);
00431     varDeepCopyField(sn, copies);
00432     varDeepCopyField(sum_nll, copies);
00433     varDeepCopyField(min_sig, copies);
00434     varDeepCopyField(init_sig, copies);
00435     deepCopyField(predictor, copies);
00436     deepCopyField(U_temp,copies);
00437     deepCopyField(F, copies);
00438     deepCopyField(distances,copies);
00439     deepCopyField(mu_temp,copies);
00440     deepCopyField(sm_temp,copies);
00441     deepCopyField(sn_temp,copies);
00442     deepCopyField(diff,copies);
00443     deepCopyField(z,copies);
00444     deepCopyField(x_minus_neighbor,copies);
00445     deepCopyField(t_row,copies);
00446     deepCopyField(neighbor_row,copies);
00447     deepCopyField(log_gauss,copies);
00448     deepCopyField(t_dist,copies);
00449     deepCopyField(t_nn,copies);
00450     deepCopyField(Ut_svd, copies);
00451     deepCopyField(V_svd, copies);
00452     deepCopyField(S_svd, copies);
00453     deepCopyField(mus, copies);
00454     deepCopyField(sns, copies);
00455     deepCopyField(sms, copies);
00456     deepCopyField(Fs, copies);
00457     deepCopyField(train_set_with_targets, copies);
00458     deepCopyField(targets_vmat, copies);
00459     varDeepCopyField(totalcost, copies);
00460     deepCopyField(paramsvalues, copies);
00461     
00462     // Public
00463 
00464     deepCopyField(parameters, copies);    
00465     deepCopyField(reference_set,copies);
00466     deepCopyField(optimizer, copies);
00467 
00468 }
00469 
00470 
00471 void NonLocalManifoldParzen::forget()
00472 {
00473     inherited::forget();
00474     if (train_set) initializeParams();
00475     if(optimizer) optimizer->reset();
00476     stage = 0;
00477 }
00478 
00479 void NonLocalManifoldParzen::train()
00480 {
00481     // Check whether gradient descent is going to be done
00482     // If not, then we don't need to store the parameters,
00483     // except for sn...
00484     bool flag = (nstages == stage);
00485 
00486     // Update sigma_min, in case it was changed,
00487     // e.g. using an HyperLearner
00488     min_sig->value[0] = sigma_min;
00489 
00490     // Set train_stats if not already done.
00491     if (!train_stats)
00492         train_stats = new VecStatsCollector();
00493 
00494     if (!cost_of_one_example)
00495         PLERROR("NonLocalManifoldParzen::train: build has not been run after setTrainingSet!");
00496 
00497     if(stage == 0)
00498     {
00499         targets_vmat = append_neighbors(
00500             train_set, nneighbors, true);
00501         nsamples = batch_size>0 ? batch_size : train_set->length();
00502 
00503         totalcost = meanOf(train_set_with_targets, cost_of_one_example, nsamples);
00504 
00505         if(optimizer)
00506         {
00507             optimizer->setToOptimize(parameters, totalcost);
00508             optimizer->build();
00509         }
00510         else PLERROR("NonLocalManifoldParzen::train can't train without setting an optimizer first!");
00511     }
00512 
00513     int optstage_per_lstage = train_set->length()/nsamples;
00514 
00515     PP<ProgressBar> pb;
00516     if(report_progress>0)
00517         pb = new ProgressBar("Training NonLocalManifoldParzen from stage " + tostring(stage) + " to " + tostring(nstages), nstages-stage);
00518 
00519     t_row.resize(train_set.width());
00520 
00521     int initial_stage = stage;
00522     bool early_stop=false;
00523     while(stage<nstages && !early_stop)
00524     {
00525         optimizer->nstages = optstage_per_lstage;
00526         train_stats->forget();
00527         optimizer->early_stop = false;
00528         optimizer->optimizeN(*train_stats);
00529         train_stats->finalize();
00530         if(verbosity>2)
00531             cout << "Epoch " << stage << " train objective: " << train_stats->getMean() << endl;
00532         ++stage;
00533         if(pb)
00534             pb->update(stage-initial_stage);
00535     }
00536     if(verbosity>1)
00537         cout << "EPOCH " << stage << " train objective: " << train_stats->getMean() << endl;
00538 
00539     if(store_prediction && !flag)
00540     {
00541         for(int t=0; t<L;t++)
00542         {
00543             reference_set->getRow(t,neighbor_row);
00544             predictor->fprop(neighbor_row, F.toVec() & mus(t) & sns.subVec(t,1));
00545             // N.B. this is the SVD of F'
00546             lapackSVD(F, Ut_svd, S_svd, V_svd,'A',1.5);
00547             for (int k=0;k<ncomponents;k++)
00548             {
00549                 sms(t,k) = mypow(S_svd[k],2);
00550                 Fs[t](k) << Ut_svd(k);
00551             }
00552             sns[t] += sigma_min - min_sig->value[0];
00553         }
00554     }
00555 }
00556 
00558 // initializeParams //
00560 void NonLocalManifoldParzen::initializeParams()
00561 {
00562     real delta = 1.0 / sqrt(real(inputsize_));
00563     random_gen->fill_random_uniform(W->value, -delta, delta);
00564     delta = 1.0 / real(nhidden);
00565     random_gen->fill_random_uniform(V->matValue, -delta, delta);
00566     random_gen->fill_random_uniform(snV->matValue, -delta, delta);
00567     random_gen->fill_random_uniform(muV->matValue, -delta, delta);
00568     W->matValue(0).clear();
00569     V->matValue(0).clear();
00570     muV->matValue(0).clear();
00571     snV->matValue(0).clear();
00572     init_sig->value[0] = sqrt(sigma_init);
00573 }
00574 
00576 // log_density //
00578 real NonLocalManifoldParzen::log_density(const Vec& x) const {
00579     // Compute log-density.
00580     real ret = 0;
00581     t_row << x;
00582     real mahal = 0;
00583     real norm_term = 0;
00584 
00585     // Update sigma_min, in case it was changed,
00586     // e.g. using an HyperLearner
00587 
00588     if(store_prediction && min_sig->value[0] != sigma_min)
00589     {
00590         for(int i=0; i<L; i++)
00591         {
00592             sns[i] += sigma_min - min_sig->value[0];
00593         }
00594     }
00595 
00596     min_sig->value[0] = sigma_min;
00597 
00598     if(nneighbors_density != L)
00599     {
00600         // Fetching nearest neighbors for density estimation.
00601         knn(reference_set,x,nneighbors_density,t_nn,0);
00602         log_gauss.resize(t_nn.length());
00603         for(int neighbor=0; neighbor<t_nn.length(); neighbor++)
00604         {
00605             reference_set->getRow(t_nn[neighbor],neighbor_row);
00606             if(!store_prediction)
00607             {
00608                 predictor->fprop(neighbor_row, F.toVec() & mu_temp & sn_temp);
00609                 // N.B. this is the SVD of F'
00610                 lapackSVD(F, Ut_svd, S_svd, V_svd,'A',1.5);
00611                 for (int k=0;k<ncomponents;k++)
00612                 {
00613                     sm_temp[k] = mypow(S_svd[k],2);
00614                     U_temp(k) << Ut_svd(k);
00615                 }
00616             }
00617             else
00618             {
00619                 if(learn_mu)
00620                     mu_temp << mus(t_nn[neighbor]);
00621                 sn_temp[0] = sns[t_nn[neighbor]];
00622                 sm_temp << sms(t_nn[neighbor]);
00623                 U_temp << Fs[t_nn[neighbor]];
00624             }
00625             if(learn_mu)
00626             {
00627                 substract(t_row,neighbor_row,x_minus_neighbor);
00628                 substract(x_minus_neighbor,mu_temp,z);
00629             }
00630             else
00631                 substract(t_row,neighbor_row,z);
00632                 
00633             mahal = -0.5*pownorm(z)/sn_temp[0];
00634             norm_term = - inputsize_/2.0 * Log2Pi 
00635                 - log_L - 0.5*(inputsize_-ncomponents)*pl_log(sn_temp[0]);
00636 
00637 
00638             for(int k=0; k<ncomponents; k++)
00639             {
00640                 mahal -= square(dot(z,U_temp(k)))*(0.5/(sm_temp[k]+sn_temp[0]) 
00641                                                    - 0.5/sn_temp[0]);
00642                 norm_term -= 0.5*pl_log(sm_temp[k]+sn_temp[0]);
00643             }
00644 
00645             log_gauss[neighbor] = mahal + norm_term;
00646         }
00647     }
00648     else
00649     {
00650         // Fetching nearest neighbors for density estimation.
00651         log_gauss.resize(L);
00652         for(int t=0; t<L;t++)
00653         {
00654             reference_set->getRow(t,neighbor_row);
00655             if(!store_prediction)
00656             {
00657                 predictor->fprop(neighbor_row, F.toVec() & mu_temp & sn_temp);
00658 
00659                 // N.B. this is the SVD of F'
00660                 lapackSVD(F, Ut_svd, S_svd, V_svd,'A',1.5);
00661                 for (int k=0;k<ncomponents;k++)
00662                 {
00663                     sm_temp[k] = mypow(S_svd[k],2);
00664                     U_temp(k) << Ut_svd(k);
00665                 }
00666             }
00667             else
00668             {
00669                 if(learn_mu)
00670                     mu_temp << mus(t);
00671                 sn_temp[0] = sns[t];
00672                 sm_temp << sms(t);
00673                 U_temp << Fs[t];
00674             }
00675 
00676             if(learn_mu)
00677             {
00678                 substract(t_row,neighbor_row,x_minus_neighbor);
00679                 substract(x_minus_neighbor,mu_temp,z);
00680             }
00681             else
00682                 substract(t_row,neighbor_row,z);
00683 
00684             mahal = -0.5*pownorm(z)/sn_temp[0];
00685             norm_term = - inputsize_/2.0 * Log2Pi - log_L 
00686                 - 0.5*(inputsize_-ncomponents)*pl_log(sn_temp[0]);
00687 
00688             for(int k=0; k<ncomponents; k++)
00689             {
00690                 mahal -= square(dot(z,U_temp(k)))*(0.5/(sm_temp[k]+sn_temp[0]) 
00691                                                    - 0.5/sn_temp[0]);
00692                 norm_term -= 0.5*pl_log(sm_temp[k]+sn_temp[0]);
00693             }
00694 
00695             log_gauss[t] = mahal + norm_term;
00696         }
00697     }
00698     ret = logadd(log_gauss);
00699 
00700     return ret;
00701 }
00702 
00703 
00705 // computeOutput //
00707 void NonLocalManifoldParzen::computeOutput(const Vec& input, Vec& output) const
00708 {
00709     switch(outputs_def[0])
00710     {
00711         /*
00712     case 'r':
00713     {
00714         string fsave = "";
00715         VMat temp;
00716         real step_size = rw_size_step;
00717         real dp;
00718         t_row << input;
00719         Vec last_F(inputsize());
00720         for(int s=0; s<rw_n_step;s++)
00721         {
00722             if(s == 0)
00723             {
00724                 predictor->fprop(t_row, F.toVec() & mu_temp & sn_temp);
00725                 last_F << F(rw_ith_component);
00726             }
00727             predictor->fprop(t_row, F.toVec() & mu_temp & sn_temp);
00728 
00729             // N.B. this is the SVD of F'
00730             lapackSVD(F, Ut_svd, S_svd, V_svd,'A',1.5);
00731             F(rw_ith_component) << Ut_svd(rw_ith_component);
00732 
00733             if(s % rw_save_every == 0)
00734             {
00735                 fsave = rw_file_name + tostring(s) + ".amat";
00736                 temp = new MemoryVMatrix(t_row.toMat(1,t_row.length()));
00737                 temp->saveAMAT(fsave,false,true);
00738                 //PLearn::save(fsave,t_row);
00739             }
00740             dp = dot(last_F,F(rw_ith_component));
00741             if(dp>0) dp = 1;
00742             else dp = -1;
00743             t_row += step_size*F(rw_ith_component)*abs(S_svd[rw_ith_component])*dp;
00744             last_F << dp*F(rw_ith_component);
00745         }
00746         output << t_row;
00747 
00748         t_row << input;
00749         for(int s=0; s<rw_n_step;s++)
00750         {
00751             if(s == 0)
00752             {
00753                 predictor->fprop(t_row, F.toVec() & mu_temp & sn_temp);
00754                 last_F << (-1.0)*F(rw_ith_component);
00755             }
00756 
00757 
00758             predictor->fprop(t_row, F.toVec() & mu_temp & sn_temp);
00759 
00760             // N.B. this is the SVD of F'
00761             lapackSVD(F, Ut_svd, S_svd, V_svd,'A',1.5);
00762             F(rw_ith_component) << Ut_svd(rw_ith_component);
00763 
00764             if(s % rw_save_every == 0)
00765             {
00766                 fsave = rw_file_name + tostring(-s) + ".amat";
00767                 temp = new MemoryVMatrix(t_row.toMat(1,t_row.length()));
00768                 temp->saveAMAT(fsave,false,true);
00769                 //PLearn::save(fsave,t_row);
00770             }
00771             dp = dot(last_F,F(rw_ith_component));
00772             if(dp>0) dp = 1;
00773             else dp = -1;
00774             t_row += step_size*F(rw_ith_component)*abs(S_svd[rw_ith_component])*dp;
00775             last_F << dp*F(rw_ith_component);
00776         }
00777         break;
00778     }
00779     case 't':
00780     {
00781         predictor->fprop(input, F.toVec() & mu_temp & sn_temp);
00782         output << F.toVec();
00783         break;
00784     }
00785         */
00786     default:
00787         
00788         inherited::computeOutput(input,output);
00789     }
00790 }
00791 
00793 // outputsize //
00795 int NonLocalManifoldParzen::outputsize() const
00796 {
00797     switch(outputs_def[0])
00798     {
00799         /*
00800     case 'm':
00801         return ncomponents;
00802         break;
00803     case 'r':
00804         return n;
00805     case 't':
00806         return ncomponents*n;
00807         */
00808     default:
00809         return inherited::outputsize();
00810     }
00811 }
00812 
00813 } // end of namespace PLearn
00814 
00815 
00816 /*
00817   Local Variables:
00818   mode:c++
00819   c-basic-offset:4
00820   c-file-style:"stroustrup"
00821   c-file-offsets:((innamespace . 0)(inline-open . 0))
00822   indent-tabs-mode:nil
00823   fill-column:79
00824   End:
00825 */
00826 // 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