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 <boost/assert.hpp> |
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 get(distance, u) == get(distance, v) + 1; |
79 | } |
80 | inline bool is_residual_edge(edge_descriptor a) { |
81 | return 0 < get(residual_capacity, a); |
82 | } |
83 | inline bool is_saturated(edge_descriptor a) { |
84 | return get(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(get(distance, u), max_active); |
97 | min_active = min BOOST_PREVENT_MACRO_SUBSTITUTION(get(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[get(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[get(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_data(num_vertices(g_)), |
124 | excess_flow(excess_flow_data.begin(), idx), |
125 | current_data(num_vertices(g_), out_edges(*vertices(g_).first, g_)), |
126 | current(current_data.begin(), idx), |
127 | distance_data(num_vertices(g_)), |
128 | distance(distance_data.begin(), idx), |
129 | color_data(num_vertices(g_)), |
130 | color(color_data.begin(), idx), |
131 | reverse_edge(rev), |
132 | residual_capacity(res), |
133 | layers(num_vertices(g_)), |
134 | layer_list_ptr_data(num_vertices(g_), |
135 | layers.front().inactive_vertices.end()), |
136 | layer_list_ptr(layer_list_ptr_data.begin(), idx), |
137 | push_count(0), update_count(0), relabel_count(0), |
138 | gap_count(0), gap_node_count(0), |
139 | work_since_last_update(0) |
140 | { |
141 | vertex_iterator u_iter, u_end; |
142 | // Don't count the reverse edges |
143 | edges_size_type m = num_edges(g) / 2; |
144 | nm = alpha() * n + m; |
145 | |
146 | // Initialize flow to zero which means initializing |
147 | // the residual capacity to equal the capacity. |
148 | out_edge_iterator ei, e_end; |
149 | for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) |
150 | for (boost::tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) { |
151 | put(residual_capacity, *ei, get(capacity, *ei)); |
152 | } |
153 | |
154 | for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
155 | vertex_descriptor u = *u_iter; |
156 | put(excess_flow, u, 0); |
157 | current[u] = out_edges(u, g); |
158 | } |
159 | |
160 | bool overflow_detected = false; |
161 | FlowValue test_excess = 0; |
162 | |
163 | out_edge_iterator a_iter, a_end; |
164 | for (boost::tie(a_iter, a_end) = out_edges(src, g); a_iter != a_end; ++a_iter) |
165 | if (target(*a_iter, g) != src) |
166 | test_excess += get(residual_capacity, *a_iter); |
167 | if (test_excess > (std::numeric_limits<FlowValue>::max)()) |
168 | overflow_detected = true; |
169 | |
170 | if (overflow_detected) |
171 | put(excess_flow, src, (std::numeric_limits<FlowValue>::max)()); |
172 | else { |
173 | put(excess_flow, src, 0); |
174 | for (boost::tie(a_iter, a_end) = out_edges(src, g); |
175 | a_iter != a_end; ++a_iter) { |
176 | edge_descriptor a = *a_iter; |
177 | vertex_descriptor tgt = target(a, g); |
178 | if (tgt != src) { |
179 | ++push_count; |
180 | FlowValue delta = get(residual_capacity, a); |
181 | put(residual_capacity, a, get(residual_capacity, a) - delta); |
182 | edge_descriptor rev = get(reverse_edge, a); |
183 | put(residual_capacity, rev, get(residual_capacity, rev) + delta); |
184 | put(excess_flow, tgt, get(excess_flow, tgt) + delta); |
185 | } |
186 | } |
187 | } |
188 | max_distance = num_vertices(g) - 1; |
189 | max_active = 0; |
190 | min_active = n; |
191 | |
192 | for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
193 | vertex_descriptor u = *u_iter; |
194 | if (u == sink) { |
195 | put(distance, u, 0); |
196 | continue; |
197 | } else if (u == src && !overflow_detected) |
198 | put(distance, u, n); |
199 | else |
200 | put(distance, u, 1); |
201 | |
202 | if (get(excess_flow, u) > 0) |
203 | add_to_active_list(u, layer&: layers[1]); |
204 | else if (get(distance, u) < n) |
205 | add_to_inactive_list(u, layer&: layers[1]); |
206 | } |
207 | |
208 | } // push_relabel constructor |
209 | |
210 | //======================================================================= |
211 | // This is a breadth-first search over the residual graph |
212 | // (well, actually the reverse of the residual graph). |
213 | // Would be cool to have a graph view adaptor for hiding certain |
214 | // edges, like the saturated (non-residual) edges in this case. |
215 | // Goldberg's implementation abused "distance" for the coloring. |
216 | void global_distance_update() |
217 | { |
218 | BOOST_USING_STD_MAX(); |
219 | ++update_count; |
220 | vertex_iterator u_iter, u_end; |
221 | for (boost::tie(u_iter,u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
222 | put(color, *u_iter, ColorTraits::white()); |
223 | put(distance, *u_iter, n); |
224 | } |
225 | put(color, sink, ColorTraits::gray()); |
226 | put(distance, sink, 0); |
227 | |
228 | for (distance_size_type l = 0; l <= max_distance; ++l) { |
229 | layers[l].active_vertices.clear(); |
230 | layers[l].inactive_vertices.clear(); |
231 | } |
232 | |
233 | max_distance = max_active = 0; |
234 | min_active = n; |
235 | |
236 | Q.push(sink); |
237 | while (! Q.empty()) { |
238 | vertex_descriptor u = Q.top(); |
239 | Q.pop(); |
240 | distance_size_type d_v = get(distance, u) + 1; |
241 | |
242 | out_edge_iterator ai, a_end; |
243 | for (boost::tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) { |
244 | edge_descriptor a = *ai; |
245 | vertex_descriptor v = target(a, g); |
246 | if (get(color, v) == ColorTraits::white() |
247 | && is_residual_edge(a: get(reverse_edge, a))) { |
248 | put(distance, v, d_v); |
249 | put(color, v, ColorTraits::gray()); |
250 | current[v] = out_edges(v, g); |
251 | max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(d_v, max_distance); |
252 | |
253 | if (get(excess_flow, v) > 0) |
254 | add_to_active_list(u: v, layer&: layers[d_v]); |
255 | else |
256 | add_to_inactive_list(u: v, layer&: layers[d_v]); |
257 | |
258 | Q.push(v); |
259 | } |
260 | } |
261 | } |
262 | } // global_distance_update() |
263 | |
264 | //======================================================================= |
265 | // This function is called "push" in Goldberg's h_prf implementation, |
266 | // but it is called "discharge" in the paper and in hi_pr.c. |
267 | void discharge(vertex_descriptor u) |
268 | { |
269 | BOOST_ASSERT(get(excess_flow, u) > 0); |
270 | while (1) { |
271 | out_edge_iterator ai, ai_end; |
272 | for (boost::tie(ai, ai_end) = current[u]; ai != ai_end; ++ai) { |
273 | edge_descriptor a = *ai; |
274 | if (is_residual_edge(a)) { |
275 | vertex_descriptor v = target(a, g); |
276 | if (is_admissible(u, v)) { |
277 | ++push_count; |
278 | if (v != sink && get(excess_flow, v) == 0) { |
279 | remove_from_inactive_list(u: v); |
280 | add_to_active_list(u: v, layer&: layers[get(distance, v)]); |
281 | } |
282 | push_flow(u_v: a); |
283 | if (get(excess_flow, u) == 0) |
284 | break; |
285 | } |
286 | } |
287 | } // for out_edges of i starting from current |
288 | |
289 | Layer& layer = layers[get(distance, u)]; |
290 | distance_size_type du = get(distance, u); |
291 | |
292 | if (ai == ai_end) { // i must be relabeled |
293 | relabel_distance(u); |
294 | if (layer.active_vertices.empty() |
295 | && layer.inactive_vertices.empty()) |
296 | gap(empty_distance: du); |
297 | if (get(distance, u) == n) |
298 | break; |
299 | } else { // i is no longer active |
300 | current[u].first = ai; |
301 | add_to_inactive_list(u, layer); |
302 | break; |
303 | } |
304 | } // while (1) |
305 | } // discharge() |
306 | |
307 | //======================================================================= |
308 | // This corresponds to the "push" update operation of the paper, |
309 | // not the "push" function in Goldberg's h_prf.c implementation. |
310 | // The idea is to push the excess flow from from vertex u to v. |
311 | void push_flow(edge_descriptor u_v) |
312 | { |
313 | vertex_descriptor |
314 | u = source(u_v, g), |
315 | v = target(u_v, g); |
316 | |
317 | BOOST_USING_STD_MIN(); |
318 | FlowValue flow_delta |
319 | = min BOOST_PREVENT_MACRO_SUBSTITUTION(get(excess_flow, u), get(residual_capacity, u_v)); |
320 | |
321 | put(residual_capacity, u_v, get(residual_capacity, u_v) - flow_delta); |
322 | edge_descriptor rev = get(reverse_edge, u_v); |
323 | put(residual_capacity, rev, get(residual_capacity, rev) + flow_delta); |
324 | |
325 | put(excess_flow, u, get(excess_flow, u) - flow_delta); |
326 | put(excess_flow, v, get(excess_flow, v) + flow_delta); |
327 | } // push_flow() |
328 | |
329 | //======================================================================= |
330 | // The main purpose of this routine is to set distance[v] |
331 | // to the smallest value allowed by the valid labeling constraints, |
332 | // which are: |
333 | // distance[t] = 0 |
334 | // distance[u] <= distance[v] + 1 for every residual edge (u,v) |
335 | // |
336 | distance_size_type relabel_distance(vertex_descriptor u) |
337 | { |
338 | BOOST_USING_STD_MAX(); |
339 | ++relabel_count; |
340 | work_since_last_update += beta(); |
341 | |
342 | distance_size_type min_distance = num_vertices(g); |
343 | put(distance, u, min_distance); |
344 | |
345 | // Examine the residual out-edges of vertex i, choosing the |
346 | // edge whose target vertex has the minimal distance. |
347 | out_edge_iterator ai, a_end, min_edge_iter; |
348 | for (boost::tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) { |
349 | ++work_since_last_update; |
350 | edge_descriptor a = *ai; |
351 | vertex_descriptor v = target(a, g); |
352 | if (is_residual_edge(a) && get(distance, v) < min_distance) { |
353 | min_distance = get(distance, v); |
354 | min_edge_iter = ai; |
355 | } |
356 | } |
357 | ++min_distance; |
358 | if (min_distance < n) { |
359 | put(distance, u, min_distance); // this is the main action |
360 | current[u].first = min_edge_iter; |
361 | max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(min_distance, max_distance); |
362 | } |
363 | return min_distance; |
364 | } // relabel_distance() |
365 | |
366 | //======================================================================= |
367 | // cleanup beyond the gap |
368 | void gap(distance_size_type empty_distance) |
369 | { |
370 | ++gap_count; |
371 | |
372 | distance_size_type r; // distance of layer before the current layer |
373 | r = empty_distance - 1; |
374 | |
375 | // Set the distance for the vertices beyond the gap to "infinity". |
376 | for (layer_iterator l = layers.begin() + empty_distance + 1; |
377 | l < layers.begin() + max_distance; ++l) { |
378 | list_iterator i; |
379 | for (i = l->inactive_vertices.begin(); |
380 | i != l->inactive_vertices.end(); ++i) { |
381 | put(distance, *i, n); |
382 | ++gap_node_count; |
383 | } |
384 | l->inactive_vertices.clear(); |
385 | } |
386 | max_distance = r; |
387 | max_active = r; |
388 | } |
389 | |
390 | //======================================================================= |
391 | // This is the core part of the algorithm, "phase one". |
392 | FlowValue maximum_preflow() |
393 | { |
394 | work_since_last_update = 0; |
395 | |
396 | while (max_active >= min_active) { // "main" loop |
397 | |
398 | Layer& layer = layers[max_active]; |
399 | list_iterator u_iter = layer.active_vertices.begin(); |
400 | |
401 | if (u_iter == layer.active_vertices.end()) |
402 | --max_active; |
403 | else { |
404 | vertex_descriptor u = *u_iter; |
405 | remove_from_active_list(u); |
406 | |
407 | discharge(u); |
408 | |
409 | if (work_since_last_update * global_update_frequency() > nm) { |
410 | global_distance_update(); |
411 | work_since_last_update = 0; |
412 | } |
413 | } |
414 | } // while (max_active >= min_active) |
415 | |
416 | return get(excess_flow, sink); |
417 | } // maximum_preflow() |
418 | |
419 | //======================================================================= |
420 | // remove excess flow, the "second phase" |
421 | // This does a DFS on the reverse flow graph of nodes with excess flow. |
422 | // If a cycle is found, cancel it. |
423 | // Return the nodes with excess flow in topological order. |
424 | // |
425 | // Unlike the prefl_to_flow() implementation, we use |
426 | // "color" instead of "distance" for the DFS labels |
427 | // "parent" instead of nl_prev for the DFS tree |
428 | // "topo_next" instead of nl_next for the topological ordering |
429 | void convert_preflow_to_flow() |
430 | { |
431 | vertex_iterator u_iter, u_end; |
432 | out_edge_iterator ai, a_end; |
433 | |
434 | vertex_descriptor r, restart, u; |
435 | |
436 | std::vector<vertex_descriptor> parent(n); |
437 | std::vector<vertex_descriptor> topo_next(n); |
438 | |
439 | vertex_descriptor tos(parent[0]), |
440 | bos(parent[0]); // bogus initialization, just to avoid warning |
441 | bool bos_null = true; |
442 | |
443 | // handle self-loops |
444 | for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) |
445 | for (boost::tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) |
446 | if (target(*ai, g) == *u_iter) |
447 | put(residual_capacity, *ai, get(capacity, *ai)); |
448 | |
449 | // initialize |
450 | for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
451 | u = *u_iter; |
452 | put(color, u, ColorTraits::white()); |
453 | parent[get(index, u)] = u; |
454 | current[u] = out_edges(u, g); |
455 | } |
456 | // eliminate flow cycles and topologically order the vertices |
457 | for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
458 | u = *u_iter; |
459 | if (get(color, u) == ColorTraits::white() |
460 | && get(excess_flow, u) > 0 |
461 | && u != src && u != sink ) { |
462 | r = u; |
463 | put(color, r, ColorTraits::gray()); |
464 | while (1) { |
465 | for (; current[u].first != current[u].second; ++current[u].first) { |
466 | edge_descriptor a = *current[u].first; |
467 | if (get(capacity, a) == 0 && is_residual_edge(a)) { |
468 | vertex_descriptor v = target(a, g); |
469 | if (get(color, v) == ColorTraits::white()) { |
470 | put(color, v, ColorTraits::gray()); |
471 | parent[get(index, v)] = u; |
472 | u = v; |
473 | break; |
474 | } else if (get(color, v) == ColorTraits::gray()) { |
475 | // find minimum flow on the cycle |
476 | FlowValue delta = get(residual_capacity, a); |
477 | while (1) { |
478 | BOOST_USING_STD_MIN(); |
479 | delta = min BOOST_PREVENT_MACRO_SUBSTITUTION(delta, get(residual_capacity, *current[v].first)); |
480 | if (v == u) |
481 | break; |
482 | else |
483 | v = target(*current[v].first, g); |
484 | } |
485 | // remove delta flow units |
486 | v = u; |
487 | while (1) { |
488 | a = *current[v].first; |
489 | put(residual_capacity, a, get(residual_capacity, a) - delta); |
490 | edge_descriptor rev = get(reverse_edge, a); |
491 | put(residual_capacity, rev, get(residual_capacity, rev) + delta); |
492 | v = target(a, g); |
493 | if (v == u) |
494 | break; |
495 | } |
496 | |
497 | // back-out of DFS to the first saturated edge |
498 | restart = u; |
499 | for (v = target(*current[u].first, g); v != u; v = target(a, g)){ |
500 | a = *current[v].first; |
501 | if (get(color, v) == ColorTraits::white() |
502 | || is_saturated(a)) { |
503 | put(color, target(*current[v].first, g), ColorTraits::white()); |
504 | if (get(color, v) != ColorTraits::white()) |
505 | restart = v; |
506 | } |
507 | } |
508 | if (restart != u) { |
509 | u = restart; |
510 | ++current[u].first; |
511 | break; |
512 | } |
513 | } // else if (color[v] == ColorTraits::gray()) |
514 | } // if (get(capacity, a) == 0 ... |
515 | } // for out_edges(u, g) (though "u" changes during loop) |
516 | |
517 | if ( current[u].first == current[u].second ) { |
518 | // scan of i is complete |
519 | put(color, u, ColorTraits::black()); |
520 | if (u != src) { |
521 | if (bos_null) { |
522 | bos = u; |
523 | bos_null = false; |
524 | tos = u; |
525 | } else { |
526 | topo_next[get(index, u)] = tos; |
527 | tos = u; |
528 | } |
529 | } |
530 | if (u != r) { |
531 | u = parent[get(index, u)]; |
532 | ++current[u].first; |
533 | } else |
534 | break; |
535 | } |
536 | } // while (1) |
537 | } // if (color[u] == white && excess_flow[u] > 0 & ...) |
538 | } // for all vertices in g |
539 | |
540 | // return excess flows |
541 | // note that the sink is not on the stack |
542 | if (! bos_null) { |
543 | for (u = tos; u != bos; u = topo_next[get(index, u)]) { |
544 | boost::tie(ai, a_end) = out_edges(u, g); |
545 | while (get(excess_flow, u) > 0 && ai != a_end) { |
546 | if (get(capacity, *ai) == 0 && is_residual_edge(a: *ai)) |
547 | push_flow(u_v: *ai); |
548 | ++ai; |
549 | } |
550 | } |
551 | // do the bottom |
552 | u = bos; |
553 | boost::tie(ai, a_end) = out_edges(u, g); |
554 | while (get(excess_flow, u) > 0 && ai != a_end) { |
555 | if (get(capacity, *ai) == 0 && is_residual_edge(a: *ai)) |
556 | push_flow(u_v: *ai); |
557 | ++ai; |
558 | } |
559 | } |
560 | |
561 | } // convert_preflow_to_flow() |
562 | |
563 | //======================================================================= |
564 | inline bool is_flow() |
565 | { |
566 | vertex_iterator u_iter, u_end; |
567 | out_edge_iterator ai, a_end; |
568 | |
569 | // check edge flow values |
570 | for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
571 | for (boost::tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) { |
572 | edge_descriptor a = *ai; |
573 | if (get(capacity, a) > 0) |
574 | if ((get(residual_capacity, a) + get(residual_capacity, get(reverse_edge, a)) |
575 | != get(capacity, a) + get(capacity, get(reverse_edge, a))) |
576 | || (get(residual_capacity, a) < 0) |
577 | || (get(residual_capacity, get(reverse_edge, a)) < 0)) |
578 | return false; |
579 | } |
580 | } |
581 | |
582 | // check conservation |
583 | FlowValue sum; |
584 | for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { |
585 | vertex_descriptor u = *u_iter; |
586 | if (u != src && u != sink) { |
587 | if (get(excess_flow, u) != 0) |
588 | return false; |
589 | sum = 0; |
590 | for (boost::tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) |
591 | if (get(capacity, *ai) > 0) |
592 | sum -= get(capacity, *ai) - get(residual_capacity, *ai); |
593 | else |
594 | sum += get(residual_capacity, *ai); |
595 | |
596 | if (get(excess_flow, u) != sum) |
597 | return false; |
598 | } |
599 | } |
600 | |
601 | return true; |
602 | } // is_flow() |
603 | |
604 | bool is_optimal() { |
605 | // check if mincut is saturated... |
606 | global_distance_update(); |
607 | return get(distance, src) >= n; |
608 | } |
609 | |
610 | void print_statistics(std::ostream& os) const { |
611 | os << "pushes: " << push_count << std::endl |
612 | << "relabels: " << relabel_count << std::endl |
613 | << "updates: " << update_count << std::endl |
614 | << "gaps: " << gap_count << std::endl |
615 | << "gap nodes: " << gap_node_count << std::endl |
616 | << std::endl; |
617 | } |
618 | |
619 | void print_flow_values(std::ostream& os) const { |
620 | os << "flow values" << std::endl; |
621 | vertex_iterator u_iter, u_end; |
622 | out_edge_iterator ei, e_end; |
623 | for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) |
624 | for (boost::tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) |
625 | if (get(capacity, *ei) > 0) |
626 | os << *u_iter << " " << target(*ei, g) << " " |
627 | << (get(capacity, *ei) - get(residual_capacity, *ei)) << std::endl; |
628 | os << std::endl; |
629 | } |
630 | |
631 | //======================================================================= |
632 | |
633 | Graph& g; |
634 | vertices_size_type n; |
635 | vertices_size_type nm; |
636 | EdgeCapacityMap capacity; |
637 | vertex_descriptor src; |
638 | vertex_descriptor sink; |
639 | VertexIndexMap index; |
640 | |
641 | // will need to use random_access_property_map with these |
642 | std::vector< FlowValue > excess_flow_data; |
643 | iterator_property_map<typename std::vector<FlowValue>::iterator, VertexIndexMap> excess_flow; |
644 | std::vector< std::pair<out_edge_iterator, out_edge_iterator> > current_data; |
645 | iterator_property_map< |
646 | typename std::vector< std::pair<out_edge_iterator, out_edge_iterator> >::iterator, |
647 | VertexIndexMap> current; |
648 | std::vector< distance_size_type > distance_data; |
649 | iterator_property_map< |
650 | typename std::vector< distance_size_type >::iterator, |
651 | VertexIndexMap> distance; |
652 | std::vector< default_color_type > color_data; |
653 | iterator_property_map< |
654 | std::vector< default_color_type >::iterator, |
655 | VertexIndexMap> color; |
656 | |
657 | // Edge Property Maps that must be interior to the graph |
658 | ReverseEdgeMap reverse_edge; |
659 | ResidualCapacityEdgeMap residual_capacity; |
660 | |
661 | LayerArray layers; |
662 | std::vector< list_iterator > layer_list_ptr_data; |
663 | iterator_property_map<typename std::vector< list_iterator >::iterator, VertexIndexMap> layer_list_ptr; |
664 | distance_size_type max_distance; // maximal distance |
665 | distance_size_type max_active; // maximal distance with active node |
666 | distance_size_type min_active; // minimal distance with active node |
667 | boost::queue<vertex_descriptor> Q; |
668 | |
669 | // Statistics counters |
670 | long push_count; |
671 | long update_count; |
672 | long relabel_count; |
673 | long gap_count; |
674 | long gap_node_count; |
675 | |
676 | inline double global_update_frequency() { return 0.5; } |
677 | inline vertices_size_type alpha() { return 6; } |
678 | inline long beta() { return 12; } |
679 | |
680 | long work_since_last_update; |
681 | }; |
682 | |
683 | } // namespace detail |
684 | |
685 | template <class Graph, |
686 | class CapacityEdgeMap, class ResidualCapacityEdgeMap, |
687 | class ReverseEdgeMap, class VertexIndexMap> |
688 | typename property_traits<CapacityEdgeMap>::value_type |
689 | push_relabel_max_flow |
690 | (Graph& g, |
691 | typename graph_traits<Graph>::vertex_descriptor src, |
692 | typename graph_traits<Graph>::vertex_descriptor sink, |
693 | CapacityEdgeMap cap, ResidualCapacityEdgeMap res, |
694 | ReverseEdgeMap rev, VertexIndexMap index_map) |
695 | { |
696 | typedef typename property_traits<CapacityEdgeMap>::value_type FlowValue; |
697 | |
698 | detail::push_relabel<Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, |
699 | ReverseEdgeMap, VertexIndexMap, FlowValue> |
700 | algo(g, cap, res, rev, src, sink, index_map); |
701 | |
702 | FlowValue flow = algo.maximum_preflow(); |
703 | |
704 | algo.convert_preflow_to_flow(); |
705 | |
706 | BOOST_ASSERT(algo.is_flow()); |
707 | BOOST_ASSERT(algo.is_optimal()); |
708 | |
709 | return flow; |
710 | } // push_relabel_max_flow() |
711 | |
712 | template <class Graph, class P, class T, class R> |
713 | typename detail::edge_capacity_value<Graph, P, T, R>::type |
714 | push_relabel_max_flow |
715 | (Graph& g, |
716 | typename graph_traits<Graph>::vertex_descriptor src, |
717 | typename graph_traits<Graph>::vertex_descriptor sink, |
718 | const bgl_named_params<P, T, R>& params) |
719 | { |
720 | return push_relabel_max_flow |
721 | (g, src, sink, |
722 | choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity), |
723 | choose_pmap(get_param(params, edge_residual_capacity), |
724 | g, edge_residual_capacity), |
725 | choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse), |
726 | choose_const_pmap(get_param(params, vertex_index), g, vertex_index) |
727 | ); |
728 | } |
729 | |
730 | template <class Graph> |
731 | typename property_traits< |
732 | typename property_map<Graph, edge_capacity_t>::const_type |
733 | >::value_type |
734 | push_relabel_max_flow |
735 | (Graph& g, |
736 | typename graph_traits<Graph>::vertex_descriptor src, |
737 | typename graph_traits<Graph>::vertex_descriptor sink) |
738 | { |
739 | bgl_named_params<int, buffer_param_t> params(0); // bogus empty param |
740 | return push_relabel_max_flow(g, src, sink, params); |
741 | } |
742 | |
743 | } // namespace boost |
744 | |
745 | #endif // BOOST_PUSH_RELABEL_MAX_FLOW_HPP |
746 | |
747 | |