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