boost/graph/push_relabel_max_flow.hpp
//======================================================================= // Copyright 2000 University of Notre Dame. // Authors: Jeremy G. Siek, Andrew Lumsdaine, Lie-Quan Lee // // Distributed under the Boost Software License, Version 1.0. (See // accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) //======================================================================= #ifndef BOOST_PUSH_RELABEL_MAX_FLOW_HPP #define BOOST_PUSH_RELABEL_MAX_FLOW_HPP #include <boost/config.hpp> #include <cassert> #include <vector> #include <list> #include <iosfwd> #include <algorithm> // for std::min and std::max #include <boost/pending/queue.hpp> #include <boost/limits.hpp> #include <boost/graph/graph_concepts.hpp> #include <boost/graph/named_function_params.hpp> namespace boost { namespace detail { // This implementation is based on Goldberg's // "On Implementing Push-Relabel Method for the Maximum Flow Problem" // by B.V. Cherkassky and A.V. Goldberg, IPCO '95, pp. 157--171 // and on the h_prf.c and hi_pr.c code written by the above authors. // This implements the highest-label version of the push-relabel method // with the global relabeling and gap relabeling heuristics. // The terms "rank", "distance", "height" are synonyms in // Goldberg's implementation, paper and in the CLR. A "layer" is a // group of vertices with the same distance. The vertices in each // layer are categorized as active or inactive. An active vertex // has positive excess flow and its distance is less than n (it is // not blocked). template <class Vertex> struct preflow_layer { std::list<Vertex> active_vertices; std::list<Vertex> inactive_vertices; }; template <class Graph, class EdgeCapacityMap, // integer value type class ResidualCapacityEdgeMap, class ReverseEdgeMap, class VertexIndexMap, // vertex_descriptor -> integer class FlowValue> class push_relabel { public: typedef graph_traits<Graph> Traits; typedef typename Traits::vertex_descriptor vertex_descriptor; typedef typename Traits::edge_descriptor edge_descriptor; typedef typename Traits::vertex_iterator vertex_iterator; typedef typename Traits::out_edge_iterator out_edge_iterator; typedef typename Traits::vertices_size_type vertices_size_type; typedef typename Traits::edges_size_type edges_size_type; typedef preflow_layer<vertex_descriptor> Layer; typedef std::vector< Layer > LayerArray; typedef typename LayerArray::iterator layer_iterator; typedef typename LayerArray::size_type distance_size_type; typedef color_traits<default_color_type> ColorTraits; //======================================================================= // Some helper predicates inline bool is_admissible(vertex_descriptor u, vertex_descriptor v) { return distance[u] == distance[v] + 1; } inline bool is_residual_edge(edge_descriptor a) { return 0 < residual_capacity[a]; } inline bool is_saturated(edge_descriptor a) { return residual_capacity[a] == 0; } //======================================================================= // Layer List Management Functions typedef typename std::list<vertex_descriptor>::iterator list_iterator; void add_to_active_list(vertex_descriptor u, Layer& layer) { BOOST_USING_STD_MIN(); BOOST_USING_STD_MAX(); layer.active_vertices.push_front(u); max_active = max BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], max_active); min_active = min BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], min_active); layer_list_ptr[u] = layer.active_vertices.begin(); } void remove_from_active_list(vertex_descriptor u) { layers[distance[u]].active_vertices.erase(layer_list_ptr[u]); } void add_to_inactive_list(vertex_descriptor u, Layer& layer) { layer.inactive_vertices.push_front(u); layer_list_ptr[u] = layer.inactive_vertices.begin(); } void remove_from_inactive_list(vertex_descriptor u) { layers[distance[u]].inactive_vertices.erase(layer_list_ptr[u]); } //======================================================================= // initialization push_relabel(Graph& g_, EdgeCapacityMap cap, ResidualCapacityEdgeMap res, ReverseEdgeMap rev, vertex_descriptor src_, vertex_descriptor sink_, VertexIndexMap idx) : g(g_), n(num_vertices(g_)), capacity(cap), src(src_), sink(sink_), index(idx), excess_flow(num_vertices(g_)), current(num_vertices(g_), out_edges(*vertices(g_).first, g_).second), distance(num_vertices(g_)), color(num_vertices(g_)), reverse_edge(rev), residual_capacity(res), layers(num_vertices(g_)), layer_list_ptr(num_vertices(g_), layers.front().inactive_vertices.end()), push_count(0), update_count(0), relabel_count(0), gap_count(0), gap_node_count(0), work_since_last_update(0) { vertex_iterator u_iter, u_end; // Don't count the reverse edges edges_size_type m = num_edges(g) / 2; nm = alpha() * n + m; // Initialize flow to zero which means initializing // the residual capacity to equal the capacity. out_edge_iterator ei, e_end; for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) { residual_capacity[*ei] = capacity[*ei]; } for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { vertex_descriptor u = *u_iter; excess_flow[u] = 0; current[u] = out_edges(u, g).first; } bool overflow_detected = false; FlowValue test_excess = 0; out_edge_iterator a_iter, a_end; for (tie(a_iter, a_end) = out_edges(src, g); a_iter != a_end; ++a_iter) if (target(*a_iter, g) != src) test_excess += residual_capacity[*a_iter]; if (test_excess > (std::numeric_limits<FlowValue>::max)()) overflow_detected = true; if (overflow_detected) excess_flow[src] = (std::numeric_limits<FlowValue>::max)(); else { excess_flow[src] = 0; for (tie(a_iter, a_end) = out_edges(src, g); a_iter != a_end; ++a_iter) { edge_descriptor a = *a_iter; if (target(a, g) != src) { ++push_count; FlowValue delta = residual_capacity[a]; residual_capacity[a] -= delta; residual_capacity[reverse_edge[a]] += delta; excess_flow[target(a, g)] += delta; } } } max_distance = num_vertices(g) - 1; max_active = 0; min_active = n; for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { vertex_descriptor u = *u_iter; if (u == sink) { distance[u] = 0; continue; } else if (u == src && !overflow_detected) distance[u] = n; else distance[u] = 1; if (excess_flow[u] > 0) add_to_active_list(u, layers[1]); else if (distance[u] < n) add_to_inactive_list(u, layers[1]); } } // push_relabel constructor //======================================================================= // This is a breadth-first search over the residual graph // (well, actually the reverse of the residual graph). // Would be cool to have a graph view adaptor for hiding certain // edges, like the saturated (non-residual) edges in this case. // Goldberg's implementation abused "distance" for the coloring. void global_distance_update() { BOOST_USING_STD_MAX(); ++update_count; vertex_iterator u_iter, u_end; for (tie(u_iter,u_end) = vertices(g); u_iter != u_end; ++u_iter) { color[*u_iter] = ColorTraits::white(); distance[*u_iter] = n; } color[sink] = ColorTraits::gray(); distance[sink] = 0; for (distance_size_type l = 0; l <= max_distance; ++l) { layers[l].active_vertices.clear(); layers[l].inactive_vertices.clear(); } max_distance = max_active = 0; min_active = n; Q.push(sink); while (! Q.empty()) { vertex_descriptor u = Q.top(); Q.pop(); distance_size_type d_v = distance[u] + 1; out_edge_iterator ai, a_end; for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) { edge_descriptor a = *ai; vertex_descriptor v = target(a, g); if (color[v] == ColorTraits::white() && is_residual_edge(reverse_edge[a])) { distance[v] = d_v; color[v] = ColorTraits::gray(); current[v] = out_edges(v, g).first; max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(d_v, max_distance); if (excess_flow[v] > 0) add_to_active_list(v, layers[d_v]); else add_to_inactive_list(v, layers[d_v]); Q.push(v); } } } } // global_distance_update() //======================================================================= // This function is called "push" in Goldberg's h_prf implementation, // but it is called "discharge" in the paper and in hi_pr.c. void discharge(vertex_descriptor u) { assert(excess_flow[u] > 0); while (1) { out_edge_iterator ai, ai_end; for (ai = current[u], ai_end = out_edges(u, g).second; ai != ai_end; ++ai) { edge_descriptor a = *ai; if (is_residual_edge(a)) { vertex_descriptor v = target(a, g); if (is_admissible(u, v)) { ++push_count; if (v != sink && excess_flow[v] == 0) { remove_from_inactive_list(v); add_to_active_list(v, layers[distance[v]]); } push_flow(a); if (excess_flow[u] == 0) break; } } } // for out_edges of i starting from current Layer& layer = layers[distance[u]]; distance_size_type du = distance[u]; if (ai == ai_end) { // i must be relabeled relabel_distance(u); if (layer.active_vertices.empty() && layer.inactive_vertices.empty()) gap(du); if (distance[u] == n) break; } else { // i is no longer active current[u] = ai; add_to_inactive_list(u, layer); break; } } // while (1) } // discharge() //======================================================================= // This corresponds to the "push" update operation of the paper, // not the "push" function in Goldberg's h_prf.c implementation. // The idea is to push the excess flow from from vertex u to v. void push_flow(edge_descriptor u_v) { vertex_descriptor u = source(u_v, g), v = target(u_v, g); BOOST_USING_STD_MIN(); FlowValue flow_delta = min BOOST_PREVENT_MACRO_SUBSTITUTION(excess_flow[u], residual_capacity[u_v]); residual_capacity[u_v] -= flow_delta; residual_capacity[reverse_edge[u_v]] += flow_delta; excess_flow[u] -= flow_delta; excess_flow[v] += flow_delta; } // push_flow() //======================================================================= // The main purpose of this routine is to set distance[v] // to the smallest value allowed by the valid labeling constraints, // which are: // distance[t] = 0 // distance[u] <= distance[v] + 1 for every residual edge (u,v) // distance_size_type relabel_distance(vertex_descriptor u) { BOOST_USING_STD_MAX(); ++relabel_count; work_since_last_update += beta(); distance_size_type min_distance = num_vertices(g); distance[u] = min_distance; // Examine the residual out-edges of vertex i, choosing the // edge whose target vertex has the minimal distance. out_edge_iterator ai, a_end, min_edge_iter; for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) { ++work_since_last_update; edge_descriptor a = *ai; vertex_descriptor v = target(a, g); if (is_residual_edge(a) && distance[v] < min_distance) { min_distance = distance[v]; min_edge_iter = ai; } } ++min_distance; if (min_distance < n) { distance[u] = min_distance; // this is the main action current[u] = min_edge_iter; max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(min_distance, max_distance); } return min_distance; } // relabel_distance() //======================================================================= // cleanup beyond the gap void gap(distance_size_type empty_distance) { ++gap_count; distance_size_type r; // distance of layer before the current layer r = empty_distance - 1; // Set the distance for the vertices beyond the gap to "infinity". for (layer_iterator l = layers.begin() + empty_distance + 1; l < layers.begin() + max_distance; ++l) { list_iterator i; for (i = l->inactive_vertices.begin(); i != l->inactive_vertices.end(); ++i) { distance[*i] = n; ++gap_node_count; } l->inactive_vertices.clear(); } max_distance = r; max_active = r; } //======================================================================= // This is the core part of the algorithm, "phase one". FlowValue maximum_preflow() { work_since_last_update = 0; while (max_active >= min_active) { // "main" loop Layer& layer = layers[max_active]; list_iterator u_iter = layer.active_vertices.begin(); if (u_iter == layer.active_vertices.end()) --max_active; else { vertex_descriptor u = *u_iter; remove_from_active_list(u); discharge(u); if (work_since_last_update * global_update_frequency() > nm) { global_distance_update(); work_since_last_update = 0; } } } // while (max_active >= min_active) return excess_flow[sink]; } // maximum_preflow() //======================================================================= // remove excess flow, the "second phase" // This does a DFS on the reverse flow graph of nodes with excess flow. // If a cycle is found, cancel it. // Return the nodes with excess flow in topological order. // // Unlike the prefl_to_flow() implementation, we use // "color" instead of "distance" for the DFS labels // "parent" instead of nl_prev for the DFS tree // "topo_next" instead of nl_next for the topological ordering void convert_preflow_to_flow() { vertex_iterator u_iter, u_end; out_edge_iterator ai, a_end; vertex_descriptor r, restart, u; std::vector<vertex_descriptor> parent(n); std::vector<vertex_descriptor> topo_next(n); vertex_descriptor tos(parent[0]), bos(parent[0]); // bogus initialization, just to avoid warning bool bos_null = true; // handle self-loops for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) if (target(*ai, g) == *u_iter) residual_capacity[*ai] = capacity[*ai]; // initialize for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { u = *u_iter; color[u] = ColorTraits::white(); parent[u] = u; current[u] = out_edges(u, g).first; } // eliminate flow cycles and topologically order the vertices for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { u = *u_iter; if (color[u] == ColorTraits::white() && excess_flow[u] > 0 && u != src && u != sink ) { r = u; color[r] = ColorTraits::gray(); while (1) { for (; current[u] != out_edges(u, g).second; ++current[u]) { edge_descriptor a = *current[u]; if (capacity[a] == 0 && is_residual_edge(a)) { vertex_descriptor v = target(a, g); if (color[v] == ColorTraits::white()) { color[v] = ColorTraits::gray(); parent[v] = u; u = v; break; } else if (color[v] == ColorTraits::gray()) { // find minimum flow on the cycle FlowValue delta = residual_capacity[a]; while (1) { BOOST_USING_STD_MIN(); delta = min BOOST_PREVENT_MACRO_SUBSTITUTION(delta, residual_capacity[*current[v]]); if (v == u) break; else v = target(*current[v], g); } // remove delta flow units v = u; while (1) { a = *current[v]; residual_capacity[a] -= delta; residual_capacity[reverse_edge[a]] += delta; v = target(a, g); if (v == u) break; } // back-out of DFS to the first saturated edge restart = u; for (v = target(*current[u], g); v != u; v = target(a, g)){ a = *current[v]; if (color[v] == ColorTraits::white() || is_saturated(a)) { color[target(*current[v], g)] = ColorTraits::white(); if (color[v] != ColorTraits::white()) restart = v; } } if (restart != u) { u = restart; ++current[u]; break; } } // else if (color[v] == ColorTraits::gray()) } // if (capacity[a] == 0 ... } // for out_edges(u, g) (though "u" changes during loop) if (current[u] == out_edges(u, g).second) { // scan of i is complete color[u] = ColorTraits::black(); if (u != src) { if (bos_null) { bos = u; bos_null = false; tos = u; } else { topo_next[u] = tos; tos = u; } } if (u != r) { u = parent[u]; ++current[u]; } else break; } } // while (1) } // if (color[u] == white && excess_flow[u] > 0 & ...) } // for all vertices in g // return excess flows // note that the sink is not on the stack if (! bos_null) { for (u = tos; u != bos; u = topo_next[u]) { ai = out_edges(u, g).first; while (excess_flow[u] > 0 && ai != out_edges(u, g).second) { if (capacity[*ai] == 0 && is_residual_edge(*ai)) push_flow(*ai); ++ai; } } // do the bottom u = bos; ai = out_edges(u, g).first; while (excess_flow[u] > 0) { if (capacity[*ai] == 0 && is_residual_edge(*ai)) push_flow(*ai); ++ai; } } } // convert_preflow_to_flow() //======================================================================= inline bool is_flow() { vertex_iterator u_iter, u_end; out_edge_iterator ai, a_end; // check edge flow values for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) { edge_descriptor a = *ai; if (capacity[a] > 0) if ((residual_capacity[a] + residual_capacity[reverse_edge[a]] != capacity[a] + capacity[reverse_edge[a]]) || (residual_capacity[a] < 0) || (residual_capacity[reverse_edge[a]] < 0)) return false; } } // check conservation FlowValue sum; for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) { vertex_descriptor u = *u_iter; if (u != src && u != sink) { if (excess_flow[u] != 0) return false; sum = 0; for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) if (capacity[*ai] > 0) sum -= capacity[*ai] - residual_capacity[*ai]; else sum += residual_capacity[*ai]; if (excess_flow[u] != sum) return false; } } return true; } // is_flow() bool is_optimal() { // check if mincut is saturated... global_distance_update(); return distance[src] >= n; } void print_statistics(std::ostream& os) const { os << "pushes: " << push_count << std::endl << "relabels: " << relabel_count << std::endl << "updates: " << update_count << std::endl << "gaps: " << gap_count << std::endl << "gap nodes: " << gap_node_count << std::endl << std::endl; } void print_flow_values(std::ostream& os) const { os << "flow values" << std::endl; vertex_iterator u_iter, u_end; out_edge_iterator ei, e_end; for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) if (capacity[*ei] > 0) os << *u_iter << " " << target(*ei, g) << " " << (capacity[*ei] - residual_capacity[*ei]) << std::endl; os << std::endl; } //======================================================================= Graph& g; vertices_size_type n; vertices_size_type nm; EdgeCapacityMap capacity; vertex_descriptor src; vertex_descriptor sink; VertexIndexMap index; // will need to use random_access_property_map with these std::vector< FlowValue > excess_flow; std::vector< out_edge_iterator > current; std::vector< distance_size_type > distance; std::vector< default_color_type > color; // Edge Property Maps that must be interior to the graph ReverseEdgeMap reverse_edge; ResidualCapacityEdgeMap residual_capacity; LayerArray layers; std::vector< list_iterator > layer_list_ptr; distance_size_type max_distance; // maximal distance distance_size_type max_active; // maximal distance with active node distance_size_type min_active; // minimal distance with active node boost::queue<vertex_descriptor> Q; // Statistics counters long push_count; long update_count; long relabel_count; long gap_count; long gap_node_count; inline double global_update_frequency() { return 0.5; } inline vertices_size_type alpha() { return 6; } inline long beta() { return 12; } long work_since_last_update; }; } // namespace detail template <class Graph, class CapacityEdgeMap, class ResidualCapacityEdgeMap, class ReverseEdgeMap, class VertexIndexMap> typename property_traits<CapacityEdgeMap>::value_type 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) { typedef typename property_traits<CapacityEdgeMap>::value_type FlowValue; detail::push_relabel<Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, ReverseEdgeMap, VertexIndexMap, FlowValue> algo(g, cap, res, rev, src, sink, index_map); FlowValue flow = algo.maximum_preflow(); algo.convert_preflow_to_flow(); assert(algo.is_flow()); assert(algo.is_optimal()); return flow; } // push_relabel_max_flow() template <class Graph, class P, class T, class R> 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) { return push_relabel_max_flow (g, src, sink, choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity), choose_pmap(get_param(params, edge_residual_capacity), g, edge_residual_capacity), choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse), choose_const_pmap(get_param(params, vertex_index), g, vertex_index) ); } template <class Graph> typename property_traits< 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) { bgl_named_params<int, buffer_param_t> params(0); // bogus empty param return push_relabel_max_flow(g, src, sink, params); } } // namespace boost #endif // BOOST_PUSH_RELABEL_MAX_FLOW_HPP