source: lemon-0.x/src/work/marci/augmenting_flow.h @ 862:732f2acb7239

// -*- C++ -*-
#ifndef HUGO_AUGMENTING_FLOW_H
#define HUGO_AUGMENTING_FLOW_H

#include <vector>
//#include <queue>
//#include <stack>
#include <iostream>

#include <hugo/graph_wrapper.h>
#include <bfs_dfs.h>
#include <hugo/invalid.h>
#include <hugo/maps.h>

/// \file
/// \brief Maximum flow algorithms.
/// \ingroup galgs

namespace hugo {

  /// \brief A map for filtering the edge-set to those edges 
  /// which are tight w.r.t. some node_potential map and 
  /// edge_distance map.
  ///
  /// A node-map node_potential is said to be a potential w.r.t. 
  /// an edge-map edge_distance 
  /// if and only if for each edge e, node_potential[g.head(e)] 
  /// <= edge_distance[e]+node_potential[g.tail(e)] 
  /// (or the reverse inequality holds for each edge).
  /// An edge is said to be tight if this inequality holds with equality, 
  /// and the map returns true exactly for those edges.
  /// To avoid rounding errors, it is recommended to use this class with exact 
  /// types, e.g. with int.
  template<typename Graph, 
           typename NodePotentialMap, typename EdgeDistanceMap>
  class TightEdgeFilterMap {
  protected:
    const Graph* g;
    NodePotentialMap* node_potential;
    EdgeDistanceMap* edge_distance;
  public:
    TightEdgeFilterMap(Graph& _g, NodePotentialMap& _node_potential, 
                       EdgeDistanceMap& _edge_distance) : 
      g(&_g), node_potential(&_node_potential), 
      edge_distance(&_edge_distance) { }
//     void set(const typename Graph::Node& n, int a) {
//       pot->set(n, a);
//     }
//     int operator[](const typename Graph::Node& n) const { 
//       return (*node_potential)[n]; 
//     }
    bool operator[](const typename Graph::Edge& e) const {
      return ((*node_potential)[g->head(e)] == 
              (*edge_distance)[e]+(*node_potential)[g->tail(e)]);
    }
  };

  /// \addtogroup galgs
  /// @{                                                                                                                                        
  /// Class for augmenting path flow algorithms.

  /// This class provides various algorithms for finding a flow of
  /// maximum value in a directed graph. The \e source node, the \e
  /// target node, the \e capacity of the edges and the \e starting \e
  /// flow value of the edges should be passed to the algorithm through the
  /// constructor. 
//   /// It is possible to change these quantities using the
//   /// functions \ref resetSource, \ref resetTarget, \ref resetCap and
//   /// \ref resetFlow. Before any subsequent runs of any algorithm of
//   /// the class \ref resetFlow should be called. 

  /// After running an algorithm of the class, the actual flow value 
  /// can be obtained by calling \ref flowValue(). The minimum
  /// value cut can be written into a \c node map of \c bools by
  /// calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
  /// the inclusionwise minimum and maximum of the minimum value
  /// cuts, resp.)                                                                                                                               
  ///\param Graph The directed graph type the algorithm runs on.
  ///\param Num The number type of the capacities and the flow values.
  ///\param CapMap The capacity map type.
  ///\param FlowMap The flow map type.                                                                                                           
  ///\author Marton Makai
  template <typename Graph, typename Num,
            typename CapMap=typename Graph::template EdgeMap<Num>,
            typename FlowMap=typename Graph::template EdgeMap<Num> >
  class AugmentingFlow {
  protected:
    typedef typename Graph::Node Node;
    typedef typename Graph::NodeIt NodeIt;
    typedef typename Graph::EdgeIt EdgeIt;
    typedef typename Graph::OutEdgeIt OutEdgeIt;
    typedef typename Graph::InEdgeIt InEdgeIt;

