[857] | 1 | //=======================================================================
|
---|
| 2 | // Copyright 2000 University of Notre Dame.
|
---|
| 3 | // Authors: Jeremy G. Siek, Andrew Lumsdaine, Lie-Quan Lee
|
---|
| 4 | //
|
---|
| 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 | //=======================================================================
|
---|
| 9 |
|
---|
| 10 | #ifndef BOOST_PUSH_RELABEL_MAX_FLOW_HPP
|
---|
| 11 | #define BOOST_PUSH_RELABEL_MAX_FLOW_HPP
|
---|
| 12 |
|
---|
| 13 | #include <boost/config.hpp>
|
---|
| 14 | #include <cassert>
|
---|
| 15 | #include <vector>
|
---|
| 16 | #include <list>
|
---|
| 17 | #include <iosfwd>
|
---|
| 18 | #include <algorithm> // for std::min and std::max
|
---|
| 19 |
|
---|
| 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>
|
---|
| 24 |
|
---|
| 25 | namespace boost {
|
---|
| 26 |
|
---|
| 27 | namespace detail {
|
---|
| 28 |
|
---|
| 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.
|
---|
| 33 |
|
---|
| 34 | // This implements the highest-label version of the push-relabel method
|
---|
| 35 | // with the global relabeling and gap relabeling heuristics.
|
---|
| 36 |
|
---|
| 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
|
---|
| 42 | // not blocked).
|
---|
| 43 |
|
---|
| 44 | template <class Vertex>
|
---|
| 45 | struct preflow_layer {
|
---|
| 46 | std::list<Vertex> active_vertices;
|
---|
| 47 | std::list<Vertex> inactive_vertices;
|
---|
| 48 | };
|
---|
| 49 |
|
---|
| 50 | template <class Graph,
|
---|
| 51 | class EdgeCapacityMap, // integer value type
|
---|
| 52 | class ResidualCapacityEdgeMap,
|
---|
| 53 | class ReverseEdgeMap,
|
---|
| 54 | class VertexIndexMap, // vertex_descriptor -> integer
|
---|
| 55 | class FlowValue>
|
---|
| 56 | class push_relabel
|
---|
| 57 | {
|
---|
| 58 | public:
|
---|
| 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;
|
---|
| 66 |
|
---|
| 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;
|
---|
| 71 |
|
---|
| 72 | typedef color_traits<default_color_type> ColorTraits;
|
---|
| 73 |
|
---|
| 74 | //=======================================================================
|
---|
| 75 | // Some helper predicates
|
---|
| 76 |
|
---|
| 77 | inline bool is_admissible(vertex_descriptor u, vertex_descriptor v) {
|
---|
| 78 | return distance[u] == distance[v] + 1;
|
---|
| 79 | }
|
---|
| 80 | inline bool is_residual_edge(edge_descriptor a) {
|
---|
| 81 | return 0 < residual_capacity[a];
|
---|
| 82 | }
|
---|
| 83 | inline bool is_saturated(edge_descriptor a) {
|
---|
| 84 | return residual_capacity[a] == 0;
|
---|
| 85 | }
|
---|
| 86 |
|
---|
| 87 | //=======================================================================
|
---|
| 88 | // Layer List Management Functions
|
---|
| 89 |
|
---|
| 90 | typedef typename std::list<vertex_descriptor>::iterator list_iterator;
|
---|
| 91 |
|
---|
| 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(distance[u], max_active);
|
---|
| 97 | min_active = min BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], min_active);
|
---|
| 98 | layer_list_ptr[u] = layer.active_vertices.begin();
|
---|
| 99 | }
|
---|
| 100 | void remove_from_active_list(vertex_descriptor u) {
|
---|
| 101 | layers[distance[u]].active_vertices.erase(layer_list_ptr[u]);
|
---|
| 102 | }
|
---|
| 103 |
|
---|
| 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();
|
---|
| 107 | }
|
---|
| 108 | void remove_from_inactive_list(vertex_descriptor u) {
|
---|
| 109 | layers[distance[u]].inactive_vertices.erase(layer_list_ptr[u]);
|
---|
| 110 | }
|
---|
| 111 |
|
---|
| 112 | //=======================================================================
|
---|
| 113 | // initialization
|
---|
| 114 | push_relabel(Graph& g_,
|
---|
| 115 | EdgeCapacityMap cap,
|
---|
| 116 | ResidualCapacityEdgeMap res,
|
---|
| 117 | ReverseEdgeMap rev,
|
---|
| 118 | vertex_descriptor src_,
|
---|
| 119 | vertex_descriptor sink_,
|
---|
| 120 | VertexIndexMap idx)
|
---|
| 121 | : g(g_), n(num_vertices(g_)), capacity(cap), src(src_), sink(sink_),
|
---|
| 122 | index(idx),
|
---|
| 123 | excess_flow(num_vertices(g_)),
|
---|
| 124 | current(num_vertices(g_), out_edges(*vertices(g_).first, g_).second),
|
---|
| 125 | distance(num_vertices(g_)),
|
---|
| 126 | color(num_vertices(g_)),
|
---|
| 127 | reverse_edge(rev),
|
---|
| 128 | residual_capacity(res),
|
---|
| 129 | layers(num_vertices(g_)),
|
---|
| 130 | layer_list_ptr(num_vertices(g_),
|
---|
| 131 | layers.front().inactive_vertices.end()),
|
---|
| 132 | push_count(0), update_count(0), relabel_count(0),
|
---|
| 133 | gap_count(0), gap_node_count(0),
|
---|
| 134 | work_since_last_update(0)
|
---|
| 135 | {
|
---|
| 136 | vertex_iterator u_iter, u_end;
|
---|
| 137 | // Don't count the reverse edges
|
---|
| 138 | edges_size_type m = num_edges(g) / 2;
|
---|
| 139 | nm = alpha() * n + m;
|
---|
| 140 |
|
---|
| 141 | // Initialize flow to zero which means initializing
|
---|
| 142 | // the residual capacity to equal the capacity.
|
---|
| 143 | out_edge_iterator ei, e_end;
|
---|
| 144 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
|
---|
| 145 | for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) {
|
---|
| 146 | residual_capacity[*ei] = capacity[*ei];
|
---|
| 147 | }
|
---|
| 148 |
|
---|
| 149 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
|
---|
| 150 | vertex_descriptor u = *u_iter;
|
---|
| 151 | excess_flow[u] = 0;
|
---|
| 152 | current[u] = out_edges(u, g).first;
|
---|
| 153 | }
|
---|
| 154 |
|
---|
| 155 | bool overflow_detected = false;
|
---|
| 156 | FlowValue test_excess = 0;
|
---|
| 157 |
|
---|
| 158 | out_edge_iterator a_iter, a_end;
|
---|
| 159 | for (tie(a_iter, a_end) = out_edges(src, g); a_iter != a_end; ++a_iter)
|
---|
| 160 | if (target(*a_iter, g) != src)
|
---|
| 161 | test_excess += residual_capacity[*a_iter];
|
---|
| 162 | if (test_excess > (std::numeric_limits<FlowValue>::max)())
|
---|
| 163 | overflow_detected = true;
|
---|
| 164 |
|
---|
| 165 | if (overflow_detected)
|
---|
| 166 | excess_flow[src] = (std::numeric_limits<FlowValue>::max)();
|
---|
| 167 | else {
|
---|
| 168 | excess_flow[src] = 0;
|
---|
| 169 | for (tie(a_iter, a_end) = out_edges(src, g);
|
---|
| 170 | a_iter != a_end; ++a_iter) {
|
---|
| 171 | edge_descriptor a = *a_iter;
|
---|
| 172 | if (target(a, g) != src) {
|
---|
| 173 | ++push_count;
|
---|
| 174 | FlowValue delta = residual_capacity[a];
|
---|
| 175 | residual_capacity[a] -= delta;
|
---|
| 176 | residual_capacity[reverse_edge[a]] += delta;
|
---|
| 177 | excess_flow[target(a, g)] += delta;
|
---|
| 178 | }
|
---|
| 179 | }
|
---|
| 180 | }
|
---|
| 181 | max_distance = num_vertices(g) - 1;
|
---|
| 182 | max_active = 0;
|
---|
| 183 | min_active = n;
|
---|
| 184 |
|
---|
| 185 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
|
---|
| 186 | vertex_descriptor u = *u_iter;
|
---|
| 187 | if (u == sink) {
|
---|
| 188 | distance[u] = 0;
|
---|
| 189 | continue;
|
---|
| 190 | } else if (u == src && !overflow_detected)
|
---|
| 191 | distance[u] = n;
|
---|
| 192 | else
|
---|
| 193 | distance[u] = 1;
|
---|
| 194 |
|
---|
| 195 | if (excess_flow[u] > 0)
|
---|
| 196 | add_to_active_list(u, layers[1]);
|
---|
| 197 | else if (distance[u] < n)
|
---|
| 198 | add_to_inactive_list(u, layers[1]);
|
---|
| 199 | }
|
---|
| 200 |
|
---|
| 201 | } // push_relabel constructor
|
---|
| 202 |
|
---|
| 203 | //=======================================================================
|
---|
| 204 | // This is a breadth-first search over the residual graph
|
---|
| 205 | // (well, actually the reverse of the residual graph).
|
---|
| 206 | // Would be cool to have a graph view adaptor for hiding certain
|
---|
| 207 | // edges, like the saturated (non-residual) edges in this case.
|
---|
| 208 | // Goldberg's implementation abused "distance" for the coloring.
|
---|
| 209 | void global_distance_update()
|
---|
| 210 | {
|
---|
| 211 | BOOST_USING_STD_MAX();
|
---|
| 212 | ++update_count;
|
---|
| 213 | vertex_iterator u_iter, u_end;
|
---|
| 214 | for (tie(u_iter,u_end) = vertices(g); u_iter != u_end; ++u_iter) {
|
---|
| 215 | color[*u_iter] = ColorTraits::white();
|
---|
| 216 | distance[*u_iter] = n;
|
---|
| 217 | }
|
---|
| 218 | color[sink] = ColorTraits::gray();
|
---|
| 219 | distance[sink] = 0;
|
---|
| 220 |
|
---|
| 221 | for (distance_size_type l = 0; l <= max_distance; ++l) {
|
---|
| 222 | layers[l].active_vertices.clear();
|
---|
| 223 | layers[l].inactive_vertices.clear();
|
---|
| 224 | }
|
---|
| 225 |
|
---|
| 226 | max_distance = max_active = 0;
|
---|
| 227 | min_active = n;
|
---|
| 228 |
|
---|
| 229 | Q.push(sink);
|
---|
| 230 | while (! Q.empty()) {
|
---|
| 231 | vertex_descriptor u = Q.top();
|
---|
| 232 | Q.pop();
|
---|
| 233 | distance_size_type d_v = distance[u] + 1;
|
---|
| 234 |
|
---|
| 235 | out_edge_iterator ai, a_end;
|
---|
| 236 | for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) {
|
---|
| 237 | edge_descriptor a = *ai;
|
---|
| 238 | vertex_descriptor v = target(a, g);
|
---|
| 239 | if (color[v] == ColorTraits::white()
|
---|
| 240 | && is_residual_edge(reverse_edge[a])) {
|
---|
| 241 | distance[v] = d_v;
|
---|
| 242 | color[v] = ColorTraits::gray();
|
---|
| 243 | current[v] = out_edges(v, g).first;
|
---|
| 244 | max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(d_v, max_distance);
|
---|
| 245 |
|
---|
| 246 | if (excess_flow[v] > 0)
|
---|
| 247 | add_to_active_list(v, layers[d_v]);
|
---|
| 248 | else
|
---|
| 249 | add_to_inactive_list(v, layers[d_v]);
|
---|
| 250 |
|
---|
| 251 | Q.push(v);
|
---|
| 252 | }
|
---|
| 253 | }
|
---|
| 254 | }
|
---|
| 255 | } // global_distance_update()
|
---|
| 256 |
|
---|
| 257 | //=======================================================================
|
---|
| 258 | // This function is called "push" in Goldberg's h_prf implementation,
|
---|
| 259 | // but it is called "discharge" in the paper and in hi_pr.c.
|
---|
| 260 | void discharge(vertex_descriptor u)
|
---|
| 261 | {
|
---|
| 262 | assert(excess_flow[u] > 0);
|
---|
| 263 | while (1) {
|
---|
| 264 | out_edge_iterator ai, ai_end;
|
---|
| 265 | for (ai = current[u], ai_end = out_edges(u, g).second;
|
---|
| 266 | ai != ai_end; ++ai) {
|
---|
| 267 | edge_descriptor a = *ai;
|
---|
| 268 | if (is_residual_edge(a)) {
|
---|
| 269 | vertex_descriptor v = target(a, g);
|
---|
| 270 | if (is_admissible(u, v)) {
|
---|
| 271 | ++push_count;
|
---|
| 272 | if (v != sink && excess_flow[v] == 0) {
|
---|
| 273 | remove_from_inactive_list(v);
|
---|
| 274 | add_to_active_list(v, layers[distance[v]]);
|
---|
| 275 | }
|
---|
| 276 | push_flow(a);
|
---|
| 277 | if (excess_flow[u] == 0)
|
---|
| 278 | break;
|
---|
| 279 | }
|
---|
| 280 | }
|
---|
| 281 | } // for out_edges of i starting from current
|
---|
| 282 |
|
---|
| 283 | Layer& layer = layers[distance[u]];
|
---|
| 284 | distance_size_type du = distance[u];
|
---|
| 285 |
|
---|
| 286 | if (ai == ai_end) { // i must be relabeled
|
---|
| 287 | relabel_distance(u);
|
---|
| 288 | if (layer.active_vertices.empty()
|
---|
| 289 | && layer.inactive_vertices.empty())
|
---|
| 290 | gap(du);
|
---|
| 291 | if (distance[u] == n)
|
---|
| 292 | break;
|
---|
| 293 | } else { // i is no longer active
|
---|
| 294 | current[u] = ai;
|
---|
| 295 | add_to_inactive_list(u, layer);
|
---|
| 296 | break;
|
---|
| 297 | }
|
---|
| 298 | } // while (1)
|
---|
| 299 | } // discharge()
|
---|
| 300 |
|
---|
| 301 | //=======================================================================
|
---|
| 302 | // This corresponds to the "push" update operation of the paper,
|
---|
| 303 | // not the "push" function in Goldberg's h_prf.c implementation.
|
---|
| 304 | // The idea is to push the excess flow from from vertex u to v.
|
---|
| 305 | void push_flow(edge_descriptor u_v)
|
---|
| 306 | {
|
---|
| 307 | vertex_descriptor
|
---|
| 308 | u = source(u_v, g),
|
---|
| 309 | v = target(u_v, g);
|
---|
| 310 |
|
---|
| 311 | BOOST_USING_STD_MIN();
|
---|
| 312 | FlowValue flow_delta
|
---|
| 313 | = min BOOST_PREVENT_MACRO_SUBSTITUTION(excess_flow[u], residual_capacity[u_v]);
|
---|
| 314 |
|
---|
| 315 | residual_capacity[u_v] -= flow_delta;
|
---|
| 316 | residual_capacity[reverse_edge[u_v]] += flow_delta;
|
---|
| 317 |
|
---|
| 318 | excess_flow[u] -= flow_delta;
|
---|
| 319 | excess_flow[v] += flow_delta;
|
---|
| 320 | } // push_flow()
|
---|
| 321 |
|
---|
| 322 | //=======================================================================
|
---|
| 323 | // The main purpose of this routine is to set distance[v]
|
---|
| 324 | // to the smallest value allowed by the valid labeling constraints,
|
---|
| 325 | // which are:
|
---|
| 326 | // distance[t] = 0
|
---|
| 327 | // distance[u] <= distance[v] + 1 for every residual edge (u,v)
|
---|
| 328 | //
|
---|
| 329 | distance_size_type relabel_distance(vertex_descriptor u)
|
---|
| 330 | {
|
---|
| 331 | BOOST_USING_STD_MAX();
|
---|
| 332 | ++relabel_count;
|
---|
| 333 | work_since_last_update += beta();
|
---|
| 334 |
|
---|
| 335 | distance_size_type min_distance = num_vertices(g);
|
---|
| 336 | distance[u] = min_distance;
|
---|
| 337 |
|
---|
| 338 | // Examine the residual out-edges of vertex i, choosing the
|
---|
| 339 | // edge whose target vertex has the minimal distance.
|
---|
| 340 | out_edge_iterator ai, a_end, min_edge_iter;
|
---|
| 341 | for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) {
|
---|
| 342 | ++work_since_last_update;
|
---|
| 343 | edge_descriptor a = *ai;
|
---|
| 344 | vertex_descriptor v = target(a, g);
|
---|
| 345 | if (is_residual_edge(a) && distance[v] < min_distance) {
|
---|
| 346 | min_distance = distance[v];
|
---|
| 347 | min_edge_iter = ai;
|
---|
| 348 | }
|
---|
| 349 | }
|
---|
| 350 | ++min_distance;
|
---|
| 351 | if (min_distance < n) {
|
---|
| 352 | distance[u] = min_distance; // this is the main action
|
---|
| 353 | current[u] = min_edge_iter;
|
---|
| 354 | max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(min_distance, max_distance);
|
---|
| 355 | }
|
---|
| 356 | return min_distance;
|
---|
| 357 | } // relabel_distance()
|
---|
| 358 |
|
---|
| 359 | //=======================================================================
|
---|
| 360 | // cleanup beyond the gap
|
---|
| 361 | void gap(distance_size_type empty_distance)
|
---|
| 362 | {
|
---|
| 363 | ++gap_count;
|
---|
| 364 |
|
---|
| 365 | distance_size_type r; // distance of layer before the current layer
|
---|
| 366 | r = empty_distance - 1;
|
---|
| 367 |
|
---|
| 368 | // Set the distance for the vertices beyond the gap to "infinity".
|
---|
| 369 | for (layer_iterator l = layers.begin() + empty_distance + 1;
|
---|
| 370 | l < layers.begin() + max_distance; ++l) {
|
---|
| 371 | list_iterator i;
|
---|
| 372 | for (i = l->inactive_vertices.begin();
|
---|
| 373 | i != l->inactive_vertices.end(); ++i) {
|
---|
| 374 | distance[*i] = n;
|
---|
| 375 | ++gap_node_count;
|
---|
| 376 | }
|
---|
| 377 | l->inactive_vertices.clear();
|
---|
| 378 | }
|
---|
| 379 | max_distance = r;
|
---|
| 380 | max_active = r;
|
---|
| 381 | }
|
---|
| 382 |
|
---|
| 383 | //=======================================================================
|
---|
| 384 | // This is the core part of the algorithm, "phase one".
|
---|
| 385 | FlowValue maximum_preflow()
|
---|
| 386 | {
|
---|
| 387 | work_since_last_update = 0;
|
---|
| 388 |
|
---|
| 389 | while (max_active >= min_active) { // "main" loop
|
---|
| 390 |
|
---|
| 391 | Layer& layer = layers[max_active];
|
---|
| 392 | list_iterator u_iter = layer.active_vertices.begin();
|
---|
| 393 |
|
---|
| 394 | if (u_iter == layer.active_vertices.end())
|
---|
| 395 | --max_active;
|
---|
| 396 | else {
|
---|
| 397 | vertex_descriptor u = *u_iter;
|
---|
| 398 | remove_from_active_list(u);
|
---|
| 399 |
|
---|
| 400 | discharge(u);
|
---|
| 401 |
|
---|
| 402 | if (work_since_last_update * global_update_frequency() > nm) {
|
---|
| 403 | global_distance_update();
|
---|
| 404 | work_since_last_update = 0;
|
---|
| 405 | }
|
---|
| 406 | }
|
---|
| 407 | } // while (max_active >= min_active)
|
---|
| 408 |
|
---|
| 409 | return excess_flow[sink];
|
---|
| 410 | } // maximum_preflow()
|
---|
| 411 |
|
---|
| 412 | //=======================================================================
|
---|
| 413 | // remove excess flow, the "second phase"
|
---|
| 414 | // This does a DFS on the reverse flow graph of nodes with excess flow.
|
---|
| 415 | // If a cycle is found, cancel it.
|
---|
| 416 | // Return the nodes with excess flow in topological order.
|
---|
| 417 | //
|
---|
| 418 | // Unlike the prefl_to_flow() implementation, we use
|
---|
| 419 | // "color" instead of "distance" for the DFS labels
|
---|
| 420 | // "parent" instead of nl_prev for the DFS tree
|
---|
| 421 | // "topo_next" instead of nl_next for the topological ordering
|
---|
| 422 | void convert_preflow_to_flow()
|
---|
| 423 | {
|
---|
| 424 | vertex_iterator u_iter, u_end;
|
---|
| 425 | out_edge_iterator ai, a_end;
|
---|
| 426 |
|
---|
| 427 | vertex_descriptor r, restart, u;
|
---|
| 428 |
|
---|
| 429 | std::vector<vertex_descriptor> parent(n);
|
---|
| 430 | std::vector<vertex_descriptor> topo_next(n);
|
---|
| 431 |
|
---|
| 432 | vertex_descriptor tos(parent[0]),
|
---|
| 433 | bos(parent[0]); // bogus initialization, just to avoid warning
|
---|
| 434 | bool bos_null = true;
|
---|
| 435 |
|
---|
| 436 | // handle self-loops
|
---|
| 437 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
|
---|
| 438 | for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai)
|
---|
| 439 | if (target(*ai, g) == *u_iter)
|
---|
| 440 | residual_capacity[*ai] = capacity[*ai];
|
---|
| 441 |
|
---|
| 442 | // initialize
|
---|
| 443 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
|
---|
| 444 | u = *u_iter;
|
---|
| 445 | color[u] = ColorTraits::white();
|
---|
| 446 | parent[u] = u;
|
---|
| 447 | current[u] = out_edges(u, g).first;
|
---|
| 448 | }
|
---|
| 449 | // eliminate flow cycles and topologically order the vertices
|
---|
| 450 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
|
---|
| 451 | u = *u_iter;
|
---|
| 452 | if (color[u] == ColorTraits::white()
|
---|
| 453 | && excess_flow[u] > 0
|
---|
| 454 | && u != src && u != sink ) {
|
---|
| 455 | r = u;
|
---|
| 456 | color[r] = ColorTraits::gray();
|
---|
| 457 | while (1) {
|
---|
| 458 | for (; current[u] != out_edges(u, g).second; ++current[u]) {
|
---|
| 459 | edge_descriptor a = *current[u];
|
---|
| 460 | if (capacity[a] == 0 && is_residual_edge(a)) {
|
---|
| 461 | vertex_descriptor v = target(a, g);
|
---|
| 462 | if (color[v] == ColorTraits::white()) {
|
---|
| 463 | color[v] = ColorTraits::gray();
|
---|
| 464 | parent[v] = u;
|
---|
| 465 | u = v;
|
---|
| 466 | break;
|
---|
| 467 | } else if (color[v] == ColorTraits::gray()) {
|
---|
| 468 | // find minimum flow on the cycle
|
---|
| 469 | FlowValue delta = residual_capacity[a];
|
---|
| 470 | while (1) {
|
---|
| 471 | BOOST_USING_STD_MIN();
|
---|
| 472 | delta = min BOOST_PREVENT_MACRO_SUBSTITUTION(delta, residual_capacity[*current[v]]);
|
---|
| 473 | if (v == u)
|
---|
| 474 | break;
|
---|
| 475 | else
|
---|
| 476 | v = target(*current[v], g);
|
---|
| 477 | }
|
---|
| 478 | // remove delta flow units
|
---|
| 479 | v = u;
|
---|
| 480 | while (1) {
|
---|
| 481 | a = *current[v];
|
---|
| 482 | residual_capacity[a] -= delta;
|
---|
| 483 | residual_capacity[reverse_edge[a]] += delta;
|
---|
| 484 | v = target(a, g);
|
---|
| 485 | if (v == u)
|
---|
| 486 | break;
|
---|
| 487 | }
|
---|
| 488 |
|
---|
| 489 | // back-out of DFS to the first saturated edge
|
---|
| 490 | restart = u;
|
---|
| 491 | for (v = target(*current[u], g); v != u; v = target(a, g)){
|
---|
| 492 | a = *current[v];
|
---|
| 493 | if (color[v] == ColorTraits::white()
|
---|
| 494 | || is_saturated(a)) {
|
---|
| 495 | color[target(*current[v], g)] = ColorTraits::white();
|
---|
| 496 | if (color[v] != ColorTraits::white())
|
---|
| 497 | restart = v;
|
---|
| 498 | }
|
---|
| 499 | }
|
---|
| 500 | if (restart != u) {
|
---|
| 501 | u = restart;
|
---|
| 502 | ++current[u];
|
---|
| 503 | break;
|
---|
| 504 | }
|
---|
| 505 | } // else if (color[v] == ColorTraits::gray())
|
---|
| 506 | } // if (capacity[a] == 0 ...
|
---|
| 507 | } // for out_edges(u, g) (though "u" changes during loop)
|
---|
| 508 |
|
---|
| 509 | if (current[u] == out_edges(u, g).second) {
|
---|
| 510 | // scan of i is complete
|
---|
| 511 | color[u] = ColorTraits::black();
|
---|
| 512 | if (u != src) {
|
---|
| 513 | if (bos_null) {
|
---|
| 514 | bos = u;
|
---|
| 515 | bos_null = false;
|
---|
| 516 | tos = u;
|
---|
| 517 | } else {
|
---|
| 518 | topo_next[u] = tos;
|
---|
| 519 | tos = u;
|
---|
| 520 | }
|
---|
| 521 | }
|
---|
| 522 | if (u != r) {
|
---|
| 523 | u = parent[u];
|
---|
| 524 | ++current[u];
|
---|
| 525 | } else
|
---|
| 526 | break;
|
---|
| 527 | }
|
---|
| 528 | } // while (1)
|
---|
| 529 | } // if (color[u] == white && excess_flow[u] > 0 & ...)
|
---|
| 530 | } // for all vertices in g
|
---|
| 531 |
|
---|
| 532 | // return excess flows
|
---|
| 533 | // note that the sink is not on the stack
|
---|
| 534 | if (! bos_null) {
|
---|
| 535 | for (u = tos; u != bos; u = topo_next[u]) {
|
---|
| 536 | ai = out_edges(u, g).first;
|
---|
| 537 | while (excess_flow[u] > 0 && ai != out_edges(u, g).second) {
|
---|
| 538 | if (capacity[*ai] == 0 && is_residual_edge(*ai))
|
---|
| 539 | push_flow(*ai);
|
---|
| 540 | ++ai;
|
---|
| 541 | }
|
---|
| 542 | }
|
---|
| 543 | // do the bottom
|
---|
| 544 | u = bos;
|
---|
| 545 | ai = out_edges(u, g).first;
|
---|
| 546 | while (excess_flow[u] > 0) {
|
---|
| 547 | if (capacity[*ai] == 0 && is_residual_edge(*ai))
|
---|
| 548 | push_flow(*ai);
|
---|
| 549 | ++ai;
|
---|
| 550 | }
|
---|
| 551 | }
|
---|
| 552 |
|
---|
| 553 | } // convert_preflow_to_flow()
|
---|
| 554 |
|
---|
| 555 | //=======================================================================
|
---|
| 556 | inline bool is_flow()
|
---|
| 557 | {
|
---|
| 558 | vertex_iterator u_iter, u_end;
|
---|
| 559 | out_edge_iterator ai, a_end;
|
---|
| 560 |
|
---|
| 561 | // check edge flow values
|
---|
| 562 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
|
---|
| 563 | for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) {
|
---|
| 564 | edge_descriptor a = *ai;
|
---|
| 565 | if (capacity[a] > 0)
|
---|
| 566 | if ((residual_capacity[a] + residual_capacity[reverse_edge[a]]
|
---|
| 567 | != capacity[a] + capacity[reverse_edge[a]])
|
---|
| 568 | || (residual_capacity[a] < 0)
|
---|
| 569 | || (residual_capacity[reverse_edge[a]] < 0))
|
---|
| 570 | return false;
|
---|
| 571 | }
|
---|
| 572 | }
|
---|
| 573 |
|
---|
| 574 | // check conservation
|
---|
| 575 | FlowValue sum;
|
---|
| 576 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
|
---|
| 577 | vertex_descriptor u = *u_iter;
|
---|
| 578 | if (u != src && u != sink) {
|
---|
| 579 | if (excess_flow[u] != 0)
|
---|
| 580 | return false;
|
---|
| 581 | sum = 0;
|
---|
| 582 | for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai)
|
---|
| 583 | if (capacity[*ai] > 0)
|
---|
| 584 | sum -= capacity[*ai] - residual_capacity[*ai];
|
---|
| 585 | else
|
---|
| 586 | sum += residual_capacity[*ai];
|
---|
| 587 |
|
---|
| 588 | if (excess_flow[u] != sum)
|
---|
| 589 | return false;
|
---|
| 590 | }
|
---|
| 591 | }
|
---|
| 592 |
|
---|
| 593 | return true;
|
---|
| 594 | } // is_flow()
|
---|
| 595 |
|
---|
| 596 | bool is_optimal() {
|
---|
| 597 | // check if mincut is saturated...
|
---|
| 598 | global_distance_update();
|
---|
| 599 | return distance[src] >= n;
|
---|
| 600 | }
|
---|
| 601 |
|
---|
| 602 | void print_statistics(std::ostream& os) const {
|
---|
| 603 | os << "pushes: " << push_count << std::endl
|
---|
| 604 | << "relabels: " << relabel_count << std::endl
|
---|
| 605 | << "updates: " << update_count << std::endl
|
---|
| 606 | << "gaps: " << gap_count << std::endl
|
---|
| 607 | << "gap nodes: " << gap_node_count << std::endl
|
---|
| 608 | << std::endl;
|
---|
| 609 | }
|
---|
| 610 |
|
---|
| 611 | void print_flow_values(std::ostream& os) const {
|
---|
| 612 | os << "flow values" << std::endl;
|
---|
| 613 | vertex_iterator u_iter, u_end;
|
---|
| 614 | out_edge_iterator ei, e_end;
|
---|
| 615 | for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
|
---|
| 616 | for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei)
|
---|
| 617 | if (capacity[*ei] > 0)
|
---|
| 618 | os << *u_iter << " " << target(*ei, g) << " "
|
---|
| 619 | << (capacity[*ei] - residual_capacity[*ei]) << std::endl;
|
---|
| 620 | os << std::endl;
|
---|
| 621 | }
|
---|
| 622 |
|
---|
| 623 | //=======================================================================
|
---|
| 624 |
|
---|
| 625 | Graph& g;
|
---|
| 626 | vertices_size_type n;
|
---|
| 627 | vertices_size_type nm;
|
---|
| 628 | EdgeCapacityMap capacity;
|
---|
| 629 | vertex_descriptor src;
|
---|
| 630 | vertex_descriptor sink;
|
---|
| 631 | VertexIndexMap index;
|
---|
| 632 |
|
---|
| 633 | // will need to use random_access_property_map with these
|
---|
| 634 | std::vector< FlowValue > excess_flow;
|
---|
| 635 | std::vector< out_edge_iterator > current;
|
---|
| 636 | std::vector< distance_size_type > distance;
|
---|
| 637 | std::vector< default_color_type > color;
|
---|
| 638 |
|
---|
| 639 | // Edge Property Maps that must be interior to the graph
|
---|
| 640 | ReverseEdgeMap reverse_edge;
|
---|
| 641 | ResidualCapacityEdgeMap residual_capacity;
|
---|
| 642 |
|
---|
| 643 | LayerArray layers;
|
---|
| 644 | std::vector< list_iterator > layer_list_ptr;
|
---|
| 645 | distance_size_type max_distance; // maximal distance
|
---|
| 646 | distance_size_type max_active; // maximal distance with active node
|
---|
| 647 | distance_size_type min_active; // minimal distance with active node
|
---|
| 648 | boost::queue<vertex_descriptor> Q;
|
---|
| 649 |
|
---|
| 650 | // Statistics counters
|
---|
| 651 | long push_count;
|
---|
| 652 | long update_count;
|
---|
| 653 | long relabel_count;
|
---|
| 654 | long gap_count;
|
---|
| 655 | long gap_node_count;
|
---|
| 656 |
|
---|
| 657 | inline double global_update_frequency() { return 0.5; }
|
---|
| 658 | inline vertices_size_type alpha() { return 6; }
|
---|
| 659 | inline long beta() { return 12; }
|
---|
| 660 |
|
---|
| 661 | long work_since_last_update;
|
---|
| 662 | };
|
---|
| 663 |
|
---|
| 664 | } // namespace detail
|
---|
| 665 |
|
---|
| 666 | template <class Graph,
|
---|
| 667 | class CapacityEdgeMap, class ResidualCapacityEdgeMap,
|
---|
| 668 | class ReverseEdgeMap, class VertexIndexMap>
|
---|
| 669 | typename property_traits<CapacityEdgeMap>::value_type
|
---|
| 670 | push_relabel_max_flow
|
---|
| 671 | (Graph& g,
|
---|
| 672 | typename graph_traits<Graph>::vertex_descriptor src,
|
---|
| 673 | typename graph_traits<Graph>::vertex_descriptor sink,
|
---|
| 674 | CapacityEdgeMap cap, ResidualCapacityEdgeMap res,
|
---|
| 675 | ReverseEdgeMap rev, VertexIndexMap index_map)
|
---|
| 676 | {
|
---|
| 677 | typedef typename property_traits<CapacityEdgeMap>::value_type FlowValue;
|
---|
| 678 |
|
---|
| 679 | detail::push_relabel<Graph, CapacityEdgeMap, ResidualCapacityEdgeMap,
|
---|
| 680 | ReverseEdgeMap, VertexIndexMap, FlowValue>
|
---|
| 681 | algo(g, cap, res, rev, src, sink, index_map);
|
---|
| 682 |
|
---|
| 683 | FlowValue flow = algo.maximum_preflow();
|
---|
| 684 |
|
---|
| 685 | algo.convert_preflow_to_flow();
|
---|
| 686 |
|
---|
| 687 | assert(algo.is_flow());
|
---|
| 688 | assert(algo.is_optimal());
|
---|
| 689 |
|
---|
| 690 | return flow;
|
---|
| 691 | } // push_relabel_max_flow()
|
---|
| 692 |
|
---|
| 693 | template <class Graph, class P, class T, class R>
|
---|
| 694 | typename detail::edge_capacity_value<Graph, P, T, R>::type
|
---|
| 695 | push_relabel_max_flow
|
---|
| 696 | (Graph& g,
|
---|
| 697 | typename graph_traits<Graph>::vertex_descriptor src,
|
---|
| 698 | typename graph_traits<Graph>::vertex_descriptor sink,
|
---|
| 699 | const bgl_named_params<P, T, R>& params)
|
---|
| 700 | {
|
---|
| 701 | return push_relabel_max_flow
|
---|
| 702 | (g, src, sink,
|
---|
| 703 | choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity),
|
---|
| 704 | choose_pmap(get_param(params, edge_residual_capacity),
|
---|
| 705 | g, edge_residual_capacity),
|
---|
| 706 | choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse),
|
---|
| 707 | choose_const_pmap(get_param(params, vertex_index), g, vertex_index)
|
---|
| 708 | );
|
---|
| 709 | }
|
---|
| 710 |
|
---|
| 711 | template <class Graph>
|
---|
| 712 | typename property_traits<
|
---|
| 713 | typename property_map<Graph, edge_capacity_t>::const_type
|
---|
| 714 | >::value_type
|
---|
| 715 | push_relabel_max_flow
|
---|
| 716 | (Graph& g,
|
---|
| 717 | typename graph_traits<Graph>::vertex_descriptor src,
|
---|
| 718 | typename graph_traits<Graph>::vertex_descriptor sink)
|
---|
| 719 | {
|
---|
| 720 | bgl_named_params<int, buffer_param_t> params(0); // bogus empty param
|
---|
| 721 | return push_relabel_max_flow(g, src, sink, params);
|
---|
| 722 | }
|
---|
| 723 |
|
---|
| 724 | } // namespace boost
|
---|
| 725 |
|
---|
| 726 | #endif // BOOST_PUSH_RELABEL_MAX_FLOW_HPP
|
---|
| 727 |
|
---|