1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/jacint/preflowproba.h Thu Apr 22 13:51:25 2004 +0000
1.3 @@ -0,0 +1,685 @@
1.4 +// -*- C++ -*-
1.5 +
1.6 +//run gyorsan tudna adni a minmincutot a 2 fazis elejen , ne vegyuk be konstruktorba egy cutmapet?
1.7 +//constzero jo igy?
1.8 +
1.9 +//majd marci megmondja betegyem-e bfs-t meg resgraphot
1.10 +
1.11 +/*
1.12 +Heuristics:
1.13 + 2 phase
1.14 + gap
1.15 + list 'level_list' on the nodes on level i implemented by hand
1.16 + stack 'active' on the active nodes on level i implemented by hand
1.17 + runs heuristic 'highest label' for H1*n relabels
1.18 + runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
1.19 +
1.20 +Parameters H0 and H1 are initialized to 20 and 10.
1.21 +
1.22 +Constructors:
1.23 +
1.24 +Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if
1.25 + FlowMap is not constant zero, and should be true if it is
1.26 +
1.27 +Members:
1.28 +
1.29 +void run()
1.30 +
1.31 +T flowValue() : returns the value of a maximum flow
1.32 +
1.33 +void minMinCut(CutMap& M) : sets M to the characteristic vector of the
1.34 + minimum min cut. M should be a map of bools initialized to false.
1.35 +
1.36 +void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
1.37 + maximum min cut. M should be a map of bools initialized to false.
1.38 +
1.39 +void minCut(CutMap& M) : sets M to the characteristic vector of
1.40 + a min cut. M should be a map of bools initialized to false.
1.41 +
1.42 +FIXME reset
1.43 +
1.44 +*/
1.45 +
1.46 +#ifndef HUGO_PREFLOW_H
1.47 +#define HUGO_PREFLOW_H
1.48 +
1.49 +#define H0 20
1.50 +#define H1 1
1.51 +
1.52 +#include <vector>
1.53 +#include <queue>
1.54 +#include<graph_wrapper.h>
1.55 +
1.56 +namespace hugo {
1.57 +
1.58 + template <typename Graph, typename T,
1.59 + typename CapMap=typename Graph::EdgeMap<T>,
1.60 + typename FlowMap=typename Graph::EdgeMap<T> >
1.61 + class Preflow {
1.62 +
1.63 + typedef typename Graph::Node Node;
1.64 + typedef typename Graph::Edge Edge;
1.65 + typedef typename Graph::NodeIt NodeIt;
1.66 + typedef typename Graph::OutEdgeIt OutEdgeIt;
1.67 + typedef typename Graph::InEdgeIt InEdgeIt;
1.68 +
1.69 + const Graph& G;
1.70 + Node s;
1.71 + Node t;
1.72 + const CapMap& capacity;
1.73 + FlowMap& flow;
1.74 + T value;
1.75 + bool constzero;
1.76 +
1.77 + typedef ResGraphWrapper<const Graph, T, CapMap, FlowMap> ResGW;
1.78 + typedef typename ResGW::OutEdgeIt ResOutEdgeIt;
1.79 + typedef typename ResGW::InEdgeIt ResInEdgeIt;
1.80 + typedef typename ResGW::Edge ResEdge;
1.81 +
1.82 + public:
1.83 + Preflow(Graph& _G, Node _s, Node _t, CapMap& _capacity,
1.84 + FlowMap& _flow, bool _constzero ) :
1.85 + G(_G), s(_s), t(_t), capacity(_capacity), flow(_flow), constzero(_constzero) {}
1.86 +
1.87 +
1.88 + void run() {
1.89 +
1.90 + ResGW res_graph(G, capacity, flow);
1.91 +
1.92 + value=0; //for the subsequent runs
1.93 +
1.94 + bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd
1.95 + int n=G.nodeNum();
1.96 + int heur0=(int)(H0*n); //time while running 'bound decrease'
1.97 + int heur1=(int)(H1*n); //time while running 'highest label'
1.98 + int heur=heur1; //starting time interval (#of relabels)
1.99 + bool what_heur=1;
1.100 + /*
1.101 + what_heur is 0 in case 'bound decrease'
1.102 + and 1 in case 'highest label'
1.103 + */
1.104 + bool end=false;
1.105 + /*
1.106 + Needed for 'bound decrease', 'true'
1.107 + means no active nodes are above bound b.
1.108 + */
1.109 + int relabel=0;
1.110 + int k=n-2; //bound on the highest level under n containing a node
1.111 + int b=k; //bound on the highest level under n of an active node
1.112 +
1.113 + typename Graph::NodeMap<int> level(G,n);
1.114 + typename Graph::NodeMap<T> excess(G);
1.115 +
1.116 + std::vector<Node> active(n-1,INVALID);
1.117 + typename Graph::NodeMap<Node> next(G,INVALID);
1.118 + //Stack of the active nodes in level i < n.
1.119 + //We use it in both phases.
1.120 +
1.121 + typename Graph::NodeMap<Node> left(G,INVALID);
1.122 + typename Graph::NodeMap<Node> right(G,INVALID);
1.123 + std::vector<Node> level_list(n,INVALID);
1.124 + /*
1.125 + List of the nodes in level i<n.
1.126 + */
1.127 +
1.128 +
1.129 + if ( constzero ) {
1.130 +
1.131 + /*Reverse_bfs from t, to find the starting level.*/
1.132 + level.set(t,0);
1.133 + std::queue<Node> bfs_queue;
1.134 + bfs_queue.push(t);
1.135 +
1.136 + while (!bfs_queue.empty()) {
1.137 +
1.138 + Node v=bfs_queue.front();
1.139 + bfs_queue.pop();
1.140 + int l=level[v]+1;
1.141 +
1.142 + InEdgeIt e;
1.143 + for(G.first(e,v); G.valid(e); G.next(e)) {
1.144 + Node w=G.tail(e);
1.145 + if ( level[w] == n && w != s ) {
1.146 + bfs_queue.push(w);
1.147 + Node first=level_list[l];
1.148 + if ( G.valid(first) ) left.set(first,w);
1.149 + right.set(w,first);
1.150 + level_list[l]=w;
1.151 + level.set(w, l);
1.152 + }
1.153 + }
1.154 + }
1.155 +
1.156 + //the starting flow
1.157 + OutEdgeIt e;
1.158 + for(G.first(e,s); G.valid(e); G.next(e))
1.159 + {
1.160 + T c=capacity[e];
1.161 + if ( c == 0 ) continue;
1.162 + Node w=G.head(e);
1.163 + if ( level[w] < n ) {
1.164 + if ( excess[w] == 0 && w!=t ) {
1.165 + next.set(w,active[level[w]]);
1.166 + active[level[w]]=w;
1.167 + }
1.168 + flow.set(e, c);
1.169 + excess.set(w, excess[w]+c);
1.170 + }
1.171 + }
1.172 + }
1.173 + else
1.174 + {
1.175 +
1.176 + /*
1.177 + Reverse_bfs from t in the residual graph,
1.178 + to find the starting level.
1.179 + */
1.180 + level.set(t,0);
1.181 + std::queue<Node> bfs_queue;
1.182 + bfs_queue.push(t);
1.183 +
1.184 + while (!bfs_queue.empty()) {
1.185 +
1.186 + Node v=bfs_queue.front();
1.187 + bfs_queue.pop();
1.188 + int l=level[v]+1;
1.189 +
1.190 + InEdgeIt e;
1.191 + for(G.first(e,v); G.valid(e); G.next(e)) {
1.192 + if ( capacity[e] == flow[e] ) continue;
1.193 + Node w=G.tail(e);
1.194 + if ( level[w] == n && w != s ) {
1.195 + bfs_queue.push(w);
1.196 + Node first=level_list[l];
1.197 + if ( G.valid(first) ) left.set(first,w);
1.198 + right.set(w,first);
1.199 + level_list[l]=w;
1.200 + level.set(w, l);
1.201 + }
1.202 + }
1.203 +
1.204 + OutEdgeIt f;
1.205 + for(G.first(f,v); G.valid(f); G.next(f)) {
1.206 + if ( 0 == flow[f] ) continue;
1.207 + Node w=G.head(f);
1.208 + if ( level[w] == n && w != s ) {
1.209 + bfs_queue.push(w);
1.210 + Node first=level_list[l];
1.211 + if ( G.valid(first) ) left.set(first,w);
1.212 + right.set(w,first);
1.213 + level_list[l]=w;
1.214 + level.set(w, l);
1.215 + }
1.216 + }
1.217 + }
1.218 +
1.219 +
1.220 + /*
1.221 + Counting the excess
1.222 + */
1.223 + NodeIt v;
1.224 + for(G.first(v); G.valid(v); G.next(v)) {
1.225 + T exc=0;
1.226 +
1.227 + InEdgeIt e;
1.228 + for(G.first(e,v); G.valid(e); G.next(e)) exc+=flow[e];
1.229 + OutEdgeIt f;
1.230 + for(G.first(f,v); G.valid(f); G.next(f)) exc-=flow[e];
1.231 +
1.232 + excess.set(v,exc);
1.233 +
1.234 + //putting the active nodes into the stack
1.235 + int lev=level[v];
1.236 + if ( exc > 0 && lev < n ) {
1.237 + next.set(v,active[lev]);
1.238 + active[lev]=v;
1.239 + }
1.240 + }
1.241 +
1.242 +
1.243 + //the starting flow
1.244 + OutEdgeIt e;
1.245 + for(G.first(e,s); G.valid(e); G.next(e))
1.246 + {
1.247 + T rem=capacity[e]-flow[e];
1.248 + if ( rem == 0 ) continue;
1.249 + Node w=G.head(e);
1.250 + if ( level[w] < n ) {
1.251 + if ( excess[w] == 0 && w!=t ) {
1.252 + next.set(w,active[level[w]]);
1.253 + active[level[w]]=w;
1.254 + }
1.255 + flow.set(e, capacity[e]);
1.256 + excess.set(w, excess[w]+rem);
1.257 + }
1.258 + }
1.259 +
1.260 + InEdgeIt f;
1.261 + for(G.first(f,s); G.valid(f); G.next(f))
1.262 + {
1.263 + if ( flow[f] == 0 ) continue;
1.264 + Node w=G.head(f);
1.265 + if ( level[w] < n ) {
1.266 + if ( excess[w] == 0 && w!=t ) {
1.267 + next.set(w,active[level[w]]);
1.268 + active[level[w]]=w;
1.269 + }
1.270 + excess.set(w, excess[w]+flow[f]);
1.271 + flow.set(f, 0);
1.272 + }
1.273 + }
1.274 + }
1.275 +
1.276 +
1.277 +
1.278 +
1.279 + /*
1.280 + End of preprocessing
1.281 + */
1.282 +
1.283 +
1.284 +
1.285 + /*
1.286 + Push/relabel on the highest level active nodes.
1.287 + */
1.288 + while ( true ) {
1.289 +
1.290 + if ( b == 0 ) {
1.291 + if ( phase ) break;
1.292 +
1.293 + if ( !what_heur && !end && k > 0 ) {
1.294 + b=k;
1.295 + end=true;
1.296 + } else {
1.297 + phase=1;
1.298 + level.set(s,0);
1.299 + std::queue<Node> bfs_queue;
1.300 + bfs_queue.push(s);
1.301 +
1.302 + while (!bfs_queue.empty()) {
1.303 +
1.304 + Node v=bfs_queue.front();
1.305 + bfs_queue.pop();
1.306 + int l=level[v]+1;
1.307 +
1.308 + ResInEdgeIt e;
1.309 + for(res_graph.first(e,s); res_graph.valid(e);
1.310 + res_graph.next(e)) {
1.311 + Node u=res_graph.tail(e);
1.312 + if ( level[u] >= n ) {
1.313 + bfs_queue.push(u);
1.314 + level.set(u, l);
1.315 + if ( excess[u] > 0 ) {
1.316 + next.set(u,active[l]);
1.317 + active[l]=u;
1.318 + }
1.319 + }
1.320 + }
1.321 + /* InEdgeIt e;
1.322 + for(G.first(e,v); G.valid(e); G.next(e)) {
1.323 + if ( capacity[e] == flow[e] ) continue;
1.324 + Node u=G.tail(e);
1.325 + if ( level[u] >= n ) {
1.326 + bfs_queue.push(u);
1.327 + level.set(u, l);
1.328 + if ( excess[u] > 0 ) {
1.329 + next.set(u,active[l]);
1.330 + active[l]=u;
1.331 + }
1.332 + }
1.333 + }
1.334 +
1.335 + OutEdgeIt f;
1.336 + for(G.first(f,v); G.valid(f); G.next(f)) {
1.337 + if ( 0 == flow[f] ) continue;
1.338 + Node u=G.head(f);
1.339 + if ( level[u] >= n ) {
1.340 + bfs_queue.push(u);
1.341 + level.set(u, l);
1.342 + if ( excess[u] > 0 ) {
1.343 + next.set(u,active[l]);
1.344 + active[l]=u;
1.345 + }
1.346 + }
1.347 + }*/
1.348 + }
1.349 + b=n-2;
1.350 + }
1.351 +
1.352 + }
1.353 +
1.354 +
1.355 + if ( !G.valid(active[b]) ) --b;
1.356 + else {
1.357 + end=false;
1.358 +
1.359 + Node w=active[b];
1.360 + active[b]=next[w];
1.361 + int lev=level[w];
1.362 + T exc=excess[w];
1.363 + int newlevel=n; //bound on the next level of w
1.364 +
1.365 + OutEdgeIt e;
1.366 + for(G.first(e,w); G.valid(e); G.next(e)) {
1.367 +
1.368 + if ( flow[e] == capacity[e] ) continue;
1.369 + Node v=G.head(e);
1.370 + //e=wv
1.371 +
1.372 + if( lev > level[v] ) {
1.373 + /*Push is allowed now*/
1.374 +
1.375 + if ( excess[v]==0 && v!=t && v!=s ) {
1.376 + int lev_v=level[v];
1.377 + next.set(v,active[lev_v]);
1.378 + active[lev_v]=v;
1.379 + }
1.380 +
1.381 + T cap=capacity[e];
1.382 + T flo=flow[e];
1.383 + T remcap=cap-flo;
1.384 +
1.385 + if ( remcap >= exc ) {
1.386 + /*A nonsaturating push.*/
1.387 +
1.388 + flow.set(e, flo+exc);
1.389 + excess.set(v, excess[v]+exc);
1.390 + exc=0;
1.391 + break;
1.392 +
1.393 + } else {
1.394 + /*A saturating push.*/
1.395 +
1.396 + flow.set(e, cap);
1.397 + excess.set(v, excess[v]+remcap);
1.398 + exc-=remcap;
1.399 + }
1.400 + } else if ( newlevel > level[v] ){
1.401 + newlevel = level[v];
1.402 + }
1.403 +
1.404 + } //for out edges wv
1.405 +
1.406 +
1.407 + if ( exc > 0 ) {
1.408 + InEdgeIt e;
1.409 + for(G.first(e,w); G.valid(e); G.next(e)) {
1.410 +
1.411 + if( flow[e] == 0 ) continue;
1.412 + Node v=G.tail(e);
1.413 + //e=vw
1.414 +
1.415 + if( lev > level[v] ) {
1.416 + /*Push is allowed now*/
1.417 +
1.418 + if ( excess[v]==0 && v!=t && v!=s ) {
1.419 + int lev_v=level[v];
1.420 + next.set(v,active[lev_v]);
1.421 + active[lev_v]=v;
1.422 + }
1.423 +
1.424 + T flo=flow[e];
1.425 +
1.426 + if ( flo >= exc ) {
1.427 + /*A nonsaturating push.*/
1.428 +
1.429 + flow.set(e, flo-exc);
1.430 + excess.set(v, excess[v]+exc);
1.431 + exc=0;
1.432 + break;
1.433 + } else {
1.434 + /*A saturating push.*/
1.435 +
1.436 + excess.set(v, excess[v]+flo);
1.437 + exc-=flo;
1.438 + flow.set(e,0);
1.439 + }
1.440 + } else if ( newlevel > level[v] ) {
1.441 + newlevel = level[v];
1.442 + }
1.443 + } //for in edges vw
1.444 +
1.445 + } // if w still has excess after the out edge for cycle
1.446 +
1.447 + excess.set(w, exc);
1.448 +
1.449 + /*
1.450 + Relabel
1.451 + */
1.452 +
1.453 +
1.454 + if ( exc > 0 ) {
1.455 + //now 'lev' is the old level of w
1.456 +
1.457 + if ( phase ) {
1.458 + level.set(w,++newlevel);
1.459 + next.set(w,active[newlevel]);
1.460 + active[newlevel]=w;
1.461 + b=newlevel;
1.462 + } else {
1.463 + //unlacing starts
1.464 + Node right_n=right[w];
1.465 + Node left_n=left[w];
1.466 +
1.467 + if ( G.valid(right_n) ) {
1.468 + if ( G.valid(left_n) ) {
1.469 + right.set(left_n, right_n);
1.470 + left.set(right_n, left_n);
1.471 + } else {
1.472 + level_list[lev]=right_n;
1.473 + left.set(right_n, INVALID);
1.474 + }
1.475 + } else {
1.476 + if ( G.valid(left_n) ) {
1.477 + right.set(left_n, INVALID);
1.478 + } else {
1.479 + level_list[lev]=INVALID;
1.480 + }
1.481 + }
1.482 + //unlacing ends
1.483 +
1.484 + if ( !G.valid(level_list[lev]) ) {
1.485 +
1.486 + //gapping starts
1.487 + for (int i=lev; i!=k ; ) {
1.488 + Node v=level_list[++i];
1.489 + while ( G.valid(v) ) {
1.490 + level.set(v,n);
1.491 + v=right[v];
1.492 + }
1.493 + level_list[i]=INVALID;
1.494 + if ( !what_heur ) active[i]=INVALID;
1.495 + }
1.496 +
1.497 + level.set(w,n);
1.498 + b=lev-1;
1.499 + k=b;
1.500 + //gapping ends
1.501 +
1.502 + } else {
1.503 +
1.504 + if ( newlevel == n ) level.set(w,n);
1.505 + else {
1.506 + level.set(w,++newlevel);
1.507 + next.set(w,active[newlevel]);
1.508 + active[newlevel]=w;
1.509 + if ( what_heur ) b=newlevel;
1.510 + if ( k < newlevel ) ++k; //now k=newlevel
1.511 + Node first=level_list[newlevel];
1.512 + if ( G.valid(first) ) left.set(first,w);
1.513 + right.set(w,first);
1.514 + left.set(w,INVALID);
1.515 + level_list[newlevel]=w;
1.516 + }
1.517 + }
1.518 +
1.519 +
1.520 + ++relabel;
1.521 + if ( relabel >= heur ) {
1.522 + relabel=0;
1.523 + if ( what_heur ) {
1.524 + what_heur=0;
1.525 + heur=heur0;
1.526 + end=false;
1.527 + } else {
1.528 + what_heur=1;
1.529 + heur=heur1;
1.530 + b=k;
1.531 + }
1.532 + }
1.533 + } //phase 0
1.534 +
1.535 +
1.536 + } // if ( exc > 0 )
1.537 +
1.538 +
1.539 + } // if stack[b] is nonempty
1.540 +
1.541 + } // while(true)
1.542 +
1.543 +
1.544 + value = excess[t];
1.545 + /*Max flow value.*/
1.546 +
1.547 + } //void run()
1.548 +
1.549 +
1.550 +
1.551 +
1.552 +
1.553 + /*
1.554 + Returns the maximum value of a flow.
1.555 + */
1.556 +
1.557 + T flowValue() {
1.558 + return value;
1.559 + }
1.560 +
1.561 +
1.562 + FlowMap Flow() {
1.563 + return flow;
1.564 + }
1.565 +
1.566 +
1.567 +
1.568 + void Flow(FlowMap& _flow ) {
1.569 + NodeIt v;
1.570 + for(G.first(v) ; G.valid(v); G.next(v))
1.571 + _flow.set(v,flow[v]);
1.572 + }
1.573 +
1.574 +
1.575 +
1.576 + /*
1.577 + Returns the minimum min cut, by a bfs from s in the residual graph.
1.578 + */
1.579 +
1.580 + template<typename _CutMap>
1.581 + void minMinCut(_CutMap& M) {
1.582 +
1.583 + std::queue<Node> queue;
1.584 +
1.585 + M.set(s,true);
1.586 + queue.push(s);
1.587 +
1.588 + while (!queue.empty()) {
1.589 + Node w=queue.front();
1.590 + queue.pop();
1.591 +
1.592 + OutEdgeIt e;
1.593 + for(G.first(e,w) ; G.valid(e); G.next(e)) {
1.594 + Node v=G.head(e);
1.595 + if (!M[v] && flow[e] < capacity[e] ) {
1.596 + queue.push(v);
1.597 + M.set(v, true);
1.598 + }
1.599 + }
1.600 +
1.601 + InEdgeIt f;
1.602 + for(G.first(f,w) ; G.valid(f); G.next(f)) {
1.603 + Node v=G.tail(f);
1.604 + if (!M[v] && flow[f] > 0 ) {
1.605 + queue.push(v);
1.606 + M.set(v, true);
1.607 + }
1.608 + }
1.609 + }
1.610 + }
1.611 +
1.612 +
1.613 +
1.614 + /*
1.615 + Returns the maximum min cut, by a reverse bfs
1.616 + from t in the residual graph.
1.617 + */
1.618 +
1.619 + template<typename _CutMap>
1.620 + void maxMinCut(_CutMap& M) {
1.621 +
1.622 + std::queue<Node> queue;
1.623 +
1.624 + M.set(t,true);
1.625 + queue.push(t);
1.626 +
1.627 + while (!queue.empty()) {
1.628 + Node w=queue.front();
1.629 + queue.pop();
1.630 +
1.631 +
1.632 + InEdgeIt e;
1.633 + for(G.first(e,w) ; G.valid(e); G.next(e)) {
1.634 + Node v=G.tail(e);
1.635 + if (!M[v] && flow[e] < capacity[e] ) {
1.636 + queue.push(v);
1.637 + M.set(v, true);
1.638 + }
1.639 + }
1.640 +
1.641 + OutEdgeIt f;
1.642 + for(G.first(f,w) ; G.valid(f); G.next(f)) {
1.643 + Node v=G.head(f);
1.644 + if (!M[v] && flow[f] > 0 ) {
1.645 + queue.push(v);
1.646 + M.set(v, true);
1.647 + }
1.648 + }
1.649 + }
1.650 +
1.651 + NodeIt v;
1.652 + for(G.first(v) ; G.valid(v); G.next(v)) {
1.653 + M.set(v, !M[v]);
1.654 + }
1.655 +
1.656 + }
1.657 +
1.658 +
1.659 +
1.660 + template<typename CutMap>
1.661 + void minCut(CutMap& M) {
1.662 + minMinCut(M);
1.663 + }
1.664 +
1.665 +
1.666 + void reset_target (Node _t) {t=_t;}
1.667 + void reset_source (Node _s) {s=_s;}
1.668 +
1.669 + template<typename _CapMap>
1.670 + void reset_cap (_CapMap _cap) {capacity=_cap;}
1.671 +
1.672 + template<typename _FlowMap>
1.673 + void reset_cap (_FlowMap _flow, bool _constzero) {
1.674 + flow=_flow;
1.675 + constzero=_constzero;
1.676 + }
1.677 +
1.678 +
1.679 +
1.680 + };
1.681 +
1.682 +} //namespace hugo
1.683 +
1.684 +#endif //PREFLOW_H
1.685 +
1.686 +
1.687 +
1.688 +