Index: src/work/edmonds_karp.h
===================================================================
--- src/work/edmonds_karp.h	(revision 268)
+++ src/work/edmonds_karp.h	(revision 269)
@@ -250,55 +250,54 @@
   template <typename GraphWrapper, typename Number, typename FlowMap, typename CapacityMap>
   class MaxFlow {
-  public:
-    typedef typename GraphWrapper::Node Node;
-    typedef typename GraphWrapper::Edge Edge;
-    typedef typename GraphWrapper::EdgeIt EdgeIt;
-    typedef typename GraphWrapper::OutEdgeIt OutEdgeIt;
-    typedef typename GraphWrapper::InEdgeIt InEdgeIt;
-
-  private:
+  protected:
+    typedef GraphWrapper GW;
+    typedef typename GW::Node Node;
+    typedef typename GW::Edge Edge;
+    typedef typename GW::EdgeIt EdgeIt;
+    typedef typename GW::OutEdgeIt OutEdgeIt;
+    typedef typename GW::InEdgeIt InEdgeIt;
     //const Graph* G;
-    GraphWrapper gw;
+    GW gw;
     Node s;
     Node t;
     FlowMap* flow;
     const CapacityMap* capacity;
-    typedef ResGraphWrapper<GraphWrapper, Number, FlowMap, CapacityMap > 
-    AugGraph;
-    typedef typename AugGraph::OutEdgeIt AugOutEdgeIt;
-    typedef typename AugGraph::Edge AugEdge;
-
+    typedef ResGraphWrapper<GW, Number, FlowMap, CapacityMap > ResGW;
+    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
+    typedef typename ResGW::Edge ResGWEdge;
   public:
-    MaxFlow(const GraphWrapper& _gw, Node _s, Node _t, FlowMap& _flow, const CapacityMap& _capacity) : 
+
+    MaxFlow(const GW& _gw, Node _s, Node _t, FlowMap& _flow, const CapacityMap& _capacity) : 
       gw(_gw), s(_s), t(_t), flow(&_flow), capacity(&_capacity) { }
+
     bool augmentOnShortestPath() {
-      AugGraph res_graph(gw, *flow, *capacity);
+      ResGW res_graph(gw, *flow, *capacity);
       bool _augment=false;
       
-      typedef typename AugGraph::NodeMap<bool> ReachedMap;
-      BfsIterator5< AugGraph, /*AugOutEdgeIt,*/ ReachedMap > res_bfs(res_graph);
-      res_bfs.pushAndSetReached(s);
+      typedef typename ResGW::NodeMap<bool> ReachedMap;
+      BfsIterator5< ResGW, ReachedMap > bfs(res_graph);
+      bfs.pushAndSetReached(s);
 	
-      typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
-      pred.set(s, AugEdge(INVALID));
+      typename ResGW::NodeMap<ResGWEdge> pred(res_graph); 
+      pred.set(s, INVALID);
       
-      typename AugGraph::NodeMap<Number> free(res_graph);
+      typename ResGW::NodeMap<Number> free(res_graph);
 	
       //searching for augmenting path
-      while ( !res_bfs.finished() ) { 
-	AugOutEdgeIt e=res_bfs;
-	if (res_graph.valid(e) && res_bfs.isBNodeNewlyReached()) {
+      while ( !bfs.finished() ) { 
+	ResGWOutEdgeIt e=bfs;
+	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
 	  Node v=res_graph.tail(e);
 	  Node w=res_graph.head(e);
 	  pred.set(w, e);
 	  if (res_graph.valid(pred.get(v))) {
-	    free.set(w, std::min(free.get(v), res_graph.free(e)));
+	    free.set(w, std::min(free.get(v), res_graph.resCap(e)));
 	  } else {
-	    free.set(w, res_graph.free(e)); 
+	    free.set(w, res_graph.resCap(e)); 
 	  }
 	  if (res_graph.head(e)==t) { _augment=true; break; }
 	}
 	
-	++res_bfs;
+	++bfs;
       } //end of searching augmenting path
 
@@ -307,5 +306,5 @@
 	Number augment_value=free.get(t);
 	while (res_graph.valid(pred.get(n))) { 
-	  AugEdge e=pred.get(n);
+	  ResGWEdge e=pred.get(n);
 	  res_graph.augment(e, augment_value); 
 	  n=res_graph.tail(e);
@@ -331,16 +330,17 @@
 
     template<typename MutableGraph> bool augmentOnBlockingFlow() {      
+      typedef MutableGraph MG;
       bool _augment=false;
 
-      AugGraph res_graph(gw, *flow, *capacity);
-
-      typedef typename AugGraph::NodeMap<bool> ReachedMap;
-      BfsIterator5< AugGraph, ReachedMap > bfs(res_graph);
+      ResGW res_graph(gw, *flow, *capacity);
+
+      typedef typename ResGW::NodeMap<bool> ReachedMap;
+      BfsIterator5< ResGW, ReachedMap > bfs(res_graph);
 
       bfs.pushAndSetReached(s);
-      //typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
-      DistanceMap<AugGraph> dist(res_graph);
+      //typename ResGW::NodeMap<int> dist(res_graph); //filled up with 0's
+      DistanceMap<ResGW> dist(res_graph);
       while ( !bfs.finished() ) { 
-	AugOutEdgeIt e=bfs;
+	ResGWOutEdgeIt e=bfs;
 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
@@ -349,28 +349,26 @@
       } //computing distances from s in the residual graph
 
-      MutableGraph F;
-      typedef SubGraphWrapper<AugGraph, DistanceMap<AugGraph> > FilterResGraph;
-      FilterResGraph filter_res_graph(res_graph, dist);
-      typename AugGraph::NodeMap<typename MutableGraph::Node> 
-	res_graph_to_F(res_graph);
-      for(typename AugGraph::NodeIt n=res_graph.template first<typename AugGraph::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
+      MG F;
+      typedef SubGraphWrapper<ResGW, DistanceMap<ResGW> > FilterResGW;
+      FilterResGW filter_res_graph(res_graph, dist);
+      typename ResGW::NodeMap<typename MG::Node> res_graph_to_F(res_graph);
+      for(typename ResGW::NodeIt n=res_graph.template first<typename ResGW::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
 	res_graph_to_F.set(n, F.addNode());
       }
-      
-      typename MutableGraph::Node sF=res_graph_to_F.get(s);
-      typename MutableGraph::Node tF=res_graph_to_F.get(t);
-
-      typename MutableGraph::EdgeMap<AugEdge> original_edge(F);
-      typename MutableGraph::EdgeMap<Number> residual_capacity(F);
+
+      typename MG::Node sF=res_graph_to_F.get(s);
+      typename MG::Node tF=res_graph_to_F.get(t);
+      typename MG::EdgeMap<ResGWEdge> original_edge(F);
+      typename MG::EdgeMap<Number> residual_capacity(F);
 
       //Making F to the graph containing the edges of the residual graph 
       //which are in some shortest paths
-      for(typename FilterResGraph::EdgeIt e=filter_res_graph.template first<typename FilterResGraph::EdgeIt>(); filter_res_graph.valid(e); filter_res_graph.next(e)) {
+      for(typename FilterResGW::EdgeIt e=filter_res_graph.template first<typename FilterResGW::EdgeIt>(); filter_res_graph.valid(e); filter_res_graph.next(e)) {
 	//if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
-	  typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
+	  typename MG::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
 	  original_edge.update();
 	  original_edge.set(f, e);
 	  residual_capacity.update();
-	  residual_capacity.set(f, res_graph.free(e));
+	  residual_capacity.set(f, res_graph.resCap(e));
 	  //} 
       }
@@ -381,19 +379,19 @@
 	__augment=false;
 	//computing blocking flow with dfs
-	typedef typename TrivGraphWrapper<MutableGraph>::NodeMap<bool> BlockingReachedMap;
-	DfsIterator5< TrivGraphWrapper<MutableGraph>, BlockingReachedMap > dfs(F);
-	typename MutableGraph::NodeMap<typename MutableGraph::Edge> pred(F);
-	pred.set(sF, typename MutableGraph::Edge(INVALID));
+	typedef typename TrivGraphWrapper<MG>::NodeMap<bool> BlockingReachedMap;
+	DfsIterator5< TrivGraphWrapper<MG>, BlockingReachedMap > dfs(F);
+	typename MG::NodeMap<typename MG::Edge> pred(F);
+	pred.set(sF, INVALID);
 	//invalid iterators for sources
 
-	typename MutableGraph::NodeMap<Number> free(F);
+	typename MG::NodeMap<Number> free(F);
 
 	dfs.pushAndSetReached(sF);      
 	while (!dfs.finished()) {
 	  ++dfs;
-	  if (F.valid(typename MutableGraph::OutEdgeIt(dfs))) {
+	  if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
 	    if (dfs.isBNodeNewlyReached()) {
-	      typename MutableGraph::Node v=F.aNode(dfs);
-	      typename MutableGraph::Node w=F.bNode(dfs);
+	      typename MG::Node v=F.aNode(dfs);
+	      typename MG::Node w=F.bNode(dfs);
 	      pred.set(w, dfs);
 	      if (F.valid(pred.get(v))) {
@@ -409,5 +407,5 @@
 	      
 	    } else {
-	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
+	      F.erase(/*typename MG::OutEdgeIt*/(dfs));
 	    }
 	  } 
@@ -415,8 +413,8 @@
 
 	if (__augment) {
-	  typename MutableGraph::Node n=tF;
+	  typename MG::Node n=tF;
 	  Number augment_value=free.get(tF);
 	  while (F.valid(pred.get(n))) { 
-	    typename MutableGraph::Edge e=pred.get(n);
+	    typename MG::Edge e=pred.get(n);
 	    res_graph.augment(original_edge.get(e), augment_value); 
 	    n=F.tail(e);
@@ -434,44 +432,45 @@
 
     template<typename MutableGraph> bool augmentOnBlockingFlow1() {      
+      typedef MutableGraph MG;
       bool _augment=false;
 
-      AugGraph res_graph(gw, *flow, *capacity);
+      ResGW res_graph(gw, *flow, *capacity);
 
       //bfs for distances on the residual graph
-      typedef typename AugGraph::NodeMap<bool> ReachedMap;
-      BfsIterator5< AugGraph, ReachedMap > bfs(res_graph);
+      typedef typename ResGW::NodeMap<bool> ReachedMap;
+      BfsIterator5< ResGW, ReachedMap > bfs(res_graph);
       bfs.pushAndSetReached(s);
-      typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
+      typename ResGW::NodeMap<int> dist(res_graph); //filled up with 0's
 
       //F will contain the physical copy of the residual graph
-      //with the set of edges which areon shortest paths
-      MutableGraph F;
-      typename AugGraph::NodeMap<typename MutableGraph::Node> 
-	res_graph_to_F(res_graph);
-      for(typename AugGraph::NodeIt n=res_graph.template first<typename AugGraph::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
+      //with the set of edges which are on shortest paths
+      MG F;
+      typename ResGW::NodeMap<typename MG::Node> res_graph_to_F(res_graph);
+      for(typename ResGW::NodeIt n=res_graph.template first<typename ResGW::NodeIt>(); res_graph.valid(n); res_graph.next(n)) {
 	res_graph_to_F.set(n, F.addNode());
       }
-      typename MutableGraph::Node sF=res_graph_to_F.get(s);
-      typename MutableGraph::Node tF=res_graph_to_F.get(t);
-      typename MutableGraph::EdgeMap<AugEdge> original_edge(F);
-      typename MutableGraph::EdgeMap<Number> residual_capacity(F);
+
+      typename MG::Node sF=res_graph_to_F.get(s);
+      typename MG::Node tF=res_graph_to_F.get(t);
+      typename MG::EdgeMap<ResGWEdge> original_edge(F);
+      typename MG::EdgeMap<Number> residual_capacity(F);
 
       while ( !bfs.finished() ) { 
-	AugOutEdgeIt e=bfs;
+	ResGWOutEdgeIt e=bfs;
 	if (res_graph.valid(e)) {
 	  if (bfs.isBNodeNewlyReached()) {
 	    dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
-	    typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
+	    typename MG::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
 	    original_edge.update();
 	    original_edge.set(f, e);
 	    residual_capacity.update();
-	    residual_capacity.set(f, res_graph.free(e));
+	    residual_capacity.set(f, res_graph.resCap(e));
 	  } else {
 	    if (dist.get(res_graph.head(e))==(dist.get(res_graph.tail(e))+1)) {
-	      typename MutableGraph::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
+	      typename MG::Edge f=F.addEdge(res_graph_to_F.get(res_graph.tail(e)), res_graph_to_F.get(res_graph.head(e)));
 	      original_edge.update();
 	      original_edge.set(f, e);
 	      residual_capacity.update();
-	      residual_capacity.set(f, res_graph.free(e));
+	      residual_capacity.set(f, res_graph.resCap(e));
 	    }
 	  }
@@ -485,19 +484,19 @@
 	__augment=false;
 	//computing blocking flow with dfs
-	typedef typename TrivGraphWrapper<MutableGraph>::NodeMap<bool> BlockingReachedMap;
-	DfsIterator5< TrivGraphWrapper<MutableGraph>/*, typename MutableGraph::OutEdgeIt*/, BlockingReachedMap > dfs(F);
-	typename MutableGraph::NodeMap<typename MutableGraph::Edge> pred(F);
-	pred.set(sF, typename MutableGraph::Edge(INVALID));
+	typedef typename TrivGraphWrapper<MG>::NodeMap<bool> BlockingReachedMap;
+	DfsIterator5< TrivGraphWrapper<MG>, BlockingReachedMap > dfs(F);
+	typename MG::NodeMap<typename MG::Edge> pred(F);
+	pred.set(sF, INVALID);
 	//invalid iterators for sources
 
-	typename MutableGraph::NodeMap<Number> free(F);
+	typename MG::NodeMap<Number> free(F);
 
 	dfs.pushAndSetReached(sF);      
 	while (!dfs.finished()) {
 	  ++dfs;
-	  if (F.valid(typename MutableGraph::OutEdgeIt(dfs))) {
+	  if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
 	    if (dfs.isBNodeNewlyReached()) {
-	      typename MutableGraph::Node v=F.aNode(dfs);
-	      typename MutableGraph::Node w=F.bNode(dfs);
+	      typename MG::Node v=F.aNode(dfs);
+	      typename MG::Node w=F.bNode(dfs);
 	      pred.set(w, dfs);
 	      if (F.valid(pred.get(v))) {
@@ -513,5 +512,5 @@
 	      
 	    } else {
-	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
+	      F.erase(/*typename MG::OutEdgeIt*/(dfs));
 	    }
 	  } 
@@ -519,8 +518,8 @@
 
 	if (__augment) {
-	  typename MutableGraph::Node n=tF;
+	  typename MG::Node n=tF;
 	  Number augment_value=free.get(tF);
 	  while (F.valid(pred.get(n))) { 
-	    typename MutableGraph::Edge e=pred.get(n);
+	    typename MG::Edge e=pred.get(n);
 	    res_graph.augment(original_edge.get(e), augment_value); 
 	    n=F.tail(e);
@@ -533,4 +532,104 @@
 	
       }
+            
+      return _augment;
+    }
+
+    bool augmentOnBlockingFlow2() {
+      bool _augment=false;
+
+      ResGW res_graph(gw, *flow, *capacity);
+
+      typedef typename ResGW::NodeMap<bool> ReachedMap;
+      BfsIterator5< ResGW, ReachedMap > bfs(res_graph);
+
+      bfs.pushAndSetReached(s);
+      DistanceMap<ResGW> dist(res_graph);
+      while ( !bfs.finished() ) { 
+ 	ResGWOutEdgeIt e=bfs;
+ 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
+ 	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
+ 	}
+	++bfs;
+      } //computing distances from s in the residual graph
+
+      //Subgraph containing the edges on some shortest paths
+      typedef SubGraphWrapper<ResGW, DistanceMap<ResGW> > FilterResGW;
+      FilterResGW filter_res_graph(res_graph, dist);
+
+      //Subgraph, which is able to delete edges which are already 
+      //met by the dfs
+      typename FilterResGW::NodeMap<typename FilterResGW::OutEdgeIt> 
+ 	first_out_edges(filter_res_graph);
+      typename FilterResGW::NodeIt v;
+      for(filter_res_graph.first(v); filter_res_graph.valid(v); 
+ 	  filter_res_graph.next(v)) 
+      {
+ 	typename FilterResGW::OutEdgeIt e;
+ 	filter_res_graph.first(e, v);
+ 	first_out_edges.set(v, e);
+      }
+      typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
+	NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
+      ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
+
+      bool __augment=true;
+
+      while (__augment) {
+
+ 	__augment=false;
+ 	//computing blocking flow with dfs
+	typedef typename ErasingResGW::NodeMap<bool> BlockingReachedMap;
+ 	DfsIterator5< ErasingResGW, BlockingReachedMap > 
+ 	  dfs(erasing_res_graph);
+ 	typename ErasingResGW::NodeMap<typename ErasingResGW::OutEdgeIt> 
+ 	  pred(erasing_res_graph); 
+ 	pred.set(s, INVALID);
+ 	//invalid iterators for sources
+
+ 	typename ErasingResGW::NodeMap<Number> free(erasing_res_graph);
+
+ 	dfs.pushAndSetReached(s);
+ 	while (!dfs.finished()) {
+ 	  ++dfs;
+ 	  if (erasing_res_graph.valid(
+ 		/*typename ErasingResGW::OutEdgeIt*/(dfs))) 
+ 	  { 
+ 	    if (dfs.isBNodeNewlyReached()) {
+	  
+ 	      typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
+ 	      typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
+
+ 	      pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
+ 	      if (erasing_res_graph.valid(pred.get(v))) {
+ 		free.set(w, std::min(free.get(v), res_graph.resCap(dfs)));
+ 	      } else {
+ 		free.set(w, res_graph.resCap(dfs)); 
+ 	      }
+	      
+ 	      if (w==t) { 
+ 		__augment=true; 
+ 		_augment=true;
+ 		break; 
+ 	      }
+	    } else {
+	      erasing_res_graph.erase(dfs);
+	    }
+	  }
+	}	
+
+ 	if (__augment) {
+ 	  typename ErasingResGW::Node n=t;
+ 	  Number augment_value=free.get(n);
+ 	  while (erasing_res_graph.valid(pred.get(n))) { 
+ 	    typename ErasingResGW::OutEdgeIt e=pred.get(n);
+ 	    res_graph.augment(e, augment_value);
+ 	    n=erasing_res_graph.tail(e);
+ 	    if (res_graph.resCap(e)==0)
+ 	      erasing_res_graph.erase(e);
+ 	  }
+ 	}
+      
+      } //while (__augment) 
             
       return _augment;
@@ -625,20 +724,23 @@
 //       return _augment;
 //     }
-//     void run() {
-//       //int num_of_augmentations=0;
-//       while (augmentOnShortestPath()) { 
-// 	//while (augmentOnBlockingFlow<MutableGraph>()) { 
-// 	//std::cout << ++num_of_augmentations << " ";
-// 	//std::cout<<std::endl;
-//       } 
-//     }
-//     template<typename MutableGraph> void run() {
-//       //int num_of_augmentations=0;
-//       //while (augmentOnShortestPath()) { 
-// 	while (augmentOnBlockingFlow<MutableGraph>()) { 
-// 	//std::cout << ++num_of_augmentations << " ";
-// 	//std::cout<<std::endl;
-//       } 
-//     }
+
+    void run() {
+      //int num_of_augmentations=0;
+      while (augmentOnShortestPath()) { 
+	//while (augmentOnBlockingFlow<MutableGraph>()) { 
+	//std::cout << ++num_of_augmentations << " ";
+	//std::cout<<std::endl;
+      } 
+    }
+
+    template<typename MutableGraph> void run() {
+      //int num_of_augmentations=0;
+      //while (augmentOnShortestPath()) { 
+	while (augmentOnBlockingFlow<MutableGraph>()) { 
+	//std::cout << ++num_of_augmentations << " ";
+	//std::cout<<std::endl;
+      } 
+    }
+
     Number flowValue() { 
       Number a=0;
@@ -649,4 +751,5 @@
       return a;
     }
+
   };
 
@@ -685,5 +788,5 @@
       
 //       typedef typename AugGraph::NodeMap<bool> ReachedMap;
-//       BfsIterator5< AugGraph, /*AugOutEdgeIt,*/ ReachedMap > res_bfs(res_graph);
+//       BfsIterator5< AugGraph, /*AugOutEdgeIt,*/ ReachedMap > bfs(res_graph);
 //       typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
 //       for(NodeIt s=G->template first<NodeIt>(); G->valid(s); G->next(s)) {
@@ -693,5 +796,5 @@
 // 	  //u+=flow->get(e);
 // 	  //if (u<1) {
-// 	    res_bfs.pushAndSetReached(s);
+// 	    bfs.pushAndSetReached(s);
 // 	    pred.set(s, AugEdge(INVALID));
 // 	    //}
@@ -703,7 +806,7 @@
 //       Node n;
 //       //searching for augmenting path
-//       while ( !res_bfs.finished() ) { 
-// 	AugOutEdgeIt e=res_bfs;
-// 	if (res_graph.valid(e) && res_bfs.isBNodeNewlyReached()) {
+//       while ( !bfs.finished() ) { 
+// 	AugOutEdgeIt e=bfs;
+// 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
 // 	  Node v=res_graph.tail(e);
 // 	  Node w=res_graph.head(e);
@@ -726,5 +829,5 @@
 // 	}
 	
-// 	++res_bfs;
+// 	++bfs;
 //       } //end of searching augmenting path
 
@@ -763,5 +866,5 @@
 // // 	    u+=flow->get(e);
 // // 	  if (u<1) {
-// // 	    res_bfs.pushAndSetReached(s);
+// // 	    bfs.pushAndSetReached(s);
 // // 	    //pred.set(s, AugEdge(INVALID));
 // // 	  }
@@ -1057,10 +1160,10 @@
       
 // //       typedef typename AugGraph::NodeMap<bool> ReachedMap;
-// //       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > res_bfs(res_graph);
+// //       BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
 // //       for(typename std::list<Node>::const_iterator i=S.begin(); 
 // // 	  i!=S.end(); ++i) {
-// // 	res_bfs.pushAndSetReached(*i);
+// // 	bfs.pushAndSetReached(*i);
 // //       }
-// //       //res_bfs.pushAndSetReached(s);
+// //       //bfs.pushAndSetReached(s);
 	
 // //       typename AugGraph::NodeMap<AugEdge> pred(res_graph); 
@@ -1070,7 +1173,7 @@
 	
 // //       //searching for augmenting path
-// //       while ( !res_bfs.finished() ) { 
-// // 	AugOutEdgeIt e=/*AugOutEdgeIt*/(res_bfs);
-// // 	if (e.valid() && res_bfs.isBNodeNewlyReached()) {
+// //       while ( !bfs.finished() ) { 
+// // 	AugOutEdgeIt e=/*AugOutEdgeIt*/(bfs);
+// // 	if (e.valid() && bfs.isBNodeNewlyReached()) {
 // // 	  Node v=res_graph.tail(e);
 // // 	  Node w=res_graph.head(e);
@@ -1088,5 +1191,5 @@
 // // 	}
 	
-// // 	++res_bfs;
+// // 	++bfs;
 // //       } //end of searching augmenting path
 
Index: src/work/marci/edmonds_karp_demo.cc
===================================================================
--- src/work/marci/edmonds_karp_demo.cc	(revision 268)
+++ src/work/marci/edmonds_karp_demo.cc	(revision 269)
@@ -91,5 +91,4 @@
 
   {
-    //std::cout << "SmartGraph..." << std::endl;
     typedef TrivGraphWrapper<const Graph> GW;
     GW gw(G);
@@ -123,5 +122,4 @@
 
   {
-    //std::cout << "SmartGraph..." << std::endl;
     typedef TrivGraphWrapper<const Graph> GW;
     GW gw(G);
@@ -154,30 +152,34 @@
   }
 
-//   {
-//     std::cout << "edmonds karp demo (on-the-fly blocking flow augmentation)..." << std::endl;
-//     Graph::EdgeMap<int> flow(G); //0 flow
-
-//     Timer ts;
-//     ts.reset();
-
-//     MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > max_flow_test(G, s, t, flow, cap);
-//     int i=0;
-//     while (max_flow_test.augmentOnBlockingFlow2()) { 
-// //     for(EdgeIt e=G.template first<EdgeIt>(); e.valid(); ++e) { 
-// //       std::cout<<"("<<G.tail(e)<< "-"<<flow.get(e)<<"->"<<G.head(e)<<") ";
-// //     }
-// //     std::cout<<std::endl;
-//       ++i; 
-//     }
-
-// //   std::cout << "maximum flow: "<< std::endl;
-// //   for(EdgeIt e=G.first<EdgeIt>(); e.valid(); ++e) { 
-// //     std::cout<<"("<<G.tail(e)<< "-"<<flow.get(e)<<"->"<<G.head(e)<<") ";
-// //   }
-// //   std::cout<<std::endl;
-//     std::cout << "elapsed time: " << ts << std::endl;
-//     std::cout << "number of augmentation phases: " << i << std::endl; 
-//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
-//   }
+  {
+    typedef TrivGraphWrapper<const Graph> GW;
+    GW gw(G);
+    std::cout << "edmonds karp demo (on-the-fly blocking flow augmentation)..." << std::endl;
+    GW::EdgeMap<int> flow(G); //0 flow
+
+    Timer ts;
+    ts.reset();
+
+    typedef GW::EdgeMapWrapper< Graph::EdgeMap<int>, int > EMW;
+    EMW cw(cap);
+    MaxFlow<GW, int, GW::EdgeMap<int>, EMW > max_flow_test(gw, s, t, flow, cw);
+    int i=0;
+    while (max_flow_test.augmentOnBlockingFlow2()) { 
+//     for(EdgeIt e=G.template first<EdgeIt>(); e.valid(); ++e) { 
+//       std::cout<<"("<<G.tail(e)<< "-"<<flow.get(e)<<"->"<<G.head(e)<<") ";
+//     }
+//     std::cout<<std::endl;
+      ++i; 
+    }
+
+//   std::cout << "maximum flow: "<< std::endl;
+//   for(EdgeIt e=G.first<EdgeIt>(); e.valid(); ++e) { 
+//     std::cout<<"("<<G.tail(e)<< "-"<<flow.get(e)<<"->"<<G.head(e)<<") ";
+//   }
+//   std::cout<<std::endl;
+    std::cout << "elapsed time: " << ts << std::endl;
+    std::cout << "number of augmentation phases: " << i << std::endl; 
+    std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
+  }
 
   {
@@ -190,5 +192,4 @@
     ts.reset();
 
-    //CM cm;
     typedef GW::EdgeMapWrapper< Graph::EdgeMap<int>, int > EMW;
     EMW cw(cap);
Index: src/work/marci/graph_wrapper.h
===================================================================
--- src/work/marci/graph_wrapper.h	(revision 266)
+++ src/work/marci/graph_wrapper.h	(revision 269)
@@ -1000,9 +1000,9 @@
       OutEdgeIt(const ResGraphWrapper<GraphWrapper, Number, FlowMap, CapacityMap>& resG, Node v) : Edge() { 
 	resG.gw.first(out, v);
-	while( resG.gw.valid(out) && !(resG.free(out)>0) ) { resG.gw.next(out); }
+	while( resG.gw.valid(out) && !(resG.resCap(out)>0) ) { resG.gw.next(out); }
 	if (!resG.gw.valid(out)) {
 	  out_or_in=0;
 	  resG.gw.first(in, v);
-	  while( resG.gw.valid(in) && !(resG.free(in)>0) ) { resG.gw.next(in); }
+	  while( resG.gw.valid(in) && !(resG.resCap(in)>0) ) { resG.gw.next(in); }
 	}
       }
@@ -1012,13 +1012,13 @@
 // 	  Node v=/*resG->*/G->aNode(out);
 // 	  ++out;
-// 	  while( out.valid() && !(Edge::free()>0) ) { ++out; }
+// 	  while( out.valid() && !(Edge::resCap()>0) ) { ++out; }
 // 	  if (!out.valid()) {
 // 	    out_or_in=0;
 // 	    G->first(in, v); 
-// 	    while( in.valid() && !(Edge::free()>0) ) { ++in; }
+// 	    while( in.valid() && !(Edge::resCap()>0) ) { ++in; }
 // 	  }
 // 	} else {
 // 	  ++in;
-// 	  while( in.valid() && !(Edge::free()>0) ) { ++in; } 
+// 	  while( in.valid() && !(Edge::resCap()>0) ) { ++in; } 
 // 	}
 // 	return *this; 
@@ -1038,20 +1038,20 @@
       EdgeIt(const ResGraphWrapper<GraphWrapper, Number, FlowMap, CapacityMap>& resG) : Edge() { 
 	resG.gw.first(v);
-	if (resG.gw.valid(v)) resG.gw.first(out, v); else out=/*OldOutEdgeIt*/(INVALID);
-	while (resG.gw.valid(out) && !(resG.free(out)>0) ) { resG.gw.next(out); }
+	if (resG.gw.valid(v)) resG.gw.first(out, v); else out=INVALID;
+	while (resG.gw.valid(out) && !(resG.resCap(out)>0) ) { resG.gw.next(out); }
 	while (resG.gw.valid(v) && !resG.gw.valid(out)) { 
 	  resG.gw.next(v); 
 	  if (resG.gw.valid(v)) resG.gw.first(out, v); 
-	  while (resG.gw.valid(out) && !(resG.free(out)>0) ) { resG.gw.next(out); }
+	  while (resG.gw.valid(out) && !(resG.resCap(out)>0) ) { resG.gw.next(out); }
 	}
 	if (!resG.gw.valid(out)) {
 	  out_or_in=0;
 	  resG.gw.first(v);
-	  if (resG.gw.valid(v)) resG.gw.first(in, v); else in=/*OldInEdgeIt*/(INVALID);
-	  while (resG.gw.valid(in) && !(resG.free(in)>0) ) { resG.gw.next(in); }
+	  if (resG.gw.valid(v)) resG.gw.first(in, v); else in=INVALID;
+	  while (resG.gw.valid(in) && !(resG.resCap(in)>0) ) { resG.gw.next(in); }
 	  while (resG.gw.valid(v) && !resG.gw.valid(in)) { 
 	    resG.gw.next(v); 
 	    if (resG.gw.valid(v)) resG.gw.first(in, v); 
-	    while (resG.gw.valid(in) && !(resG.free(in)>0) ) { resG.gw.next(in); }
+	    while (resG.gw.valid(in) && !(resG.resCap(in)>0) ) { resG.gw.next(in); }
 	  }
 	}
@@ -1060,9 +1060,9 @@
 // 	if (out_or_in) {
 // 	  ++out;
-// 	  while (out.valid() && !(Edge::free()>0) ) { ++out; }
+// 	  while (out.valid() && !(Edge::resCap()>0) ) { ++out; }
 // 	  while (v.valid() && !out.valid()) { 
 // 	    ++v; 
 // 	    if (v.valid()) G->first(out, v); 
-// 	    while (out.valid() && !(Edge::free()>0) ) { ++out; }
+// 	    while (out.valid() && !(Edge::resCap()>0) ) { ++out; }
 // 	  }
 // 	  if (!out.valid()) {
@@ -1070,18 +1070,18 @@
 // 	    G->first(v);
 // 	    if (v.valid()) G->first(in, v); else in=OldInEdgeIt();
-// 	    while (in.valid() && !(Edge::free()>0) ) { ++in; }
+// 	    while (in.valid() && !(Edge::resCap()>0) ) { ++in; }
 // 	    while (v.valid() && !in.valid()) { 
 // 	      ++v; 
 // 	      if (v.valid()) G->first(in, v); 
-// 	      while (in.valid() && !(Edge::free()>0) ) { ++in; }
+// 	      while (in.valid() && !(Edge::resCap()>0) ) { ++in; }
 // 	    }  
 // 	  }
 // 	} else {
 // 	  ++in;
-// 	  while (in.valid() && !(Edge::free()>0) ) { ++in; }
+// 	  while (in.valid() && !(Edge::resCap()>0) ) { ++in; }
 // 	  while (v.valid() && !in.valid()) { 
 // 	    ++v; 
 // 	    if (v.valid()) G->first(in, v); 
-// 	    while (in.valid() && !(Edge::free()>0) ) { ++in; }
+// 	    while (in.valid() && !(Edge::resCap()>0) ) { ++in; }
 // 	  }
 // 	}
@@ -1106,13 +1106,13 @@
 	Node v=gw.aNode(e.out);
 	gw.next(e.out);
-	while( gw.valid(e.out) && !(free(e.out)>0) ) { gw.next(e.out); }
+	while( gw.valid(e.out) && !(resCap(e.out)>0) ) { gw.next(e.out); }
 	if (!gw.valid(e.out)) {
 	  e.out_or_in=0;
 	  gw.first(e.in, v); 
-	  while( gw.valid(e.in) && !(free(e.in)>0) ) { gw.next(e.in); }
+	  while( gw.valid(e.in) && !(resCap(e.in)>0) ) { gw.next(e.in); }
 	}
       } else {
 	gw.next(e.in);
-	while( gw.valid(e.in) && !(free(e.in)>0) ) { gw.next(e.in); } 
+	while( gw.valid(e.in) && !(resCap(e.in)>0) ) { gw.next(e.in); } 
       }
       return e;
@@ -1122,28 +1122,28 @@
       if (e.out_or_in) {
 	gw.next(e.out);
-	while (gw.valid(e.out) && !(free(e.out)>0) ) { gw.next(e.out); }
+	while (gw.valid(e.out) && !(resCap(e.out)>0) ) { gw.next(e.out); }
 	  while (gw.valid(e.v) && !gw.valid(e.out)) { 
 	    gw.next(e.v); 
 	    if (gw.valid(e.v)) gw.first(e.out, e.v); 
-	    while (gw.valid(e.out) && !(free(e.out)>0) ) { gw.next(e.out); }
+	    while (gw.valid(e.out) && !(resCap(e.out)>0) ) { gw.next(e.out); }
 	  }
 	  if (!gw.valid(e.out)) {
 	    e.out_or_in=0;
 	    gw.first(e.v);
-	    if (gw.valid(e.v)) gw.first(e.in, e.v); else e.in=/*OldInEdgeIt*/(INVALID);
-	    while (gw.valid(e.in) && !(free(e.in)>0) ) { gw.next(e.in); }
+	    if (gw.valid(e.v)) gw.first(e.in, e.v); else e.in=INVALID;
+	    while (gw.valid(e.in) && !(resCap(e.in)>0) ) { gw.next(e.in); }
 	    while (gw.valid(e.v) && !gw.valid(e.in)) { 
 	      gw.next(e.v); 
 	      if (gw.valid(e.v)) gw.first(e.in, e.v); 
-	      while (gw.valid(e.in) && !(free(e.in)>0) ) { gw.next(e.in); }
+	      while (gw.valid(e.in) && !(resCap(e.in)>0) ) { gw.next(e.in); }
 	    }  
 	  }
 	} else {
 	  gw.next(e.in);
-	  while (gw.valid(e.in) && !(free(e.in)>0) ) { gw.next(e.in); }
+	  while (gw.valid(e.in) && !(resCap(e.in)>0) ) { gw.next(e.in); }
 	  while (gw.valid(e.v) && !gw.valid(e.in)) { 
 	    gw.next(e.v); 
 	    if (gw.valid(e.v)) gw.first(e.in, e.v); 
-	    while (gw.valid(e.in) && !(free(e.in)>0) ) { gw.next(e.in); }
+	    while (gw.valid(e.in) && !(resCap(e.in)>0) ) { gw.next(e.in); }
 	  }
 	}
@@ -1194,5 +1194,5 @@
     }
 
-    Number free(const Edge& e) const { 
+    Number resCap(const Edge& e) const { 
       if (e.out_or_in) 
 	return (capacity->get(e.out)-flow->get(e.out)); 
@@ -1201,9 +1201,9 @@
     }
 
-    Number free(OldOutEdgeIt out) const { 
+    Number resCap(OldOutEdgeIt out) const { 
       return (capacity->get(out)-flow->get(out)); 
     }
     
-    Number free(OldInEdgeIt in) const { 
+    Number resCap(OldInEdgeIt in) const { 
       return (flow->get(in)); 
     }
@@ -1246,4 +1246,57 @@
       }
     };
+  };
+
+  //Subgraph on the same node-set and partial edge-set
+  template<typename GraphWrapper, typename FirstOutEdgesMap>
+  class ErasingFirstGraphWrapper : public GraphWrapperSkeleton<GraphWrapper> {
+  protected:
+    FirstOutEdgesMap* first_out_edges;
+  public:
+    typedef typename GraphWrapperSkeleton<GraphWrapper>::Node Node;
+    typedef typename GraphWrapperSkeleton<GraphWrapper>::NodeIt NodeIt;
+    typedef typename GraphWrapperSkeleton<GraphWrapper>::Edge Edge;
+    typedef typename GraphWrapperSkeleton<GraphWrapper>::EdgeIt EdgeIt;
+    typedef typename GraphWrapperSkeleton<GraphWrapper>::InEdgeIt InEdgeIt;
+    typedef typename GraphWrapperSkeleton<GraphWrapper>::OutEdgeIt OutEdgeIt;
+
+    ErasingFirstGraphWrapper(GraphWrapper _gw, FirstOutEdgesMap& _first_out_edges) : 
+      GraphWrapperSkeleton<GraphWrapper>(_gw), first_out_edges(&_first_out_edges) { }  
+
+    template<typename I> I& first(I& i) const { 
+      gw.first(i); 
+      //while (gw.valid(i) && !filter_map->get(i)) { gw.next(i); }
+      return i;
+    }
+    OutEdgeIt& first(OutEdgeIt& e, const Node& n) const {
+      e=first_out_edges->get(n);
+      return e;
+    }
+    template<typename I, typename P> I& first(I& i, const P& p) const { 
+      gw.first(i, p); 
+      //while (gw.valid(i) && !filter_map->get(i)) { gw.next(i); }
+      return i;
+    }
+    
+    //template<typename I> I getNext(const I& i) const { 
+    //  return gw.getNext(i); 
+    //}
+    template<typename I> I& next(I &i) const { 
+      gw.next(i); 
+      //while (gw.valid(i) && !filter_map->get(i)) { gw.next(i); }
+      return i;
+    }
+    
+    template< typename It > It first() const { 
+      It e; this->first(e); return e; }
+    
+    template< typename It > It first(const Node& v) const { 
+      It e; this->first(e, v); return e; }
+
+    void erase(const OutEdgeIt& e) const {
+      OutEdgeIt f=e;
+      this->next(f);
+      first_out_edges->set(this->tail(e), f);
+    }
   };
 
