// -*- mode:c++ -*-

#ifndef EDMONDS_KARP_HH
#define EDMONDS_KARP_HH

#include <queue>

//#include <marci_property_vector.hh>

#include <algorithm>

//#include <bfs_iterator.hh>

namespace hugo {
  template <typename Graph, typename FlowMap, typename CapacityMap>
  typename FlowMap::ValueType maxFlow(Graph &G,
				      FlowMap &f,
				      CapacityMap &c,
				      typename Graph::NodeIt s,
				      typename Graph::NodeIt t)
  { 
    typedef typename Graph::NodeIt NodeIt;
    typedef typename Graph::EachNodeIt EachNodeIt;
    typedef typename Graph::EdgeIt EdgeIt;
    typedef typename Graph::EachEdgeIt EachEdgeIt;
    typedef typename Graph::OutEdgeIt OutEdgeIt;
    typedef typename Graph::InEdgeIt InEdgeIt;
    typedef typename FlowMap::ValueType T;
    
    T flow_val = 0;
    T aug_val;

    for(EachEdgeIt e(G);G.valid(e);G.goNext(e))
      f.set(e,0);
    
    std::queue<NodeIt> bfs_queue;
    typename Graph::NodeMap<int> visited(G); //0: unvisited,
                                             //1: reached by a forward edge
                                             //2: reached by a backward edge
                                             //3: it is node s
    typename Graph::NodeMap<EdgeIt> tree(G);
    
    NodeIt gn;  //FIXME: it might be too global for some people...
    
  augment:
    
    for(EachNodeIt n(G);G.valid(n);G.goNext(n))
      visited.set(n,0);
    
    visited.set(s,3);
    
    //There must be a better way to do this:
    while(!bfs_queue.empty()) bfs_queue.pop();
    
    bfs_queue.push(s);
    
    while(!bfs_queue.empty() && !visited.get(t))
      {
	NodeIt n(bfs_queue.front());
	bfs_queue.pop();
	for(OutEdgeIt e(G,n);G.valid(e);G.goNext(e))
	  if(f.get(e)<c.get(e) &&   //FIXME: <
	     !visited.get(G.bNode(e))) 
	    {
	      NodeIt nn(G.bNode(e)); //It improves nothing
	      tree.set(nn,e);
	      visited.set(nn,1);
	      bfs_queue.push(nn);
	    }
	for(InEdgeIt e(G,n);G.valid(e);G.goNext(e))
	  if(f.get(e)>0 &&   //FIXME: >
	     !visited.get(G.bNode(e))) 
	    {
	      NodeIt nn(G.bNode(e));
	      tree.set(nn,e);
	      visited.set(nn,2);
	      bfs_queue.push(nn);
	    }
      }
    
    if(!visited.get(t)) return flow_val;

    // Augmenting value computation
    aug_val = visited.get(t)==1 ?
      c.get(tree.get(t))-f.get(tree.get(t)) : f.get(tree.get(t));
    //FIXME: I would need 'G.opposite(e,n)'
    gn = visited.get(t)==1 ? G.tail(tree.get(t)) : G.head(tree.get(t));
    while(gn!=s) if(visited.get(gn)==1)
      {
	//FIXME: nonstandard gcc extension!
	aug_val <?= c.get(tree.get(gn))-f.get(tree.get(gn));
	gn=G.tail(tree.get(gn));
      }
    else {
      //FIXME: nonstandard gcc extension!
      aug_val <?= f.get(tree.get(gn));
      gn=G.head(tree.get(gn));
    }
	
    // The augmentation itself
    gn = t;
    while(gn!=s) if(visited.get(gn)==1)
      {
	f.set(tree.get(gn),f.get(tree.get(gn))+aug_val);
	gn=G.tail(tree.get(gn));
      }
    else {
      f.set(tree.get(gn),f.get(tree.get(gn))-aug_val);
      gn=G.head(tree.get(gn));
    }

    flow_val+=aug_val;

    goto augment;   // Vivat goto forever!
  }
  
} // namespace hugo

#endif //EDMONDS_KARP_HH