X-Git-Url: https://fleuret.org/cgi-bin/gitweb/gitweb.cgi?a=blobdiff_plain;f=miniksp.cc;h=92e9de54a50adac278b322621a027caab6ae6dba;hb=58f6c3e7578dcef2007819bbdad95d144e547d44;hp=c6e83542004b49fc542b5928e8bbf32d660e8965;hpb=0428b7ae98e49c2ed98a24868f4ca2da8f3ab6e6;p=mtp.git diff --git a/miniksp.cc b/miniksp.cc index c6e8354..92e9de5 100644 --- a/miniksp.cc +++ b/miniksp.cc @@ -1,12 +1,26 @@ -//////////////////////////////////////////////////////////////////// -// START_IP_HEADER // -// // -// Written by Francois Fleuret // -// Contact 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 . // +// // +// Written by and Copyright (C) Francois Fleuret // +// Contact for comments & bug reports // +// // +// END_IP_HEADER // +/////////////////////////////////////////////////////////////////////////// + +// #define VERBOSE #include #include @@ -28,78 +42,207 @@ typedef float scalar_t; #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 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 *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; + 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); do { - find_shortest(nb_vertices, nb_edges, es, ea, eb, source, sink, dist, pred); - add_dpsi_es(nb_edges, es, ea, eb, dist); - for(int v = 0; v < nb_edges; v++) { - + 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) { + + 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(); + } while(s < 0); + delete[] positive_el; delete[] dist; - delete[] pred; + delete[] edge_back; } int main(int argc, char **argv) { - int nb_time_steps = 4; - int nb_locations = 10; - int nb_neighbors = 3; - // Add the source and sink - int nb_vertices = nb_time_steps * nb_locations + 2; - int nb_edges = nb_locations + (nb_time_steps - 1) * nb_locations + + if(argc < 2) { + cerr << argv[0] << " " << 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; + + 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); + + } + + delete file; + exit(EXIT_SUCCESS); }