// -*- 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 <bfs_mm.h>

 #include <lemon/invalid.h>

 #include <lemon/maps.h>

 #include <demo/tight_edge_filter_map.h>

 /// \file

 /// \brief Maximum flow algorithms.

 /// \ingroup galgs

 namespace lemon {

   using lemon::marci::BfsIterator;

   using lemon::marci::DfsIterator;

   /// \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:

       typedef Node Key;

       typedef bool Value;

       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->target(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->source(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.source(e);

 	Node w=res_graph.target(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.target(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.source(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.source(e);

 	Node w=res_graph.target(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.target(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.source(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);

     {

       for(typename ResGW::NodeIt n(res_graph); 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.target(e), dist[res_graph.source(e)]+1);

 	  typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.source(e)],

 					res_graph_to_F[res_graph.target(e)]);

 	  //original_edge.update();

 	  original_edge.set(f, e);

 	  //residual_capacity.update();

 	  residual_capacity.set(f, res_cap[e]);

 	} else {

 	  if (dist[res_graph.target(e)]==(dist[res_graph.source(e)]+1)) {

 	    typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.source(e)],

 					  res_graph_to_F[res_graph.target(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.source(dfs);

 	    typename MG::Node w=F.target(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.source(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.target(e), potential[res_graph.source(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.source(dfs);

 	    typename ErasingResGW::Node w=erasing_res_graph.target(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.source(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