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