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