1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/epoc32/include/stdapis/boost/graph/push_relabel_max_flow.hpp Tue Mar 16 16:12:26 2010 +0000
1.3 @@ -0,0 +1,727 @@
1.4 +//=======================================================================
1.5 +// Copyright 2000 University of Notre Dame.
1.6 +// Authors: Jeremy G. Siek, Andrew Lumsdaine, Lie-Quan Lee
1.7 +//
1.8 +// Distributed under the Boost Software License, Version 1.0. (See
1.9 +// accompanying file LICENSE_1_0.txt or copy at
1.10 +// http://www.boost.org/LICENSE_1_0.txt)
1.11 +//=======================================================================
1.12 +
1.13 +#ifndef BOOST_PUSH_RELABEL_MAX_FLOW_HPP
1.14 +#define BOOST_PUSH_RELABEL_MAX_FLOW_HPP
1.15 +
1.16 +#include <boost/config.hpp>
1.17 +#include <cassert>
1.18 +#include <vector>
1.19 +#include <list>
1.20 +#include <iosfwd>
1.21 +#include <algorithm> // for std::min and std::max
1.22 +
1.23 +#include <boost/pending/queue.hpp>
1.24 +#include <boost/limits.hpp>
1.25 +#include <boost/graph/graph_concepts.hpp>
1.26 +#include <boost/graph/named_function_params.hpp>
1.27 +
1.28 +namespace boost {
1.29 +
1.30 + namespace detail {
1.31 +
1.32 + // This implementation is based on Goldberg's
1.33 + // "On Implementing Push-Relabel Method for the Maximum Flow Problem"
1.34 + // by B.V. Cherkassky and A.V. Goldberg, IPCO '95, pp. 157--171
1.35 + // and on the h_prf.c and hi_pr.c code written by the above authors.
1.36 +
1.37 + // This implements the highest-label version of the push-relabel method
1.38 + // with the global relabeling and gap relabeling heuristics.
1.39 +
1.40 + // The terms "rank", "distance", "height" are synonyms in
1.41 + // Goldberg's implementation, paper and in the CLR. A "layer" is a
1.42 + // group of vertices with the same distance. The vertices in each
1.43 + // layer are categorized as active or inactive. An active vertex
1.44 + // has positive excess flow and its distance is less than n (it is
1.45 + // not blocked).
1.46 +
1.47 + template <class Vertex>
1.48 + struct preflow_layer {
1.49 + std::list<Vertex> active_vertices;
1.50 + std::list<Vertex> inactive_vertices;
1.51 + };
1.52 +
1.53 + template <class Graph,
1.54 + class EdgeCapacityMap, // integer value type
1.55 + class ResidualCapacityEdgeMap,
1.56 + class ReverseEdgeMap,
1.57 + class VertexIndexMap, // vertex_descriptor -> integer
1.58 + class FlowValue>
1.59 + class push_relabel
1.60 + {
1.61 + public:
1.62 + typedef graph_traits<Graph> Traits;
1.63 + typedef typename Traits::vertex_descriptor vertex_descriptor;
1.64 + typedef typename Traits::edge_descriptor edge_descriptor;
1.65 + typedef typename Traits::vertex_iterator vertex_iterator;
1.66 + typedef typename Traits::out_edge_iterator out_edge_iterator;
1.67 + typedef typename Traits::vertices_size_type vertices_size_type;
1.68 + typedef typename Traits::edges_size_type edges_size_type;
1.69 +
1.70 + typedef preflow_layer<vertex_descriptor> Layer;
1.71 + typedef std::vector< Layer > LayerArray;
1.72 + typedef typename LayerArray::iterator layer_iterator;
1.73 + typedef typename LayerArray::size_type distance_size_type;
1.74 +
1.75 + typedef color_traits<default_color_type> ColorTraits;
1.76 +
1.77 + //=======================================================================
1.78 + // Some helper predicates
1.79 +
1.80 + inline bool is_admissible(vertex_descriptor u, vertex_descriptor v) {
1.81 + return distance[u] == distance[v] + 1;
1.82 + }
1.83 + inline bool is_residual_edge(edge_descriptor a) {
1.84 + return 0 < residual_capacity[a];
1.85 + }
1.86 + inline bool is_saturated(edge_descriptor a) {
1.87 + return residual_capacity[a] == 0;
1.88 + }
1.89 +
1.90 + //=======================================================================
1.91 + // Layer List Management Functions
1.92 +
1.93 + typedef typename std::list<vertex_descriptor>::iterator list_iterator;
1.94 +
1.95 + void add_to_active_list(vertex_descriptor u, Layer& layer) {
1.96 + BOOST_USING_STD_MIN();
1.97 + BOOST_USING_STD_MAX();
1.98 + layer.active_vertices.push_front(u);
1.99 + max_active = max BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], max_active);
1.100 + min_active = min BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], min_active);
1.101 + layer_list_ptr[u] = layer.active_vertices.begin();
1.102 + }
1.103 + void remove_from_active_list(vertex_descriptor u) {
1.104 + layers[distance[u]].active_vertices.erase(layer_list_ptr[u]);
1.105 + }
1.106 +
1.107 + void add_to_inactive_list(vertex_descriptor u, Layer& layer) {
1.108 + layer.inactive_vertices.push_front(u);
1.109 + layer_list_ptr[u] = layer.inactive_vertices.begin();
1.110 + }
1.111 + void remove_from_inactive_list(vertex_descriptor u) {
1.112 + layers[distance[u]].inactive_vertices.erase(layer_list_ptr[u]);
1.113 + }
1.114 +
1.115 + //=======================================================================
1.116 + // initialization
1.117 + push_relabel(Graph& g_,
1.118 + EdgeCapacityMap cap,
1.119 + ResidualCapacityEdgeMap res,
1.120 + ReverseEdgeMap rev,
1.121 + vertex_descriptor src_,
1.122 + vertex_descriptor sink_,
1.123 + VertexIndexMap idx)
1.124 + : g(g_), n(num_vertices(g_)), capacity(cap), src(src_), sink(sink_),
1.125 + index(idx),
1.126 + excess_flow(num_vertices(g_)),
1.127 + current(num_vertices(g_), out_edges(*vertices(g_).first, g_).second),
1.128 + distance(num_vertices(g_)),
1.129 + color(num_vertices(g_)),
1.130 + reverse_edge(rev),
1.131 + residual_capacity(res),
1.132 + layers(num_vertices(g_)),
1.133 + layer_list_ptr(num_vertices(g_),
1.134 + layers.front().inactive_vertices.end()),
1.135 + push_count(0), update_count(0), relabel_count(0),
1.136 + gap_count(0), gap_node_count(0),
1.137 + work_since_last_update(0)
1.138 + {
1.139 + vertex_iterator u_iter, u_end;
1.140 + // Don't count the reverse edges
1.141 + edges_size_type m = num_edges(g) / 2;
1.142 + nm = alpha() * n + m;
1.143 +
1.144 + // Initialize flow to zero which means initializing
1.145 + // the residual capacity to equal the capacity.
1.146 + out_edge_iterator ei, e_end;
1.147 + for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
1.148 + for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) {
1.149 + residual_capacity[*ei] = capacity[*ei];
1.150 + }
1.151 +
1.152 + for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
1.153 + vertex_descriptor u = *u_iter;
1.154 + excess_flow[u] = 0;
1.155 + current[u] = out_edges(u, g).first;
1.156 + }
1.157 +
1.158 + bool overflow_detected = false;
1.159 + FlowValue test_excess = 0;
1.160 +
1.161 + out_edge_iterator a_iter, a_end;
1.162 + for (tie(a_iter, a_end) = out_edges(src, g); a_iter != a_end; ++a_iter)
1.163 + if (target(*a_iter, g) != src)
1.164 + test_excess += residual_capacity[*a_iter];
1.165 + if (test_excess > (std::numeric_limits<FlowValue>::max)())
1.166 + overflow_detected = true;
1.167 +
1.168 + if (overflow_detected)
1.169 + excess_flow[src] = (std::numeric_limits<FlowValue>::max)();
1.170 + else {
1.171 + excess_flow[src] = 0;
1.172 + for (tie(a_iter, a_end) = out_edges(src, g);
1.173 + a_iter != a_end; ++a_iter) {
1.174 + edge_descriptor a = *a_iter;
1.175 + if (target(a, g) != src) {
1.176 + ++push_count;
1.177 + FlowValue delta = residual_capacity[a];
1.178 + residual_capacity[a] -= delta;
1.179 + residual_capacity[reverse_edge[a]] += delta;
1.180 + excess_flow[target(a, g)] += delta;
1.181 + }
1.182 + }
1.183 + }
1.184 + max_distance = num_vertices(g) - 1;
1.185 + max_active = 0;
1.186 + min_active = n;
1.187 +
1.188 + for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
1.189 + vertex_descriptor u = *u_iter;
1.190 + if (u == sink) {
1.191 + distance[u] = 0;
1.192 + continue;
1.193 + } else if (u == src && !overflow_detected)
1.194 + distance[u] = n;
1.195 + else
1.196 + distance[u] = 1;
1.197 +
1.198 + if (excess_flow[u] > 0)
1.199 + add_to_active_list(u, layers[1]);
1.200 + else if (distance[u] < n)
1.201 + add_to_inactive_list(u, layers[1]);
1.202 + }
1.203 +
1.204 + } // push_relabel constructor
1.205 +
1.206 + //=======================================================================
1.207 + // This is a breadth-first search over the residual graph
1.208 + // (well, actually the reverse of the residual graph).
1.209 + // Would be cool to have a graph view adaptor for hiding certain
1.210 + // edges, like the saturated (non-residual) edges in this case.
1.211 + // Goldberg's implementation abused "distance" for the coloring.
1.212 + void global_distance_update()
1.213 + {
1.214 + BOOST_USING_STD_MAX();
1.215 + ++update_count;
1.216 + vertex_iterator u_iter, u_end;
1.217 + for (tie(u_iter,u_end) = vertices(g); u_iter != u_end; ++u_iter) {
1.218 + color[*u_iter] = ColorTraits::white();
1.219 + distance[*u_iter] = n;
1.220 + }
1.221 + color[sink] = ColorTraits::gray();
1.222 + distance[sink] = 0;
1.223 +
1.224 + for (distance_size_type l = 0; l <= max_distance; ++l) {
1.225 + layers[l].active_vertices.clear();
1.226 + layers[l].inactive_vertices.clear();
1.227 + }
1.228 +
1.229 + max_distance = max_active = 0;
1.230 + min_active = n;
1.231 +
1.232 + Q.push(sink);
1.233 + while (! Q.empty()) {
1.234 + vertex_descriptor u = Q.top();
1.235 + Q.pop();
1.236 + distance_size_type d_v = distance[u] + 1;
1.237 +
1.238 + out_edge_iterator ai, a_end;
1.239 + for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) {
1.240 + edge_descriptor a = *ai;
1.241 + vertex_descriptor v = target(a, g);
1.242 + if (color[v] == ColorTraits::white()
1.243 + && is_residual_edge(reverse_edge[a])) {
1.244 + distance[v] = d_v;
1.245 + color[v] = ColorTraits::gray();
1.246 + current[v] = out_edges(v, g).first;
1.247 + max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(d_v, max_distance);
1.248 +
1.249 + if (excess_flow[v] > 0)
1.250 + add_to_active_list(v, layers[d_v]);
1.251 + else
1.252 + add_to_inactive_list(v, layers[d_v]);
1.253 +
1.254 + Q.push(v);
1.255 + }
1.256 + }
1.257 + }
1.258 + } // global_distance_update()
1.259 +
1.260 + //=======================================================================
1.261 + // This function is called "push" in Goldberg's h_prf implementation,
1.262 + // but it is called "discharge" in the paper and in hi_pr.c.
1.263 + void discharge(vertex_descriptor u)
1.264 + {
1.265 + assert(excess_flow[u] > 0);
1.266 + while (1) {
1.267 + out_edge_iterator ai, ai_end;
1.268 + for (ai = current[u], ai_end = out_edges(u, g).second;
1.269 + ai != ai_end; ++ai) {
1.270 + edge_descriptor a = *ai;
1.271 + if (is_residual_edge(a)) {
1.272 + vertex_descriptor v = target(a, g);
1.273 + if (is_admissible(u, v)) {
1.274 + ++push_count;
1.275 + if (v != sink && excess_flow[v] == 0) {
1.276 + remove_from_inactive_list(v);
1.277 + add_to_active_list(v, layers[distance[v]]);
1.278 + }
1.279 + push_flow(a);
1.280 + if (excess_flow[u] == 0)
1.281 + break;
1.282 + }
1.283 + }
1.284 + } // for out_edges of i starting from current
1.285 +
1.286 + Layer& layer = layers[distance[u]];
1.287 + distance_size_type du = distance[u];
1.288 +
1.289 + if (ai == ai_end) { // i must be relabeled
1.290 + relabel_distance(u);
1.291 + if (layer.active_vertices.empty()
1.292 + && layer.inactive_vertices.empty())
1.293 + gap(du);
1.294 + if (distance[u] == n)
1.295 + break;
1.296 + } else { // i is no longer active
1.297 + current[u] = ai;
1.298 + add_to_inactive_list(u, layer);
1.299 + break;
1.300 + }
1.301 + } // while (1)
1.302 + } // discharge()
1.303 +
1.304 + //=======================================================================
1.305 + // This corresponds to the "push" update operation of the paper,
1.306 + // not the "push" function in Goldberg's h_prf.c implementation.
1.307 + // The idea is to push the excess flow from from vertex u to v.
1.308 + void push_flow(edge_descriptor u_v)
1.309 + {
1.310 + vertex_descriptor
1.311 + u = source(u_v, g),
1.312 + v = target(u_v, g);
1.313 +
1.314 + BOOST_USING_STD_MIN();
1.315 + FlowValue flow_delta
1.316 + = min BOOST_PREVENT_MACRO_SUBSTITUTION(excess_flow[u], residual_capacity[u_v]);
1.317 +
1.318 + residual_capacity[u_v] -= flow_delta;
1.319 + residual_capacity[reverse_edge[u_v]] += flow_delta;
1.320 +
1.321 + excess_flow[u] -= flow_delta;
1.322 + excess_flow[v] += flow_delta;
1.323 + } // push_flow()
1.324 +
1.325 + //=======================================================================
1.326 + // The main purpose of this routine is to set distance[v]
1.327 + // to the smallest value allowed by the valid labeling constraints,
1.328 + // which are:
1.329 + // distance[t] = 0
1.330 + // distance[u] <= distance[v] + 1 for every residual edge (u,v)
1.331 + //
1.332 + distance_size_type relabel_distance(vertex_descriptor u)
1.333 + {
1.334 + BOOST_USING_STD_MAX();
1.335 + ++relabel_count;
1.336 + work_since_last_update += beta();
1.337 +
1.338 + distance_size_type min_distance = num_vertices(g);
1.339 + distance[u] = min_distance;
1.340 +
1.341 + // Examine the residual out-edges of vertex i, choosing the
1.342 + // edge whose target vertex has the minimal distance.
1.343 + out_edge_iterator ai, a_end, min_edge_iter;
1.344 + for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) {
1.345 + ++work_since_last_update;
1.346 + edge_descriptor a = *ai;
1.347 + vertex_descriptor v = target(a, g);
1.348 + if (is_residual_edge(a) && distance[v] < min_distance) {
1.349 + min_distance = distance[v];
1.350 + min_edge_iter = ai;
1.351 + }
1.352 + }
1.353 + ++min_distance;
1.354 + if (min_distance < n) {
1.355 + distance[u] = min_distance; // this is the main action
1.356 + current[u] = min_edge_iter;
1.357 + max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(min_distance, max_distance);
1.358 + }
1.359 + return min_distance;
1.360 + } // relabel_distance()
1.361 +
1.362 + //=======================================================================
1.363 + // cleanup beyond the gap
1.364 + void gap(distance_size_type empty_distance)
1.365 + {
1.366 + ++gap_count;
1.367 +
1.368 + distance_size_type r; // distance of layer before the current layer
1.369 + r = empty_distance - 1;
1.370 +
1.371 + // Set the distance for the vertices beyond the gap to "infinity".
1.372 + for (layer_iterator l = layers.begin() + empty_distance + 1;
1.373 + l < layers.begin() + max_distance; ++l) {
1.374 + list_iterator i;
1.375 + for (i = l->inactive_vertices.begin();
1.376 + i != l->inactive_vertices.end(); ++i) {
1.377 + distance[*i] = n;
1.378 + ++gap_node_count;
1.379 + }
1.380 + l->inactive_vertices.clear();
1.381 + }
1.382 + max_distance = r;
1.383 + max_active = r;
1.384 + }
1.385 +
1.386 + //=======================================================================
1.387 + // This is the core part of the algorithm, "phase one".
1.388 + FlowValue maximum_preflow()
1.389 + {
1.390 + work_since_last_update = 0;
1.391 +
1.392 + while (max_active >= min_active) { // "main" loop
1.393 +
1.394 + Layer& layer = layers[max_active];
1.395 + list_iterator u_iter = layer.active_vertices.begin();
1.396 +
1.397 + if (u_iter == layer.active_vertices.end())
1.398 + --max_active;
1.399 + else {
1.400 + vertex_descriptor u = *u_iter;
1.401 + remove_from_active_list(u);
1.402 +
1.403 + discharge(u);
1.404 +
1.405 + if (work_since_last_update * global_update_frequency() > nm) {
1.406 + global_distance_update();
1.407 + work_since_last_update = 0;
1.408 + }
1.409 + }
1.410 + } // while (max_active >= min_active)
1.411 +
1.412 + return excess_flow[sink];
1.413 + } // maximum_preflow()
1.414 +
1.415 + //=======================================================================
1.416 + // remove excess flow, the "second phase"
1.417 + // This does a DFS on the reverse flow graph of nodes with excess flow.
1.418 + // If a cycle is found, cancel it.
1.419 + // Return the nodes with excess flow in topological order.
1.420 + //
1.421 + // Unlike the prefl_to_flow() implementation, we use
1.422 + // "color" instead of "distance" for the DFS labels
1.423 + // "parent" instead of nl_prev for the DFS tree
1.424 + // "topo_next" instead of nl_next for the topological ordering
1.425 + void convert_preflow_to_flow()
1.426 + {
1.427 + vertex_iterator u_iter, u_end;
1.428 + out_edge_iterator ai, a_end;
1.429 +
1.430 + vertex_descriptor r, restart, u;
1.431 +
1.432 + std::vector<vertex_descriptor> parent(n);
1.433 + std::vector<vertex_descriptor> topo_next(n);
1.434 +
1.435 + vertex_descriptor tos(parent[0]),
1.436 + bos(parent[0]); // bogus initialization, just to avoid warning
1.437 + bool bos_null = true;
1.438 +
1.439 + // handle self-loops
1.440 + for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
1.441 + for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai)
1.442 + if (target(*ai, g) == *u_iter)
1.443 + residual_capacity[*ai] = capacity[*ai];
1.444 +
1.445 + // initialize
1.446 + for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
1.447 + u = *u_iter;
1.448 + color[u] = ColorTraits::white();
1.449 + parent[u] = u;
1.450 + current[u] = out_edges(u, g).first;
1.451 + }
1.452 + // eliminate flow cycles and topologically order the vertices
1.453 + for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
1.454 + u = *u_iter;
1.455 + if (color[u] == ColorTraits::white()
1.456 + && excess_flow[u] > 0
1.457 + && u != src && u != sink ) {
1.458 + r = u;
1.459 + color[r] = ColorTraits::gray();
1.460 + while (1) {
1.461 + for (; current[u] != out_edges(u, g).second; ++current[u]) {
1.462 + edge_descriptor a = *current[u];
1.463 + if (capacity[a] == 0 && is_residual_edge(a)) {
1.464 + vertex_descriptor v = target(a, g);
1.465 + if (color[v] == ColorTraits::white()) {
1.466 + color[v] = ColorTraits::gray();
1.467 + parent[v] = u;
1.468 + u = v;
1.469 + break;
1.470 + } else if (color[v] == ColorTraits::gray()) {
1.471 + // find minimum flow on the cycle
1.472 + FlowValue delta = residual_capacity[a];
1.473 + while (1) {
1.474 + BOOST_USING_STD_MIN();
1.475 + delta = min BOOST_PREVENT_MACRO_SUBSTITUTION(delta, residual_capacity[*current[v]]);
1.476 + if (v == u)
1.477 + break;
1.478 + else
1.479 + v = target(*current[v], g);
1.480 + }
1.481 + // remove delta flow units
1.482 + v = u;
1.483 + while (1) {
1.484 + a = *current[v];
1.485 + residual_capacity[a] -= delta;
1.486 + residual_capacity[reverse_edge[a]] += delta;
1.487 + v = target(a, g);
1.488 + if (v == u)
1.489 + break;
1.490 + }
1.491 +
1.492 + // back-out of DFS to the first saturated edge
1.493 + restart = u;
1.494 + for (v = target(*current[u], g); v != u; v = target(a, g)){
1.495 + a = *current[v];
1.496 + if (color[v] == ColorTraits::white()
1.497 + || is_saturated(a)) {
1.498 + color[target(*current[v], g)] = ColorTraits::white();
1.499 + if (color[v] != ColorTraits::white())
1.500 + restart = v;
1.501 + }
1.502 + }
1.503 + if (restart != u) {
1.504 + u = restart;
1.505 + ++current[u];
1.506 + break;
1.507 + }
1.508 + } // else if (color[v] == ColorTraits::gray())
1.509 + } // if (capacity[a] == 0 ...
1.510 + } // for out_edges(u, g) (though "u" changes during loop)
1.511 +
1.512 + if (current[u] == out_edges(u, g).second) {
1.513 + // scan of i is complete
1.514 + color[u] = ColorTraits::black();
1.515 + if (u != src) {
1.516 + if (bos_null) {
1.517 + bos = u;
1.518 + bos_null = false;
1.519 + tos = u;
1.520 + } else {
1.521 + topo_next[u] = tos;
1.522 + tos = u;
1.523 + }
1.524 + }
1.525 + if (u != r) {
1.526 + u = parent[u];
1.527 + ++current[u];
1.528 + } else
1.529 + break;
1.530 + }
1.531 + } // while (1)
1.532 + } // if (color[u] == white && excess_flow[u] > 0 & ...)
1.533 + } // for all vertices in g
1.534 +
1.535 + // return excess flows
1.536 + // note that the sink is not on the stack
1.537 + if (! bos_null) {
1.538 + for (u = tos; u != bos; u = topo_next[u]) {
1.539 + ai = out_edges(u, g).first;
1.540 + while (excess_flow[u] > 0 && ai != out_edges(u, g).second) {
1.541 + if (capacity[*ai] == 0 && is_residual_edge(*ai))
1.542 + push_flow(*ai);
1.543 + ++ai;
1.544 + }
1.545 + }
1.546 + // do the bottom
1.547 + u = bos;
1.548 + ai = out_edges(u, g).first;
1.549 + while (excess_flow[u] > 0) {
1.550 + if (capacity[*ai] == 0 && is_residual_edge(*ai))
1.551 + push_flow(*ai);
1.552 + ++ai;
1.553 + }
1.554 + }
1.555 +
1.556 + } // convert_preflow_to_flow()
1.557 +
1.558 + //=======================================================================
1.559 + inline bool is_flow()
1.560 + {
1.561 + vertex_iterator u_iter, u_end;
1.562 + out_edge_iterator ai, a_end;
1.563 +
1.564 + // check edge flow values
1.565 + for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
1.566 + for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) {
1.567 + edge_descriptor a = *ai;
1.568 + if (capacity[a] > 0)
1.569 + if ((residual_capacity[a] + residual_capacity[reverse_edge[a]]
1.570 + != capacity[a] + capacity[reverse_edge[a]])
1.571 + || (residual_capacity[a] < 0)
1.572 + || (residual_capacity[reverse_edge[a]] < 0))
1.573 + return false;
1.574 + }
1.575 + }
1.576 +
1.577 + // check conservation
1.578 + FlowValue sum;
1.579 + for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
1.580 + vertex_descriptor u = *u_iter;
1.581 + if (u != src && u != sink) {
1.582 + if (excess_flow[u] != 0)
1.583 + return false;
1.584 + sum = 0;
1.585 + for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai)
1.586 + if (capacity[*ai] > 0)
1.587 + sum -= capacity[*ai] - residual_capacity[*ai];
1.588 + else
1.589 + sum += residual_capacity[*ai];
1.590 +
1.591 + if (excess_flow[u] != sum)
1.592 + return false;
1.593 + }
1.594 + }
1.595 +
1.596 + return true;
1.597 + } // is_flow()
1.598 +
1.599 + bool is_optimal() {
1.600 + // check if mincut is saturated...
1.601 + global_distance_update();
1.602 + return distance[src] >= n;
1.603 + }
1.604 +
1.605 + void print_statistics(std::ostream& os) const {
1.606 + os << "pushes: " << push_count << std::endl
1.607 + << "relabels: " << relabel_count << std::endl
1.608 + << "updates: " << update_count << std::endl
1.609 + << "gaps: " << gap_count << std::endl
1.610 + << "gap nodes: " << gap_node_count << std::endl
1.611 + << std::endl;
1.612 + }
1.613 +
1.614 + void print_flow_values(std::ostream& os) const {
1.615 + os << "flow values" << std::endl;
1.616 + vertex_iterator u_iter, u_end;
1.617 + out_edge_iterator ei, e_end;
1.618 + for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
1.619 + for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei)
1.620 + if (capacity[*ei] > 0)
1.621 + os << *u_iter << " " << target(*ei, g) << " "
1.622 + << (capacity[*ei] - residual_capacity[*ei]) << std::endl;
1.623 + os << std::endl;
1.624 + }
1.625 +
1.626 + //=======================================================================
1.627 +
1.628 + Graph& g;
1.629 + vertices_size_type n;
1.630 + vertices_size_type nm;
1.631 + EdgeCapacityMap capacity;
1.632 + vertex_descriptor src;
1.633 + vertex_descriptor sink;
1.634 + VertexIndexMap index;
1.635 +
1.636 + // will need to use random_access_property_map with these
1.637 + std::vector< FlowValue > excess_flow;
1.638 + std::vector< out_edge_iterator > current;
1.639 + std::vector< distance_size_type > distance;
1.640 + std::vector< default_color_type > color;
1.641 +
1.642 + // Edge Property Maps that must be interior to the graph
1.643 + ReverseEdgeMap reverse_edge;
1.644 + ResidualCapacityEdgeMap residual_capacity;
1.645 +
1.646 + LayerArray layers;
1.647 + std::vector< list_iterator > layer_list_ptr;
1.648 + distance_size_type max_distance; // maximal distance
1.649 + distance_size_type max_active; // maximal distance with active node
1.650 + distance_size_type min_active; // minimal distance with active node
1.651 + boost::queue<vertex_descriptor> Q;
1.652 +
1.653 + // Statistics counters
1.654 + long push_count;
1.655 + long update_count;
1.656 + long relabel_count;
1.657 + long gap_count;
1.658 + long gap_node_count;
1.659 +
1.660 + inline double global_update_frequency() { return 0.5; }
1.661 + inline vertices_size_type alpha() { return 6; }
1.662 + inline long beta() { return 12; }
1.663 +
1.664 + long work_since_last_update;
1.665 + };
1.666 +
1.667 + } // namespace detail
1.668 +
1.669 + template <class Graph,
1.670 + class CapacityEdgeMap, class ResidualCapacityEdgeMap,
1.671 + class ReverseEdgeMap, class VertexIndexMap>
1.672 + typename property_traits<CapacityEdgeMap>::value_type
1.673 + push_relabel_max_flow
1.674 + (Graph& g,
1.675 + typename graph_traits<Graph>::vertex_descriptor src,
1.676 + typename graph_traits<Graph>::vertex_descriptor sink,
1.677 + CapacityEdgeMap cap, ResidualCapacityEdgeMap res,
1.678 + ReverseEdgeMap rev, VertexIndexMap index_map)
1.679 + {
1.680 + typedef typename property_traits<CapacityEdgeMap>::value_type FlowValue;
1.681 +
1.682 + detail::push_relabel<Graph, CapacityEdgeMap, ResidualCapacityEdgeMap,
1.683 + ReverseEdgeMap, VertexIndexMap, FlowValue>
1.684 + algo(g, cap, res, rev, src, sink, index_map);
1.685 +
1.686 + FlowValue flow = algo.maximum_preflow();
1.687 +
1.688 + algo.convert_preflow_to_flow();
1.689 +
1.690 + assert(algo.is_flow());
1.691 + assert(algo.is_optimal());
1.692 +
1.693 + return flow;
1.694 + } // push_relabel_max_flow()
1.695 +
1.696 + template <class Graph, class P, class T, class R>
1.697 + typename detail::edge_capacity_value<Graph, P, T, R>::type
1.698 + push_relabel_max_flow
1.699 + (Graph& g,
1.700 + typename graph_traits<Graph>::vertex_descriptor src,
1.701 + typename graph_traits<Graph>::vertex_descriptor sink,
1.702 + const bgl_named_params<P, T, R>& params)
1.703 + {
1.704 + return push_relabel_max_flow
1.705 + (g, src, sink,
1.706 + choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity),
1.707 + choose_pmap(get_param(params, edge_residual_capacity),
1.708 + g, edge_residual_capacity),
1.709 + choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse),
1.710 + choose_const_pmap(get_param(params, vertex_index), g, vertex_index)
1.711 + );
1.712 + }
1.713 +
1.714 + template <class Graph>
1.715 + typename property_traits<
1.716 + typename property_map<Graph, edge_capacity_t>::const_type
1.717 + >::value_type
1.718 + push_relabel_max_flow
1.719 + (Graph& g,
1.720 + typename graph_traits<Graph>::vertex_descriptor src,
1.721 + typename graph_traits<Graph>::vertex_descriptor sink)
1.722 + {
1.723 + bgl_named_params<int, buffer_param_t> params(0); // bogus empty param
1.724 + return push_relabel_max_flow(g, src, sink, params);
1.725 + }
1.726 +
1.727 +} // namespace boost
1.728 +
1.729 +#endif // BOOST_PUSH_RELABEL_MAX_FLOW_HPP
1.730 +