]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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, layers[1]); | |
204 | else if (get(distance, u) < n) | |
205 | add_to_inactive_list(u, 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(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(v, layers[d_v]); | |
255 | else | |
256 | add_to_inactive_list(v, 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(v); | |
280 | add_to_active_list(v, layers[get(distance, v)]); | |
281 | } | |
282 | push_flow(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(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(*ai)) | |
547 | push_flow(*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(*ai)) | |
556 | push_flow(*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 |