1.1 --- a/src/work/jacint/preflow.h Thu Apr 29 16:26:01 2004 +0000
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,1016 +0,0 @@
1.4 -// -*- C++ -*-
1.5 -
1.6 -/*
1.7 -Heuristics:
1.8 - 2 phase
1.9 - gap
1.10 - list 'level_list' on the nodes on level i implemented by hand
1.11 - stack 'active' on the active nodes on level i
1.12 - runs heuristic 'highest label' for H1*n relabels
1.13 - runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
1.14 -
1.15 -Parameters H0 and H1 are initialized to 20 and 1.
1.16 -
1.17 -Constructors:
1.18 -
1.19 -Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if
1.20 - FlowMap is not constant zero, and should be true if it is
1.21 -
1.22 -Members:
1.23 -
1.24 -void run()
1.25 -
1.26 -Num flowValue() : returns the value of a maximum flow
1.27 -
1.28 -void minMinCut(CutMap& M) : sets M to the characteristic vector of the
1.29 - minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
1.30 -
1.31 -void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
1.32 - maximum min cut. M should be a map of bools initialized to false.
1.33 -
1.34 -void minCut(CutMap& M) : sets M to the characteristic vector of
1.35 - a min cut. M should be a map of bools initialized to false.
1.36 -
1.37 -*/
1.38 -
1.39 -#ifndef HUGO_PREFLOW_H
1.40 -#define HUGO_PREFLOW_H
1.41 -
1.42 -#define H0 20
1.43 -#define H1 1
1.44 -
1.45 -#include <vector>
1.46 -#include <queue>
1.47 -#include <stack>
1.48 -
1.49 -#include <graph_wrapper.h>
1.50 -#include <bfs_iterator.h>
1.51 -#include <invalid.h>
1.52 -#include <maps.h>
1.53 -#include <for_each_macros.h>
1.54 -
1.55 -
1.56 -namespace hugo {
1.57 -
1.58 - template <typename Graph, typename Num,
1.59 - typename CapMap=typename Graph::template EdgeMap<Num>,
1.60 - typename FlowMap=typename Graph::template EdgeMap<Num> >
1.61 - class MaxFlow {
1.62 -
1.63 - typedef typename Graph::Node Node;
1.64 - typedef typename Graph::NodeIt NodeIt;
1.65 - typedef typename Graph::OutEdgeIt OutEdgeIt;
1.66 - typedef typename Graph::InEdgeIt InEdgeIt;
1.67 -
1.68 - typedef typename std::vector<std::stack<Node> > VecStack;
1.69 - typedef typename Graph::template NodeMap<Node> NNMap;
1.70 - typedef typename std::vector<Node> VecNode;
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 - int n; //the number of nodes of G
1.78 - typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
1.79 - typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
1.80 - typedef typename ResGW::Edge ResGWEdge;
1.81 - //typedef typename ResGW::template NodeMap<bool> ReachedMap;
1.82 - typedef typename Graph::template NodeMap<int> ReachedMap;
1.83 - ReachedMap level;
1.84 - //level works as a bool map in augmenting path algorithms
1.85 - //and is used by bfs for storing reached information.
1.86 - //In preflow, it shows levels of nodes.
1.87 - //typename Graph::template NodeMap<int> level;
1.88 - typename Graph::template NodeMap<Num> excess;
1.89 -
1.90 - public:
1.91 -
1.92 - enum flowEnum{
1.93 - ZERO_FLOW=0,
1.94 - GEN_FLOW=1,
1.95 - PREFLOW=2
1.96 - };
1.97 -
1.98 - MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
1.99 - FlowMap& _flow) :
1.100 - g(&_G), s(_s), t(_t), capacity(&_capacity),
1.101 - flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
1.102 -
1.103 - void run() {
1.104 - preflow( ZERO_FLOW );
1.105 - }
1.106 -
1.107 - void preflow( flowEnum fe ) {
1.108 - preflowPhase0(fe);
1.109 - preflowPhase1();
1.110 - }
1.111 -
1.112 - void preflowPhase0( flowEnum fe );
1.113 -
1.114 - void preflowPhase1();
1.115 -
1.116 - bool augmentOnShortestPath();
1.117 -
1.118 - template<typename MutableGraph> bool augmentOnBlockingFlow();
1.119 -
1.120 - bool augmentOnBlockingFlow2();
1.121 -
1.122 - /// Returns the actual flow value.
1.123 - /// More precisely, it returns the negative excess of s, thus
1.124 - /// this works also for preflows.
1.125 - Num flowValue() {
1.126 - Num a=0;
1.127 - FOR_EACH_INC_LOC(OutEdgeIt, e, *g, s) a+=(*flow)[e];
1.128 - FOR_EACH_INC_LOC(InEdgeIt, e, *g, s) a-=(*flow)[e];
1.129 - return a;
1.130 - }
1.131 -
1.132 - //should be used only between preflowPhase0 and preflowPhase1
1.133 - template<typename _CutMap>
1.134 - void actMinCut(_CutMap& M) {
1.135 - NodeIt v;
1.136 - for(g->first(v); g->valid(v); g->next(v))
1.137 - if ( level[v] < n ) {
1.138 - M.set(v,false);
1.139 - } else {
1.140 - M.set(v,true);
1.141 - }
1.142 - }
1.143 -
1.144 -
1.145 -
1.146 - /*
1.147 - Returns the minimum min cut, by a bfs from s in the residual graph.
1.148 - */
1.149 - template<typename _CutMap>
1.150 - void minMinCut(_CutMap& M) {
1.151 -
1.152 - std::queue<Node> queue;
1.153 -
1.154 - M.set(s,true);
1.155 - queue.push(s);
1.156 -
1.157 - while (!queue.empty()) {
1.158 - Node w=queue.front();
1.159 - queue.pop();
1.160 -
1.161 - OutEdgeIt e;
1.162 - for(g->first(e,w) ; g->valid(e); g->next(e)) {
1.163 - Node v=g->head(e);
1.164 - if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
1.165 - queue.push(v);
1.166 - M.set(v, true);
1.167 - }
1.168 - }
1.169 -
1.170 - InEdgeIt f;
1.171 - for(g->first(f,w) ; g->valid(f); g->next(f)) {
1.172 - Node v=g->tail(f);
1.173 - if (!M[v] && (*flow)[f] > 0 ) {
1.174 - queue.push(v);
1.175 - M.set(v, true);
1.176 - }
1.177 - }
1.178 - }
1.179 - }
1.180 -
1.181 -
1.182 -
1.183 - /*
1.184 - Returns the maximum min cut, by a reverse bfs
1.185 - from t in the residual graph.
1.186 - */
1.187 -
1.188 - template<typename _CutMap>
1.189 - void maxMinCut(_CutMap& M) {
1.190 -
1.191 - NodeIt v;
1.192 - for(g->first(v) ; g->valid(v); g->next(v)) {
1.193 - M.set(v, true);
1.194 - }
1.195 -
1.196 - std::queue<Node> queue;
1.197 -
1.198 - M.set(t,false);
1.199 - queue.push(t);
1.200 -
1.201 - while (!queue.empty()) {
1.202 - Node w=queue.front();
1.203 - queue.pop();
1.204 -
1.205 -
1.206 - InEdgeIt e;
1.207 - for(g->first(e,w) ; g->valid(e); g->next(e)) {
1.208 - Node v=g->tail(e);
1.209 - if (M[v] && (*flow)[e] < (*capacity)[e] ) {
1.210 - queue.push(v);
1.211 - M.set(v, false);
1.212 - }
1.213 - }
1.214 -
1.215 - OutEdgeIt f;
1.216 - for(g->first(f,w) ; g->valid(f); g->next(f)) {
1.217 - Node v=g->head(f);
1.218 - if (M[v] && (*flow)[f] > 0 ) {
1.219 - queue.push(v);
1.220 - M.set(v, false);
1.221 - }
1.222 - }
1.223 - }
1.224 - }
1.225 -
1.226 -
1.227 - template<typename CutMap>
1.228 - void minCut(CutMap& M) {
1.229 - minMinCut(M);
1.230 - }
1.231 -
1.232 - void resetTarget(Node _t) {t=_t;}
1.233 - void resetSource(Node _s) {s=_s;}
1.234 -
1.235 - void resetCap(const CapMap& _cap) {
1.236 - capacity=&_cap;
1.237 - }
1.238 -
1.239 - void resetFlow(FlowMap& _flow) {
1.240 - flow=&_flow;
1.241 - }
1.242 -
1.243 -
1.244 - private:
1.245 -
1.246 - int push(Node w, VecStack& active) {
1.247 -
1.248 - int lev=level[w];
1.249 - Num exc=excess[w];
1.250 - int newlevel=n; //bound on the next level of w
1.251 -
1.252 - OutEdgeIt e;
1.253 - for(g->first(e,w); g->valid(e); g->next(e)) {
1.254 -
1.255 - if ( (*flow)[e] >= (*capacity)[e] ) continue;
1.256 - Node v=g->head(e);
1.257 -
1.258 - if( lev > level[v] ) { //Push is allowed now
1.259 -
1.260 - if ( excess[v]<=0 && v!=t && v!=s ) {
1.261 - int lev_v=level[v];
1.262 - active[lev_v].push(v);
1.263 - }
1.264 -
1.265 - Num cap=(*capacity)[e];
1.266 - Num flo=(*flow)[e];
1.267 - Num remcap=cap-flo;
1.268 -
1.269 - if ( remcap >= exc ) { //A nonsaturating push.
1.270 -
1.271 - flow->set(e, flo+exc);
1.272 - excess.set(v, excess[v]+exc);
1.273 - exc=0;
1.274 - break;
1.275 -
1.276 - } else { //A saturating push.
1.277 - flow->set(e, cap);
1.278 - excess.set(v, excess[v]+remcap);
1.279 - exc-=remcap;
1.280 - }
1.281 - } else if ( newlevel > level[v] ) newlevel = level[v];
1.282 - } //for out edges wv
1.283 -
1.284 - if ( exc > 0 ) {
1.285 - InEdgeIt e;
1.286 - for(g->first(e,w); g->valid(e); g->next(e)) {
1.287 -
1.288 - if( (*flow)[e] <= 0 ) continue;
1.289 - Node v=g->tail(e);
1.290 -
1.291 - if( lev > level[v] ) { //Push is allowed now
1.292 -
1.293 - if ( excess[v]<=0 && v!=t && v!=s ) {
1.294 - int lev_v=level[v];
1.295 - active[lev_v].push(v);
1.296 - }
1.297 -
1.298 - Num flo=(*flow)[e];
1.299 -
1.300 - if ( flo >= exc ) { //A nonsaturating push.
1.301 -
1.302 - flow->set(e, flo-exc);
1.303 - excess.set(v, excess[v]+exc);
1.304 - exc=0;
1.305 - break;
1.306 - } else { //A saturating push.
1.307 -
1.308 - excess.set(v, excess[v]+flo);
1.309 - exc-=flo;
1.310 - flow->set(e,0);
1.311 - }
1.312 - } else if ( newlevel > level[v] ) newlevel = level[v];
1.313 - } //for in edges vw
1.314 -
1.315 - } // if w still has excess after the out edge for cycle
1.316 -
1.317 - excess.set(w, exc);
1.318 -
1.319 - return newlevel;
1.320 - }
1.321 -
1.322 -
1.323 - void preflowPreproc ( flowEnum fe, VecStack& active,
1.324 - VecNode& level_list, NNMap& left, NNMap& right ) {
1.325 -
1.326 - std::queue<Node> bfs_queue;
1.327 -
1.328 - switch ( fe ) {
1.329 - case ZERO_FLOW:
1.330 - {
1.331 - //Reverse_bfs from t, to find the starting level.
1.332 - level.set(t,0);
1.333 - bfs_queue.push(t);
1.334 -
1.335 - while (!bfs_queue.empty()) {
1.336 -
1.337 - Node v=bfs_queue.front();
1.338 - bfs_queue.pop();
1.339 - int l=level[v]+1;
1.340 -
1.341 - InEdgeIt e;
1.342 - for(g->first(e,v); g->valid(e); g->next(e)) {
1.343 - Node w=g->tail(e);
1.344 - if ( level[w] == n && w != s ) {
1.345 - bfs_queue.push(w);
1.346 - Node first=level_list[l];
1.347 - if ( g->valid(first) ) left.set(first,w);
1.348 - right.set(w,first);
1.349 - level_list[l]=w;
1.350 - level.set(w, l);
1.351 - }
1.352 - }
1.353 - }
1.354 -
1.355 - //the starting flow
1.356 - OutEdgeIt e;
1.357 - for(g->first(e,s); g->valid(e); g->next(e))
1.358 - {
1.359 - Num c=(*capacity)[e];
1.360 - if ( c <= 0 ) continue;
1.361 - Node w=g->head(e);
1.362 - if ( level[w] < n ) {
1.363 - if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
1.364 - flow->set(e, c);
1.365 - excess.set(w, excess[w]+c);
1.366 - }
1.367 - }
1.368 - break;
1.369 - }
1.370 -
1.371 - case GEN_FLOW:
1.372 - case PREFLOW:
1.373 - {
1.374 - //Reverse_bfs from t in the residual graph,
1.375 - //to find the starting level.
1.376 - level.set(t,0);
1.377 - bfs_queue.push(t);
1.378 -
1.379 - while (!bfs_queue.empty()) {
1.380 -
1.381 - Node v=bfs_queue.front();
1.382 - bfs_queue.pop();
1.383 - int l=level[v]+1;
1.384 -
1.385 - InEdgeIt e;
1.386 - for(g->first(e,v); g->valid(e); g->next(e)) {
1.387 - if ( (*capacity)[e] <= (*flow)[e] ) continue;
1.388 - Node w=g->tail(e);
1.389 - if ( level[w] == n && w != s ) {
1.390 - bfs_queue.push(w);
1.391 - Node first=level_list[l];
1.392 - if ( g->valid(first) ) left.set(first,w);
1.393 - right.set(w,first);
1.394 - level_list[l]=w;
1.395 - level.set(w, l);
1.396 - }
1.397 - }
1.398 -
1.399 - OutEdgeIt f;
1.400 - for(g->first(f,v); g->valid(f); g->next(f)) {
1.401 - if ( 0 >= (*flow)[f] ) continue;
1.402 - Node w=g->head(f);
1.403 - if ( level[w] == n && w != s ) {
1.404 - bfs_queue.push(w);
1.405 - Node first=level_list[l];
1.406 - if ( g->valid(first) ) left.set(first,w);
1.407 - right.set(w,first);
1.408 - level_list[l]=w;
1.409 - level.set(w, l);
1.410 - }
1.411 - }
1.412 - }
1.413 -
1.414 -
1.415 - //the starting flow
1.416 - OutEdgeIt e;
1.417 - for(g->first(e,s); g->valid(e); g->next(e))
1.418 - {
1.419 - Num rem=(*capacity)[e]-(*flow)[e];
1.420 - if ( rem <= 0 ) continue;
1.421 - Node w=g->head(e);
1.422 - if ( level[w] < n ) {
1.423 - if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
1.424 - flow->set(e, (*capacity)[e]);
1.425 - excess.set(w, excess[w]+rem);
1.426 - }
1.427 - }
1.428 -
1.429 - InEdgeIt f;
1.430 - for(g->first(f,s); g->valid(f); g->next(f))
1.431 - {
1.432 - if ( (*flow)[f] <= 0 ) continue;
1.433 - Node w=g->tail(f);
1.434 - if ( level[w] < n ) {
1.435 - if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
1.436 - excess.set(w, excess[w]+(*flow)[f]);
1.437 - flow->set(f, 0);
1.438 - }
1.439 - }
1.440 - break;
1.441 - } //case PREFLOW
1.442 - }
1.443 - } //preflowPreproc
1.444 -
1.445 -
1.446 -
1.447 - void relabel(Node w, int newlevel, VecStack& active,
1.448 - VecNode& level_list, NNMap& left,
1.449 - NNMap& right, int& b, int& k, bool what_heur )
1.450 - {
1.451 -
1.452 - Num lev=level[w];
1.453 -
1.454 - Node right_n=right[w];
1.455 - Node left_n=left[w];
1.456 -
1.457 - //unlacing starts
1.458 - if ( g->valid(right_n) ) {
1.459 - if ( g->valid(left_n) ) {
1.460 - right.set(left_n, right_n);
1.461 - left.set(right_n, left_n);
1.462 - } else {
1.463 - level_list[lev]=right_n;
1.464 - left.set(right_n, INVALID);
1.465 - }
1.466 - } else {
1.467 - if ( g->valid(left_n) ) {
1.468 - right.set(left_n, INVALID);
1.469 - } else {
1.470 - level_list[lev]=INVALID;
1.471 - }
1.472 - }
1.473 - //unlacing ends
1.474 -
1.475 - if ( !g->valid(level_list[lev]) ) {
1.476 -
1.477 - //gapping starts
1.478 - for (int i=lev; i!=k ; ) {
1.479 - Node v=level_list[++i];
1.480 - while ( g->valid(v) ) {
1.481 - level.set(v,n);
1.482 - v=right[v];
1.483 - }
1.484 - level_list[i]=INVALID;
1.485 - if ( !what_heur ) {
1.486 - while ( !active[i].empty() ) {
1.487 - active[i].pop(); //FIXME: ezt szebben kene
1.488 - }
1.489 - }
1.490 - }
1.491 -
1.492 - level.set(w,n);
1.493 - b=lev-1;
1.494 - k=b;
1.495 - //gapping ends
1.496 -
1.497 - } else {
1.498 -
1.499 - if ( newlevel == n ) level.set(w,n);
1.500 - else {
1.501 - level.set(w,++newlevel);
1.502 - active[newlevel].push(w);
1.503 - if ( what_heur ) b=newlevel;
1.504 - if ( k < newlevel ) ++k; //now k=newlevel
1.505 - Node first=level_list[newlevel];
1.506 - if ( g->valid(first) ) left.set(first,w);
1.507 - right.set(w,first);
1.508 - left.set(w,INVALID);
1.509 - level_list[newlevel]=w;
1.510 - }
1.511 - }
1.512 -
1.513 - } //relabel
1.514 -
1.515 -
1.516 - template<typename MapGraphWrapper>
1.517 - class DistanceMap {
1.518 - protected:
1.519 - const MapGraphWrapper* g;
1.520 - typename MapGraphWrapper::template NodeMap<int> dist;
1.521 - public:
1.522 - DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
1.523 - void set(const typename MapGraphWrapper::Node& n, int a) {
1.524 - dist.set(n, a);
1.525 - }
1.526 - int operator[](const typename MapGraphWrapper::Node& n)
1.527 - { return dist[n]; }
1.528 -// int get(const typename MapGraphWrapper::Node& n) const {
1.529 -// return dist[n]; }
1.530 -// bool get(const typename MapGraphWrapper::Edge& e) const {
1.531 -// return (dist.get(g->tail(e))<dist.get(g->head(e))); }
1.532 - bool operator[](const typename MapGraphWrapper::Edge& e) const {
1.533 - return (dist[g->tail(e)]<dist[g->head(e)]);
1.534 - }
1.535 - };
1.536 -
1.537 - };
1.538 -
1.539 -
1.540 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.541 - void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe )
1.542 - {
1.543 -
1.544 - int heur0=(int)(H0*n); //time while running 'bound decrease'
1.545 - int heur1=(int)(H1*n); //time while running 'highest label'
1.546 - int heur=heur1; //starting time interval (#of relabels)
1.547 - int numrelabel=0;
1.548 -
1.549 - bool what_heur=1;
1.550 - //It is 0 in case 'bound decrease' and 1 in case 'highest label'
1.551 -
1.552 - bool end=false;
1.553 - //Needed for 'bound decrease', true means no active nodes are above bound b.
1.554 -
1.555 - int k=n-2; //bound on the highest level under n containing a node
1.556 - int b=k; //bound on the highest level under n of an active node
1.557 -
1.558 - VecStack active(n);
1.559 -
1.560 - NNMap left(*g, INVALID);
1.561 - NNMap right(*g, INVALID);
1.562 - VecNode level_list(n,INVALID);
1.563 - //List of the nodes in level i<n, set to n.
1.564 -
1.565 - NodeIt v;
1.566 - for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
1.567 - //setting each node to level n
1.568 -
1.569 - switch ( fe ) {
1.570 - case PREFLOW:
1.571 - {
1.572 - //counting the excess
1.573 - NodeIt v;
1.574 - for(g->first(v); g->valid(v); g->next(v)) {
1.575 - Num exc=0;
1.576 -
1.577 - InEdgeIt e;
1.578 - for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
1.579 - OutEdgeIt f;
1.580 - for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
1.581 -
1.582 - excess.set(v,exc);
1.583 -
1.584 - //putting the active nodes into the stack
1.585 - int lev=level[v];
1.586 - if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
1.587 - }
1.588 - break;
1.589 - }
1.590 - case GEN_FLOW:
1.591 - {
1.592 - //Counting the excess of t
1.593 - Num exc=0;
1.594 -
1.595 - InEdgeIt e;
1.596 - for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
1.597 - OutEdgeIt f;
1.598 - for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
1.599 -
1.600 - excess.set(t,exc);
1.601 -
1.602 - break;
1.603 - }
1.604 - default:
1.605 - break;
1.606 - }
1.607 -
1.608 - preflowPreproc( fe, active, level_list, left, right );
1.609 - //End of preprocessing
1.610 -
1.611 -
1.612 - //Push/relabel on the highest level active nodes.
1.613 - while ( true ) {
1.614 - if ( b == 0 ) {
1.615 - if ( !what_heur && !end && k > 0 ) {
1.616 - b=k;
1.617 - end=true;
1.618 - } else break;
1.619 - }
1.620 -
1.621 - if ( active[b].empty() ) --b;
1.622 - else {
1.623 - end=false;
1.624 - Node w=active[b].top();
1.625 - active[b].pop();
1.626 - int newlevel=push(w,active);
1.627 - if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
1.628 - left, right, b, k, what_heur);
1.629 -
1.630 - ++numrelabel;
1.631 - if ( numrelabel >= heur ) {
1.632 - numrelabel=0;
1.633 - if ( what_heur ) {
1.634 - what_heur=0;
1.635 - heur=heur0;
1.636 - end=false;
1.637 - } else {
1.638 - what_heur=1;
1.639 - heur=heur1;
1.640 - b=k;
1.641 - }
1.642 - }
1.643 - }
1.644 - }
1.645 - }
1.646 -
1.647 -
1.648 -
1.649 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.650 - void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1()
1.651 - {
1.652 -
1.653 - int k=n-2; //bound on the highest level under n containing a node
1.654 - int b=k; //bound on the highest level under n of an active node
1.655 -
1.656 - VecStack active(n);
1.657 - level.set(s,0);
1.658 - std::queue<Node> bfs_queue;
1.659 - bfs_queue.push(s);
1.660 -
1.661 - while (!bfs_queue.empty()) {
1.662 -
1.663 - Node v=bfs_queue.front();
1.664 - bfs_queue.pop();
1.665 - int l=level[v]+1;
1.666 -
1.667 - InEdgeIt e;
1.668 - for(g->first(e,v); g->valid(e); g->next(e)) {
1.669 - if ( (*capacity)[e] <= (*flow)[e] ) continue;
1.670 - Node u=g->tail(e);
1.671 - if ( level[u] >= n ) {
1.672 - bfs_queue.push(u);
1.673 - level.set(u, l);
1.674 - if ( excess[u] > 0 ) active[l].push(u);
1.675 - }
1.676 - }
1.677 -
1.678 - OutEdgeIt f;
1.679 - for(g->first(f,v); g->valid(f); g->next(f)) {
1.680 - if ( 0 >= (*flow)[f] ) continue;
1.681 - Node u=g->head(f);
1.682 - if ( level[u] >= n ) {
1.683 - bfs_queue.push(u);
1.684 - level.set(u, l);
1.685 - if ( excess[u] > 0 ) active[l].push(u);
1.686 - }
1.687 - }
1.688 - }
1.689 - b=n-2;
1.690 -
1.691 - while ( true ) {
1.692 -
1.693 - if ( b == 0 ) break;
1.694 -
1.695 - if ( active[b].empty() ) --b;
1.696 - else {
1.697 - Node w=active[b].top();
1.698 - active[b].pop();
1.699 - int newlevel=push(w,active);
1.700 -
1.701 - //relabel
1.702 - if ( excess[w] > 0 ) {
1.703 - level.set(w,++newlevel);
1.704 - active[newlevel].push(w);
1.705 - b=newlevel;
1.706 - }
1.707 - } // if stack[b] is nonempty
1.708 - } // while(true)
1.709 - }
1.710 -
1.711 -
1.712 -
1.713 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.714 - bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
1.715 - {
1.716 - ResGW res_graph(*g, *capacity, *flow);
1.717 - bool _augment=false;
1.718 -
1.719 - //ReachedMap level(res_graph);
1.720 - FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1.721 - BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
1.722 - bfs.pushAndSetReached(s);
1.723 -
1.724 - typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
1.725 - pred.set(s, INVALID);
1.726 -
1.727 - typename ResGW::template NodeMap<Num> free(res_graph);
1.728 -
1.729 - //searching for augmenting path
1.730 - while ( !bfs.finished() ) {
1.731 - ResGWOutEdgeIt e=bfs;
1.732 - if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
1.733 - Node v=res_graph.tail(e);
1.734 - Node w=res_graph.head(e);
1.735 - pred.set(w, e);
1.736 - if (res_graph.valid(pred[v])) {
1.737 - free.set(w, std::min(free[v], res_graph.resCap(e)));
1.738 - } else {
1.739 - free.set(w, res_graph.resCap(e));
1.740 - }
1.741 - if (res_graph.head(e)==t) { _augment=true; break; }
1.742 - }
1.743 -
1.744 - ++bfs;
1.745 - } //end of searching augmenting path
1.746 -
1.747 - if (_augment) {
1.748 - Node n=t;
1.749 - Num augment_value=free[t];
1.750 - while (res_graph.valid(pred[n])) {
1.751 - ResGWEdge e=pred[n];
1.752 - res_graph.augment(e, augment_value);
1.753 - n=res_graph.tail(e);
1.754 - }
1.755 - }
1.756 -
1.757 - return _augment;
1.758 - }
1.759 -
1.760 -
1.761 -
1.762 -
1.763 -
1.764 -
1.765 -
1.766 -
1.767 -
1.768 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.769 - template<typename MutableGraph>
1.770 - bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
1.771 - {
1.772 - typedef MutableGraph MG;
1.773 - bool _augment=false;
1.774 -
1.775 - ResGW res_graph(*g, *capacity, *flow);
1.776 -
1.777 - //bfs for distances on the residual graph
1.778 - //ReachedMap level(res_graph);
1.779 - FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1.780 - BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
1.781 - bfs.pushAndSetReached(s);
1.782 - typename ResGW::template NodeMap<int>
1.783 - dist(res_graph); //filled up with 0's
1.784 -
1.785 - //F will contain the physical copy of the residual graph
1.786 - //with the set of edges which are on shortest paths
1.787 - MG F;
1.788 - typename ResGW::template NodeMap<typename MG::Node>
1.789 - res_graph_to_F(res_graph);
1.790 - {
1.791 - typename ResGW::NodeIt n;
1.792 - for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
1.793 - res_graph_to_F.set(n, F.addNode());
1.794 - }
1.795 - }
1.796 -
1.797 - typename MG::Node sF=res_graph_to_F[s];
1.798 - typename MG::Node tF=res_graph_to_F[t];
1.799 - typename MG::template EdgeMap<ResGWEdge> original_edge(F);
1.800 - typename MG::template EdgeMap<Num> residual_capacity(F);
1.801 -
1.802 - while ( !bfs.finished() ) {
1.803 - ResGWOutEdgeIt e=bfs;
1.804 - if (res_graph.valid(e)) {
1.805 - if (bfs.isBNodeNewlyReached()) {
1.806 - dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
1.807 - typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
1.808 - original_edge.update();
1.809 - original_edge.set(f, e);
1.810 - residual_capacity.update();
1.811 - residual_capacity.set(f, res_graph.resCap(e));
1.812 - } else {
1.813 - if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
1.814 - typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
1.815 - original_edge.update();
1.816 - original_edge.set(f, e);
1.817 - residual_capacity.update();
1.818 - residual_capacity.set(f, res_graph.resCap(e));
1.819 - }
1.820 - }
1.821 - }
1.822 - ++bfs;
1.823 - } //computing distances from s in the residual graph
1.824 -
1.825 - bool __augment=true;
1.826 -
1.827 - while (__augment) {
1.828 - __augment=false;
1.829 - //computing blocking flow with dfs
1.830 - DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
1.831 - typename MG::template NodeMap<typename MG::Edge> pred(F);
1.832 - pred.set(sF, INVALID);
1.833 - //invalid iterators for sources
1.834 -
1.835 - typename MG::template NodeMap<Num> free(F);
1.836 -
1.837 - dfs.pushAndSetReached(sF);
1.838 - while (!dfs.finished()) {
1.839 - ++dfs;
1.840 - if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
1.841 - if (dfs.isBNodeNewlyReached()) {
1.842 - typename MG::Node v=F.aNode(dfs);
1.843 - typename MG::Node w=F.bNode(dfs);
1.844 - pred.set(w, dfs);
1.845 - if (F.valid(pred[v])) {
1.846 - free.set(w, std::min(free[v], residual_capacity[dfs]));
1.847 - } else {
1.848 - free.set(w, residual_capacity[dfs]);
1.849 - }
1.850 - if (w==tF) {
1.851 - __augment=true;
1.852 - _augment=true;
1.853 - break;
1.854 - }
1.855 -
1.856 - } else {
1.857 - F.erase(/*typename MG::OutEdgeIt*/(dfs));
1.858 - }
1.859 - }
1.860 - }
1.861 -
1.862 - if (__augment) {
1.863 - typename MG::Node n=tF;
1.864 - Num augment_value=free[tF];
1.865 - while (F.valid(pred[n])) {
1.866 - typename MG::Edge e=pred[n];
1.867 - res_graph.augment(original_edge[e], augment_value);
1.868 - n=F.tail(e);
1.869 - if (residual_capacity[e]==augment_value)
1.870 - F.erase(e);
1.871 - else
1.872 - residual_capacity.set(e, residual_capacity[e]-augment_value);
1.873 - }
1.874 - }
1.875 -
1.876 - }
1.877 -
1.878 - return _augment;
1.879 - }
1.880 -
1.881 -
1.882 -
1.883 -
1.884 -
1.885 -
1.886 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
1.887 - bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
1.888 - {
1.889 - bool _augment=false;
1.890 -
1.891 - ResGW res_graph(*g, *capacity, *flow);
1.892 -
1.893 - //ReachedMap level(res_graph);
1.894 - FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
1.895 - BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
1.896 -
1.897 - bfs.pushAndSetReached(s);
1.898 - DistanceMap<ResGW> dist(res_graph);
1.899 - while ( !bfs.finished() ) {
1.900 - ResGWOutEdgeIt e=bfs;
1.901 - if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
1.902 - dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
1.903 - }
1.904 - ++bfs;
1.905 - } //computing distances from s in the residual graph
1.906 -
1.907 - //Subgraph containing the edges on some shortest paths
1.908 - ConstMap<typename ResGW::Node, bool> true_map(true);
1.909 - typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
1.910 - DistanceMap<ResGW> > FilterResGW;
1.911 - FilterResGW filter_res_graph(res_graph, true_map, dist);
1.912 -
1.913 - //Subgraph, which is able to delete edges which are already
1.914 - //met by the dfs
1.915 - typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
1.916 - first_out_edges(filter_res_graph);
1.917 - typename FilterResGW::NodeIt v;
1.918 - for(filter_res_graph.first(v); filter_res_graph.valid(v);
1.919 - filter_res_graph.next(v))
1.920 - {
1.921 - typename FilterResGW::OutEdgeIt e;
1.922 - filter_res_graph.first(e, v);
1.923 - first_out_edges.set(v, e);
1.924 - }
1.925 - typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
1.926 - template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
1.927 - ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
1.928 -
1.929 - bool __augment=true;
1.930 -
1.931 - while (__augment) {
1.932 -
1.933 - __augment=false;
1.934 - //computing blocking flow with dfs
1.935 - DfsIterator< ErasingResGW,
1.936 - typename ErasingResGW::template NodeMap<bool> >
1.937 - dfs(erasing_res_graph);
1.938 - typename ErasingResGW::
1.939 - template NodeMap<typename ErasingResGW::OutEdgeIt>
1.940 - pred(erasing_res_graph);
1.941 - pred.set(s, INVALID);
1.942 - //invalid iterators for sources
1.943 -
1.944 - typename ErasingResGW::template NodeMap<Num>
1.945 - free1(erasing_res_graph);
1.946 -
1.947 - dfs.pushAndSetReached(
1.948 - typename ErasingResGW::Node(
1.949 - typename FilterResGW::Node(
1.950 - typename ResGW::Node(s)
1.951 - )
1.952 - )
1.953 - );
1.954 - while (!dfs.finished()) {
1.955 - ++dfs;
1.956 - if (erasing_res_graph.valid(
1.957 - typename ErasingResGW::OutEdgeIt(dfs)))
1.958 - {
1.959 - if (dfs.isBNodeNewlyReached()) {
1.960 -
1.961 - typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
1.962 - typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
1.963 -
1.964 - pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
1.965 - if (erasing_res_graph.valid(pred[v])) {
1.966 - free1.set(w, std::min(free1[v], res_graph.resCap(
1.967 - typename ErasingResGW::OutEdgeIt(dfs))));
1.968 - } else {
1.969 - free1.set(w, res_graph.resCap(
1.970 - typename ErasingResGW::OutEdgeIt(dfs)));
1.971 - }
1.972 -
1.973 - if (w==t) {
1.974 - __augment=true;
1.975 - _augment=true;
1.976 - break;
1.977 - }
1.978 - } else {
1.979 - erasing_res_graph.erase(dfs);
1.980 - }
1.981 - }
1.982 - }
1.983 -
1.984 - if (__augment) {
1.985 - typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
1.986 -// typename ResGW::NodeMap<Num> a(res_graph);
1.987 -// typename ResGW::Node b;
1.988 -// Num j=a[b];
1.989 -// typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
1.990 -// typename FilterResGW::Node b1;
1.991 -// Num j1=a1[b1];
1.992 -// typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
1.993 -// typename ErasingResGW::Node b2;
1.994 -// Num j2=a2[b2];
1.995 - Num augment_value=free1[n];
1.996 - while (erasing_res_graph.valid(pred[n])) {
1.997 - typename ErasingResGW::OutEdgeIt e=pred[n];
1.998 - res_graph.augment(e, augment_value);
1.999 - n=erasing_res_graph.tail(e);
1.1000 - if (res_graph.resCap(e)==0)
1.1001 - erasing_res_graph.erase(e);
1.1002 - }
1.1003 - }
1.1004 -
1.1005 - } //while (__augment)
1.1006 -
1.1007 - return _augment;
1.1008 - }
1.1009 -
1.1010 -
1.1011 -
1.1012 -
1.1013 -} //namespace hugo
1.1014 -
1.1015 -#endif //HUGO_PREFLOW_H
1.1016 -
1.1017 -
1.1018 -
1.1019 -