miniksp.cc

   1
   2 ///////////////////////////////////////////////////////////////////////////
   3 // START_IP_HEADER                                                       //
   4 //                                                                       //
   5 // This program is free software: you can redistribute it and/or modify  //
   6 // it under the terms of the version 3 of the GNU General Public License //
   7 // as published by the Free Software Foundation.                         //
   8 //                                                                       //
   9 // This program is distributed in the hope that it will be useful, but   //
  10 // WITHOUT ANY WARRANTY; without even the implied warranty of            //
  11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU      //
  12 // General Public License for more details.                              //
  13 //                                                                       //
  14 // You should have received a copy of the GNU General Public License     //
  15 // along with this program. If not, see <http://www.gnu.org/licenses/>.  //
  16 //                                                                       //
  17 // Written by and Copyright (C) Francois Fleuret                         //
  18 // Contact <francois.fleuret@idiap.ch> for comments & bug reports        //
  19 //                                                                       //
  20 // END_IP_HEADER                                                         //
  21 ///////////////////////////////////////////////////////////////////////////
  22
  23 // #define VERBOSE
  24
  25 #include <iostream>
  26 #include <fstream>
  27 #include <cmath>
  28 #include <stdio.h>
  29 #include <stdlib.h>
  30 #include <float.h>
  31
  32 using namespace std;
  33
  34 typedef float scalar_t;
  35
  36 #ifdef DEBUG
  37 #define ASSERT(x) if(!(x)) { \
  38   std::cerr << "ASSERT FAILED IN " << __FILE__ << ":" << __LINE__ << endl; \
  39   abort(); \
  40 }
  41 #else
  42 #define ASSERT(x)
  43 #endif
  44
  45 // Adds to all the edge length the length of the shortest (which can
  46 // be negative)
  47 void raise_es(int nb_edges, scalar_t *el) {
  48   scalar_t min_es = el[0];
  49   for(int e = 1; e < nb_edges; e++) {
  50     min_es = min(min_es, el[e]);
  51   }
  52   for(int e = 0; e < nb_edges; e++) {
  53     el[e] -= min_es;
  54   }
  55 }
  56
  57 // Add to every edge length the differential of the psi function on it
  58 void add_dpsi_es(int nb_edges, scalar_t *el, int *ea, int *eb, scalar_t *psi) {
  59   for(int e = 0; e < nb_edges; e++) {
  60     el[e] += psi[ea[e]] - psi[eb[e]];
  61   }
  62 }
  63
  64 // Find the shortest path in the graph and return in return_edge_back,
  65 // for each vertex the edge to follow back from it to reach the source
  66 // with the shortest path, and in return_dist the distance to the
  67 // source. The edge lengths have to be positive.
  68 void find_shortest(int nb_vertices,
  69                    int nb_edges, scalar_t *el, int *ea, int *eb,
  70                    int source, int sink,
  71                    int *return_edge_back, scalar_t *return_dist) {
  72   for(int v = 0; v < nb_vertices; v++) {
  73     dist[v] = FLT_MAX;
  74     edge_back[v] = -1;
  75   }
  76
  77   dist[source] = 0;
  78
  79 #ifdef DEBUG
  80   for(int e = 0; e < nb_edges; e++) {
  81     if(el[e] < 0) abort();
  82   }
  83 #endif
  84
  85   int nb_changes;
  86   scalar_t d;
  87   do {
  88     nb_changes = 0;
  89     for(int e = 0; e < nb_edges; e++) {
  90       d = dist[ea[e]] + el[e];
  91       if(d < dist[eb[e]]) {
  92         nb_changes++;
  93         dist[eb[e]] = d;
  94         edge_back[eb[e]] = e;
  95       }
  96     }
  97   } while(nb_changes > 0);
  98 }
  99
 100 void find_best_paths(int nb_vertices,
 101                      int nb_edges, scalar_t *el, int *ea, int *eb,
 102                      int source, int sink,
 103                      int *edge_occupation) {
 104   scalar_t *dist = new scalar_t[nb_vertices];
 105   int *edge_back = new int[nb_vertices];
 106   scalar_t *positive_el = new scalar_t[nb_edges];
 107   scalar_t s;
 108   int v;
 109
 110   for(int e = 0; e < nb_edges; e++) {
 111     positive_el[e] = el[e];
 112     edge_occupation[e] = 0;
 113   }
 114
 115   raise_es(nb_edges, positive_el);
 116
 117   do {
 118     find_shortest(nb_vertices, nb_edges, positive_el, ea, eb, source, sink, edge_back, dist);
 119     add_dpsi_es(nb_edges, positive_el, ea, eb, dist);
 120     s = 0.0;
 121
 122     // If the new path reaches the sink, we will backtrack on it to
 123     // compute its score and invert edges
 124
 125     if(edge_back[sink] >= 0) {
 126
 127       v = sink;
 128       while(v != source) {
 129         int e = edge_back[v];
 130         ASSERT(eb[e] == v);
 131         v = ea[e];
 132         s += el[e];
 133       }
 134
 135       // We found a good path (score < 0), we revert the edges along
 136       // the path and note that they are occupied
 137
 138       if(s < 0) {
 139         v = sink;
 140 #ifdef VERBOSE
 141         cout << "FOUND A PATH OF LENGTH " << s << endl;
 142 #endif
 143         while(v != source) {
 144           int e = edge_back[v];
 145           ASSERT(eb[e] == v);
 146           v = ea[e];
 147 #ifdef VERBOSE
 148           cout << "INVERTING " << ea[e] << " -> " << eb[e] << " [" << el[e] << "]" << endl;
 149 #endif
 150           int k = eb[e];
 151           eb[e] = ea[e];
 152           ea[e] = k;
 153           positive_el[e] = - positive_el[e];
 154           el[e] = - el[e];
 155           edge_occupation[e] = 1 - edge_occupation[e];
 156         }
 157       }
 158     }
 159   } while(s < 0);
 160
 161   delete[] positive_el;
 162   delete[] dist;
 163   delete[] edge_back;
 164 }
 165
 166 int main(int argc, char **argv) {
 167
 168   if(argc < 2) {
 169     cerr << argv[0] << " <graph file>" << endl;
 170     exit(EXIT_FAILURE);
 171   }
 172
 173   ifstream *file = new ifstream(argv[1]);
 174
 175   int nb_edges, nb_vertices;
 176   int source, sink;
 177
 178   if(file->good()) {
 179
 180     (*file) >> nb_vertices >> nb_edges;
 181     (*file) >> source >> sink;
 182
 183     cout << "INPUT nb_edges " << nb_edges << endl;
 184     cout << "INPUT nb_vertices " << nb_vertices << endl;
 185     cout << "INPUT source " << source << endl;
 186     cout << "INPUT sink " << sink << endl;
 187
 188     scalar_t *el = new scalar_t[nb_edges];
 189     int *ea = new int[nb_edges];
 190     int *eb = new int[nb_edges];
 191     int *edge_occupation = new int[nb_edges];
 192
 193     for(int e = 0; e < nb_edges; e++) {
 194       (*file) >> ea[e] >> eb[e] >> el[e];
 195       cout << "INPUT_EDGE " << ea[e] << " " << eb[e] << " " << el[e] << endl;
 196     }
 197
 198     find_best_paths(nb_vertices, nb_edges, el, ea, eb, source, sink,
 199                     edge_occupation);
 200
 201     for(int e = 0; e < nb_edges; e++) {
 202       if(edge_occupation[e]) {
 203         cout << "RESULT_OCCUPIED_EDGE " << ea[e] << " " << eb[e] << endl;
 204       }
 205     }
 206
 207     delete[] edge_occupation;
 208     delete[] el;
 209     delete[] ea;
 210     delete[] eb;
 211
 212   } else {
 213
 214     cerr << "Can not open " << argv[1] << endl;
 215
 216     delete file;
 217     exit(EXIT_FAILURE);
 218
 219   }
 220
 221   delete file;
 222   exit(EXIT_SUCCESS);
 223 }