PLearn 0.1
GeodesicDistanceKernel.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // GeodesicDistanceKernel.cc
00004 //
00005 // Copyright (C) 2004 Olivier Delalleau 
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: GeodesicDistanceKernel.cc 9117 2008-06-11 15:29:02Z tihocan $ 
00037  ******************************************************* */
00038 
00039 // Authors: Olivier Delalleau
00040 
00044 #include "DistanceKernel.h"
00045 #include <plearn/vmat/FileVMatrix.h>
00046 #include "GeodesicDistanceKernel.h"
00047 
00048 namespace PLearn {
00049 using namespace std;
00050 
00052 // GeodesicDistanceKernel //
00054 GeodesicDistanceKernel::GeodesicDistanceKernel() 
00055     : geodesic_file(""),
00056       knn(10),
00057       pow_distance(false),
00058       shortest_algo("floyd")
00059 {
00060     distance_kernel = new DistanceKernel(2);
00061 }
00062 
00063 GeodesicDistanceKernel::GeodesicDistanceKernel(
00064         Ker the_distance_kernel, int the_knn,
00065         const PPath& the_geodesic_file, bool the_pow_distance,
00066         const string& the_method,
00067         bool call_build_):
00068     inherited(true, call_build_),
00069     geodesic_file(the_geodesic_file),
00070     knn(the_knn),
00071     pow_distance(the_pow_distance),
00072     shortest_algo(the_method)
00073 {
00074     distance_kernel = the_distance_kernel;
00075     if (call_build_)
00076         build_();
00077 }
00078 
00079 PLEARN_IMPLEMENT_OBJECT(GeodesicDistanceKernel,
00080                         "Computes the geodesic distance based on k nearest neighbors.",
00081                         ""
00082     );
00083 
00085 // declareOptions //
00087 void GeodesicDistanceKernel::declareOptions(OptionList& ol)
00088 {
00089     // Build options.
00090 
00091     declareOption(ol, "knn", &GeodesicDistanceKernel::knn, OptionBase::buildoption,
00092                   "The number of nearest neighbors considered.");
00093 
00094     declareOption(ol, "distance_kernel", &GeodesicDistanceKernel::distance_kernel, OptionBase::buildoption,
00095                   "The kernel giving the distance between two points.");
00096 
00097     declareOption(ol, "pow_distance", &GeodesicDistanceKernel::pow_distance, OptionBase::buildoption,
00098                   "If set to 1, then it will compute the squared geodesic distance.");
00099 
00100     declareOption(ol, "geodesic_file", &GeodesicDistanceKernel::geodesic_file, OptionBase::buildoption,
00101                   "If provided, the geodesic distances will be saved in this file in binary format.");
00102 
00103     declareOption(ol, "shortest_algo", &GeodesicDistanceKernel::shortest_algo, OptionBase::buildoption,
00104                   "The algorithm used to compute the geodesic distances:\n"
00105                   " - floyd     : Floyd's algorithm\n"
00106                   " - dijkstra  : Dijkstra's algorithm");
00107 
00108     // Learnt options.
00109 
00110     declareOption(ol, "geo_distances", &GeodesicDistanceKernel::geo_distances, OptionBase::learntoption,
00111                   "The geodesic distances between training points.");
00112 
00113     // Now call the parent class' declareOptions
00114     inherited::declareOptions(ol);
00115 }
00116 
00118 // build //
00120 void GeodesicDistanceKernel::build()
00121 {
00122     inherited::build();
00123     build_();
00124 }
00125 
00127 // build_ //
00129 void GeodesicDistanceKernel::build_()
00130 {
00131 }
00132 
00134 // computeNearestGeodesicNeighbour //
00136 int GeodesicDistanceKernel::computeNearestGeodesicNeighbour(int i, const Mat& distances_xi_x_sorted, real* dist_i) const {
00137     real min = distances_xi_x_sorted(0,0) + geo_distances->get(i, int(distances_xi_x_sorted(0,1)));
00138     real dist;
00139     int indice = 0;
00140     for (int j = 1; j < knn; j++) {
00141         dist = distances_xi_x_sorted(j,0) + geo_distances->get(i, int(distances_xi_x_sorted(j,1)));
00142         if (dist < min) {
00143             min = dist;
00144             indice = j;
00145         }
00146     }
00147     if (dist_i)
00148         *dist_i = min;
00149     return int(distances_xi_x_sorted(indice,1));
00150 }
00151 
00153 // computeShortestDistance //
00155 real GeodesicDistanceKernel::computeShortestDistance(int i, const Mat& distances_xi_x_sorted) const {
00156     static real result;
00157     computeNearestGeodesicNeighbour(i, distances_xi_x_sorted, &result);
00158     return result;
00159 }
00160 
00162 // evaluate //
00164 real GeodesicDistanceKernel::evaluate(const Vec& x1, const Vec& x2) const {
00165     distance_kernel->computeNearestNeighbors(x1, dist_xi_x_sorted1, knn);
00166     distance_kernel->computeNearestNeighbors(x2, dist_xi_x_sorted2, knn);
00167     real min = REAL_MAX;
00168     real dist;
00169     for (int j = 0; j < knn; j++) {
00170         for (int k = 0; k < knn; k++) {
00171             dist = dist_xi_x_sorted1(j,0) + dist_xi_x_sorted2(k,0)
00172                 + geo_distances->get(int(dist_xi_x_sorted1(j,1)), int(dist_xi_x_sorted2(k,1)));
00173             if (dist < min) {
00174                 min = dist;
00175             }
00176         }
00177     }
00178     if (pow_distance) {
00179         return square(min);
00180     } else {
00181         return min;
00182     }
00183 }
00184 
00186 // evaluate_i_j //
00188 real GeodesicDistanceKernel::evaluate_i_j(int i, int j) const {
00189     if (pow_distance) {
00190         return square(geo_distances->get(i,j));
00191     } else {
00192         return geo_distances->get(i,j);
00193     }
00194 }
00195 
00197 // evaluate_i_x //
00199 real GeodesicDistanceKernel::evaluate_i_x(int i, const Vec& x, real squared_norm_of_x) const {
00200     return evaluate_i_x_again(i, x, squared_norm_of_x, true);
00201 }
00202 
00204 // evaluate_i_x_from_distances //
00206 real GeodesicDistanceKernel::evaluate_i_x_from_distances(int i, const Mat& distances_xi_x_sorted) const {
00207     if (pow_distance) {
00208         return square(computeShortestDistance(i, distances_xi_x_sorted));
00209     } else {
00210         return computeShortestDistance(i, distances_xi_x_sorted);
00211     }
00212 }
00213 
00215 // evaluate_i_x_again //
00217 real GeodesicDistanceKernel::evaluate_i_x_again(int i, const Vec& x, real squared_norm_of_x, bool first_time) const {
00218     if (first_time) {
00219         distance_kernel->computeNearestNeighbors(x, dist_xi_x_sorted, knn);
00220     }
00221     if (pow_distance) {
00222         return square(computeShortestDistance(i, dist_xi_x_sorted));
00223     } else {
00224         return computeShortestDistance(i, dist_xi_x_sorted);
00225     }
00226 }
00227 
00229 // makeDeepCopyFromShallowCopy //
00231 void GeodesicDistanceKernel::makeDeepCopyFromShallowCopy(CopiesMap& copies)
00232 {
00233     inherited::makeDeepCopyFromShallowCopy(copies);
00234     deepCopyField(dist_xi_x_sorted1, copies);
00235     deepCopyField(dist_xi_x_sorted2, copies);
00236     deepCopyField(dist_xi_x_sorted, copies);
00237     deepCopyField(distance_kernel, copies);
00238     deepCopyField(geo_distances, copies);
00239 }
00240 
00242 // setDataForKernelMatrix //
00244 void GeodesicDistanceKernel::setDataForKernelMatrix(VMat the_data) {
00245     inherited::setDataForKernelMatrix(the_data);
00246     distance_kernel->setDataForKernelMatrix(the_data);
00247     int n = n_examples;
00248     // Check whether we have already compute the geodesic distances.
00249     if (geo_distances && geo_distances->length() == n && geo_distances->width() == n) {
00250         return;
00251     }
00252     // Compute pair distances.
00253     Mat distances(n,n);
00254     distance_kernel->computeGramMatrix(distances);
00255     // Compute knn - nearest neighbors.
00256     TMat<int> neighborhoods =
00257         Kernel::computeKNNeighbourMatrixFromDistanceMatrix(
00258             distances, knn, true, report_progress != 0);
00259     // Compute geodesic distance by Floyd or Dijkstra's algorithm.
00260     Mat geodesic(n,n);
00261     real big_value = REAL_MAX / 3.0; // To make sure no overflow.
00262     PP<ProgressBar> pb;
00263     if (report_progress)
00264         pb = new ProgressBar("Computing geodesic distances", n);
00265     if (shortest_algo == "floyd") {
00266         // First initialize the geodesic distances matrix.
00267         geodesic.fill(big_value);
00268         int neighbor;
00269         real d;
00270         for (int i = 0; i < n; i++) {
00271             geodesic(i,i) = 0;
00272             for (int j = 1; j < knn; j++) {
00273                 neighbor = neighborhoods(i,j);
00274                 d = distances(i, neighbor);
00275                 geodesic(i, neighbor) = d;
00276                 geodesic(neighbor, i) = d;
00277             }
00278         }
00279         // And iterate to find geodesic distances.
00280         real dist;
00281         for (int k = 0; k < n; k++) {
00282             for (int i = 0; i < n; i++) {
00283                 for (int j = 0; j < n; j++) {
00284                     dist = geodesic(i,k) + geodesic(k,j);
00285                     if (geodesic(i,j) > dist) {
00286                         geodesic(i,j) = dist;
00287                     }
00288                 }
00289             }
00290             if (report_progress)
00291                 pb->update(k + 1);
00292         }
00293     } else if (shortest_algo == "dijkstra") {
00294         // First build a symmetric neighborhoods matrix
00295         // (j is a neighbor of i if it was already a neighbor, or if i was a
00296         // neighbor of j).
00297         TVec< TVec<int> > sym_neighborhoods(n);
00298         int neighb, i;
00299         for (i = 0; i < n; i++) {
00300             for (int j = 1; j < knn; j++) {
00301                 neighb = neighborhoods(i, j);
00302                 sym_neighborhoods[i].append(neighb);
00303                 sym_neighborhoods[neighb].append(i);
00304             }
00305         }
00306         Vec d;
00307         TVec<bool> T(n);
00308         int t, min, j, m, k;
00309         real dist;
00310         for (k = 0; k < n; k++) {
00311             d = geodesic(k);
00312             d.fill(big_value);
00313             d[k] = 0;
00314             T.fill(true);
00315             for (i = 0; i < n; i++) {
00316                 min = 0;
00317                 while (!T[min])
00318                     min++;
00319                 for (m = min + 1; m < n; m++) {
00320                     if (T[m] && d[m] < d[min]) {
00321                         min = m;
00322                     }
00323                 }
00324                 for (j = 0; j < sym_neighborhoods[min].length(); j++) {
00325                     t = sym_neighborhoods[min][j];
00326                     if (T[t]) {
00327                         dist = d[min] + distances(min, t);
00328                         if (d[t] > dist) {
00329                             d[t] = dist;
00330                         }
00331                     }
00332                 }
00333                 T[min] = false;
00334             }
00335             if (report_progress)
00336                 pb->update(k+1);
00337         }
00338     } else {
00339         PLERROR("In GeodesicDistanceKernel::setDataForKernelMatrix - Unknown "
00340                 "value for 'shortest_algo': %s",
00341                 shortest_algo.c_str());
00342     }
00343     // Save the result in geo_distances.
00344     if (geodesic_file.isEmpty()) {
00345         geo_distances = VMat(geodesic);
00346     } else {
00347         // Use a FileVMatrix to save on disk.
00348         geo_distances = new FileVMatrix(geodesic_file, n, n);
00349         geo_distances->putMat(0, 0, geodesic);
00350     }
00351 }
00352 
00353 } // end of namespace PLearn
00354 
00355 
00356 /*
00357   Local Variables:
00358   mode:c++
00359   c-basic-offset:4
00360   c-file-style:"stroustrup"
00361   c-file-offsets:((innamespace . 0)(inline-open . 0))
00362   indent-tabs-mode:nil
00363   fill-column:79
00364   End:
00365 */
00366 // 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