#ifndef PREFLOW_PUSH_HH

 #define PREFLOW_PUSH_HH

 #include <algorithm>

 #include <list>

 #include <vector>

 //#include "pf_hiba.hh"

 //#include <marci_list_graph.hh>

 #include <marci_graph_traits.hh>

 #include "reverse_bfs.hh"

 using namespace std;

 namespace marci {

   template <typename graph_type, typename T>

   class preflow_push {

     //Hasznos typedef-ek

     typedef graph_traits<graph_type>::node_iterator node_iterator;

     typedef graph_traits<graph_type>::edge_iterator edge_iterator;

     typedef graph_traits<graph_type>::each_node_iterator each_node_iterator;

     typedef graph_traits<graph_type>::each_edge_iterator each_edge_iterator;

     typedef graph_traits<graph_type>::out_edge_iterator out_edge_iterator;

     typedef graph_traits<graph_type>::in_edge_iterator in_edge_iterator;

     typedef graph_traits<graph_type>::sym_edge_iterator sym_edge_iterator;

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

     //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_type& G;

     node_iterator s;

     node_iterator t;

     edge_property_vector<graph_type, T> &capacity;

     //output

     edge_property_vector<graph_type, T> preflow;

     T maxflow_value;

     //auxiliary variables for computation

     int number_of_nodes;

     node_property_vector<graph_type, int> level;

     node_property_vector<graph_type, T> excess;

     //Number of nodes on each level

     vector<int> num_of_nodes_on_level;

     //For the FIFO implementation

     list<node_iterator> fifo_nodes;

     //For 'highest label' implementation

     int highest_active;

     //int second_highest_active;

     vector< list<node_iterator> > active_nodes;

   public:

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

     preflow_push(

 		      graph_type& _G, 

 		      node_iterator _s, 

 		      node_iterator _t, 

 		      edge_property_vector<graph_type, T>& _capacity)

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

 	capacity(_capacity), 

 	preflow(_G),

 	//Counting the number of nodes

 	number_of_nodes(number_of(G.first_node())),

 	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.resize(2*number_of_nodes-1);

 	active_nodes.clear();

 	break;

       }

       default:

 	break;

       }

     }

     //Returns the value of a maximal flow 

     T run();

     edge_property_vector<graph_type, T> getmaxflow(){

       return preflow;

     }

   private:

     //For testing purposes only

     //Lists the node_properties

     void write_property_vector(node_property_vector<graph_type, T> a, 

 			       char* prop_name="property"){

       for(each_node_iterator i=G.first_node(); i.valid(); ++i) {

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

       }

       cout<<endl;

     }

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

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

 	T old_value=excess.get(a);

 	excess.put(a,old_value+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_iterator j, const T& v){

       //Auxiliary variable

       T old_value;

       //Modifiyng the edge

       old_value=preflow.get(j);

       preflow.put(j,old_value+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_iterator get_active_node(){

       //cout<<highest_active<<endl;

       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_iterator a=active_nodes[highest_active].front();

 	  active_nodes[highest_active].pop_front();

 	  return a;

 	}

 	else {

 	  return node_iterator();

 	}

 	break;

       }

       case impl_fifo : {

 	if( ! fifo_nodes.empty() ) {

 	  node_iterator a=fifo_nodes.front();

 	  fifo_nodes.pop_front();

 	  return a;

 	}

 	else {

 	  return node_iterator();

 	}

 	break;

       }

       }

       //

       return node_iterator();

     }

     //Puts node 'a' among the active nodes

     void make_active(const node_iterator& a){

       //s and t never become active

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

 	switch(implementation){

 	case impl_highest_label :

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

 	  break;

 	case impl_fifo :

 	  fifo_nodes.push_back(a);

 	  break;

 	}

       }

       //Update highest_active label

       if (highest_active<level.get(a)){

 	highest_active=level.get(a);

       }

     }

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

     void change_level_to(node_iterator a, int new_value){

       int seged = level.get(a);

       level.put(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(each_node_iterator i=G.first_node(); i.valid(); ++i) {

         level.put(i,0);

         excess.put(i,0);

 	for(out_edge_iterator j=G.first_out_edge(i); j.valid(); ++j) 

 	  preflow.put(j, 0);

       }

       num_of_nodes_on_level[0]=number_of_nodes;

       highest_active=0;

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

       //Initialize parameters

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

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

       //This is the only part that uses BFS

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

       //Setting starting level values using reverse bfs

       reverse_bfs<graph_type> rev_bfs(G,t);

       rev_bfs.run();

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

       for(each_node_iterator i=G.first_node(); i.valid(); ++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(out_edge_iterator j=G.first_out_edge(s); j.valid(); ++j){ 

 	modify_preflow(j,capacity.get(j) );

 	make_active(G.head(j));

 	int lev=level.get(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(out_edge_iterator j, int& new_level){

       if (capacity.get(j)>preflow.get(j)){

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

 	  return true;

 	}

 	else{

 	  if (level.get(G.head(j)) < new_level)

 	    new_level=level.get(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(in_edge_iterator j, int& new_level){

       if (0<preflow.get(j)){

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

 	  return true;

 	}

 	else{

 	  if (level.get(G.tail(j)) < new_level)

 	    new_level=level.get(G.tail(j));

 	}

       }

       return false;

     }

   };  //class preflow_push  

   template<typename graph_type, typename T>

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

     preprocess();

     T e,v;

     node_iterator a;

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

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

       //write_property_vector(excess,"excess");

       //write_property_vector(level,"level");

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

       int new_level=2*number_of_nodes;

       bool go_to_next_node=false;

       e = excess.get(a);

       while (!go_to_next_node){

 	//Out edges from node a

 	{

 	  out_edge_iterator j=G.first_out_edge(a);

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

 	    if (is_admissible_forward_edge(j,new_level)){

 	      v=min(e,capacity.get(j) - preflow.get(j));

 	      e -= v;

 	      //New node might become active

 	      if (excess.get(G.head(j))==0){

 		make_active(G.head(j));

 	      }

 	      modify_preflow(j,v);

 	    }

 	    ++j;

 	  }

 	}

 	//In edges to node a

 	{

 	  in_edge_iterator j=G.first_in_edge(a);

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

 	    if (is_admissible_backward_edge(j,new_level)){

 	      v=min(e,preflow.get(j));

 	      e -= v;

 	      //New node might become active

 	      if (excess.get(G.tail(j))==0){

 		make_active(G.tail(j));

 	      }

 	      modify_preflow(j,-v);

 	    }

 	    ++j;

 	  }

 	}

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

 	  //Level remains empty

 	  if (num_of_nodes_on_level[level.get(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.get(t);

     return maxflow_value;

   }//run

 }//namespace marci

 #endif //PREFLOW_PUSH_HH