PLearn 0.1
|
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 :