Update.

diff --git a/miniksp.cc b/miniksp.cc
index 513da0e..f562423 100644 (file)
--- a/miniksp.cc
+++ b/miniksp.cc
@@ -1,12 +1,26 @@
 
-////////////////////////////////////////////////////////////////////
-// START_IP_HEADER                                                //
-//                                                                //
-// Written by Francois Fleuret                                    //
-// Contact <francois.fleuret@idiap.ch> for comments & bug reports //
-//                                                                //
-// END_IP_HEADER                                                  //
-////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////
+// START_IP_HEADER                                                       //
+//                                                                       //
+// This program is free software: you can redistribute it and/or modify  //
+// it under the terms of the version 3 of the GNU General Public License //
+// as published by the Free Software Foundation.                         //
+//                                                                       //
+// This program is distributed in the hope that it will be useful, but   //
+// WITHOUT ANY WARRANTY; without even the implied warranty of            //
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU      //
+// General Public License for more details.                              //
+//                                                                       //
+// You should have received a copy of the GNU General Public License     //
+// along with this program. If not, see <http://www.gnu.org/licenses/>.  //
+//                                                                       //
+// Written by and Copyright (C) Francois Fleuret                         //
+// Contact <francois.fleuret@idiap.ch> for comments & bug reports        //
+//                                                                       //
+// END_IP_HEADER                                                         //
+///////////////////////////////////////////////////////////////////////////
+
+#define VERBOSE
 
 #include <iostream>
 #include <fstream>
@@ -28,35 +42,39 @@ typedef float scalar_t;
 #define ASSERT(x)
 #endif
 
-void raise_es(int nb_edges, scalar_t *es) {
-  scalar_t min_es = es[0];
+void raise_es(int nb_edges, scalar_t *el) {
+  scalar_t min_es = el[0];
   for(int e = 1; e < nb_edges; e++) {
-    min_es = min(min_es, es[e]);
+    min_es = min(min_es, el[e]);
   }
   for(int e = 0; e < nb_edges; e++) {
-    es[e] -= min_es;
+    el[e] -= min_es;
   }
 }
 
-void add_dpsi_es(int nb_edges, scalar_t *es, int *ea, int *eb, scalar_t *psi) {
+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++) {
-    es[e] += psi[ea[e]] - psi[eb[e]];
+    el[e] += psi[ea[e]] - psi[eb[e]];
   }
 }
 
 void find_shortest(int nb_vertices,
-                   int nb_edges, scalar_t *es, int *ea, int *eb,
+                   int nb_edges, scalar_t *el, int *ea, int *eb,
                    int source, int sink,
-                   int *pred, scalar_t *dist) {
+                   int *edge_back, scalar_t *dist) {
   for(int v = 0; v < nb_vertices; v++) {
-    dist[v] = FLT_MAX;
+    // dist[v] = FLT_MAX;
+    dist[v] = 999;
   }
 
   dist[source] = 0;
 
+  for(int v = 0; v < nb_vertices; v++) {
+    edge_back[v] = -1;
+  }
+
   for(int e = 0; e < nb_edges; e++) {
-    pred[e] = -1;
-    ASSERT(se[e] >= 0);
+    ASSERT(el[e] >= 0);
   }
 
   int nb_changes;
@@ -64,71 +82,116 @@ void find_shortest(int nb_vertices,
   do {
     nb_changes = 0;
     for(int e = 0; e < nb_edges; e++) {
-      d = dist[ea[e]] + es[e];
+      d = dist[ea[e]] + el[e];
       if(d < dist[eb[e]]) {
         nb_changes++;
         dist[eb[e]] = d;
-        pred[eb[e]] = ea[e];
+        edge_back[eb[e]] = e;
       }
     }
   } while(nb_changes > 0);
-
-  ASSERT(pred[sink] >= 0);
 }
 
 void find_best_paths(int nb_vertices,
-                     int nb_edges, scalar_t *es, int *ea, int *eb,
-                     int source, int sink) {
+                     int nb_edges, scalar_t *el, int *ea, int *eb,
+                     int source, int sink,
+                     int *edge_occupation) {
   scalar_t *dist = new scalar_t[nb_vertices];
-  int *pred = new int[nb_vertices];
+  int *edge_back = new int[nb_vertices];
+  scalar_t *tmp_el = new scalar_t[nb_edges];
+
+  for(int e = 0; e < nb_edges; e++) {
+    tmp_el[e] = el[e];
+    edge_occupation[e] = 0;
+  }
 
-  raise_es(nb_edges, es);
+  raise_es(nb_edges, tmp_el);
 
   scalar_t s;
   do {
-    find_shortest(nb_vertices, nb_edges, es, ea, eb, source, sink, pred, dist);
-    add_dpsi_es(nb_edges, es, ea, eb, dist);
+    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);
     s = 0.0;
-    for(int e = 0; e < nb_edges; e++) {
-      if(pred[eb[e]] == ea[e]) {
-        s += es[e];
+    if(edge_back[sink] >= 0) {
+      int 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;
-        es[e] = - es[e];
+        tmp_el[e] = - tmp_el[e];
+        edge_occupation[e] = 1 - edge_occupation[e];
       }
+#ifdef VERBOSE
+      cout << "FOUND A PATH OF LENGTH " << s << endl;
+#endif
     }
   } while(s < 0);
 
+  delete[] tmp_el;
   delete[] dist;
-  delete[] pred;
+  delete[] edge_back;
 }
 
 int main(int argc, char **argv) {
+
+  if(argc < 2) {
+    cerr << argv[0] << " <graph file>" << endl;
+    exit(EXIT_FAILURE);
+  }
+
   ifstream file(argv[1]);
+
   int nb_edges, nb_vertices;
   int source, sink;
 
   if(file.good()) {
     file >> nb_vertices >> nb_edges;
     file >> source >> sink;
-  }
 
-  cout << "nb_edges = " << nb_edges << endl;
-  cout << "nb_vertices = " << nb_vertices << endl;
-  cout << "source = " << source << endl;
-  cout << "sink = " << sink << endl;
+    cout << "nb_edges = " << nb_edges << endl;
+    cout << "nb_vertices = " << nb_vertices << endl;
+    cout << "source = " << source << endl;
+    cout << "sink = " << sink << endl;
 
-  scalar_t *es = new scalar_t[nb_edges];
-  int *ea = new int[nb_edges];
-  int *eb = new int[nb_edges];
+    scalar_t *el = new scalar_t[nb_edges];
+    int *ea = new int[nb_edges];
+    int *eb = new int[nb_edges];
+    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
+    }
+
+    find_best_paths(nb_vertices, nb_edges, el, ea, eb, source, sink,
+                    edge_occupation);
+
+    for(int e = 0; e < nb_edges; e++) {
+      if(edge_occupation[e]) {
+        cout << ea[e] << " " << eb[e] << endl;
+      }
+    }
+
+    delete[] edge_occupation;
+    delete[] el;
+    delete[] ea;
+    delete[] eb;
+
+  } else {
+
+    cerr << "Can not open " << argv[1] << endl;
+    exit(EXIT_FAILURE);
 
-  for(int e = 0; e < nb_edges; e++) {
-    file >> ea[e] >> eb[e] >> es[e];
-    cout << ea[e] << " -> " << eb[e] << " [" << es[e] << "]" << endl;
   }
 
-  delete[] es;
-  delete[] ea;
-  delete[] eb;
+  exit(EXIT_SUCCESS);
 }