Google Git
Sign in
android / platform / external / boost / ac861f8c0f33538060790a8e50701464ca9982d3 / . / boost / graph / push_relabel_max_flow.hpp
blob: 2b65470979d3fa0a8ba4e6752f415bc8df260604 [file] [log] [blame]
//=======================================================================
// 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_)),
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);
}
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);
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 (tie(ai, ai_end) = current[u]; 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].first = 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].first = 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);
}
// 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].first != current[u].second; ++current[u].first) {
edge_descriptor a = *current[u].first;
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].first]);
if (v == u)
break;
else
v = target(*current[v].first, g);
}
// remove delta flow units
v = u;
while (1) {
a = *current[v].first;
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].first, g); v != u; v = target(a, g)){
a = *current[v].first;
if (color[v] == ColorTraits::white()
|| is_saturated(a)) {
color[target(*current[v].first, g)] = ColorTraits::white();
if (color[v] != ColorTraits::white())
restart = v;
}
}
if (restart != u) {
u = restart;
++current[u].first;
break;
}
} // else if (color[v] == ColorTraits::gray())
} // if (capacity[a] == 0 ...
} // for out_edges(u, g) (though "u" changes during loop)
if ( current[u].first == current[u].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].first;
} 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]) {
tie(ai, a_end) = out_edges(u, g);
while (excess_flow[u] > 0 && ai != a_end) {
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< std::pair<out_edge_iterator, 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