# HG changeset patch
# User marci
# Date 1077885555 0
# Node ID 0631992fe7a19c0767457700e5c8edf24a84cbcc
# Parent  1ac27e476e252aab8e9438eb159f4a404e0acd05
Dinits blocking flow added to edmonds_karp_demo.hh.

diff -r 1ac27e476e25 -r 0631992fe7a1 src/work/bfs_iterator.hh
--- a/src/work/bfs_iterator.hh	Thu Feb 26 16:07:40 2004 +0000
+++ b/src/work/bfs_iterator.hh	Fri Feb 27 12:39:15 2004 +0000
@@ -644,6 +644,7 @@
       //}
     }
     void pushAndSetReached(NodeIt s) { 
+      actual_node=s;
       reached.set(s, true);
       dfs_stack.push(G.template first<OutEdgeIt>(s)); 
     }
@@ -659,6 +660,7 @@
 	  reached.set(w, true);
 	  b_node_newly_reached=true;
 	} else {
+	  actual_node=G.aNode(actual_edge);
 	  ++(dfs_stack.top());
 	  b_node_newly_reached=false;
 	}
@@ -672,7 +674,7 @@
     operator OutEdgeIt () const { return actual_edge; }
     bool isBNodeNewlyReached() const { return b_node_newly_reached; }
     bool isANodeExamined() const { return !(actual_edge.valid()); }
-    NodeIt aNode() const { return actual_node; }
+    NodeIt aNode() const { return actual_node; /*FIXME*/}
     NodeIt bNode() const { return G.bNode(actual_edge); }
     const ReachedMap& getReachedMap() const { return reached; }
     const std::stack<OutEdgeIt>& getDfsStack() const { return dfs_stack; }
diff -r 1ac27e476e25 -r 0631992fe7a1 src/work/edmonds_karp.hh
--- a/src/work/edmonds_karp.hh	Thu Feb 26 16:07:40 2004 +0000
+++ b/src/work/edmonds_karp.hh	Fri Feb 27 12:39:15 2004 +0000
@@ -6,7 +6,7 @@
 #include <iterator>
 
 #include <bfs_iterator.hh>
-#include <time_measure.h>
+//#include <time_measure.h>
 
 namespace hugo {
 
@@ -281,6 +281,7 @@
       bool out_or_in; //1, iff out
     public:
       EdgeIt() : out_or_in(1) { } 
+      //EdgeIt(const EdgeIt& e) : G(e.G), flow(e.flow), capacity(e.capacity), out(e.out), in(e.in), out_or_in(e.out_or_in) { }
       Number free() const { 
 	if (out_or_in) { 
 	  return (/*resG->*/capacity->get(out)-/*resG->*/flow->get(out)); 
@@ -297,6 +298,23 @@
 	  /*resG->*/flow->set(in, /*resG->*/flow->get(in)-a);
 	}
       }
+      void print() { 
+	if (out_or_in) {
+	  std::cout << "out "; 
+	  if (out.valid()) 
+	    std::cout << G->id(G->tail(out)) << "--"<< G->id(out) <<"->"<< G->id(G->head(out)); 
+	  else 
+	    std::cout << "invalid"; 
+	}
+	else {
+	  std::cout << "in "; 
+	  if (in.valid()) 
+	    std::cout << G->id(G->head(in)) << "<-"<< G->id(in) <<"--"<< G->id(G->tail(in)); 
+	  else 
+	    std::cout << "invalid"; 
+	}
+	std::cout << std::endl;
+      }
     };
 
     class OutEdgeIt : public EdgeIt {
@@ -308,11 +326,13 @@
       	G=&_G;
 	flow=&_flow;
 	capacity=&_capacity;
-	out=/*resG->*/G->template first<OldOutEdgeIt>(v);
+	//out=/*resG->*/G->template first<OldOutEdgeIt>(v);
+	G->getFirst(out, v);
 	while( out.valid() && !(free()>0) ) { ++out; }
 	if (!out.valid()) {
 	  out_or_in=0;
-	  in=/*resG->*/G->template first<OldInEdgeIt>(v);
+	  //in=/*resG->*/G->template first<OldInEdgeIt>(v);
+	  G->getFirst(in, v);
 	  while( in.valid() && !(free()>0) ) { ++in; }
 	}
       }
