Index: src/work/edmonds_karp.hh
===================================================================
--- src/work/edmonds_karp.hh	(revision 69)
+++ src/work/edmonds_karp.hh	(revision 75)
@@ -3,4 +3,6 @@
 
 #include <algorithm>
+#include <list>
+#include <iterator>
 
 #include <bfs_iterator.hh>
@@ -8,5 +10,5 @@
 namespace marci {
 
-  template<typename Graph, typename T, typename FlowMap, typename CapacityMap>
+  template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
   class ResGraph {
     typedef typename Graph::NodeIt NodeIt;
@@ -27,12 +29,12 @@
 
     class EdgeIt {
-      friend class ResGraph<Graph, T, FlowMap, CapacityMap>;
+      friend class ResGraph<Graph, Number, FlowMap, CapacityMap>;
     protected:
-      const ResGraph<Graph, T, FlowMap, CapacityMap>* resG;
+      const ResGraph<Graph, Number, FlowMap, CapacityMap>* resG;
       OldSymEdgeIt sym;
     public:
       EdgeIt() { } 
       //EdgeIt(const EdgeIt& e) : resG(e.resG), sym(e.sym) { }
-      T free() const { 
+      Number free() const { 
 	if (resG->G.aNode(sym)==resG->G.tail(sym)) { 
 	  return (resG->capacity.get(sym)-resG->flow.get(sym)); 
@@ -42,9 +44,11 @@
       }
       bool valid() const { return sym.valid(); }
-      void augment(T a) const {
+      void augment(Number a) const {
 	if (resG->G.aNode(sym)==resG->G.tail(sym)) { 
 	  resG->flow.set(sym, resG->flow.get(sym)+a);
+	  //resG->flow[sym]+=a;
 	} else { 
 	  resG->flow.set(sym, resG->flow.get(sym)-a);
+	  //resG->flow[sym]-=a;
 	}
       }
@@ -52,10 +56,10 @@
 
     class OutEdgeIt : public EdgeIt {
-      friend class ResGraph<Graph, T, FlowMap, CapacityMap>;
+      friend class ResGraph<Graph, Number, FlowMap, CapacityMap>;
     public:
       OutEdgeIt() { }
       //OutEdgeIt(const OutEdgeIt& e) { resG=e.resG; sym=e.sym; }
     private:
-      OutEdgeIt(const ResGraph<Graph, T, FlowMap, CapacityMap>& _resG, NodeIt v) { 
+      OutEdgeIt(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _resG, NodeIt v) { 
       	resG=&_resG;
 	sym=resG->G.template first<OldSymEdgeIt>(v);
@@ -66,4 +70,129 @@
 	++sym; 
 	while( sym.valid() && !(free()>0) ) { ++sym; }
+	return *this; 
+      }
+    };
+
+    void getFirst(OutEdgeIt& e, NodeIt v) const { 
+      e=OutEdgeIt(*this, v); 
+    }
+    void getFirst(EachNodeIt& v) const { G.getFirst(v); }
+    
+    template< typename It >
+    It first() const { 
+      It e;      
+      getFirst(e);
+      return e; 
+    }
+
+    template< typename It >
+    It first(NodeIt v) const { 
+      It e;
+      getFirst(e, v);
+      return e; 
+    }
+
+    NodeIt tail(EdgeIt e) const { return G.aNode(e.sym); }
+    NodeIt head(EdgeIt e) const { return G.bNode(e.sym); }
+
+    NodeIt aNode(OutEdgeIt e) const { return G.aNode(e.sym); }
+    NodeIt bNode(OutEdgeIt e) const { return G.bNode(e.sym); }
+
+    int id(NodeIt v) const { return G.id(v); }
+
+    template <typename S>
+    class NodeMap {
+      typename Graph::NodeMap<S> node_map; 
+    public:
+      NodeMap(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _G) : node_map(_G.G) { }
+      NodeMap(const ResGraph<Graph, Number, FlowMap, CapacityMap>& _G, S a) : node_map(_G.G, a) { }
+      void set(NodeIt nit, S a) { node_map.set(nit, a); }
+      S get(NodeIt nit) const { return node_map.get(nit); }
+      S& operator[](NodeIt nit) { return node_map[nit]; } 
+      const S& operator[](NodeIt nit) const { return node_map[nit]; } 
+    };
+
+  };
+
+
+  template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
+  class ResGraph2 {
+    typedef typename Graph::NodeIt NodeIt;
+    typedef typename Graph::EachNodeIt EachNodeIt;
+    //typedef typename Graph::SymEdgeIt OldSymEdgeIt;
+    typedef typename Graph::OutEdgeIt OldOutEdgeIt;
+    typedef typename Graph::InEdgeIt OldInEdgeIt;
+    
+    const Graph& G;
+    FlowMap& flow;
+    const CapacityMap& capacity;
+  public:
+    ResGraph2(const Graph& _G, FlowMap& _flow, 
+	     const CapacityMap& _capacity) : 
+      G(_G), flow(_flow), capacity(_capacity) { }
+
+    class EdgeIt; 
+    class OutEdgeIt; 
+    friend class EdgeIt; 
+    friend class OutEdgeIt; 
+
+    class EdgeIt {
+      friend class ResGraph2<Graph, Number, FlowMap, CapacityMap>;
+    protected:
+      const ResGraph2<Graph, Number, FlowMap, CapacityMap>* resG;
+      //OldSymEdgeIt sym;
+      OldOutEdgeIt out;
+      OldInEdgeIt in;
+      bool out_or_in; //1, iff out
+    public:
+      EdgeIt() : out_or_in(1) { } 
+      Number free() const { 
+	if (out_or_in) { 
+	  return (resG->capacity.get(out)-resG->flow.get(out)); 
+	} else { 
+	  return (resG->flow.get(in)); 
+	}
+      }
+      bool valid() const { 
+	return out_or_in && out.valid() || in.valid(); }
+      void augment(Number a) const {
+	if (out_or_in) { 
+	  resG->flow.set(out, resG->flow.get(out)+a);
+	} else { 
+	  resG->flow.set(in, resG->flow.get(in)-a);
+	}
+      }
+    };
+
+    class OutEdgeIt : public EdgeIt {
+      friend class ResGraph2<Graph, Number, FlowMap, CapacityMap>;
+    public:
+      OutEdgeIt() { }
+    private:
+      OutEdgeIt(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _resG, NodeIt v) { 
+      	resG=&_resG;
+	out=resG->G.template first<OldOutEdgeIt>(v);
+	while( out.valid() && !(free()>0) ) { ++out; }
+	if (!out.valid()) {
+	  out_or_in=0;
+	  in=resG->G.template first<OldInEdgeIt>(v);
+	  while( in.valid() && !(free()>0) ) { ++in; }
+	}
+      }
+    public:
+      OutEdgeIt& operator++() { 
+	if (out_or_in) {
+	  NodeIt v=resG->G.aNode(out);
+	  ++out;
+	  while( out.valid() && !(free()>0) ) { ++out; }
+	  if (!out.valid()) {
+	    out_or_in=0;
+	    in=resG->G.template first<OldInEdgeIt>(v);
+	    while( in.valid() && !(free()>0) ) { ++in; }
+	  }
+	} else {
+	  ++in;
+	  while( in.valid() && !(free()>0) ) { ++in; } 
+	}
 	return *this; 
       }
@@ -89,9 +218,13 @@
     }
 
-    NodeIt tail(EdgeIt e) const { return G.aNode(e.sym); }
-    NodeIt head(EdgeIt e) const { return G.bNode(e.sym); }
-
-    NodeIt aNode(OutEdgeIt e) const { return G.aNode(e.sym); }
-    NodeIt bNode(OutEdgeIt e) const { return G.bNode(e.sym); }
+    NodeIt tail(EdgeIt e) const { 
+      return ((e.out_or_in) ? G.aNode(e.out) : G.aNode(e.in)); }
+    NodeIt head(EdgeIt e) const { 
+      return ((e.out_or_in) ? G.bNode(e.out) : G.bNode(e.in)); }
+
+    NodeIt aNode(OutEdgeIt e) const { 
+      return ((e.out_or_in) ? G.aNode(e.out) : G.aNode(e.in)); }
+    NodeIt bNode(OutEdgeIt e) const { 
+      return ((e.out_or_in) ? G.bNode(e.out) : G.bNode(e.in)); }
 
     int id(NodeIt v) const { return G.id(v); }
@@ -101,13 +234,144 @@
       typename Graph::NodeMap<S> node_map; 
     public:
-      NodeMap(const ResGraph<Graph, T, FlowMap, CapacityMap>& _G) : node_map(_G.G) { }
-      NodeMap(const ResGraph<Graph, T, FlowMap, CapacityMap>& _G, S a) : node_map(_G.G, a) { }
+      NodeMap(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _G) : node_map(_G.G) { }
+      NodeMap(const ResGraph2<Graph, Number, FlowMap, CapacityMap>& _G, S a) : node_map(_G.G, a) { }
       void set(NodeIt nit, S a) { node_map.set(nit, a); }
       S get(NodeIt nit) const { return node_map.get(nit); }
     };
-
   };
 
-  template <typename Graph, typename T, typename FlowMap, typename CapacityMap>
+  template<typename Graph, typename Number, typename FlowMap, typename CapacityMap>
+  class ResGraph3 {
+    typedef typename Graph::NodeIt NodeIt;
+    typedef typename Graph::EachNodeIt EachNodeIt;
+    //typedef typename Graph::SymEdgeIt OldSymEdgeIt;
+    typedef typename Graph::OutEdgeIt OldOutEdgeIt;
+    typedef typename Graph::InEdgeIt OldInEdgeIt;
+    
+    const Graph& G;
+    FlowMap& flow;
+    const CapacityMap& capacity;
+  public:
+    ResGraph3(const Graph& _G, FlowMap& _flow, 
+	     const CapacityMap& _capacity) : 
+      G(_G), flow(_flow), capacity(_capacity) { }
+
+    class EdgeIt; 
+    class OutEdgeIt; 
+    friend class EdgeIt; 
+    friend class OutEdgeIt; 
+
+    class EdgeIt {
+      friend class ResGraph3<Graph, Number, FlowMap, CapacityMap>;
+    protected:
+      //const ResGraph3<Graph, Number, FlowMap, CapacityMap>* resG;
+      const Graph* G;
+      FlowMap* flow;
+      const CapacityMap* capacity;
+      //OldSymEdgeIt sym;
+      OldOutEdgeIt out;
+      OldInEdgeIt in;
+      bool out_or_in; //1, iff out
+    public:
+      EdgeIt() : out_or_in(1) { } 
+      Number free() const { 
+	if (out_or_in) { 
+	  return (/*resG->*/capacity->get(out)-/*resG->*/flow->get(out)); 
+	} else { 
+	  return (/*resG->*/flow->get(in)); 
+	}
+      }
+      bool valid() const { 
+	return out_or_in && out.valid() || in.valid(); }
+      void augment(Number a) const {
+	if (out_or_in) { 
+	  /*resG->*/flow->set(out, /*resG->*/flow->get(out)+a);
+	} else { 
+	  /*resG->*/flow->set(in, /*resG->*/flow->get(in)-a);
+	}
+      }
+    };
+
+    class OutEdgeIt : public EdgeIt {
+      friend class ResGraph3<Graph, Number, FlowMap, CapacityMap>;
+    public:
+      OutEdgeIt() { }
+    private:
+      OutEdgeIt(const Graph& _G, NodeIt v, FlowMap& _flow, const CapacityMap& _capacity) { 
+      	G=&_G;
+	flow=&_flow;
+	capacity=&_capacity;
+	out=/*resG->*/G->template first<OldOutEdgeIt>(v);
+	while( out.valid() && !(free()>0) ) { ++out; }
+	if (!out.valid()) {
+	  out_or_in=0;
+	  in=/*resG->*/G->template first<OldInEdgeIt>(v);
+	  while( in.valid() && !(free()>0) ) { ++in; }
+	}
+      }
+    public:
+      OutEdgeIt& operator++() { 
+	if (out_or_in) {
+	  NodeIt v=/*resG->*/G->aNode(out);
+	  ++out;
+	  while( out.valid() && !(free()>0) ) { ++out; }
+	  if (!out.valid()) {
+	    out_or_in=0;
+	    in=/*resG->*/G->template first<OldInEdgeIt>(v);
+	    while( in.valid() && !(free()>0) ) { ++in; }
+	  }
+	} else {
+	  ++in;
+	  while( in.valid() && !(free()>0) ) { ++in; } 
+	}
+	return *this; 
+      }
+    };
+
+    void getFirst(OutEdgeIt& e, NodeIt v) const { 
+      e=OutEdgeIt(G, v, flow, capacity); 
+    }
+    void getFirst(EachNodeIt& v) const { G.getFirst(v); }
+    
+    template< typename It >
+    It first() const { 
+      It e;
+      getFirst(e);
+      return e; 
+    }
+
+    template< typename It >
+    It first(NodeIt v) const { 
+      It e;
+      getFirst(e, v);
+      return e; 
+    }
+
+    NodeIt tail(EdgeIt e) const { 
+      return ((e.out_or_in) ? G.aNode(e.out) : G.aNode(e.in)); }
+    NodeIt head(EdgeIt e) const { 
+      return ((e.out_or_in) ? G.bNode(e.out) : G.bNode(e.in)); }
+
+    NodeIt aNode(OutEdgeIt e) const { 
+      return ((e.out_or_in) ? G.aNode(e.out) : G.aNode(e.in)); }
+    NodeIt bNode(OutEdgeIt e) const { 
+      return ((e.out_or_in) ? G.bNode(e.out) : G.bNode(e.in)); }
+
+    int id(NodeIt v) const { return G.id(v); }
+
+    template <typename S>
+    class NodeMap {
+      typename Graph::NodeMap<S> node_map; 
+    public:
+      NodeMap(const ResGraph3<Graph, Number, FlowMap, CapacityMap>& _G) : node_map(_G.G) { }
+      NodeMap(const ResGraph3<Graph, Number, FlowMap, CapacityMap>& _G, S a) : node_map(_G.G, a) { }
+      void set(NodeIt nit, S a) { node_map.set(nit, a); }
+      S get(NodeIt nit) const { return node_map.get(nit); }
+    };
+
+  };
+
+
+  template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
   class MaxFlow {
     typedef typename Graph::NodeIt NodeIt;
@@ -121,9 +385,17 @@
     FlowMap& flow;
     const CapacityMap& capacity;
-    typedef ResGraph<Graph, T, FlowMap, CapacityMap > AugGraph;
+    typedef ResGraph3<Graph, Number, FlowMap, CapacityMap > AugGraph;
     typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
     typedef typename AugGraph::EdgeIt AugEdgeIt;
+
+    //AugGraph res_graph;    
+    //typedef typename AugGraph::NodeMap<bool> ReachedMap;
+    //typename AugGraph::NodeMap<AugEdgeIt> pred; 
+    //typename AugGraph::NodeMap<int> free;
   public:
-    MaxFlow(const Graph& _G, NodeIt _s, NodeIt _t, FlowMap& _flow, const CapacityMap& _capacity) : G(_G), s(_s), t(_t), flow(_flow), capacity(_capacity) { }
+    MaxFlow(const Graph& _G, NodeIt _s, NodeIt _t, FlowMap& _flow, const CapacityMap& _capacity) : 
+      G(_G), s(_s), t(_t), flow(_flow), capacity(_capacity) //,  
+      //res_graph(G, flow, capacity), pred(res_graph), free(res_graph) 
+      { }
     bool augment() {
       AugGraph res_graph(G, flow, capacity);
@@ -131,9 +403,10 @@
       
       typedef typename AugGraph::NodeMap<bool> ReachedMap;
-      BfsIterator2< AugGraph, AugOutEdgeIt, ReachedMap > res_bfs(res_graph);
+      BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > res_bfs(res_graph);
       res_bfs.pushAndSetReached(s);
 	
       typename AugGraph::NodeMap<AugEdgeIt> pred(res_graph); 
-      //filled with invalid iterators
+      //filled up with invalid iterators
+      //pred.set(s, AugEdgeIt());
       
       typename AugGraph::NodeMap<int> free(res_graph);
@@ -159,5 +432,5 @@
       if (_augment) {
 	NodeIt n=t;
-	T augment_value=free.get(t);
+	Number augment_value=free.get(t);
 	while (pred.get(n).valid()) { 
 	  AugEdgeIt e=pred.get(n);
@@ -172,6 +445,6 @@
       while (augment()) { } 
     }
-    T flowValue() { 
-      T a=0;
+    Number flowValue() { 
+      Number a=0;
       for(OutEdgeIt i=G.template first<OutEdgeIt>(s); i.valid(); ++i) {
 	a+=flow.get(i);
@@ -181,4 +454,151 @@
   };
 
+
+/*
+  template <typename Graph>
+  class IteratorBoolNodeMap {
+    typedef typename Graph::NodeIt NodeIt;
+    typedef typename Graph::EachNodeIt EachNodeIt;
+    class BoolItem {
+    public:
+      bool value;
+      NodeIt prev;
+      NodeIt next;
+      BoolItem() : value(false), prev(), next() {}
+    };
+    NodeIt first_true;
+    //NodeIt last_true;
+    NodeIt first_false;
+    //NodeIt last_false;
+    const Graph& G; 
+    typename Graph::NodeMap<BoolItem> container;
+  public:
+    typedef bool ValueType;
+    typedef NodeIt KeyType;
+    IteratorBoolNodeMap(const Graph& _G) : G(_G), container(G), first_true(NodeIt()) { 
+      //for (EachNodeIt e=G.template first<EacNodeIt>(); e.valid(); ++e) {
+      //BoolItem b=container.get(e);
+      //b.me=e;
+      //container.set(b);
+      //} 
+      G.getFirst(first_false);
+      NodeIt e_prev;
+      for (EachNodeIt e=G.template first<EacNodeIt>(); e.valid(); ++e) {
+	container[e].prev=e_prev;
+	if (e_prev.valid()) container[e_prev].next=e;
+	e_prev=e;
+      }
+    }
+    //NodeMap(const Graph& _G, T a) : 
+    //  G(_G), container(G.node_id, a) { }
+    //FIXME
+    void set(NodeIt nit, T a) { container[G.id(nit)]=a; }
+    T get(NodeIt nit) const { return container[G.id(nit)]; }
+    //void resize() { container.resize(G.node_id); }
+    //void resize(T a) { container.resize(G.node_id, a); }
+  };
+*/
+
+  
+  template <typename Graph, typename Number, typename FlowMap, typename CapacityMap>
+  class MaxFlow2 {
+    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;
+    const Graph& G;
+    std::list<NodeIt>& S;
+    std::list<NodeIt>& T;
+    FlowMap& flow;
+    const CapacityMap& capacity;
+    typedef ResGraph<Graph, Number, FlowMap, CapacityMap > AugGraph;
+    typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
+    typedef typename AugGraph::EdgeIt AugEdgeIt;
+    typename Graph::NodeMap<bool> SMap;
+    typename Graph::NodeMap<bool> TMap;
+  public:
+    MaxFlow2(const Graph& _G, std::list<NodeIt>& _S, std::list<NodeIt>& _T, FlowMap& _flow, const CapacityMap& _capacity) : G(_G), S(_S), T(_T), flow(_flow), capacity(_capacity), SMap(_G), TMap(_G) { 
+      for(typename std::list<NodeIt>::const_iterator i=S.begin(); 
+	  i!=S.end(); ++i) { 
+	SMap.set(*i, true); 
+      }
+      for (typename std::list<NodeIt>::const_iterator i=T.begin(); 
+	   i!=T.end(); ++i) { 
+	TMap.set(*i, true); 
+      }
+    }
+    bool augment() {
+      AugGraph res_graph(G, flow, capacity);
+      bool _augment=false;
+      NodeIt reached_t_node;
+      
+      typedef typename AugGraph::NodeMap<bool> ReachedMap;
+      BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > res_bfs(res_graph);
+      for(typename std::list<NodeIt>::const_iterator i=S.begin(); 
+	  i!=S.end(); ++i) {
+	res_bfs.pushAndSetReached(*i);
+      }
+      //res_bfs.pushAndSetReached(s);
+	
+      typename AugGraph::NodeMap<AugEdgeIt> pred(res_graph); 
+      //filled up with invalid iterators
+      
+      typename AugGraph::NodeMap<int> free(res_graph);
+	
+      //searching for augmenting path
+      while ( !res_bfs.finished() ) { 
+	AugOutEdgeIt e=AugOutEdgeIt(res_bfs);
+	if (e.valid() && res_bfs.isBNodeNewlyReached()) {
+	  NodeIt v=res_graph.tail(e);
+	  NodeIt w=res_graph.head(e);
+	  pred.set(w, e);
+	  if (pred.get(v).valid()) {
+	    free.set(w, std::min(free.get(v), e.free()));
+	  } else {
+	    free.set(w, e.free()); 
+	  }
+	  if (TMap.get(res_graph.head(e))) { 
+	    _augment=true; 
+	    reached_t_node=res_graph.head(e);
+	    break; 
+	  }
+	}
+	
+	++res_bfs;
+      } //end of searching augmenting path
+
+      if (_augment) {
+	NodeIt n=reached_t_node;
+	Number augment_value=free.get(reached_t_node);
+	while (pred.get(n).valid()) { 
+	  AugEdgeIt e=pred.get(n);
+	  e.augment(augment_value); 
+	  n=res_graph.tail(e);
+	}
+      }
+
+      return _augment;
+    }
+    void run() {
+      while (augment()) { } 
+    }
+    Number flowValue() { 
+      Number a=0;
+      for(typename std::list<NodeIt>::const_iterator i=S.begin(); 
+	  i!=S.end(); ++i) { 
+	for(OutEdgeIt e=G.template first<OutEdgeIt>(*i); e.valid(); ++e) {
+	  a+=flow.get(e);
+	}
+	for(InEdgeIt e=G.template first<InEdgeIt>(*i); e.valid(); ++e) {
+	  a-=flow.get(e);
+	}
+      }
+      return a;
+    }
+  };
+
+
+
 } // namespace marci
 
