// -*- C++ -*-

#ifndef LEMON_AUGMENTING_FLOW_H

#define LEMON_AUGMENTING_FLOW_H

#include <vector>

#include <iostream>

#include <lemon/graph_wrapper.h>

#include <bfs_dfs.h>

#include <lemon/invalid.h>

#include <lemon/maps.h>

#include <lemon/tight_edge_filter_map.h>

/// \file

/// \brief Maximum flow algorithms.

/// \ingroup galgs

namespace lemon {

  /// \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);

    typename ResGW::ResCap res_cap(res_graph);

    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_cap[e]));

	} else {

	  free.set(w, res_cap[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);

    typename ResGW::ResCap res_cap(res_graph);

    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_cap[e]));

	} else {

	  free.set(w, res_cap[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);

    typename ResGW::ResCap res_cap(res_graph);

    //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_cap[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_cap[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);

    typename ResGW::ResCap res_cap(res_graph);

    //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 lemon 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_cap

			     [typename ErasingResGW::Edge(dfs)]));

	    } else {

	      free1.set

		(w, res_cap

		 [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_cap[e]==0)

	    erasing_res_graph.erase(e);

	}

      }

    } //while (__augment)

    status=AFTER_AUGMENTING;

    return _augment;

  }

} //namespace lemon

#endif //LEMON_AUGMENTING_FLOW_H