@@ -324,7 +344,7 @@
 	  while( out.valid() && !(free()>0) ) { ++out; }
 	  if (!out.valid()) {
 	    out_or_in=0;
-	    in=/*resG->*/G->template first<OldInEdgeIt>(v);
+	    G->getFirst(in, v); //=/*resG->*/G->template first<OldInEdgeIt>(v);
 	    while( in.valid() && !(free()>0) ) { ++in; }
 	  }
 	} else {
@@ -335,9 +355,75 @@
       }
     };
 
+    class EachEdgeIt : public EdgeIt {
+      typename Graph::EachNodeIt v;
+    public:
+      EachEdgeIt() { }
+      //EachEdgeIt(const EachEdgeIt& e) : EdgeIt(e), v(e.v) { }
+      EachEdgeIt(const Graph& _G, FlowMap& _flow, const CapacityMap& _capacity) { 
+	G=&_G;
+	flow=&_flow;
+	capacity=&_capacity;
+	out_or_in=1;
+	G->getFirst(v);
+	if (v.valid()) G->getFirst(out, v); else out=OldOutEdgeIt();
+	while (out.valid() && !(free()>0) ) { ++out; }
+	while (v.valid() && !out.valid()) { 
+	  ++v; 
+	  if (v.valid()) G->getFirst(out, v); 
+	  while (out.valid() && !(free()>0) ) { ++out; }
+	}
+	if (!out.valid()) {
+	  out_or_in=0;
+	  G->getFirst(v);
+	  if (v.valid()) G->getFirst(in, v); else in=OldInEdgeIt();
+	  while (in.valid() && !(free()>0) ) { ++in; }
+	  while (v.valid() && !in.valid()) { 
+	    ++v; 
+	    if (v.valid()) G->getFirst(in, v); 
+	    while (in.valid() && !(free()>0) ) { ++in; }
+	  }
+	}
+      }
+      EachEdgeIt& operator++() { 
+	if (out_or_in) {
+	  ++out;
+	  while (out.valid() && !(free()>0) ) { ++out; }
+	  while (v.valid() && !out.valid()) { 
+	    ++v; 
+	    if (v.valid()) G->getFirst(out, v); 
+	    while (out.valid() && !(free()>0) ) { ++out; }
+	  }
+	  if (!out.valid()) {
+	    out_or_in=0;
+	    G->getFirst(v);
+	    if (v.valid()) G->getFirst(in, v); else in=OldInEdgeIt();
+	    while (in.valid() && !(free()>0) ) { ++in; }
+	    while (v.valid() && !in.valid()) { 
+	      ++v; 
+	      if (v.valid()) G->getFirst(in, v); 
+	      while (in.valid() && !(free()>0) ) { ++in; }
+	    }  
+	  }
+	} else {
+	  ++in;
+	  while (in.valid() && !(free()>0) ) { ++in; }
+	  while (v.valid() && !in.valid()) { 
+	    ++v; 
+	    if (v.valid()) G->getFirst(in, v); 
+	    while (in.valid() && !(free()>0) ) { ++in; }
+	  }
+	}
+	return *this;
+      }
+    };
+
     void getFirst(OutEdgeIt& e, NodeIt v) const { 
       e=OutEdgeIt(G, v, flow, capacity); 
     }
+    void getFirst(EachEdgeIt& e) const { 
+      e=EachEdgeIt(G, flow, capacity); 
+    }
     void getFirst(EachNodeIt& v) const { G.getFirst(v); }
     
     template< typename It >
@@ -401,13 +487,13 @@
     //AugGraph res_graph;    
     //typedef typename AugGraph::NodeMap<bool> ReachedMap;
     //typename AugGraph::NodeMap<AugEdgeIt> pred; 
-    //typename AugGraph::NodeMap<int> free;
+    //typename AugGraph::NodeMap<Number> 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) //,  
       //res_graph(G, flow, capacity), pred(res_graph), free(res_graph) 
       { }
