src/work/jacint/preflow.h
changeset 472 052af4060f3e
parent 471 a40985a922d0
child 476 cfe550761745
     1.1 --- a/src/work/jacint/preflow.h	Thu Apr 29 15:01:52 2004 +0000
     1.2 +++ b/src/work/jacint/preflow.h	Thu Apr 29 15:58:34 2004 +0000
     1.3 @@ -44,6 +44,11 @@
     1.4  #include <stack>
     1.5  
     1.6  #include <graph_wrapper.h>
     1.7 +#include <bfs_iterator.h>
     1.8 +#include <invalid.h>
     1.9 +#include <maps.h>
    1.10 +#include <for_each_macros.h>
    1.11 +
    1.12  
    1.13  namespace hugo {
    1.14  
    1.15 @@ -111,9 +116,14 @@
    1.16  
    1.17      bool augmentOnBlockingFlow2();
    1.18  
    1.19 -    //Returns the maximum value of a flow.
    1.20 -    Num flowValue() {
    1.21 -      return excess[t];
    1.22 +    /// Returns the actual flow value.
    1.23 +    /// More precisely, it returns the negative excess of s, thus 
    1.24 +    /// this works also for preflows.
    1.25 +    Num flowValue() { 
    1.26 +      Num a=0;
    1.27 +      FOR_EACH_INC_LOC(OutEdgeIt, e, *g, s) a+=(*flow)[e];
    1.28 +      FOR_EACH_INC_LOC(InEdgeIt, e, *g, s) a-=(*flow)[e];
    1.29 +      return a;
    1.30      }
    1.31  
    1.32      //should be used only between preflowPhase0 and preflowPhase1
    1.33 @@ -433,7 +443,8 @@
    1.34  
    1.35      void relabel(Node w, int newlevel, VecStack& active,  
    1.36  		 VecNode& level_list, NNMap& left, 
    1.37 -		 NNMap& right, int& b, int& k, bool what_heur ) {
    1.38 +		 NNMap& right, int& b, int& k, bool what_heur ) 
    1.39 +    {
    1.40  
    1.41        Num lev=level[w];	
    1.42        
    1.43 @@ -497,6 +508,28 @@
    1.44        }
    1.45        
    1.46      } //relabel
    1.47 +
    1.48 +
    1.49 +    template<typename MapGraphWrapper> 
    1.50 +    class DistanceMap {
    1.51 +    protected:
    1.52 +      const MapGraphWrapper* g;
    1.53 +      typename MapGraphWrapper::template NodeMap<int> dist; 
    1.54 +    public:
    1.55 +      DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
    1.56 +      void set(const typename MapGraphWrapper::Node& n, int a) { 
    1.57 +	dist.set(n, a); 
    1.58 +      }
    1.59 +      int operator[](const typename MapGraphWrapper::Node& n) 
    1.60 +	{ return dist[n]; }
    1.61 +//       int get(const typename MapGraphWrapper::Node& n) const { 
    1.62 +// 	return dist[n]; }
    1.63 +//       bool get(const typename MapGraphWrapper::Edge& e) const { 
    1.64 +// 	return (dist.get(g->tail(e))<dist.get(g->head(e))); }
    1.65 +      bool operator[](const typename MapGraphWrapper::Edge& e) const { 
    1.66 +	return (dist[g->tail(e)]<dist[g->head(e)]); 
    1.67 +      }
    1.68 +    };
    1.69      
    1.70    };
    1.71  
    1.72 @@ -674,8 +707,52 @@
    1.73  
    1.74  
    1.75  
    1.76 +  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
    1.77 +  bool Preflow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath() 
    1.78 +  {
    1.79 +      ResGW res_graph(*g, *capacity, *flow);
    1.80 +      bool _augment=false;
    1.81 +      
    1.82 +      //ReachedMap level(res_graph);
    1.83 +      FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
    1.84 +      BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
    1.85 +      bfs.pushAndSetReached(s);
    1.86 +	
    1.87 +      typename ResGW::template NodeMap<ResGWEdge> pred(res_graph); 
    1.88 +      pred.set(s, INVALID);
    1.89 +      
    1.90 +      typename ResGW::template NodeMap<Num> free(res_graph);
    1.91 +	
    1.92 +      //searching for augmenting path
    1.93 +      while ( !bfs.finished() ) { 
    1.94 +	ResGWOutEdgeIt e=bfs;
    1.95 +	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
    1.96 +	  Node v=res_graph.tail(e);
    1.97 +	  Node w=res_graph.head(e);
    1.98 +	  pred.set(w, e);
    1.99 +	  if (res_graph.valid(pred[v])) {
   1.100 +	    free.set(w, std::min(free[v], res_graph.resCap(e)));
   1.101 +	  } else {
   1.102 +	    free.set(w, res_graph.resCap(e)); 
   1.103 +	  }
   1.104 +	  if (res_graph.head(e)==t) { _augment=true; break; }
   1.105 +	}
   1.106 +	
   1.107 +	++bfs;
   1.108 +      } //end of searching augmenting path
   1.109  
   1.110 +      if (_augment) {
   1.111 +	Node n=t;
   1.112 +	Num augment_value=free[t];
   1.113 +	while (res_graph.valid(pred[n])) { 
   1.114 +	  ResGWEdge e=pred[n];
   1.115 +	  res_graph.augment(e, augment_value); 
   1.116 +	  n=res_graph.tail(e);
   1.117 +	}
   1.118 +      }
   1.119  
   1.120 +      return _augment;
   1.121 +    }
   1.122  
   1.123  
   1.124  
   1.125 @@ -683,6 +760,253 @@
   1.126  
   1.127  
   1.128  
   1.129 +
   1.130 +
   1.131 +  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   1.132 +  template<typename MutableGraph> 
   1.133 +  bool Preflow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow() 
   1.134 +  {      
   1.135 +      typedef MutableGraph MG;
   1.136 +      bool _augment=false;
   1.137 +
   1.138 +      ResGW res_graph(*g, *capacity, *flow);
   1.139 +
   1.140 +      //bfs for distances on the residual graph
   1.141 +      //ReachedMap level(res_graph);
   1.142 +      FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   1.143 +      BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   1.144 +      bfs.pushAndSetReached(s);
   1.145 +      typename ResGW::template NodeMap<int> 
   1.146 +	dist(res_graph); //filled up with 0's
   1.147 +
   1.148 +      //F will contain the physical copy of the residual graph
   1.149 +      //with the set of edges which are on shortest paths
   1.150 +      MG F;
   1.151 +      typename ResGW::template NodeMap<typename MG::Node> 
   1.152 +	res_graph_to_F(res_graph);
   1.153 +      {
   1.154 +	typename ResGW::NodeIt n;
   1.155 +	for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
   1.156 +	  res_graph_to_F.set(n, F.addNode());
   1.157 +	}
   1.158 +      }
   1.159 +
   1.160 +      typename MG::Node sF=res_graph_to_F[s];
   1.161 +      typename MG::Node tF=res_graph_to_F[t];
   1.162 +      typename MG::template EdgeMap<ResGWEdge> original_edge(F);
   1.163 +      typename MG::template EdgeMap<Num> residual_capacity(F);
   1.164 +
   1.165 +      while ( !bfs.finished() ) { 
   1.166 +	ResGWOutEdgeIt e=bfs;
   1.167 +	if (res_graph.valid(e)) {
   1.168 +	  if (bfs.isBNodeNewlyReached()) {
   1.169 +	    dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
   1.170 +	    typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
   1.171 +	    original_edge.update();
   1.172 +	    original_edge.set(f, e);
   1.173 +	    residual_capacity.update();
   1.174 +	    residual_capacity.set(f, res_graph.resCap(e));
   1.175 +	  } else {
   1.176 +	    if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
   1.177 +	      typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
   1.178 +	      original_edge.update();
   1.179 +	      original_edge.set(f, e);
   1.180 +	      residual_capacity.update();
   1.181 +	      residual_capacity.set(f, res_graph.resCap(e));
   1.182 +	    }
   1.183 +	  }
   1.184 +	}
   1.185 +	++bfs;
   1.186 +      } //computing distances from s in the residual graph
   1.187 +
   1.188 +      bool __augment=true;
   1.189 +
   1.190 +      while (__augment) {
   1.191 +	__augment=false;
   1.192 +	//computing blocking flow with dfs
   1.193 +	DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
   1.194 +	typename MG::template NodeMap<typename MG::Edge> pred(F);
   1.195 +	pred.set(sF, INVALID);
   1.196 +	//invalid iterators for sources
   1.197 +
   1.198 +	typename MG::template NodeMap<Num> free(F);
   1.199 +
   1.200 +	dfs.pushAndSetReached(sF);      
   1.201 +	while (!dfs.finished()) {
   1.202 +	  ++dfs;
   1.203 +	  if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
   1.204 +	    if (dfs.isBNodeNewlyReached()) {
   1.205 +	      typename MG::Node v=F.aNode(dfs);
   1.206 +	      typename MG::Node w=F.bNode(dfs);
   1.207 +	      pred.set(w, dfs);
   1.208 +	      if (F.valid(pred[v])) {
   1.209 +		free.set(w, std::min(free[v], residual_capacity[dfs]));
   1.210 +	      } else {
   1.211 +		free.set(w, residual_capacity[dfs]); 
   1.212 +	      }
   1.213 +	      if (w==tF) { 
   1.214 +		__augment=true; 
   1.215 +		_augment=true;
   1.216 +		break; 
   1.217 +	      }
   1.218 +	      
   1.219 +	    } else {
   1.220 +	      F.erase(/*typename MG::OutEdgeIt*/(dfs));
   1.221 +	    }
   1.222 +	  } 
   1.223 +	}
   1.224 +
   1.225 +	if (__augment) {
   1.226 +	  typename MG::Node n=tF;
   1.227 +	  Num augment_value=free[tF];
   1.228 +	  while (F.valid(pred[n])) { 
   1.229 +	    typename MG::Edge e=pred[n];
   1.230 +	    res_graph.augment(original_edge[e], augment_value); 
   1.231 +	    n=F.tail(e);
   1.232 +	    if (residual_capacity[e]==augment_value) 
   1.233 +	      F.erase(e); 
   1.234 +	    else 
   1.235 +	      residual_capacity.set(e, residual_capacity[e]-augment_value);
   1.236 +	  }
   1.237 +	}
   1.238 +	
   1.239 +      }
   1.240 +            
   1.241 +      return _augment;
   1.242 +    }
   1.243 +
   1.244 +
   1.245 +
   1.246 +
   1.247 +
   1.248 +
   1.249 +  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
   1.250 +  bool Preflow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2() 
   1.251 +  {
   1.252 +      bool _augment=false;
   1.253 +
   1.254 +      ResGW res_graph(*g, *capacity, *flow);
   1.255 +      
   1.256 +      //ReachedMap level(res_graph);
   1.257 +      FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
   1.258 +      BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
   1.259 +
   1.260 +      bfs.pushAndSetReached(s);
   1.261 +      DistanceMap<ResGW> dist(res_graph);
   1.262 +      while ( !bfs.finished() ) { 
   1.263 + 	ResGWOutEdgeIt e=bfs;
   1.264 + 	if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
   1.265 + 	  dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
   1.266 + 	}
   1.267 +	++bfs;
   1.268 +      } //computing distances from s in the residual graph
   1.269 +
   1.270 +      //Subgraph containing the edges on some shortest paths
   1.271 +      ConstMap<typename ResGW::Node, bool> true_map(true);
   1.272 +      typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>, 
   1.273 +	DistanceMap<ResGW> > FilterResGW;
   1.274 +      FilterResGW filter_res_graph(res_graph, true_map, dist);
   1.275 +
   1.276 +      //Subgraph, which is able to delete edges which are already 
   1.277 +      //met by the dfs
   1.278 +      typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt> 
   1.279 + 	first_out_edges(filter_res_graph);
   1.280 +      typename FilterResGW::NodeIt v;
   1.281 +      for(filter_res_graph.first(v); filter_res_graph.valid(v); 
   1.282 + 	  filter_res_graph.next(v)) 
   1.283 +      {
   1.284 + 	typename FilterResGW::OutEdgeIt e;
   1.285 + 	filter_res_graph.first(e, v);
   1.286 + 	first_out_edges.set(v, e);
   1.287 +      }
   1.288 +      typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
   1.289 +	template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
   1.290 +      ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
   1.291 +
   1.292 +      bool __augment=true;
   1.293 +
   1.294 +      while (__augment) {
   1.295 +
   1.296 + 	__augment=false;
   1.297 +  	//computing blocking flow with dfs
   1.298 +  	DfsIterator< ErasingResGW, 
   1.299 +	  typename ErasingResGW::template NodeMap<bool> > 
   1.300 +  	  dfs(erasing_res_graph);
   1.301 + 	typename ErasingResGW::
   1.302 +	  template NodeMap<typename ErasingResGW::OutEdgeIt> 
   1.303 +	  pred(erasing_res_graph); 
   1.304 + 	pred.set(s, INVALID);
   1.305 +  	//invalid iterators for sources
   1.306 +
   1.307 +  	typename ErasingResGW::template NodeMap<Num> 
   1.308 +	  free1(erasing_res_graph);
   1.309 +
   1.310 + 	dfs.pushAndSetReached(
   1.311 +	  typename ErasingResGW::Node(
   1.312 +	    typename FilterResGW::Node(
   1.313 +	      typename ResGW::Node(s)
   1.314 +	      )
   1.315 +	    )
   1.316 +	  );
   1.317 + 	while (!dfs.finished()) {
   1.318 + 	  ++dfs;
   1.319 + 	  if (erasing_res_graph.valid(
   1.320 + 		typename ErasingResGW::OutEdgeIt(dfs))) 
   1.321 + 	  { 
   1.322 +  	    if (dfs.isBNodeNewlyReached()) {
   1.323 +	  
   1.324 + 	      typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
   1.325 + 	      typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
   1.326 +
   1.327 + 	      pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
   1.328 + 	      if (erasing_res_graph.valid(pred[v])) {
   1.329 + 		free1.set(w, std::min(free1[v], res_graph.resCap(
   1.330 +				       typename ErasingResGW::OutEdgeIt(dfs))));
   1.331 + 	      } else {
   1.332 + 		free1.set(w, res_graph.resCap(
   1.333 +			   typename ErasingResGW::OutEdgeIt(dfs))); 
   1.334 + 	      }
   1.335 +	      
   1.336 + 	      if (w==t) { 
   1.337 + 		__augment=true; 
   1.338 + 		_augment=true;
   1.339 + 		break; 
   1.340 + 	      }
   1.341 + 	    } else {
   1.342 + 	      erasing_res_graph.erase(dfs);
   1.343 +	    }
   1.344 +	  }
   1.345 +	}	
   1.346 +
   1.347 +  	if (__augment) {
   1.348 +   	  typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
   1.349 +// 	  typename ResGW::NodeMap<Num> a(res_graph);
   1.350 +// 	  typename ResGW::Node b;
   1.351 +// 	  Num j=a[b];
   1.352 +// 	  typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
   1.353 +// 	  typename FilterResGW::Node b1;
   1.354 +// 	  Num j1=a1[b1];
   1.355 +// 	  typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
   1.356 +// 	  typename ErasingResGW::Node b2;
   1.357 +// 	  Num j2=a2[b2];
   1.358 + 	  Num augment_value=free1[n];
   1.359 + 	  while (erasing_res_graph.valid(pred[n])) { 
   1.360 + 	    typename ErasingResGW::OutEdgeIt e=pred[n];
   1.361 + 	    res_graph.augment(e, augment_value);
   1.362 + 	    n=erasing_res_graph.tail(e);
   1.363 + 	    if (res_graph.resCap(e)==0)
   1.364 + 	      erasing_res_graph.erase(e);
   1.365 +	}
   1.366 +      }
   1.367 +      
   1.368 +      } //while (__augment) 
   1.369 +            
   1.370 +      return _augment;
   1.371 +    }
   1.372 +
   1.373 +
   1.374 +
   1.375 +
   1.376  } //namespace hugo
   1.377  
   1.378  #endif //HUGO_PREFLOW_H