2 ///////////////////////////////////////////////////////////////////////////
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. //
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. //
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/>. //
17 // Written by and Copyright (C) Francois Fleuret //
18 // Contact <francois.fleuret@idiap.ch> for comments & bug reports //
21 ///////////////////////////////////////////////////////////////////////////
34 typedef float scalar_t;
37 #define ASSERT(x) if(!(x)) { \
38 std::cerr << "ASSERT FAILED IN " << __FILE__ << ":" << __LINE__ << endl; \
47 // * el[e] is the length of edge e
48 // * ea[e] is its starting node
49 // * eb[e] is its destination node.
51 // Adds to all the edge length the length of the shortest (which can
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]);
58 for(int e = 0; e < nb_edges; e++) {
63 // Adds to every edge length the differential of the psi function on
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]];
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,
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;
84 result_dist[source] = 0;
87 for(int e = 0; e < nb_edges; e++) {
88 if(el[e] < 0) abort();
96 for(int e = 0; e < nb_edges; e++) {
97 d = result_dist[ea[e]] + el[e];
98 if(d < result_dist[eb[e]]) {
100 result_dist[eb[e]] = d;
101 result_edge_back[eb[e]] = e;
104 } while(nb_changes > 0);
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];
120 for(int e = 0; e < nb_edges; e++) {
121 positive_el[e] = el[e];
122 result_edge_occupation[e] = 0;
125 raise_es(nb_edges, positive_el);
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);
132 // If the new path reaches the sink, we will backtrack on it to
133 // compute its score and invert edges
135 if(edge_back[sink] >= 0) {
139 int e = edge_back[v];
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)
152 cout << "FOUND A PATH OF LENGTH " << s << endl;
155 int e = edge_back[v];
159 cout << "INVERTING " << ea[e] << " -> " << eb[e] << " [" << el[e] << "]" << endl;
164 positive_el[e] = - positive_el[e];
166 result_edge_occupation[e] = 1 - result_edge_occupation[e];
172 delete[] positive_el;
177 int main(int argc, char **argv) {
180 cerr << argv[0] << " <graph file>" << endl;
184 ifstream *file = new ifstream(argv[1]);
186 int nb_edges, nb_vertices;
191 (*file) >> nb_vertices >> nb_edges;
192 (*file) >> source >> sink;
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;
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];
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;
209 find_best_paths(nb_vertices, nb_edges, el, ea, eb, source, sink,
213 // Sanity check on the overall resulting score (the edge lengths
214 // have been changed, hence should be the opposite of the sum of
217 for(int e = 0; e < nb_edges; e++) {
218 if(edge_occupation[e]) s += el[e];
220 cout << "RESULT_SANITY_CHECK_SCORE " << s << endl;
223 for(int e = 0; e < nb_edges; e++) {
224 if(edge_occupation[e]) {
225 cout << "RESULT_OCCUPIED_EDGE " << ea[e] << " " << eb[e] << endl;
229 delete[] edge_occupation;
236 cerr << "Can not open " << argv[1] << endl;