1 //=======================================================================
2 // Copyright 2000 University of Notre Dame.
3 // Authors: Jeremy G. Siek, Andrew Lumsdaine, Lie-Quan Lee
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 //=======================================================================
10 #ifndef BOOST_PUSH_RELABEL_MAX_FLOW_HPP
11 #define BOOST_PUSH_RELABEL_MAX_FLOW_HPP
13 #include <boost/config.hpp>
14 #include <boost/assert.hpp>
18 #include <algorithm> // for std::min and std::max
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>
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.
34 // This implements the highest-label version of the push-relabel method
35 // with the global relabeling and gap relabeling heuristics.
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
44 template <class Vertex>
45 struct preflow_layer {
46 std::list<Vertex> active_vertices;
47 std::list<Vertex> inactive_vertices;
50 template <class Graph,
51 class EdgeCapacityMap, // integer value type
52 class ResidualCapacityEdgeMap,
54 class VertexIndexMap, // vertex_descriptor -> integer
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;
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;
72 typedef color_traits<default_color_type> ColorTraits;
74 //=======================================================================
75 // Some helper predicates
77 inline bool is_admissible(vertex_descriptor u, vertex_descriptor v) {
78 return get(distance, u) == get(distance, v) + 1;
80 inline bool is_residual_edge(edge_descriptor a) {
81 return 0 < get(residual_capacity, a);
83 inline bool is_saturated(edge_descriptor a) {
84 return get(residual_capacity, a) == 0;
87 //=======================================================================
88 // Layer List Management Functions
90 typedef typename std::list<vertex_descriptor>::iterator list_iterator;
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();
100 void remove_from_active_list(vertex_descriptor u) {
101 layers[get(distance, u)].active_vertices.erase(layer_list_ptr[u]);
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();
108 void remove_from_inactive_list(vertex_descriptor u) {
109 layers[get(distance, u)].inactive_vertices.erase(layer_list_ptr[u]);
112 //=======================================================================
114 push_relabel(Graph& g_,
116 ResidualCapacityEdgeMap res,
118 vertex_descriptor src_,
119 vertex_descriptor sink_,
121 : g(g_), n(num_vertices(g_)), capacity(cap), src(src_), sink(sink_),
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),
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)
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;
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));
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);
160 bool overflow_detected = false;
161 FlowValue test_excess = 0;
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;
170 if (overflow_detected)
171 put(excess_flow, src, (std::numeric_limits<FlowValue>::max)());
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);
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);
188 max_distance = num_vertices(g) - 1;
192 for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
193 vertex_descriptor u = *u_iter;
197 } else if (u == src && !overflow_detected)
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]);
208 } // push_relabel constructor
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()
218 BOOST_USING_STD_MAX();
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);
225 put(color, sink, ColorTraits::gray());
226 put(distance, sink, 0);
228 for (distance_size_type l = 0; l <= max_distance; ++l) {
229 layers[l].active_vertices.clear();
230 layers[l].inactive_vertices.clear();
233 max_distance = max_active = 0;
237 while (! Q.empty()) {
238 vertex_descriptor u = Q.top();
240 distance_size_type d_v = get(distance, u) + 1;
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);
253 if (get(excess_flow, v) > 0)
254 add_to_active_list(v, layers[d_v]);
256 add_to_inactive_list(v, layers[d_v]);
262 } // global_distance_update()
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)
269 BOOST_ASSERT(get(excess_flow, u) > 0);
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)) {
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)]);
283 if (get(excess_flow, u) == 0)
287 } // for out_edges of i starting from current
289 Layer& layer = layers[get(distance, u)];
290 distance_size_type du = get(distance, u);
292 if (ai == ai_end) { // i must be relabeled
294 if (layer.active_vertices.empty()
295 && layer.inactive_vertices.empty())
297 if (get(distance, u) == n)
299 } else { // i is no longer active
300 current[u].first = ai;
301 add_to_inactive_list(u, layer);
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)
317 BOOST_USING_STD_MIN();
319 = min BOOST_PREVENT_MACRO_SUBSTITUTION(get(excess_flow, u), get(residual_capacity, u_v));
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);
325 put(excess_flow, u, get(excess_flow, u) - flow_delta);
326 put(excess_flow, v, get(excess_flow, v) + flow_delta);
329 //=======================================================================
330 // The main purpose of this routine is to set distance[v]
331 // to the smallest value allowed by the valid labeling constraints,
334 // distance[u] <= distance[v] + 1 for every residual edge (u,v)
336 distance_size_type relabel_distance(vertex_descriptor u)
338 BOOST_USING_STD_MAX();
340 work_since_last_update += beta();
342 distance_size_type min_distance = num_vertices(g);
343 put(distance, u, min_distance);
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);
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);
364 } // relabel_distance()
366 //=======================================================================
367 // cleanup beyond the gap
368 void gap(distance_size_type empty_distance)
372 distance_size_type r; // distance of layer before the current layer
373 r = empty_distance - 1;
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) {
379 for (i = l->inactive_vertices.begin();
380 i != l->inactive_vertices.end(); ++i) {
381 put(distance, *i, n);
384 l->inactive_vertices.clear();
390 //=======================================================================
391 // This is the core part of the algorithm, "phase one".
392 FlowValue maximum_preflow()
394 work_since_last_update = 0;
396 while (max_active >= min_active) { // "main" loop
398 Layer& layer = layers[max_active];
399 list_iterator u_iter = layer.active_vertices.begin();
401 if (u_iter == layer.active_vertices.end())
404 vertex_descriptor u = *u_iter;
405 remove_from_active_list(u);
409 if (work_since_last_update * global_update_frequency() > nm) {
410 global_distance_update();
411 work_since_last_update = 0;
414 } // while (max_active >= min_active)
416 return get(excess_flow, sink);
417 } // maximum_preflow()
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.
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()
431 vertex_iterator u_iter, u_end;
432 out_edge_iterator ai, a_end;
434 vertex_descriptor r, restart, u;
436 std::vector<vertex_descriptor> parent(n);
437 std::vector<vertex_descriptor> topo_next(n);
439 vertex_descriptor tos(parent[0]),
440 bos(parent[0]); // bogus initialization, just to avoid warning
441 bool bos_null = true;
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));
450 for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
452 put(color, u, ColorTraits::white());
453 parent[get(index, u)] = u;
454 current[u] = out_edges(u, g);
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) {
459 if (get(color, u) == ColorTraits::white()
460 && get(excess_flow, u) > 0
461 && u != src && u != sink ) {
463 put(color, r, ColorTraits::gray());
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;
474 } else if (get(color, v) == ColorTraits::gray()) {
475 // find minimum flow on the cycle
476 FlowValue delta = get(residual_capacity, a);
478 BOOST_USING_STD_MIN();
479 delta = min BOOST_PREVENT_MACRO_SUBSTITUTION(delta, get(residual_capacity, *current[v].first));
483 v = target(*current[v].first, g);
485 // remove delta flow units
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);
497 // back-out of DFS to the first saturated edge
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())
513 } // else if (color[v] == ColorTraits::gray())
514 } // if (get(capacity, a) == 0 ...
515 } // for out_edges(u, g) (though "u" changes during loop)
517 if ( current[u].first == current[u].second ) {
518 // scan of i is complete
519 put(color, u, ColorTraits::black());
526 topo_next[get(index, u)] = tos;
531 u = parent[get(index, u)];
537 } // if (color[u] == white && excess_flow[u] > 0 & ...)
538 } // for all vertices in g
540 // return excess flows
541 // note that the sink is not on the stack
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))
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))
561 } // convert_preflow_to_flow()
563 //=======================================================================
564 inline bool is_flow()
566 vertex_iterator u_iter, u_end;
567 out_edge_iterator ai, a_end;
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))
582 // check conservation
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)
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);
594 sum += get(residual_capacity, *ai);
596 if (get(excess_flow, u) != sum)
605 // check if mincut is saturated...
606 global_distance_update();
607 return get(distance, src) >= n;
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
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;
631 //=======================================================================
634 vertices_size_type n;
635 vertices_size_type nm;
636 EdgeCapacityMap capacity;
637 vertex_descriptor src;
638 vertex_descriptor sink;
639 VertexIndexMap index;
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;
657 // Edge Property Maps that must be interior to the graph
658 ReverseEdgeMap reverse_edge;
659 ResidualCapacityEdgeMap residual_capacity;
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;
669 // Statistics counters
676 inline double global_update_frequency() { return 0.5; }
677 inline vertices_size_type alpha() { return 6; }
678 inline long beta() { return 12; }
680 long work_since_last_update;
683 } // namespace detail
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
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)
696 typedef typename property_traits<CapacityEdgeMap>::value_type FlowValue;
698 detail::push_relabel<Graph, CapacityEdgeMap, ResidualCapacityEdgeMap,
699 ReverseEdgeMap, VertexIndexMap, FlowValue>
700 algo(g, cap, res, rev, src, sink, index_map);
702 FlowValue flow = algo.maximum_preflow();
704 algo.convert_preflow_to_flow();
706 BOOST_ASSERT(algo.is_flow());
707 BOOST_ASSERT(algo.is_optimal());
710 } // push_relabel_max_flow()
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
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)
720 return push_relabel_max_flow
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)
730 template <class Graph>
731 typename property_traits<
732 typename property_map<Graph, edge_capacity_t>::const_type
734 push_relabel_max_flow
736 typename graph_traits<Graph>::vertex_descriptor src,
737 typename graph_traits<Graph>::vertex_descriptor sink)
739 bgl_named_params<int, buffer_param_t> params(0); // bogus empty param
740 return push_relabel_max_flow(g, src, sink, params);
745 #endif // BOOST_PUSH_RELABEL_MAX_FLOW_HPP