# HG changeset patch
# User marci
# Date 1083254314 0
# Node ID 052af4060f3ebab3621c4add3b4f128ecc0a3f63
# Parent  a40985a922d038d8863581ef19ed48ecad09317f
preflow, maxflow

diff -r a40985a922d0 -r 052af4060f3e src/work/jacint/preflow.h
--- a/src/work/jacint/preflow.h	Thu Apr 29 15:01:52 2004 +0000
+++ b/src/work/jacint/preflow.h	Thu Apr 29 15:58:34 2004 +0000
@@ -44,6 +44,11 @@
 #include <stack>
 
 #include <graph_wrapper.h>
+#include <bfs_iterator.h>
+#include <invalid.h>
+#include <maps.h>
+#include <for_each_macros.h>
+
 
 namespace hugo {
 
@@ -111,9 +116,14 @@
 
     bool augmentOnBlockingFlow2();
 
-    //Returns the maximum value of a flow.
-    Num flowValue() {
-      return excess[t];
+    /// Returns the actual flow value.
+    /// More precisely, it returns the negative excess of s, thus 
+    /// this works also for preflows.
+    Num flowValue() { 
+      Num a=0;
+      FOR_EACH_INC_LOC(OutEdgeIt, e, *g, s) a+=(*flow)[e];
+      FOR_EACH_INC_LOC(InEdgeIt, e, *g, s) a-=(*flow)[e];
+      return a;
     }
 
     //should be used only between preflowPhase0 and preflowPhase1
@@ -433,7 +443,8 @@
 
     void relabel(Node w, int newlevel, VecStack& active,  
 		 VecNode& level_list, NNMap& left, 
-		 NNMap& right, int& b, int& k, bool what_heur ) {
+		 NNMap& right, int& b, int& k, bool what_heur ) 
+    {
 
       Num lev=level[w];	
       
@@ -497,6 +508,28 @@
       }
       
     } //relabel
+
+
+    template<typename MapGraphWrapper> 
+    class DistanceMap {
+    protected:
+      const MapGraphWrapper* g;
+      typename MapGraphWrapper::template NodeMap<int> dist; 
+    public:
+      DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
+      void set(const typename MapGraphWrapper::Node& n, int a) { 
+	dist.set(n, a); 
+      }
+      int operator[](const typename MapGraphWrapper::Node& n) 
+	{ return dist[n]; }
+//       int get(const typename MapGraphWrapper::Node& n) const { 
+// 	return dist[n]; }
+//       bool get(const typename MapGraphWrapper::Edge& e) const { 
+// 	return (dist.get(g->tail(e))<dist.get(g->head(e))); }
+      bool operator[](const typename MapGraphWrapper::Edge& e) const { 
+	return (dist[g->tail(e)]<dist[g->head(e)]); 
+      }
+    };
     
   };
 
@@ -674,8 +707,52 @@
 
 
 
+  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
+  bool Preflow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath() 
+  {
+      ResGW res_graph(*g, *capacity, *flow);
+      bool _augment=false;
+      
+      //ReachedMap level(res_graph);
+      FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
+      BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
+      bfs.pushAndSetReached(s);
+	
+      typename ResGW::template NodeMap<ResGWEdge> pred(res_graph); 
+      pred.set(s, INVALID);
+      
+      typename ResGW::template NodeMap<Num> free(res_graph);
+	
+      //searching for augmenting path
+      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[v])) {
+	    free.set(w, std::min(free[v], res_graph.resCap(e)));
+	  } else {
+	    free.set(w, res_graph.resCap(e)); 
+	  }
+	  if (res_graph.head(e)==t) { _augment=true; break; }
+	}
+	
+	++bfs;
+      } //end of searching augmenting path
 
+      if (_augment) {
+	Node n=t;
+	Num augment_value=free[t];
+	while (res_graph.valid(pred[n])) { 
+	  ResGWEdge e=pred[n];
+	  res_graph.augment(e, augment_value); 
+	  n=res_graph.tail(e);
+	}
+      }
 
+      return _augment;
+    }
 
 
 
@@ -683,6 +760,253 @@
 
 
 
