//=======================================================================

 // Copyright 2000 University of Notre Dame.

 // Authors: Jeremy G. Siek, Andrew Lumsdaine, Lie-Quan Lee

 //

 // Distributed under the Boost Software License, Version 1.0. (See

 // accompanying file LICENSE_1_0.txt or copy at

 // http://www.boost.org/LICENSE_1_0.txt)

 //=======================================================================

 #ifndef BOOST_PUSH_RELABEL_MAX_FLOW_HPP

 #define BOOST_PUSH_RELABEL_MAX_FLOW_HPP

 #include <boost/config.hpp>

 #include <cassert>

 #include <vector>

 #include <list>

 #include <iosfwd>

 #include <algorithm> // for std::min and std::max

 #include <boost/pending/queue.hpp>

 #include <boost/limits.hpp>

 #include <boost/graph/graph_concepts.hpp>

 #include <boost/graph/named_function_params.hpp>

 namespace boost {

   namespace detail {

    // This implementation is based on Goldberg's 

    // "On Implementing Push-Relabel Method for the Maximum Flow Problem"

    // by B.V. Cherkassky and A.V. Goldberg, IPCO '95, pp. 157--171

    // and on the h_prf.c and hi_pr.c code written by the above authors.

    // This implements the highest-label version of the push-relabel method

    // with the global relabeling and gap relabeling heuristics.

    // The terms "rank", "distance", "height" are synonyms in

    // Goldberg's implementation, paper and in the CLR.  A "layer" is a

    // group of vertices with the same distance. The vertices in each

    // layer are categorized as active or inactive.  An active vertex

    // has positive excess flow and its distance is less than n (it is

    // not blocked).

     template <class Vertex>

     struct preflow_layer {

       std::list<Vertex> active_vertices;

       std::list<Vertex> inactive_vertices;

     };

     template <class Graph, 

               class EdgeCapacityMap,    // integer value type

               class ResidualCapacityEdgeMap,

               class ReverseEdgeMap,

               class VertexIndexMap,     // vertex_descriptor -> integer

               class FlowValue>

     class push_relabel

     {

     public:

       typedef graph_traits<Graph> Traits;

       typedef typename Traits::vertex_descriptor vertex_descriptor;

       typedef typename Traits::edge_descriptor edge_descriptor;

       typedef typename Traits::vertex_iterator vertex_iterator;

       typedef typename Traits::out_edge_iterator out_edge_iterator;

       typedef typename Traits::vertices_size_type vertices_size_type;

       typedef typename Traits::edges_size_type edges_size_type;

       typedef preflow_layer<vertex_descriptor> Layer;

       typedef std::vector< Layer > LayerArray;

       typedef typename LayerArray::iterator layer_iterator;

       typedef typename LayerArray::size_type distance_size_type;

       typedef color_traits<default_color_type> ColorTraits;

       //=======================================================================

       // Some helper predicates

       inline bool is_admissible(vertex_descriptor u, vertex_descriptor v) {

         return distance[u] == distance[v] + 1;

       }

       inline bool is_residual_edge(edge_descriptor a) {

         return 0 < residual_capacity[a];

       }

       inline bool is_saturated(edge_descriptor a) {

         return residual_capacity[a] == 0;

       }

       //=======================================================================

       // Layer List Management Functions

       typedef typename std::list<vertex_descriptor>::iterator list_iterator;

       void add_to_active_list(vertex_descriptor u, Layer& layer) {

         BOOST_USING_STD_MIN();

         BOOST_USING_STD_MAX();

         layer.active_vertices.push_front(u);

         max_active = max BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], max_active);

         min_active = min BOOST_PREVENT_MACRO_SUBSTITUTION(distance[u], min_active);

         layer_list_ptr[u] = layer.active_vertices.begin();

       }

       void remove_from_active_list(vertex_descriptor u) {

         layers[distance[u]].active_vertices.erase(layer_list_ptr[u]);    

       }

       void add_to_inactive_list(vertex_descriptor u, Layer& layer) {

         layer.inactive_vertices.push_front(u);

         layer_list_ptr[u] = layer.inactive_vertices.begin();

       }

       void remove_from_inactive_list(vertex_descriptor u) {

         layers[distance[u]].inactive_vertices.erase(layer_list_ptr[u]);    

       }

       //=======================================================================

       // initialization

       push_relabel(Graph& g_, 

                    EdgeCapacityMap cap,

                    ResidualCapacityEdgeMap res,

                    ReverseEdgeMap rev,

                    vertex_descriptor src_, 

                    vertex_descriptor sink_,

                    VertexIndexMap idx)

         : g(g_), n(num_vertices(g_)), capacity(cap), src(src_), sink(sink_), 

           index(idx),

           excess_flow(num_vertices(g_)),

           current(num_vertices(g_), out_edges(*vertices(g_).first, g_).second),

           distance(num_vertices(g_)),

           color(num_vertices(g_)),

           reverse_edge(rev),

           residual_capacity(res),

           layers(num_vertices(g_)),

           layer_list_ptr(num_vertices(g_), 

                          layers.front().inactive_vertices.end()),

           push_count(0), update_count(0), relabel_count(0), 

           gap_count(0), gap_node_count(0),

           work_since_last_update(0)

       {

         vertex_iterator u_iter, u_end;

         // Don't count the reverse edges

         edges_size_type m = num_edges(g) / 2;

         nm = alpha() * n + m;

         // Initialize flow to zero which means initializing

         // the residual capacity to equal the capacity.

         out_edge_iterator ei, e_end;

         for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)

           for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) {

             residual_capacity[*ei] = capacity[*ei];

           }

         for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {

           vertex_descriptor u = *u_iter;

           excess_flow[u] = 0;

           current[u] = out_edges(u, g).first;

         }

         bool overflow_detected = false;

         FlowValue test_excess = 0;

         out_edge_iterator a_iter, a_end;

         for (tie(a_iter, a_end) = out_edges(src, g); a_iter != a_end; ++a_iter)

           if (target(*a_iter, g) != src)

             test_excess += residual_capacity[*a_iter];

         if (test_excess > (std::numeric_limits<FlowValue>::max)())

           overflow_detected = true;

         if (overflow_detected)

           excess_flow[src] = (std::numeric_limits<FlowValue>::max)();

         else {

           excess_flow[src] = 0;

           for (tie(a_iter, a_end) = out_edges(src, g); 

                a_iter != a_end; ++a_iter) {

             edge_descriptor a = *a_iter;

             if (target(a, g) != src) {

               ++push_count;

               FlowValue delta = residual_capacity[a];

               residual_capacity[a] -= delta;

               residual_capacity[reverse_edge[a]] += delta;

               excess_flow[target(a, g)] += delta;

             }

           }

         }

         max_distance = num_vertices(g) - 1;

         max_active = 0;

         min_active = n;

         for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {

           vertex_descriptor u = *u_iter;

           if (u == sink) {

             distance[u] = 0;

             continue;

           } else if (u == src && !overflow_detected)

             distance[u] = n;

           else

             distance[u] = 1;

           if (excess_flow[u] > 0)

             add_to_active_list(u, layers[1]);

           else if (distance[u] < n)

             add_to_inactive_list(u, layers[1]);

         }       

       } // push_relabel constructor

       //=======================================================================

       // This is a breadth-first search over the residual graph

       // (well, actually the reverse of the residual graph).

       // Would be cool to have a graph view adaptor for hiding certain

       // edges, like the saturated (non-residual) edges in this case.

       // Goldberg's implementation abused "distance" for the coloring.

       void global_distance_update()

       {

         BOOST_USING_STD_MAX();

         ++update_count;

         vertex_iterator u_iter, u_end;

         for (tie(u_iter,u_end) = vertices(g); u_iter != u_end; ++u_iter) {

           color[*u_iter] = ColorTraits::white();

           distance[*u_iter] = n;

         }

         color[sink] = ColorTraits::gray();

         distance[sink] = 0;

         for (distance_size_type l = 0; l <= max_distance; ++l) {

           layers[l].active_vertices.clear();

           layers[l].inactive_vertices.clear();

         }

         max_distance = max_active = 0;

         min_active = n;

         Q.push(sink);

         while (! Q.empty()) {

           vertex_descriptor u = Q.top();

           Q.pop();

           distance_size_type d_v = distance[u] + 1;

           out_edge_iterator ai, a_end;

           for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) {

             edge_descriptor a = *ai;

             vertex_descriptor v = target(a, g);

             if (color[v] == ColorTraits::white()

                 && is_residual_edge(reverse_edge[a])) {

               distance[v] = d_v;

               color[v] = ColorTraits::gray();

               current[v] = out_edges(v, g).first;

               max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(d_v, max_distance);

               if (excess_flow[v] > 0)

                 add_to_active_list(v, layers[d_v]);

               else

                 add_to_inactive_list(v, layers[d_v]);

               Q.push(v);

             }

           }

         }

       } // global_distance_update()

       //=======================================================================

       // This function is called "push" in Goldberg's h_prf implementation,

       // but it is called "discharge" in the paper and in hi_pr.c.

       void discharge(vertex_descriptor u)

       {

         assert(excess_flow[u] > 0);

         while (1) {

           out_edge_iterator ai, ai_end;

           for (ai = current[u], ai_end = out_edges(u, g).second;

                ai != ai_end; ++ai) {

             edge_descriptor a = *ai;

             if (is_residual_edge(a)) {

               vertex_descriptor v = target(a, g);

               if (is_admissible(u, v)) {

                 ++push_count;

                 if (v != sink && excess_flow[v] == 0) {

                   remove_from_inactive_list(v);

                   add_to_active_list(v, layers[distance[v]]);

                 }

                 push_flow(a);

                 if (excess_flow[u] == 0)

                   break;

               } 

             } 

           } // for out_edges of i starting from current

           Layer& layer = layers[distance[u]];

           distance_size_type du = distance[u];

           if (ai == ai_end) {   // i must be relabeled

             relabel_distance(u);

             if (layer.active_vertices.empty()

                 && layer.inactive_vertices.empty())

               gap(du);

             if (distance[u] == n)

               break;

           } else {              // i is no longer active

             current[u] = ai;

             add_to_inactive_list(u, layer);

             break;

           }

         } // while (1)

       } // discharge()

       //=======================================================================

       // This corresponds to the "push" update operation of the paper,

       // not the "push" function in Goldberg's h_prf.c implementation.

       // The idea is to push the excess flow from from vertex u to v.

       void push_flow(edge_descriptor u_v)

       {

         vertex_descriptor

           u = source(u_v, g),

           v = target(u_v, g);

         BOOST_USING_STD_MIN();

         FlowValue flow_delta

           = min BOOST_PREVENT_MACRO_SUBSTITUTION(excess_flow[u], residual_capacity[u_v]);

         residual_capacity[u_v] -= flow_delta;

         residual_capacity[reverse_edge[u_v]] += flow_delta;

         excess_flow[u] -= flow_delta;

         excess_flow[v] += flow_delta;

       } // push_flow()

       //=======================================================================

       // The main purpose of this routine is to set distance[v]

       // to the smallest value allowed by the valid labeling constraints,

       // which are:

       // distance[t] = 0

       // distance[u] <= distance[v] + 1   for every residual edge (u,v)

       //

       distance_size_type relabel_distance(vertex_descriptor u)

       {

         BOOST_USING_STD_MAX();

         ++relabel_count;

         work_since_last_update += beta();

         distance_size_type min_distance = num_vertices(g);

         distance[u] = min_distance;

         // Examine the residual out-edges of vertex i, choosing the

         // edge whose target vertex has the minimal distance.

         out_edge_iterator ai, a_end, min_edge_iter;

         for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) {

           ++work_since_last_update;

           edge_descriptor a = *ai;

           vertex_descriptor v = target(a, g);

           if (is_residual_edge(a) && distance[v] < min_distance) {

             min_distance = distance[v];

             min_edge_iter = ai;

           }

         }

         ++min_distance;

         if (min_distance < n) {

           distance[u] = min_distance;     // this is the main action

           current[u] = min_edge_iter;

           max_distance = max BOOST_PREVENT_MACRO_SUBSTITUTION(min_distance, max_distance);

         }

         return min_distance;

       } // relabel_distance()

       //=======================================================================

       // cleanup beyond the gap

       void gap(distance_size_type empty_distance)

       {

         ++gap_count;

         distance_size_type r; // distance of layer before the current layer

         r = empty_distance - 1;

         // Set the distance for the vertices beyond the gap to "infinity".

         for (layer_iterator l = layers.begin() + empty_distance + 1;

              l < layers.begin() + max_distance; ++l) {

           list_iterator i;

           for (i = l->inactive_vertices.begin(); 

                i != l->inactive_vertices.end(); ++i) {

             distance[*i] = n;

             ++gap_node_count;

           }

           l->inactive_vertices.clear();

         }

         max_distance = r;

         max_active = r;

       }

       //=======================================================================

       // This is the core part of the algorithm, "phase one".

       FlowValue maximum_preflow()

       {

         work_since_last_update = 0;

         while (max_active >= min_active) { // "main" loop

           Layer& layer = layers[max_active];

           list_iterator u_iter = layer.active_vertices.begin();

           if (u_iter == layer.active_vertices.end())

             --max_active;

           else {

             vertex_descriptor u = *u_iter;

             remove_from_active_list(u);

             discharge(u);

             if (work_since_last_update * global_update_frequency() > nm) {

               global_distance_update();

               work_since_last_update = 0;

             }

           }

         } // while (max_active >= min_active)

         return excess_flow[sink];

       } // maximum_preflow()

       //=======================================================================

       // remove excess flow, the "second phase"

       // This does a DFS on the reverse flow graph of nodes with excess flow.

       // If a cycle is found, cancel it.

       // Return the nodes with excess flow in topological order.

       //

       // Unlike the prefl_to_flow() implementation, we use

       //   "color" instead of "distance" for the DFS labels

       //   "parent" instead of nl_prev for the DFS tree

       //   "topo_next" instead of nl_next for the topological ordering

       void convert_preflow_to_flow()

       {

         vertex_iterator u_iter, u_end;

         out_edge_iterator ai, a_end;

         vertex_descriptor r, restart, u;

         std::vector<vertex_descriptor> parent(n);

         std::vector<vertex_descriptor> topo_next(n);

         vertex_descriptor tos(parent[0]), 

           bos(parent[0]); // bogus initialization, just to avoid warning

         bool bos_null = true;

         // handle self-loops

         for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)

           for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai)

             if (target(*ai, g) == *u_iter)

               residual_capacity[*ai] = capacity[*ai];

         // initialize

         for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {

           u = *u_iter;

           color[u] = ColorTraits::white();

           parent[u] = u;

           current[u] = out_edges(u, g).first;

         }

         // eliminate flow cycles and topologically order the vertices

         for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {

           u = *u_iter;

           if (color[u] == ColorTraits::white() 

               && excess_flow[u] > 0

               && u != src && u != sink ) {

             r = u;

             color[r] = ColorTraits::gray();

             while (1) {

               for (; current[u] != out_edges(u, g).second; ++current[u]) {

                 edge_descriptor a = *current[u];

                 if (capacity[a] == 0 && is_residual_edge(a)) {

                   vertex_descriptor v = target(a, g);

                   if (color[v] == ColorTraits::white()) {

                     color[v] = ColorTraits::gray();

                     parent[v] = u;

                     u = v;

                     break;

                   } else if (color[v] == ColorTraits::gray()) {

                     // find minimum flow on the cycle

                     FlowValue delta = residual_capacity[a];

                     while (1) {

                       BOOST_USING_STD_MIN();

                       delta = min BOOST_PREVENT_MACRO_SUBSTITUTION(delta, residual_capacity[*current[v]]);

                       if (v == u)

                         break;

                       else

                         v = target(*current[v], g);

                     }

                     // remove delta flow units

                     v = u;

                     while (1) {

                       a = *current[v];

                       residual_capacity[a] -= delta;

                       residual_capacity[reverse_edge[a]] += delta;

                       v = target(a, g);

                       if (v == u)

                         break;

                     }

                     // back-out of DFS to the first saturated edge

                     restart = u;

                     for (v = target(*current[u], g); v != u; v = target(a, g)){

                       a = *current[v];

                       if (color[v] == ColorTraits::white() 

                           || is_saturated(a)) {

                         color[target(*current[v], g)] = ColorTraits::white();

                         if (color[v] != ColorTraits::white())

                           restart = v;

                       }

                     }

                     if (restart != u) {

                       u = restart;

                       ++current[u];

                       break;

                     }

                   } // else if (color[v] == ColorTraits::gray())

                 } // if (capacity[a] == 0 ...

               } // for out_edges(u, g)  (though "u" changes during loop)

               if (current[u] == out_edges(u, g).second) {

                 // scan of i is complete

                 color[u] = ColorTraits::black();

                 if (u != src) {

                   if (bos_null) {

                     bos = u;

                     bos_null = false;

                     tos = u;

                   } else {

                     topo_next[u] = tos;

                     tos = u;

                   }

                 }

                 if (u != r) {

                   u = parent[u];

                   ++current[u];

                 } else

                   break;

               }

             } // while (1)

           } // if (color[u] == white && excess_flow[u] > 0 & ...)

         } // for all vertices in g

         // return excess flows

         // note that the sink is not on the stack

         if (! bos_null) {

           for (u = tos; u != bos; u = topo_next[u]) {

             ai = out_edges(u, g).first;

             while (excess_flow[u] > 0 && ai != out_edges(u, g).second) {

               if (capacity[*ai] == 0 && is_residual_edge(*ai))

                 push_flow(*ai);

               ++ai;

             }

           }

           // do the bottom

           u = bos;

           ai = out_edges(u, g).first;

           while (excess_flow[u] > 0) {

             if (capacity[*ai] == 0 && is_residual_edge(*ai))

               push_flow(*ai);

             ++ai;

           }

         }

       } // convert_preflow_to_flow()

       //=======================================================================

       inline bool is_flow()

       {

         vertex_iterator u_iter, u_end;

         out_edge_iterator ai, a_end;

         // check edge flow values

         for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {

           for (tie(ai, a_end) = out_edges(*u_iter, g); ai != a_end; ++ai) {

             edge_descriptor a = *ai;

             if (capacity[a] > 0)

               if ((residual_capacity[a] + residual_capacity[reverse_edge[a]]

                    != capacity[a] + capacity[reverse_edge[a]])

                   || (residual_capacity[a] < 0)

                   || (residual_capacity[reverse_edge[a]] < 0))

               return false;

           }

         }

         // check conservation

         FlowValue sum;  

         for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) {

           vertex_descriptor u = *u_iter;

           if (u != src && u != sink) {

             if (excess_flow[u] != 0)

               return false;

             sum = 0;

             for (tie(ai, a_end) = out_edges(u, g); ai != a_end; ++ai) 

               if (capacity[*ai] > 0)

                 sum -= capacity[*ai] - residual_capacity[*ai];

               else

                 sum += residual_capacity[*ai];

             if (excess_flow[u] != sum)

               return false;

           }

         }

         return true;

       } // is_flow()

       bool is_optimal() {

         // check if mincut is saturated...

         global_distance_update();

         return distance[src] >= n;

       }

       void print_statistics(std::ostream& os) const {

         os << "pushes:     " << push_count << std::endl

            << "relabels:   " << relabel_count << std::endl

            << "updates:    " << update_count << std::endl

            << "gaps:       " << gap_count << std::endl

            << "gap nodes:  " << gap_node_count << std::endl

            << std::endl;

       }

       void print_flow_values(std::ostream& os) const {

         os << "flow values" << std::endl;

         vertex_iterator u_iter, u_end;

         out_edge_iterator ei, e_end;

         for (tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter)

           for (tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei)

             if (capacity[*ei] > 0)

               os << *u_iter << " " << target(*ei, g) << " " 

                  << (capacity[*ei] - residual_capacity[*ei]) << std::endl;

         os << std::endl;

       }

       //=======================================================================

       Graph& g;

       vertices_size_type n;

       vertices_size_type nm;

       EdgeCapacityMap capacity;

       vertex_descriptor src;

       vertex_descriptor sink;

       VertexIndexMap index;

       // will need to use random_access_property_map with these

       std::vector< FlowValue > excess_flow;

       std::vector< out_edge_iterator > current;

       std::vector< distance_size_type > distance;

       std::vector< default_color_type > color;

       // Edge Property Maps that must be interior to the graph

       ReverseEdgeMap reverse_edge;

       ResidualCapacityEdgeMap residual_capacity;

       LayerArray layers;

       std::vector< list_iterator > layer_list_ptr;

       distance_size_type max_distance;  // maximal distance

       distance_size_type max_active;    // maximal distance with active node

       distance_size_type min_active;    // minimal distance with active node

       boost::queue<vertex_descriptor> Q;

       // Statistics counters

       long push_count;

       long update_count;

       long relabel_count;

       long gap_count;

       long gap_node_count;

       inline double global_update_frequency() { return 0.5; }

       inline vertices_size_type alpha() { return 6; }

       inline long beta() { return 12; }

       long work_since_last_update;

     };

   } // namespace detail

   template <class Graph, 

             class CapacityEdgeMap, class ResidualCapacityEdgeMap,

             class ReverseEdgeMap, class VertexIndexMap>

   typename property_traits<CapacityEdgeMap>::value_type

   push_relabel_max_flow

     (Graph& g, 

      typename graph_traits<Graph>::vertex_descriptor src,

      typename graph_traits<Graph>::vertex_descriptor sink,

      CapacityEdgeMap cap, ResidualCapacityEdgeMap res,

      ReverseEdgeMap rev, VertexIndexMap index_map)

   {

     typedef typename property_traits<CapacityEdgeMap>::value_type FlowValue;

     detail::push_relabel<Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, 

       ReverseEdgeMap, VertexIndexMap, FlowValue>

       algo(g, cap, res, rev, src, sink, index_map);

     FlowValue flow = algo.maximum_preflow();

     algo.convert_preflow_to_flow();

     assert(algo.is_flow());

     assert(algo.is_optimal());

     return flow;

   } // push_relabel_max_flow()

   template <class Graph, class P, class T, class R>

   typename detail::edge_capacity_value<Graph, P, T, R>::type

   push_relabel_max_flow

     (Graph& g, 

      typename graph_traits<Graph>::vertex_descriptor src,

      typename graph_traits<Graph>::vertex_descriptor sink,

      const bgl_named_params<P, T, R>& params)

   {

     return push_relabel_max_flow

       (g, src, sink,

        choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity),

        choose_pmap(get_param(params, edge_residual_capacity), 

                    g, edge_residual_capacity),

        choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse),

        choose_const_pmap(get_param(params, vertex_index), g, vertex_index)

        );

   }

   template <class Graph>

   typename property_traits<

     typename property_map<Graph, edge_capacity_t>::const_type

   >::value_type

   push_relabel_max_flow

     (Graph& g, 

      typename graph_traits<Graph>::vertex_descriptor src,

      typename graph_traits<Graph>::vertex_descriptor sink)

   {

     bgl_named_params<int, buffer_param_t> params(0); // bogus empty param

     return push_relabel_max_flow(g, src, sink, params);

   }

 } // namespace boost

 #endif // BOOST_PUSH_RELABEL_MAX_FLOW_HPP