lemon-0.x: comparison src/work/jacint/preflow.h

equal deleted inserted replaced

-:a8c0ea17699b
+:03a22be75ac4
 #include <vector>
 #include <queue>
 #include <stack>
 #include <graph_wrapper.h>
+#include <bfs_iterator.h>
+#include <invalid.h>
+#include <maps.h>
+#include <for_each_macros.h>
 namespace hugo {
 template <typename Graph, typename Num,
 	    typename CapMap=typename Graph::template EdgeMap<Num>,
 template<typename MutableGraph> bool augmentOnBlockingFlow();
 bool augmentOnBlockingFlow2();
-//Returns the maximum value of a flow.
+/// Returns the actual flow value.
-Num flowValue() {
+/// More precisely, it returns the negative excess of s, thus
-return excess[t];
+/// 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
 template<typename _CutMap>
 void actMinCut(_CutMap& M) {
 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];
 Node right_n=right[w];
 Node left_n=left[w];
 	  level_list[newlevel]=w;
 	}
 }
 } //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)]);
+}
+};
 };
 template <typename Graph, typename Num, typename CapMap, typename FlowMap>
 } // while(true)
 }
+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;
+}
+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

changeset 473	2cef25dcde3f
parent 471	a40985a922d0
child 476	cfe550761745