    const Graph* g;
    Node s;
    Node t;
    const CapMap* capacity;
    FlowMap* flow;
//    int n;      //the number of nodes of G
    typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
    //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    typedef typename ResGW::Edge ResGWEdge;
    //typedef typename ResGW::template NodeMap<bool> ReachedMap;
    typedef typename Graph::template NodeMap<int> ReachedMap;

    //level works as a bool map in augmenting path algorithms and is
    //used by bfs for storing reached information.  In preflow, it
    //shows the levels of nodes.     
    ReachedMap level;

  public:
    ///Indicates the property of the starting flow.

    ///Indicates the property of the starting flow. The meanings are as follows:
    ///- \c ZERO_FLOW: constant zero flow
    ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
    ///the sum of the out-flows in every node except the \e source and
    ///the \e target.
    ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at 
    ///least the sum of the out-flows in every node except the \e source.
    ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be 
    ///set to the constant zero flow in the beginning of the algorithm in this case.
    enum FlowEnum{
      ZERO_FLOW,
      GEN_FLOW,
      PRE_FLOW,
      NO_FLOW
    };

    enum StatusEnum {
      AFTER_NOTHING,
      AFTER_AUGMENTING,
      AFTER_FAST_AUGMENTING, 
      AFTER_PRE_FLOW_PHASE_1,      
      AFTER_PRE_FLOW_PHASE_2
    };

    /// Don not needle this flag only if necessary.
    StatusEnum status;
    int number_of_augmentations;


    template<typename IntMap>
    class TrickyReachedMap {
    protected:
      IntMap* map;
      int* number_of_augmentations;
    public:
      TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) : 
        map(&_map), number_of_augmentations(&_number_of_augmentations) { }
      void set(const Node& n, bool b) {
        if (b)
          map->set(n, *number_of_augmentations);
        else 
          map->set(n, *number_of_augmentations-1);
      }
      bool operator[](const Node& n) const { 
        return (*map)[n]==*number_of_augmentations; 
      }
    };
    
    AugmentingFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
                   FlowMap& _flow) :
      g(&_G), s(_s), t(_t), capacity(&_capacity),
      flow(&_flow), //n(_G.nodeNum()), 
      level(_G), //excess(_G,0), 
      status(AFTER_NOTHING), number_of_augmentations(0) { }

    /// Starting from a flow, this method searches for an augmenting path
    /// according to the Edmonds-Karp algorithm
    /// and augments the flow on if any.
    /// The return value shows if the augmentation was succesful.
    bool augmentOnShortestPath();
    bool augmentOnShortestPath2();

    /// Starting from a flow, this method searches for an augmenting blocking
    /// flow according to Dinits' algorithm and augments the flow on if any.
    /// The blocking flow is computed in a physically constructed
    /// residual graph of type \c Mutablegraph.
    /// The return value show sif the augmentation was succesful.
    template<typename MutableGraph> bool augmentOnBlockingFlow();

    /// The same as \c augmentOnBlockingFlow<MutableGraph> but the
    /// residual graph is not constructed physically.
    /// The return value shows if the augmentation was succesful.
    bool augmentOnBlockingFlow2();

    template<typename _CutMap>
    void actMinCut(_CutMap& M) const {
      NodeIt v;
      switch (status) {
        case AFTER_PRE_FLOW_PHASE_1:
//      std::cout << "AFTER_PRE_FLOW_PHASE_1" << std::endl;
//      for(g->first(v); g->valid(v); g->next(v)) {
//        if (level[v] < n) {
//          M.set(v, false);
//        } else {
//          M.set(v, true);
//        }
//      }
        break;
      case AFTER_PRE_FLOW_PHASE_2:
//      std::cout << "AFTER_PRE_FLOW_PHASE_2" << std::endl;
        break;
      case AFTER_NOTHING:
//      std::cout << "AFTER_NOTHING" << std::endl;
        minMinCut(M);
        break;
      case AFTER_AUGMENTING:
//      std::cout << "AFTER_AUGMENTING" << std::endl;
        for(g->first(v); v!=INVALID; ++v) {
          if (level[v]) {
            M.set(v, true);
          } else {
            M.set(v, false);
          }
        }
        break;
      case AFTER_FAST_AUGMENTING:
//      std::cout << "AFTER_FAST_AUGMENTING" << std::endl;
        for(g->first(v); v!=INVALID; ++v) {
          if (level[v]==number_of_augmentations) {
            M.set(v, true);
          } else {
            M.set(v, false);
          }
        }
        break;
      }
    }

    template<typename _CutMap>
    void minMinCut(_CutMap& M) const {
      std::queue<Node> queue;

      M.set(s,true);
      queue.push(s);

      while (!queue.empty()) {
        Node w=queue.front();
        queue.pop();

        OutEdgeIt e;
        for(g->first(e,w) ; e!=INVALID; ++e) {
          Node v=g->head(e);
          if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
            queue.push(v);
            M.set(v, true);
          }
        }

        InEdgeIt f;
        for(g->first(f,w) ; f!=INVALID; ++f) {
          Node v=g->tail(f);
          if (!M[v] && (*flow)[f] > 0 ) {
            queue.push(v);
            M.set(v, true);
          }
        }
      }
    }

    template<typename _CutMap>
    void minMinCut2(_CutMap& M) const {
      ResGW res_graph(*g, *capacity, *flow);
      BfsIterator<ResGW, _CutMap> bfs(res_graph, M);
      bfs.pushAndSetReached(s);
      while (!bfs.finished()) ++bfs;
    }

    Num flowValue() const {
      Num a=0;
      for (InEdgeIt e(*g, t); e!=INVALID; ++e) a+=(*flow)[e];
      for (OutEdgeIt e(*g, t); e!=INVALID; ++e) a-=(*flow)[e];
      return a;
      //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
    }

  };



  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
  bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
  {
    ResGW res_graph(*g, *capacity, *flow);
    bool _augment=false;

    //ReachedMap level(res_graph);
    for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);
    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
    bfs.pushAndSetReached(s);

    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
    pred.set(s, INVALID);

    typename ResGW::template NodeMap<Num> free(res_graph);

    //searching for augmenting path
    while ( !bfs.finished() ) {
      ResGWEdge e=bfs;
      if (e!=INVALID && bfs.isBNodeNewlyReached()) {
        Node v=res_graph.tail(e);
        Node w=res_graph.head(e);
        pred.set(w, e);
        if (pred[v]!=INVALID) {
          free.set(w, std::min(free[v], res_graph.resCap(e)));
        } else {
          free.set(w, res_graph.resCap(e));
        }
        if (res_graph.head(e)==t) { _augment=true; break; }
      }

      ++bfs;
    } //end of searching augmenting path

    if (_augment) {
      Node n=t;
      Num augment_value=free[t];
      while (pred[n]!=INVALID) {
        ResGWEdge e=pred[n];
        res_graph.augment(e, augment_value);
        n=res_graph.tail(e);
      }
    }

    status=AFTER_AUGMENTING;
    return _augment;
  }

  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
  bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath2()
  {
    ResGW res_graph(*g, *capacity, *flow);
    bool _augment=false;

    if (status!=AFTER_FAST_AUGMENTING) {
      for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0); 
      number_of_augmentations=1;
    } else {
      ++number_of_augmentations;
    }
    TrickyReachedMap<ReachedMap> 
      tricky_reached_map(level, number_of_augmentations);
    //ReachedMap level(res_graph);
