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

#ifndef EDMONDS_KARP_HH

#define EDMONDS_KARP_HH

#include <queue>

//#include <marci_property_vector.hh>

#include <algorithm>

//#include <bfs_iterator.hh>

namespace marci {

  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::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.next(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(NodeIt n(G);G.valid(n);G.next(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());

	for(OutEdgeIt e(G,n);G.valid(e);G.next(e))

	  if(f.get(e)<c.get(e) &&   //FIXME: <

	     !visited.get(G.bNode(e))) 

	    {

	      tree.set(G.bNode(e),e);

	      visited.set(G.bNode(e),1);

	    }

	for(InEdgeIt e(G,n);G.valid(e);G.next(e))

	  if(f.get(e)>0 &&   //FIXME: >

	     !visited.get(G.bNode(e))) 

	    {

	      tree.set(G.bNode(e),e);

	      visited.set(G.bNode(e),2);

	    }

      }

    if(!visited.get(t)) return flow_val;

    // Augmenting value computation

    aug_val = visited.get(t)==2 ?

      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)==2 ? G.from(tree.get(t)) : G.to(tree.get(t));

    while(gn!=s) if(visited.get(gn)==2)

      {

	//FIXME: nonstandars. gcc extension!

	aug_val <?= c.get(tree.get(gn))-f.get(tree.get(gn));

	gn=G.from(tree.get(gn));

      }

    else {

      //FIXME: nonstandars. gcc extension!

      aug_val <?= f.get(tree.get(gn));

      gn=G.from(tree.get(gn));

    }

    // The augmentation itself

    gn = t;

    while(gn!=s) if(visited.get(gn)==2)

      {

	f.set(tree.get(gn),f.get(tree.get(gn))+aug_val);

	gn=G.from(tree.get(gn));

      }

    else {

      f.set(tree.get(gn),f.get(tree.get(gn))-aug_val);

      gn=G.from(tree.get(gn));

    }

    flow_val+=aug_val;

    goto augment;   // Vivat goto forever!

  }

} // namespace marci

#endif //EDMONDS_KARP_HH