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 // In all the code:
  46 //
  47 //  * el[e] is the length of edge e
  48 //  * ea[e] is its starting node
  49 //  * eb[e] is its destination node.
  50
  51 // Adds to all the edge length the length of the shortest (which can
  52 // be negative)
  53 void raise_es(int nb_edges, scalar_t *el) {
  54   scalar_t min_es = el[0];
  55   for(int e = 1; e < nb_edges; e++) {
  56     min_es = min(min_es, el[e]);
  57   }
  58   for(int e = 0; e < nb_edges; e++) {
  59     el[e] -= min_es;
  60   }
  61 }
  62
  63 // Adds to every edge length the differential of the psi function on
  64 // it
  65 void add_dpsi_es(int nb_edges, scalar_t *el, int *ea, int *eb, scalar_t *psi) {
  66   for(int e = 0; e < nb_edges; e++) {
  67     el[e] += psi[ea[e]] - psi[eb[e]];
  68   }
  69 }
  70
  71 // Finds the shortest path in the graph and return in
  72 // result_edge_back, for each vertex the edge to follow back from it
  73 // to reach the source with the shortest path, and in result_dist the
  74 // distance to the source. The edge lengths have to be positive.
  75 void find_shortest(int nb_vertices,
  76                    int nb_edges, scalar_t *el, int *ea, int *eb,
  77                    int source, int sink,
  78                    int *result_edge_back, scalar_t *result_dist) {
  79   for(int v = 0; v < nb_vertices; v++) {
  80     result_dist[v] = FLT_MAX;
  81     result_edge_back[v] = -1;
  82   }
  83
  84   result_dist[source] = 0;
  85
  86 #ifdef DEBUG
  87   for(int e = 0; e < nb_edges; e++) {
  88     if(el[e] < 0) abort();
  89   }
  90 #endif
  91
  92   int nb_changes;
  93   scalar_t d;
  94   do {
  95     nb_changes = 0;
  96     for(int e = 0; e < nb_edges; e++) {
  97       d = result_dist[ea[e]] + el[e];
  98       if(d < result_dist[eb[e]]) {
  99         nb_changes++;
 100         result_dist[eb[e]] = d;
 101         result_edge_back[eb[e]] = e;
 102       }
 103     }
 104   } while(nb_changes > 0);
 105 }
 106
 107 // Iterates find_shortest as long as it finds paths of negative
 108 // lengths. Notes which edges are occupied by the superposition of
 109 // paths in result_edge_occupation
 110 void find_best_paths(int nb_vertices,
 111                      int nb_edges, scalar_t *el, int *ea, int *eb,
 112                      int source, int sink,
 113                      int *result_edge_occupation) {
 114   scalar_t *dist = new scalar_t[nb_vertices];
 115   int *edge_back = new int[nb_vertices];
 116   scalar_t *positive_el = new scalar_t[nb_edges];
 117   scalar_t s;
 118   int v;
 119
 120   for(int e = 0; e < nb_edges; e++) {
 121     positive_el[e] = el[e];
 122     result_edge_occupation[e] = 0;
 123   }
 124
 125   raise_es(nb_edges, positive_el);
 126
 127   do {
 128     find_shortest(nb_vertices, nb_edges, positive_el, ea, eb, source, sink, edge_back, dist);
 129     add_dpsi_es(nb_edges, positive_el, ea, eb, dist);
 130     s = 0.0;
 131
 132     // If the new path reaches the sink, we will backtrack on it to
 133     // compute its score and invert edges
 134
 135     if(edge_back[sink] >= 0) {
 136
 137       v = sink;
 138       while(v != source) {
 139         int e = edge_back[v];
 140         ASSERT(eb[e] == v);
 141         v = ea[e];
 142         s += el[e];
 143       }
 144
 145       // We found a good path (score < 0), we revert the edges along
 146       // the path and invert their occupation (since adding a path on
 147       // a path already occupied means removing flow on it)
 148
 149       if(s < 0) {
 150         v = sink;
 151 #ifdef VERBOSE
 152         cout << "FOUND A PATH OF LENGTH " << s << endl;
 153 #endif
 154         while(v != source) {
 155           int e = edge_back[v];
 156           ASSERT(eb[e] == v);
 157           v = ea[e];
 158 #ifdef VERBOSE
 159           cout << "INVERTING " << ea[e] << " -> " << eb[e] << " [" << el[e] << "]" << endl;
 160 #endif
 161           int k = eb[e];
 162           eb[e] = ea[e];
 163           ea[e] = k;
 164           positive_el[e] = - positive_el[e];
 165           el[e] = - el[e];
 166           result_edge_occupation[e] = 1 - result_edge_occupation[e];
 167         }
 168       }
 169     }
 170   } while(s < 0);
 171
 172   delete[] positive_el;
 173   delete[] dist;
 174   delete[] edge_back;
 175 }
 176
 177 int main(int argc, char **argv) {
 178
 179   if(argc < 2) {
 180     cerr << argv[0] << " <graph file>" << endl;
 181     exit(EXIT_FAILURE);
 182   }
 183
 184   ifstream *file = new ifstream(argv[1]);
 185
 186   int nb_edges, nb_vertices;
 187   int source, sink;
 188
 189   if(file->good()) {
 190
 191     (*file) >> nb_vertices >> nb_edges;
 192     (*file) >> source >> sink;
 193
 194     cout << "INPUT nb_edges " << nb_edges << endl;
 195     cout << "INPUT nb_vertices " << nb_vertices << endl;
 196     cout << "INPUT source " << source << endl;
 197     cout << "INPUT sink " << sink << endl;
 198
 199     scalar_t *el = new scalar_t[nb_edges];
 200     int *ea = new int[nb_edges];
 201     int *eb = new int[nb_edges];
 202     int *edge_occupation = new int[nb_edges];
 203
 204     for(int e = 0; e < nb_edges; e++) {
 205       (*file) >> ea[e] >> eb[e] >> el[e];
 206       cout << "INPUT_EDGE " << ea[e] << " " << eb[e] << " " << el[e] << endl;
 207     }
 208
 209     find_best_paths(nb_vertices, nb_edges, el, ea, eb, source, sink,
 210                     edge_occupation);
 211
 212 #ifdef VERBOSE
 213     // Sanity check on the overall resulting score (the edge lengths
 214     // have been changed, hence should be the opposite of the sum of
 215     // the path lengths)
 216     scalar_t s = 0;
 217     for(int e = 0; e < nb_edges; e++) {
 218       if(edge_occupation[e]) s += el[e];
 219     }
 220     cout << "RESULT_SANITY_CHECK_SCORE " << s << endl;
 221 #endif
 222
 223     for(int e = 0; e < nb_edges; e++) {
 224       if(edge_occupation[e]) {
 225         cout << "RESULT_OCCUPIED_EDGE " << ea[e] << " " << eb[e] << endl;
 226       }
 227     }
 228
 229     delete[] edge_occupation;
 230     delete[] el;
 231     delete[] ea;
 232     delete[] eb;
 233
 234   } else {
 235
 236     cerr << "Can not open " << argv[1] << endl;
 237
 238     delete file;
 239     exit(EXIT_FAILURE);
 240
 241   }
 242
 243   delete file;
 244   exit(EXIT_SUCCESS);
 245 }