#ifndef HUGO_PREFLOW_PUSH_HH

 #define HUGO_PREFLOW_PUSH_HH

 //#include <algorithm>

 #include <list>

 #include <vector>

 #include <queue>

 //#include "pf_hiba.hh"

 //#include <marci_list_graph.hh>

 //#include <marci_graph_traits.hh>

 #include <invalid.h>

 //#include <reverse_bfs.hh>

 using namespace std;

 namespace hugo {

   template <typename Graph, typename T>

   class preflow_push {

     //Useful typedefs

     typedef typename Graph::Node Node;

     typedef typename Graph::NodeIt NodeIt;

     typedef typename Graph::Edge Edge;

     typedef typename Graph::OutEdgeIt OutEdgeIt;

     typedef typename Graph::InEdgeIt InEdgeIt;

     //---------------------------------------------

     //Parameters of the algorithm

     //---------------------------------------------

     //Fully examine an active node until excess becomes 0

     enum node_examination_t {examine_full, examine_to_relabel};

     //No more implemented yet:, examine_only_one_edge};

     node_examination_t node_examination;

     //Which implementation to be used

     enum implementation_t {impl_fifo, impl_highest_label};

     //No more implemented yet:};

     implementation_t implementation;

     //---------------------------------------------

     //Parameters of the algorithm

     //---------------------------------------------

   private:

     //input

     Graph& G;

     Node s;

     Node t;

     typename Graph::EdgeMap<T> &capacity;

     //output

     typename Graph::EdgeMap<T> preflow;

     T maxflow_value;

     //auxiliary variables for computation

     //The number of the nodes

     int number_of_nodes;

     //A nodemap for the level

     typename Graph::NodeMap<int> level;

     //A nodemap for the excess

     typename Graph::NodeMap<T> excess;

     //Number of nodes on each level

     vector<int> num_of_nodes_on_level;

     //For the FIFO implementation

     list<Node> fifo_nodes;

     //For 'highest label' implementation

     int highest_active;

     //int second_highest_active;

     vector< list<Node> > active_nodes;

   public:

     //Constructing the object using the graph, source, sink and capacity vector

     preflow_push(

 		      Graph& _G, 

 		      Node _s, 

 		      Node _t, 

 		      typename Graph::EdgeMap<T> & _capacity)

       : G(_G), s(_s), t(_t), 

 	capacity(_capacity), 

 	preflow(_G),

 	//Counting the number of nodes

 	//number_of_nodes(count(G.first<EachNodeIt>())),

 	number_of_nodes(G.nodeNum()),

 	level(_G),

 	excess(_G)//,

         // Default constructor: active_nodes()

     { 

       //Simplest parameter settings

       node_examination = examine_full;//examine_to_relabel;//

       //Which implementation to be usedexamine_full

       implementation = impl_highest_label;//impl_fifo;

       //

       num_of_nodes_on_level.resize(2*number_of_nodes-1);

       num_of_nodes_on_level.clear();

       switch(implementation){

       case impl_highest_label :{

 	active_nodes.clear();

 	active_nodes.resize(2*number_of_nodes-1);

 	break;

       }

       default:

 	break;

       }

     }

     //Returns the value of a maximal flow 

     T run();

     typename Graph::EdgeMap<T>  getmaxflow(){

       return preflow;

     }

   private:

     //For testing purposes only

     //Lists the node_properties

     void write_property_vector(typename Graph::NodeMap<T> a,

 			       //node_property_vector<Graph, T> a, 

 			       char* prop_name="property"){

       for(NodeIt i=G.template first<NodeIt>(); G.valid(i); G.next(i)) {

 	cout<<"Node id.: "<<G.id(i)<<", "<<prop_name<<" value: "<<a[i]<<endl;

       }

       cout<<endl;

     }

     //Modifies the excess of the node and makes sufficient changes

     void modify_excess(const Node& a ,T v){

       //T old_value=excess[a];

       excess[a] += v;

     }

     //This private procedure is supposed to modify the preflow on edge j

     //by value v (which can be positive or negative as well) 

     //and maintain the excess on the head and tail

     //Here we do not check whether this is possible or not

     void modify_preflow(Edge j, const T& v){

       //Modifiyng the edge

       preflow[j] += v;

       //Modifiyng the head

       modify_excess(G.head(j),v);

       //Modifiyng the tail

       modify_excess(G.tail(j),-v);

     }

     //Gives the active node to work with 

     //(depending on the implementation to be used)

     Node get_active_node(){

       switch(implementation) {

       case impl_highest_label : {

 	//First need to find the highest label for which there's an active node

 	while( highest_active>=0 && active_nodes[highest_active].empty() ){ 

 	  --highest_active;

 	}

 	if( highest_active>=0) {

 	  Node a=active_nodes[highest_active].front();

 	  active_nodes[highest_active].pop_front();

 	  return a;

 	}

 	else {

 	  return INVALID;

 	}

 	break;

       }

       case impl_fifo : {

 	if( ! fifo_nodes.empty() ) {

 	  Node a=fifo_nodes.front();

 	  fifo_nodes.pop_front();

 	  return a;

 	}

 	else {

 	  return INVALID;

 	}

 	break;

       }

       }

       //

       return INVALID;

     }

     //Puts node 'a' among the active nodes

     void make_active(const Node& a){

       //s and t never become active

       if (a!=s && a!= t){

 	switch(implementation){

 	case impl_highest_label :

 	  active_nodes[level[a]].push_back(a);

 	  break;

 	case impl_fifo :

 	  fifo_nodes.push_back(a);

 	  break;

 	}

       }

       //Update highest_active label

       if (highest_active<level[a]){

 	highest_active=level[a];

       }

     }

     //Changes the level of node a and make sufficent changes

     void change_level_to(Node a, int new_value){

       int seged = level[a];

       level.set(a,new_value);

       --num_of_nodes_on_level[seged];

       ++num_of_nodes_on_level[new_value];

     }

     //Collection of things useful (or necessary) to do before running

     void preprocess(){

       //---------------------------------------

       //Initialize parameters

       //---------------------------------------

       //Setting starting preflow, level and excess values to zero

       //This can be important, if the algorithm is run more then once

       for(NodeIt i=G.template first<NodeIt>(); G.valid(i); G.next(i)) {

         level.set(i,0);

         excess.set(i,0);

 	for(OutEdgeIt j=G.template first<OutEdgeIt>(i); G.valid(j); G.next(j)) 

 	  preflow.set(j, 0);

       }

       num_of_nodes_on_level[0]=number_of_nodes;

       highest_active=0;

       //---------------------------------------

       //Initialize parameters

       //---------------------------------------

       //------------------------------------

       //This is the only part that uses BFS

       //------------------------------------

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

       //Copyright: Jacint

       change_level_to(t,0);

       std::queue<Node> bfs_queue;

       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] == number_of_nodes && 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;

 	    change_level_to(w, l);

 	    //level.set(w, l);

 	  }

 	}

       }

       change_level_to(s,number_of_nodes);

       //level.set(s,number_of_nodes);

       /*

       //Setting starting level values using reverse bfs

       reverse_bfs<Graph> rev_bfs(G,t);

       rev_bfs.run();

       //write_property_vector(rev_bfs.dist,"rev_bfs");

       for(NodeIt i=G.template first<NodeIt>(); G.valid(i); G.next(i)) {

         change_level_to(i,rev_bfs.dist(i));

 	//level.put(i,rev_bfs.dist.get(i));

       }

       */

       //------------------------------------

       //This is the only part that uses BFS

       //------------------------------------

       //Starting level of s

       change_level_to(s,number_of_nodes);

       //level.put(s,number_of_nodes);

       //we push as much preflow from s as possible to start with

       for(OutEdgeIt j=G.template first<OutEdgeIt>(s); G.valid(j); G.next(j)){ 

 	modify_preflow(j,capacity[j] );

 	make_active(G.head(j));

 	int lev=level[G.head(j)];

 	if(highest_active<lev){

 	  highest_active=lev;

 	}

       }

       //cout<<highest_active<<endl;

     } 

     //If the preflow is less than the capacity on the given edge

     //then it is an edge in the residual graph

     bool is_admissible_forward_edge(Edge j, int& new_level){

       if (capacity[j]>preflow[j]){

 	if(level[G.tail(j)]==level[G.head(j)]+1){

 	  return true;

 	}

 	else{

 	  if (level[G.head(j)] < new_level)

 	    new_level=level[G.head(j)];

 	}

       }

       return false;

     }

     //If the preflow is greater than 0 on the given edge

     //then the edge reversd is an edge in the residual graph

     bool is_admissible_backward_edge(Edge j, int& new_level){

       if (0<preflow[j]){

 	if(level[G.tail(j)]==level[G.head(j)]-1){

 	  return true;

 	}

 	else{

 	  if (level[G.tail(j)] < new_level)

 	    new_level=level[G.tail(j)];

 	}

       }

       return false;

     }

   };  //class preflow_push  

   template<typename Graph, typename T>

     T preflow_push<Graph, T>::run() {

     preprocess();

     //write_property_vector(level,"level");

     T e,v;

     Node a;

     while (a=get_active_node(), G.valid(a)){

       //cout<<G.id(a)<<endl;

       //write_property_vector(excess,"excess");

       //write_property_vector(level,"level");

       bool go_to_next_node=false;

       e = excess[a];

       while (!go_to_next_node){

 	//Initial value for the new level for the active node we are dealing with

 	int new_level=2*number_of_nodes;

 	//write_property_vector(excess,"excess");

 	//write_property_vector(level,"level");

 	//cout<<G.id(a)<<endl;

 	//Out edges from node a

 	{

 	  OutEdgeIt j=G.template first<OutEdgeIt>(a);

 	  while (G.valid(j) && e){

 	    if (is_admissible_forward_edge(j,new_level)){

 	      v=min(e,capacity[j] - preflow[j]);

 	      e -= v;

 	      //New node might become active

 	      if (excess[G.head(j)]==0){

 		make_active(G.head(j));

 	      }

 	      modify_preflow(j,v);

 	    }

 	    G.next(j);

 	  }

 	}

 	//In edges to node a

 	{

 	  InEdgeIt j=G.template first<InEdgeIt>(a);

 	  while (G.valid(j) && e){

 	    if (is_admissible_backward_edge(j,new_level)){

 	      v=min(e,preflow[j]);

 	      e -= v;

 	      //New node might become active

 	      if (excess[G.tail(j)]==0){

 		make_active(G.tail(j));

 	      }

 	      modify_preflow(j,-v);

 	    }

 	    G.next(j);

 	  }

 	}

 	//if (G.id(a)==999)

 	//cout<<new_level<<" e: "<<e<<endl;

 	//cout<<G.id(a)<<" "<<new_level<<endl;

 	if (0==e){

 	  //Saturating push

 	  go_to_next_node=true;

 	}

 	else{//If there is still excess in node a

 	  //change_level_to(a,new_level+1);

 	  //Level remains empty

 	  if (num_of_nodes_on_level[level[a]]==1){

 	    change_level_to(a,number_of_nodes);

 	    //go_to_next_node=True;

 	  }

 	  else{

 	    change_level_to(a,new_level+1);

 	    //increase_level(a);

 	  }

 	  switch(node_examination){

 	  case examine_to_relabel:

 	    make_active(a);

 	    go_to_next_node = true;

 	    break;

 	  default:

 	    break;

 	  }

 	}//if (0==e)

       }

     }

     maxflow_value = excess[t];

     return maxflow_value;

   }//run

 }//namespace hugo

 #endif //PREFLOW_PUSH_HH