-////////////////////////////////////////////////////////////////////
-// 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>
#define ASSERT(x)
#endif
-void raise_es(int nb_edges, scalar_t *es) {
- scalar_t min_es = es[0];
+// 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++) {
- 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) {
+// 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++) {
- es[e] += psi[ea[e]] - psi[eb[e]];
+ el[e] += psi[ea[e]] - psi[eb[e]];
}
}
+// Finds the shortest path in the graph and return 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 *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 *result_edge_back, scalar_t *result_dist) {
for(int v = 0; v < nb_vertices; v++) {
- dist[v] = FLT_MAX;
+ result_dist[v] = FLT_MAX;
+ result_edge_back[v] = -1;
}
- dist[source] = 0;
+ result_dist[source] = 0;
+#ifdef DEBUG
for(int e = 0; e < nb_edges; e++) {
- pred[e] = -1;
- ASSERT(se[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]] + es[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;
- pred[eb[e]] = ea[e];
+ result_dist[eb[e]] = d;
+ result_edge_back[eb[e]] = e;
}
}
} while(nb_changes > 0);
-
- ASSERT(pred[sink] >= 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 *es, int *ea, int *eb,
- int source, int sink) {
+ int nb_edges, scalar_t *el, int *ea, int *eb,
+ int source, int sink,
+ int *result_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 *positive_el = new scalar_t[nb_edges];
+ scalar_t s;
+ int v;
- raise_es(nb_edges, es);
+ for(int e = 0; e < nb_edges; e++) {
+ positive_el[e] = el[e];
+ result_edge_occupation[e] = 0;
+ }
+
+ raise_es(nb_edges, positive_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, positive_el, ea, eb, source, sink, edge_back, dist);
+ add_dpsi_es(nb_edges, positive_el, ea, eb, dist);
s = 0.0;
- for(int e = 0; e < nb_edges; e++) {
- if(pred[eb[e]] == ea[e]) {
- s += es[e];
- int k = eb[e];
- eb[e] = ea[e];
- ea[e] = k;
- es[e] = - es[e];
+
+ // If the new path reaches the sink, we will backtrack on it to
+ // compute its score and invert edges
+
+ if(edge_back[sink] >= 0) {
+
+ v = sink;
+ while(v != source) {
+ int e = edge_back[v];
+ ASSERT(eb[e] == v);
+ v = ea[e];
+ s += el[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 << "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[] positive_el;
delete[] dist;
- delete[] pred;
+ delete[] edge_back;
}
int main(int argc, char **argv) {
- ifstream file(argv[1]);
+
+ if(argc < 2) {
+ cerr << argv[0] << " <graph file>" << endl;
+ exit(EXIT_FAILURE);
+ }
+
+ 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;
- scalar_t *es = new scalar_t[nb_edges];
- int *ea = new int[nb_edges];
- int *eb = new int[nb_edges];
+ 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];
+ 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];
+ 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 << "RESULT_OCCUPIED_EDGE " << ea[e] << " " << eb[e] << endl;
+ }
+ }
+
+ delete[] edge_occupation;
+ delete[] el;
+ delete[] ea;
+ delete[] eb;
+
+ } else {
+
+ cerr << "Can not open " << argv[1] << endl;
+
+ delete file;
+ 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;
+ delete file;
+ exit(EXIT_SUCCESS);
}
projects / mtp.git / blobdiff