// -*- C++ -*-

 /*

 preflow_h5.h

 by jacint. 

 Heuristics: 

  2 phase

  gap

  list 'level_list' on the nodes on level i implemented by hand

  highest label

  relevel: in phase 0, after BFS*n relabels, it runs a reverse 

    bfs from t in the res graph to relevel the nodes reachable from t.

    BFS is initialized to 20

 Due to the last heuristic, this algorithm is quite fast on very 

 sparse graphs, but relatively bad on even the dense graphs.

 'NodeMap<bool> cut' is a member, in this way we can count it fast, after

 the algorithm was run.

 The constructor runs the algorithm.

 Members:

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

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

 FlowMap Flow() : returns the fixed maximum flow x

 void Flow(FlowMap& _flow ) : returns the fixed maximum flow x

 void minMinCut(CutMap& M) : sets M to the characteristic vector of the 

      minimum min cut. M should be a map of bools initialized to false.

 void maxMinCut(CutMap& M) : sets M to the characteristic vector of the 

      maximum min cut. M should be a map of bools initialized to false.

 void minCut(CutMap& M) : fast function, sets M to the characteristic 

      vector of a minimum cut. 

 Different member from the other preflow_hl-s (here we have a member 

 'NodeMap<bool> cut').

 CutMap minCut() : fast function, giving the characteristic 

      vector of a minimum cut.

 */

 #ifndef PREFLOW_HL4_H

 #define PREFLOW_HL4_H

 #define BFS 20

 #include <vector>

 #include <queue>

 #include <time_measure.h>  //for test

 namespace hugo {

   template <typename Graph, typename T, 

     typename FlowMap=typename Graph::EdgeMap<T>, 

     typename CutMap=typename Graph::NodeMap<bool>,

     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;  

     CutMap cut;

     T value;

   public:

     double time;

     preflow_hl4(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :

       G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), 

       cut(G, false) {

       bool phase=0;   //phase 0 is the 1st phase, phase 1 is the 2nd

       int n=G.nodeNum(); 

       int relabel=0;

       int heur=(int)BFS*n;

       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<NodeIt> active(n);

       typename Graph::NodeMap<NodeIt> next(G);

       //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 && w !=s ) {

 	    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 ) {

 	      next.set(w,active[level.get(w)]);

 	      active[level.get(w)]=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;

 	  //Now have a min cut.

 	  for( EachNodeIt v=G.template first<EachNodeIt>(); 

 	       v.valid(); ++v) 

 	    if (level.get(v) >= n ) cut.set(v,true);

 	  time=currTime();

 	  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 ) {

 		    next.set(u,active[l]);

 		    active[l]=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 ) {

 		    next.set(u,active[l]);

 		    active[l]=u;

 		}

 	      }

 	    }

 	  }

 	  b=n-2;

 	}

 	if ( active[b] == 0 ) --b;

 	else {

 	  NodeIt w=active[b];

 	  active[b]=next.get(w);

 	  int lev=level.get(w);

 	  T exc=excess.get(w);

 	  int newlevel=n;          //bound on 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 )  {

 		int lev_v=level.get(v);

 		next.set(v,active[lev_v]);

 		active[lev_v]=v;

 	      }

 	      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 ) {

 		int lev_v=level.get(v);

 		next.set(v,active[lev_v]);

 		active[lev_v]=v;

 	      }

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

 	    next.set(w,active[newlevel]);

 	    active[newlevel]=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;

 	      } 

 	    }

 	    //unlacing ends

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

 		next.set(w,active[newlevel]);

 		active[newlevel]=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;

 	      }

 	    }

 	    ++relabel;

 	    if ( relabel >= heur ) {

 	      relabel=0;

 	      b=n-2;

 	      k=b;

 	      for ( int i=1; i!=n; ++i ) { 

 		active[i]=0;

 		level_list[i]=0;

 	      }

 	      //bfs from t in the res graph to relevel the nodes

 	      for( EachNodeIt v=G.template first<EachNodeIt>(); 

 		   v.valid(); ++v) level.set(v,n);

 	      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) {

 		  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 ) {

 		      next.set(u,active[l]);

 		      active[l]=u;

 		    }

 		    NodeIt first=level_list[l];

 		    if ( first != 0 ) left.set(first,w);

 		    right.set(w,first);

 		    left.set(w,0);

 		    level_list[l]=w;

 		  }

 		}

 		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 ) {

 		      next.set(u,active[l]);

 		      active[l]=u;

 		    }

 		    NodeIt first=level_list[l];

 		    if ( first != 0 ) left.set(first,w);

 		    right.set(w,first);

 		    left.set(w,0);

 		    level_list[l]=w;

 		  }

 		}

 	      }

 	      level.set(s,n);

 	    }

 	  } //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 Flow() {

       return flow;

     }

     void Flow(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 minMinCut(_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 maxMinCut(_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 minCut(_CutMap& M) {

       for( EachNodeIt v=G.template first<EachNodeIt>(); 

 	   v.valid(); ++v) 

 	M.set(v, cut.get(v));

     }

     CutMap minCut() {

       return cut;

     }

   };

 }//namespace marci

 #endif