Chris@16: //=======================================================================
Chris@16: // Copyright 2000 University of Notre Dame.
Chris@16: // Authors: Jeremy G. Siek, Andrew Lumsdaine, Lie-Quan Lee
Chris@16: //
Chris@16: // Distributed under the Boost Software License, Version 1.0. (See
Chris@16: // accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: //=======================================================================
Chris@16: 
Chris@16: #ifndef BOOST_PUSH_RELABEL_MAX_FLOW_HPP
Chris@16: #define BOOST_PUSH_RELABEL_MAX_FLOW_HPP
Chris@16: 
Chris@16: #include <boost/config.hpp>
Chris@16: #include <boost/assert.hpp>
Chris@16: #include <vector>
Chris@16: #include <list>
Chris@16: #include <iosfwd>
Chris@16: #include <algorithm> // for std::min and std::max
Chris@16: 
Chris@16: #include <boost/pending/queue.hpp>
Chris@16: #include <boost/limits.hpp>
Chris@16: #include <boost/graph/graph_concepts.hpp>
Chris@16: #include <boost/graph/named_function_params.hpp>
Chris@16: 
Chris@16: namespace boost {
Chris@16: 
Chris@16:   namespace detail {
Chris@16:     
Chris@16:    // This implementation is based on Goldberg's 
Chris@16:    // "On Implementing Push-Relabel Method for the Maximum Flow Problem"
Chris@16:    // by B.V. Cherkassky and A.V. Goldberg, IPCO '95, pp. 157--171
Chris@16:    // and on the h_prf.c and hi_pr.c code written by the above authors.
Chris@16: 
Chris@16:    // This implements the highest-label version of the push-relabel method
Chris@16:    // with the global relabeling and gap relabeling heuristics.
Chris@16: 
Chris@16:    // The terms "rank", "distance", "height" are synonyms in
Chris@16:    // Goldberg's implementation, paper and in the CLR.  A "layer" is a
Chris@16:    // group of vertices with the same distance. The vertices in each
Chris@16:    // layer are categorized as active or inactive.  An active vertex
Chris@16:    // has positive excess flow and its distance is less than n (it is
Chris@16:    // not blocked).
Chris@16: 
Chris@16:     template <class Vertex>
Chris@16:     struct preflow_layer {
Chris@16:       std::list<Vertex> active_vertices;
Chris@16:       std::list<Vertex> inactive_vertices;
Chris@16:     };
Chris@16: 
Chris@16:     template <class Graph, 
Chris@16:               class EdgeCapacityMap,    // integer value type
Chris@16:               class ResidualCapacityEdgeMap,
Chris@16:               class ReverseEdgeMap,
Chris@16:               class VertexIndexMap,     // vertex_descriptor -> integer
Chris@16:               class FlowValue>
Chris@16:     class push_relabel
Chris@16:     {
Chris@16:     public:
Chris@16:       typedef graph_traits<Graph> Traits;
Chris@16:       typedef typename Traits::vertex_descriptor vertex_descriptor;
Chris@16:       typedef typename Traits::edge_descriptor edge_descriptor;
Chris@16:       typedef typename Traits::vertex_iterator vertex_iterator;
Chris@16:       typedef typename Traits::out_edge_iterator out_edge_iterator;
Chris@16:       typedef typename Traits::vertices_size_type vertices_size_type;
Chris@16:       typedef typename Traits::edges_size_type edges_size_type;
Chris@16: 
Chris@16:       typedef preflow_layer<vertex_descriptor> Layer;
Chris@16:       typedef std::vector< Layer > LayerArray;
Chris@16:       typedef typename LayerArray::iterator layer_iterator;
Chris@16:       typedef typename LayerArray::size_type distance_size_type;
Chris@16: 
Chris@16:       typedef color_traits<default_color_type> ColorTraits;
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       // Some helper predicates
Chris@16: 
Chris@16:       inline bool is_admissible(vertex_descriptor u, vertex_descriptor v) {
Chris@16:         return get(distance, u) == get(distance, v) + 1;
Chris@16:       }
Chris@16:       inline bool is_residual_edge(edge_descriptor a) {
Chris@16:         return 0 < get(residual_capacity, a);
Chris@16:       }
Chris@16:       inline bool is_saturated(edge_descriptor a) {
Chris@16:         return get(residual_capacity, a) == 0;
Chris@16:       }
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       // Layer List Management Functions
Chris@16: 
Chris@16:       typedef typename std::list<vertex_descriptor>::iterator list_iterator;
Chris@16: 
Chris@16:       void add_to_active_list(vertex_descriptor u, Layer& layer) {
Chris@16:         BOOST_USING_STD_MIN();
Chris@16:         BOOST_USING_STD_MAX();
Chris@16:         layer.active_vertices.push_front(u);
Chris@16:         max_active = max BOOST_PREVENT_MACRO_SUBSTITUTION(get(distance, u), max_active);
Chris@16:         min_active = min BOOST_PREVENT_MACRO_SUBSTITUTION(get(distance, u), min_active);
Chris@16:         layer_list_ptr[u] = layer.active_vertices.begin();
Chris@16:       }
Chris@16:       void remove_from_active_list(vertex_descriptor u) {
Chris@16:         layers[get(distance, u)].active_vertices.erase(layer_list_ptr[u]);    
Chris@16:       }
Chris@16: 
Chris@16:       void add_to_inactive_list(vertex_descriptor u, Layer& layer) {
Chris@16:         layer.inactive_vertices.push_front(u);
Chris@16:         layer_list_ptr[u] = layer.inactive_vertices.begin();
Chris@16:       }
Chris@16:       void remove_from_inactive_list(vertex_descriptor u) {
Chris@16:         layers[get(distance, u)].inactive_vertices.erase(layer_list_ptr[u]);    
Chris@16:       }
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       // initialization
Chris@16:       push_relabel(Graph& g_, 
Chris@16:                    EdgeCapacityMap cap,
Chris@16:                    ResidualCapacityEdgeMap res,
Chris@16:                    ReverseEdgeMap rev,
Chris@16:                    vertex_descriptor src_, 
Chris@16:                    vertex_descriptor sink_,
Chris@16:                    VertexIndexMap idx)
Chris@16:         : g(g_), n(num_vertices(g_)), capacity(cap), src(src_), sink(sink_), 
Chris@16:           index(idx),
Chris@16:           excess_flow_data(num_vertices(g_)),
Chris@16:           excess_flow(excess_flow_data.begin(), idx),
Chris@16:           current_data(num_vertices(g_), out_edges(*vertices(g_).first, g_)),
Chris@16:           current(current_data.begin(), idx),
Chris@16:           distance_data(num_vertices(g_)),
Chris@16:           distance(distance_data.begin(), idx),
Chris@16:           color_data(num_vertices(g_)),
Chris@16:           color(color_data.begin(), idx),
Chris@16:           reverse_edge(rev),
Chris@16:           residual_capacity(res),
Chris@16:           layers(num_vertices(g_)),
Chris@16:           layer_list_ptr_data(num_vertices(g_), 
Chris@16:                               layers.front().inactive_vertices.end()),
Chris@16:           layer_list_ptr(layer_list_ptr_data.begin(), idx),
Chris@16:           push_count(0), update_count(0), relabel_count(0), 
Chris@16:           gap_count(0), gap_node_count(0),
Chris@16:           work_since_last_update(0)
Chris@16:       {
Chris@16:         vertex_iterator u_iter, u_end;
Chris@16:         // Don't count the reverse edges
Chris@16:         edges_size_type m = num_edges(g) / 2;
Chris@16:         nm = alpha() * n + m;
Chris@16: 
Chris@16:         // Initialize flow to zero which means initializing
Chris@16:         // the residual capacity to equal the capacity.
Chris@16:         out_edge_iterator ei, e_end;
Chris@16:         for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
Chris@16:           for (boost::tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) {
Chris@16:             put(residual_capacity, *ei, get(capacity, *ei));
Chris@16:           }
Chris@16: 
Chris@16:         for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
Chris@16:           vertex_descriptor u = *u_iter;
Chris@16:           put(excess_flow, u, 0);
Chris@16:           current[u] = out_edges(u, g);
Chris@16:         }
Chris@16: 
Chris@16:         bool overflow_detected = false;
Chris@16:         FlowValue test_excess = 0;
Chris@16: 
Chris@16:         out_edge_iterator a_iter, a_end;
Chris@16:         for (boost::tie(a_iter, a_end) = out_edges(src, g); a_iter != a_end; ++a_iter)
Chris@16:           if (target(*a_iter, g) != src)
Chris@16:             test_excess += get(residual_capacity, *a_iter);
Chris@16:         if (test_excess > (std::numeric_limits<FlowValue>::max)())
Chris@16:           overflow_detected = true;
Chris@16: 
Chris@16:         if (overflow_detected)
Chris@16:           put(excess_flow, src, (std::numeric_limits<FlowValue>::max)());
Chris@16:         else {
Chris@16:           put(excess_flow, src, 0);
Chris@16:           for (boost::tie(a_iter, a_end) = out_edges(src, g); 
Chris@16:                a_iter != a_end; ++a_iter) {
Chris@16:             edge_descriptor a = *a_iter;
Chris@16:             vertex_descriptor tgt = target(a, g);
Chris@16:             if (tgt != src) {
Chris@16:               ++push_count;
Chris@16:               FlowValue delta = get(residual_capacity, a);
Chris@16:               put(residual_capacity, a, get(residual_capacity, a) - delta);
Chris@16:               edge_descriptor rev = get(reverse_edge, a);
Chris@16:               put(residual_capacity, rev, get(residual_capacity, rev) + delta);
Chris@16:               put(excess_flow, tgt, get(excess_flow, tgt) + delta);
Chris@16:             }
Chris@16:           }
Chris@16:         }
Chris@16:         max_distance = num_vertices(g) - 1;
Chris@16:         max_active = 0;
Chris@16:         min_active = n;
Chris@16: 
Chris@16:         for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
Chris@16:           vertex_descriptor u = *u_iter;
Chris@16:           if (u == sink) {
Chris@16:             put(distance, u, 0);
Chris@16:             continue;
Chris@16:           } else if (u == src && !overflow_detected)
Chris@16:             put(distance, u, n);
Chris@16:           else
Chris@16:             put(distance, u, 1);
Chris@16: 
Chris@16:           if (get(excess_flow, u) > 0)
Chris@16:             add_to_active_list(u, layers[1]);
Chris@16:           else if (get(distance, u) < n)
Chris@16:             add_to_inactive_list(u, layers[1]);
Chris@16:         }       
Chris@16: 
Chris@16:       } // push_relabel constructor
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       // This is a breadth-first search over the residual graph
Chris@16:       // (well, actually the reverse of the residual graph).
Chris@16:       // Would be cool to have a graph view adaptor for hiding certain
Chris@16:       // edges, like the saturated (non-residual) edges in this case.
Chris@16:       // Goldberg's implementation abused "distance" for the coloring.
Chris@16:       void global_distance_update()
Chris@16:       {
Chris@16:         BOOST_USING_STD_MAX();
Chris@16:         ++update_count;
Chris@16:         vertex_iterator u_iter, u_end;
Chris@16:         for (boost::tie(u_iter,u_end) = vertices(g); u_iter != u_end; ++u_iter) {
Chris@16:           put(color, *u_iter, ColorTraits::white());
Chris@16:           put(distance, *u_iter, n);
Chris@16:         }
Chris@16:         put(color, sink, ColorTraits::gray());
Chris@16:         put(distance, sink, 0);
Chris@16:         
Chris@16:         for (distance_size_type l = 0; l <= max_distance; ++l) {
Chris@16:           layers[l].active_vertices.clear();
Chris@16:           layers[l].inactive_vertices.clear();
Chris@16:         }
Chris@16:         
Chris@16:         max_distance = max_active = 0;
Chris@16:         min_active = n;
Chris@16: 
Chris@16:         Q.push(sink);
Chris@16:         while (! Q.empty()) {
Chris@16:           vertex_descriptor u = Q.top();
Chris@16:           Q.pop();
Chris@16:           distance_size_type d_v = get(distance, u) + 1;
Chris@16: 
Chris@16:           out_edge_iterator ai, a_end;
Chris@16:           for (boost::tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) {
Chris@16:             edge_descriptor a = *ai;
Chris@16:             vertex_descriptor v = target(a, g);
Chris@16:             if (get(color, v) == ColorTraits::white()
Chris@16:                 && is_residual_edge(get(reverse_edge, a))) {
Chris@16:               put(distance, v, d_v);
Chris@16:               put(color, v, ColorTraits::gray());
Chris@16:               current[v] = out_edges(v, g);
Chris@16:               max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(d_v, max_distance);
Chris@16: 
Chris@16:               if (get(excess_flow, v) > 0)
Chris@16:                 add_to_active_list(v, layers[d_v]);
Chris@16:               else
Chris@16:                 add_to_inactive_list(v, layers[d_v]);
Chris@16: 
Chris@16:               Q.push(v);
Chris@16:             }
Chris@16:           }
Chris@16:         }
Chris@16:       } // global_distance_update()
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       // This function is called "push" in Goldberg's h_prf implementation,
Chris@16:       // but it is called "discharge" in the paper and in hi_pr.c.
Chris@16:       void discharge(vertex_descriptor u)
Chris@16:       {
Chris@16:         BOOST_ASSERT(get(excess_flow, u) > 0);
Chris@16:         while (1) {
Chris@16:           out_edge_iterator ai, ai_end;
Chris@16:           for (boost::tie(ai, ai_end) = current[u]; ai != ai_end; ++ai) {
Chris@16:             edge_descriptor a = *ai;
Chris@16:             if (is_residual_edge(a)) {
Chris@16:               vertex_descriptor v = target(a, g);
Chris@16:               if (is_admissible(u, v)) {
Chris@16:                 ++push_count;
Chris@16:                 if (v != sink && get(excess_flow, v) == 0) {
Chris@16:                   remove_from_inactive_list(v);
Chris@16:                   add_to_active_list(v, layers[get(distance, v)]);
Chris@16:                 }
Chris@16:                 push_flow(a);
Chris@16:                 if (get(excess_flow, u) == 0)
Chris@16:                   break;
Chris@16:               } 
Chris@16:             } 
Chris@16:           } // for out_edges of i starting from current
Chris@16: 
Chris@16:           Layer& layer = layers[get(distance, u)];
Chris@16:           distance_size_type du = get(distance, u);
Chris@16: 
Chris@16:           if (ai == ai_end) {   // i must be relabeled
Chris@16:             relabel_distance(u);
Chris@16:             if (layer.active_vertices.empty()
Chris@16:                 && layer.inactive_vertices.empty())
Chris@16:               gap(du);
Chris@16:             if (get(distance, u) == n)
Chris@16:               break;
Chris@16:           } else {              // i is no longer active
Chris@16:             current[u].first = ai;
Chris@16:             add_to_inactive_list(u, layer);
Chris@16:             break;
Chris@16:           }
Chris@16:         } // while (1)
Chris@16:       } // discharge()
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       // This corresponds to the "push" update operation of the paper,
Chris@16:       // not the "push" function in Goldberg's h_prf.c implementation.
Chris@16:       // The idea is to push the excess flow from from vertex u to v.
Chris@16:       void push_flow(edge_descriptor u_v)
Chris@16:       {
Chris@16:         vertex_descriptor
Chris@16:           u = source(u_v, g),
Chris@16:           v = target(u_v, g);
Chris@16:         
Chris@16:         BOOST_USING_STD_MIN();
Chris@16:         FlowValue flow_delta
Chris@16:           = min BOOST_PREVENT_MACRO_SUBSTITUTION(get(excess_flow, u), get(residual_capacity, u_v));
Chris@16: 
Chris@16:         put(residual_capacity, u_v, get(residual_capacity, u_v) - flow_delta);
Chris@16:         edge_descriptor rev = get(reverse_edge, u_v);
Chris@16:         put(residual_capacity, rev, get(residual_capacity, rev) + flow_delta);
Chris@16: 
Chris@16:         put(excess_flow, u, get(excess_flow, u) - flow_delta);
Chris@16:         put(excess_flow, v, get(excess_flow, v) + flow_delta);
Chris@16:       } // push_flow()
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       // The main purpose of this routine is to set distance[v]
Chris@16:       // to the smallest value allowed by the valid labeling constraints,
Chris@16:       // which are:
Chris@16:       // distance[t] = 0
Chris@16:       // distance[u] <= distance[v] + 1   for every residual edge (u,v)
Chris@16:       //
Chris@16:       distance_size_type relabel_distance(vertex_descriptor u)
Chris@16:       {
Chris@16:         BOOST_USING_STD_MAX();
Chris@16:         ++relabel_count;
Chris@16:         work_since_last_update += beta();
Chris@16: 
Chris@16:         distance_size_type min_distance = num_vertices(g);
Chris@16:         put(distance, u, min_distance);
Chris@16: 
Chris@16:         // Examine the residual out-edges of vertex i, choosing the
Chris@16:         // edge whose target vertex has the minimal distance.
Chris@16:         out_edge_iterator ai, a_end, min_edge_iter;
Chris@16:         for (boost::tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) {
Chris@16:           ++work_since_last_update;
Chris@16:           edge_descriptor a = *ai;
Chris@16:           vertex_descriptor v = target(a, g);
Chris@16:           if (is_residual_edge(a) && get(distance, v) < min_distance) {
Chris@16:             min_distance = get(distance, v);
Chris@16:             min_edge_iter = ai;
Chris@16:           }
Chris@16:         }
Chris@16:         ++min_distance;
Chris@16:         if (min_distance < n) {
Chris@16:           put(distance, u, min_distance);     // this is the main action
Chris@16:           current[u].first = min_edge_iter;
Chris@16:           max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(min_distance, max_distance);
Chris@16:         }
Chris@16:         return min_distance;
Chris@16:       } // relabel_distance()
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       // cleanup beyond the gap
Chris@16:       void gap(distance_size_type empty_distance)
Chris@16:       {
Chris@16:         ++gap_count;
Chris@16: 
Chris@16:         distance_size_type r; // distance of layer before the current layer
Chris@16:         r = empty_distance - 1;
Chris@16: 
Chris@16:         // Set the distance for the vertices beyond the gap to "infinity".
Chris@16:         for (layer_iterator l = layers.begin() + empty_distance + 1;
Chris@16:              l < layers.begin() + max_distance; ++l) {
Chris@16:           list_iterator i;
Chris@16:           for (i = l->inactive_vertices.begin(); 
Chris@16:                i != l->inactive_vertices.end(); ++i) {
Chris@16:             put(distance, *i, n);
Chris@16:             ++gap_node_count;
Chris@16:           }
Chris@16:           l->inactive_vertices.clear();
Chris@16:         }
Chris@16:         max_distance = r;
Chris@16:         max_active = r;
Chris@16:       }
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       // This is the core part of the algorithm, "phase one".
Chris@16:       FlowValue maximum_preflow()
Chris@16:       {
Chris@16:         work_since_last_update = 0;
Chris@16: 
Chris@16:         while (max_active >= min_active) { // "main" loop
Chris@16: 
Chris@16:           Layer& layer = layers[max_active];
Chris@16:           list_iterator u_iter = layer.active_vertices.begin();
Chris@16: 
Chris@16:           if (u_iter == layer.active_vertices.end())
Chris@16:             --max_active;
Chris@16:           else {
Chris@16:             vertex_descriptor u = *u_iter;
Chris@16:             remove_from_active_list(u);
Chris@16:             
Chris@16:             discharge(u);
Chris@16: 
Chris@16:             if (work_since_last_update * global_update_frequency() > nm) {
Chris@16:               global_distance_update();
Chris@16:               work_since_last_update = 0;
Chris@16:             }
Chris@16:           }
Chris@16:         } // while (max_active >= min_active)
Chris@16: 
Chris@16:         return get(excess_flow, sink);
Chris@16:       } // maximum_preflow()
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       // remove excess flow, the "second phase"
Chris@16:       // This does a DFS on the reverse flow graph of nodes with excess flow.
Chris@16:       // If a cycle is found, cancel it.
Chris@16:       // Return the nodes with excess flow in topological order.
Chris@16:       //
Chris@16:       // Unlike the prefl_to_flow() implementation, we use
Chris@16:       //   "color" instead of "distance" for the DFS labels
Chris@16:       //   "parent" instead of nl_prev for the DFS tree
Chris@16:       //   "topo_next" instead of nl_next for the topological ordering
Chris@16:       void convert_preflow_to_flow()
Chris@16:       {
Chris@16:         vertex_iterator u_iter, u_end;
Chris@16:         out_edge_iterator ai, a_end;
Chris@16: 
Chris@16:         vertex_descriptor r, restart, u;
Chris@16: 
Chris@16:         std::vector<vertex_descriptor> parent(n);
Chris@16:         std::vector<vertex_descriptor> topo_next(n);
Chris@16: 
Chris@16:         vertex_descriptor tos(parent[0]), 
Chris@16:           bos(parent[0]); // bogus initialization, just to avoid warning
Chris@16:         bool bos_null = true;
Chris@16: 
Chris@16:         // handle self-loops
Chris@16:         for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
Chris@16:           for (boost::tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai)
Chris@16:             if (target(*ai, g) == *u_iter)
Chris@16:               put(residual_capacity, *ai, get(capacity, *ai));
Chris@16: 
Chris@16:         // initialize
Chris@16:         for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
Chris@16:           u = *u_iter;
Chris@16:           put(color, u, ColorTraits::white());
Chris@16:           parent[get(index, u)] = u;
Chris@16:           current[u] = out_edges(u, g);
Chris@16:         }
Chris@16:         // eliminate flow cycles and topologically order the vertices
Chris@16:         for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
Chris@16:           u = *u_iter;
Chris@16:           if (get(color, u) == ColorTraits::white() 
Chris@16:               && get(excess_flow, u) > 0
Chris@16:               && u != src && u != sink ) {
Chris@16:             r = u;
Chris@16:             put(color, r, ColorTraits::gray());
Chris@16:             while (1) {
Chris@16:               for (; current[u].first != current[u].second; ++current[u].first) {
Chris@16:                 edge_descriptor a = *current[u].first;
Chris@16:                 if (get(capacity, a) == 0 && is_residual_edge(a)) {
Chris@16:                   vertex_descriptor v = target(a, g);
Chris@16:                   if (get(color, v) == ColorTraits::white()) {
Chris@16:                     put(color, v, ColorTraits::gray());
Chris@16:                     parent[get(index, v)] = u;
Chris@16:                     u = v;
Chris@16:                     break;
Chris@16:                   } else if (get(color, v) == ColorTraits::gray()) {
Chris@16:                     // find minimum flow on the cycle
Chris@16:                     FlowValue delta = get(residual_capacity, a);
Chris@16:                     while (1) {
Chris@16:                       BOOST_USING_STD_MIN();
Chris@16:                       delta = min BOOST_PREVENT_MACRO_SUBSTITUTION(delta, get(residual_capacity, *current[v].first));
Chris@16:                       if (v == u)
Chris@16:                         break;
Chris@16:                       else
Chris@16:                         v = target(*current[v].first, g);
Chris@16:                     }
Chris@16:                     // remove delta flow units
Chris@16:                     v = u;
Chris@16:                     while (1) {
Chris@16:                       a = *current[v].first;
Chris@16:                       put(residual_capacity, a, get(residual_capacity, a) - delta);
Chris@16:                       edge_descriptor rev = get(reverse_edge, a);
Chris@16:                       put(residual_capacity, rev, get(residual_capacity, rev) + delta);
Chris@16:                       v = target(a, g);
Chris@16:                       if (v == u)
Chris@16:                         break;
Chris@16:                     }
Chris@16: 
Chris@16:                     // back-out of DFS to the first saturated edge
Chris@16:                     restart = u;
Chris@16:                     for (v = target(*current[u].first, g); v != u; v = target(a, g)){
Chris@16:                       a = *current[v].first;
Chris@16:                       if (get(color, v) == ColorTraits::white() 
Chris@16:                           || is_saturated(a)) {
Chris@16:                         put(color, target(*current[v].first, g), ColorTraits::white());
Chris@16:                         if (get(color, v) != ColorTraits::white())
Chris@16:                           restart = v;
Chris@16:                       }
Chris@16:                     }
Chris@16:                     if (restart != u) {
Chris@16:                       u = restart;
Chris@16:                       ++current[u].first;
Chris@16:                       break;
Chris@16:                     }
Chris@16:                   } // else if (color[v] == ColorTraits::gray())
Chris@16:                 } // if (get(capacity, a) == 0 ...
Chris@16:               } // for out_edges(u, g)  (though "u" changes during loop)
Chris@16:               
Chris@16:               if ( current[u].first == current[u].second ) {
Chris@16:                 // scan of i is complete
Chris@16:                 put(color, u, ColorTraits::black());
Chris@16:                 if (u != src) {
Chris@16:                   if (bos_null) {
Chris@16:                     bos = u;
Chris@16:                     bos_null = false;
Chris@16:                     tos = u;
Chris@16:                   } else {
Chris@16:                     topo_next[get(index, u)] = tos;
Chris@16:                     tos = u;
Chris@16:                   }
Chris@16:                 }
Chris@16:                 if (u != r) {
Chris@16:                   u = parent[get(index, u)];
Chris@16:                   ++current[u].first;
Chris@16:                 } else
Chris@16:                   break;
Chris@16:               }
Chris@16:             } // while (1)
Chris@16:           } // if (color[u] == white && excess_flow[u] > 0 & ...)
Chris@16:         } // for all vertices in g
Chris@16: 
Chris@16:         // return excess flows
Chris@16:         // note that the sink is not on the stack
Chris@16:         if (! bos_null) {
Chris@16:           for (u = tos; u != bos; u = topo_next[get(index, u)]) {
Chris@16:             boost::tie(ai, a_end) = out_edges(u, g);
Chris@16:             while (get(excess_flow, u) > 0 && ai != a_end) {
Chris@16:               if (get(capacity, *ai) == 0 && is_residual_edge(*ai))
Chris@16:                 push_flow(*ai);
Chris@16:               ++ai;
Chris@16:             }
Chris@16:           }
Chris@16:           // do the bottom
Chris@16:           u = bos;
Chris@16:           boost::tie(ai, a_end) = out_edges(u, g);
Chris@16:           while (get(excess_flow, u) > 0 && ai != a_end) {
Chris@16:             if (get(capacity, *ai) == 0 && is_residual_edge(*ai))
Chris@16:               push_flow(*ai);
Chris@16:             ++ai;
Chris@16:           }
Chris@16:         }
Chris@16:         
Chris@16:       } // convert_preflow_to_flow()
Chris@16: 
Chris@16:       //=======================================================================
Chris@16:       inline bool is_flow()
Chris@16:       {
Chris@16:         vertex_iterator u_iter, u_end;
Chris@16:         out_edge_iterator ai, a_end;
Chris@16: 
Chris@16:         // check edge flow values
Chris@16:         for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
Chris@16:           for (boost::tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) {
Chris@16:             edge_descriptor a = *ai;
Chris@16:             if (get(capacity, a) > 0)
Chris@16:               if ((get(residual_capacity, a) + get(residual_capacity, get(reverse_edge, a))
Chris@16:                    != get(capacity, a) + get(capacity, get(reverse_edge, a)))
Chris@16:                   || (get(residual_capacity, a) < 0)
Chris@16:                   || (get(residual_capacity, get(reverse_edge, a)) < 0))
Chris@16:               return false;
Chris@16:           }
Chris@16:         }
Chris@16:         
Chris@16:         // check conservation
Chris@16:         FlowValue sum;  
Chris@16:         for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {
Chris@16:           vertex_descriptor u = *u_iter;
Chris@16:           if (u != src && u != sink) {
Chris@16:             if (get(excess_flow, u) != 0)
Chris@16:               return false;
Chris@16:             sum = 0;
Chris@16:             for (boost::tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) 
Chris@16:               if (get(capacity, *ai) > 0)
Chris@16:                 sum -= get(capacity, *ai) - get(residual_capacity, *ai);
Chris@16:               else
Chris@16:                 sum += get(residual_capacity, *ai);
Chris@16: 
Chris@16:             if (get(excess_flow, u) != sum)
Chris@16:               return false;
Chris@16:           }
Chris@16:         }
Chris@16: 
Chris@16:         return true;
Chris@16:       } // is_flow()
Chris@16: 
Chris@16:       bool is_optimal() {
Chris@16:         // check if mincut is saturated...
Chris@16:         global_distance_update();
Chris@16:         return get(distance, src) >= n;
Chris@16:       }
Chris@16: 
Chris@16:       void print_statistics(std::ostream& os) const {
Chris@16:         os << "pushes:     " << push_count << std::endl
Chris@16:            << "relabels:   " << relabel_count << std::endl
Chris@16:            << "updates:    " << update_count << std::endl
Chris@16:            << "gaps:       " << gap_count << std::endl
Chris@16:            << "gap nodes:  " << gap_node_count << std::endl
Chris@16:            << std::endl;
Chris@16:       }
Chris@16: 
Chris@16:       void print_flow_values(std::ostream& os) const {
Chris@16:         os << "flow values" << std::endl;
Chris@16:         vertex_iterator u_iter, u_end;
Chris@16:         out_edge_iterator ei, e_end;
Chris@16:         for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)
Chris@16:           for (boost::tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei)
Chris@16:             if (get(capacity, *ei) > 0)
Chris@16:               os << *u_iter << " " << target(*ei, g) << " " 
Chris@16:                  << (get(capacity, *ei) - get(residual_capacity, *ei)) << std::endl;
Chris@16:         os << std::endl;
Chris@16:       }
Chris@16: 
Chris@16:       //=======================================================================
Chris@16: 
Chris@16:       Graph& g;
Chris@16:       vertices_size_type n;
Chris@16:       vertices_size_type nm;
Chris@16:       EdgeCapacityMap capacity;
Chris@16:       vertex_descriptor src;
Chris@16:       vertex_descriptor sink;
Chris@16:       VertexIndexMap index;
Chris@16: 
Chris@16:       // will need to use random_access_property_map with these
Chris@16:       std::vector< FlowValue > excess_flow_data;
Chris@16:       iterator_property_map<typename std::vector<FlowValue>::iterator, VertexIndexMap> excess_flow;
Chris@16:       std::vector< std::pair<out_edge_iterator, out_edge_iterator> > current_data;
Chris@16:       iterator_property_map<
Chris@16:         typename std::vector< std::pair<out_edge_iterator, out_edge_iterator> >::iterator,
Chris@16:         VertexIndexMap> current;
Chris@16:       std::vector< distance_size_type > distance_data;
Chris@16:       iterator_property_map<
Chris@16:         typename std::vector< distance_size_type >::iterator,
Chris@16:         VertexIndexMap> distance;
Chris@16:       std::vector< default_color_type > color_data;
Chris@16:       iterator_property_map<
Chris@16:         std::vector< default_color_type >::iterator,
Chris@16:         VertexIndexMap> color;
Chris@16: 
Chris@16:       // Edge Property Maps that must be interior to the graph
Chris@16:       ReverseEdgeMap reverse_edge;
Chris@16:       ResidualCapacityEdgeMap residual_capacity;
Chris@16: 
Chris@16:       LayerArray layers;
Chris@16:       std::vector< list_iterator > layer_list_ptr_data;
Chris@16:       iterator_property_map<typename std::vector< list_iterator >::iterator, VertexIndexMap> layer_list_ptr;
Chris@16:       distance_size_type max_distance;  // maximal distance
Chris@16:       distance_size_type max_active;    // maximal distance with active node
Chris@16:       distance_size_type min_active;    // minimal distance with active node
Chris@16:       boost::queue<vertex_descriptor> Q;
Chris@16: 
Chris@16:       // Statistics counters
Chris@16:       long push_count;
Chris@16:       long update_count;
Chris@16:       long relabel_count;
Chris@16:       long gap_count;
Chris@16:       long gap_node_count;
Chris@16: 
Chris@16:       inline double global_update_frequency() { return 0.5; }
Chris@16:       inline vertices_size_type alpha() { return 6; }
Chris@16:       inline long beta() { return 12; }
Chris@16: 
Chris@16:       long work_since_last_update;
Chris@16:     };
Chris@16: 
Chris@16:   } // namespace detail
Chris@16:   
Chris@16:   template <class Graph, 
Chris@16:             class CapacityEdgeMap, class ResidualCapacityEdgeMap,
Chris@16:             class ReverseEdgeMap, class VertexIndexMap>
Chris@16:   typename property_traits<CapacityEdgeMap>::value_type
Chris@16:   push_relabel_max_flow
Chris@16:     (Graph& g, 
Chris@16:      typename graph_traits<Graph>::vertex_descriptor src,
Chris@16:      typename graph_traits<Graph>::vertex_descriptor sink,
Chris@16:      CapacityEdgeMap cap, ResidualCapacityEdgeMap res,
Chris@16:      ReverseEdgeMap rev, VertexIndexMap index_map)
Chris@16:   {
Chris@16:     typedef typename property_traits<CapacityEdgeMap>::value_type FlowValue;
Chris@16:     
Chris@16:     detail::push_relabel<Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, 
Chris@16:       ReverseEdgeMap, VertexIndexMap, FlowValue>
Chris@16:       algo(g, cap, res, rev, src, sink, index_map);
Chris@16:     
Chris@16:     FlowValue flow = algo.maximum_preflow();
Chris@16:     
Chris@16:     algo.convert_preflow_to_flow();
Chris@16:     
Chris@16:     BOOST_ASSERT(algo.is_flow());
Chris@16:     BOOST_ASSERT(algo.is_optimal());
Chris@16:     
Chris@16:     return flow;
Chris@16:   } // push_relabel_max_flow()
Chris@16:   
Chris@16:   template <class Graph, class P, class T, class R>
Chris@16:   typename detail::edge_capacity_value<Graph, P, T, R>::type
Chris@16:   push_relabel_max_flow
Chris@16:     (Graph& g, 
Chris@16:      typename graph_traits<Graph>::vertex_descriptor src,
Chris@16:      typename graph_traits<Graph>::vertex_descriptor sink,
Chris@16:      const bgl_named_params<P, T, R>& params)
Chris@16:   {
Chris@16:     return push_relabel_max_flow
Chris@16:       (g, src, sink,
Chris@16:        choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity),
Chris@16:        choose_pmap(get_param(params, edge_residual_capacity), 
Chris@16:                    g, edge_residual_capacity),
Chris@16:        choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse),
Chris@16:        choose_const_pmap(get_param(params, vertex_index), g, vertex_index)
Chris@16:        );
Chris@16:   }
Chris@16: 
Chris@16:   template <class Graph>
Chris@16:   typename property_traits<
Chris@16:     typename property_map<Graph, edge_capacity_t>::const_type
Chris@16:   >::value_type
Chris@16:   push_relabel_max_flow
Chris@16:     (Graph& g, 
Chris@16:      typename graph_traits<Graph>::vertex_descriptor src,
Chris@16:      typename graph_traits<Graph>::vertex_descriptor sink)
Chris@16:   {
Chris@16:     bgl_named_params<int, buffer_param_t> params(0); // bogus empty param
Chris@16:     return push_relabel_max_flow(g, src, sink, params);
Chris@16:   }
Chris@16: 
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_PUSH_RELABEL_MAX_FLOW_HPP
Chris@16: