source: NonGTP/Boost/boost/graph/push_relabel_max_flow.hpp @ 857

//=======================================================================
// 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