1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/jacint/max_flow.h Thu Apr 29 16:26:01 2004 +0000
1.3 @@ -0,0 +1,1016 @@
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 +
2.1 --- a/src/work/jacint/preflow.h Thu Apr 29 16:26:01 2004 +0000
2.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
2.3 @@ -1,1016 +0,0 @@
2.4 -// -*- C++ -*-
2.5 -
2.6 -/*
2.7 -Heuristics:
2.8 - 2 phase
2.9 - gap
2.10 - list 'level_list' on the nodes on level i implemented by hand
2.11 - stack 'active' on the active nodes on level i
2.12 - runs heuristic 'highest label' for H1*n relabels
2.13 - runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
2.14 -
2.15 -Parameters H0 and H1 are initialized to 20 and 1.
2.16 -
2.17 -Constructors:
2.18 -
2.19 -Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if
2.20 - FlowMap is not constant zero, and should be true if it is
2.21 -
2.22 -Members:
2.23 -
2.24 -void run()
2.25 -
2.26 -Num flowValue() : returns the value of a maximum flow
2.27 -
2.28 -void minMinCut(CutMap& M) : sets M to the characteristic vector of the
2.29 - minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
2.30 -
2.31 -void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
2.32 - maximum min cut. M should be a map of bools initialized to false.
2.33 -
2.34 -void minCut(CutMap& M) : sets M to the characteristic vector of
2.35 - a min cut. M should be a map of bools initialized to false.
2.36 -
2.37 -*/
2.38 -
2.39 -#ifndef HUGO_PREFLOW_H
2.40 -#define HUGO_PREFLOW_H
2.41 -
2.42 -#define H0 20
2.43 -#define H1 1
2.44 -
2.45 -#include <vector>
2.46 -#include <queue>
2.47 -#include <stack>
2.48 -
2.49 -#include <graph_wrapper.h>
2.50 -#include <bfs_iterator.h>
2.51 -#include <invalid.h>
2.52 -#include <maps.h>
2.53 -#include <for_each_macros.h>
2.54 -
2.55 -
2.56 -namespace hugo {
2.57 -
2.58 - template <typename Graph, typename Num,
2.59 - typename CapMap=typename Graph::template EdgeMap<Num>,
2.60 - typename FlowMap=typename Graph::template EdgeMap<Num> >
2.61 - class MaxFlow {
2.62 -
2.63 - typedef typename Graph::Node Node;
2.64 - typedef typename Graph::NodeIt NodeIt;
2.65 - typedef typename Graph::OutEdgeIt OutEdgeIt;
2.66 - typedef typename Graph::InEdgeIt InEdgeIt;
2.67 -
2.68 - typedef typename std::vector<std::stack<Node> > VecStack;
2.69 - typedef typename Graph::template NodeMap<Node> NNMap;
2.70 - typedef typename std::vector<Node> VecNode;
2.71 -
2.72 - const Graph* g;
2.73 - Node s;
2.74 - Node t;
2.75 - const CapMap* capacity;
2.76 - FlowMap* flow;
2.77 - int n; //the number of nodes of G
2.78 - typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
2.79 - typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
2.80 - typedef typename ResGW::Edge ResGWEdge;
2.81 - //typedef typename ResGW::template NodeMap<bool> ReachedMap;
2.82 - typedef typename Graph::template NodeMap<int> ReachedMap;
2.83 - ReachedMap level;
2.84 - //level works as a bool map in augmenting path algorithms
2.85 - //and is used by bfs for storing reached information.
2.86 - //In preflow, it shows levels of nodes.
2.87 - //typename Graph::template NodeMap<int> level;
2.88 - typename Graph::template NodeMap<Num> excess;
2.89 -
2.90 - public:
2.91 -
2.92 - enum flowEnum{
2.93 - ZERO_FLOW=0,
2.94 - GEN_FLOW=1,
2.95 - PREFLOW=2
2.96 - };
2.97 -
2.98 - MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
2.99 - FlowMap& _flow) :
2.100 - g(&_G), s(_s), t(_t), capacity(&_capacity),
2.101 - flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
2.102 -
2.103 - void run() {
2.104 - preflow( ZERO_FLOW );
2.105 - }
2.106 -
2.107 - void preflow( flowEnum fe ) {
2.108 - preflowPhase0(fe);
2.109 - preflowPhase1();
2.110 - }
2.111 -
2.112 - void preflowPhase0( flowEnum fe );
2.113 -
2.114 - void preflowPhase1();
2.115 -
2.116 - bool augmentOnShortestPath();
2.117 -
2.118 - template<typename MutableGraph> bool augmentOnBlockingFlow();
2.119 -
2.120 - bool augmentOnBlockingFlow2();
2.121 -
2.122 - /// Returns the actual flow value.
2.123 - /// More precisely, it returns the negative excess of s, thus
2.124 - /// this works also for preflows.
2.125 - Num flowValue() {
2.126 - Num a=0;
2.127 - FOR_EACH_INC_LOC(OutEdgeIt, e, *g, s) a+=(*flow)[e];
2.128 - FOR_EACH_INC_LOC(InEdgeIt, e, *g, s) a-=(*flow)[e];
2.129 - return a;
2.130 - }
2.131 -
2.132 - //should be used only between preflowPhase0 and preflowPhase1
2.133 - template<typename _CutMap>
2.134 - void actMinCut(_CutMap& M) {
2.135 - NodeIt v;
2.136 - for(g->first(v); g->valid(v); g->next(v))
2.137 - if ( level[v] < n ) {
2.138 - M.set(v,false);
2.139 - } else {
2.140 - M.set(v,true);
2.141 - }
2.142 - }
2.143 -
2.144 -
2.145 -
2.146 - /*
2.147 - Returns the minimum min cut, by a bfs from s in the residual graph.
2.148 - */
2.149 - template<typename _CutMap>
2.150 - void minMinCut(_CutMap& M) {
2.151 -
2.152 - std::queue<Node> queue;
2.153 -
2.154 - M.set(s,true);
2.155 - queue.push(s);
2.156 -
2.157 - while (!queue.empty()) {
2.158 - Node w=queue.front();
2.159 - queue.pop();
2.160 -
2.161 - OutEdgeIt e;
2.162 - for(g->first(e,w) ; g->valid(e); g->next(e)) {
2.163 - Node v=g->head(e);
2.164 - if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
2.165 - queue.push(v);
2.166 - M.set(v, true);
2.167 - }
2.168 - }
2.169 -
2.170 - InEdgeIt f;
2.171 - for(g->first(f,w) ; g->valid(f); g->next(f)) {
2.172 - Node v=g->tail(f);
2.173 - if (!M[v] && (*flow)[f] > 0 ) {
2.174 - queue.push(v);
2.175 - M.set(v, true);
2.176 - }
2.177 - }
2.178 - }
2.179 - }
2.180 -
2.181 -
2.182 -
2.183 - /*
2.184 - Returns the maximum min cut, by a reverse bfs
2.185 - from t in the residual graph.
2.186 - */
2.187 -
2.188 - template<typename _CutMap>
2.189 - void maxMinCut(_CutMap& M) {
2.190 -
2.191 - NodeIt v;
2.192 - for(g->first(v) ; g->valid(v); g->next(v)) {
2.193 - M.set(v, true);
2.194 - }
2.195 -
2.196 - std::queue<Node> queue;
2.197 -
2.198 - M.set(t,false);
2.199 - queue.push(t);
2.200 -
2.201 - while (!queue.empty()) {
2.202 - Node w=queue.front();
2.203 - queue.pop();
2.204 -
2.205 -
2.206 - InEdgeIt e;
2.207 - for(g->first(e,w) ; g->valid(e); g->next(e)) {
2.208 - Node v=g->tail(e);
2.209 - if (M[v] && (*flow)[e] < (*capacity)[e] ) {
2.210 - queue.push(v);
2.211 - M.set(v, false);
2.212 - }
2.213 - }
2.214 -
2.215 - OutEdgeIt f;
2.216 - for(g->first(f,w) ; g->valid(f); g->next(f)) {
2.217 - Node v=g->head(f);
2.218 - if (M[v] && (*flow)[f] > 0 ) {
2.219 - queue.push(v);
2.220 - M.set(v, false);
2.221 - }
2.222 - }
2.223 - }
2.224 - }
2.225 -
2.226 -
2.227 - template<typename CutMap>
2.228 - void minCut(CutMap& M) {
2.229 - minMinCut(M);
2.230 - }
2.231 -
2.232 - void resetTarget(Node _t) {t=_t;}
2.233 - void resetSource(Node _s) {s=_s;}
2.234 -
2.235 - void resetCap(const CapMap& _cap) {
2.236 - capacity=&_cap;
2.237 - }
2.238 -
2.239 - void resetFlow(FlowMap& _flow) {
2.240 - flow=&_flow;
2.241 - }
2.242 -
2.243 -
2.244 - private:
2.245 -
2.246 - int push(Node w, VecStack& active) {
2.247 -
2.248 - int lev=level[w];
2.249 - Num exc=excess[w];
2.250 - int newlevel=n; //bound on the next level of w
2.251 -
2.252 - OutEdgeIt e;
2.253 - for(g->first(e,w); g->valid(e); g->next(e)) {
2.254 -
2.255 - if ( (*flow)[e] >= (*capacity)[e] ) continue;
2.256 - Node v=g->head(e);
2.257 -
2.258 - if( lev > level[v] ) { //Push is allowed now
2.259 -
2.260 - if ( excess[v]<=0 && v!=t && v!=s ) {
2.261 - int lev_v=level[v];
2.262 - active[lev_v].push(v);
2.263 - }
2.264 -
2.265 - Num cap=(*capacity)[e];
2.266 - Num flo=(*flow)[e];
2.267 - Num remcap=cap-flo;
2.268 -
2.269 - if ( remcap >= exc ) { //A nonsaturating push.
2.270 -
2.271 - flow->set(e, flo+exc);
2.272 - excess.set(v, excess[v]+exc);
2.273 - exc=0;
2.274 - break;
2.275 -
2.276 - } else { //A saturating push.
2.277 - flow->set(e, cap);
2.278 - excess.set(v, excess[v]+remcap);
2.279 - exc-=remcap;
2.280 - }
2.281 - } else if ( newlevel > level[v] ) newlevel = level[v];
2.282 - } //for out edges wv
2.283 -
2.284 - if ( exc > 0 ) {
2.285 - InEdgeIt e;
2.286 - for(g->first(e,w); g->valid(e); g->next(e)) {
2.287 -
2.288 - if( (*flow)[e] <= 0 ) continue;
2.289 - Node v=g->tail(e);
2.290 -
2.291 - if( lev > level[v] ) { //Push is allowed now
2.292 -
2.293 - if ( excess[v]<=0 && v!=t && v!=s ) {
2.294 - int lev_v=level[v];
2.295 - active[lev_v].push(v);
2.296 - }
2.297 -
2.298 - Num flo=(*flow)[e];
2.299 -
2.300 - if ( flo >= exc ) { //A nonsaturating push.
2.301 -
2.302 - flow->set(e, flo-exc);
2.303 - excess.set(v, excess[v]+exc);
2.304 - exc=0;
2.305 - break;
2.306 - } else { //A saturating push.
2.307 -
2.308 - excess.set(v, excess[v]+flo);
2.309 - exc-=flo;
2.310 - flow->set(e,0);
2.311 - }
2.312 - } else if ( newlevel > level[v] ) newlevel = level[v];
2.313 - } //for in edges vw
2.314 -
2.315 - } // if w still has excess after the out edge for cycle
2.316 -
2.317 - excess.set(w, exc);
2.318 -
2.319 - return newlevel;
2.320 - }
2.321 -
2.322 -
2.323 - void preflowPreproc ( flowEnum fe, VecStack& active,
2.324 - VecNode& level_list, NNMap& left, NNMap& right ) {
2.325 -
2.326 - std::queue<Node> bfs_queue;
2.327 -
2.328 - switch ( fe ) {
2.329 - case ZERO_FLOW:
2.330 - {
2.331 - //Reverse_bfs from t, to find the starting level.
2.332 - level.set(t,0);
2.333 - bfs_queue.push(t);
2.334 -
2.335 - while (!bfs_queue.empty()) {
2.336 -
2.337 - Node v=bfs_queue.front();
2.338 - bfs_queue.pop();
2.339 - int l=level[v]+1;
2.340 -
2.341 - InEdgeIt e;
2.342 - for(g->first(e,v); g->valid(e); g->next(e)) {
2.343 - Node w=g->tail(e);
2.344 - if ( level[w] == n && w != s ) {
2.345 - bfs_queue.push(w);
2.346 - Node first=level_list[l];
2.347 - if ( g->valid(first) ) left.set(first,w);
2.348 - right.set(w,first);
2.349 - level_list[l]=w;
2.350 - level.set(w, l);
2.351 - }
2.352 - }
2.353 - }
2.354 -
2.355 - //the starting flow
2.356 - OutEdgeIt e;
2.357 - for(g->first(e,s); g->valid(e); g->next(e))
2.358 - {
2.359 - Num c=(*capacity)[e];
2.360 - if ( c <= 0 ) continue;
2.361 - Node w=g->head(e);
2.362 - if ( level[w] < n ) {
2.363 - if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
2.364 - flow->set(e, c);
2.365 - excess.set(w, excess[w]+c);
2.366 - }
2.367 - }
2.368 - break;
2.369 - }
2.370 -
2.371 - case GEN_FLOW:
2.372 - case PREFLOW:
2.373 - {
2.374 - //Reverse_bfs from t in the residual graph,
2.375 - //to find the starting level.
2.376 - level.set(t,0);
2.377 - bfs_queue.push(t);
2.378 -
2.379 - while (!bfs_queue.empty()) {
2.380 -
2.381 - Node v=bfs_queue.front();
2.382 - bfs_queue.pop();
2.383 - int l=level[v]+1;
2.384 -
2.385 - InEdgeIt e;
2.386 - for(g->first(e,v); g->valid(e); g->next(e)) {
2.387 - if ( (*capacity)[e] <= (*flow)[e] ) continue;
2.388 - Node w=g->tail(e);
2.389 - if ( level[w] == n && w != s ) {
2.390 - bfs_queue.push(w);
2.391 - Node first=level_list[l];
2.392 - if ( g->valid(first) ) left.set(first,w);
2.393 - right.set(w,first);
2.394 - level_list[l]=w;
2.395 - level.set(w, l);
2.396 - }
2.397 - }
2.398 -
2.399 - OutEdgeIt f;
2.400 - for(g->first(f,v); g->valid(f); g->next(f)) {
2.401 - if ( 0 >= (*flow)[f] ) continue;
2.402 - Node w=g->head(f);
2.403 - if ( level[w] == n && w != s ) {
2.404 - bfs_queue.push(w);
2.405 - Node first=level_list[l];
2.406 - if ( g->valid(first) ) left.set(first,w);
2.407 - right.set(w,first);
2.408 - level_list[l]=w;
2.409 - level.set(w, l);
2.410 - }
2.411 - }
2.412 - }
2.413 -
2.414 -
2.415 - //the starting flow
2.416 - OutEdgeIt e;
2.417 - for(g->first(e,s); g->valid(e); g->next(e))
2.418 - {
2.419 - Num rem=(*capacity)[e]-(*flow)[e];
2.420 - if ( rem <= 0 ) continue;
2.421 - Node w=g->head(e);
2.422 - if ( level[w] < n ) {
2.423 - if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
2.424 - flow->set(e, (*capacity)[e]);
2.425 - excess.set(w, excess[w]+rem);
2.426 - }
2.427 - }
2.428 -
2.429 - InEdgeIt f;
2.430 - for(g->first(f,s); g->valid(f); g->next(f))
2.431 - {
2.432 - if ( (*flow)[f] <= 0 ) continue;
2.433 - Node w=g->tail(f);
2.434 - if ( level[w] < n ) {
2.435 - if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w);
2.436 - excess.set(w, excess[w]+(*flow)[f]);
2.437 - flow->set(f, 0);
2.438 - }
2.439 - }
2.440 - break;
2.441 - } //case PREFLOW
2.442 - }
2.443 - } //preflowPreproc
2.444 -
2.445 -
2.446 -
2.447 - void relabel(Node w, int newlevel, VecStack& active,
2.448 - VecNode& level_list, NNMap& left,
2.449 - NNMap& right, int& b, int& k, bool what_heur )
2.450 - {
2.451 -
2.452 - Num lev=level[w];
2.453 -
2.454 - Node right_n=right[w];
2.455 - Node left_n=left[w];
2.456 -
2.457 - //unlacing starts
2.458 - if ( g->valid(right_n) ) {
2.459 - if ( g->valid(left_n) ) {
2.460 - right.set(left_n, right_n);
2.461 - left.set(right_n, left_n);
2.462 - } else {
2.463 - level_list[lev]=right_n;
2.464 - left.set(right_n, INVALID);
2.465 - }
2.466 - } else {
2.467 - if ( g->valid(left_n) ) {
2.468 - right.set(left_n, INVALID);
2.469 - } else {
2.470 - level_list[lev]=INVALID;
2.471 - }
2.472 - }
2.473 - //unlacing ends
2.474 -
2.475 - if ( !g->valid(level_list[lev]) ) {
2.476 -
2.477 - //gapping starts
2.478 - for (int i=lev; i!=k ; ) {
2.479 - Node v=level_list[++i];
2.480 - while ( g->valid(v) ) {
2.481 - level.set(v,n);
2.482 - v=right[v];
2.483 - }
2.484 - level_list[i]=INVALID;
2.485 - if ( !what_heur ) {
2.486 - while ( !active[i].empty() ) {
2.487 - active[i].pop(); //FIXME: ezt szebben kene
2.488 - }
2.489 - }
2.490 - }
2.491 -
2.492 - level.set(w,n);
2.493 - b=lev-1;
2.494 - k=b;
2.495 - //gapping ends
2.496 -
2.497 - } else {
2.498 -
2.499 - if ( newlevel == n ) level.set(w,n);
2.500 - else {
2.501 - level.set(w,++newlevel);
2.502 - active[newlevel].push(w);
2.503 - if ( what_heur ) b=newlevel;
2.504 - if ( k < newlevel ) ++k; //now k=newlevel
2.505 - Node first=level_list[newlevel];
2.506 - if ( g->valid(first) ) left.set(first,w);
2.507 - right.set(w,first);
2.508 - left.set(w,INVALID);
2.509 - level_list[newlevel]=w;
2.510 - }
2.511 - }
2.512 -
2.513 - } //relabel
2.514 -
2.515 -
2.516 - template<typename MapGraphWrapper>
2.517 - class DistanceMap {
2.518 - protected:
2.519 - const MapGraphWrapper* g;
2.520 - typename MapGraphWrapper::template NodeMap<int> dist;
2.521 - public:
2.522 - DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
2.523 - void set(const typename MapGraphWrapper::Node& n, int a) {
2.524 - dist.set(n, a);
2.525 - }
2.526 - int operator[](const typename MapGraphWrapper::Node& n)
2.527 - { return dist[n]; }
2.528 -// int get(const typename MapGraphWrapper::Node& n) const {
2.529 -// return dist[n]; }
2.530 -// bool get(const typename MapGraphWrapper::Edge& e) const {
2.531 -// return (dist.get(g->tail(e))<dist.get(g->head(e))); }
2.532 - bool operator[](const typename MapGraphWrapper::Edge& e) const {
2.533 - return (dist[g->tail(e)]<dist[g->head(e)]);
2.534 - }
2.535 - };
2.536 -
2.537 - };
2.538 -
2.539 -
2.540 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.541 - void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe )
2.542 - {
2.543 -
2.544 - int heur0=(int)(H0*n); //time while running 'bound decrease'
2.545 - int heur1=(int)(H1*n); //time while running 'highest label'
2.546 - int heur=heur1; //starting time interval (#of relabels)
2.547 - int numrelabel=0;
2.548 -
2.549 - bool what_heur=1;
2.550 - //It is 0 in case 'bound decrease' and 1 in case 'highest label'
2.551 -
2.552 - bool end=false;
2.553 - //Needed for 'bound decrease', true means no active nodes are above bound b.
2.554 -
2.555 - int k=n-2; //bound on the highest level under n containing a node
2.556 - int b=k; //bound on the highest level under n of an active node
2.557 -
2.558 - VecStack active(n);
2.559 -
2.560 - NNMap left(*g, INVALID);
2.561 - NNMap right(*g, INVALID);
2.562 - VecNode level_list(n,INVALID);
2.563 - //List of the nodes in level i<n, set to n.
2.564 -
2.565 - NodeIt v;
2.566 - for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
2.567 - //setting each node to level n
2.568 -
2.569 - switch ( fe ) {
2.570 - case PREFLOW:
2.571 - {
2.572 - //counting the excess
2.573 - NodeIt v;
2.574 - for(g->first(v); g->valid(v); g->next(v)) {
2.575 - Num exc=0;
2.576 -
2.577 - InEdgeIt e;
2.578 - for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
2.579 - OutEdgeIt f;
2.580 - for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
2.581 -
2.582 - excess.set(v,exc);
2.583 -
2.584 - //putting the active nodes into the stack
2.585 - int lev=level[v];
2.586 - if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
2.587 - }
2.588 - break;
2.589 - }
2.590 - case GEN_FLOW:
2.591 - {
2.592 - //Counting the excess of t
2.593 - Num exc=0;
2.594 -
2.595 - InEdgeIt e;
2.596 - for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
2.597 - OutEdgeIt f;
2.598 - for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
2.599 -
2.600 - excess.set(t,exc);
2.601 -
2.602 - break;
2.603 - }
2.604 - default:
2.605 - break;
2.606 - }
2.607 -
2.608 - preflowPreproc( fe, active, level_list, left, right );
2.609 - //End of preprocessing
2.610 -
2.611 -
2.612 - //Push/relabel on the highest level active nodes.
2.613 - while ( true ) {
2.614 - if ( b == 0 ) {
2.615 - if ( !what_heur && !end && k > 0 ) {
2.616 - b=k;
2.617 - end=true;
2.618 - } else break;
2.619 - }
2.620 -
2.621 - if ( active[b].empty() ) --b;
2.622 - else {
2.623 - end=false;
2.624 - Node w=active[b].top();
2.625 - active[b].pop();
2.626 - int newlevel=push(w,active);
2.627 - if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list,
2.628 - left, right, b, k, what_heur);
2.629 -
2.630 - ++numrelabel;
2.631 - if ( numrelabel >= heur ) {
2.632 - numrelabel=0;
2.633 - if ( what_heur ) {
2.634 - what_heur=0;
2.635 - heur=heur0;
2.636 - end=false;
2.637 - } else {
2.638 - what_heur=1;
2.639 - heur=heur1;
2.640 - b=k;
2.641 - }
2.642 - }
2.643 - }
2.644 - }
2.645 - }
2.646 -
2.647 -
2.648 -
2.649 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.650 - void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1()
2.651 - {
2.652 -
2.653 - int k=n-2; //bound on the highest level under n containing a node
2.654 - int b=k; //bound on the highest level under n of an active node
2.655 -
2.656 - VecStack active(n);
2.657 - level.set(s,0);
2.658 - std::queue<Node> bfs_queue;
2.659 - bfs_queue.push(s);
2.660 -
2.661 - while (!bfs_queue.empty()) {
2.662 -
2.663 - Node v=bfs_queue.front();
2.664 - bfs_queue.pop();
2.665 - int l=level[v]+1;
2.666 -
2.667 - InEdgeIt e;
2.668 - for(g->first(e,v); g->valid(e); g->next(e)) {
2.669 - if ( (*capacity)[e] <= (*flow)[e] ) continue;
2.670 - Node u=g->tail(e);
2.671 - if ( level[u] >= n ) {
2.672 - bfs_queue.push(u);
2.673 - level.set(u, l);
2.674 - if ( excess[u] > 0 ) active[l].push(u);
2.675 - }
2.676 - }
2.677 -
2.678 - OutEdgeIt f;
2.679 - for(g->first(f,v); g->valid(f); g->next(f)) {
2.680 - if ( 0 >= (*flow)[f] ) continue;
2.681 - Node u=g->head(f);
2.682 - if ( level[u] >= n ) {
2.683 - bfs_queue.push(u);
2.684 - level.set(u, l);
2.685 - if ( excess[u] > 0 ) active[l].push(u);
2.686 - }
2.687 - }
2.688 - }
2.689 - b=n-2;
2.690 -
2.691 - while ( true ) {
2.692 -
2.693 - if ( b == 0 ) break;
2.694 -
2.695 - if ( active[b].empty() ) --b;
2.696 - else {
2.697 - Node w=active[b].top();
2.698 - active[b].pop();
2.699 - int newlevel=push(w,active);
2.700 -
2.701 - //relabel
2.702 - if ( excess[w] > 0 ) {
2.703 - level.set(w,++newlevel);
2.704 - active[newlevel].push(w);
2.705 - b=newlevel;
2.706 - }
2.707 - } // if stack[b] is nonempty
2.708 - } // while(true)
2.709 - }
2.710 -
2.711 -
2.712 -
2.713 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.714 - bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
2.715 - {
2.716 - ResGW res_graph(*g, *capacity, *flow);
2.717 - bool _augment=false;
2.718 -
2.719 - //ReachedMap level(res_graph);
2.720 - FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
2.721 - BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
2.722 - bfs.pushAndSetReached(s);
2.723 -
2.724 - typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
2.725 - pred.set(s, INVALID);
2.726 -
2.727 - typename ResGW::template NodeMap<Num> free(res_graph);
2.728 -
2.729 - //searching for augmenting path
2.730 - while ( !bfs.finished() ) {
2.731 - ResGWOutEdgeIt e=bfs;
2.732 - if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
2.733 - Node v=res_graph.tail(e);
2.734 - Node w=res_graph.head(e);
2.735 - pred.set(w, e);
2.736 - if (res_graph.valid(pred[v])) {
2.737 - free.set(w, std::min(free[v], res_graph.resCap(e)));
2.738 - } else {
2.739 - free.set(w, res_graph.resCap(e));
2.740 - }
2.741 - if (res_graph.head(e)==t) { _augment=true; break; }
2.742 - }
2.743 -
2.744 - ++bfs;
2.745 - } //end of searching augmenting path
2.746 -
2.747 - if (_augment) {
2.748 - Node n=t;
2.749 - Num augment_value=free[t];
2.750 - while (res_graph.valid(pred[n])) {
2.751 - ResGWEdge e=pred[n];
2.752 - res_graph.augment(e, augment_value);
2.753 - n=res_graph.tail(e);
2.754 - }
2.755 - }
2.756 -
2.757 - return _augment;
2.758 - }
2.759 -
2.760 -
2.761 -
2.762 -
2.763 -
2.764 -
2.765 -
2.766 -
2.767 -
2.768 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.769 - template<typename MutableGraph>
2.770 - bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
2.771 - {
2.772 - typedef MutableGraph MG;
2.773 - bool _augment=false;
2.774 -
2.775 - ResGW res_graph(*g, *capacity, *flow);
2.776 -
2.777 - //bfs for distances on the residual graph
2.778 - //ReachedMap level(res_graph);
2.779 - FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
2.780 - BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
2.781 - bfs.pushAndSetReached(s);
2.782 - typename ResGW::template NodeMap<int>
2.783 - dist(res_graph); //filled up with 0's
2.784 -
2.785 - //F will contain the physical copy of the residual graph
2.786 - //with the set of edges which are on shortest paths
2.787 - MG F;
2.788 - typename ResGW::template NodeMap<typename MG::Node>
2.789 - res_graph_to_F(res_graph);
2.790 - {
2.791 - typename ResGW::NodeIt n;
2.792 - for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
2.793 - res_graph_to_F.set(n, F.addNode());
2.794 - }
2.795 - }
2.796 -
2.797 - typename MG::Node sF=res_graph_to_F[s];
2.798 - typename MG::Node tF=res_graph_to_F[t];
2.799 - typename MG::template EdgeMap<ResGWEdge> original_edge(F);
2.800 - typename MG::template EdgeMap<Num> residual_capacity(F);
2.801 -
2.802 - while ( !bfs.finished() ) {
2.803 - ResGWOutEdgeIt e=bfs;
2.804 - if (res_graph.valid(e)) {
2.805 - if (bfs.isBNodeNewlyReached()) {
2.806 - dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
2.807 - typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
2.808 - original_edge.update();
2.809 - original_edge.set(f, e);
2.810 - residual_capacity.update();
2.811 - residual_capacity.set(f, res_graph.resCap(e));
2.812 - } else {
2.813 - if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
2.814 - typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]);
2.815 - original_edge.update();
2.816 - original_edge.set(f, e);
2.817 - residual_capacity.update();
2.818 - residual_capacity.set(f, res_graph.resCap(e));
2.819 - }
2.820 - }
2.821 - }
2.822 - ++bfs;
2.823 - } //computing distances from s in the residual graph
2.824 -
2.825 - bool __augment=true;
2.826 -
2.827 - while (__augment) {
2.828 - __augment=false;
2.829 - //computing blocking flow with dfs
2.830 - DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
2.831 - typename MG::template NodeMap<typename MG::Edge> pred(F);
2.832 - pred.set(sF, INVALID);
2.833 - //invalid iterators for sources
2.834 -
2.835 - typename MG::template NodeMap<Num> free(F);
2.836 -
2.837 - dfs.pushAndSetReached(sF);
2.838 - while (!dfs.finished()) {
2.839 - ++dfs;
2.840 - if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
2.841 - if (dfs.isBNodeNewlyReached()) {
2.842 - typename MG::Node v=F.aNode(dfs);
2.843 - typename MG::Node w=F.bNode(dfs);
2.844 - pred.set(w, dfs);
2.845 - if (F.valid(pred[v])) {
2.846 - free.set(w, std::min(free[v], residual_capacity[dfs]));
2.847 - } else {
2.848 - free.set(w, residual_capacity[dfs]);
2.849 - }
2.850 - if (w==tF) {
2.851 - __augment=true;
2.852 - _augment=true;
2.853 - break;
2.854 - }
2.855 -
2.856 - } else {
2.857 - F.erase(/*typename MG::OutEdgeIt*/(dfs));
2.858 - }
2.859 - }
2.860 - }
2.861 -
2.862 - if (__augment) {
2.863 - typename MG::Node n=tF;
2.864 - Num augment_value=free[tF];
2.865 - while (F.valid(pred[n])) {
2.866 - typename MG::Edge e=pred[n];
2.867 - res_graph.augment(original_edge[e], augment_value);
2.868 - n=F.tail(e);
2.869 - if (residual_capacity[e]==augment_value)
2.870 - F.erase(e);
2.871 - else
2.872 - residual_capacity.set(e, residual_capacity[e]-augment_value);
2.873 - }
2.874 - }
2.875 -
2.876 - }
2.877 -
2.878 - return _augment;
2.879 - }
2.880 -
2.881 -
2.882 -
2.883 -
2.884 -
2.885 -
2.886 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.887 - bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
2.888 - {
2.889 - bool _augment=false;
2.890 -
2.891 - ResGW res_graph(*g, *capacity, *flow);
2.892 -
2.893 - //ReachedMap level(res_graph);
2.894 - FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
2.895 - BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
2.896 -
2.897 - bfs.pushAndSetReached(s);
2.898 - DistanceMap<ResGW> dist(res_graph);
2.899 - while ( !bfs.finished() ) {
2.900 - ResGWOutEdgeIt e=bfs;
2.901 - if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
2.902 - dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
2.903 - }
2.904 - ++bfs;
2.905 - } //computing distances from s in the residual graph
2.906 -
2.907 - //Subgraph containing the edges on some shortest paths
2.908 - ConstMap<typename ResGW::Node, bool> true_map(true);
2.909 - typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
2.910 - DistanceMap<ResGW> > FilterResGW;
2.911 - FilterResGW filter_res_graph(res_graph, true_map, dist);
2.912 -
2.913 - //Subgraph, which is able to delete edges which are already
2.914 - //met by the dfs
2.915 - typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
2.916 - first_out_edges(filter_res_graph);
2.917 - typename FilterResGW::NodeIt v;
2.918 - for(filter_res_graph.first(v); filter_res_graph.valid(v);
2.919 - filter_res_graph.next(v))
2.920 - {
2.921 - typename FilterResGW::OutEdgeIt e;
2.922 - filter_res_graph.first(e, v);
2.923 - first_out_edges.set(v, e);
2.924 - }
2.925 - typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
2.926 - template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
2.927 - ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
2.928 -
2.929 - bool __augment=true;
2.930 -
2.931 - while (__augment) {
2.932 -
2.933 - __augment=false;
2.934 - //computing blocking flow with dfs
2.935 - DfsIterator< ErasingResGW,
2.936 - typename ErasingResGW::template NodeMap<bool> >
2.937 - dfs(erasing_res_graph);
2.938 - typename ErasingResGW::
2.939 - template NodeMap<typename ErasingResGW::OutEdgeIt>
2.940 - pred(erasing_res_graph);
2.941 - pred.set(s, INVALID);
2.942 - //invalid iterators for sources
2.943 -
2.944 - typename ErasingResGW::template NodeMap<Num>
2.945 - free1(erasing_res_graph);
2.946 -
2.947 - dfs.pushAndSetReached(
2.948 - typename ErasingResGW::Node(
2.949 - typename FilterResGW::Node(
2.950 - typename ResGW::Node(s)
2.951 - )
2.952 - )
2.953 - );
2.954 - while (!dfs.finished()) {
2.955 - ++dfs;
2.956 - if (erasing_res_graph.valid(
2.957 - typename ErasingResGW::OutEdgeIt(dfs)))
2.958 - {
2.959 - if (dfs.isBNodeNewlyReached()) {
2.960 -
2.961 - typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
2.962 - typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
2.963 -
2.964 - pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
2.965 - if (erasing_res_graph.valid(pred[v])) {
2.966 - free1.set(w, std::min(free1[v], res_graph.resCap(
2.967 - typename ErasingResGW::OutEdgeIt(dfs))));
2.968 - } else {
2.969 - free1.set(w, res_graph.resCap(
2.970 - typename ErasingResGW::OutEdgeIt(dfs)));
2.971 - }
2.972 -
2.973 - if (w==t) {
2.974 - __augment=true;
2.975 - _augment=true;
2.976 - break;
2.977 - }
2.978 - } else {
2.979 - erasing_res_graph.erase(dfs);
2.980 - }
2.981 - }
2.982 - }
2.983 -
2.984 - if (__augment) {
2.985 - typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t));
2.986 -// typename ResGW::NodeMap<Num> a(res_graph);
2.987 -// typename ResGW::Node b;
2.988 -// Num j=a[b];
2.989 -// typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
2.990 -// typename FilterResGW::Node b1;
2.991 -// Num j1=a1[b1];
2.992 -// typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
2.993 -// typename ErasingResGW::Node b2;
2.994 -// Num j2=a2[b2];
2.995 - Num augment_value=free1[n];
2.996 - while (erasing_res_graph.valid(pred[n])) {
2.997 - typename ErasingResGW::OutEdgeIt e=pred[n];
2.998 - res_graph.augment(e, augment_value);
2.999 - n=erasing_res_graph.tail(e);
2.1000 - if (res_graph.resCap(e)==0)
2.1001 - erasing_res_graph.erase(e);
2.1002 - }
2.1003 - }
2.1004 -
2.1005 - } //while (__augment)
2.1006 -
2.1007 - return _augment;
2.1008 - }
2.1009 -
2.1010 -
2.1011 -
2.1012 -
2.1013 -} //namespace hugo
2.1014 -
2.1015 -#endif //HUGO_PREFLOW_H
2.1016 -
2.1017 -
2.1018 -
2.1019 -