git.sesse.net Git - casparcg/blob - dependencies64/boost/boost/graph/push_relabel_max_flow.hpp

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