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