// -*- C++ -*-

/*

Heuristics: 

 2 phase

 gap

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

 stack 'active' on the active nodes on level i

 runs heuristic 'highest label' for H1*n relabels

 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'

Parameters H0 and H1 are initialized to 20 and 1.

Constructors:

Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if 

     FlowMap is not constant zero, and should be true if it is

Members:

void run()

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

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. ??Is it OK?

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) : sets M to the characteristic vector of 

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

*/

#ifndef HUGO_PREFLOW_H

#define HUGO_PREFLOW_H

#define H0 20

#define H1 1

#include <vector>

#include <queue>

#include <stack>

namespace hugo {

  template <typename Graph, typename T, 

	    typename CapMap=typename Graph::template EdgeMap<T>, 

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

  class Preflow {

    typedef typename Graph::Node Node;

    typedef typename Graph::NodeIt NodeIt;

    typedef typename Graph::OutEdgeIt OutEdgeIt;

    typedef typename Graph::InEdgeIt InEdgeIt;

    typedef typename std::vector<std::stack<Node> > VecStack;

    typedef typename Graph::template NodeMap<Node> NNMap;

    typedef typename std::vector<Node> VecNode;

    const Graph& G;

    Node s;

    Node t;

    CapMap* capacity;  

    FlowMap* flow;

    int n;      //the number of nodes of G

    typename Graph::template NodeMap<int> level;      

    typename Graph::template NodeMap<T> excess; 

  public:

    enum flowEnum{

      ZERO_FLOW=0,

      GEN_FLOW=1,

      PREFLOW=2

    };

    Preflow(Graph& _G, Node _s, Node _t, CapMap& _capacity, 

	    FlowMap& _flow) :

      G(_G), s(_s), t(_t), capacity(&_capacity), 

      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}

    void run() {

      preflow( ZERO_FLOW );

    }

    void preflow( flowEnum fe ) {

      preflowPhase0(fe);

      preflowPhase1();

    }

