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