//    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
    BfsIterator<ResGW, TrickyReachedMap<ReachedMap> > 
      bfs(res_graph, tricky_reached_map);
    bfs.pushAndSetReached(s);

    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
    pred.set(s, INVALID);

    typename ResGW::template NodeMap<Num> free(res_graph);

    //searching for augmenting path
    while ( !bfs.finished() ) {
      ResGWEdge e=bfs;
      if (e!=INVALID && bfs.isBNodeNewlyReached()) {
        Node v=res_graph.tail(e);
        Node w=res_graph.head(e);
        pred.set(w, e);
        if (pred[v]!=INVALID) {
          free.set(w, std::min(free[v], res_graph.resCap(e)));
        } else {
          free.set(w, res_graph.resCap(e));
        }
        if (res_graph.head(e)==t) { _augment=true; break; }
      }

      ++bfs;
    } //end of searching augmenting path

    if (_augment) {
      Node n=t;
      Num augment_value=free[t];
      while (pred[n]!=INVALID) {
        ResGWEdge e=pred[n];
        res_graph.augment(e, augment_value);
        n=res_graph.tail(e);
      }
    }

    status=AFTER_FAST_AUGMENTING;
    return _augment;
  }


  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
  template<typename MutableGraph>
  bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
  {
    typedef MutableGraph MG;
    bool _augment=false;

    ResGW res_graph(*g, *capacity, *flow);

    //bfs for distances on the residual graph
    //ReachedMap level(res_graph);
    for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);
    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
    bfs.pushAndSetReached(s);
    typename ResGW::template NodeMap<int>
      dist(res_graph); //filled up with 0's

    //F will contain the physical copy of the residual graph
    //with the set of edges which are on shortest paths
    MG F;
    typename ResGW::template NodeMap<typename MG::Node>
      res_graph_to_F(res_graph);
    {
      typename ResGW::NodeIt n;
      for(res_graph.first(n); n!=INVALID; ++n) 
        res_graph_to_F.set(n, F.addNode());
    }

    typename MG::Node sF=res_graph_to_F[s];
    typename MG::Node tF=res_graph_to_F[t];
    typename MG::template EdgeMap<ResGWEdge> original_edge(F);
    typename MG::template EdgeMap<Num> residual_capacity(F);

    while ( !bfs.finished() ) {
      ResGWEdge e=bfs;
      if (e!=INVALID) {
        if (bfs.isBNodeNewlyReached()) {
          dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
          typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
                                        res_graph_to_F[res_graph.head(e)]);
          //original_edge.update();
          original_edge.set(f, e);
          //residual_capacity.update();
          residual_capacity.set(f, res_graph.resCap(e));
        } else {
          if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
            typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
                                          res_graph_to_F[res_graph.head(e)]);
            //original_edge.update();
            original_edge.set(f, e);
            //residual_capacity.update();
            residual_capacity.set(f, res_graph.resCap(e));
          }
        }
      }
      ++bfs;
    } //computing distances from s in the residual graph

    bool __augment=true;

    while (__augment) {
      __augment=false;
      //computing blocking flow with dfs
      DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
      typename MG::template NodeMap<typename MG::Edge> pred(F);
      pred.set(sF, INVALID);
      //invalid iterators for sources

      typename MG::template NodeMap<Num> free(F);

      dfs.pushAndSetReached(sF);
      while (!dfs.finished()) {
        ++dfs;
        if (typename MG::Edge(dfs)!=INVALID) {
          if (dfs.isBNodeNewlyReached()) {
            typename MG::Node v=F.tail(dfs);
            typename MG::Node w=F.head(dfs);
            pred.set(w, dfs);
            if (pred[v]!=INVALID) {
              free.set(w, std::min(free[v], residual_capacity[dfs]));
            } else {
              free.set(w, residual_capacity[dfs]);
            }
            if (w==tF) {
              __augment=true;
              _augment=true;
              break;
            }

          } else {
            F.erase(typename MG::Edge(dfs));
          }
        }
      }

      if (__augment) {
        typename MG::Node n=tF;
        Num augment_value=free[tF];
        while (pred[n]!=INVALID) {
          typename MG::Edge e=pred[n];
          res_graph.augment(original_edge[e], augment_value);
          n=F.tail(e);
          if (residual_capacity[e]==augment_value)
            F.erase(e);
          else
            residual_capacity.set(e, residual_capacity[e]-augment_value);
        }
      }

    }

    status=AFTER_AUGMENTING;
    return _augment;
  }

  /// Blocking flow augmentation without constructing the layered 
  /// graph physically in which the blocking flow is computed.
  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
  bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
  {
    bool _augment=false;

    ResGW res_graph(*g, *capacity, *flow);

    //Potential map, for distances from s
    typename ResGW::template NodeMap<int> potential(res_graph, 0);
    typedef ConstMap<typename ResGW::Edge, int> Const1Map; 
    Const1Map const_1_map(1);
    TightEdgeFilterMap<ResGW, typename ResGW::template NodeMap<int>,
      Const1Map> tight_edge_filter(res_graph, potential, const_1_map);

    for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);
    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
    bfs.pushAndSetReached(s);

    //computing distances from s in the residual graph
    while ( !bfs.finished() ) {
      ResGWEdge e=bfs;
      if (e!=INVALID && bfs.isBNodeNewlyReached())
        potential.set(res_graph.head(e), potential[res_graph.tail(e)]+1);
      ++bfs;
    } 

    //Subgraph containing the edges on some shortest paths 
    //(i.e. tight edges)
    ConstMap<typename ResGW::Node, bool> true_map(true);
    typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
      TightEdgeFilterMap<ResGW, typename ResGW::template NodeMap<int>, 
      Const1Map> > FilterResGW;
    FilterResGW filter_res_graph(res_graph, true_map, tight_edge_filter);

    //Subgraph, which is able to delete edges which are already
    //met by the dfs
    typename FilterResGW::template NodeMap<typename FilterResGW::Edge>
      first_out_edges(filter_res_graph);
    for (typename FilterResGW::NodeIt v(filter_res_graph); v!=INVALID; ++v)
      first_out_edges.set
        (v, typename FilterResGW::OutEdgeIt(filter_res_graph, v));

    typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
      template NodeMap<typename FilterResGW::Edge> > ErasingResGW;
    ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);

    bool __augment=true;

    while (__augment) {

      __augment=false;
      //computing blocking flow with dfs
      DfsIterator< ErasingResGW,
        typename ErasingResGW::template NodeMap<bool> >
        dfs(erasing_res_graph);
      typename ErasingResGW::
        template NodeMap<typename ErasingResGW::Edge> pred(erasing_res_graph);
      pred.set(s, INVALID);
      //invalid iterators for sources

      typename ErasingResGW::template NodeMap<Num>
        free1(erasing_res_graph);

      dfs.pushAndSetReached
        /// \bug hugo 0.2
        (typename ErasingResGW::Node
         (typename FilterResGW::Node
          (typename ResGW::Node(s)
           )
          )
         );
        
      while (!dfs.finished()) {
        ++dfs;
        if (typename ErasingResGW::Edge(dfs)!=INVALID) {
          if (dfs.isBNodeNewlyReached()) {
            
            typename ErasingResGW::Node v=erasing_res_graph.tail(dfs);
            typename ErasingResGW::Node w=erasing_res_graph.head(dfs);

            pred.set(w, typename ErasingResGW::Edge(dfs));
            if (pred[v]!=INVALID) {
              free1.set
                (w, std::min(free1[v], res_graph.resCap
                             (typename ErasingResGW::Edge(dfs))));
            } else {
              free1.set
                (w, res_graph.resCap
                 (typename ErasingResGW::Edge(dfs)));
            }

            if (w==t) {
              __augment=true;
              _augment=true;
              break;
            }
          } else {
            erasing_res_graph.erase(dfs);
          }
        }
      }

      if (__augment) {
        typename ErasingResGW::Node
          n=typename FilterResGW::Node(typename ResGW::Node(t));
        Num augment_value=free1[n];
        while (pred[n]!=INVALID) {
          typename ErasingResGW::Edge e=pred[n];
          res_graph.augment(e, augment_value);
          n=erasing_res_graph.tail(e);
          if (res_graph.resCap(e)==0)
            erasing_res_graph.erase(e);
        }
      }

    } //while (__augment)

    status=AFTER_AUGMENTING;
    return _augment;
  }


} //namespace hugo

#endif //HUGO_AUGMENTING_FLOW_H