1 | //======================================================================= |
---|
2 | // Copyright 2000 University of Notre Dame. |
---|
3 | // Authors: Jeremy G. Siek, Andrew Lumsdaine, Lie-Quan Lee |
---|
4 | // |
---|
5 | // Distributed under the Boost Software License, Version 1.0. (See |
---|
6 | // accompanying file LICENSE_1_0.txt or copy at |
---|
7 | // http://www.boost.org/LICENSE_1_0.txt) |
---|
8 | //======================================================================= |
---|
9 | |
---|
10 | #ifndef BOOST_PUSH_RELABEL_MAX_FLOW_HPP |
---|
11 | #define BOOST_PUSH_RELABEL_MAX_FLOW_HPP |
---|
12 | |
---|
13 | #include <boost/config.hpp> |
---|
14 | #include <cassert> |
---|
15 | #include <vector> |
---|
16 | #include <list> |
---|
17 | #include <iosfwd> |
---|
18 | #include <algorithm> // for std::min and std::max |
---|
19 | |
---|
20 | #include <boost/pending/queue.hpp> |
---|
21 | #include <boost/limits.hpp> |
---|
22 | #include <boost/graph/graph_concepts.hpp> |
---|
23 | #include <boost/graph/named_function_params.hpp> |
---|
24 | |
---|
25 | namespace boost { |
---|
26 | |
---|
27 | namespace detail { |
---|
28 | |
---|
29 | // This implementation is based on Goldberg's |
---|
30 | // "On Implementing Push-Relabel Method for the Maximum Flow Problem" |
---|
31 | // by B.V. Cherkassky and A.V. Goldberg, IPCO '95, pp. 157--171 |
---|
32 | // and on the h_prf.c and hi_pr.c code written by the above authors. |
---|
33 | |
---|
34 | // This implements the highest-label version of the push-relabel method |
---|
35 | // with the global relabeling and gap relabeling heuristics. |
---|
36 | |
---|
37 | // The terms "rank", "distance", "height" are synonyms in |
---|
38 | // Goldberg's implementation, paper and in the CLR. A "layer" is a |
---|
39 | // group of vertices with the same distance. The vertices in each |
---|
40 | // layer are categorized as active or inactive. An active vertex |
---|
41 | // has positive excess flow and its distance is less than n (it is |
---|
42 | // not blocked). |
---|
43 | |
---|
44 | template <class Vertex> |
---|
45 | struct preflow_layer { |
---|
46 | std::list<Vertex> active_vertices; |
---|
47 | std::list<Vertex> inactive_vertices; |
---|
48 | }; |
---|
49 | |
---|
50 | template <class Graph, |
---|
51 | class EdgeCapacityMap, // integer value type |
---|
52 | class ResidualCapacityEdgeMap, |
---|
53 | class ReverseEdgeMap, |
---|
54 | class VertexIndexMap, // vertex_descriptor -> integer |
---|
55 | class FlowValue> |
---|
56 | class push_relabel |
---|
57 | { |
---|
58 | public: |
---|
59 | typedef graph_traits<Graph> Traits; |
---|
60 | typedef typename Traits::vertex_descriptor vertex_descriptor; |
---|
61 | typedef typename Traits::edge_descriptor edge_descriptor; |
---|
62 | typedef typename Traits::vertex_iterator vertex_iterator; |
---|
63 | typedef typename Traits::out_edge_iterator out_edge_iterator; |
---|
64 | typedef typename Traits::vertices_size_type vertices_size_type; |
---|
65 | typedef typename Traits::edges_size_type edges_size_type; |
---|
66 | |
---|
67 | typedef preflow_layer<vertex_descriptor> Layer; |
---|
68 | typedef std::vector< Layer > LayerArray; |
---|
69 | typedef typename LayerArray::iterator layer_iterator; |
---|
70 | typedef typename LayerArray::size_type distance_size_type; |
---|
71 | |
---|
72 | typedef color_traits<default_color_type> ColorTraits; |
---|
73 | |
---|
74 | //======================================================================= |
---|
75 | // Some helper predicates |
---|
76 | |
---|
77 | inline bool is_admissible(vertex_descriptor u, vertex_descriptor v) { |
---|
78 | return distance[u] == distance[v] + 1; |
---|
79 | } |
---|
80 | inline bool is_residual_edge(edge_descriptor a) { |
---|
81 | return 0 < residual_capacity[a]; |
---|
82 | } |
---|
83 | inline bool is_saturated(edge_descriptor a) { |
---|
84 | return residual_capacity[a] == 0; |
---|
85 | } |
---|
86 | |
---|
87 | //======================================================================= |
---|
88 | // Layer List Management Functions |
---|
89 | |
---|
90 | typedef typename std::list<vertex_descriptor>::iterator list_iterator; |
---|
91 | |
---|
92 | void add_to_active_list(vertex_descriptor u, Layer& layer) { |
---|
93 | BOOST_USING_STD_MIN(); |
---|
94 | BOOST_USING_STD_MAX(); |
---|
95 | layer.active_vertices.push_front(u); |
---|
96 | max_active = max BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], max_active); |
---|
97 | min_active = min BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], min_active); |
---|
98 | layer_list_ptr[u] = layer.active_vertices.begin(); |
---|
99 | } |
---|
100 | void remove_from_active_list(vertex_descriptor u) { |
---|
101 | layers[distance[u]].active_vertices.erase(layer_list_ptr[u]); |
---|
102 | } |
---|
103 | |
---|
104 | void add_to_inactive_list(vertex_descriptor u, Layer& layer) { |
---|
105 | layer.inactive_vertices.push_front(u); |
---|
106 | layer_list_ptr[u] = layer.inactive_vertices.begin(); |
---|
107 | } |
---|
108 | void remove_from_inactive_list(vertex_descriptor u) { |
---|
109 | layers[distance[u]].inactive_vertices.erase(layer_list_ptr[u]); |
---|
110 | } |
---|
111 | |
---|
112 | //======================================================================= |
---|
113 | // initialization |
---|
114 | push_relabel(Graph& g_, |
---|
115 | EdgeCapacityMap cap, |
---|
116 | ResidualCapacityEdgeMap res, |
---|
117 | ReverseEdgeMap rev, |
---|
118 | vertex_descriptor src_, |
---|
119 | vertex_descriptor sink_, |
---|
120 | VertexIndexMap idx) |
---|
121 | : g(g_), n(num_vertices(g_)), capacity(cap), src(src_), sink(sink_), |
---|
122 | index(idx), |
---|
123 | excess_flow(num_vertices(g_)), |
---|
124 | current(num_vertices(g_), out_edges(*vertices(g_).first, g_).second), |
---|
125 | distance(num_vertices(g_)), |
---|
126 | color(num_vertices(g_)), |
---|
127 | reverse_edge(rev), |
---|
128 | residual_capacity(res), |
---|
129 | layers(num_vertices(g_)), |
---|
130 | layer_list_ptr(num_vertices(g_), |
---|
131 | layers.front().inactive_vertices.end()), |
---|
132 | push_count(0), update_count(0), relabel_count(0), |
---|
133 | gap_count(0), gap_node_count(0), |
---|
134 | work_since_last_update(0) |
---|
135 | { |
---|
136 | vertex_iterator u_iter, u_end; |
---|
137 | // Don't count the reverse edges |
---|
138 | edges_size_type m = num_edges(g) / 2; |
---|
139 | nm = alpha() * n + m; |
---|
140 | |
---|
141 | // Initialize flow to zero which means initializing |
---|
142 | // the residual capacity to equal the capacity. |
---|
143 | out_edge_iterator ei, e_end; |
---|
144 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) |
---|
145 | for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) { |
---|
146 | residual_capacity[*ei] = capacity[*ei]; |
---|
147 | } |
---|
148 | |
---|
149 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
---|
150 | vertex_descriptor u = *u_iter; |
---|
151 | excess_flow[u] = 0; |
---|
152 | current[u] = out_edges(u, g).first; |
---|
153 | } |
---|
154 | |
---|
155 | bool overflow_detected = false; |
---|
156 | FlowValue test_excess = 0; |
---|
157 | |
---|
158 | out_edge_iterator a_iter, a_end; |
---|
159 | for (tie(a_iter, a_end) = out_edges(src, g); a_iter != a_end; ++a_iter) |
---|
160 | if (target(*a_iter, g) != src) |
---|
161 | test_excess += residual_capacity[*a_iter]; |
---|
162 | if (test_excess > (std::numeric_limits<FlowValue>::max)()) |
---|
163 | overflow_detected = true; |
---|
164 | |
---|
165 | if (overflow_detected) |
---|
166 | excess_flow[src] = (std::numeric_limits<FlowValue>::max)(); |
---|
167 | else { |
---|
168 | excess_flow[src] = 0; |
---|
169 | for (tie(a_iter, a_end) = out_edges(src, g); |
---|
170 | a_iter != a_end; ++a_iter) { |
---|
171 | edge_descriptor a = *a_iter; |
---|
172 | if (target(a, g) != src) { |
---|
173 | ++push_count; |
---|
174 | FlowValue delta = residual_capacity[a]; |
---|
175 | residual_capacity[a] -= delta; |
---|
176 | residual_capacity[reverse_edge[a]] += delta; |
---|
177 | excess_flow[target(a, g)] += delta; |
---|
178 | } |
---|
179 | } |
---|
180 | } |
---|
181 | max_distance = num_vertices(g) - 1; |
---|
182 | max_active = 0; |
---|
183 | min_active = n; |
---|
184 | |
---|
185 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
---|
186 | vertex_descriptor u = *u_iter; |
---|
187 | if (u == sink) { |
---|
188 | distance[u] = 0; |
---|
189 | continue; |
---|
190 | } else if (u == src && !overflow_detected) |
---|
191 | distance[u] = n; |
---|
192 | else |
---|
193 | distance[u] = 1; |
---|
194 | |
---|
195 | if (excess_flow[u] > 0) |
---|
196 | add_to_active_list(u, layers[1]); |
---|
197 | else if (distance[u] < n) |
---|
198 | add_to_inactive_list(u, layers[1]); |
---|
199 | } |
---|
200 | |
---|
201 | } // push_relabel constructor |
---|
202 | |
---|
203 | //======================================================================= |
---|
204 | // This is a breadth-first search over the residual graph |
---|
205 | // (well, actually the reverse of the residual graph). |
---|
206 | // Would be cool to have a graph view adaptor for hiding certain |
---|
207 | // edges, like the saturated (non-residual) edges in this case. |
---|
208 | // Goldberg's implementation abused "distance" for the coloring. |
---|
209 | void global_distance_update() |
---|
210 | { |
---|
211 | BOOST_USING_STD_MAX(); |
---|
212 | ++update_count; |
---|
213 | vertex_iterator u_iter, u_end; |
---|
214 | for (tie(u_iter,u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
---|
215 | color[*u_iter] = ColorTraits::white(); |
---|
216 | distance[*u_iter] = n; |
---|
217 | } |
---|
218 | color[sink] = ColorTraits::gray(); |
---|
219 | distance[sink] = 0; |
---|
220 | |
---|
221 | for (distance_size_type l = 0; l <= max_distance; ++l) { |
---|
222 | layers[l].active_vertices.clear(); |
---|
223 | layers[l].inactive_vertices.clear(); |
---|
224 | } |
---|
225 | |
---|
226 | max_distance = max_active = 0; |
---|
227 | min_active = n; |
---|
228 | |
---|
229 | Q.push(sink); |
---|
230 | while (! Q.empty()) { |
---|
231 | vertex_descriptor u = Q.top(); |
---|
232 | Q.pop(); |
---|
233 | distance_size_type d_v = distance[u] + 1; |
---|
234 | |
---|
235 | out_edge_iterator ai, a_end; |
---|
236 | for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) { |
---|
237 | edge_descriptor a = *ai; |
---|
238 | vertex_descriptor v = target(a, g); |
---|
239 | if (color[v] == ColorTraits::white() |
---|
240 | && is_residual_edge(reverse_edge[a])) { |
---|
241 | distance[v] = d_v; |
---|
242 | color[v] = ColorTraits::gray(); |
---|
243 | current[v] = out_edges(v, g).first; |
---|
244 | max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(d_v, max_distance); |
---|
245 | |
---|
246 | if (excess_flow[v] > 0) |
---|
247 | add_to_active_list(v, layers[d_v]); |
---|
248 | else |
---|
249 | add_to_inactive_list(v, layers[d_v]); |
---|
250 | |
---|
251 | Q.push(v); |
---|
252 | } |
---|
253 | } |
---|
254 | } |
---|
255 | } // global_distance_update() |
---|
256 | |
---|
257 | //======================================================================= |
---|
258 | // This function is called "push" in Goldberg's h_prf implementation, |
---|
259 | // but it is called "discharge" in the paper and in hi_pr.c. |
---|
260 | void discharge(vertex_descriptor u) |
---|
261 | { |
---|
262 | assert(excess_flow[u] > 0); |
---|
263 | while (1) { |
---|
264 | out_edge_iterator ai, ai_end; |
---|
265 | for (ai = current[u], ai_end = out_edges(u, g).second; |
---|
266 | ai != ai_end; ++ai) { |
---|
267 | edge_descriptor a = *ai; |
---|
268 | if (is_residual_edge(a)) { |
---|
269 | vertex_descriptor v = target(a, g); |
---|
270 | if (is_admissible(u, v)) { |
---|
271 | ++push_count; |
---|
272 | if (v != sink && excess_flow[v] == 0) { |
---|
273 | remove_from_inactive_list(v); |
---|
274 | add_to_active_list(v, layers[distance[v]]); |
---|
275 | } |
---|
276 | push_flow(a); |
---|
277 | if (excess_flow[u] == 0) |
---|
278 | break; |
---|
279 | } |
---|
280 | } |
---|
281 | } // for out_edges of i starting from current |
---|
282 | |
---|
283 | Layer& layer = layers[distance[u]]; |
---|
284 | distance_size_type du = distance[u]; |
---|
285 | |
---|
286 | if (ai == ai_end) { // i must be relabeled |
---|
287 | relabel_distance(u); |
---|
288 | if (layer.active_vertices.empty() |
---|
289 | && layer.inactive_vertices.empty()) |
---|
290 | gap(du); |
---|
291 | if (distance[u] == n) |
---|
292 | break; |
---|
293 | } else { // i is no longer active |
---|
294 | current[u] = ai; |
---|
295 | add_to_inactive_list(u, layer); |
---|
296 | break; |
---|
297 | } |
---|
298 | } // while (1) |
---|
299 | } // discharge() |
---|
300 | |
---|
301 | //======================================================================= |
---|
302 | // This corresponds to the "push" update operation of the paper, |
---|
303 | // not the "push" function in Goldberg's h_prf.c implementation. |
---|
304 | // The idea is to push the excess flow from from vertex u to v. |
---|
305 | void push_flow(edge_descriptor u_v) |
---|
306 | { |
---|
307 | vertex_descriptor |
---|
308 | u = source(u_v, g), |
---|
309 | v = target(u_v, g); |
---|
310 | |
---|
311 | BOOST_USING_STD_MIN(); |
---|
312 | FlowValue flow_delta |
---|
313 | = min BOOST_PREVENT_MACRO_SUBSTITUTION(excess_flow[u], residual_capacity[u_v]); |
---|
314 | |
---|
315 | residual_capacity[u_v] -= flow_delta; |
---|
316 | residual_capacity[reverse_edge[u_v]] += flow_delta; |
---|
317 | |
---|
318 | excess_flow[u] -= flow_delta; |
---|
319 | excess_flow[v] += flow_delta; |
---|
320 | } // push_flow() |
---|
321 | |
---|
322 | //======================================================================= |
---|
323 | // The main purpose of this routine is to set distance[v] |
---|
324 | // to the smallest value allowed by the valid labeling constraints, |
---|
325 | // which are: |
---|
326 | // distance[t] = 0 |
---|
327 | // distance[u] <= distance[v] + 1 for every residual edge (u,v) |
---|
328 | // |
---|
329 | distance_size_type relabel_distance(vertex_descriptor u) |
---|
330 | { |
---|
331 | BOOST_USING_STD_MAX(); |
---|
332 | ++relabel_count; |
---|
333 | work_since_last_update += beta(); |
---|
334 | |
---|
335 | distance_size_type min_distance = num_vertices(g); |
---|
336 | distance[u] = min_distance; |
---|
337 | |
---|
338 | // Examine the residual out-edges of vertex i, choosing the |
---|
339 | // edge whose target vertex has the minimal distance. |
---|
340 | out_edge_iterator ai, a_end, min_edge_iter; |
---|
341 | for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) { |
---|
342 | ++work_since_last_update; |
---|
343 | edge_descriptor a = *ai; |
---|
344 | vertex_descriptor v = target(a, g); |
---|
345 | if (is_residual_edge(a) && distance[v] < min_distance) { |
---|
346 | min_distance = distance[v]; |
---|
347 | min_edge_iter = ai; |
---|
348 | } |
---|
349 | } |
---|
350 | ++min_distance; |
---|
351 | if (min_distance < n) { |
---|
352 | distance[u] = min_distance; // this is the main action |
---|
353 | current[u] = min_edge_iter; |
---|
354 | max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(min_distance, max_distance); |
---|
355 | } |
---|
356 | return min_distance; |
---|
357 | } // relabel_distance() |
---|
358 | |
---|
359 | //======================================================================= |
---|
360 | // cleanup beyond the gap |
---|
361 | void gap(distance_size_type empty_distance) |
---|
362 | { |
---|
363 | ++gap_count; |
---|
364 | |
---|
365 | distance_size_type r; // distance of layer before the current layer |
---|
366 | r = empty_distance - 1; |
---|
367 | |
---|
368 | // Set the distance for the vertices beyond the gap to "infinity". |
---|
369 | for (layer_iterator l = layers.begin() + empty_distance + 1; |
---|
370 | l < layers.begin() + max_distance; ++l) { |
---|
371 | list_iterator i; |
---|
372 | for (i = l->inactive_vertices.begin(); |
---|
373 | i != l->inactive_vertices.end(); ++i) { |
---|
374 | distance[*i] = n; |
---|
375 | ++gap_node_count; |
---|
376 | } |
---|
377 | l->inactive_vertices.clear(); |
---|
378 | } |
---|
379 | max_distance = r; |
---|
380 | max_active = r; |
---|
381 | } |
---|
382 | |
---|
383 | //======================================================================= |
---|
384 | // This is the core part of the algorithm, "phase one". |
---|
385 | FlowValue maximum_preflow() |
---|
386 | { |
---|
387 | work_since_last_update = 0; |
---|
388 | |
---|
389 | while (max_active >= min_active) { // "main" loop |
---|
390 | |
---|
391 | Layer& layer = layers[max_active]; |
---|
392 | list_iterator u_iter = layer.active_vertices.begin(); |
---|
393 | |
---|
394 | if (u_iter == layer.active_vertices.end()) |
---|
395 | --max_active; |
---|
396 | else { |
---|
397 | vertex_descriptor u = *u_iter; |
---|
398 | remove_from_active_list(u); |
---|
399 | |
---|
400 | discharge(u); |
---|
401 | |
---|
402 | if (work_since_last_update * global_update_frequency() > nm) { |
---|
403 | global_distance_update(); |
---|
404 | work_since_last_update = 0; |
---|
405 | } |
---|
406 | } |
---|
407 | } // while (max_active >= min_active) |
---|
408 | |
---|
409 | return excess_flow[sink]; |
---|
410 | } // maximum_preflow() |
---|
411 | |
---|
412 | //======================================================================= |
---|
413 | // remove excess flow, the "second phase" |
---|
414 | // This does a DFS on the reverse flow graph of nodes with excess flow. |
---|
415 | // If a cycle is found, cancel it. |
---|
416 | // Return the nodes with excess flow in topological order. |
---|
417 | // |
---|
418 | // Unlike the prefl_to_flow() implementation, we use |
---|
419 | // "color" instead of "distance" for the DFS labels |
---|
420 | // "parent" instead of nl_prev for the DFS tree |
---|
421 | // "topo_next" instead of nl_next for the topological ordering |
---|
422 | void convert_preflow_to_flow() |
---|
423 | { |
---|
424 | vertex_iterator u_iter, u_end; |
---|
425 | out_edge_iterator ai, a_end; |
---|
426 | |
---|
427 | vertex_descriptor r, restart, u; |
---|
428 | |
---|
429 | std::vector<vertex_descriptor> parent(n); |
---|
430 | std::vector<vertex_descriptor> topo_next(n); |
---|
431 | |
---|
432 | vertex_descriptor tos(parent[0]), |
---|
433 | bos(parent[0]); // bogus initialization, just to avoid warning |
---|
434 | bool bos_null = true; |
---|
435 | |
---|
436 | // handle self-loops |
---|
437 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) |
---|
438 | for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) |
---|
439 | if (target(*ai, g) == *u_iter) |
---|
440 | residual_capacity[*ai] = capacity[*ai]; |
---|
441 | |
---|
442 | // initialize |
---|
443 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
---|
444 | u = *u_iter; |
---|
445 | color[u] = ColorTraits::white(); |
---|
446 | parent[u] = u; |
---|
447 | current[u] = out_edges(u, g).first; |
---|
448 | } |
---|
449 | // eliminate flow cycles and topologically order the vertices |
---|
450 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
---|
451 | u = *u_iter; |
---|
452 | if (color[u] == ColorTraits::white() |
---|
453 | && excess_flow[u] > 0 |
---|
454 | && u != src && u != sink ) { |
---|
455 | r = u; |
---|
456 | color[r] = ColorTraits::gray(); |
---|
457 | while (1) { |
---|
458 | for (; current[u] != out_edges(u, g).second; ++current[u]) { |
---|
459 | edge_descriptor a = *current[u]; |
---|
460 | if (capacity[a] == 0 && is_residual_edge(a)) { |
---|
461 | vertex_descriptor v = target(a, g); |
---|
462 | if (color[v] == ColorTraits::white()) { |
---|
463 | color[v] = ColorTraits::gray(); |
---|
464 | parent[v] = u; |
---|
465 | u = v; |
---|
466 | break; |
---|
467 | } else if (color[v] == ColorTraits::gray()) { |
---|
468 | // find minimum flow on the cycle |
---|
469 | FlowValue delta = residual_capacity[a]; |
---|
470 | while (1) { |
---|
471 | BOOST_USING_STD_MIN(); |
---|
472 | delta = min BOOST_PREVENT_MACRO_SUBSTITUTION(delta, residual_capacity[*current[v]]); |
---|
473 | if (v == u) |
---|
474 | break; |
---|
475 | else |
---|
476 | v = target(*current[v], g); |
---|
477 | } |
---|
478 | // remove delta flow units |
---|
479 | v = u; |
---|
480 | while (1) { |
---|
481 | a = *current[v]; |
---|
482 | residual_capacity[a] -= delta; |
---|
483 | residual_capacity[reverse_edge[a]] += delta; |
---|
484 | v = target(a, g); |
---|
485 | if (v == u) |
---|
486 | break; |
---|
487 | } |
---|
488 | |
---|
489 | // back-out of DFS to the first saturated edge |
---|
490 | restart = u; |
---|
491 | for (v = target(*current[u], g); v != u; v = target(a, g)){ |
---|
492 | a = *current[v]; |
---|
493 | if (color[v] == ColorTraits::white() |
---|
494 | || is_saturated(a)) { |
---|
495 | color[target(*current[v], g)] = ColorTraits::white(); |
---|
496 | if (color[v] != ColorTraits::white()) |
---|
497 | restart = v; |
---|
498 | } |
---|
499 | } |
---|
500 | if (restart != u) { |
---|
501 | u = restart; |
---|
502 | ++current[u]; |
---|
503 | break; |
---|
504 | } |
---|
505 | } // else if (color[v] == ColorTraits::gray()) |
---|
506 | } // if (capacity[a] == 0 ... |
---|
507 | } // for out_edges(u, g) (though "u" changes during loop) |
---|
508 | |
---|
509 | if (current[u] == out_edges(u, g).second) { |
---|
510 | // scan of i is complete |
---|
511 | color[u] = ColorTraits::black(); |
---|
512 | if (u != src) { |
---|
513 | if (bos_null) { |
---|
514 | bos = u; |
---|
515 | bos_null = false; |
---|
516 | tos = u; |
---|
517 | } else { |
---|
518 | topo_next[u] = tos; |
---|
519 | tos = u; |
---|
520 | } |
---|
521 | } |
---|
522 | if (u != r) { |
---|
523 | u = parent[u]; |
---|
524 | ++current[u]; |
---|
525 | } else |
---|
526 | break; |
---|
527 | } |
---|
528 | } // while (1) |
---|
529 | } // if (color[u] == white && excess_flow[u] > 0 & ...) |
---|
530 | } // for all vertices in g |
---|
531 | |
---|
532 | // return excess flows |
---|
533 | // note that the sink is not on the stack |
---|
534 | if (! bos_null) { |
---|
535 | for (u = tos; u != bos; u = topo_next[u]) { |
---|
536 | ai = out_edges(u, g).first; |
---|
537 | while (excess_flow[u] > 0 && ai != out_edges(u, g).second) { |
---|
538 | if (capacity[*ai] == 0 && is_residual_edge(*ai)) |
---|
539 | push_flow(*ai); |
---|
540 | ++ai; |
---|
541 | } |
---|
542 | } |
---|
543 | // do the bottom |
---|
544 | u = bos; |
---|
545 | ai = out_edges(u, g).first; |
---|
546 | while (excess_flow[u] > 0) { |
---|
547 | if (capacity[*ai] == 0 && is_residual_edge(*ai)) |
---|
548 | push_flow(*ai); |
---|
549 | ++ai; |
---|
550 | } |
---|
551 | } |
---|
552 | |
---|
553 | } // convert_preflow_to_flow() |
---|
554 | |
---|
555 | //======================================================================= |
---|
556 | inline bool is_flow() |
---|
557 | { |
---|
558 | vertex_iterator u_iter, u_end; |
---|
559 | out_edge_iterator ai, a_end; |
---|
560 | |
---|
561 | // check edge flow values |
---|
562 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
---|
563 | for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) { |
---|
564 | edge_descriptor a = *ai; |
---|
565 | if (capacity[a] > 0) |
---|
566 | if ((residual_capacity[a] + residual_capacity[reverse_edge[a]] |
---|
567 | != capacity[a] + capacity[reverse_edge[a]]) |
---|
568 | || (residual_capacity[a] < 0) |
---|
569 | || (residual_capacity[reverse_edge[a]] < 0)) |
---|
570 | return false; |
---|
571 | } |
---|
572 | } |
---|
573 | |
---|
574 | // check conservation |
---|
575 | FlowValue sum; |
---|
576 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
---|
577 | vertex_descriptor u = *u_iter; |
---|
578 | if (u != src && u != sink) { |
---|
579 | if (excess_flow[u] != 0) |
---|
580 | return false; |
---|
581 | sum = 0; |
---|
582 | for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) |
---|
583 | if (capacity[*ai] > 0) |
---|
584 | sum -= capacity[*ai] - residual_capacity[*ai]; |
---|
585 | else |
---|
586 | sum += residual_capacity[*ai]; |
---|
587 | |
---|
588 | if (excess_flow[u] != sum) |
---|
589 | return false; |
---|
590 | } |
---|
591 | } |
---|
592 | |
---|
593 | return true; |
---|
594 | } // is_flow() |
---|
595 | |
---|
596 | bool is_optimal() { |
---|
597 | // check if mincut is saturated... |
---|
598 | global_distance_update(); |
---|
599 | return distance[src] >= n; |
---|
600 | } |
---|
601 | |
---|
602 | void print_statistics(std::ostream& os) const { |
---|
603 | os << "pushes: " << push_count << std::endl |
---|
604 | << "relabels: " << relabel_count << std::endl |
---|
605 | << "updates: " << update_count << std::endl |
---|
606 | << "gaps: " << gap_count << std::endl |
---|
607 | << "gap nodes: " << gap_node_count << std::endl |
---|
608 | << std::endl; |
---|
609 | } |
---|
610 | |
---|
611 | void print_flow_values(std::ostream& os) const { |
---|
612 | os << "flow values" << std::endl; |
---|
613 | vertex_iterator u_iter, u_end; |
---|
614 | out_edge_iterator ei, e_end; |
---|
615 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) |
---|
616 | for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) |
---|
617 | if (capacity[*ei] > 0) |
---|
618 | os << *u_iter << " " << target(*ei, g) << " " |
---|
619 | << (capacity[*ei] - residual_capacity[*ei]) << std::endl; |
---|
620 | os << std::endl; |
---|
621 | } |
---|
622 | |
---|
623 | //======================================================================= |
---|
624 | |
---|
625 | Graph& g; |
---|
626 | vertices_size_type n; |
---|
627 | vertices_size_type nm; |
---|
628 | EdgeCapacityMap capacity; |
---|
629 | vertex_descriptor src; |
---|
630 | vertex_descriptor sink; |
---|
631 | VertexIndexMap index; |
---|
632 | |
---|
633 | // will need to use random_access_property_map with these |
---|
634 | std::vector< FlowValue > excess_flow; |
---|
635 | std::vector< out_edge_iterator > current; |
---|
636 | std::vector< distance_size_type > distance; |
---|
637 | std::vector< default_color_type > color; |
---|
638 | |
---|
639 | // Edge Property Maps that must be interior to the graph |
---|
640 | ReverseEdgeMap reverse_edge; |
---|
641 | ResidualCapacityEdgeMap residual_capacity; |
---|
642 | |
---|
643 | LayerArray layers; |
---|
644 | std::vector< list_iterator > layer_list_ptr; |
---|
645 | distance_size_type max_distance; // maximal distance |
---|
646 | distance_size_type max_active; // maximal distance with active node |
---|
647 | distance_size_type min_active; // minimal distance with active node |
---|
648 | boost::queue<vertex_descriptor> Q; |
---|
649 | |
---|
650 | // Statistics counters |
---|
651 | long push_count; |
---|
652 | long update_count; |
---|
653 | long relabel_count; |
---|
654 | long gap_count; |
---|
655 | long gap_node_count; |
---|
656 | |
---|
657 | inline double global_update_frequency() { return 0.5; } |
---|
658 | inline vertices_size_type alpha() { return 6; } |
---|
659 | inline long beta() { return 12; } |
---|
660 | |
---|
661 | long work_since_last_update; |
---|
662 | }; |
---|
663 | |
---|
664 | } // namespace detail |
---|
665 | |
---|
666 | template <class Graph, |
---|
667 | class CapacityEdgeMap, class ResidualCapacityEdgeMap, |
---|
668 | class ReverseEdgeMap, class VertexIndexMap> |
---|
669 | typename property_traits<CapacityEdgeMap>::value_type |
---|
670 | push_relabel_max_flow |
---|
671 | (Graph& g, |
---|
672 | typename graph_traits<Graph>::vertex_descriptor src, |
---|
673 | typename graph_traits<Graph>::vertex_descriptor sink, |
---|
674 | CapacityEdgeMap cap, ResidualCapacityEdgeMap res, |
---|
675 | ReverseEdgeMap rev, VertexIndexMap index_map) |
---|
676 | { |
---|
677 | typedef typename property_traits<CapacityEdgeMap>::value_type FlowValue; |
---|
678 | |
---|
679 | detail::push_relabel<Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, |
---|
680 | ReverseEdgeMap, VertexIndexMap, FlowValue> |
---|
681 | algo(g, cap, res, rev, src, sink, index_map); |
---|
682 | |
---|
683 | FlowValue flow = algo.maximum_preflow(); |
---|
684 | |
---|
685 | algo.convert_preflow_to_flow(); |
---|
686 | |
---|
687 | assert(algo.is_flow()); |
---|
688 | assert(algo.is_optimal()); |
---|
689 | |
---|
690 | return flow; |
---|
691 | } // push_relabel_max_flow() |
---|
692 | |
---|
693 | template <class Graph, class P, class T, class R> |
---|
694 | typename detail::edge_capacity_value<Graph, P, T, R>::type |
---|
695 | push_relabel_max_flow |
---|
696 | (Graph& g, |
---|
697 | typename graph_traits<Graph>::vertex_descriptor src, |
---|
698 | typename graph_traits<Graph>::vertex_descriptor sink, |
---|
699 | const bgl_named_params<P, T, R>& params) |
---|
700 | { |
---|
701 | return push_relabel_max_flow |
---|
702 | (g, src, sink, |
---|
703 | choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity), |
---|
704 | choose_pmap(get_param(params, edge_residual_capacity), |
---|
705 | g, edge_residual_capacity), |
---|
706 | choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse), |
---|
707 | choose_const_pmap(get_param(params, vertex_index), g, vertex_index) |
---|
708 | ); |
---|
709 | } |
---|
710 | |
---|
711 | template <class Graph> |
---|
712 | typename property_traits< |
---|
713 | typename property_map<Graph, edge_capacity_t>::const_type |
---|
714 | >::value_type |
---|
715 | push_relabel_max_flow |
---|
716 | (Graph& g, |
---|
717 | typename graph_traits<Graph>::vertex_descriptor src, |
---|
718 | typename graph_traits<Graph>::vertex_descriptor sink) |
---|
719 | { |
---|
720 | bgl_named_params<int, buffer_param_t> params(0); // bogus empty param |
---|
721 | return push_relabel_max_flow(g, src, sink, params); |
---|
722 | } |
---|
723 | |
---|
724 | } // namespace boost |
---|
725 | |
---|
726 | #endif // BOOST_PUSH_RELABEL_MAX_FLOW_HPP |
---|
727 | |
---|