1.1 --- a/src/work/jacint/preflow_push_max_flow.h Wed Feb 18 13:06:41 2004 +0000
1.2 +++ b/src/work/jacint/preflow_push_max_flow.h Wed Feb 18 14:42:38 2004 +0000
1.3 @@ -1,3 +1,4 @@
1.4 +// -*- C++ -*-
1.5 /*
1.6 preflow_push_max_flow_h
1.7 by jacint.
1.8 @@ -15,13 +16,16 @@
1.9
1.10 T maxflow() : returns the value of a maximum flow
1.11
1.12 -NodeMap<bool> mincut(): returns a
1.13 - characteristic vector of a minimum cut.
1.14 +void mincut(CutMap& M) : sets M to the characteristic vector of a
1.15 + minimum cut. M should be a map of bools initialized to false.
1.16 +
1.17 */
1.18
1.19 #ifndef PREFLOW_PUSH_MAX_FLOW_H
1.20 #define PREFLOW_PUSH_MAX_FLOW_H
1.21
1.22 +#define A 1
1.23 +
1.24 #include <algorithm>
1.25 #include <vector>
1.26 #include <stack>
1.27 @@ -31,7 +35,9 @@
1.28
1.29 namespace marci {
1.30
1.31 - template <typename Graph, typename T>
1.32 + template <typename Graph, typename T,
1.33 + typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>,
1.34 + typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
1.35 class preflow_push_max_flow {
1.36
1.37 typedef typename Graph::NodeIt NodeIt;
1.38 @@ -42,17 +48,15 @@
1.39 Graph& G;
1.40 NodeIt s;
1.41 NodeIt t;
1.42 - typename Graph::EdgeMap<T>& capacity;
1.43 - T value;
1.44 - typename Graph::NodeMap<bool> mincutvector;
1.45 -
1.46 -
1.47 -
1.48 + IntMap level;
1.49 + CapMap& capacity;
1.50 + int empty_level; //an empty level in the end of run()
1.51 + T value;
1.52 +
1.53 public:
1.54 -
1.55 - preflow_push_max_flow ( Graph& _G, NodeIt _s, NodeIt _t,
1.56 - typename Graph::EdgeMap<T>& _capacity) :
1.57 - G(_G), s(_s), t(_t), capacity(_capacity), mincutvector(_G, false) { }
1.58 +
1.59 + preflow_push_max_flow(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
1.60 + G(_G), s(_s), t(_t), level(_G), capacity(_capacity) { }
1.61
1.62
1.63 /*
1.64 @@ -62,223 +66,200 @@
1.65 */
1.66 void run() {
1.67
1.68 - typename Graph::EdgeMap<T> flow(G, 0);
1.69 - typename Graph::NodeMap<int> level(G);
1.70 - typename Graph::NodeMap<T> excess(G);
1.71 -
1.72 - int n=G.nodeNum();
1.73 + int n=G.nodeNum();
1.74 int b=n-2;
1.75 /*
1.76 - b is a bound on the highest level of an active Node.
1.77 - In the beginning it is at most n-2.
1.78 + b is a bound on the highest level of an active node.
1.79 */
1.80 -
1.81 - std::vector<int> numb(n); //The number of Nodes on level i < n.
1.82 - std::vector<std::stack<NodeIt> > stack(2*n-1);
1.83 - //Stack of the active Nodes in level i.
1.84 +
1.85 + IntMap level(G,n);
1.86 + TMap excess(G);
1.87 + FlowMap flow(G,0);
1.88 +
1.89 + std::vector<int> numb(n);
1.90 + /*
1.91 + The number of nodes on level i < n. It is
1.92 + initialized to n+1, because of the reverse_bfs-part.
1.93 + */
1.94 +
1.95 + std::vector<std::stack<NodeIt> > stack(n);
1.96 + //Stack of the active nodes in level i.
1.97 +
1.98
1.99 /*Reverse_bfs from t, to find the starting level.*/
1.100 - reverse_bfs<Graph> bfs(G, t);
1.101 - bfs.run();
1.102 - for(EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v)
1.103 - {
1.104 - int dist=bfs.dist(v);
1.105 - level.set(v, dist);
1.106 - ++numb[dist];
1.107 + level.set(t,0);
1.108 + std::queue<NodeIt> bfs_queue;
1.109 + bfs_queue.push(t);
1.110 +
1.111 + while (!bfs_queue.empty()) {
1.112 +
1.113 + NodeIt v=bfs_queue.front();
1.114 + bfs_queue.pop();
1.115 + int l=level.get(v)+1;
1.116 +
1.117 + for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
1.118 + NodeIt w=G.tail(e);
1.119 + if ( level.get(w) == n ) {
1.120 + bfs_queue.push(w);
1.121 + ++numb[l];
1.122 + level.set(w, l);
1.123 + }
1.124 }
1.125 -
1.126 + }
1.127 +
1.128 level.set(s,n);
1.129
1.130 - /* Starting flow. It is everywhere 0 at the moment. */
1.131 +
1.132 +
1.133 + /* Starting flow. It is everywhere 0 at the moment. */
1.134 for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e)
1.135 {
1.136 - if ( capacity.get(e) > 0 ) {
1.137 - NodeIt w=G.head(e);
1.138 + if ( capacity.get(e) == 0 ) continue;
1.139 + NodeIt w=G.head(e);
1.140 + if ( level.get(w) < n ) {
1.141 + if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w);
1.142 flow.set(e, capacity.get(e));
1.143 - stack[level.get(w)].push(w);
1.144 excess.set(w, excess.get(w)+capacity.get(e));
1.145 }
1.146 }
1.147 -
1.148 +
1.149 /*
1.150 End of preprocessing
1.151 */
1.152
1.153
1.154 + /*
1.155 + Push/relabel on the highest level active nodes.
1.156 + */
1.157 + /*While there exists an active node.*/
1.158 + while (b) {
1.159 + if ( stack[b].empty() ) {
1.160 + --b;
1.161 + continue;
1.162 + }
1.163 +
1.164 + NodeIt w=stack[b].top(); //w is a highest label active node.
1.165 + stack[b].pop();
1.166 + int lev=level.get(w);
1.167 + int exc=excess.get(w);
1.168 + int newlevel=2*n-2; //In newlevel we bound the next level of w.
1.169 +
1.170 + // if ( level.get(w) < n ) { //Nem tudom ez mukodik-e
1.171 + for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
1.172 +
1.173 + if ( flow.get(e) == capacity.get(e) ) continue;
1.174 + NodeIt v=G.head(e);
1.175 + //e=wv
1.176 +
1.177 + if( lev > level.get(v) ) {
1.178 + /*Push is allowed now*/
1.179 +
1.180 + if ( excess.get(v)==0 && v != s && v !=t )
1.181 + stack[level.get(v)].push(v);
1.182 + /*v becomes active.*/
1.183 +
1.184 + int cap=capacity.get(e);
1.185 + int flo=flow.get(e);
1.186 + int remcap=cap-flo;
1.187 +
1.188 + if ( remcap >= exc ) {
1.189 + /*A nonsaturating push.*/
1.190 +
1.191 + flow.set(e, flo+exc);
1.192 + excess.set(v, excess.get(v)+exc);
1.193 + exc=0;
1.194 + break;
1.195 +
1.196 + } else {
1.197 + /*A saturating push.*/
1.198 +
1.199 + flow.set(e, cap );
1.200 + excess.set(v, excess.get(v)+remcap);
1.201 + exc-=remcap;
1.202 + }
1.203 + } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
1.204 +
1.205 + } //for out edges wv
1.206 +
1.207 +
1.208 + if ( exc > 0 ) {
1.209 + for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
1.210 +
1.211 + if( flow.get(e) == 0 ) continue;
1.212 + NodeIt v=G.tail(e);
1.213 + //e=vw
1.214 +
1.215 + if( lev > level.get(v) ) {
1.216 + /*Push is allowed now*/
1.217 +
1.218 + if ( excess.get(v)==0 && v != s && v !=t)
1.219 + stack[level.get(v)].push(v);
1.220 + /*v becomes active.*/
1.221 +
1.222 + int flo=flow.get(e);
1.223 +
1.224 + if ( flo >= exc ) {
1.225 + /*A nonsaturating push.*/
1.226 +
1.227 + flow.set(e, flo-exc);
1.228 + excess.set(v, excess.get(v)+exc);
1.229 + exc=0;
1.230 + break;
1.231 + } else {
1.232 + /*A saturating push.*/
1.233 +
1.234 + excess.set(v, excess.get(v)+flo);
1.235 + exc-=flo;
1.236 + flow.set(e,0);
1.237 + }
1.238 + } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
1.239 +
1.240 + } //for in edges vw
1.241 +
1.242 + } // if w still has excess after the out edge for cycle
1.243 +
1.244 + excess.set(w, exc);
1.245 +
1.246 +
1.247 + /*
1.248 + Relabel
1.249 + */
1.250 +
1.251 + if ( exc > 0 ) {
1.252 + //now 'lev' is the old level of w
1.253 + level.set(w,++newlevel);
1.254 + --numb[lev];
1.255 +
1.256 + if ( !numb[lev] && lev < A*n ) { //If the level of w gets empty.
1.257 +
1.258 + for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
1.259 + if (level.get(v) > lev ) level.set(v,n);
1.260 + }
1.261 + for (int i=lev+1 ; i!=n ; ++i) numb[i]=0;
1.262 + if ( newlevel < n ) newlevel=n;
1.263 + } else if ( newlevel < n ) {
1.264 + ++numb[newlevel];
1.265 + stack[newlevel].push(w);
1.266 + b=newlevel;
1.267 + }
1.268 + }
1.269
1.270 - /*
1.271 - Push/relabel on the highest level active Nodes.
1.272 - */
1.273 -
1.274 - /*While there exists an active Node.*/
1.275 - while (b) {
1.276
1.277 - /*We decrease the bound if there is no active node of level b.*/
1.278 - if (stack[b].empty()) {
1.279 - --b;
1.280 - } else {
1.281
1.282 - NodeIt w=stack[b].top(); //w is the highest label active Node.
1.283 - stack[b].pop(); //We delete w from the stack.
1.284 -
1.285 - int newlevel=2*n-2; //In newlevel we maintain the next level of w.
1.286 -
1.287 - for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
1.288 - NodeIt v=G.head(e);
1.289 - /*e is the Edge wv.*/
1.290 -
1.291 - if (flow.get(e)<capacity.get(e)) {
1.292 - /*e is an Edge of the residual graph */
1.293 -
1.294 - if(level.get(w)==level.get(v)+1) {
1.295 - /*Push is allowed now*/
1.296 -
1.297 - if (capacity.get(e)-flow.get(e) > excess.get(w)) {
1.298 - /*A nonsaturating push.*/
1.299 -
1.300 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
1.301 - /*v becomes active.*/
1.302 -
1.303 - flow.set(e, flow.get(e)+excess.get(w));
1.304 - excess.set(v, excess.get(v)+excess.get(w));
1.305 - excess.set(w,0);
1.306 - //std::cout << w << " " << v <<" elore elen nonsat pump " << std::endl;
1.307 - break;
1.308 - } else {
1.309 - /*A saturating push.*/
1.310 -
1.311 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
1.312 - /*v becomes active.*/
1.313 -
1.314 - excess.set(v, excess.get(v)+capacity.get(e)-flow.get(e));
1.315 - excess.set(w, excess.get(w)-capacity.get(e)+flow.get(e));
1.316 - flow.set(e, capacity.get(e));
1.317 - //std::cout << w <<" " << v <<" elore elen sat pump " << std::endl;
1.318 - if (excess.get(w)==0) break;
1.319 - /*If w is not active any more, then we go on to the next Node.*/
1.320 -
1.321 - } // if (capacity.get(e)-flow.get(e) > excess.get(w))
1.322 - } // if (level.get(w)==level.get(v)+1)
1.323 -
1.324 - else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
1.325 -
1.326 - } //if (flow.get(e)<capacity.get(e))
1.327 -
1.328 - } //for(OutEdgeIt e=G.first_OutEdge(w); e.valid(); ++e)
1.329 -
1.330 -
1.331 -
1.332 - for(InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
1.333 - NodeIt v=G.tail(e);
1.334 - /*e is the Edge vw.*/
1.335 -
1.336 - if (excess.get(w)==0) break;
1.337 - /*It may happen, that w became inactive in the first 'for' cycle.*/
1.338 -
1.339 - if(flow.get(e)>0) {
1.340 - /*e is an Edge of the residual graph */
1.341 -
1.342 - if(level.get(w)==level.get(v)+1) {
1.343 - /*Push is allowed now*/
1.344 -
1.345 - if (flow.get(e) > excess.get(w)) {
1.346 - /*A nonsaturating push.*/
1.347 -
1.348 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
1.349 - /*v becomes active.*/
1.350 -
1.351 - flow.set(e, flow.get(e)-excess.get(w));
1.352 - excess.set(v, excess.get(v)+excess.get(w));
1.353 - excess.set(w,0);
1.354 - //std::cout << v << " " << w << " vissza elen nonsat pump " << std::endl;
1.355 - break;
1.356 - } else {
1.357 - /*A saturating push.*/
1.358 -
1.359 - if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v);
1.360 - /*v becomes active.*/
1.361 -
1.362 - flow.set(e,0);
1.363 - excess.set(v, excess.get(v)+flow.get(e));
1.364 - excess.set(w, excess.get(w)-flow.get(e));
1.365 - //std::cout << v <<" " << w << " vissza elen sat pump " << std::endl;
1.366 - if (excess.get(w)==0) { break;}
1.367 - } //if (flow.get(e) > excess.get(v))
1.368 - } //if(level.get(w)==level.get(v)+1)
1.369 -
1.370 - else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
1.371 - //std::cout << "Leveldecrease of Node " << w << " to " << newlevel << std::endl;
1.372 -
1.373 - } //if (flow.get(e)>0)
1.374 -
1.375 - } //for in-Edge
1.376 -
1.377 -
1.378 -
1.379 -
1.380 - /*
1.381 - Relabel
1.382 - */
1.383 - if (excess.get(w)>0) {
1.384 - /*Now newlevel <= n*/
1.385 -
1.386 - int l=level.get(w); //l is the old level of w.
1.387 - --numb[l];
1.388 -
1.389 - if (newlevel == n) {
1.390 - level.set(w,n);
1.391 -
1.392 - } else {
1.393 -
1.394 - if (numb[l]) {
1.395 - /*If the level of w remains nonempty.*/
1.396 -
1.397 - level.set(w,++newlevel);
1.398 - ++numb[newlevel];
1.399 - stack[newlevel].push(w);
1.400 - b=newlevel;
1.401 - } else {
1.402 - /*If the level of w gets empty.*/
1.403 -
1.404 - for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
1.405 - if (level.get(v) >= l ) {
1.406 - level.set(v,n);
1.407 - }
1.408 - }
1.409 -
1.410 - for (int i=l+1 ; i!=n ; ++i) numb[i]=0;
1.411 - } //if (numb[l])
1.412 -
1.413 - } // if (newlevel = n)
1.414 -
1.415 - } // if (excess.get(w)>0)
1.416 -
1.417 -
1.418 - } //else
1.419 -
1.420 } //while(b)
1.421
1.422 value=excess.get(t);
1.423 /*Max flow value.*/
1.424
1.425
1.426 -
1.427 - /*
1.428 - We find an empty level, e. The Nodes above this level give
1.429 - a minimum cut.
1.430 + /*
1.431 + We count empty_level. The nodes above this level is a mincut.
1.432 */
1.433 -
1.434 - int e=1;
1.435 -
1.436 - while(e) {
1.437 - if(numb[e]) ++e;
1.438 + while(true) {
1.439 + if(numb[empty_level]) ++empty_level;
1.440 else break;
1.441 }
1.442 - for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) {
1.443 - if (level.get(v) > e) mincutvector.set(v, true);
1.444 - }
1.445
1.446 -
1.447 } // void run()
1.448
1.449
1.450 @@ -295,12 +276,15 @@
1.451
1.452 /*
1.453 Returns a minimum cut.
1.454 - */
1.455 -
1.456 - typename Graph::NodeMap<bool> mincut() {
1.457 - return mincutvector;
1.458 + */
1.459 + template<typename CutMap>
1.460 + void mincut(CutMap& M) {
1.461 +
1.462 + for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) {
1.463 + if ( level.get(v) > empty_level ) M.set(v, true);
1.464 + }
1.465 }
1.466 -
1.467 +
1.468
1.469 };
1.470 }//namespace marci