X-Git-Url: https://fleuret.org/cgi-bin/gitweb/gitweb.cgi?a=blobdiff_plain;f=miniksp.cc;h=92e9de54a50adac278b322621a027caab6ae6dba;hb=3d87c6606fb8750621ace760032473d6a7a5d82c;hp=f562423c0aea3844b2fb8af2410023f562f6b37f;hpb=65a716e6dbb3a6b8f061033a2ee3a95ba7bb7a5c;p=mtp.git

diff --git a/miniksp.cc b/miniksp.cc
index f562423..92e9de5 100644
--- a/miniksp.cc
+++ b/miniksp.cc
@@ -20,7 +20,7 @@
 // END_IP_HEADER                                                         //
 ///////////////////////////////////////////////////////////////////////////
 
-#define VERBOSE
+// #define VERBOSE
 
 #include <iostream>
 #include <fstream>
@@ -42,6 +42,14 @@ typedef float scalar_t;
 #define ASSERT(x)
 #endif
 
+// In all the code:
+//
+//  * el[e] is the length of edge e
+//  * ea[e] is its starting node
+//  * eb[e] is its destination node.
+
+// Adds to all the edge length the length of the shortest (which can
+// be negative)
 void raise_es(int nb_edges, scalar_t *el) {
   scalar_t min_es = el[0];
   for(int e = 1; e < nb_edges; e++) {
@@ -52,89 +60,119 @@ void raise_es(int nb_edges, scalar_t *el) {
   }
 }
 
+// Adds to every edge length the differential of the psi function on
+// it
 void add_dpsi_es(int nb_edges, scalar_t *el, int *ea, int *eb, scalar_t *psi) {
   for(int e = 0; e < nb_edges; e++) {
     el[e] += psi[ea[e]] - psi[eb[e]];
   }
 }
 
+// Finds the shortest path in the graph and returns in
+// result_edge_back, for each vertex, the edge to follow back from it
+// to reach the source with the shortest path, and in result_dist the
+// distance to the source. The edge lengths have to be positive.
 void find_shortest(int nb_vertices,
                    int nb_edges, scalar_t *el, int *ea, int *eb,
                    int source, int sink,
-                   int *edge_back, scalar_t *dist) {
+                   int *result_edge_back, scalar_t *result_dist) {
   for(int v = 0; v < nb_vertices; v++) {
-    // dist[v] = FLT_MAX;
-    dist[v] = 999;
+    result_dist[v] = FLT_MAX;
+    result_edge_back[v] = -1;
   }
 
-  dist[source] = 0;
-
-  for(int v = 0; v < nb_vertices; v++) {
-    edge_back[v] = -1;
-  }
+  result_dist[source] = 0;
 
+#ifdef DEBUG
   for(int e = 0; e < nb_edges; e++) {
-    ASSERT(el[e] >= 0);
+    if(el[e] < 0) abort();
   }
+#endif
 
   int nb_changes;
   scalar_t d;
   do {
     nb_changes = 0;
     for(int e = 0; e < nb_edges; e++) {
-      d = dist[ea[e]] + el[e];
-      if(d < dist[eb[e]]) {
+      d = result_dist[ea[e]] + el[e];
+      if(d < result_dist[eb[e]]) {
         nb_changes++;
-        dist[eb[e]] = d;
-        edge_back[eb[e]] = e;
+        result_dist[eb[e]] = d;
+        result_edge_back[eb[e]] = e;
       }
     }
   } while(nb_changes > 0);
 }
 
+// Iterates find_shortest as long as it finds paths of negative
+// lengths. Returns which edges are occupied by the superposition of
+// paths in result_edge_occupation.
+//
+// **WARNING** this routine changes the values of el, ea, and eb
+// (i.e. the occupied edges are inverted).
 void find_best_paths(int nb_vertices,
                      int nb_edges, scalar_t *el, int *ea, int *eb,
                      int source, int sink,
-                     int *edge_occupation) {
+                     int *result_edge_occupation) {
   scalar_t *dist = new scalar_t[nb_vertices];
   int *edge_back = new int[nb_vertices];
-  scalar_t *tmp_el = new scalar_t[nb_edges];
+  scalar_t *positive_el = new scalar_t[nb_edges];
+  scalar_t s;
+  int v;
 
   for(int e = 0; e < nb_edges; e++) {
-    tmp_el[e] = el[e];
-    edge_occupation[e] = 0;
+    positive_el[e] = el[e];
+    result_edge_occupation[e] = 0;
   }
 
-  raise_es(nb_edges, tmp_el);
+  raise_es(nb_edges, positive_el);
 
-  scalar_t s;
   do {
-    find_shortest(nb_vertices, nb_edges, tmp_el, ea, eb, source, sink, edge_back, dist);
-    add_dpsi_es(nb_edges, tmp_el, ea, eb, dist);
+    find_shortest(nb_vertices, nb_edges, positive_el, ea, eb, source, sink, edge_back, dist);
+    add_dpsi_es(nb_edges, positive_el, ea, eb, dist);
     s = 0.0;
+
+    // If the new path reaches the sink, we will backtrack on it to
+    // compute its score and invert edges
+
     if(edge_back[sink] >= 0) {
-      int v = sink;
+
+      v = sink;
       while(v != source) {
         int e = edge_back[v];
         ASSERT(eb[e] == v);
         v = ea[e];
-#ifdef VERBOSE
-        cout << "INVERTING " << ea[e] << " -> " << eb[e] << " [" << el[e] << "]" << endl;
-#endif
         s += el[e];
-        int k = eb[e];
-        eb[e] = ea[e];
-        ea[e] = k;
-        tmp_el[e] = - tmp_el[e];
-        edge_occupation[e] = 1 - edge_occupation[e];
       }
+
+      // We found a good path (score < 0), we revert the edges along
+      // the path and invert their occupation (since adding a path on
+      // a path already occupied means removing flow on it)
+
+      if(s < 0) {
+        v = sink;
+#ifdef VERBOSE
+        cout << "FOUND A PATH OF LENGTH " << s << endl;
+#endif
+        while(v != source) {
+          int e = edge_back[v];
+          ASSERT(eb[e] == v);
+          v = ea[e];
 #ifdef VERBOSE
-      cout << "FOUND A PATH OF LENGTH " << s << endl;
+          cout << "INVERTING " << ea[e] << " -> " << eb[e] << " [" << el[e] << "]" << endl;
 #endif
+          int k = eb[e];
+          eb[e] = ea[e];
+          ea[e] = k;
+          positive_el[e] = - positive_el[e];
+          el[e] = - el[e];
+          result_edge_occupation[e] = 1 - result_edge_occupation[e];
+        }
+      }
     }
   } while(s < 0);
 
-  delete[] tmp_el;
+  delete[] positive_el;
   delete[] dist;
   delete[] edge_back;
 }
@@ -146,19 +184,20 @@ int main(int argc, char **argv) {
     exit(EXIT_FAILURE);
   }
 
-  ifstream file(argv[1]);
+  ifstream *file = new ifstream(argv[1]);
 
   int nb_edges, nb_vertices;
   int source, sink;
 
-  if(file.good()) {
-    file >> nb_vertices >> nb_edges;
-    file >> source >> sink;
+  if(file->good()) {
 
-    cout << "nb_edges = " << nb_edges << endl;
-    cout << "nb_vertices = " << nb_vertices << endl;
-    cout << "source = " << source << endl;
-    cout << "sink = " << sink << endl;
+    (*file) >> nb_vertices >> nb_edges;
+    (*file) >> source >> sink;
+
+    cout << "INPUT nb_edges " << nb_edges << endl;
+    cout << "INPUT nb_vertices " << nb_vertices << endl;
+    cout << "INPUT source " << source << endl;
+    cout << "INPUT sink " << sink << endl;
 
     scalar_t *el = new scalar_t[nb_edges];
     int *ea = new int[nb_edges];
@@ -166,18 +205,27 @@ int main(int argc, char **argv) {
     int *edge_occupation = new int[nb_edges];
 
     for(int e = 0; e < nb_edges; e++) {
-      file >> ea[e] >> eb[e] >> el[e];
-#ifdef VERBOSE
-      cout << ea[e] << " -> " << eb[e] << " [" << el[e] << "]" << endl;
-#endif
+      (*file) >> ea[e] >> eb[e] >> el[e];
+      cout << "INPUT_EDGE " << ea[e] << " " << eb[e] << " " << el[e] << endl;
     }
 
     find_best_paths(nb_vertices, nb_edges, el, ea, eb, source, sink,
                     edge_occupation);
 
+#ifdef VERBOSE
+    // Sanity check on the overall resulting score (the edge lengths
+    // have been changed, hence should be the opposite of the sum of
+    // the path lengths)
+    scalar_t s = 0;
+    for(int e = 0; e < nb_edges; e++) {
+      if(edge_occupation[e]) s += el[e];
+    }
+    cout << "RESULT_SANITY_CHECK_SCORE " << s << endl;
+#endif
+
     for(int e = 0; e < nb_edges; e++) {
       if(edge_occupation[e]) {
-        cout << ea[e] << " " << eb[e] << endl;
+        cout << "RESULT_OCCUPIED_EDGE " << ea[e] << " " << eb[e] << endl;
       }
     }
 
@@ -189,9 +237,12 @@ int main(int argc, char **argv) {
   } else {
 
     cerr << "Can not open " << argv[1] << endl;
+
+    delete file;
     exit(EXIT_FAILURE);
 
   }
 
+  delete file;
   exit(EXIT_SUCCESS);
 }