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