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