// -*- C++ -*-
/*
preflow_hl4.h
by jacint. 
Runs the two phase highest label preflow push algorithm. In phase 0
we maintain in a list the nodes in level i < n, and we maintain a 
bound k on the max level i < n containing a node, so we can do
the gap heuristic fast. Phase 1 is the same. (The algorithm is the 
same as preflow.hl3, the only diff is that here we use the gap
heuristic with the list of the nodes on level i, and not a bfs form the
upgraded node.)

In phase 1 we shift everything downwards by n.

Member functions:

void run() : runs the algorithm

 The following functions should be used after run() was already run.

T maxflow() : returns the value of a maximum flow

T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e) 

FlowMap allflow() : returns a maximum flow

void allflow(FlowMap& _flow ) : returns a maximum flow

void mincut(CutMap& M) : sets M to the characteristic vector of a 
     minimum cut. M should be a map of bools initialized to false.

void min_mincut(CutMap& M) : sets M to the characteristic vector of the 
     minimum min cut. M should be a map of bools initialized to false.

void max_mincut(CutMap& M) : sets M to the characteristic vector of the 
     maximum min cut. M should be a map of bools initialized to false.

*/

#ifndef PREFLOW_HL4_H
#define PREFLOW_HL4_H

#include <vector>
#include <stack>
#include <queue>

namespace hugo {

  template <typename Graph, typename T, 
    typename FlowMap=typename Graph::EdgeMap<T>, 
    typename CapMap=typename Graph::EdgeMap<T> >
  class preflow_hl4 {
    
    typedef typename Graph::NodeIt NodeIt;
    typedef typename Graph::EdgeIt EdgeIt;
    typedef typename Graph::EachNodeIt EachNodeIt;
    typedef typename Graph::OutEdgeIt OutEdgeIt;
    typedef typename Graph::InEdgeIt InEdgeIt;
    
    Graph& G;
    NodeIt s;
    NodeIt t;
    FlowMap flow;
    CapMap& capacity;  
    T value;
    
  public:

    preflow_hl4(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
      G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }


    void run() {
 
      bool phase=0;
      int n=G.nodeNum(); 
      int k=n-2;
      int b=k;
      /*
	b is a bound on the highest level of the stack. 
	k is a bound on the highest nonempty level i < n.
      */

      typename Graph::NodeMap<int> level(G,n);      
      typename Graph::NodeMap<T> excess(G); 
      std::vector<std::stack<NodeIt> > stack(n);    
      //Stack of the active nodes in level i < n.
      //We use it in both phases.

      typename Graph::NodeMap<NodeIt> left(G);
      typename Graph::NodeMap<NodeIt> right(G);
      std::vector<NodeIt> level_list(n);
      /*
	Needed for the list of the nodes in level i.
      */

      /*Reverse_bfs from t, to find the starting level.*/
      level.set(t,0);
      std::queue<NodeIt> bfs_queue;
      bfs_queue.push(t);

      while (!bfs_queue.empty()) {

	NodeIt v=bfs_queue.front();	
	bfs_queue.pop();
	int l=level.get(v)+1;

	for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
	  NodeIt w=G.tail(e);
	  if ( level.get(w) == n ) {
	    bfs_queue.push(w);
	    NodeIt first=level_list[l];
	    if ( first != 0 ) left.set(first,w);
	    right.set(w,first);
	    level_list[l]=w;
	    level.set(w, l);
	  }
	}
      }
      
      level.set(s,n);


