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