    void preflowPhase0( flowEnum fe ) {

      int heur0=(int)(H0*n);  //time while running 'bound decrease' 

      int heur1=(int)(H1*n);  //time while running 'highest label'

      int heur=heur1;         //starting time interval (#of relabels)

      int numrelabel=0;

      bool what_heur=1;       

      //It is 0 in case 'bound decrease' and 1 in case 'highest label'

      bool end=false;     

      //Needed for 'bound decrease', true means no active nodes are above bound b.

      int k=n-2;  //bound on the highest level under n containing a node

      int b=k;    //bound on the highest level under n of an active node

      VecStack active(n);

      NNMap left(G,INVALID);

      NNMap right(G,INVALID);

      VecNode level_list(n,INVALID);

      //List of the nodes in level i<n, set to n.

      NodeIt v;

      for(G.first(v); G.valid(v); G.next(v)) level.set(v,n);

      //setting each node to level n

      switch ( fe ) {

      case PREFLOW:

	{

	  //counting the excess

	  NodeIt v;

	  for(G.first(v); G.valid(v); G.next(v)) {

	    T exc=0;

	    InEdgeIt e;

	    for(G.first(e,v); G.valid(e); G.next(e)) exc+=(*flow)[e];

	    OutEdgeIt f;

	    for(G.first(f,v); G.valid(f); G.next(f)) exc-=(*flow)[f];

	    excess.set(v,exc);	  

	    //putting the active nodes into the stack

	    int lev=level[v];

	    if ( exc > 0 && lev < n && v != t ) active[lev].push(v);

	  }

	  break;

	}

      case GEN_FLOW:

	{

	  //Counting the excess of t

	  T exc=0;

	  InEdgeIt e;

	  for(G.first(e,t); G.valid(e); G.next(e)) exc+=(*flow)[e];

	  OutEdgeIt f;

	  for(G.first(f,t); G.valid(f); G.next(f)) exc-=(*flow)[f];

	  excess.set(t,exc);	

	  break;

	}

      }

      preflowPreproc( fe, active, level_list, left, right );

      //End of preprocessing 

      //Push/relabel on the highest level active nodes.

      while ( true ) {

	if ( b == 0 ) {

	  if ( !what_heur && !end && k > 0 ) {

	    b=k;

	    end=true;

	  } else break;

	}

	if ( active[b].empty() ) --b; 

	else {

	  end=false;  

	  Node w=active[b].top();

	  active[b].pop();

	  int newlevel=push(w,active);

	  if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, 

				       left, right, b, k, what_heur);

	  ++numrelabel; 

	  if ( numrelabel >= heur ) {

	    numrelabel=0;

	    if ( what_heur ) {

	      what_heur=0;

	      heur=heur0;

	      end=false;

	    } else {

	      what_heur=1;

	      heur=heur1;

	      b=k; 

	    }

	  }

	} 

      } 

    }

    void preflowPhase1() {

      int k=n-2;  //bound on the highest level under n containing a node

      int b=k;    //bound on the highest level under n of an active node

      VecStack active(n);

      level.set(s,0);

      std::queue<Node> bfs_queue;

      bfs_queue.push(s);

      while (!bfs_queue.empty()) {

	Node v=bfs_queue.front();	

	bfs_queue.pop();

	int l=level[v]+1;

	InEdgeIt e;

	for(G.first(e,v); G.valid(e); G.next(e)) {

	  if ( (*capacity)[e] == (*flow)[e] ) continue;

	  Node u=G.tail(e);

	  if ( level[u] >= n ) { 

	    bfs_queue.push(u);

	    level.set(u, l);

	    if ( excess[u] > 0 ) active[l].push(u);

	  }

	}

	OutEdgeIt f;

	for(G.first(f,v); G.valid(f); G.next(f)) {

	  if ( 0 == (*flow)[f] ) continue;

	  Node u=G.head(f);

	  if ( level[u] >= n ) { 

	    bfs_queue.push(u);

	    level.set(u, l);

	    if ( excess[u] > 0 ) active[l].push(u);

	  }

	}

      }

      b=n-2;

      while ( true ) {

	if ( b == 0 ) break;

	if ( active[b].empty() ) --b; 

	else {

	  Node w=active[b].top();

	  active[b].pop();

	  int newlevel=push(w,active);	  

	  //relabel

	  if ( excess[w] > 0 ) {

	    level.set(w,++newlevel);

	    active[newlevel].push(w);

	    b=newlevel;

	  }

	}  // if stack[b] is nonempty

      } // while(true)

    }

    //Returns the maximum value of a flow.

    T flowValue() {

      return excess[t];

    }

    //should be used only between preflowPhase0 and preflowPhase1

    template<typename _CutMap>

    void actMinCut(_CutMap& M) {

      NodeIt v;

      for(G.first(v); G.valid(v); G.next(v)) 

	if ( level[v] < n ) M.set(v,false);

	else M.set(v,true);

    }

    /*

      Returns the minimum min cut, by a bfs from s in the residual graph.

    */

    template<typename _CutMap>

    void minMinCut(_CutMap& M) {

      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) ; G.valid(e); G.next(e)) {

	  Node v=G.head(e);

	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {

	    queue.push(v);

	    M.set(v, true);

	  }

	} 

	InEdgeIt f;

	for(G.first(f,w) ; G.valid(f); G.next(f)) {

	  Node v=G.tail(f);

	  if (!M[v] && (*flow)[f] > 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) {

      NodeIt v;

      for(G.first(v) ; G.valid(v); G.next(v)) {

	M.set(v, true);

      }

      std::queue<Node> queue;

      M.set(t,false);        

      queue.push(t);

      while (!queue.empty()) {

        Node w=queue.front();

	queue.pop();

	InEdgeIt e;

	for(G.first(e,w) ; G.valid(e); G.next(e)) {

	  Node v=G.tail(e);

	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {

	    queue.push(v);

	    M.set(v, false);

	  }

	}

	OutEdgeIt f;

	for(G.first(f,w) ; G.valid(f); G.next(f)) {

	  Node v=G.head(f);

	  if (M[v] && (*flow)[f] > 0 ) {

	    queue.push(v);

	    M.set(v, false);

	  }

	}

      }

    }

    template<typename CutMap>

    void minCut(CutMap& M) {

      minMinCut(M);

    }

    void resetTarget (const Node _t) {t=_t;}

    void resetSource (const Node _s) {s=_s;}

    void resetCap (const CapMap& _cap) {

      capacity=&_cap;

    }

    void resetFlow (FlowMap& _flow) {

      flow=&_flow;

    }

  private:

    int push(const Node w, VecStack& active) {

      int lev=level[w];

      T exc=excess[w];

      int newlevel=n;       //bound on the next level of w

      OutEdgeIt e;

      for(G.first(e,w); G.valid(e); G.next(e)) {

	if ( (*flow)[e] == (*capacity)[e] ) continue; 

	Node v=G.head(e);            

	if( lev > level[v] ) { //Push is allowed now

	  if ( excess[v]==0 && v!=t && v!=s ) {

	    int lev_v=level[v];

	    active[lev_v].push(v);

	  }

	  T cap=(*capacity)[e];

	  T flo=(*flow)[e];

	  T remcap=cap-flo;

	  if ( remcap >= exc ) { //A nonsaturating push.

	    flow->set(e, flo+exc);

	    excess.set(v, excess[v]+exc);

	    exc=0;

	    break; 

	  } else { //A saturating push.

	    flow->set(e, cap);

	    excess.set(v, excess[v]+remcap);

	    exc-=remcap;

	  }

	} else if ( newlevel > level[v] ) newlevel = level[v];

      } //for out edges wv 

      if ( exc > 0 ) {	

	InEdgeIt e;

	for(G.first(e,w); G.valid(e); G.next(e)) {

	  if( (*flow)[e] == 0 ) continue; 

	  Node v=G.tail(e); 

	  if( lev > level[v] ) { //Push is allowed now

	    if ( excess[v]==0 && v!=t && v!=s ) {

	      int lev_v=level[v];

	      active[lev_v].push(v);

	    }

	    T flo=(*flow)[e];

	    if ( flo >= exc ) { //A nonsaturating push.

	      flow->set(e, flo-exc);

	      excess.set(v, excess[v]+exc);

	      exc=0;

	      break; 

	    } else {  //A saturating push.

	      excess.set(v, excess[v]+flo);

	      exc-=flo;

	      flow->set(e,0);

	    }  

	  } else if ( newlevel > level[v] ) newlevel = level[v];

	} //for in edges vw

      } // if w still has excess after the out edge for cycle

      excess.set(w, exc);

      return newlevel;

     }

    void preflowPreproc ( flowEnum fe, VecStack& active, 

			  VecNode& level_list, NNMap& left, NNMap& right ) {

      std::queue<Node> bfs_queue;

      switch ( fe ) {

      case ZERO_FLOW: 

	{

	  //Reverse_bfs from t, to find the starting level.

	  level.set(t,0);

	  bfs_queue.push(t);

	  while (!bfs_queue.empty()) {

	    Node v=bfs_queue.front();	

	    bfs_queue.pop();

	    int l=level[v]+1;

	    InEdgeIt e;

	    for(G.first(e,v); G.valid(e); G.next(e)) {

	      Node w=G.tail(e);

	      if ( level[w] == n && w != s ) {

		bfs_queue.push(w);

		Node first=level_list[l];

		if ( G.valid(first) ) left.set(first,w);

		right.set(w,first);

		level_list[l]=w;

		level.set(w, l);

	      }

	    }

	  }

	  //the starting flow

	  OutEdgeIt e;

	  for(G.first(e,s); G.valid(e); G.next(e)) 

	    {

	      T c=(*capacity)[e];

	      if ( c == 0 ) continue;

	      Node w=G.head(e);

	      if ( level[w] < n ) {	  

		if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);

		flow->set(e, c); 

		excess.set(w, excess[w]+c);

	      }

	    }

	  break;

	}

      case GEN_FLOW:

      case PREFLOW:

	{

	  //Reverse_bfs from t in the residual graph, 

	  //to find the starting level.

	  level.set(t,0);

	  bfs_queue.push(t);

	  while (!bfs_queue.empty()) {

	    Node v=bfs_queue.front();	

	    bfs_queue.pop();

	    int l=level[v]+1;

	    InEdgeIt e;

	    for(G.first(e,v); G.valid(e); G.next(e)) {

	      if ( (*capacity)[e] == (*flow)[e] ) continue;

	      Node w=G.tail(e);

	      if ( level[w] == n && w != s ) {

		bfs_queue.push(w);

		Node first=level_list[l];

		if ( G.valid(first) ) left.set(first,w);

		right.set(w,first);

		level_list[l]=w;

		level.set(w, l);

	      }

	    }

	    OutEdgeIt f;

	    for(G.first(f,v); G.valid(f); G.next(f)) {

	      if ( 0 == (*flow)[f] ) continue;

	      Node w=G.head(f);

	      if ( level[w] == n && w != s ) {

		bfs_queue.push(w);

		Node first=level_list[l];

		if ( G.valid(first) ) left.set(first,w);

		right.set(w,first);

		level_list[l]=w;

		level.set(w, l);

	      }

	    }

	  }

	  //the starting flow

	  OutEdgeIt e;

	  for(G.first(e,s); G.valid(e); G.next(e)) 

	    {

	      T rem=(*capacity)[e]-(*flow)[e];

	      if ( rem == 0 ) continue;

	      Node w=G.head(e);

	      if ( level[w] < n ) {	  

		if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);

		flow->set(e, (*capacity)[e]); 

		excess.set(w, excess[w]+rem);

	      }

	    }

	  InEdgeIt f;

	  for(G.first(f,s); G.valid(f); G.next(f)) 

	    {

	      if ( (*flow)[f] == 0 ) continue;

	      Node w=G.tail(f);

	      if ( level[w] < n ) {	  

		if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);

		excess.set(w, excess[w]+(*flow)[f]);

		flow->set(f, 0); 

	      }

	    }  

	  break;

	} //case PREFLOW

      }

    } //preflowPreproc

    void relabel( const Node w, int newlevel, VecStack& active,  

		  VecNode& level_list, NNMap& left, 

		  NNMap& right, int& b, int& k, const bool what_heur ) {

      T lev=level[w];	

      Node right_n=right[w];

      Node left_n=left[w];

      //unlacing starts

      if ( G.valid(right_n) ) {

	if ( G.valid(left_n) ) {

	  right.set(left_n, right_n);

	  left.set(right_n, left_n);

	} else {

	  level_list[lev]=right_n;   

	  left.set(right_n, INVALID);

	} 

      } else {

	if ( G.valid(left_n) ) {

	  right.set(left_n, INVALID);

	} else { 

	  level_list[lev]=INVALID;   

	} 

      } 

      //unlacing ends

      if ( !G.valid(level_list[lev]) ) {

	//gapping starts

	for (int i=lev; i!=k ; ) {

	  Node v=level_list[++i];

	  while ( G.valid(v) ) {

	    level.set(v,n);

	    v=right[v];

	  }

	  level_list[i]=INVALID;

	  if ( !what_heur ) {

	    while ( !active[i].empty() ) {

	      active[i].pop();    //FIXME: ezt szebben kene

	    }

	  }	     

	}

	level.set(w,n);

	b=lev-1;

	k=b;

	//gapping ends

      } else {

	if ( newlevel == n ) level.set(w,n); 

	else {

	  level.set(w,++newlevel);

	  active[newlevel].push(w);

	  if ( what_heur ) b=newlevel;

	  if ( k < newlevel ) ++k;      //now k=newlevel

	  Node first=level_list[newlevel];

	  if ( G.valid(first) ) left.set(first,w);

	  right.set(w,first);

	  left.set(w,INVALID);

	  level_list[newlevel]=w;

	}

      }

    } //relabel

  };

} //namespace hugo

#endif //HUGO_PREFLOW_H