1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/hugo/max_flow.h Thu Jul 22 14:19:23 2004 +0000
1.3 @@ -0,0 +1,892 @@
1.4 +// -*- C++ -*-
1.5 +#ifndef HUGO_MAX_FLOW_NO_STACK_H
1.6 +#define HUGO_MAX_FLOW_NO_STACK_H
1.7 +
1.8 +#include <vector>
1.9 +#include <queue>
1.10 +//#include <stack>
1.11 +
1.12 +#include <hugo/graph_wrapper.h>
1.13 +#include <hugo/invalid.h>
1.14 +#include <hugo/maps.h>
1.15 +
1.16 +/// \file
1.17 +/// \brief The same as max_flow.h, but without using stl stack for the active nodes. Only for test.
1.18 +/// \ingroup galgs
1.19 +
1.20 +namespace hugo {
1.21 +
1.22 + /// \addtogroup galgs
1.23 + /// @{
1.24 + ///Maximum flow algorithms class.
1.25 +
1.26 + ///This class provides various algorithms for finding a flow of
1.27 + ///maximum value in a directed graph. The \e source node, the \e
1.28 + ///target node, the \e capacity of the edges and the \e starting \e
1.29 + ///flow value of the edges should be passed to the algorithm through the
1.30 + ///constructor. It is possible to change these quantities using the
1.31 + ///functions \ref resetSource, \ref resetTarget, \ref resetCap and
1.32 + ///\ref resetFlow. Before any subsequent runs of any algorithm of
1.33 + ///the class \ref resetFlow should be called.
1.34 +
1.35 + ///After running an algorithm of the class, the actual flow value
1.36 + ///can be obtained by calling \ref flowValue(). The minimum
1.37 + ///value cut can be written into a \c node map of \c bools by
1.38 + ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
1.39 + ///the inclusionwise minimum and maximum of the minimum value
1.40 + ///cuts, resp.)
1.41 + ///\param Graph The directed graph type the algorithm runs on.
1.42 + ///\param Num The number type of the capacities and the flow values.
1.43 + ///\param CapMap The capacity map type.
1.44 + ///\param FlowMap The flow map type.
1.45 + ///\author Marton Makai, Jacint Szabo
1.46 + template <typename Graph, typename Num,
1.47 + typename CapMap=typename Graph::template EdgeMap<Num>,
1.48 + typename FlowMap=typename Graph::template EdgeMap<Num> >
1.49 + class MaxFlow {
1.50 + protected:
1.51 + typedef typename Graph::Node Node;
1.52 + typedef typename Graph::NodeIt NodeIt;
1.53 + typedef typename Graph::EdgeIt EdgeIt;
1.54 + typedef typename Graph::OutEdgeIt OutEdgeIt;
1.55 + typedef typename Graph::InEdgeIt InEdgeIt;
1.56 +
1.57 + // typedef typename std::vector<std::stack<Node> > VecStack;
1.58 + typedef typename std::vector<Node> VecFirst;
1.59 + typedef typename Graph::template NodeMap<Node> NNMap;
1.60 + typedef typename std::vector<Node> VecNode;
1.61 +
1.62 + const Graph* g;
1.63 + Node s;
1.64 + Node t;
1.65 + const CapMap* capacity;
1.66 + FlowMap* flow;
1.67 + int n; //the number of nodes of G
1.68 + typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
1.69 + //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
1.70 + typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
1.71 + typedef typename ResGW::Edge ResGWEdge;
1.72 + //typedef typename ResGW::template NodeMap<bool> ReachedMap;
1.73 + typedef typename Graph::template NodeMap<int> ReachedMap;
1.74 +
1.75 +
1.76 + //level works as a bool map in augmenting path algorithms and is
1.77 + //used by bfs for storing reached information. In preflow, it
1.78 + //shows the levels of nodes.
1.79 + ReachedMap level;
1.80 +
1.81 + //excess is needed only in preflow
1.82 + typename Graph::template NodeMap<Num> excess;
1.83 +
1.84 + // constants used for heuristics
1.85 + static const int H0=20;
1.86 + static const int H1=1;
1.87 +
1.88 + public:
1.89 +
1.90 + ///Indicates the property of the starting flow.
1.91 +
1.92 + ///Indicates the property of the starting flow. The meanings are as follows:
1.93 + ///- \c ZERO_FLOW: constant zero flow
1.94 + ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
1.95 + ///the sum of the out-flows in every node except the \e source and
1.96 + ///the \e target.
1.97 + ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
1.98 + ///least the sum of the out-flows in every node except the \e source.
1.99 + ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be
1.100 + ///set to the constant zero flow in the beginning of the algorithm in this case.
1.101 + enum FlowEnum{
1.102 + ZERO_FLOW,
1.103 + GEN_FLOW,
1.104 + PRE_FLOW,
1.105 + NO_FLOW
1.106 + };
1.107 +
1.108 + enum StatusEnum {
1.109 + AFTER_NOTHING,
1.110 + AFTER_AUGMENTING,
1.111 + AFTER_FAST_AUGMENTING,
1.112 + AFTER_PRE_FLOW_PHASE_1,
1.113 + AFTER_PRE_FLOW_PHASE_2
1.114 + };
1.115 +
1.116 + /// Don not needle this flag only if necessary.
1.117 + StatusEnum status;
1.118 +
1.119 +// int number_of_augmentations;
1.120 +
1.121 +
1.122 +// template<typename IntMap>
1.123 +// class TrickyReachedMap {
1.124 +// protected:
1.125 +// IntMap* map;
1.126 +// int* number_of_augmentations;
1.127 +// public:
1.128 +// TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) :
1.129 +// map(&_map), number_of_augmentations(&_number_of_augmentations) { }
1.130 +// void set(const Node& n, bool b) {
1.131 +// if (b)
1.132 +// map->set(n, *number_of_augmentations);
1.133 +// else
1.134 +// map->set(n, *number_of_augmentations-1);
1.135 +// }
1.136 +// bool operator[](const Node& n) const {
1.137 +// return (*map)[n]==*number_of_augmentations;
1.138 +// }
1.139 +// };
1.140 +
1.141 + ///Constructor
1.142 +
1.143 + ///\todo Document, please.
1.144 + ///
1.145 + MaxFlow(const Graph& _G, Node _s, Node _t,
1.146 + const CapMap& _capacity, FlowMap& _flow) :
1.147 + g(&_G), s(_s), t(_t), capacity(&_capacity),
1.148 + flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0),
1.149 + status(AFTER_NOTHING) { }
1.150 +
1.151 + ///Runs a maximum flow algorithm.
1.152 +
1.153 + ///Runs a preflow algorithm, which is the fastest maximum flow
1.154 + ///algorithm up-to-date. The default for \c fe is ZERO_FLOW.
1.155 + ///\pre The starting flow must be
1.156 + /// - a constant zero flow if \c fe is \c ZERO_FLOW,
1.157 + /// - an arbitary flow if \c fe is \c GEN_FLOW,
1.158 + /// - an arbitary preflow if \c fe is \c PRE_FLOW,
1.159 + /// - any map if \c fe is NO_FLOW.
1.160 + void run(FlowEnum fe=ZERO_FLOW) {
1.161 + preflow(fe);
1.162 + }
1.163 +
1.164 +
1.165 + ///Runs a preflow algorithm.
1.166 +
1.167 + ///Runs a preflow algorithm. The preflow algorithms provide the
1.168 + ///fastest way to compute a maximum flow in a directed graph.
1.169 + ///\pre The starting flow must be
1.170 + /// - a constant zero flow if \c fe is \c ZERO_FLOW,
1.171 + /// - an arbitary flow if \c fe is \c GEN_FLOW,
1.172 + /// - an arbitary preflow if \c fe is \c PRE_FLOW,
1.173 + /// - any map if \c fe is NO_FLOW.
1.174 + ///
1.175 + ///\todo NO_FLOW should be the default flow.
1.176 + void preflow(FlowEnum fe) {
1.177 + preflowPhase1(fe);
1.178 + preflowPhase2();
1.179 + }
1.180 + // Heuristics:
1.181 + // 2 phase
1.182 + // gap
1.183 + // list 'level_list' on the nodes on level i implemented by hand
1.184 + // stack 'active' on the active nodes on level i
1.185 + // runs heuristic 'highest label' for H1*n relabels
1.186 + // runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
1.187 + // Parameters H0 and H1 are initialized to 20 and 1.
1.188 +
1.189 + ///Runs the first phase of the preflow algorithm.
1.190 +
1.191 + ///The preflow algorithm consists of two phases, this method runs the
1.192 + ///first phase. After the first phase the maximum flow value and a
1.193 + ///minimum value cut can already be computed, though a maximum flow
1.194 + ///is net yet obtained. So after calling this method \ref flowValue
1.195 + ///and \ref actMinCut gives proper results.
1.196 + ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
1.197 + ///give minimum value cuts unless calling \ref preflowPhase2.
1.198 + ///\pre The starting flow must be
1.199 + /// - a constant zero flow if \c fe is \c ZERO_FLOW,
1.200 + /// - an arbitary flow if \c fe is \c GEN_FLOW,
1.201 + /// - an arbitary preflow if \c fe is \c PRE_FLOW,
1.202 + /// - any map if \c fe is NO_FLOW.
1.203 + void preflowPhase1(FlowEnum fe)
1.204 + {
1.205 +
1.206 + int heur0=(int)(H0*n); //time while running 'bound decrease'
1.207 + int heur1=(int)(H1*n); //time while running 'highest label'
1.208 + int heur=heur1; //starting time interval (#of relabels)
1.209 + int numrelabel=0;
1.210 +
1.211 + bool what_heur=1;
1.212 + //It is 0 in case 'bound decrease' and 1 in case 'highest label'
1.213 +
1.214 + bool end=false;
1.215 + //Needed for 'bound decrease', true means no active nodes are above bound
1.216 + //b.
1.217 +
1.218 + int k=n-2; //bound on the highest level under n containing a node
1.219 + int b=k; //bound on the highest level under n of an active node
1.220 +
1.221 + VecFirst first(n, INVALID);
1.222 + NNMap next(*g, INVALID); //maybe INVALID is not needed
1.223 + // VecStack active(n);
1.224 +
1.225 + NNMap left(*g, INVALID);
1.226 + NNMap right(*g, INVALID);
1.227 + VecNode level_list(n,INVALID);
1.228 + //List of the nodes in level i<n, set to n.
1.229 +
1.230 + NodeIt v;
1.231 + for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
1.232 + //setting each node to level n
1.233 +
1.234 + if ( fe == NO_FLOW ) {
1.235 + EdgeIt e;
1.236 + for(g->first(e); g->valid(e); g->next(e)) flow->set(e,0);
1.237 + }
1.238 +
1.239 + switch (fe) { //computing the excess
1.240 + case PRE_FLOW:
1.241 + {
1.242 + NodeIt v;
1.243 + for(g->first(v); g->valid(v); g->next(v)) {
1.244 + Num exc=0;
1.245 +
1.246 + InEdgeIt e;
1.247 + for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
1.248 + OutEdgeIt f;
1.249 + for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
1.250 +
1.251 + excess.set(v,exc);
1.252 +
1.253 + //putting the active nodes into the stack
1.254 + int lev=level[v];
1.255 + if ( exc > 0 && lev < n && v != t )
1.256 + {
1.257 + next.set(v,first[lev]);
1.258 + first[lev]=v;
1.259 + }
1.260 + // active[lev].push(v);
1.261 + }
1.262 + break;
1.263 + }
1.264 + case GEN_FLOW:
1.265 + {
1.266 + NodeIt v;
1.267 + for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
1.268 +
1.269 + Num exc=0;
1.270 + InEdgeIt e;
1.271 + for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
1.272 + OutEdgeIt f;
1.273 + for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
1.274 + excess.set(t,exc);
1.275 + break;
1.276 + }
1.277 + case ZERO_FLOW:
1.278 + case NO_FLOW:
1.279 + {
1.280 + NodeIt v;
1.281 + for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
1.282 + break;
1.283 + }
1.284 + }
1.285 +
1.286 + preflowPreproc(fe, next, first,/*active*/ level_list, left, right);
1.287 + //End of preprocessing
1.288 +
1.289 +
1.290 + //Push/relabel on the highest level active nodes.
1.291 + while ( true ) {
1.292 + if ( b == 0 ) {
1.293 + if ( !what_heur && !end && k > 0 ) {
1.294 + b=k;
1.295 + end=true;
1.296 + } else break;
1.297 + }
1.298 +
1.299 + if ( !g->valid(first[b])/*active[b].empty()*/ ) --b;
1.300 + else {
1.301 + end=false;
1.302 + Node w=first[b];
1.303 + first[b]=next[w];
1.304 + /* Node w=active[b].top();
1.305 + active[b].pop();*/
1.306 + int newlevel=push(w,/*active*/next, first);
1.307 + if ( excess[w] > 0 ) relabel(w, newlevel, /*active*/next, first, level_list,
1.308 + left, right, b, k, what_heur);
1.309 +
1.310 + ++numrelabel;
1.311 + if ( numrelabel >= heur ) {
1.312 + numrelabel=0;
1.313 + if ( what_heur ) {
1.314 + what_heur=0;
1.315 + heur=heur0;
1.316 + end=false;
1.317 + } else {
1.318 + what_heur=1;
1.319 + heur=heur1;
1.320 + b=k;
1.321 + }
1.322 + }
1.323 + }
1.324 + }
1.325 +
1.326 + status=AFTER_PRE_FLOW_PHASE_1;
1.327 + }
1.328 +
1.329 +
1.330 + ///Runs the second phase of the preflow algorithm.
1.331 +
1.332 + ///The preflow algorithm consists of two phases, this method runs
1.333 + ///the second phase. After calling \ref preflowPhase1 and then
1.334 + ///\ref preflowPhase2 the methods \ref flowValue, \ref minCut,
1.335 + ///\ref minMinCut and \ref maxMinCut give proper results.
1.336 + ///\pre \ref preflowPhase1 must be called before.
1.337 + void preflowPhase2()
1.338 + {
1.339 +
1.340 + int k=n-2; //bound on the highest level under n containing a node
1.341 + int b=k; //bound on the highest level under n of an active node
1.342 +
1.343 +
1.344 + VecFirst first(n, INVALID);
1.345 + NNMap next(*g, INVALID); //maybe INVALID is not needed
1.346 + // VecStack active(n);
1.347 + level.set(s,0);
1.348 + std::queue<Node> bfs_queue;
1.349 + bfs_queue.push(s);
1.350 +
1.351 + while (!bfs_queue.empty()) {
1.352 +
1.353 + Node v=bfs_queue.front();
1.354 + bfs_queue.pop();
1.355 + int l=level[v]+1;
1.356 +
1.357 + InEdgeIt e;
1.358 + for(g->first(e,v); g->valid(e); g->next(e)) {
1.359 + if ( (*capacity)[e] <= (*flow)[e] ) continue;
1.360 + Node u=g->tail(e);
1.361 + if ( level[u] >= n ) {
1.362 + bfs_queue.push(u);
1.363 + level.set(u, l);
1.364 + if ( excess[u] > 0 ) {
1.365 + next.set(u,first[l]);
1.366 + first[l]=u;
1.367 + //active[l].push(u);
1.368 + }
1.369 + }
1.370 + }
1.371 +
1.372 + OutEdgeIt f;
1.373 + for(g->first(f,v); g->valid(f); g->next(f)) {
1.374 + if ( 0 >= (*flow)[f] ) continue;
1.375 + Node u=g->head(f);
1.376 + if ( level[u] >= n ) {
1.377 + bfs_queue.push(u);
1.378 + level.set(u, l);
1.379 + if ( excess[u] > 0 ) {
1.380 + next.set(u,first[l]);
1.381 + first[l]=u;
1.382 + //active[l].push(u);
1.383 + }
1.384 + }
1.385 + }
1.386 + }
1.387 + b=n-2;
1.388 +
1.389 + while ( true ) {
1.390 +
1.391 + if ( b == 0 ) break;
1.392 +
1.393 + if ( !g->valid(first[b])/*active[b].empty()*/ ) --b;
1.394 + else {
1.395 +
1.396 + Node w=first[b];
1.397 + first[b]=next[w];
1.398 + /* Node w=active[b].top();
1.399 + active[b].pop();*/
1.400 + int newlevel=push(w,next, first/*active*/);
1.401 +
1.402 + //relabel
1.403 + if ( excess[w] > 0 ) {
1.404 + level.set(w,++newlevel);
1.405 + next.set(w,first[newlevel]);
1.406 + first[newlevel]=w;
1.407 + //active[newlevel].push(w);
1.408 + b=newlevel;
1.409 + }
1.410 + } // if stack[b] is nonempty
1.411 + } // while(true)
1.412 +
1.413 + status=AFTER_PRE_FLOW_PHASE_2;
1.414 + }
1.415 +
1.416 +
1.417 + /// Returns the maximum value of a flow.
1.418 +
1.419 + /// Returns the maximum value of a flow, by counting the
1.420 + /// over-flow of the target node \ref t.
1.421 + /// It can be called already after running \ref preflowPhase1.
1.422 + Num flowValue() const {
1.423 + Num a=0;
1.424 + for(InEdgeIt e(*g,t);g->valid(e);G.next(e)) a+=(*flow)[e];
1.425 + for(OutEdgeIt e(*g,t);g->valid(e);G.next(e)) a-=(*flow)[e];
1.426 +
1.427 + //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan
1.428 + }
1.429 +
1.430 + ///Returns a minimum value cut after calling \ref preflowPhase1.
1.431 +
1.432 + ///After the first phase of the preflow algorithm the maximum flow
1.433 + ///value and a minimum value cut can already be computed. This
1.434 + ///method can be called after running \ref preflowPhase1 for
1.435 + ///obtaining a minimum value cut.
1.436 + /// \warning Gives proper result only right after calling \ref
1.437 + /// preflowPhase1.
1.438 + /// \todo We have to make some status variable which shows the
1.439 + /// actual state
1.440 + /// of the class. This enables us to determine which methods are valid
1.441 + /// for MinCut computation
1.442 + template<typename _CutMap>
1.443 + void actMinCut(_CutMap& M) const {
1.444 + NodeIt v;
1.445 + switch (status) {
1.446 + case AFTER_PRE_FLOW_PHASE_1:
1.447 + for(g->first(v); g->valid(v); g->next(v)) {
1.448 + if (level[v] < n) {
1.449 + M.set(v, false);
1.450 + } else {
1.451 + M.set(v, true);
1.452 + }
1.453 + }
1.454 + break;
1.455 + case AFTER_PRE_FLOW_PHASE_2:
1.456 + case AFTER_NOTHING:
1.457 + minMinCut(M);
1.458 + break;
1.459 + case AFTER_AUGMENTING:
1.460 + for(g->first(v); g->valid(v); g->next(v)) {
1.461 + if (level[v]) {
1.462 + M.set(v, true);
1.463 + } else {
1.464 + M.set(v, false);
1.465 + }
1.466 + }
1.467 + break;
1.468 + case AFTER_FAST_AUGMENTING:
1.469 + for(g->first(v); g->valid(v); g->next(v)) {
1.470 + if (level[v]==number_of_augmentations) {
1.471 + M.set(v, true);
1.472 + } else {
1.473 + M.set(v, false);
1.474 + }
1.475 + }
1.476 + break;
1.477 + }
1.478 + }
1.479 +
1.480 + ///Returns the inclusionwise minimum of the minimum value cuts.
1.481 +
1.482 + ///Sets \c M to the characteristic vector of the minimum value cut
1.483 + ///which is inclusionwise minimum. It is computed by processing
1.484 + ///a bfs from the source node \c s in the residual graph.
1.485 + ///\pre M should be a node map of bools initialized to false.
1.486 + ///\pre \c flow must be a maximum flow.
1.487 + template<typename _CutMap>
1.488 + void minMinCut(_CutMap& M) const {
1.489 + std::queue<Node> queue;
1.490 +
1.491 + M.set(s,true);
1.492 + queue.push(s);
1.493 +
1.494 + while (!queue.empty()) {
1.495 + Node w=queue.front();
1.496 + queue.pop();
1.497 +
1.498 + OutEdgeIt e;
1.499 + for(g->first(e,w) ; g->valid(e); g->next(e)) {
1.500 + Node v=g->head(e);
1.501 + if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
1.502 + queue.push(v);
1.503 + M.set(v, true);
1.504 + }
1.505 + }
1.506 +
1.507 + InEdgeIt f;
1.508 + for(g->first(f,w) ; g->valid(f); g->next(f)) {
1.509 + Node v=g->tail(f);
1.510 + if (!M[v] && (*flow)[f] > 0 ) {
1.511 + queue.push(v);
1.512 + M.set(v, true);
1.513 + }
1.514 + }
1.515 + }
1.516 + }
1.517 +
1.518 + ///Returns the inclusionwise maximum of the minimum value cuts.
1.519 +
1.520 + ///Sets \c M to the characteristic vector of the minimum value cut
1.521 + ///which is inclusionwise maximum. It is computed by processing a
1.522 + ///backward bfs from the target node \c t in the residual graph.
1.523 + ///\pre M should be a node map of bools initialized to false.
1.524 + ///\pre \c flow must be a maximum flow.
1.525 + template<typename _CutMap>
1.526 + void maxMinCut(_CutMap& M) const {
1.527 +
1.528 + NodeIt v;
1.529 + for(g->first(v) ; g->valid(v); g->next(v)) {
1.530 + M.set(v, true);
1.531 + }
1.532 +
1.533 + std::queue<Node> queue;
1.534 +
1.535 + M.set(t,false);
1.536 + queue.push(t);
1.537 +
1.538 + while (!queue.empty()) {
1.539 + Node w=queue.front();
1.540 + queue.pop();
1.541 +
1.542 + InEdgeIt e;
1.543 + for(g->first(e,w) ; g->valid(e); g->next(e)) {
1.544 + Node v=g->tail(e);
1.545 + if (M[v] && (*flow)[e] < (*capacity)[e] ) {
1.546 + queue.push(v);
1.547 + M.set(v, false);
1.548 + }
1.549 + }
1.550 +
1.551 + OutEdgeIt f;
1.552 + for(g->first(f,w) ; g->valid(f); g->next(f)) {
1.553 + Node v=g->head(f);
1.554 + if (M[v] && (*flow)[f] > 0 ) {
1.555 + queue.push(v);
1.556 + M.set(v, false);
1.557 + }
1.558 + }
1.559 + }
1.560 + }
1.561 +
1.562 + ///Returns a minimum value cut.
1.563 +
1.564 + ///Sets \c M to the characteristic vector of a minimum value cut.
1.565 + ///\pre M should be a node map of bools initialized to false.
1.566 + ///\pre \c flow must be a maximum flow.
1.567 + template<typename CutMap>
1.568 + void minCut(CutMap& M) const { minMinCut(M); }
1.569 +
1.570 + ///Resets the source node to \c _s.
1.571 +
1.572 + ///Resets the source node to \c _s.
1.573 + ///
1.574 + void resetSource(Node _s) { s=_s; status=AFTER_NOTHING; }
1.575 +
1.576 + ///Resets the target node to \c _t.
1.577 +
1.578 + ///Resets the target node to \c _t.
1.579 + ///
1.580 + void resetTarget(Node _t) { t=_t; status=AFTER_NOTHING; }
1.581 +
1.582 + /// Resets the edge map of the capacities to _cap.
1.583 +
1.584 + /// Resets the edge map of the capacities to _cap.
1.585 + ///
1.586 + void resetCap(const CapMap& _cap)
1.587 + { capacity=&_cap; status=AFTER_NOTHING; }
1.588 +
1.589 + /// Resets the edge map of the flows to _flow.
1.590 +
1.591 + /// Resets the edge map of the flows to _flow.
1.592 + ///
1.593 + void resetFlow(FlowMap& _flow) { flow=&_flow; status=AFTER_NOTHING; }
1.594 +
1.595 +
1.596 + private:
1.597 +
1.598 + int push(Node w, NNMap& next, VecFirst& first) {
1.599 +
1.600 + int lev=level[w];
1.601 + Num exc=excess[w];
1.602 + int newlevel=n; //bound on the next level of w
1.603 +
1.604 + OutEdgeIt e;
1.605 + for(g->first(e,w); g->valid(e); g->next(e)) {
1.606 +
1.607 + if ( (*flow)[e] >= (*capacity)[e] ) continue;
1.608 + Node v=g->head(e);
1.609 +
1.610 + if( lev > level[v] ) { //Push is allowed now
1.611 +
1.612 + if ( excess[v]<=0 && v!=t && v!=s ) {
1.613 + next.set(v,first[level[v]]);
1.614 + first[level[v]]=v;
1.615 + // int lev_v=level[v];
1.616 + //active[lev_v].push(v);
1.617 + }
1.618 +
1.619 + Num cap=(*capacity)[e];
1.620 + Num flo=(*flow)[e];
1.621 + Num remcap=cap-flo;
1.622 +
1.623 + if ( remcap >= exc ) { //A nonsaturating push.
1.624 +
1.625 + flow->set(e, flo+exc);
1.626 + excess.set(v, excess[v]+exc);
1.627 + exc=0;
1.628 + break;
1.629 +
1.630 + } else { //A saturating push.
1.631 + flow->set(e, cap);
1.632 + excess.set(v, excess[v]+remcap);
1.633 + exc-=remcap;
1.634 + }
1.635 + } else if ( newlevel > level[v] ) newlevel = level[v];
1.636 + } //for out edges wv
1.637 +
1.638 + if ( exc > 0 ) {
1.639 + InEdgeIt e;
1.640 + for(g->first(e,w); g->valid(e); g->next(e)) {
1.641 +
1.642 + if( (*flow)[e] <= 0 ) continue;
1.643 + Node v=g->tail(e);
1.644 +
1.645 + if( lev > level[v] ) { //Push is allowed now
1.646 +
1.647 + if ( excess[v]<=0 && v!=t && v!=s ) {
1.648 + next.set(v,first[level[v]]);
1.649 + first[level[v]]=v;
1.650 + //int lev_v=level[v];
1.651 + //active[lev_v].push(v);
1.652 + }
1.653 +
1.654 + Num flo=(*flow)[e];
1.655 +
1.656 + if ( flo >= exc ) { //A nonsaturating push.
1.657 +
1.658 + flow->set(e, flo-exc);
1.659 + excess.set(v, excess[v]+exc);
1.660 + exc=0;
1.661 + break;
1.662 + } else { //A saturating push.
1.663 +
1.664 + excess.set(v, excess[v]+flo);
1.665 + exc-=flo;
1.666 + flow->set(e,0);
1.667 + }
1.668 + } else if ( newlevel > level[v] ) newlevel = level[v];
1.669 + } //for in edges vw
1.670 +
1.671 + } // if w still has excess after the out edge for cycle
1.672 +
1.673 + excess.set(w, exc);
1.674 +
1.675 + return newlevel;
1.676 + }
1.677 +
1.678 +
1.679 + void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
1.680 + VecNode& level_list, NNMap& left, NNMap& right)
1.681 + {
1.682 + std::queue<Node> bfs_queue;
1.683 +
1.684 + switch (fe) {
1.685 + case NO_FLOW: //flow is already set to const zero in this case
1.686 + case ZERO_FLOW:
1.687 + {
1.688 + //Reverse_bfs from t, to find the starting level.
1.689 + level.set(t,0);
1.690 + bfs_queue.push(t);
1.691 +
1.692 + while (!bfs_queue.empty()) {
1.693 +
1.694 + Node v=bfs_queue.front();
1.695 + bfs_queue.pop();
1.696 + int l=level[v]+1;
1.697 +
1.698 + InEdgeIt e;
1.699 + for(g->first(e,v); g->valid(e); g->next(e)) {
1.700 + Node w=g->tail(e);
1.701 + if ( level[w] == n && w != s ) {
1.702 + bfs_queue.push(w);
1.703 + Node z=level_list[l];
1.704 + if ( g->valid(z) ) left.set(z,w);
1.705 + right.set(w,z);
1.706 + level_list[l]=w;
1.707 + level.set(w, l);
1.708 + }
1.709 + }
1.710 + }
1.711 +
1.712 + //the starting flow
1.713 + OutEdgeIt e;
1.714 + for(g->first(e,s); g->valid(e); g->next(e))
1.715 + {
1.716 + Num c=(*capacity)[e];
1.717 + if ( c <= 0 ) continue;
1.718 + Node w=g->head(e);
1.719 + if ( level[w] < n ) {
1.720 + if ( excess[w] <= 0 && w!=t )
1.721 + {
1.722 + next.set(w,first[level[w]]);
1.723 + first[level[w]]=w;
1.724 + //active[level[w]].push(w);
1.725 + }
1.726 + flow->set(e, c);
1.727 + excess.set(w, excess[w]+c);
1.728 + }
1.729 + }
1.730 + break;
1.731 + }
1.732 +
1.733 + case GEN_FLOW:
1.734 + case PRE_FLOW:
1.735 + {
1.736 + //Reverse_bfs from t in the residual graph,
1.737 + //to find the starting level.
1.738 + level.set(t,0);
1.739 + bfs_queue.push(t);
1.740 +
1.741 + while (!bfs_queue.empty()) {
1.742 +
1.743 + Node v=bfs_queue.front();
1.744 + bfs_queue.pop();
1.745 + int l=level[v]+1;
1.746 +
1.747 + InEdgeIt e;
1.748 + for(g->first(e,v); g->valid(e); g->next(e)) {
1.749 + if ( (*capacity)[e] <= (*flow)[e] ) continue;
1.750 + Node w=g->tail(e);
1.751 + if ( level[w] == n && w != s ) {
1.752 + bfs_queue.push(w);
1.753 + Node z=level_list[l];
1.754 + if ( g->valid(z) ) left.set(z,w);
1.755 + right.set(w,z);
1.756 + level_list[l]=w;
1.757 + level.set(w, l);
1.758 + }
1.759 + }
1.760 +
1.761 + OutEdgeIt f;
1.762 + for(g->first(f,v); g->valid(f); g->next(f)) {
1.763 + if ( 0 >= (*flow)[f] ) continue;
1.764 + Node w=g->head(f);
1.765 + if ( level[w] == n && w != s ) {
1.766 + bfs_queue.push(w);
1.767 + Node z=level_list[l];
1.768 + if ( g->valid(z) ) left.set(z,w);
1.769 + right.set(w,z);
1.770 + level_list[l]=w;
1.771 + level.set(w, l);
1.772 + }
1.773 + }
1.774 + }
1.775 +
1.776 +
1.777 + //the starting flow
1.778 + OutEdgeIt e;
1.779 + for(g->first(e,s); g->valid(e); g->next(e))
1.780 + {
1.781 + Num rem=(*capacity)[e]-(*flow)[e];
1.782 + if ( rem <= 0 ) continue;
1.783 + Node w=g->head(e);
1.784 + if ( level[w] < n ) {
1.785 + if ( excess[w] <= 0 && w!=t )
1.786 + {
1.787 + next.set(w,first[level[w]]);
1.788 + first[level[w]]=w;
1.789 + //active[level[w]].push(w);
1.790 + }
1.791 + flow->set(e, (*capacity)[e]);
1.792 + excess.set(w, excess[w]+rem);
1.793 + }
1.794 + }
1.795 +
1.796 + InEdgeIt f;
1.797 + for(g->first(f,s); g->valid(f); g->next(f))
1.798 + {
1.799 + if ( (*flow)[f] <= 0 ) continue;
1.800 + Node w=g->tail(f);
1.801 + if ( level[w] < n ) {
1.802 + if ( excess[w] <= 0 && w!=t )
1.803 + {
1.804 + next.set(w,first[level[w]]);
1.805 + first[level[w]]=w;
1.806 + //active[level[w]].push(w);
1.807 + }
1.808 + excess.set(w, excess[w]+(*flow)[f]);
1.809 + flow->set(f, 0);
1.810 + }
1.811 + }
1.812 + break;
1.813 + } //case PRE_FLOW
1.814 + }
1.815 + } //preflowPreproc
1.816 +
1.817 +
1.818 +
1.819 + void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
1.820 + VecNode& level_list, NNMap& left,
1.821 + NNMap& right, int& b, int& k, bool what_heur )
1.822 + {
1.823 +
1.824 + Num lev=level[w];
1.825 +
1.826 + Node right_n=right[w];
1.827 + Node left_n=left[w];
1.828 +
1.829 + //unlacing starts
1.830 + if ( g->valid(right_n) ) {
1.831 + if ( g->valid(left_n) ) {
1.832 + right.set(left_n, right_n);
1.833 + left.set(right_n, left_n);
1.834 + } else {
1.835 + level_list[lev]=right_n;
1.836 + left.set(right_n, INVALID);
1.837 + }
1.838 + } else {
1.839 + if ( g->valid(left_n) ) {
1.840 + right.set(left_n, INVALID);
1.841 + } else {
1.842 + level_list[lev]=INVALID;
1.843 + }
1.844 + }
1.845 + //unlacing ends
1.846 +
1.847 + if ( !g->valid(level_list[lev]) ) {
1.848 +
1.849 + //gapping starts
1.850 + for (int i=lev; i!=k ; ) {
1.851 + Node v=level_list[++i];
1.852 + while ( g->valid(v) ) {
1.853 + level.set(v,n);
1.854 + v=right[v];
1.855 + }
1.856 + level_list[i]=INVALID;
1.857 + if ( !what_heur ) first[i]=INVALID;
1.858 + /*{
1.859 + while ( !active[i].empty() ) {
1.860 + active[i].pop(); //FIXME: ezt szebben kene
1.861 + }
1.862 + }*/
1.863 + }
1.864 +
1.865 + level.set(w,n);
1.866 + b=lev-1;
1.867 + k=b;
1.868 + //gapping ends
1.869 +
1.870 + } else {
1.871 +
1.872 + if ( newlevel == n ) level.set(w,n);
1.873 + else {
1.874 + level.set(w,++newlevel);
1.875 + next.set(w,first[newlevel]);
1.876 + first[newlevel]=w;
1.877 + // active[newlevel].push(w);
1.878 + if ( what_heur ) b=newlevel;
1.879 + if ( k < newlevel ) ++k; //now k=newlevel
1.880 + Node z=level_list[newlevel];
1.881 + if ( g->valid(z) ) left.set(z,w);
1.882 + right.set(w,z);
1.883 + left.set(w,INVALID);
1.884 + level_list[newlevel]=w;
1.885 + }
1.886 + }
1.887 + } //relabel
1.888 + }; //class MaxFlow
1.889 +} //namespace hugo
1.890 +
1.891 +#endif //HUGO_MAX_FLOW_H
1.892 +
1.893 +
1.894 +
1.895 +
2.1 --- a/src/work/jacint/max_flow_no_stack.h Thu Jul 22 14:09:21 2004 +0000
2.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
2.3 @@ -1,1318 +0,0 @@
2.4 -// -*- C++ -*-
2.5 -#ifndef HUGO_MAX_FLOW_NO_STACK_H
2.6 -#define HUGO_MAX_FLOW_NO_STACK_H
2.7 -
2.8 -#include <vector>
2.9 -#include <queue>
2.10 -//#include <stack>
2.11 -
2.12 -#include <hugo/graph_wrapper.h>
2.13 -#include <bfs_dfs.h>
2.14 -#include <hugo/invalid.h>
2.15 -#include <hugo/maps.h>
2.16 -#include <hugo/for_each_macros.h>
2.17 -
2.18 -/// \file
2.19 -/// \brief The same as max_flow.h, but without using stl stack for the active nodes. Only for test.
2.20 -/// \ingroup galgs
2.21 -
2.22 -namespace hugo {
2.23 -
2.24 - /// \addtogroup galgs
2.25 - /// @{
2.26 - ///Maximum flow algorithms class.
2.27 -
2.28 - ///This class provides various algorithms for finding a flow of
2.29 - ///maximum value in a directed graph. The \e source node, the \e
2.30 - ///target node, the \e capacity of the edges and the \e starting \e
2.31 - ///flow value of the edges should be passed to the algorithm through the
2.32 - ///constructor. It is possible to change these quantities using the
2.33 - ///functions \ref resetSource, \ref resetTarget, \ref resetCap and
2.34 - ///\ref resetFlow. Before any subsequent runs of any algorithm of
2.35 - ///the class \ref resetFlow should be called.
2.36 -
2.37 - ///After running an algorithm of the class, the actual flow value
2.38 - ///can be obtained by calling \ref flowValue(). The minimum
2.39 - ///value cut can be written into a \c node map of \c bools by
2.40 - ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
2.41 - ///the inclusionwise minimum and maximum of the minimum value
2.42 - ///cuts, resp.)
2.43 - ///\param Graph The directed graph type the algorithm runs on.
2.44 - ///\param Num The number type of the capacities and the flow values.
2.45 - ///\param CapMap The capacity map type.
2.46 - ///\param FlowMap The flow map type.
2.47 - ///\author Marton Makai, Jacint Szabo
2.48 - template <typename Graph, typename Num,
2.49 - typename CapMap=typename Graph::template EdgeMap<Num>,
2.50 - typename FlowMap=typename Graph::template EdgeMap<Num> >
2.51 - class MaxFlowNoStack {
2.52 - protected:
2.53 - typedef typename Graph::Node Node;
2.54 - typedef typename Graph::NodeIt NodeIt;
2.55 - typedef typename Graph::EdgeIt EdgeIt;
2.56 - typedef typename Graph::OutEdgeIt OutEdgeIt;
2.57 - typedef typename Graph::InEdgeIt InEdgeIt;
2.58 -
2.59 - // typedef typename std::vector<std::stack<Node> > VecStack;
2.60 - typedef typename std::vector<Node> VecFirst;
2.61 - typedef typename Graph::template NodeMap<Node> NNMap;
2.62 - typedef typename std::vector<Node> VecNode;
2.63 -
2.64 - const Graph* g;
2.65 - Node s;
2.66 - Node t;
2.67 - const CapMap* capacity;
2.68 - FlowMap* flow;
2.69 - int n; //the number of nodes of G
2.70 - typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
2.71 - //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
2.72 - typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
2.73 - typedef typename ResGW::Edge ResGWEdge;
2.74 - //typedef typename ResGW::template NodeMap<bool> ReachedMap;
2.75 - typedef typename Graph::template NodeMap<int> ReachedMap;
2.76 -
2.77 -
2.78 - //level works as a bool map in augmenting path algorithms and is
2.79 - //used by bfs for storing reached information. In preflow, it
2.80 - //shows the levels of nodes.
2.81 - ReachedMap level;
2.82 -
2.83 - //excess is needed only in preflow
2.84 - typename Graph::template NodeMap<Num> excess;
2.85 -
2.86 - //fixme
2.87 -// protected:
2.88 - // MaxFlow() { }
2.89 - // void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
2.90 - // FlowMap& _flow)
2.91 - // {
2.92 - // g=&_G;
2.93 - // s=_s;
2.94 - // t=_t;
2.95 - // capacity=&_capacity;
2.96 - // flow=&_flow;
2.97 - // n=_G.nodeNum;
2.98 - // level.set (_G); //kellene vmi ilyesmi fv
2.99 - // excess(_G,0); //itt is
2.100 - // }
2.101 -
2.102 - // constants used for heuristics
2.103 - static const int H0=20;
2.104 - static const int H1=1;
2.105 -
2.106 - public:
2.107 -
2.108 - ///Indicates the property of the starting flow.
2.109 -
2.110 - ///Indicates the property of the starting flow. The meanings are as follows:
2.111 - ///- \c ZERO_FLOW: constant zero flow
2.112 - ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
2.113 - ///the sum of the out-flows in every node except the \e source and
2.114 - ///the \e target.
2.115 - ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
2.116 - ///least the sum of the out-flows in every node except the \e source.
2.117 - ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be
2.118 - ///set to the constant zero flow in the beginning of the algorithm in this case.
2.119 - enum FlowEnum{
2.120 - ZERO_FLOW,
2.121 - GEN_FLOW,
2.122 - PRE_FLOW,
2.123 - NO_FLOW
2.124 - };
2.125 -
2.126 - enum StatusEnum {
2.127 - AFTER_NOTHING,
2.128 - AFTER_AUGMENTING,
2.129 - AFTER_FAST_AUGMENTING,
2.130 - AFTER_PRE_FLOW_PHASE_1,
2.131 - AFTER_PRE_FLOW_PHASE_2
2.132 - };
2.133 -
2.134 - /// Don not needle this flag only if necessary.
2.135 - StatusEnum status;
2.136 - int number_of_augmentations;
2.137 -
2.138 -
2.139 - template<typename IntMap>
2.140 - class TrickyReachedMap {
2.141 - protected:
2.142 - IntMap* map;
2.143 - int* number_of_augmentations;
2.144 - public:
2.145 - TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) :
2.146 - map(&_map), number_of_augmentations(&_number_of_augmentations) { }
2.147 - void set(const Node& n, bool b) {
2.148 - if (b)
2.149 - map->set(n, *number_of_augmentations);
2.150 - else
2.151 - map->set(n, *number_of_augmentations-1);
2.152 - }
2.153 - bool operator[](const Node& n) const {
2.154 - return (*map)[n]==*number_of_augmentations;
2.155 - }
2.156 - };
2.157 -
2.158 - ///Constructor
2.159 -
2.160 - ///\todo Document, please.
2.161 - ///
2.162 - MaxFlowNoStack(const Graph& _G, Node _s, Node _t,
2.163 - const CapMap& _capacity, FlowMap& _flow) :
2.164 - g(&_G), s(_s), t(_t), capacity(&_capacity),
2.165 - flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0),
2.166 - status(AFTER_NOTHING), number_of_augmentations(0) { }
2.167 -
2.168 - ///Runs a maximum flow algorithm.
2.169 -
2.170 - ///Runs a preflow algorithm, which is the fastest maximum flow
2.171 - ///algorithm up-to-date. The default for \c fe is ZERO_FLOW.
2.172 - ///\pre The starting flow must be
2.173 - /// - a constant zero flow if \c fe is \c ZERO_FLOW,
2.174 - /// - an arbitary flow if \c fe is \c GEN_FLOW,
2.175 - /// - an arbitary preflow if \c fe is \c PRE_FLOW,
2.176 - /// - any map if \c fe is NO_FLOW.
2.177 - void run(FlowEnum fe=ZERO_FLOW) {
2.178 - preflow(fe);
2.179 - }
2.180 -
2.181 -
2.182 - ///Runs a preflow algorithm.
2.183 -
2.184 - ///Runs a preflow algorithm. The preflow algorithms provide the
2.185 - ///fastest way to compute a maximum flow in a directed graph.
2.186 - ///\pre The starting flow must be
2.187 - /// - a constant zero flow if \c fe is \c ZERO_FLOW,
2.188 - /// - an arbitary flow if \c fe is \c GEN_FLOW,
2.189 - /// - an arbitary preflow if \c fe is \c PRE_FLOW,
2.190 - /// - any map if \c fe is NO_FLOW.
2.191 - ///
2.192 - ///\todo NO_FLOW should be the default flow.
2.193 - void preflow(FlowEnum fe) {
2.194 - preflowPhase1(fe);
2.195 - preflowPhase2();
2.196 - }
2.197 - // Heuristics:
2.198 - // 2 phase
2.199 - // gap
2.200 - // list 'level_list' on the nodes on level i implemented by hand
2.201 - // stack 'active' on the active nodes on level i
2.202 - // runs heuristic 'highest label' for H1*n relabels
2.203 - // runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
2.204 - // Parameters H0 and H1 are initialized to 20 and 1.
2.205 -
2.206 - ///Runs the first phase of the preflow algorithm.
2.207 -
2.208 - ///The preflow algorithm consists of two phases, this method runs the
2.209 - ///first phase. After the first phase the maximum flow value and a
2.210 - ///minimum value cut can already be computed, though a maximum flow
2.211 - ///is net yet obtained. So after calling this method \ref flowValue
2.212 - ///and \ref actMinCut gives proper results.
2.213 - ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
2.214 - ///give minimum value cuts unless calling \ref preflowPhase2.
2.215 - ///\pre The starting flow must be
2.216 - /// - a constant zero flow if \c fe is \c ZERO_FLOW,
2.217 - /// - an arbitary flow if \c fe is \c GEN_FLOW,
2.218 - /// - an arbitary preflow if \c fe is \c PRE_FLOW,
2.219 - /// - any map if \c fe is NO_FLOW.
2.220 - void preflowPhase1(FlowEnum fe);
2.221 -
2.222 - ///Runs the second phase of the preflow algorithm.
2.223 -
2.224 - ///The preflow algorithm consists of two phases, this method runs
2.225 - ///the second phase. After calling \ref preflowPhase1 and then
2.226 - ///\ref preflowPhase2 the methods \ref flowValue, \ref minCut,
2.227 - ///\ref minMinCut and \ref maxMinCut give proper results.
2.228 - ///\pre \ref preflowPhase1 must be called before.
2.229 - void preflowPhase2();
2.230 -
2.231 - /// Starting from a flow, this method searches for an augmenting path
2.232 - /// according to the Edmonds-Karp algorithm
2.233 - /// and augments the flow on if any.
2.234 - /// The return value shows if the augmentation was succesful.
2.235 - bool augmentOnShortestPath();
2.236 - bool augmentOnShortestPath2();
2.237 -
2.238 - /// Starting from a flow, this method searches for an augmenting blocking
2.239 - /// flow according to Dinits' algorithm and augments the flow on if any.
2.240 - /// The blocking flow is computed in a physically constructed
2.241 - /// residual graph of type \c Mutablegraph.
2.242 - /// The return value show sif the augmentation was succesful.
2.243 - template<typename MutableGraph> bool augmentOnBlockingFlow();
2.244 -
2.245 - /// The same as \c augmentOnBlockingFlow<MutableGraph> but the
2.246 - /// residual graph is not constructed physically.
2.247 - /// The return value shows if the augmentation was succesful.
2.248 - bool augmentOnBlockingFlow2();
2.249 -
2.250 - /// Returns the maximum value of a flow.
2.251 -
2.252 - /// Returns the maximum value of a flow, by counting the
2.253 - /// over-flow of the target node \ref t.
2.254 - /// It can be called already after running \ref preflowPhase1.
2.255 - Num flowValue() const {
2.256 - Num a=0;
2.257 - FOR_EACH_INC_LOC(InEdgeIt, e, *g, t) a+=(*flow)[e];
2.258 - FOR_EACH_INC_LOC(OutEdgeIt, e, *g, t) a-=(*flow)[e];
2.259 - return a;
2.260 - //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan
2.261 - }
2.262 -
2.263 - ///Returns a minimum value cut after calling \ref preflowPhase1.
2.264 -
2.265 - ///After the first phase of the preflow algorithm the maximum flow
2.266 - ///value and a minimum value cut can already be computed. This
2.267 - ///method can be called after running \ref preflowPhase1 for
2.268 - ///obtaining a minimum value cut.
2.269 - /// \warning Gives proper result only right after calling \ref
2.270 - /// preflowPhase1.
2.271 - /// \todo We have to make some status variable which shows the
2.272 - /// actual state
2.273 - /// of the class. This enables us to determine which methods are valid
2.274 - /// for MinCut computation
2.275 - template<typename _CutMap>
2.276 - void actMinCut(_CutMap& M) const {
2.277 - NodeIt v;
2.278 - switch (status) {
2.279 - case AFTER_PRE_FLOW_PHASE_1:
2.280 - for(g->first(v); g->valid(v); g->next(v)) {
2.281 - if (level[v] < n) {
2.282 - M.set(v, false);
2.283 - } else {
2.284 - M.set(v, true);
2.285 - }
2.286 - }
2.287 - break;
2.288 - case AFTER_PRE_FLOW_PHASE_2:
2.289 - case AFTER_NOTHING:
2.290 - minMinCut(M);
2.291 - break;
2.292 - case AFTER_AUGMENTING:
2.293 - for(g->first(v); g->valid(v); g->next(v)) {
2.294 - if (level[v]) {
2.295 - M.set(v, true);
2.296 - } else {
2.297 - M.set(v, false);
2.298 - }
2.299 - }
2.300 - break;
2.301 - case AFTER_FAST_AUGMENTING:
2.302 - for(g->first(v); g->valid(v); g->next(v)) {
2.303 - if (level[v]==number_of_augmentations) {
2.304 - M.set(v, true);
2.305 - } else {
2.306 - M.set(v, false);
2.307 - }
2.308 - }
2.309 - break;
2.310 - }
2.311 - }
2.312 -
2.313 - ///Returns the inclusionwise minimum of the minimum value cuts.
2.314 -
2.315 - ///Sets \c M to the characteristic vector of the minimum value cut
2.316 - ///which is inclusionwise minimum. It is computed by processing
2.317 - ///a bfs from the source node \c s in the residual graph.
2.318 - ///\pre M should be a node map of bools initialized to false.
2.319 - ///\pre \c flow must be a maximum flow.
2.320 - template<typename _CutMap>
2.321 - void minMinCut(_CutMap& M) const {
2.322 - std::queue<Node> queue;
2.323 -
2.324 - M.set(s,true);
2.325 - queue.push(s);
2.326 -
2.327 - while (!queue.empty()) {
2.328 - Node w=queue.front();
2.329 - queue.pop();
2.330 -
2.331 - OutEdgeIt e;
2.332 - for(g->first(e,w) ; g->valid(e); g->next(e)) {
2.333 - Node v=g->head(e);
2.334 - if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
2.335 - queue.push(v);
2.336 - M.set(v, true);
2.337 - }
2.338 - }
2.339 -
2.340 - InEdgeIt f;
2.341 - for(g->first(f,w) ; g->valid(f); g->next(f)) {
2.342 - Node v=g->tail(f);
2.343 - if (!M[v] && (*flow)[f] > 0 ) {
2.344 - queue.push(v);
2.345 - M.set(v, true);
2.346 - }
2.347 - }
2.348 - }
2.349 - }
2.350 -
2.351 - ///Returns the inclusionwise maximum of the minimum value cuts.
2.352 -
2.353 - ///Sets \c M to the characteristic vector of the minimum value cut
2.354 - ///which is inclusionwise maximum. It is computed by processing a
2.355 - ///backward bfs from the target node \c t in the residual graph.
2.356 - ///\pre M should be a node map of bools initialized to false.
2.357 - ///\pre \c flow must be a maximum flow.
2.358 - template<typename _CutMap>
2.359 - void maxMinCut(_CutMap& M) const {
2.360 -
2.361 - NodeIt v;
2.362 - for(g->first(v) ; g->valid(v); g->next(v)) {
2.363 - M.set(v, true);
2.364 - }
2.365 -
2.366 - std::queue<Node> queue;
2.367 -
2.368 - M.set(t,false);
2.369 - queue.push(t);
2.370 -
2.371 - while (!queue.empty()) {
2.372 - Node w=queue.front();
2.373 - queue.pop();
2.374 -
2.375 - InEdgeIt e;
2.376 - for(g->first(e,w) ; g->valid(e); g->next(e)) {
2.377 - Node v=g->tail(e);
2.378 - if (M[v] && (*flow)[e] < (*capacity)[e] ) {
2.379 - queue.push(v);
2.380 - M.set(v, false);
2.381 - }
2.382 - }
2.383 -
2.384 - OutEdgeIt f;
2.385 - for(g->first(f,w) ; g->valid(f); g->next(f)) {
2.386 - Node v=g->head(f);
2.387 - if (M[v] && (*flow)[f] > 0 ) {
2.388 - queue.push(v);
2.389 - M.set(v, false);
2.390 - }
2.391 - }
2.392 - }
2.393 - }
2.394 -
2.395 - ///Returns a minimum value cut.
2.396 -
2.397 - ///Sets \c M to the characteristic vector of a minimum value cut.
2.398 - ///\pre M should be a node map of bools initialized to false.
2.399 - ///\pre \c flow must be a maximum flow.
2.400 - template<typename CutMap>
2.401 - void minCut(CutMap& M) const { minMinCut(M); }
2.402 -
2.403 - ///Resets the source node to \c _s.
2.404 -
2.405 - ///Resets the source node to \c _s.
2.406 - ///
2.407 - void resetSource(Node _s) { s=_s; status=AFTER_NOTHING; }
2.408 -
2.409 - ///Resets the target node to \c _t.
2.410 -
2.411 - ///Resets the target node to \c _t.
2.412 - ///
2.413 - void resetTarget(Node _t) { t=_t; status=AFTER_NOTHING; }
2.414 -
2.415 - /// Resets the edge map of the capacities to _cap.
2.416 -
2.417 - /// Resets the edge map of the capacities to _cap.
2.418 - ///
2.419 - void resetCap(const CapMap& _cap)
2.420 - { capacity=&_cap; status=AFTER_NOTHING; }
2.421 -
2.422 - /// Resets the edge map of the flows to _flow.
2.423 -
2.424 - /// Resets the edge map of the flows to _flow.
2.425 - ///
2.426 - void resetFlow(FlowMap& _flow) { flow=&_flow; status=AFTER_NOTHING; }
2.427 -
2.428 -
2.429 - private:
2.430 -
2.431 - int push(Node w, NNMap& next, VecFirst& first) {
2.432 -
2.433 - int lev=level[w];
2.434 - Num exc=excess[w];
2.435 - int newlevel=n; //bound on the next level of w
2.436 -
2.437 - OutEdgeIt e;
2.438 - for(g->first(e,w); g->valid(e); g->next(e)) {
2.439 -
2.440 - if ( (*flow)[e] >= (*capacity)[e] ) continue;
2.441 - Node v=g->head(e);
2.442 -
2.443 - if( lev > level[v] ) { //Push is allowed now
2.444 -
2.445 - if ( excess[v]<=0 && v!=t && v!=s ) {
2.446 - next.set(v,first[level[v]]);
2.447 - first[level[v]]=v;
2.448 - // int lev_v=level[v];
2.449 - //active[lev_v].push(v);
2.450 - }
2.451 -
2.452 - Num cap=(*capacity)[e];
2.453 - Num flo=(*flow)[e];
2.454 - Num remcap=cap-flo;
2.455 -
2.456 - if ( remcap >= exc ) { //A nonsaturating push.
2.457 -
2.458 - flow->set(e, flo+exc);
2.459 - excess.set(v, excess[v]+exc);
2.460 - exc=0;
2.461 - break;
2.462 -
2.463 - } else { //A saturating push.
2.464 - flow->set(e, cap);
2.465 - excess.set(v, excess[v]+remcap);
2.466 - exc-=remcap;
2.467 - }
2.468 - } else if ( newlevel > level[v] ) newlevel = level[v];
2.469 - } //for out edges wv
2.470 -
2.471 - if ( exc > 0 ) {
2.472 - InEdgeIt e;
2.473 - for(g->first(e,w); g->valid(e); g->next(e)) {
2.474 -
2.475 - if( (*flow)[e] <= 0 ) continue;
2.476 - Node v=g->tail(e);
2.477 -
2.478 - if( lev > level[v] ) { //Push is allowed now
2.479 -
2.480 - if ( excess[v]<=0 && v!=t && v!=s ) {
2.481 - next.set(v,first[level[v]]);
2.482 - first[level[v]]=v;
2.483 - //int lev_v=level[v];
2.484 - //active[lev_v].push(v);
2.485 - }
2.486 -
2.487 - Num flo=(*flow)[e];
2.488 -
2.489 - if ( flo >= exc ) { //A nonsaturating push.
2.490 -
2.491 - flow->set(e, flo-exc);
2.492 - excess.set(v, excess[v]+exc);
2.493 - exc=0;
2.494 - break;
2.495 - } else { //A saturating push.
2.496 -
2.497 - excess.set(v, excess[v]+flo);
2.498 - exc-=flo;
2.499 - flow->set(e,0);
2.500 - }
2.501 - } else if ( newlevel > level[v] ) newlevel = level[v];
2.502 - } //for in edges vw
2.503 -
2.504 - } // if w still has excess after the out edge for cycle
2.505 -
2.506 - excess.set(w, exc);
2.507 -
2.508 - return newlevel;
2.509 - }
2.510 -
2.511 -
2.512 - void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
2.513 - VecNode& level_list, NNMap& left, NNMap& right)
2.514 - {
2.515 - std::queue<Node> bfs_queue;
2.516 -
2.517 - switch (fe) {
2.518 - case NO_FLOW: //flow is already set to const zero in this case
2.519 - case ZERO_FLOW:
2.520 - {
2.521 - //Reverse_bfs from t, to find the starting level.
2.522 - level.set(t,0);
2.523 - bfs_queue.push(t);
2.524 -
2.525 - while (!bfs_queue.empty()) {
2.526 -
2.527 - Node v=bfs_queue.front();
2.528 - bfs_queue.pop();
2.529 - int l=level[v]+1;
2.530 -
2.531 - InEdgeIt e;
2.532 - for(g->first(e,v); g->valid(e); g->next(e)) {
2.533 - Node w=g->tail(e);
2.534 - if ( level[w] == n && w != s ) {
2.535 - bfs_queue.push(w);
2.536 - Node z=level_list[l];
2.537 - if ( g->valid(z) ) left.set(z,w);
2.538 - right.set(w,z);
2.539 - level_list[l]=w;
2.540 - level.set(w, l);
2.541 - }
2.542 - }
2.543 - }
2.544 -
2.545 - //the starting flow
2.546 - OutEdgeIt e;
2.547 - for(g->first(e,s); g->valid(e); g->next(e))
2.548 - {
2.549 - Num c=(*capacity)[e];
2.550 - if ( c <= 0 ) continue;
2.551 - Node w=g->head(e);
2.552 - if ( level[w] < n ) {
2.553 - if ( excess[w] <= 0 && w!=t )
2.554 - {
2.555 - next.set(w,first[level[w]]);
2.556 - first[level[w]]=w;
2.557 - //active[level[w]].push(w);
2.558 - }
2.559 - flow->set(e, c);
2.560 - excess.set(w, excess[w]+c);
2.561 - }
2.562 - }
2.563 - break;
2.564 - }
2.565 -
2.566 - case GEN_FLOW:
2.567 - case PRE_FLOW:
2.568 - {
2.569 - //Reverse_bfs from t in the residual graph,
2.570 - //to find the starting level.
2.571 - level.set(t,0);
2.572 - bfs_queue.push(t);
2.573 -
2.574 - while (!bfs_queue.empty()) {
2.575 -
2.576 - Node v=bfs_queue.front();
2.577 - bfs_queue.pop();
2.578 - int l=level[v]+1;
2.579 -
2.580 - InEdgeIt e;
2.581 - for(g->first(e,v); g->valid(e); g->next(e)) {
2.582 - if ( (*capacity)[e] <= (*flow)[e] ) continue;
2.583 - Node w=g->tail(e);
2.584 - if ( level[w] == n && w != s ) {
2.585 - bfs_queue.push(w);
2.586 - Node z=level_list[l];
2.587 - if ( g->valid(z) ) left.set(z,w);
2.588 - right.set(w,z);
2.589 - level_list[l]=w;
2.590 - level.set(w, l);
2.591 - }
2.592 - }
2.593 -
2.594 - OutEdgeIt f;
2.595 - for(g->first(f,v); g->valid(f); g->next(f)) {
2.596 - if ( 0 >= (*flow)[f] ) continue;
2.597 - Node w=g->head(f);
2.598 - if ( level[w] == n && w != s ) {
2.599 - bfs_queue.push(w);
2.600 - Node z=level_list[l];
2.601 - if ( g->valid(z) ) left.set(z,w);
2.602 - right.set(w,z);
2.603 - level_list[l]=w;
2.604 - level.set(w, l);
2.605 - }
2.606 - }
2.607 - }
2.608 -
2.609 -
2.610 - //the starting flow
2.611 - OutEdgeIt e;
2.612 - for(g->first(e,s); g->valid(e); g->next(e))
2.613 - {
2.614 - Num rem=(*capacity)[e]-(*flow)[e];
2.615 - if ( rem <= 0 ) continue;
2.616 - Node w=g->head(e);
2.617 - if ( level[w] < n ) {
2.618 - if ( excess[w] <= 0 && w!=t )
2.619 - {
2.620 - next.set(w,first[level[w]]);
2.621 - first[level[w]]=w;
2.622 - //active[level[w]].push(w);
2.623 - }
2.624 - flow->set(e, (*capacity)[e]);
2.625 - excess.set(w, excess[w]+rem);
2.626 - }
2.627 - }
2.628 -
2.629 - InEdgeIt f;
2.630 - for(g->first(f,s); g->valid(f); g->next(f))
2.631 - {
2.632 - if ( (*flow)[f] <= 0 ) continue;
2.633 - Node w=g->tail(f);
2.634 - if ( level[w] < n ) {
2.635 - if ( excess[w] <= 0 && w!=t )
2.636 - {
2.637 - next.set(w,first[level[w]]);
2.638 - first[level[w]]=w;
2.639 - //active[level[w]].push(w);
2.640 - }
2.641 - excess.set(w, excess[w]+(*flow)[f]);
2.642 - flow->set(f, 0);
2.643 - }
2.644 - }
2.645 - break;
2.646 - } //case PRE_FLOW
2.647 - }
2.648 - } //preflowPreproc
2.649 -
2.650 -
2.651 -
2.652 - void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
2.653 - VecNode& level_list, NNMap& left,
2.654 - NNMap& right, int& b, int& k, bool what_heur )
2.655 - {
2.656 -
2.657 - Num lev=level[w];
2.658 -
2.659 - Node right_n=right[w];
2.660 - Node left_n=left[w];
2.661 -
2.662 - //unlacing starts
2.663 - if ( g->valid(right_n) ) {
2.664 - if ( g->valid(left_n) ) {
2.665 - right.set(left_n, right_n);
2.666 - left.set(right_n, left_n);
2.667 - } else {
2.668 - level_list[lev]=right_n;
2.669 - left.set(right_n, INVALID);
2.670 - }
2.671 - } else {
2.672 - if ( g->valid(left_n) ) {
2.673 - right.set(left_n, INVALID);
2.674 - } else {
2.675 - level_list[lev]=INVALID;
2.676 - }
2.677 - }
2.678 - //unlacing ends
2.679 -
2.680 - if ( !g->valid(level_list[lev]) ) {
2.681 -
2.682 - //gapping starts
2.683 - for (int i=lev; i!=k ; ) {
2.684 - Node v=level_list[++i];
2.685 - while ( g->valid(v) ) {
2.686 - level.set(v,n);
2.687 - v=right[v];
2.688 - }
2.689 - level_list[i]=INVALID;
2.690 - if ( !what_heur ) first[i]=INVALID;
2.691 - /*{
2.692 - while ( !active[i].empty() ) {
2.693 - active[i].pop(); //FIXME: ezt szebben kene
2.694 - }
2.695 - }*/
2.696 - }
2.697 -
2.698 - level.set(w,n);
2.699 - b=lev-1;
2.700 - k=b;
2.701 - //gapping ends
2.702 -
2.703 - } else {
2.704 -
2.705 - if ( newlevel == n ) level.set(w,n);
2.706 - else {
2.707 - level.set(w,++newlevel);
2.708 - next.set(w,first[newlevel]);
2.709 - first[newlevel]=w;
2.710 - // active[newlevel].push(w);
2.711 - if ( what_heur ) b=newlevel;
2.712 - if ( k < newlevel ) ++k; //now k=newlevel
2.713 - Node z=level_list[newlevel];
2.714 - if ( g->valid(z) ) left.set(z,w);
2.715 - right.set(w,z);
2.716 - left.set(w,INVALID);
2.717 - level_list[newlevel]=w;
2.718 - }
2.719 - }
2.720 -
2.721 - } //relabel
2.722 -
2.723 -
2.724 - template<typename MapGraphWrapper>
2.725 - class DistanceMap {
2.726 - protected:
2.727 - const MapGraphWrapper* g;
2.728 - typename MapGraphWrapper::template NodeMap<int> dist;
2.729 - public:
2.730 - DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { }
2.731 - void set(const typename MapGraphWrapper::Node& n, int a) {
2.732 - dist.set(n, a);
2.733 - }
2.734 - int operator[](const typename MapGraphWrapper::Node& n) const {
2.735 - return dist[n];
2.736 - }
2.737 - // int get(const typename MapGraphWrapper::Node& n) const {
2.738 - // return dist[n]; }
2.739 - // bool get(const typename MapGraphWrapper::Edge& e) const {
2.740 - // return (dist.get(g->tail(e))<dist.get(g->head(e))); }
2.741 - bool operator[](const typename MapGraphWrapper::Edge& e) const {
2.742 - return (dist[g->tail(e)]<dist[g->head(e)]);
2.743 - }
2.744 - };
2.745 -
2.746 - };
2.747 -
2.748 -
2.749 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.750 - void MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::preflowPhase1(FlowEnum fe)
2.751 - {
2.752 -
2.753 - int heur0=(int)(H0*n); //time while running 'bound decrease'
2.754 - int heur1=(int)(H1*n); //time while running 'highest label'
2.755 - int heur=heur1; //starting time interval (#of relabels)
2.756 - int numrelabel=0;
2.757 -
2.758 - bool what_heur=1;
2.759 - //It is 0 in case 'bound decrease' and 1 in case 'highest label'
2.760 -
2.761 - bool end=false;
2.762 - //Needed for 'bound decrease', true means no active nodes are above bound
2.763 - //b.
2.764 -
2.765 - int k=n-2; //bound on the highest level under n containing a node
2.766 - int b=k; //bound on the highest level under n of an active node
2.767 -
2.768 - VecFirst first(n, INVALID);
2.769 - NNMap next(*g, INVALID); //maybe INVALID is not needed
2.770 - // VecStack active(n);
2.771 -
2.772 - NNMap left(*g, INVALID);
2.773 - NNMap right(*g, INVALID);
2.774 - VecNode level_list(n,INVALID);
2.775 - //List of the nodes in level i<n, set to n.
2.776 -
2.777 - NodeIt v;
2.778 - for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
2.779 - //setting each node to level n
2.780 -
2.781 - if ( fe == NO_FLOW ) {
2.782 - EdgeIt e;
2.783 - for(g->first(e); g->valid(e); g->next(e)) flow->set(e,0);
2.784 - }
2.785 -
2.786 - switch (fe) { //computing the excess
2.787 - case PRE_FLOW:
2.788 - {
2.789 - NodeIt v;
2.790 - for(g->first(v); g->valid(v); g->next(v)) {
2.791 - Num exc=0;
2.792 -
2.793 - InEdgeIt e;
2.794 - for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e];
2.795 - OutEdgeIt f;
2.796 - for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f];
2.797 -
2.798 - excess.set(v,exc);
2.799 -
2.800 - //putting the active nodes into the stack
2.801 - int lev=level[v];
2.802 - if ( exc > 0 && lev < n && v != t )
2.803 - {
2.804 - next.set(v,first[lev]);
2.805 - first[lev]=v;
2.806 - }
2.807 - // active[lev].push(v);
2.808 - }
2.809 - break;
2.810 - }
2.811 - case GEN_FLOW:
2.812 - {
2.813 - NodeIt v;
2.814 - for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
2.815 -
2.816 - Num exc=0;
2.817 - InEdgeIt e;
2.818 - for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
2.819 - OutEdgeIt f;
2.820 - for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
2.821 - excess.set(t,exc);
2.822 - break;
2.823 - }
2.824 - case ZERO_FLOW:
2.825 - case NO_FLOW:
2.826 - {
2.827 - NodeIt v;
2.828 - for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
2.829 - break;
2.830 - }
2.831 - }
2.832 -
2.833 - preflowPreproc(fe, next, first,/*active*/ level_list, left, right);
2.834 - //End of preprocessing
2.835 -
2.836 -
2.837 - //Push/relabel on the highest level active nodes.
2.838 - while ( true ) {
2.839 - if ( b == 0 ) {
2.840 - if ( !what_heur && !end && k > 0 ) {
2.841 - b=k;
2.842 - end=true;
2.843 - } else break;
2.844 - }
2.845 -
2.846 - if ( !g->valid(first[b])/*active[b].empty()*/ ) --b;
2.847 - else {
2.848 - end=false;
2.849 - Node w=first[b];
2.850 - first[b]=next[w];
2.851 - /* Node w=active[b].top();
2.852 - active[b].pop();*/
2.853 - int newlevel=push(w,/*active*/next, first);
2.854 - if ( excess[w] > 0 ) relabel(w, newlevel, /*active*/next, first, level_list,
2.855 - left, right, b, k, what_heur);
2.856 -
2.857 - ++numrelabel;
2.858 - if ( numrelabel >= heur ) {
2.859 - numrelabel=0;
2.860 - if ( what_heur ) {
2.861 - what_heur=0;
2.862 - heur=heur0;
2.863 - end=false;
2.864 - } else {
2.865 - what_heur=1;
2.866 - heur=heur1;
2.867 - b=k;
2.868 - }
2.869 - }
2.870 - }
2.871 - }
2.872 -
2.873 - status=AFTER_PRE_FLOW_PHASE_1;
2.874 - }
2.875 -
2.876 -
2.877 -
2.878 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.879 - void MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::preflowPhase2()
2.880 - {
2.881 -
2.882 - int k=n-2; //bound on the highest level under n containing a node
2.883 - int b=k; //bound on the highest level under n of an active node
2.884 -
2.885 -
2.886 - VecFirst first(n, INVALID);
2.887 - NNMap next(*g, INVALID); //maybe INVALID is not needed
2.888 - // VecStack active(n);
2.889 - level.set(s,0);
2.890 - std::queue<Node> bfs_queue;
2.891 - bfs_queue.push(s);
2.892 -
2.893 - while (!bfs_queue.empty()) {
2.894 -
2.895 - Node v=bfs_queue.front();
2.896 - bfs_queue.pop();
2.897 - int l=level[v]+1;
2.898 -
2.899 - InEdgeIt e;
2.900 - for(g->first(e,v); g->valid(e); g->next(e)) {
2.901 - if ( (*capacity)[e] <= (*flow)[e] ) continue;
2.902 - Node u=g->tail(e);
2.903 - if ( level[u] >= n ) {
2.904 - bfs_queue.push(u);
2.905 - level.set(u, l);
2.906 - if ( excess[u] > 0 ) {
2.907 - next.set(u,first[l]);
2.908 - first[l]=u;
2.909 - //active[l].push(u);
2.910 - }
2.911 - }
2.912 - }
2.913 -
2.914 - OutEdgeIt f;
2.915 - for(g->first(f,v); g->valid(f); g->next(f)) {
2.916 - if ( 0 >= (*flow)[f] ) continue;
2.917 - Node u=g->head(f);
2.918 - if ( level[u] >= n ) {
2.919 - bfs_queue.push(u);
2.920 - level.set(u, l);
2.921 - if ( excess[u] > 0 ) {
2.922 - next.set(u,first[l]);
2.923 - first[l]=u;
2.924 - //active[l].push(u);
2.925 - }
2.926 - }
2.927 - }
2.928 - }
2.929 - b=n-2;
2.930 -
2.931 - while ( true ) {
2.932 -
2.933 - if ( b == 0 ) break;
2.934 -
2.935 - if ( !g->valid(first[b])/*active[b].empty()*/ ) --b;
2.936 - else {
2.937 -
2.938 - Node w=first[b];
2.939 - first[b]=next[w];
2.940 - /* Node w=active[b].top();
2.941 - active[b].pop();*/
2.942 - int newlevel=push(w,next, first/*active*/);
2.943 -
2.944 - //relabel
2.945 - if ( excess[w] > 0 ) {
2.946 - level.set(w,++newlevel);
2.947 - next.set(w,first[newlevel]);
2.948 - first[newlevel]=w;
2.949 - //active[newlevel].push(w);
2.950 - b=newlevel;
2.951 - }
2.952 - } // if stack[b] is nonempty
2.953 - } // while(true)
2.954 -
2.955 - status=AFTER_PRE_FLOW_PHASE_2;
2.956 - }
2.957 -
2.958 -
2.959 -
2.960 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.961 - bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
2.962 - {
2.963 - ResGW res_graph(*g, *capacity, *flow);
2.964 - bool _augment=false;
2.965 -
2.966 - //ReachedMap level(res_graph);
2.967 - FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
2.968 - BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
2.969 - bfs.pushAndSetReached(s);
2.970 -
2.971 - typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
2.972 - pred.set(s, INVALID);
2.973 -
2.974 - typename ResGW::template NodeMap<Num> free(res_graph);
2.975 -
2.976 - //searching for augmenting path
2.977 - while ( !bfs.finished() ) {
2.978 - ResGWOutEdgeIt e=bfs;
2.979 - if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
2.980 - Node v=res_graph.tail(e);
2.981 - Node w=res_graph.head(e);
2.982 - pred.set(w, e);
2.983 - if (res_graph.valid(pred[v])) {
2.984 - free.set(w, std::min(free[v], res_graph.resCap(e)));
2.985 - } else {
2.986 - free.set(w, res_graph.resCap(e));
2.987 - }
2.988 - if (res_graph.head(e)==t) { _augment=true; break; }
2.989 - }
2.990 -
2.991 - ++bfs;
2.992 - } //end of searching augmenting path
2.993 -
2.994 - if (_augment) {
2.995 - Node n=t;
2.996 - Num augment_value=free[t];
2.997 - while (res_graph.valid(pred[n])) {
2.998 - ResGWEdge e=pred[n];
2.999 - res_graph.augment(e, augment_value);
2.1000 - n=res_graph.tail(e);
2.1001 - }
2.1002 - }
2.1003 -
2.1004 - status=AFTER_AUGMENTING;
2.1005 - return _augment;
2.1006 - }
2.1007 -
2.1008 -
2.1009 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.1010 - bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath2()
2.1011 - {
2.1012 - ResGW res_graph(*g, *capacity, *flow);
2.1013 - bool _augment=false;
2.1014 -
2.1015 - if (status!=AFTER_FAST_AUGMENTING) {
2.1016 - FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
2.1017 - number_of_augmentations=1;
2.1018 - } else {
2.1019 - ++number_of_augmentations;
2.1020 - }
2.1021 - TrickyReachedMap<ReachedMap>
2.1022 - tricky_reached_map(level, number_of_augmentations);
2.1023 - //ReachedMap level(res_graph);
2.1024 -// FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
2.1025 - BfsIterator<ResGW, TrickyReachedMap<ReachedMap> >
2.1026 - bfs(res_graph, tricky_reached_map);
2.1027 - bfs.pushAndSetReached(s);
2.1028 -
2.1029 - typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
2.1030 - pred.set(s, INVALID);
2.1031 -
2.1032 - typename ResGW::template NodeMap<Num> free(res_graph);
2.1033 -
2.1034 - //searching for augmenting path
2.1035 - while ( !bfs.finished() ) {
2.1036 - ResGWOutEdgeIt e=bfs;
2.1037 - if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
2.1038 - Node v=res_graph.tail(e);
2.1039 - Node w=res_graph.head(e);
2.1040 - pred.set(w, e);
2.1041 - if (res_graph.valid(pred[v])) {
2.1042 - free.set(w, std::min(free[v], res_graph.resCap(e)));
2.1043 - } else {
2.1044 - free.set(w, res_graph.resCap(e));
2.1045 - }
2.1046 - if (res_graph.head(e)==t) { _augment=true; break; }
2.1047 - }
2.1048 -
2.1049 - ++bfs;
2.1050 - } //end of searching augmenting path
2.1051 -
2.1052 - if (_augment) {
2.1053 - Node n=t;
2.1054 - Num augment_value=free[t];
2.1055 - while (res_graph.valid(pred[n])) {
2.1056 - ResGWEdge e=pred[n];
2.1057 - res_graph.augment(e, augment_value);
2.1058 - n=res_graph.tail(e);
2.1059 - }
2.1060 - }
2.1061 -
2.1062 - status=AFTER_FAST_AUGMENTING;
2.1063 - return _augment;
2.1064 - }
2.1065 -
2.1066 -
2.1067 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.1068 - template<typename MutableGraph>
2.1069 - bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
2.1070 - {
2.1071 - typedef MutableGraph MG;
2.1072 - bool _augment=false;
2.1073 -
2.1074 - ResGW res_graph(*g, *capacity, *flow);
2.1075 -
2.1076 - //bfs for distances on the residual graph
2.1077 - //ReachedMap level(res_graph);
2.1078 - FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
2.1079 - BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
2.1080 - bfs.pushAndSetReached(s);
2.1081 - typename ResGW::template NodeMap<int>
2.1082 - dist(res_graph); //filled up with 0's
2.1083 -
2.1084 - //F will contain the physical copy of the residual graph
2.1085 - //with the set of edges which are on shortest paths
2.1086 - MG F;
2.1087 - typename ResGW::template NodeMap<typename MG::Node>
2.1088 - res_graph_to_F(res_graph);
2.1089 - {
2.1090 - typename ResGW::NodeIt n;
2.1091 - for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) {
2.1092 - res_graph_to_F.set(n, F.addNode());
2.1093 - }
2.1094 - }
2.1095 -
2.1096 - typename MG::Node sF=res_graph_to_F[s];
2.1097 - typename MG::Node tF=res_graph_to_F[t];
2.1098 - typename MG::template EdgeMap<ResGWEdge> original_edge(F);
2.1099 - typename MG::template EdgeMap<Num> residual_capacity(F);
2.1100 -
2.1101 - while ( !bfs.finished() ) {
2.1102 - ResGWOutEdgeIt e=bfs;
2.1103 - if (res_graph.valid(e)) {
2.1104 - if (bfs.isBNodeNewlyReached()) {
2.1105 - dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
2.1106 - typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
2.1107 - res_graph_to_F[res_graph.head(e)]);
2.1108 - original_edge.update();
2.1109 - original_edge.set(f, e);
2.1110 - residual_capacity.update();
2.1111 - residual_capacity.set(f, res_graph.resCap(e));
2.1112 - } else {
2.1113 - if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
2.1114 - typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
2.1115 - res_graph_to_F[res_graph.head(e)]);
2.1116 - original_edge.update();
2.1117 - original_edge.set(f, e);
2.1118 - residual_capacity.update();
2.1119 - residual_capacity.set(f, res_graph.resCap(e));
2.1120 - }
2.1121 - }
2.1122 - }
2.1123 - ++bfs;
2.1124 - } //computing distances from s in the residual graph
2.1125 -
2.1126 - bool __augment=true;
2.1127 -
2.1128 - while (__augment) {
2.1129 - __augment=false;
2.1130 - //computing blocking flow with dfs
2.1131 - DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
2.1132 - typename MG::template NodeMap<typename MG::Edge> pred(F);
2.1133 - pred.set(sF, INVALID);
2.1134 - //invalid iterators for sources
2.1135 -
2.1136 - typename MG::template NodeMap<Num> free(F);
2.1137 -
2.1138 - dfs.pushAndSetReached(sF);
2.1139 - while (!dfs.finished()) {
2.1140 - ++dfs;
2.1141 - if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) {
2.1142 - if (dfs.isBNodeNewlyReached()) {
2.1143 - typename MG::Node v=F.aNode(dfs);
2.1144 - typename MG::Node w=F.bNode(dfs);
2.1145 - pred.set(w, dfs);
2.1146 - if (F.valid(pred[v])) {
2.1147 - free.set(w, std::min(free[v], residual_capacity[dfs]));
2.1148 - } else {
2.1149 - free.set(w, residual_capacity[dfs]);
2.1150 - }
2.1151 - if (w==tF) {
2.1152 - __augment=true;
2.1153 - _augment=true;
2.1154 - break;
2.1155 - }
2.1156 -
2.1157 - } else {
2.1158 - F.erase(/*typename MG::OutEdgeIt*/(dfs));
2.1159 - }
2.1160 - }
2.1161 - }
2.1162 -
2.1163 - if (__augment) {
2.1164 - typename MG::Node n=tF;
2.1165 - Num augment_value=free[tF];
2.1166 - while (F.valid(pred[n])) {
2.1167 - typename MG::Edge e=pred[n];
2.1168 - res_graph.augment(original_edge[e], augment_value);
2.1169 - n=F.tail(e);
2.1170 - if (residual_capacity[e]==augment_value)
2.1171 - F.erase(e);
2.1172 - else
2.1173 - residual_capacity.set(e, residual_capacity[e]-augment_value);
2.1174 - }
2.1175 - }
2.1176 -
2.1177 - }
2.1178 -
2.1179 - status=AFTER_AUGMENTING;
2.1180 - return _augment;
2.1181 - }
2.1182 -
2.1183 -
2.1184 -
2.1185 -
2.1186 - template <typename Graph, typename Num, typename CapMap, typename FlowMap>
2.1187 - bool MaxFlowNoStack<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
2.1188 - {
2.1189 - bool _augment=false;
2.1190 -
2.1191 - ResGW res_graph(*g, *capacity, *flow);
2.1192 -
2.1193 - //ReachedMap level(res_graph);
2.1194 - FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
2.1195 - BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
2.1196 -
2.1197 - bfs.pushAndSetReached(s);
2.1198 - DistanceMap<ResGW> dist(res_graph);
2.1199 - while ( !bfs.finished() ) {
2.1200 - ResGWOutEdgeIt e=bfs;
2.1201 - if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) {
2.1202 - dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
2.1203 - }
2.1204 - ++bfs;
2.1205 - } //computing distances from s in the residual graph
2.1206 -
2.1207 - //Subgraph containing the edges on some shortest paths
2.1208 - ConstMap<typename ResGW::Node, bool> true_map(true);
2.1209 - typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
2.1210 - DistanceMap<ResGW> > FilterResGW;
2.1211 - FilterResGW filter_res_graph(res_graph, true_map, dist);
2.1212 -
2.1213 - //Subgraph, which is able to delete edges which are already
2.1214 - //met by the dfs
2.1215 - typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt>
2.1216 - first_out_edges(filter_res_graph);
2.1217 - typename FilterResGW::NodeIt v;
2.1218 - for(filter_res_graph.first(v); filter_res_graph.valid(v);
2.1219 - filter_res_graph.next(v))
2.1220 - {
2.1221 - typename FilterResGW::OutEdgeIt e;
2.1222 - filter_res_graph.first(e, v);
2.1223 - first_out_edges.set(v, e);
2.1224 - }
2.1225 - typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
2.1226 - template NodeMap<typename FilterResGW::OutEdgeIt> > ErasingResGW;
2.1227 - ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
2.1228 -
2.1229 - bool __augment=true;
2.1230 -
2.1231 - while (__augment) {
2.1232 -
2.1233 - __augment=false;
2.1234 - //computing blocking flow with dfs
2.1235 - DfsIterator< ErasingResGW,
2.1236 - typename ErasingResGW::template NodeMap<bool> >
2.1237 - dfs(erasing_res_graph);
2.1238 - typename ErasingResGW::
2.1239 - template NodeMap<typename ErasingResGW::OutEdgeIt>
2.1240 - pred(erasing_res_graph);
2.1241 - pred.set(s, INVALID);
2.1242 - //invalid iterators for sources
2.1243 -
2.1244 - typename ErasingResGW::template NodeMap<Num>
2.1245 - free1(erasing_res_graph);
2.1246 -
2.1247 - dfs.pushAndSetReached
2.1248 - ///\bug hugo 0.2
2.1249 - (typename ErasingResGW::Node
2.1250 - (typename FilterResGW::Node
2.1251 - (typename ResGW::Node(s)
2.1252 - )
2.1253 - )
2.1254 - );
2.1255 - while (!dfs.finished()) {
2.1256 - ++dfs;
2.1257 - if (erasing_res_graph.valid(typename ErasingResGW::OutEdgeIt(dfs)))
2.1258 - {
2.1259 - if (dfs.isBNodeNewlyReached()) {
2.1260 -
2.1261 - typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs);
2.1262 - typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs);
2.1263 -
2.1264 - pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs));
2.1265 - if (erasing_res_graph.valid(pred[v])) {
2.1266 - free1.set
2.1267 - (w, std::min(free1[v], res_graph.resCap
2.1268 - (typename ErasingResGW::OutEdgeIt(dfs))));
2.1269 - } else {
2.1270 - free1.set
2.1271 - (w, res_graph.resCap
2.1272 - (typename ErasingResGW::OutEdgeIt(dfs)));
2.1273 - }
2.1274 -
2.1275 - if (w==t) {
2.1276 - __augment=true;
2.1277 - _augment=true;
2.1278 - break;
2.1279 - }
2.1280 - } else {
2.1281 - erasing_res_graph.erase(dfs);
2.1282 - }
2.1283 - }
2.1284 - }
2.1285 -
2.1286 - if (__augment) {
2.1287 - typename ErasingResGW::Node
2.1288 - n=typename FilterResGW::Node(typename ResGW::Node(t));
2.1289 - // typename ResGW::NodeMap<Num> a(res_graph);
2.1290 - // typename ResGW::Node b;
2.1291 - // Num j=a[b];
2.1292 - // typename FilterResGW::NodeMap<Num> a1(filter_res_graph);
2.1293 - // typename FilterResGW::Node b1;
2.1294 - // Num j1=a1[b1];
2.1295 - // typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph);
2.1296 - // typename ErasingResGW::Node b2;
2.1297 - // Num j2=a2[b2];
2.1298 - Num augment_value=free1[n];
2.1299 - while (erasing_res_graph.valid(pred[n])) {
2.1300 - typename ErasingResGW::OutEdgeIt e=pred[n];
2.1301 - res_graph.augment(e, augment_value);
2.1302 - n=erasing_res_graph.tail(e);
2.1303 - if (res_graph.resCap(e)==0)
2.1304 - erasing_res_graph.erase(e);
2.1305 - }
2.1306 - }
2.1307 -
2.1308 - } //while (__augment)
2.1309 -
2.1310 - status=AFTER_AUGMENTING;
2.1311 - return _augment;
2.1312 - }
2.1313 -
2.1314 -
2.1315 -} //namespace hugo
2.1316 -
2.1317 -#endif //HUGO_MAX_FLOW_H
2.1318 -
2.1319 -
2.1320 -
2.1321 -