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