// -*- 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 marci {

  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 marci

#endif