      /* Starting flow. It is everywhere 0 at the moment. */     
      for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e) 
	{
	  T c=capacity.get(e);
	  if ( c == 0 ) continue;
	  NodeIt w=G.head(e);
	  if ( level.get(w) < n ) {	  
	    if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); 
	    flow.set(e, c); 
	    excess.set(w, excess.get(w)+c);
	  }
	}
      /* 
	 End of preprocessing 
      */


      /*
	Push/relabel on the highest level active nodes.
      */	
      while ( true ) {

	if ( b == 0 ) {
	  if ( phase ) break;
	  
	  /*
	    In the end of phase 0 we apply a bfs from s in
	    the residual graph.
	  */
	  phase=1;
	  level.set(s,0);
	  std::queue<NodeIt> bfs_queue;
	  bfs_queue.push(s);
	  
	  while (!bfs_queue.empty()) {
	    
	    NodeIt v=bfs_queue.front();	
	    bfs_queue.pop();
	    int l=level.get(v)+1;

	    for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
	      if ( capacity.get(e) == flow.get(e) ) continue;
	      NodeIt u=G.tail(e);
	      if ( level.get(u) >= n ) { 
		bfs_queue.push(u);
		level.set(u, l);
		if ( excess.get(u) > 0 ) stack[l].push(u);
	      }
	    }

	    for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
	      if ( 0 == flow.get(e) ) continue;
	      NodeIt u=G.head(e);
	      if ( level.get(u) >= n ) { 
		bfs_queue.push(u);
		level.set(u, l);
		if ( excess.get(u) > 0 ) stack[l].push(u);
	      }
	    }
	  }
	  b=n-2;
	}


	if ( stack[b].empty() ) --b;
	else {
	  
	  NodeIt w=stack[b].top();        //w is a highest label active node.
	  stack[b].pop();           
	  int lev=level.get(w);
	  T exc=excess.get(w);
	  int newlevel=n;          //In newlevel we bound the next level of w.
	  
	  for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
	    
	    if ( flow.get(e) == capacity.get(e) ) continue; 
	    NodeIt v=G.head(e);            
	    //e=wv	    
	    
	    if( lev > level.get(v) ) {      
	      /*Push is allowed now*/
	      
	      if ( excess.get(v)==0 && v!=t && v!=s ) 
		stack[level.get(v)].push(v); 
	      /*v becomes active.*/
	      
	      T cap=capacity.get(e);
	      T flo=flow.get(e);
	      T remcap=cap-flo;
	      
	      if ( remcap >= exc ) {       
		/*A nonsaturating push.*/
		
		flow.set(e, flo+exc);
		excess.set(v, excess.get(v)+exc);
		exc=0;
		break; 
		
	      } else { 
		/*A saturating push.*/
		
		flow.set(e, cap);
		excess.set(v, excess.get(v)+remcap);
		exc-=remcap;
	      }
	    } else if ( newlevel > level.get(v) ){
	      newlevel = level.get(v);
	    }	    
	    
	  } //for out edges wv 
	
	
	if ( exc > 0 ) {	
	  for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
	    
	    if( flow.get(e) == 0 ) continue; 
	    NodeIt v=G.tail(e);  
	    //e=vw
	    
	    if( lev > level.get(v) ) {  
	      /*Push is allowed now*/
	      
	      if ( excess.get(v)==0 && v!=t && v!=s ) 
		stack[level.get(v)].push(v); 
	      /*v becomes active.*/
	      
	      T flo=flow.get(e);
	      
	      if ( flo >= exc ) { 
		/*A nonsaturating push.*/
		
		flow.set(e, flo-exc);
		excess.set(v, excess.get(v)+exc);
		exc=0;
		break; 
	      } else {                                               
		/*A saturating push.*/
		
		excess.set(v, excess.get(v)+flo);
		exc-=flo;
		flow.set(e,0);
	      }  
	    } else if ( newlevel > level.get(v) ) {
	      newlevel = level.get(v);
	    }	    
	  } //for in edges vw
	  
	} // if w still has excess after the out edge for cycle
	
	excess.set(w, exc);
	 
	/*
	  Relabel
	*/
	
	if ( exc > 0 ) {
	  //now 'lev' is the old level of w
	
	  if ( phase ) {
	    level.set(w,++newlevel);
	    stack[newlevel].push(w);
	    b=newlevel;
	  } else {
	    //unlacing
	    NodeIt right_n=right.get(w);
	    NodeIt left_n=left.get(w);

	    if ( right_n != 0 ) {
	      if ( left_n != 0 ) {
		right.set(left_n, right_n);
		left.set(right_n, left_n);
	      } else {
		level_list[lev]=right_n;
		left.set(right_n, 0);
	      } 
	    } else {
	      if ( left_n != 0 ) {
		right.set(left_n, 0);
	      } else { 
		level_list[lev]=0;
	      } 
	    }
		

	    if ( level_list[lev]==0 ) {

	      for (int i=lev; i!=k ; ) {
		NodeIt v=level_list[++i];
		while ( v != 0 ) {
		  level.set(v,n);
		  v=right.get(v);
		}
		level_list[i]=0;
	      }	     

	      level.set(w,n);

	      b=--lev;
	      k=b;

	    } else {

	      if ( newlevel == n ) {
		level.set(w,n);
	      } else {
		
		level.set(w,++newlevel);
		stack[newlevel].push(w);
		b=newlevel;
		if ( k < newlevel ) ++k;
		NodeIt first=level_list[newlevel];
		if ( first != 0 ) left.set(first,w);
		right.set(w,first);
		left.set(w,0);
		level_list[newlevel]=w;
	      }
	    }
	  } //phase 0
	} // if ( exc > 0 )
 
	
	} // if stack[b] is nonempty
	
      } // while(true)


      value = excess.get(t);
      /*Max flow value.*/


    } //void run()





    /*
      Returns the maximum value of a flow.
     */

    T maxflow() {
      return value;
    }



    /*
      For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e). 
    */

    T flowonedge(EdgeIt e) {
      return flow.get(e);
    }



    FlowMap allflow() {
      return flow;
    }


    
    void allflow(FlowMap& _flow ) {
      for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v)
	_flow.set(v,flow.get(v));
    }



    /*
      Returns the minimum min cut, by a bfs from s in the residual graph.
    */
    
    template<typename CutMap>
    void mincut(CutMap& M) {
    
      std::queue<NodeIt> queue;
      
      M.set(s,true);      
      queue.push(s);

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

	for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
	  NodeIt v=G.head(e);
	  if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
	    queue.push(v);
	    M.set(v, true);
	  }
	} 

	for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
	  NodeIt v=G.tail(e);
	  if (!M.get(v) && flow.get(e) > 0 ) {
	    queue.push(v);
	    M.set(v, true);
	  }
	} 

      }

    }



    /*
      Returns the maximum min cut, by a reverse bfs 
      from t in the residual graph.
    */
    
    template<typename CutMap>
    void max_mincut(CutMap& M) {
    
      std::queue<NodeIt> queue;
      
      M.set(t,true);        
      queue.push(t);

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

	for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
	  NodeIt v=G.tail(e);
	  if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
	    queue.push(v);
	    M.set(v, true);
	  }
	}

	for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
	  NodeIt v=G.head(e);
	  if (!M.get(v) && flow.get(e) > 0 ) {
	    queue.push(v);
	    M.set(v, true);
	  }
	}
      }

      for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
	M.set(v, !M.get(v));
      }

    }



    template<typename CutMap>
    void min_mincut(CutMap& M) {
      mincut(M);
    }



  };
}//namespace hugo
#endif