+
+
+  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
+  template<typename MutableGraph> 
+  bool Preflow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow() 
+  {      
+      typedef MutableGraph MG;
+      bool _augment=false;
+
+      ResGW res_graph(*g, *capacity, *flow);
+
+      //bfs for distances on the residual graph
+      //ReachedMap level(res_graph);
+      FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
+      BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
+      bfs.pushAndSetReached(s);
+      typename ResGW::template 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 are on shortest paths
+      MG F;
+      typename ResGW::template NodeMap<typename MG::Node> 
+	res_graph_to_F(res_graph);
+      {
+	typename ResGW::NodeIt n;
+	for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
+	  res_graph_to_F.set(n, F.addNode());
+	}
+      }
+
+      typename MG::Node sF=res_graph_to_F[s];
+      typename MG::Node tF=res_graph_to_F[t];
+      typename MG::template EdgeMap<ResGWEdge> original_edge(F);
+      typename MG::template EdgeMap<Num> residual_capacity(F);
+
+      while ( !bfs.finished() ) { 
+	ResGWOutEdgeIt e=bfs;
+	if (res_graph.valid(e)) {
+	  if (bfs.isBNodeNewlyReached()) {
+	    dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
+	    typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
+	    original_edge.update();
+	    original_edge.set(f, e);
+	    residual_capacity.update();
+	    residual_capacity.set(f, res_graph.resCap(e));
+	  } else {
+	    if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
+	      typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
+	      original_edge.update();
+	      original_edge.set(f, e);
+	      residual_capacity.update();
+	      residual_capacity.set(f, res_graph.resCap(e));
+	    }
+	  }
+	}
+	++bfs;
+      } //computing distances from s in the residual graph
+
+      bool __augment=true;
+
+      while (__augment) {
+	__augment=false;
+	//computing blocking flow with dfs
+	DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
+	typename MG::template NodeMap<typename MG::Edge> pred(F);
+	pred.set(sF, INVALID);
+	//invalid iterators for sources
+
+	typename MG::template NodeMap<Num> free(F);
+
+	dfs.pushAndSetReached(sF);      
+	while (!dfs.finished()) {
+	  ++dfs;
+	  if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
+	    if (dfs.isBNodeNewlyReached()) {
+	      typename MG::Node v=F.aNode(dfs);
+	      typename MG::Node w=F.bNode(dfs);
+	      pred.set(w, dfs);
+	      if (F.valid(pred[v])) {
+		free.set(w, std::min(free[v], residual_capacity[dfs]));
+	      } else {
+		free.set(w, residual_capacity[dfs]); 
+	      }
+	      if (w==tF) { 
+		__augment=true; 
+		_augment=true;
+		break; 
+	      }
+	      
+	    } else {
+	      F.erase(/*typename MG::OutEdgeIt*/(dfs));
+	    }
+	  } 
+	}
+
+	if (__augment) {
+	  typename MG::Node n=tF;
+	  Num augment_value=free[tF];
+	  while (F.valid(pred[n])) { 
+	    typename MG::Edge e=pred[n];
+	    res_graph.augment(original_edge[e], augment_value); 
+	    n=F.tail(e);
+	    if (residual_capacity[e]==augment_value) 
+	      F.erase(e); 
+	    else 
+	      residual_capacity.set(e, residual_capacity[e]-augment_value);
+	  }
+	}
+	
+      }
+            
+      return _augment;
+    }
+
+
+
+
+
+
+  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
+  bool Preflow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2() 
+  {
+      bool _augment=false;
+
+      ResGW res_graph(*g, *capacity, *flow);
+      
+      //ReachedMap level(res_graph);
+      FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
+      BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
+
+      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[res_graph.tail(e)]+1);
+ 	}
+	++bfs;
+      } //computing distances from s in the residual graph
+
+      //Subgraph containing the edges on some shortest paths
+      ConstMap<typename ResGW::Node, bool> true_map(true);
+      typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>, 
+	DistanceMap<ResGW> > FilterResGW;
+      FilterResGW filter_res_graph(res_graph, true_map, dist);
+
+      //Subgraph, which is able to delete edges which are already 
+      //met by the dfs
+      typename FilterResGW::template 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::
+	template 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
+  	DfsIterator< ErasingResGW, 
+	  typename ErasingResGW::template NodeMap<bool> > 
+  	  dfs(erasing_res_graph);
+ 	typename ErasingResGW::
+	  template NodeMap<typename ErasingResGW::OutEdgeIt> 
+	  pred(erasing_res_graph); 
+ 	pred.set(s, INVALID);
+  	//invalid iterators for sources
+
+  	typename ErasingResGW::template NodeMap<Num> 
+	  free1(erasing_res_graph);
+
+ 	dfs.pushAndSetReached(
+	  typename ErasingResGW::Node(
+	    typename FilterResGW::Node(
+	      typename ResGW::Node(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[v])) {
+ 		free1.set(w, std::min(free1[v], res_graph.resCap(
+				       typename ErasingResGW::OutEdgeIt(dfs))));
+ 	      } else {
+ 		free1.set(w, res_graph.resCap(
+			   typename ErasingResGW::OutEdgeIt(dfs))); 
+ 	      }
+	      
+ 	      if (w==t) { 
+ 		__augment=true; 
+ 		_augment=true;
+ 		break; 
+ 	      }
+ 	    } else {
+ 	      erasing_res_graph.erase(dfs);
+	    }
+	  }
+	}	
+
+  	if (__augment) {
+   	  typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
+// 	  typename ResGW::NodeMap<Num> a(res_graph);
+// 	  typename ResGW::Node b;
+// 	  Num j=a[b];
+// 	  typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
+// 	  typename FilterResGW::Node b1;
+// 	  Num j1=a1[b1];
+// 	  typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
+// 	  typename ErasingResGW::Node b2;
+// 	  Num j2=a2[b2];
+ 	  Num augment_value=free1[n];
+ 	  while (erasing_res_graph.valid(pred[n])) { 
+ 	    typename ErasingResGW::OutEdgeIt e=pred[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;
+    }
+
+
+
+
 } //namespace hugo
 
 #endif //HUGO_PREFLOW_H
diff -r a40985a922d0 -r 052af4060f3e src/work/marci/edmonds_karp_demo.cc
--- a/src/work/marci/edmonds_karp_demo.cc	Thu Apr 29 15:01:52 2004 +0000
+++ b/src/work/marci/edmonds_karp_demo.cc	Thu Apr 29 15:58:34 2004 +0000
@@ -5,11 +5,11 @@
 #include <list_graph.h>
 #include <smart_graph.h>
 #include <dimacs.h>
-#include <edmonds_karp.h>
+//#include <edmonds_karp.h>
 #include <time_measure.h>
 //#include <graph_wrapper.h>
 #include <preflow.h>
-#include <preflow_res.h>
+//#include <preflow_res.h>
 #include <for_each_macros.h>
 
 using namespace hugo;
@@ -72,12 +72,12 @@
   Graph::EdgeMap<int> flow(G); //0 flow
   Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
     pre_flow_test(G, s, t, cap, flow/*, true*/);
-  Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
-    pre_flow_ize(G, s, t, cap, flow/*, false*/);
+  //  Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
+  //  pre_flow_ize(G, s, t, cap, flow/*, false*/);
 //   PreflowRes<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
 //     pre_flow_res(G, s, t, cap, flow/*, true*/);
-  MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
-    max_flow_test(G, s, t, cap, flow);
+  //MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
+  //  max_flow_test(G, s, t, cap, flow);
 
   {
     std::cout << "preflow ..." << std::endl;
@@ -91,9 +91,9 @@
     std::cout << "preflow ..." << std::endl;
     FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
     ts.reset();
-    pre_flow_ize.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
+    pre_flow_test.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
     std::cout << "elapsed time: " << ts << std::endl;
-    std::cout << "flow value: "<< pre_flow_ize.flowValue() << std::endl;
+    std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
   }
 
 //   {
@@ -110,32 +110,32 @@
     FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
     ts.reset();
     int i=0;
-    while (max_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
+    while (pre_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
     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;
+    std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
   }
 
-  {
-    std::cout << "faster physical blocking flow augmentation ..." << std::endl;
-    FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
-    ts.reset();
-    int i=0;
-    while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
-    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;
-  }
+//   {
+//     std::cout << "faster physical blocking flow augmentation ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
+//     ts.reset();
+//     int i=0;
+//     while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
+//     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;
+//   }
 
   {
     std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
     FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
     ts.reset();
     int i=0;
-    while (max_flow_test.augmentOnBlockingFlow2()) { ++i; }
+    while (pre_flow_test.augmentOnBlockingFlow2()) { ++i; }
     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;
+    std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
   }
 
   {
@@ -143,10 +143,10 @@
     FOR_EACH_LOC(Graph::EdgeIt, e, G) flow.set(e, 0);
     ts.reset();
     int i=0;
-    while (max_flow_test.augmentOnShortestPath()) { ++i; }
+    while (pre_flow_test.augmentOnShortestPath()) { ++i; }
     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;
+    std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
   }