1.1 --- a/src/work/jacint/preflow_push_hl.h Wed Feb 18 13:06:41 2004 +0000
1.2 +++ b/src/work/jacint/preflow_push_hl.h Wed Feb 18 14:42:38 2004 +0000
1.3 @@ -3,7 +3,9 @@
1.4 preflow_push_hl.h
1.5 by jacint.
1.6 Runs the highest label variant of the preflow push algorithm with
1.7 -running time O(n^2\sqrt(m)).
1.8 +running time O(n^2\sqrt(m)), and with the 'empty level' heuristic.
1.9 +
1.10 +'A' is a parameter for the empty_level heuristic
1.11
1.12 Member functions:
1.13
1.14 @@ -15,11 +17,17 @@
1.15
1.16 T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
1.17
1.18 -Graph::EdgeMap<T> allflow() : returns the fixed maximum flow x
1.19 +FlowMap allflow() : returns the fixed maximum flow x
1.20
1.21 -Graph::NodeMap<bool> mincut() : returns a
1.22 - characteristic vector of a minimum cut. (An empty level
1.23 - in the algorithm gives a minimum cut.)
1.24 +void mincut(CutMap& M) : sets M to the characteristic vector of a
1.25 + minimum cut. M should be a map of bools initialized to false.
1.26 +
1.27 +void min_mincut(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 max_mincut(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
1.35 #ifndef PREFLOW_PUSH_HL_H
1.36 @@ -29,12 +37,13 @@
1.37
1.38 #include <vector>
1.39 #include <stack>
1.40 -
1.41 -#include <reverse_bfs.h>
1.42 +#include <queue>
1.43
1.44 namespace marci {
1.45
1.46 - template <typename Graph, typename T>
1.47 + template <typename Graph, typename T,
1.48 + typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>,
1.49 + typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
1.50 class preflow_push_hl {
1.51
1.52 typedef typename Graph::NodeIt NodeIt;
1.53 @@ -46,207 +55,207 @@
1.54 Graph& G;
1.55 NodeIt s;
1.56 NodeIt t;
1.57 - typename Graph::EdgeMap<T> flow;
1.58 - typename Graph::EdgeMap<T> capacity;
1.59 + FlowMap flow;
1.60 + CapMap& capacity;
1.61 T value;
1.62 - typename Graph::NodeMap<bool> mincutvector;
1.63 -
1.64 +
1.65 public:
1.66
1.67 - preflow_push_hl(Graph& _G, NodeIt _s, NodeIt _t,
1.68 - typename Graph::EdgeMap<T>& _capacity) :
1.69 - G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity),
1.70 - mincutvector(_G, true) { }
1.71 + preflow_push_hl(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
1.72 + G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
1.73
1.74
1.75 - /*
1.76 - The run() function runs the highest label preflow-push,
1.77 - running time: O(n^2\sqrt(m))
1.78 - */
1.79 +
1.80 +
1.81 void run() {
1.82
1.83 - std::cout<<"A is "<<A<<" ";
1.84 -
1.85 - typename Graph::NodeMap<int> level(G);
1.86 - typename Graph::NodeMap<T> excess(G);
1.87 -
1.88 int n=G.nodeNum();
1.89 int b=n-2;
1.90 /*
1.91 b is a bound on the highest level of an active node.
1.92 - In the beginning it is at most n-2.
1.93 */
1.94
1.95 - std::vector<int> numb(n); //The number of nodes on level i < n.
1.96 + IntMap level(G,n);
1.97 + TMap excess(G);
1.98 +
1.99 + std::vector<int> numb(n);
1.100 + /*
1.101 + The number of nodes on level i < n. It is
1.102 + initialized to n+1, because of the reverse_bfs-part.
1.103 + */
1.104 +
1.105 std::vector<std::stack<NodeIt> > stack(2*n-1);
1.106 //Stack of the active nodes in level i.
1.107
1.108
1.109 /*Reverse_bfs from t, to find the starting level.*/
1.110 - reverse_bfs<Graph> bfs(G, t);
1.111 - bfs.run();
1.112 - for(EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v)
1.113 - {
1.114 - int dist=bfs.dist(v);
1.115 - level.set(v, dist);
1.116 - ++numb[dist];
1.117 + level.set(t,0);
1.118 + std::queue<NodeIt> bfs_queue;
1.119 + bfs_queue.push(t);
1.120 +
1.121 + while (!bfs_queue.empty()) {
1.122 +
1.123 + NodeIt v=bfs_queue.front();
1.124 + bfs_queue.pop();
1.125 + int l=level.get(v)+1;
1.126 +
1.127 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
1.128 + NodeIt w=G.tail(e);
1.129 + if ( level.get(w) == n ) {
1.130 + bfs_queue.push(w);
1.131 + ++numb[l];
1.132 + level.set(w, l);
1.133 + }
1.134 }
1.135 + }
1.136 +
1.137 + level.set(s,n);
1.138
1.139 - level.set(s,n);
1.140
1.141
1.142 /* Starting flow. It is everywhere 0 at the moment. */
1.143 for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
1.144 {
1.145 - if ( capacity.get(e) > 0 ) {
1.146 - NodeIt w=G.head(e);
1.147 - if ( w!=s ) {
1.148 - if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w);
1.149 - flow.set(e, capacity.get(e));
1.150 - excess.set(w, excess.get(w)+capacity.get(e));
1.151 - }
1.152 + if ( capacity.get(e) == 0 ) continue;
1.153 + NodeIt w=G.head(e);
1.154 + if ( w!=s ) {
1.155 + if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w);
1.156 + flow.set(e, capacity.get(e));
1.157 + excess.set(w, excess.get(w)+capacity.get(e));
1.158 }
1.159 }
1.160 -
1.161 +
1.162 /*
1.163 End of preprocessing
1.164 */
1.165
1.166
1.167 -
1.168 /*
1.169 Push/relabel on the highest level active nodes.
1.170 */
1.171 -
1.172 /*While there exists an active node.*/
1.173 while (b) {
1.174 -
1.175 - /*We decrease the bound if there is no active node of level b.*/
1.176 - if (stack[b].empty()) {
1.177 + if ( stack[b].empty() ) {
1.178 --b;
1.179 - } else {
1.180 -
1.181 - NodeIt w=stack[b].top(); //w is a highest label active node.
1.182 - stack[b].pop();
1.183 + continue;
1.184 + }
1.185
1.186 - int newlevel=2*n-2; //In newlevel we bound the next level of w.
1.187 + NodeIt w=stack[b].top(); //w is a highest label active node.
1.188 + stack[b].pop();
1.189 + int lev=level.get(w);
1.190 + int exc=excess.get(w);
1.191 + int newlevel=2*n-2; //In newlevel we bound the next level of w.
1.192
1.193 + // if ( level.get(w) < n ) { //Nem tudom ez mukodik-e
1.194 for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
1.195
1.196 - if ( flow.get(e) < capacity.get(e) ) {
1.197 - /*e is an edge of the residual graph */
1.198 + if ( flow.get(e) == capacity.get(e) ) continue;
1.199 + NodeIt v=G.head(e);
1.200 + //e=wv
1.201 +
1.202 + if( lev > level.get(v) ) {
1.203 + /*Push is allowed now*/
1.204 +
1.205 + if ( excess.get(v)==0 && v != s && v !=t )
1.206 + stack[level.get(v)].push(v);
1.207 + /*v becomes active.*/
1.208 +
1.209 + int cap=capacity.get(e);
1.210 + int flo=flow.get(e);
1.211 + int remcap=cap-flo;
1.212 +
1.213 + if ( remcap >= exc ) {
1.214 + /*A nonsaturating push.*/
1.215 +
1.216 + flow.set(e, flo+exc);
1.217 + excess.set(v, excess.get(v)+exc);
1.218 + exc=0;
1.219 + break;
1.220 +
1.221 + } else {
1.222 + /*A saturating push.*/
1.223 +
1.224 + flow.set(e, cap );
1.225 + excess.set(v, excess.get(v)+remcap);
1.226 + exc-=remcap;
1.227 + }
1.228 + } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
1.229 +
1.230 + } //for out edges wv
1.231 +
1.232 +
1.233 + if ( exc > 0 ) {
1.234 + for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
1.235 +
1.236 + if( flow.get(e) == 0 ) continue;
1.237 + NodeIt v=G.tail(e);
1.238 + //e=vw
1.239 +
1.240 + if( lev > level.get(v) ) {
1.241 + /*Push is allowed now*/
1.242 +
1.243 + if ( excess.get(v)==0 && v != s && v !=t)
1.244 + stack[level.get(v)].push(v);
1.245 + /*v becomes active.*/
1.246 +
1.247 + int flo=flow.get(e);
1.248 +
1.249 + if ( flo >= exc ) {
1.250 + /*A nonsaturating push.*/
1.251 +
1.252 + flow.set(e, flo-exc);
1.253 + excess.set(v, excess.get(v)+exc);
1.254 + exc=0;
1.255 + break;
1.256 + } else {
1.257 + /*A saturating push.*/
1.258 +
1.259 + excess.set(v, excess.get(v)+flo);
1.260 + exc-=flo;
1.261 + flow.set(e,0);
1.262 + }
1.263 + } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
1.264 +
1.265 + } //for in edges vw
1.266 +
1.267 + } // if w still has excess after the out edge for cycle
1.268 +
1.269 + excess.set(w, exc);
1.270 +
1.271
1.272 - NodeIt v=G.head(e); /*e is the edge wv.*/
1.273 +
1.274
1.275 - if( level.get(w) == level.get(v)+1 ) {
1.276 - /*Push is allowed now*/
1.277 -
1.278 - if ( excess.get(v)==0 && v != s && v !=t ) stack[level.get(v)].push(v);
1.279 - /*v becomes active.*/
1.280 -
1.281 - if ( capacity.get(e)-flow.get(e) > excess.get(w) ) {
1.282 - /*A nonsaturating push.*/
1.283 -
1.284 - flow.set(e, flow.get(e)+excess.get(w));
1.285 - excess.set(v, excess.get(v)+excess.get(w));
1.286 - excess.set(w,0);
1.287 - break;
1.288 -
1.289 - } else {
1.290 - /*A saturating push.*/
1.291 -
1.292 - excess.set(v, excess.get(v)+capacity.get(e)-flow.get(e));
1.293 - excess.set(w, excess.get(w)-capacity.get(e)+flow.get(e));
1.294 - flow.set(e, capacity.get(e));
1.295 - if ( excess.get(w)==0 ) break;
1.296 - /*If w is not active any more, then we go on to the next node.*/
1.297 -
1.298 - }
1.299 - } else {
1.300 - newlevel = newlevel < level.get(v) ? newlevel : level.get(v);
1.301 + /*
1.302 + Relabel
1.303 + */
1.304 +
1.305 + if ( exc > 0 ) {
1.306 + //now 'lev' is the old level of w
1.307 + level.set(w,++newlevel);
1.308 +
1.309 + if ( lev < n ) {
1.310 + --numb[lev];
1.311 +
1.312 + if ( !numb[lev] && lev < A*n ) { //If the level of w gets empty.
1.313 +
1.314 + for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
1.315 + if (level.get(v) > lev && level.get(v) < n ) level.set(v,n);
1.316 }
1.317 -
1.318 - } //if the out edge wv is in the res graph
1.319 -
1.320 - } //for out edges wv
1.321 + for (int i=lev+1 ; i!=n ; ++i) numb[i]=0;
1.322 + if ( newlevel < n ) newlevel=n;
1.323 + } else {
1.324 + if ( newlevel < n ) ++numb[newlevel];
1.325 + }
1.326 + }
1.327
1.328 -
1.329 - if ( excess.get(w) > 0 ) {
1.330 -
1.331 - for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
1.332 - NodeIt v=G.tail(e); /*e is the edge vw.*/
1.333 -
1.334 - if( flow.get(e) > 0 ) {
1.335 - /*e is an edge of the residual graph */
1.336 -
1.337 - if( level.get(w)==level.get(v)+1 ) {
1.338 - /*Push is allowed now*/
1.339 -
1.340 - if ( excess.get(v)==0 && v != s && v !=t) stack[level.get(v)].push(v);
1.341 - /*v becomes active.*/
1.342 -
1.343 - if ( flow.get(e) > excess.get(w) ) {
1.344 - /*A nonsaturating push.*/
1.345 -
1.346 - flow.set(e, flow.get(e)-excess.get(w));
1.347 - excess.set(v, excess.get(v)+excess.get(w));
1.348 - excess.set(w,0);
1.349 - break;
1.350 - } else {
1.351 - /*A saturating push.*/
1.352 -
1.353 - excess.set(v, excess.get(v)+flow.get(e));
1.354 - excess.set(w, excess.get(w)-flow.get(e));
1.355 - flow.set(e,0);
1.356 - if ( excess.get(w)==0 ) break;
1.357 - }
1.358 - } else {
1.359 - newlevel = newlevel < level.get(v) ? newlevel : level.get(v);
1.360 - }
1.361 -
1.362 - } //if in edge vw is in the res graph
1.363 -
1.364 - } //for in edges vw
1.365 -
1.366 - } // if w still has excess after the out edge for cycle
1.367 -
1.368 -
1.369 - /*
1.370 - Relabel
1.371 - */
1.372 + stack[newlevel].push(w);
1.373 + b=newlevel;
1.374
1.375 - if ( excess.get(w) > 0 ) {
1.376 -
1.377 - int oldlevel=level.get(w);
1.378 - level.set(w,++newlevel);
1.379 -
1.380 - if ( oldlevel < n ) {
1.381 - --numb[oldlevel];
1.382 -
1.383 - if ( !numb[oldlevel] && oldlevel < A*n ) { //If the level of w gets empty.
1.384 -
1.385 - for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
1.386 - if (level.get(v) > oldlevel && level.get(v) < n ) level.set(v,n);
1.387 - }
1.388 - for (int i=oldlevel+1 ; i!=n ; ++i) numb[i]=0;
1.389 - if ( newlevel < n ) newlevel=n;
1.390 - } else {
1.391 - if ( newlevel < n ) ++numb[newlevel];
1.392 - }
1.393 - } else {
1.394 - if ( newlevel < n ) ++numb[newlevel];
1.395 - }
1.396 -
1.397 - stack[newlevel].push(w);
1.398 - b=newlevel;
1.399 -
1.400 - }
1.401 -
1.402 - } // if stack[b] is nonempty
1.403 -
1.404 + }
1.405 +
1.406 } // while(b)
1.407 -
1.408 -
1.409 +
1.410 +
1.411 value = excess.get(t);
1.412 /*Max flow value.*/
1.413
1.414 @@ -271,7 +280,7 @@
1.415 For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e).
1.416 */
1.417
1.418 - T flowonedge(EdgeIt e) {
1.419 + T flowonedge(const EdgeIt e) {
1.420 return flow.get(e);
1.421 }
1.422
1.423 @@ -281,21 +290,61 @@
1.424 Returns the maximum flow x found by the algorithm.
1.425 */
1.426
1.427 - typename Graph::EdgeMap<T> allflow() {
1.428 + FlowMap allflow() {
1.429 return flow;
1.430 }
1.431
1.432
1.433
1.434 +
1.435 /*
1.436 - Returns a minimum cut by using a reverse bfs from t in the residual graph.
1.437 + Returns the minimum min cut, by a bfs from s in the residual graph.
1.438 */
1.439
1.440 - typename Graph::NodeMap<bool> mincut() {
1.441 + template<typename CutMap>
1.442 + void mincut(CutMap& M) {
1.443
1.444 std::queue<NodeIt> queue;
1.445
1.446 - mincutvector.set(t,false);
1.447 + M.set(s,true);
1.448 + queue.push(s);
1.449 +
1.450 + while (!queue.empty()) {
1.451 + NodeIt w=queue.front();
1.452 + queue.pop();
1.453 +
1.454 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
1.455 + NodeIt v=G.head(e);
1.456 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
1.457 + queue.push(v);
1.458 + M.set(v, true);
1.459 + }
1.460 + }
1.461 +
1.462 + for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
1.463 + NodeIt v=G.tail(e);
1.464 + if (!M.get(v) && flow.get(e) > 0 ) {
1.465 + queue.push(v);
1.466 + M.set(v, true);
1.467 + }
1.468 + }
1.469 +
1.470 + }
1.471 + }
1.472 +
1.473 +
1.474 +
1.475 + /*
1.476 + Returns the maximum min cut, by a reverse bfs
1.477 + from t in the residual graph.
1.478 + */
1.479 +
1.480 + template<typename CutMap>
1.481 + void max_mincut(CutMap& M) {
1.482 +
1.483 + std::queue<NodeIt> queue;
1.484 +
1.485 + M.set(t,true);
1.486 queue.push(t);
1.487
1.488 while (!queue.empty()) {
1.489 @@ -304,25 +353,36 @@
1.490
1.491 for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
1.492 NodeIt v=G.tail(e);
1.493 - if (mincutvector.get(v) && flow.get(e) < capacity.get(e) ) {
1.494 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
1.495 queue.push(v);
1.496 - mincutvector.set(v, false);
1.497 + M.set(v, true);
1.498 }
1.499 - } // for
1.500 + }
1.501
1.502 for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
1.503 NodeIt v=G.head(e);
1.504 - if (mincutvector.get(v) && flow.get(e) > 0 ) {
1.505 + if (!M.get(v) && flow.get(e) > 0 ) {
1.506 queue.push(v);
1.507 - mincutvector.set(v, false);
1.508 + M.set(v, true);
1.509 }
1.510 - } // for
1.511 -
1.512 + }
1.513 }
1.514
1.515 - return mincutvector;
1.516 -
1.517 + for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
1.518 + M.set(v, !M.get(v));
1.519 + }
1.520 +
1.521 }
1.522 +
1.523 +
1.524 +
1.525 + template<typename CutMap>
1.526 + void min_mincut(CutMap& M) {
1.527 + mincut(M);
1.528 + }
1.529 +
1.530 +
1.531 +
1.532 };
1.533 }//namespace marci
1.534 #endif
2.1 --- a/src/work/jacint/preflow_push_max_flow.h Wed Feb 18 13:06:41 2004 +0000
2.2 +++ b/src/work/jacint/preflow_push_max_flow.h Wed Feb 18 14:42:38 2004 +0000
2.3 @@ -1,3 +1,4 @@
2.4 +// -*- C++ -*-
2.5 /*
2.6 preflow_push_max_flow_h
2.7 by jacint.
2.8 @@ -15,13 +16,16 @@
2.9
2.10 T maxflow() : returns the value of a maximum flow
2.11
2.12 -NodeMap<bool> mincut(): returns a
2.13 - characteristic vector of a minimum cut.
2.14 +void mincut(CutMap& M) : sets M to the characteristic vector of a
2.15 + minimum cut. M should be a map of bools initialized to false.
2.16 +
2.17 */
2.18
2.19 #ifndef PREFLOW_PUSH_MAX_FLOW_H
2.20 #define PREFLOW_PUSH_MAX_FLOW_H
2.21
2.22 +#define A 1
2.23 +
2.24 #include <algorithm>
2.25 #include <vector>
2.26 #include <stack>
2.27 @@ -31,7 +35,9 @@
2.28
2.29 namespace marci {
2.30
2.31 - template <typename Graph, typename T>
2.32 + template <typename Graph, typename T,
2.33 + typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>,
2.34 + typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
2.35 class preflow_push_max_flow {
2.36
2.37 typedef typename Graph::NodeIt NodeIt;
2.38 @@ -42,17 +48,15 @@
2.39 Graph& G;
2.40 NodeIt s;
2.41 NodeIt t;
2.42 - typename Graph::EdgeMap<T>& capacity;
2.43 - T value;
2.44 - typename Graph::NodeMap<bool> mincutvector;
2.45 -
2.46 -
2.47 -
2.48 + IntMap level;
2.49 + CapMap& capacity;
2.50 + int empty_level; //an empty level in the end of run()
2.51 + T value;
2.52 +
2.53 public:
2.54 -
2.55 - preflow_push_max_flow ( Graph& _G, NodeIt _s, NodeIt _t,
2.56 - typename Graph::EdgeMap<T>& _capacity) :
2.57 - G(_G), s(_s), t(_t), capacity(_capacity), mincutvector(_G, false) { }
2.58 +
2.59 + preflow_push_max_flow(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
2.60 + G(_G), s(_s), t(_t), level(_G), capacity(_capacity) { }
2.61
2.62
2.63 /*
2.64 @@ -62,223 +66,200 @@
2.65 */
2.66 void run() {
2.67
2.68 - typename Graph::EdgeMap<T> flow(G, 0);
2.69 - typename Graph::NodeMap<int> level(G);
2.70 - typename Graph::NodeMap<T> excess(G);
2.71 -
2.72 - int n=G.nodeNum();
2.73 + int n=G.nodeNum();
2.74 int b=n-2;
2.75 /*
2.76 - b is a bound on the highest level of an active Node.
2.77 - In the beginning it is at most n-2.
2.78 + b is a bound on the highest level of an active node.
2.79 */
2.80 -
2.81 - std::vector<int> numb(n); //The number of Nodes on level i < n.
2.82 - std::vector<std::stack<NodeIt> > stack(2*n-1);
2.83 - //Stack of the active Nodes in level i.
2.84 +
2.85 + IntMap level(G,n);
2.86 + TMap excess(G);
2.87 + FlowMap flow(G,0);
2.88 +
2.89 + std::vector<int> numb(n);
2.90 + /*
2.91 + The number of nodes on level i < n. It is
2.92 + initialized to n+1, because of the reverse_bfs-part.
2.93 + */
2.94 +
2.95 + std::vector<std::stack<NodeIt> > stack(n);
2.96 + //Stack of the active nodes in level i.
2.97 +
2.98
2.99 /*Reverse_bfs from t, to find the starting level.*/
2.100 - reverse_bfs<Graph> bfs(G, t);
2.101 - bfs.run();
2.102 - for(EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v)
2.103 - {
2.104 - int dist=bfs.dist(v);
2.105 - level.set(v, dist);
2.106 - ++numb[dist];
2.107 + level.set(t,0);
2.108 + std::queue<NodeIt> bfs_queue;
2.109 + bfs_queue.push(t);
2.110 +
2.111 + while (!bfs_queue.empty()) {
2.112 +
2.113 + NodeIt v=bfs_queue.front();
2.114 + bfs_queue.pop();
2.115 + int l=level.get(v)+1;
2.116 +
2.117 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
2.118 + NodeIt w=G.tail(e);
2.119 + if ( level.get(w) == n ) {
2.120 + bfs_queue.push(w);
2.121 + ++numb[l];
2.122 + level.set(w, l);
2.123 + }
2.124 }
2.125 -
2.126 + }
2.127 +
2.128 level.set(s,n);
2.129
2.130 - /* Starting flow. It is everywhere 0 at the moment. */
2.131 +
2.132 +
2.133 + /* Starting flow. It is everywhere 0 at the moment. */
2.134 for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
2.135 {
2.136 - if ( capacity.get(e) > 0 ) {
2.137 - NodeIt w=G.head(e);
2.138 + if ( capacity.get(e) == 0 ) continue;
2.139 + NodeIt w=G.head(e);
2.140 + if ( level.get(w) < n ) {
2.141 + if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w);
2.142 flow.set(e, capacity.get(e));
2.143 - stack[level.get(w)].push(w);
2.144 excess.set(w, excess.get(w)+capacity.get(e));
2.145 }
2.146 }
2.147 -
2.148 +
2.149 /*
2.150 End of preprocessing
2.151 */
2.152
2.153
2.154 + /*
2.155 + Push/relabel on the highest level active nodes.
2.156 + */
2.157 + /*While there exists an active node.*/
2.158 + while (b) {
2.159 + if ( stack[b].empty() ) {
2.160 + --b;
2.161 + continue;
2.162 + }
2.163 +
2.164 + NodeIt w=stack[b].top(); //w is a highest label active node.
2.165 + stack[b].pop();
2.166 + int lev=level.get(w);
2.167 + int exc=excess.get(w);
2.168 + int newlevel=2*n-2; //In newlevel we bound the next level of w.
2.169 +
2.170 + // if ( level.get(w) < n ) { //Nem tudom ez mukodik-e
2.171 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
2.172 +
2.173 + if ( flow.get(e) == capacity.get(e) ) continue;
2.174 + NodeIt v=G.head(e);
2.175 + //e=wv
2.176 +
2.177 + if( lev > level.get(v) ) {
2.178 + /*Push is allowed now*/
2.179 +
2.180 + if ( excess.get(v)==0 && v != s && v !=t )
2.181 + stack[level.get(v)].push(v);
2.182 + /*v becomes active.*/
2.183 +
2.184 + int cap=capacity.get(e);
2.185 + int flo=flow.get(e);
2.186 + int remcap=cap-flo;
2.187 +
2.188 + if ( remcap >= exc ) {
2.189 + /*A nonsaturating push.*/
2.190 +
2.191 + flow.set(e, flo+exc);
2.192 + excess.set(v, excess.get(v)+exc);
2.193 + exc=0;
2.194 + break;
2.195 +
2.196 + } else {
2.197 + /*A saturating push.*/
2.198 +
2.199 + flow.set(e, cap );
2.200 + excess.set(v, excess.get(v)+remcap);
2.201 + exc-=remcap;
2.202 + }
2.203 + } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
2.204 +
2.205 + } //for out edges wv
2.206 +
2.207 +
2.208 + if ( exc > 0 ) {
2.209 + for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
2.210 +
2.211 + if( flow.get(e) == 0 ) continue;
2.212 + NodeIt v=G.tail(e);
2.213 + //e=vw
2.214 +
2.215 + if( lev > level.get(v) ) {
2.216 + /*Push is allowed now*/
2.217 +
2.218 + if ( excess.get(v)==0 && v != s && v !=t)
2.219 + stack[level.get(v)].push(v);
2.220 + /*v becomes active.*/
2.221 +
2.222 + int flo=flow.get(e);
2.223 +
2.224 + if ( flo >= exc ) {
2.225 + /*A nonsaturating push.*/
2.226 +
2.227 + flow.set(e, flo-exc);
2.228 + excess.set(v, excess.get(v)+exc);
2.229 + exc=0;
2.230 + break;
2.231 + } else {
2.232 + /*A saturating push.*/
2.233 +
2.234 + excess.set(v, excess.get(v)+flo);
2.235 + exc-=flo;
2.236 + flow.set(e,0);
2.237 + }
2.238 + } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
2.239 +
2.240 + } //for in edges vw
2.241 +
2.242 + } // if w still has excess after the out edge for cycle
2.243 +
2.244 + excess.set(w, exc);
2.245 +
2.246 +
2.247 + /*
2.248 + Relabel
2.249 + */
2.250 +
2.251 + if ( exc > 0 ) {
2.252 + //now 'lev' is the old level of w
2.253 + level.set(w,++newlevel);
2.254 + --numb[lev];
2.255 +
2.256 + if ( !numb[lev] && lev < A*n ) { //If the level of w gets empty.
2.257 +
2.258 + for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
2.259 + if (level.get(v) > lev ) level.set(v,n);
2.260 + }
2.261 + for (int i=lev+1 ; i!=n ; ++i) numb[i]=0;
2.262 + if ( newlevel < n ) newlevel=n;
2.263 + } else if ( newlevel < n ) {
2.264 + ++numb[newlevel];
2.265 + stack[newlevel].push(w);
2.266 + b=newlevel;
2.267 + }
2.268 + }
2.269
2.270 - /*
2.271 - Push/relabel on the highest level active Nodes.
2.272 - */
2.273 -
2.274 - /*While there exists an active Node.*/
2.275 - while (b) {
2.276
2.277 - /*We decrease the bound if there is no active node of level b.*/
2.278 - if (stack[b].empty()) {
2.279 - --b;
2.280 - } else {
2.281
2.282 - NodeIt w=stack[b].top(); //w is the highest label active Node.
2.283 - stack[b].pop(); //We delete w from the stack.
2.284 -
2.285 - int newlevel=2*n-2; //In newlevel we maintain the next level of w.
2.286 -
2.287 - for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
2.288 - NodeIt v=G.head(e);
2.289 - /*e is the Edge wv.*/
2.290 -
2.291 - if (flow.get(e)<capacity.get(e)) {
2.292 - /*e is an Edge of the residual graph */
2.293 -
2.294 - if(level.get(w)==level.get(v)+1) {
2.295 - /*Push is allowed now*/
2.296 -
2.297 - if (capacity.get(e)-flow.get(e) > excess.get(w)) {
2.298 - /*A nonsaturating push.*/
2.299 -
2.300 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
2.301 - /*v becomes active.*/
2.302 -
2.303 - flow.set(e, flow.get(e)+excess.get(w));
2.304 - excess.set(v, excess.get(v)+excess.get(w));
2.305 - excess.set(w,0);
2.306 - //std::cout << w << " " << v <<" elore elen nonsat pump " << std::endl;
2.307 - break;
2.308 - } else {
2.309 - /*A saturating push.*/
2.310 -
2.311 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
2.312 - /*v becomes active.*/
2.313 -
2.314 - excess.set(v, excess.get(v)+capacity.get(e)-flow.get(e));
2.315 - excess.set(w, excess.get(w)-capacity.get(e)+flow.get(e));
2.316 - flow.set(e, capacity.get(e));
2.317 - //std::cout << w <<" " << v <<" elore elen sat pump " << std::endl;
2.318 - if (excess.get(w)==0) break;
2.319 - /*If w is not active any more, then we go on to the next Node.*/
2.320 -
2.321 - } // if (capacity.get(e)-flow.get(e) > excess.get(w))
2.322 - } // if (level.get(w)==level.get(v)+1)
2.323 -
2.324 - else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
2.325 -
2.326 - } //if (flow.get(e)<capacity.get(e))
2.327 -
2.328 - } //for(OutEdgeIt e=G.first_OutEdge(w); e.valid(); ++e)
2.329 -
2.330 -
2.331 -
2.332 - for(InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
2.333 - NodeIt v=G.tail(e);
2.334 - /*e is the Edge vw.*/
2.335 -
2.336 - if (excess.get(w)==0) break;
2.337 - /*It may happen, that w became inactive in the first 'for' cycle.*/
2.338 -
2.339 - if(flow.get(e)>0) {
2.340 - /*e is an Edge of the residual graph */
2.341 -
2.342 - if(level.get(w)==level.get(v)+1) {
2.343 - /*Push is allowed now*/
2.344 -
2.345 - if (flow.get(e) > excess.get(w)) {
2.346 - /*A nonsaturating push.*/
2.347 -
2.348 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
2.349 - /*v becomes active.*/
2.350 -
2.351 - flow.set(e, flow.get(e)-excess.get(w));
2.352 - excess.set(v, excess.get(v)+excess.get(w));
2.353 - excess.set(w,0);
2.354 - //std::cout << v << " " << w << " vissza elen nonsat pump " << std::endl;
2.355 - break;
2.356 - } else {
2.357 - /*A saturating push.*/
2.358 -
2.359 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
2.360 - /*v becomes active.*/
2.361 -
2.362 - flow.set(e,0);
2.363 - excess.set(v, excess.get(v)+flow.get(e));
2.364 - excess.set(w, excess.get(w)-flow.get(e));
2.365 - //std::cout << v <<" " << w << " vissza elen sat pump " << std::endl;
2.366 - if (excess.get(w)==0) { break;}
2.367 - } //if (flow.get(e) > excess.get(v))
2.368 - } //if(level.get(w)==level.get(v)+1)
2.369 -
2.370 - else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
2.371 - //std::cout << "Leveldecrease of Node " << w << " to " << newlevel << std::endl;
2.372 -
2.373 - } //if (flow.get(e)>0)
2.374 -
2.375 - } //for in-Edge
2.376 -
2.377 -
2.378 -
2.379 -
2.380 - /*
2.381 - Relabel
2.382 - */
2.383 - if (excess.get(w)>0) {
2.384 - /*Now newlevel <= n*/
2.385 -
2.386 - int l=level.get(w); //l is the old level of w.
2.387 - --numb[l];
2.388 -
2.389 - if (newlevel == n) {
2.390 - level.set(w,n);
2.391 -
2.392 - } else {
2.393 -
2.394 - if (numb[l]) {
2.395 - /*If the level of w remains nonempty.*/
2.396 -
2.397 - level.set(w,++newlevel);
2.398 - ++numb[newlevel];
2.399 - stack[newlevel].push(w);
2.400 - b=newlevel;
2.401 - } else {
2.402 - /*If the level of w gets empty.*/
2.403 -
2.404 - for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
2.405 - if (level.get(v) >= l ) {
2.406 - level.set(v,n);
2.407 - }
2.408 - }
2.409 -
2.410 - for (int i=l+1 ; i!=n ; ++i) numb[i]=0;
2.411 - } //if (numb[l])
2.412 -
2.413 - } // if (newlevel = n)
2.414 -
2.415 - } // if (excess.get(w)>0)
2.416 -
2.417 -
2.418 - } //else
2.419 -
2.420 } //while(b)
2.421
2.422 value=excess.get(t);
2.423 /*Max flow value.*/
2.424
2.425
2.426 -
2.427 - /*
2.428 - We find an empty level, e. The Nodes above this level give
2.429 - a minimum cut.
2.430 + /*
2.431 + We count empty_level. The nodes above this level is a mincut.
2.432 */
2.433 -
2.434 - int e=1;
2.435 -
2.436 - while(e) {
2.437 - if(numb[e]) ++e;
2.438 + while(true) {
2.439 + if(numb[empty_level]) ++empty_level;
2.440 else break;
2.441 }
2.442 - for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) {
2.443 - if (level.get(v) > e) mincutvector.set(v, true);
2.444 - }
2.445
2.446 -
2.447 } // void run()
2.448
2.449
2.450 @@ -295,12 +276,15 @@
2.451
2.452 /*
2.453 Returns a minimum cut.
2.454 - */
2.455 -
2.456 - typename Graph::NodeMap<bool> mincut() {
2.457 - return mincutvector;
2.458 + */
2.459 + template<typename CutMap>
2.460 + void mincut(CutMap& M) {
2.461 +
2.462 + for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) {
2.463 + if ( level.get(v) > empty_level ) M.set(v, true);
2.464 + }
2.465 }
2.466 -
2.467 +
2.468
2.469 };
2.470 }//namespace marci