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