FCL/sf/os/ossrv: comparison ossrv_pub/boost_apis/boost/graph/push_relabel_max

equal deleted inserted replaced

--1:000000000000
+:e4d67989cc36
+//=======================================================================
+// 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