1.1 --- a/src/work/jacint/preflow_hl0.h Mon Mar 01 17:24:34 2004 +0000
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,508 +0,0 @@
1.4 -// -*- C++ -*-
1.5 -/*
1.6 -preflow_hl0.h
1.7 -by jacint.
1.8 -Heuristics:
1.9 - 2 phase
1.10 - gap
1.11 - list 'level_list' on the nodes on level i implemented by hand
1.12 - stack 'active' on the active nodes on level i implemented by hand
1.13 - bound decrease
1.14 -
1.15 -The bound decrease heuristic behaves unexpectedly well.
1.16 -
1.17 -The constructor runs the algorithm.
1.18 -
1.19 -Members:
1.20 -
1.21 -T maxFlow() : returns the value of a maximum flow
1.22 -
1.23 -T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
1.24 -
1.25 -FlowMap Flow() : returns the fixed maximum flow x
1.26 -
1.27 -void minMinCut(CutMap& M) : sets M to the characteristic vector of the
1.28 - minimum min cut. M should be a map of bools initialized to false.
1.29 -
1.30 -void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
1.31 - maximum min cut. M should be a map of bools initialized to false.
1.32 -
1.33 -
1.34 -void minCut(CutMap& M) : sets M to the characteristic vector of
1.35 - a min cut. M should be a map of bools initialized to false.
1.36 -
1.37 -*/
1.38 -
1.39 -#ifndef PREFLOW_HL0_H
1.40 -#define PREFLOW_HL0_H
1.41 -
1.42 -#include <vector>
1.43 -#include <queue>
1.44 -
1.45 -#include <time_measure.h> //for test
1.46 -
1.47 -namespace hugo {
1.48 -
1.49 - template <typename Graph, typename T,
1.50 - typename FlowMap=typename Graph::EdgeMap<T>,
1.51 - typename CapMap=typename Graph::EdgeMap<T> >
1.52 - class preflow_hl0 {
1.53 -
1.54 - typedef typename Graph::NodeIt NodeIt;
1.55 - typedef typename Graph::EdgeIt EdgeIt;
1.56 - typedef typename Graph::EachNodeIt EachNodeIt;
1.57 - typedef typename Graph::OutEdgeIt OutEdgeIt;
1.58 - typedef typename Graph::InEdgeIt InEdgeIt;
1.59 -
1.60 - Graph& G;
1.61 - NodeIt s;
1.62 - NodeIt t;
1.63 - FlowMap flow;
1.64 - CapMap& capacity;
1.65 - T value;
1.66 -
1.67 - public:
1.68 - double time;
1.69 -
1.70 - preflow_hl0(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity ) :
1.71 - G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) {
1.72 -
1.73 - bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd
1.74 - int n=G.nodeNum();
1.75 - bool end=false;
1.76 - /*
1.77 - 'true' means no active nodes are above bound b.
1.78 - */
1.79 - int k=n-2; //bound on the highest level under n containing a node
1.80 - int b=k; //bound on the highest level under n of an active node
1.81 - /*
1.82 - b is a bound on the highest level of the stack.
1.83 - k is a bound on the highest nonempty level i < n.
1.84 - */
1.85 -
1.86 - typename Graph::NodeMap<int> level(G,n);
1.87 - typename Graph::NodeMap<T> excess(G);
1.88 -
1.89 - std::vector<NodeIt> active(n);
1.90 - typename Graph::NodeMap<NodeIt> next(G);
1.91 - //Stack of the active nodes in level i < n.
1.92 - //We use it in both phases.
1.93 -
1.94 - typename Graph::NodeMap<NodeIt> left(G);
1.95 - typename Graph::NodeMap<NodeIt> right(G);
1.96 - std::vector<NodeIt> level_list(n);
1.97 - /*
1.98 - List of the nodes in level i<n.
1.99 - */
1.100 -
1.101 - /*Reverse_bfs from t, to find the starting level.*/
1.102 - level.set(t,0);
1.103 - std::queue<NodeIt> bfs_queue;
1.104 - bfs_queue.push(t);
1.105 -
1.106 - while (!bfs_queue.empty()) {
1.107 -
1.108 - NodeIt v=bfs_queue.front();
1.109 - bfs_queue.pop();
1.110 - int l=level.get(v)+1;
1.111 -
1.112 - for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
1.113 - NodeIt w=G.tail(e);
1.114 - if ( level.get(w) == n && w != s ) {
1.115 - bfs_queue.push(w);
1.116 - NodeIt first=level_list[l];
1.117 - if ( first != 0 ) left.set(first,w);
1.118 - right.set(w,first);
1.119 - level_list[l]=w;
1.120 - level.set(w, l);
1.121 - }
1.122 - }
1.123 - }
1.124 -
1.125 - level.set(s,n);
1.126 -
1.127 -
1.128 - /* Starting flow. It is everywhere 0 at the moment. */
1.129 - for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
1.130 - {
1.131 - T c=capacity.get(e);
1.132 - if ( c == 0 ) continue;
1.133 - NodeIt w=G.head(e);
1.134 - if ( level.get(w) < n ) {
1.135 - if ( excess.get(w) == 0 && w!=t ) {
1.136 - next.set(w,active[level.get(w)]);
1.137 - active[level.get(w)]=w;
1.138 - }
1.139 - flow.set(e, c);
1.140 - excess.set(w, excess.get(w)+c);
1.141 - }
1.142 - }
1.143 -
1.144 - /*
1.145 - End of preprocessing
1.146 - */
1.147 -
1.148 -
1.149 -
1.150 - /*
1.151 - Push/relabel on the highest level active nodes.
1.152 - */
1.153 - while ( true ) {
1.154 -
1.155 - if ( b == 0 ) {
1.156 - if ( phase ) break;
1.157 -
1.158 - if ( !end && k > 0 ) {
1.159 - b=k;
1.160 - end=true;
1.161 - } else {
1.162 - phase=1;
1.163 - time=currTime();
1.164 - level.set(s,0);
1.165 - std::queue<NodeIt> bfs_queue;
1.166 - bfs_queue.push(s);
1.167 -
1.168 - while (!bfs_queue.empty()) {
1.169 -
1.170 - NodeIt v=bfs_queue.front();
1.171 - bfs_queue.pop();
1.172 - int l=level.get(v)+1;
1.173 -
1.174 - for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
1.175 - if ( capacity.get(e) == flow.get(e) ) continue;
1.176 - NodeIt u=G.tail(e);
1.177 - if ( level.get(u) >= n ) {
1.178 - bfs_queue.push(u);
1.179 - level.set(u, l);
1.180 - if ( excess.get(u) > 0 ) {
1.181 - next.set(u,active[l]);
1.182 - active[l]=u;
1.183 - }
1.184 - }
1.185 - }
1.186 -
1.187 - for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
1.188 - if ( 0 == flow.get(e) ) continue;
1.189 - NodeIt u=G.head(e);
1.190 - if ( level.get(u) >= n ) {
1.191 - bfs_queue.push(u);
1.192 - level.set(u, l);
1.193 - if ( excess.get(u) > 0 ) {
1.194 - next.set(u,active[l]);
1.195 - active[l]=u;
1.196 - }
1.197 - }
1.198 - }
1.199 - }
1.200 - b=n-2;
1.201 - }
1.202 -
1.203 - }
1.204 -
1.205 -
1.206 - if ( active[b] == 0 ) --b;
1.207 - else {
1.208 - end=false;
1.209 -
1.210 - NodeIt w=active[b];
1.211 - active[b]=next.get(w);
1.212 - int lev=level.get(w);
1.213 - T exc=excess.get(w);
1.214 - int newlevel=n; //bound on the next level of w
1.215 -
1.216 - for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
1.217 -
1.218 - if ( flow.get(e) == capacity.get(e) ) continue;
1.219 - NodeIt v=G.head(e);
1.220 - //e=wv
1.221 -
1.222 - if( lev > level.get(v) ) {
1.223 - /*Push is allowed now*/
1.224 -
1.225 - if ( excess.get(v)==0 && v!=t && v!=s ) {
1.226 - int lev_v=level.get(v);
1.227 - next.set(v,active[lev_v]);
1.228 - active[lev_v]=v;
1.229 - }
1.230 -
1.231 - T cap=capacity.get(e);
1.232 - T flo=flow.get(e);
1.233 - T remcap=cap-flo;
1.234 -
1.235 - if ( remcap >= exc ) {
1.236 - /*A nonsaturating push.*/
1.237 -
1.238 - flow.set(e, flo+exc);
1.239 - excess.set(v, excess.get(v)+exc);
1.240 - exc=0;
1.241 - break;
1.242 -
1.243 - } else {
1.244 - /*A saturating push.*/
1.245 -
1.246 - flow.set(e, cap);
1.247 - excess.set(v, excess.get(v)+remcap);
1.248 - exc-=remcap;
1.249 - }
1.250 - } else if ( newlevel > level.get(v) ){
1.251 - newlevel = level.get(v);
1.252 - }
1.253 -
1.254 - } //for out edges wv
1.255 -
1.256 -
1.257 - if ( exc > 0 ) {
1.258 - for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
1.259 -
1.260 - if( flow.get(e) == 0 ) continue;
1.261 - NodeIt v=G.tail(e);
1.262 - //e=vw
1.263 -
1.264 - if( lev > level.get(v) ) {
1.265 - /*Push is allowed now*/
1.266 -
1.267 - if ( excess.get(v)==0 && v!=t && v!=s ) {
1.268 - int lev_v=level.get(v);
1.269 - next.set(v,active[lev_v]);
1.270 - active[lev_v]=v;
1.271 - }
1.272 -
1.273 - T flo=flow.get(e);
1.274 -
1.275 - if ( flo >= exc ) {
1.276 - /*A nonsaturating push.*/
1.277 -
1.278 - flow.set(e, flo-exc);
1.279 - excess.set(v, excess.get(v)+exc);
1.280 - exc=0;
1.281 - break;
1.282 - } else {
1.283 - /*A saturating push.*/
1.284 -
1.285 - excess.set(v, excess.get(v)+flo);
1.286 - exc-=flo;
1.287 - flow.set(e,0);
1.288 - }
1.289 - } else if ( newlevel > level.get(v) ) {
1.290 - newlevel = level.get(v);
1.291 - }
1.292 - } //for in edges vw
1.293 -
1.294 - } // if w still has excess after the out edge for cycle
1.295 -
1.296 - excess.set(w, exc);
1.297 -
1.298 - /*
1.299 - Relabel
1.300 - */
1.301 -
1.302 -
1.303 - if ( exc > 0 ) {
1.304 - //now 'lev' is the old level of w
1.305 -
1.306 - if ( phase ) {
1.307 - level.set(w,++newlevel);
1.308 - next.set(w,active[newlevel]);
1.309 - active[newlevel]=w;
1.310 - b=newlevel;
1.311 - } else {
1.312 - //unlacing starts
1.313 - NodeIt right_n=right.get(w);
1.314 - NodeIt left_n=left.get(w);
1.315 -
1.316 - if ( right_n != 0 ) {
1.317 - if ( left_n != 0 ) {
1.318 - right.set(left_n, right_n);
1.319 - left.set(right_n, left_n);
1.320 - } else {
1.321 - level_list[lev]=right_n;
1.322 - left.set(right_n, 0);
1.323 - }
1.324 - } else {
1.325 - if ( left_n != 0 ) {
1.326 - right.set(left_n, 0);
1.327 - } else {
1.328 - level_list[lev]=0;
1.329 -
1.330 - }
1.331 - }
1.332 - //unlacing ends
1.333 -
1.334 - //gapping starts
1.335 - if ( level_list[lev]==0 ) {
1.336 -
1.337 - for (int i=lev; i!=k ; ) {
1.338 - NodeIt v=level_list[++i];
1.339 - while ( v != 0 ) {
1.340 - level.set(v,n);
1.341 - v=right.get(v);
1.342 - }
1.343 - level_list[i]=0;
1.344 - active[i]=0;
1.345 - }
1.346 -
1.347 - level.set(w,n);
1.348 - b=lev-1;
1.349 - k=b;
1.350 - //gapping ends
1.351 - } else {
1.352 -
1.353 - if ( newlevel == n ) level.set(w,n);
1.354 - else {
1.355 - level.set(w,++newlevel);
1.356 - next.set(w,active[newlevel]);
1.357 - active[newlevel]=w;
1.358 - if ( k < newlevel ) ++k;
1.359 - NodeIt first=level_list[newlevel];
1.360 - if ( first != 0 ) left.set(first,w);
1.361 - right.set(w,first);
1.362 - left.set(w,0);
1.363 - level_list[newlevel]=w;
1.364 - }
1.365 - }
1.366 -
1.367 -
1.368 - } //phase 0
1.369 -
1.370 -
1.371 - } // if ( exc > 0 )
1.372 -
1.373 -
1.374 - } // if stack[b] is nonempty
1.375 -
1.376 - } // while(true)
1.377 -
1.378 -
1.379 - value = excess.get(t);
1.380 - /*Max flow value.*/
1.381 -
1.382 - } //void run()
1.383 -
1.384 -
1.385 -
1.386 -
1.387 -
1.388 - /*
1.389 - Returns the maximum value of a flow.
1.390 - */
1.391 -
1.392 - T maxFlow() {
1.393 - return value;
1.394 - }
1.395 -
1.396 -
1.397 -
1.398 - /*
1.399 - For the maximum flow x found by the algorithm,
1.400 - it returns the flow value on edge e, i.e. x(e).
1.401 - */
1.402 -
1.403 - T flowOnEdge(EdgeIt e) {
1.404 - return flow.get(e);
1.405 - }
1.406 -
1.407 -
1.408 -
1.409 - FlowMap Flow() {
1.410 - return flow;
1.411 - }
1.412 -
1.413 -
1.414 -
1.415 - void Flow(FlowMap& _flow ) {
1.416 - for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v)
1.417 - _flow.set(v,flow.get(v));
1.418 - }
1.419 -
1.420 -
1.421 -
1.422 - /*
1.423 - Returns the minimum min cut, by a bfs from s in the residual graph.
1.424 - */
1.425 -
1.426 - template<typename _CutMap>
1.427 - void minMinCut(_CutMap& M) {
1.428 -
1.429 - std::queue<NodeIt> queue;
1.430 -
1.431 - M.set(s,true);
1.432 - queue.push(s);
1.433 -
1.434 - while (!queue.empty()) {
1.435 - NodeIt w=queue.front();
1.436 - queue.pop();
1.437 -
1.438 - for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
1.439 - NodeIt v=G.head(e);
1.440 - if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
1.441 - queue.push(v);
1.442 - M.set(v, true);
1.443 - }
1.444 - }
1.445 -
1.446 - for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
1.447 - NodeIt v=G.tail(e);
1.448 - if (!M.get(v) && flow.get(e) > 0 ) {
1.449 - queue.push(v);
1.450 - M.set(v, true);
1.451 - }
1.452 - }
1.453 - }
1.454 - }
1.455 -
1.456 -
1.457 -
1.458 - /*
1.459 - Returns the maximum min cut, by a reverse bfs
1.460 - from t in the residual graph.
1.461 - */
1.462 -
1.463 - template<typename _CutMap>
1.464 - void maxMinCut(_CutMap& M) {
1.465 -
1.466 - std::queue<NodeIt> queue;
1.467 -
1.468 - M.set(t,true);
1.469 - queue.push(t);
1.470 -
1.471 - while (!queue.empty()) {
1.472 - NodeIt w=queue.front();
1.473 - queue.pop();
1.474 -
1.475 - for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
1.476 - NodeIt v=G.tail(e);
1.477 - if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
1.478 - queue.push(v);
1.479 - M.set(v, true);
1.480 - }
1.481 - }
1.482 -
1.483 - for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
1.484 - NodeIt v=G.head(e);
1.485 - if (!M.get(v) && flow.get(e) > 0 ) {
1.486 - queue.push(v);
1.487 - M.set(v, true);
1.488 - }
1.489 - }
1.490 - }
1.491 -
1.492 - for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
1.493 - M.set(v, !M.get(v));
1.494 - }
1.495 -
1.496 - }
1.497 -
1.498 -
1.499 -
1.500 - template<typename _CutMap>
1.501 - void minCut(_CutMap& M) {
1.502 - minMinCut(M);
1.503 - }
1.504 -
1.505 - };
1.506 -}//namespace marci
1.507 -#endif
1.508 -
1.509 -
1.510 -
1.511 -