-    bool augment() {
+    bool augmentOnShortestPath() {
       AugGraph res_graph(G, flow, capacity);
       bool _augment=false;
       
@@ -419,7 +505,7 @@
       //filled up with invalid iterators
       //pred.set(s, AugEdgeIt());
       
-      typename AugGraph::NodeMap<int> free(res_graph);
+      typename AugGraph::NodeMap<Number> free(res_graph);
 	
       //searching for augmenting path
       while ( !res_bfs.finished() ) { 
@@ -451,48 +537,130 @@
 
       return _augment;
     }
-    bool augmentWithBlockingFlow() {
-      BfsIterator4< Graph, OutEdgeIt, typename Graph::NodeMap<bool> > bfs(G);
+    template<typename MutableGraph> bool augmentOnBlockingFlow() {
+      bool _augment=false;
+
+      AugGraph res_graph(G, flow, capacity);
+
+      typedef typename AugGraph::NodeMap<bool> ReachedMap;
+      BfsIterator4< AugGraph, AugOutEdgeIt, ReachedMap > bfs(res_graph);
+
       bfs.pushAndSetReached(s);
-      typename Graph::NodeMap<int> dist(G); //filled up with 0's
+      typename AugGraph::NodeMap<int> dist(res_graph); //filled up with 0's
       while ( !bfs.finished() ) { 
-	OutEdgeIt e=OutEdgeIt(bfs);
+	AugOutEdgeIt e=AugOutEdgeIt(bfs);
 	if (e.valid() && bfs.isBNodeNewlyReached()) {
-	  dist.set(G.head(e), dist.get(G.tail(e))+1);
-	  //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 (res_graph.head(e)==t) { _augment=true; break; }
+	  dist.set(res_graph.head(e), dist.get(res_graph.tail(e))+1);
 	}
 	
 	++bfs;
-      } //end of searching augmenting path
+      } //computing distances from s in the residual graph
 
-      double pre_time_copy=currTime();
-      typedef Graph MutableGraph;
       MutableGraph F;
-      typename Graph::NodeMap<NodeIt> G_to_F(G);
-      for(typename Graph::EachNodeIt n=G.template first<typename Graph::EachNodeIt>(); n.valid(); ++n) {
-	G_to_F.set(n, F.addNode());
+      typename AugGraph::NodeMap<NodeIt> res_graph_to_F(res_graph);
+      for(typename AugGraph::EachNodeIt n=res_graph.template first<typename AugGraph::EachNodeIt>(); n.valid(); ++n) {
+	res_graph_to_F.set(n, F.addNode());
       }
-      for(typename Graph::EachEdgeIt e=G.template first<typename Graph::EachEdgeIt>(); e.valid(); ++e) {
-	if (dist.get(G.head(e))==dist.get(G.tail(e))+1) {
-	  F.addEdge(G_to_F.get(G.tail(e)), G_to_F.get(G.head(e)));
+      
+      typename MutableGraph::NodeIt sF=res_graph_to_F.get(s);
+      typename MutableGraph::NodeIt tF=res_graph_to_F.get(t);
+
+      typename MutableGraph::EdgeMap<AugEdgeIt> original_edge(F);
+      typename MutableGraph::EdgeMap<Number> free_on_edge(F);
+
+      //Making F to the graph containing the edges of the residual graph 
+      //which are in some shortest paths
+      for(typename AugGraph::EachEdgeIt e=res_graph.template first<typename AugGraph::EachEdgeIt>(); e.valid(); ++e) {
+	if (dist.get(res_graph.head(e))==dist.get(res_graph.tail(e))+1) {
+	  typename MutableGraph::EdgeIt 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);
+	  free_on_edge.update();
+	  free_on_edge.set(f, e.free());
+	} 
+      }
+
+      bool __augment=true;
+
+      while (__augment) {
+	__augment=false;
+	//computing blocking flow with dfs
+	typedef typename MutableGraph::NodeMap<bool> BlockingReachedMap;
+	DfsIterator4< MutableGraph, typename MutableGraph::OutEdgeIt, BlockingReachedMap > dfs(F);
+	typename MutableGraph::NodeMap<EdgeIt> pred(F); //invalid iterators
+	typename MutableGraph::NodeMap<Number> free(F);
+
+	dfs.pushAndSetReached(sF);      
+	while (!dfs.finished()) {
+	  ++dfs;
+	  if (typename MutableGraph::OutEdgeIt(dfs).valid()) {
+	    //std::cout << "OutEdgeIt: " << dfs; 
+	    //std::cout << " aNode: " << F.aNode(dfs); 
+	    //std::cout << " bNode: " << F.bNode(dfs) << " ";
+	  
+	    typename MutableGraph::NodeIt v=F.aNode(dfs);
+	    typename MutableGraph::NodeIt w=F.bNode(dfs);
+	    pred.set(w, dfs);
+	    if (pred.get(v).valid()) {
+	      free.set(w, std::min(free.get(v), free_on_edge.get(dfs)));
+	    } else {
+	      free.set(w, free_on_edge.get(dfs)/*original_edge.get(dfs).free()*/); 
+	    }
+	    if (w==tF) { 
+	      //std::cout << "AUGMENTATION"<<std::endl;
+	      __augment=true; 
+	      _augment=true;
+	      break; 
+	    }
+	  } else { 
+	    //std::cout << "OutEdgeIt: " << "invalid"; 
+	    //std::cout << " aNode: " << dfs.aNode(); 
+	    //std::cout << " bNode: " << "invalid" << " ";
+	  }
+	  if (dfs.isBNodeNewlyReached()) { 
+	    //std::cout << "bNodeIsNewlyReached ";
+	  } else { 
+	    //std::cout << "bNodeIsNotNewlyReached ";
+	    if (typename MutableGraph::OutEdgeIt(dfs).valid()) {
+	      //std::cout << "DELETE ";
+	      F.erase(typename MutableGraph::OutEdgeIt(dfs));
+	    }
+	  } 
+	  //if (dfs.isANodeExamined()) { 
+	    //std::cout << "aNodeIsExamined ";
+	    //} else { 
+	    //std::cout << "aNodeIsNotExamined ";
+	    //} 
+	  //std::cout<<std::endl;
 	}
+
+	if (__augment) {
+	  typename MutableGraph::NodeIt n=tF;
+	  Number augment_value=free.get(tF);
+	  while (pred.get(n).valid()) { 
+	    typename MutableGraph::EdgeIt e=pred.get(n);
+	    original_edge.get(e).augment(augment_value); 
+	    n=F.tail(e);
+	    F.erase(e);
+	  }
+	}
+
       }
-      double post_time_copy=currTime();
-      std::cout << "copy time: " << post_time_copy-pre_time_copy << " sec"<< std::endl; 
-
-      return 0;
+            
+      return _augment;
     }
     void run() {
       //int num_of_augmentations=0;
-      while (augment()) { 
+      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;
       } 
@@ -599,7 +767,7 @@
       typename AugGraph::NodeMap<AugEdgeIt> pred(res_graph); 
       //filled up with invalid iterators
       
-      typename AugGraph::NodeMap<int> free(res_graph);
+      typename AugGraph::NodeMap<Number> free(res_graph);
 	
       //searching for augmenting path
       while ( !res_bfs.finished() ) { 
diff -r 1ac27e476e25 -r 0631992fe7a1 src/work/marci/edmonds_karp_demo.cc
--- a/src/work/marci/edmonds_karp_demo.cc	Thu Feb 26 16:07:40 2004 +0000
+++ b/src/work/marci/edmonds_karp_demo.cc	Fri Feb 27 12:39:15 2004 +0000
@@ -19,29 +19,16 @@
   ListGraph::EdgeMap<int> cap(G);
   readDimacsMaxFlow(std::cin, G, s, t, cap);
 
-/*
-  double pre_time_copy=currTime();
-  ListGraph F;
-  ListGraph::NodeMap<NodeIt> G_to_F(G);
-  typedef ListGraph::EachNodeIt EachNodeIt;
-  for(EachNodeIt n=G.first<EachNodeIt>(); n.valid(); ++n) {
-    G_to_F.set(n, F.addNode());
-  }
-  for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) {
-    F.addEdge(G_to_F.get(G.tail(e)), G_to_F.get(G.head(e)));
-  }
-  double post_time_copy=currTime();
-  std::cout << "copy time: " << post_time_copy-pre_time_copy << " sec"<< std::endl; 
-*/
-
-  std::cout << "edmonds karp demo..." << std::endl;
+  {
+  std::cout << "edmonds karp demo (blocking flow augmentation)..." << std::endl;
   ListGraph::EdgeMap<int> flow(G); //0 flow
 
   double pre_time=currTime();
   MaxFlow<ListGraph, int, ListGraph::EdgeMap<int>, ListGraph::EdgeMap<int> > max_flow_test(G, s, t, flow, cap);
-  max_flow_test.augmentWithBlockingFlow();
-  max_flow_test.run();
+  //max_flow_test.augmentWithBlockingFlow<ListGraph>();
+  max_flow_test.run<ListGraph>();
   double post_time=currTime();
+
   //std::cout << "maximum flow: "<< std::endl;
   //for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) { 
   //  std::cout<<"("<<G.tail(e)<< "-"<<flow.get(e)<<"->"<<G.head(e)<<") ";
@@ -49,6 +36,26 @@
   //std::cout<<std::endl;
   std::cout << "elapsed time: " << post_time-pre_time << " sec"<< std::endl; 
   std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
+  }
+
+  {
+  std::cout << "edmonds karp demo (shortest path augmentation)..." << std::endl;
+  ListGraph::EdgeMap<int> flow(G); //0 flow
+
+  double pre_time=currTime();
+  MaxFlow<ListGraph, int, ListGraph::EdgeMap<int>, ListGraph::EdgeMap<int> > max_flow_test(G, s, t, flow, cap);
+  //max_flow_test.augmentWithBlockingFlow<ListGraph>();
+  max_flow_test.run();
+  double post_time=currTime();
+
+  //std::cout << "maximum flow: "<< std::endl;
+  //for(EachEdgeIt e=G.first<EachEdgeIt>(); e.valid(); ++e) { 
+  //  std::cout<<"("<<G.tail(e)<< "-"<<flow.get(e)<<"->"<<G.head(e)<<") ";
+  //}
+  //std::cout<<std::endl;
+  std::cout << "elapsed time: " << post_time-pre_time << " sec"<< std::endl; 
+  std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
+  }
 
   return 0;
 }
diff -r 1ac27e476e25 -r 0631992fe7a1 src/work/marci_graph_demo.cc
--- a/src/work/marci_graph_demo.cc	Thu Feb 26 16:07:40 2004 +0000
+++ b/src/work/marci_graph_demo.cc	Fri Feb 27 12:39:15 2004 +0000
@@ -225,6 +225,17 @@
   {
     ListGraph::EdgeMap<int> flow(flowG, 0);
     MaxFlow<ListGraph, int, ListGraph::EdgeMap<int>, ListGraph::EdgeMap<int> > max_flow_test(flowG, s, t, flow, cap);
+    /*
+    max_flow_test.augmentOnBlockingFlow<ListGraph>();
+    for(EachEdgeIt e=flowG.template first<EachEdgeIt>(); e.valid(); ++e) { 
+      std::cout<<"("<<flowG.tail(e)<< "-"<<flow.get(e)<<"->"<<flowG.head(e)<<") ";
+    }
+    std::cout<<std::endl;
+    max_flow_test.augmentOnBlockingFlow<ListGraph>();
+    for(EachEdgeIt e=flowG.template first<EachEdgeIt>(); e.valid(); ++e) { 
+      std::cout<<"("<<flowG.tail(e)<< "-"<<flow.get(e)<<"->"<<flowG.head(e)<<") ";
+    }
+    std::cout<<std::endl;*/
     max_flow_test.run();
     
     std::cout << "maximum flow: "<< std::endl;