1.1 --- a/src/hugo/max_flow.h Mon Sep 13 11:24:35 2004 +0000
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,903 +0,0 @@
1.4 -// -*- C++ -*-
1.5 -#ifndef HUGO_MAX_FLOW_H
1.6 -#define HUGO_MAX_FLOW_H
1.7 -
1.8 -#include <vector>
1.9 -#include <queue>
1.10 -
1.11 -//#include <hugo/graph_wrapper.h>
1.12 -#include <hugo/invalid.h>
1.13 -#include <hugo/maps.h>
1.14 -
1.15 -/// \file
1.16 -/// \ingroup flowalgs
1.17 -
1.18 -namespace hugo {
1.19 -
1.20 - /// \addtogroup flowalgs
1.21 - /// @{
1.22 -
1.23 - ///Maximum flow algorithms class.
1.24 -
1.25 - ///This class provides various algorithms for finding a flow of
1.26 - ///maximum value in a directed graph. The \e source node, the \e
1.27 - ///target node, the \e capacity of the edges and the \e starting \e
1.28 - ///flow value of the edges should be passed to the algorithm through the
1.29 - ///constructor. It is possible to change these quantities using the
1.30 - ///functions \ref setSource, \ref setTarget, \ref setCap and
1.31 - ///\ref setFlow. Before any subsequent runs of any algorithm of
1.32 - ///the class \ref setFlow should be called.
1.33 - ///
1.34 - ///After running an algorithm of the class, the actual flow value
1.35 - ///can be obtained by calling \ref flowValue(). The minimum
1.36 - ///value cut can be written into a \c node map of \c bools by
1.37 - ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
1.38 - ///the inclusionwise minimum and maximum of the minimum value
1.39 - ///cuts, resp.)
1.40 - ///
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 - ///
1.46 - ///\author Marton Makai, Jacint Szabo
1.47 - template <typename Graph, typename Num,
1.48 - typename CapMap=typename Graph::template EdgeMap<Num>,
1.49 - typename FlowMap=typename Graph::template EdgeMap<Num> >
1.50 - class MaxFlow {
1.51 - protected:
1.52 - typedef typename Graph::Node Node;
1.53 - typedef typename Graph::NodeIt NodeIt;
1.54 - typedef typename Graph::EdgeIt EdgeIt;
1.55 - typedef typename Graph::OutEdgeIt OutEdgeIt;
1.56 - typedef typename Graph::InEdgeIt InEdgeIt;
1.57 -
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 Graph::template NodeMap<int> ReachedMap;
1.73 -
1.74 -
1.75 - //level works as a bool map in augmenting path algorithms and is
1.76 - //used by bfs for storing reached information. In preflow, it
1.77 - //shows the levels of nodes.
1.78 - ReachedMap level;
1.79 -
1.80 - //excess is needed only in preflow
1.81 - typename Graph::template NodeMap<Num> excess;
1.82 -
1.83 - // constants used for heuristics
1.84 - static const int H0=20;
1.85 - static const int H1=1;
1.86 -
1.87 - public:
1.88 -
1.89 - ///Indicates the property of the starting flow.
1.90 -
1.91 - ///Indicates the property of the starting flow. The meanings are as follows:
1.92 - ///- \c ZERO_FLOW: constant zero flow
1.93 - ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
1.94 - ///the sum of the out-flows in every node except the \e source and
1.95 - ///the \e target.
1.96 - ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
1.97 - ///least the sum of the out-flows in every node except the \e source.
1.98 - ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be
1.99 - ///set to the constant zero flow in the beginning of the algorithm in this case.
1.100 - enum FlowEnum{
1.101 - ZERO_FLOW,
1.102 - GEN_FLOW,
1.103 - PRE_FLOW,
1.104 - NO_FLOW
1.105 - };
1.106 -
1.107 - enum StatusEnum {
1.108 - AFTER_NOTHING,
1.109 - AFTER_AUGMENTING,
1.110 - AFTER_FAST_AUGMENTING,
1.111 - AFTER_PRE_FLOW_PHASE_1,
1.112 - AFTER_PRE_FLOW_PHASE_2
1.113 - };
1.114 -
1.115 - /// Do not needle this flag only if necessary.
1.116 - StatusEnum status;
1.117 -
1.118 - // int number_of_augmentations;
1.119 -
1.120 -
1.121 - // template<typename IntMap>
1.122 - // class TrickyReachedMap {
1.123 - // protected:
1.124 - // IntMap* map;
1.125 - // int* number_of_augmentations;
1.126 - // public:
1.127 - // TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) :
1.128 - // map(&_map), number_of_augmentations(&_number_of_augmentations) { }
1.129 - // void set(const Node& n, bool b) {
1.130 - // if (b)
1.131 - // map->set(n, *number_of_augmentations);
1.132 - // else
1.133 - // map->set(n, *number_of_augmentations-1);
1.134 - // }
1.135 - // bool operator[](const Node& n) const {
1.136 - // return (*map)[n]==*number_of_augmentations;
1.137 - // }
1.138 - // };
1.139 -
1.140 - ///Constructor
1.141 -
1.142 - ///\todo Document, please.
1.143 - ///
1.144 - MaxFlow(const Graph& _G, Node _s, Node _t,
1.145 - const CapMap& _capacity, FlowMap& _flow) :
1.146 - g(&_G), s(_s), t(_t), capacity(&_capacity),
1.147 - flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0),
1.148 - status(AFTER_NOTHING) { }
1.149 -
1.150 - ///Runs a maximum flow algorithm.
1.151 -
1.152 - ///Runs a preflow algorithm, which is the fastest maximum flow
1.153 - ///algorithm up-to-date. The default for \c fe is ZERO_FLOW.
1.154 - ///\pre The starting flow must be
1.155 - /// - a constant zero flow if \c fe is \c ZERO_FLOW,
1.156 - /// - an arbitary flow if \c fe is \c GEN_FLOW,
1.157 - /// - an arbitary preflow if \c fe is \c PRE_FLOW,
1.158 - /// - any map if \c fe is NO_FLOW.
1.159 - void run(FlowEnum fe=ZERO_FLOW) {
1.160 - preflow(fe);
1.161 - }
1.162 -
1.163 -
1.164 - ///Runs a preflow algorithm.
1.165 -
1.166 - ///Runs a preflow algorithm. The preflow algorithms provide the
1.167 - ///fastest way to compute a maximum flow in a directed graph.
1.168 - ///\pre The starting flow must be
1.169 - /// - a constant zero flow if \c fe is \c ZERO_FLOW,
1.170 - /// - an arbitary flow if \c fe is \c GEN_FLOW,
1.171 - /// - an arbitary preflow if \c fe is \c PRE_FLOW,
1.172 - /// - any map if \c fe is NO_FLOW.
1.173 - ///
1.174 - ///\todo NO_FLOW should be the default flow.
1.175 - void preflow(FlowEnum fe) {
1.176 - preflowPhase1(fe);
1.177 - preflowPhase2();
1.178 - }
1.179 - // Heuristics:
1.180 - // 2 phase
1.181 - // gap
1.182 - // list 'level_list' on the nodes on level i implemented by hand
1.183 - // stack 'active' on the active nodes on level i
1.184 - // runs heuristic 'highest label' for H1*n relabels
1.185 - // runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
1.186 - // Parameters H0 and H1 are initialized to 20 and 1.
1.187 -
1.188 - ///Runs the first phase of the preflow algorithm.
1.189 -
1.190 - ///The preflow algorithm consists of two phases, this method runs the
1.191 - ///first phase. After the first phase the maximum flow value and a
1.192 - ///minimum value cut can already be computed, though a maximum flow
1.193 - ///is not yet obtained. So after calling this method \ref flowValue
1.194 - ///and \ref actMinCut gives proper results.
1.195 - ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
1.196 - ///give minimum value cuts unless calling \ref preflowPhase2.
1.197 - ///\pre The starting flow must be
1.198 - /// - a constant zero flow if \c fe is \c ZERO_FLOW,
1.199 - /// - an arbitary flow if \c fe is \c GEN_FLOW,
1.200 - /// - an arbitary preflow if \c fe is \c PRE_FLOW,
1.201 - /// - any map if \c fe is NO_FLOW.
1.202 - void preflowPhase1(FlowEnum fe)
1.203 - {
1.204 -
1.205 - int heur0=(int)(H0*n); //time while running 'bound decrease'
1.206 - int heur1=(int)(H1*n); //time while running 'highest label'
1.207 - int heur=heur1; //starting time interval (#of relabels)
1.208 - int numrelabel=0;
1.209 -
1.210 - bool what_heur=1;
1.211 - //It is 0 in case 'bound decrease' and 1 in case 'highest label'
1.212 -
1.213 - bool end=false;
1.214 - //Needed for 'bound decrease', true means no active nodes are above bound
1.215 - //b.
1.216 -
1.217 - int k=n-2; //bound on the highest level under n containing a node
1.218 - int b=k; //bound on the highest level under n of an active node
1.219 -
1.220 - VecFirst first(n, INVALID);
1.221 - NNMap next(*g, INVALID); //maybe INVALID is not needed
1.222 -
1.223 - NNMap left(*g, INVALID);
1.224 - NNMap right(*g, INVALID);
1.225 - VecNode level_list(n,INVALID);
1.226 - //List of the nodes in level i<n, set to n.
1.227 -
1.228 - preflowPreproc(fe, next, first, level_list, left, right);
1.229 - //End of preprocessing
1.230 -
1.231 - //Push/relabel on the highest level active nodes.
1.232 - while ( true ) {
1.233 - if ( b == 0 ) {
1.234 - if ( !what_heur && !end && k > 0 ) {
1.235 - b=k;
1.236 - end=true;
1.237 - } else break;
1.238 - }
1.239 -
1.240 - if ( first[b]==INVALID ) --b;
1.241 - else {
1.242 - end=false;
1.243 - Node w=first[b];
1.244 - first[b]=next[w];
1.245 - int newlevel=push(w, next, first);
1.246 - if ( excess[w] > 0 ) relabel(w, newlevel, next, first, level_list,
1.247 - left, right, b, k, what_heur);
1.248 -
1.249 - ++numrelabel;
1.250 - if ( numrelabel >= heur ) {
1.251 - numrelabel=0;
1.252 - if ( what_heur ) {
1.253 - what_heur=0;
1.254 - heur=heur0;
1.255 - end=false;
1.256 - } else {
1.257 - what_heur=1;
1.258 - heur=heur1;
1.259 - b=k;
1.260 - }
1.261 - }
1.262 - }
1.263 - }
1.264 -
1.265 - status=AFTER_PRE_FLOW_PHASE_1;
1.266 - }
1.267 -
1.268 -
1.269 - ///Runs the second phase of the preflow algorithm.
1.270 -
1.271 - ///The preflow algorithm consists of two phases, this method runs
1.272 - ///the second phase. After calling \ref preflowPhase1 and then
1.273 - ///\ref preflowPhase2 the methods \ref flowValue, \ref minCut,
1.274 - ///\ref minMinCut and \ref maxMinCut give proper results.
1.275 - ///\pre \ref preflowPhase1 must be called before.
1.276 - void preflowPhase2()
1.277 - {
1.278 -
1.279 - int k=n-2; //bound on the highest level under n containing a node
1.280 - int b=k; //bound on the highest level under n of an active node
1.281 -
1.282 -
1.283 - VecFirst first(n, INVALID);
1.284 - NNMap next(*g, INVALID); //maybe INVALID is not needed
1.285 - level.set(s,0);
1.286 - std::queue<Node> bfs_queue;
1.287 - bfs_queue.push(s);
1.288 -
1.289 - while (!bfs_queue.empty()) {
1.290 -
1.291 - Node v=bfs_queue.front();
1.292 - bfs_queue.pop();
1.293 - int l=level[v]+1;
1.294 -
1.295 - for(InEdgeIt e(*g,v); e!=INVALID; ++e) {
1.296 - if ( (*capacity)[e] <= (*flow)[e] ) continue;
1.297 - Node u=g->tail(e);
1.298 - if ( level[u] >= n ) {
1.299 - bfs_queue.push(u);
1.300 - level.set(u, l);
1.301 - if ( excess[u] > 0 ) {
1.302 - next.set(u,first[l]);
1.303 - first[l]=u;
1.304 - }
1.305 - }
1.306 - }
1.307 -
1.308 - for(OutEdgeIt e(*g,v); e!=INVALID; ++e) {
1.309 - if ( 0 >= (*flow)[e] ) continue;
1.310 - Node u=g->head(e);
1.311 - if ( level[u] >= n ) {
1.312 - bfs_queue.push(u);
1.313 - level.set(u, l);
1.314 - if ( excess[u] > 0 ) {
1.315 - next.set(u,first[l]);
1.316 - first[l]=u;
1.317 - }
1.318 - }
1.319 - }
1.320 - }
1.321 - b=n-2;
1.322 -
1.323 - while ( true ) {
1.324 -
1.325 - if ( b == 0 ) break;
1.326 -
1.327 - if ( first[b]==INVALID ) --b;
1.328 - else {
1.329 -
1.330 - Node w=first[b];
1.331 - first[b]=next[w];
1.332 - int newlevel=push(w,next, first/*active*/);
1.333 -
1.334 - //relabel
1.335 - if ( excess[w] > 0 ) {
1.336 - level.set(w,++newlevel);
1.337 - next.set(w,first[newlevel]);
1.338 - first[newlevel]=w;
1.339 - b=newlevel;
1.340 - }
1.341 - }
1.342 - } // while(true)
1.343 -
1.344 - status=AFTER_PRE_FLOW_PHASE_2;
1.345 - }
1.346 -
1.347 -
1.348 - /// Returns the value of the maximum flow.
1.349 -
1.350 - /// Returns the excess of the target node \ref t.
1.351 - /// After running \ref preflowPhase1, this is the value of
1.352 - /// the maximum flow.
1.353 - /// It can be called already after running \ref preflowPhase1.
1.354 - Num flowValue() const {
1.355 - // Num a=0;
1.356 - // for(InEdgeIt e(*g,t);g->valid(e);g->next(e)) a+=(*flow)[e];
1.357 - // for(OutEdgeIt e(*g,t);g->valid(e);g->next(e)) a-=(*flow)[e];
1.358 - // return a;
1.359 - return excess[t];
1.360 - //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan
1.361 - }
1.362 -
1.363 -
1.364 - ///Returns a minimum value cut after calling \ref preflowPhase1.
1.365 -
1.366 - ///After the first phase of the preflow algorithm the maximum flow
1.367 - ///value and a minimum value cut can already be computed. This
1.368 - ///method can be called after running \ref preflowPhase1 for
1.369 - ///obtaining a minimum value cut.
1.370 - /// \warning Gives proper result only right after calling \ref
1.371 - /// preflowPhase1.
1.372 - /// \todo We have to make some status variable which shows the
1.373 - /// actual state
1.374 - /// of the class. This enables us to determine which methods are valid
1.375 - /// for MinCut computation
1.376 - template<typename _CutMap>
1.377 - void actMinCut(_CutMap& M) const {
1.378 - switch (status) {
1.379 - case AFTER_PRE_FLOW_PHASE_1:
1.380 - for(NodeIt v(*g); v!=INVALID; ++v) {
1.381 - if (level[v] < n) {
1.382 - M.set(v, false);
1.383 - } else {
1.384 - M.set(v, true);
1.385 - }
1.386 - }
1.387 - break;
1.388 - case AFTER_PRE_FLOW_PHASE_2:
1.389 - case AFTER_NOTHING:
1.390 - case AFTER_AUGMENTING:
1.391 - case AFTER_FAST_AUGMENTING:
1.392 - minMinCut(M);
1.393 - break;
1.394 - }
1.395 - }
1.396 -
1.397 - ///Returns the inclusionwise minimum of the minimum value cuts.
1.398 -
1.399 - ///Sets \c M to the characteristic vector of the minimum value cut
1.400 - ///which is inclusionwise minimum. It is computed by processing
1.401 - ///a bfs from the source node \c s in the residual graph.
1.402 - ///\pre M should be a node map of bools initialized to false.
1.403 - ///\pre \c flow must be a maximum flow.
1.404 - template<typename _CutMap>
1.405 - void minMinCut(_CutMap& M) const {
1.406 - std::queue<Node> queue;
1.407 -
1.408 - M.set(s,true);
1.409 - queue.push(s);
1.410 -
1.411 - while (!queue.empty()) {
1.412 - Node w=queue.front();
1.413 - queue.pop();
1.414 -
1.415 - for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
1.416 - Node v=g->head(e);
1.417 - if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
1.418 - queue.push(v);
1.419 - M.set(v, true);
1.420 - }
1.421 - }
1.422 -
1.423 - for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
1.424 - Node v=g->tail(e);
1.425 - if (!M[v] && (*flow)[e] > 0 ) {
1.426 - queue.push(v);
1.427 - M.set(v, true);
1.428 - }
1.429 - }
1.430 - }
1.431 - }
1.432 -
1.433 - ///Returns the inclusionwise maximum of the minimum value cuts.
1.434 -
1.435 - ///Sets \c M to the characteristic vector of the minimum value cut
1.436 - ///which is inclusionwise maximum. It is computed by processing a
1.437 - ///backward bfs from the target node \c t in the residual graph.
1.438 - ///\pre M should be a node map of bools initialized to false.
1.439 - ///\pre \c flow must be a maximum flow.
1.440 - template<typename _CutMap>
1.441 - void maxMinCut(_CutMap& M) const {
1.442 -
1.443 - for(NodeIt v(*g) ; v!=INVALID; ++v) M.set(v, true);
1.444 -
1.445 - std::queue<Node> queue;
1.446 -
1.447 - M.set(t,false);
1.448 - queue.push(t);
1.449 -
1.450 - while (!queue.empty()) {
1.451 - Node w=queue.front();
1.452 - queue.pop();
1.453 -
1.454 - for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
1.455 - Node v=g->tail(e);
1.456 - if (M[v] && (*flow)[e] < (*capacity)[e] ) {
1.457 - queue.push(v);
1.458 - M.set(v, false);
1.459 - }
1.460 - }
1.461 -
1.462 - for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
1.463 - Node v=g->head(e);
1.464 - if (M[v] && (*flow)[e] > 0 ) {
1.465 - queue.push(v);
1.466 - M.set(v, false);
1.467 - }
1.468 - }
1.469 - }
1.470 - }
1.471 -
1.472 - ///Returns a minimum value cut.
1.473 -
1.474 - ///Sets \c M to the characteristic vector of a minimum value cut.
1.475 - ///\pre M should be a node map of bools initialized to false.
1.476 - ///\pre \c flow must be a maximum flow.
1.477 - template<typename CutMap>
1.478 - void minCut(CutMap& M) const { minMinCut(M); }
1.479 -
1.480 - ///Sets the source node to \c _s.
1.481 -
1.482 - ///Sets the source node to \c _s.
1.483 - ///
1.484 - void setSource(Node _s) { s=_s; status=AFTER_NOTHING; }
1.485 -
1.486 - ///Sets the target node to \c _t.
1.487 -
1.488 - ///Sets the target node to \c _t.
1.489 - ///
1.490 - void setTarget(Node _t) { t=_t; status=AFTER_NOTHING; }
1.491 -
1.492 - /// Sets the edge map of the capacities to _cap.
1.493 -
1.494 - /// Sets the edge map of the capacities to _cap.
1.495 - ///
1.496 - void setCap(const CapMap& _cap)
1.497 - { capacity=&_cap; status=AFTER_NOTHING; }
1.498 -
1.499 - /// Sets the edge map of the flows to _flow.
1.500 -
1.501 - /// Sets the edge map of the flows to _flow.
1.502 - ///
1.503 - void setFlow(FlowMap& _flow) { flow=&_flow; status=AFTER_NOTHING; }
1.504 -
1.505 -
1.506 - private:
1.507 -
1.508 - int push(Node w, NNMap& next, VecFirst& first) {
1.509 -
1.510 - int lev=level[w];
1.511 - Num exc=excess[w];
1.512 - int newlevel=n; //bound on the next level of w
1.513 -
1.514 - for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
1.515 - if ( (*flow)[e] >= (*capacity)[e] ) continue;
1.516 - Node v=g->head(e);
1.517 -
1.518 - if( lev > level[v] ) { //Push is allowed now
1.519 -
1.520 - if ( excess[v]<=0 && v!=t && v!=s ) {
1.521 - next.set(v,first[level[v]]);
1.522 - first[level[v]]=v;
1.523 - }
1.524 -
1.525 - Num cap=(*capacity)[e];
1.526 - Num flo=(*flow)[e];
1.527 - Num remcap=cap-flo;
1.528 -
1.529 - if ( remcap >= exc ) { //A nonsaturating push.
1.530 -
1.531 - flow->set(e, flo+exc);
1.532 - excess.set(v, excess[v]+exc);
1.533 - exc=0;
1.534 - break;
1.535 -
1.536 - } else { //A saturating push.
1.537 - flow->set(e, cap);
1.538 - excess.set(v, excess[v]+remcap);
1.539 - exc-=remcap;
1.540 - }
1.541 - } else if ( newlevel > level[v] ) newlevel = level[v];
1.542 - } //for out edges wv
1.543 -
1.544 - if ( exc > 0 ) {
1.545 - for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
1.546 -
1.547 - if( (*flow)[e] <= 0 ) continue;
1.548 - Node v=g->tail(e);
1.549 -
1.550 - if( lev > level[v] ) { //Push is allowed now
1.551 -
1.552 - if ( excess[v]<=0 && v!=t && v!=s ) {
1.553 - next.set(v,first[level[v]]);
1.554 - first[level[v]]=v;
1.555 - }
1.556 -
1.557 - Num flo=(*flow)[e];
1.558 -
1.559 - if ( flo >= exc ) { //A nonsaturating push.
1.560 -
1.561 - flow->set(e, flo-exc);
1.562 - excess.set(v, excess[v]+exc);
1.563 - exc=0;
1.564 - break;
1.565 - } else { //A saturating push.
1.566 -
1.567 - excess.set(v, excess[v]+flo);
1.568 - exc-=flo;
1.569 - flow->set(e,0);
1.570 - }
1.571 - } else if ( newlevel > level[v] ) newlevel = level[v];
1.572 - } //for in edges vw
1.573 -
1.574 - } // if w still has excess after the out edge for cycle
1.575 -
1.576 - excess.set(w, exc);
1.577 -
1.578 - return newlevel;
1.579 - }
1.580 -
1.581 -
1.582 -
1.583 - void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
1.584 - VecNode& level_list, NNMap& left, NNMap& right)
1.585 - {
1.586 - switch (fe) { //setting excess
1.587 - case NO_FLOW:
1.588 - for(EdgeIt e(*g); e!=INVALID; ++e) flow->set(e,0);
1.589 - for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
1.590 - break;
1.591 - case ZERO_FLOW:
1.592 - for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
1.593 - break;
1.594 - case GEN_FLOW:
1.595 - for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
1.596 - {
1.597 - Num exc=0;
1.598 - for(InEdgeIt e(*g,t) ; e!=INVALID; ++e) exc+=(*flow)[e];
1.599 - for(OutEdgeIt e(*g,t) ; e!=INVALID; ++e) exc-=(*flow)[e];
1.600 - excess.set(t,exc);
1.601 - }
1.602 - break;
1.603 - default:
1.604 - break;
1.605 - }
1.606 -
1.607 - for(NodeIt v(*g); v!=INVALID; ++v) level.set(v,n);
1.608 - //setting each node to level n
1.609 -
1.610 - std::queue<Node> bfs_queue;
1.611 -
1.612 -
1.613 - switch (fe) {
1.614 - case NO_FLOW: //flow is already set to const zero
1.615 - case ZERO_FLOW:
1.616 - //Reverse_bfs from t, to find the starting level.
1.617 - level.set(t,0);
1.618 - bfs_queue.push(t);
1.619 -
1.620 - while (!bfs_queue.empty()) {
1.621 -
1.622 - Node v=bfs_queue.front();
1.623 - bfs_queue.pop();
1.624 - int l=level[v]+1;
1.625 -
1.626 - for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
1.627 - Node w=g->tail(e);
1.628 - if ( level[w] == n && w != s ) {
1.629 - bfs_queue.push(w);
1.630 - Node z=level_list[l];
1.631 - if ( z!=INVALID ) left.set(z,w);
1.632 - right.set(w,z);
1.633 - level_list[l]=w;
1.634 - level.set(w, l);
1.635 - }
1.636 - }
1.637 - }
1.638 -
1.639 - //the starting flow
1.640 - for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e)
1.641 - {
1.642 - Num c=(*capacity)[e];
1.643 - if ( c <= 0 ) continue;
1.644 - Node w=g->head(e);
1.645 - if ( level[w] < n ) {
1.646 - if ( excess[w] <= 0 && w!=t ) //putting into the stack
1.647 - {
1.648 - next.set(w,first[level[w]]);
1.649 - first[level[w]]=w;
1.650 - }
1.651 - flow->set(e, c);
1.652 - excess.set(w, excess[w]+c);
1.653 - }
1.654 - }
1.655 - break;
1.656 - case GEN_FLOW:
1.657 - //Reverse_bfs from t in the residual graph,
1.658 - //to find the starting level.
1.659 - level.set(t,0);
1.660 - bfs_queue.push(t);
1.661 -
1.662 - while (!bfs_queue.empty()) {
1.663 -
1.664 - Node v=bfs_queue.front();
1.665 - bfs_queue.pop();
1.666 - int l=level[v]+1;
1.667 -
1.668 - for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
1.669 - if ( (*capacity)[e] <= (*flow)[e] ) continue;
1.670 - Node w=g->tail(e);
1.671 - if ( level[w] == n && w != s ) {
1.672 - bfs_queue.push(w);
1.673 - Node z=level_list[l];
1.674 - if ( z!=INVALID ) left.set(z,w);
1.675 - right.set(w,z);
1.676 - level_list[l]=w;
1.677 - level.set(w, l);
1.678 - }
1.679 - }
1.680 -
1.681 - for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) {
1.682 - if ( 0 >= (*flow)[e] ) continue;
1.683 - Node w=g->head(e);
1.684 - if ( level[w] == n && w != s ) {
1.685 - bfs_queue.push(w);
1.686 - Node z=level_list[l];
1.687 - if ( z!=INVALID ) left.set(z,w);
1.688 - right.set(w,z);
1.689 - level_list[l]=w;
1.690 - level.set(w, l);
1.691 - }
1.692 - }
1.693 - }
1.694 -
1.695 - //the starting flow
1.696 - for(OutEdgeIt e(*g,s); e!=INVALID; ++e)
1.697 - {
1.698 - Num rem=(*capacity)[e]-(*flow)[e];
1.699 - if ( rem <= 0 ) continue;
1.700 - Node w=g->head(e);
1.701 - if ( level[w] < n ) {
1.702 - if ( excess[w] <= 0 && w!=t ) //putting into the stack
1.703 - {
1.704 - next.set(w,first[level[w]]);
1.705 - first[level[w]]=w;
1.706 - }
1.707 - flow->set(e, (*capacity)[e]);
1.708 - excess.set(w, excess[w]+rem);
1.709 - }
1.710 - }
1.711 -
1.712 - for(InEdgeIt e(*g,s); e!=INVALID; ++e)
1.713 - {
1.714 - if ( (*flow)[e] <= 0 ) continue;
1.715 - Node w=g->tail(e);
1.716 - if ( level[w] < n ) {
1.717 - if ( excess[w] <= 0 && w!=t )
1.718 - {
1.719 - next.set(w,first[level[w]]);
1.720 - first[level[w]]=w;
1.721 - }
1.722 - excess.set(w, excess[w]+(*flow)[e]);
1.723 - flow->set(e, 0);
1.724 - }
1.725 - }
1.726 - break;
1.727 - case PRE_FLOW:
1.728 - //Reverse_bfs from t in the residual graph,
1.729 - //to find the starting level.
1.730 - level.set(t,0);
1.731 - bfs_queue.push(t);
1.732 -
1.733 - while (!bfs_queue.empty()) {
1.734 -
1.735 - Node v=bfs_queue.front();
1.736 - bfs_queue.pop();
1.737 - int l=level[v]+1;
1.738 -
1.739 - for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
1.740 - if ( (*capacity)[e] <= (*flow)[e] ) continue;
1.741 - Node w=g->tail(e);
1.742 - if ( level[w] == n && w != s ) {
1.743 - bfs_queue.push(w);
1.744 - Node z=level_list[l];
1.745 - if ( z!=INVALID ) left.set(z,w);
1.746 - right.set(w,z);
1.747 - level_list[l]=w;
1.748 - level.set(w, l);
1.749 - }
1.750 - }
1.751 -
1.752 - for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) {
1.753 - if ( 0 >= (*flow)[e] ) continue;
1.754 - Node w=g->head(e);
1.755 - if ( level[w] == n && w != s ) {
1.756 - bfs_queue.push(w);
1.757 - Node z=level_list[l];
1.758 - if ( z!=INVALID ) left.set(z,w);
1.759 - right.set(w,z);
1.760 - level_list[l]=w;
1.761 - level.set(w, l);
1.762 - }
1.763 - }
1.764 - }
1.765 -
1.766 -
1.767 - //the starting flow
1.768 - for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
1.769 - Num rem=(*capacity)[e]-(*flow)[e];
1.770 - if ( rem <= 0 ) continue;
1.771 - Node w=g->head(e);
1.772 - if ( level[w] < n ) {
1.773 - flow->set(e, (*capacity)[e]);
1.774 - excess.set(w, excess[w]+rem);
1.775 - }
1.776 - }
1.777 -
1.778 - for(InEdgeIt e(*g,s) ; e!=INVALID; ++e) {
1.779 - if ( (*flow)[e] <= 0 ) continue;
1.780 - Node w=g->tail(e);
1.781 - if ( level[w] < n ) {
1.782 - excess.set(w, excess[w]+(*flow)[e]);
1.783 - flow->set(e, 0);
1.784 - }
1.785 - }
1.786 -
1.787 - //computing the excess
1.788 - for(NodeIt w(*g); w!=INVALID; ++w) {
1.789 - Num exc=0;
1.790 -
1.791 - for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) exc+=(*flow)[e];
1.792 - for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) exc-=(*flow)[e];
1.793 -
1.794 - excess.set(w,exc);
1.795 -
1.796 - //putting the active nodes into the stack
1.797 - int lev=level[w];
1.798 - if ( exc > 0 && lev < n && Node(w) != t )
1.799 - ///\bug if ( exc > 0 && lev < n && w != t ) temporarily for working with wrappers.
1.800 - {
1.801 - next.set(w,first[lev]);
1.802 - first[lev]=w;
1.803 - }
1.804 - }
1.805 - break;
1.806 - } //switch
1.807 - } //preflowPreproc
1.808 -
1.809 -
1.810 - void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
1.811 - VecNode& level_list, NNMap& left,
1.812 - NNMap& right, int& b, int& k, bool what_heur )
1.813 - {
1.814 -
1.815 - int lev=level[w];
1.816 -
1.817 - Node right_n=right[w];
1.818 - Node left_n=left[w];
1.819 -
1.820 - //unlacing starts
1.821 - if ( right_n!=INVALID ) {
1.822 - if ( left_n!=INVALID ) {
1.823 - right.set(left_n, right_n);
1.824 - left.set(right_n, left_n);
1.825 - } else {
1.826 - level_list[lev]=right_n;
1.827 - left.set(right_n, INVALID);
1.828 - }
1.829 - } else {
1.830 - if ( left_n!=INVALID ) {
1.831 - right.set(left_n, INVALID);
1.832 - } else {
1.833 - level_list[lev]=INVALID;
1.834 - }
1.835 - }
1.836 - //unlacing ends
1.837 -
1.838 - if ( level_list[lev]==INVALID ) {
1.839 -
1.840 - //gapping starts
1.841 - for (int i=lev; i!=k ; ) {
1.842 - Node v=level_list[++i];
1.843 - while ( v!=INVALID ) {
1.844 - level.set(v,n);
1.845 - v=right[v];
1.846 - }
1.847 - level_list[i]=INVALID;
1.848 - if ( !what_heur ) first[i]=INVALID;
1.849 - }
1.850 -
1.851 - level.set(w,n);
1.852 - b=lev-1;
1.853 - k=b;
1.854 - //gapping ends
1.855 -
1.856 - } else {
1.857 -
1.858 - if ( newlevel == n ) level.set(w,n);
1.859 - else {
1.860 - level.set(w,++newlevel);
1.861 - next.set(w,first[newlevel]);
1.862 - first[newlevel]=w;
1.863 - if ( what_heur ) b=newlevel;
1.864 - if ( k < newlevel ) ++k; //now k=newlevel
1.865 - Node z=level_list[newlevel];
1.866 - if ( z!=INVALID ) left.set(z,w);
1.867 - right.set(w,z);
1.868 - left.set(w,INVALID);
1.869 - level_list[newlevel]=w;
1.870 - }
1.871 - }
1.872 - } //relabel
1.873 -
1.874 - void printexcess() {////
1.875 - std::cout << "Excesses:" <<std::endl;
1.876 -
1.877 - for(NodeIt v(*g); v!=INVALID ; ++v) {
1.878 - std::cout << 1+(g->id(v)) << ":" << excess[v]<<std::endl;
1.879 - }
1.880 - }
1.881 -
1.882 - void printlevel() {////
1.883 - std::cout << "Levels:" <<std::endl;
1.884 -
1.885 - for(NodeIt v(*g); v!=INVALID ; ++v) {
1.886 - std::cout << 1+(g->id(v)) << ":" << level[v]<<std::endl;
1.887 - }
1.888 - }
1.889 -
1.890 - void printactive() {////
1.891 - std::cout << "Levels:" <<std::endl;
1.892 -
1.893 - for(NodeIt v(*g); v!=INVALID ; ++v) {
1.894 - std::cout << 1+(g->id(v)) << ":" << level[v]<<std::endl;
1.895 - }
1.896 - }
1.897 -
1.898 -
1.899 - }; //class MaxFlow
1.900 -} //namespace hugo
1.901 -
1.902 -#endif //HUGO_MAX_FLOW_H
1.903 -
1.904 -
1.905 -
1.906 -
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/src/hugo/preflow.h Mon Sep 13 13:57:13 2004 +0000
2.3 @@ -0,0 +1,796 @@
2.4 +// -*- C++ -*-
2.5 +#ifndef HUGO_PREFLOW_H
2.6 +#define HUGO_PREFLOW_H
2.7 +
2.8 +#include <vector>
2.9 +#include <queue>
2.10 +
2.11 +#include <hugo/invalid.h>
2.12 +#include <hugo/maps.h>
2.13 +
2.14 +/// \file
2.15 +/// \ingroup flowalgs
2.16 +
2.17 +namespace hugo {
2.18 +
2.19 + /// \addtogroup flowalgs
2.20 + /// @{
2.21 +
2.22 + ///Preflow algorithms class.
2.23 +
2.24 + ///This class provides an implementation of the \e preflow \e
2.25 + ///algorithm producing a flow of maximum value in a directed
2.26 + ///graph. The preflow algorithms are the fastest max flow algorithms
2.27 + ///up-to-date. The \e source node, the \e target node, the \e
2.28 + ///capacity of the edges and the \e starting \e flow value of the
2.29 + ///edges should be passed to the algorithm through the
2.30 + ///constructor. It is possible to change these quantities using the
2.31 + ///functions \ref setSource, \ref setTarget, \ref setCap and \ref
2.32 + ///setFlow.
2.33 + ///
2.34 + ///After running \c phase1 or \c preflow, the actual flow
2.35 + ///value can be obtained by calling \ref flowValue(). The minimum
2.36 + ///value cut can be written into a \c node map of \c bools by
2.37 + ///calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
2.38 + ///the inclusionwise minimum and maximum of the minimum value cuts,
2.39 + ///resp.)
2.40 + ///
2.41 + ///\param Graph The directed graph type the algorithm runs on.
2.42 + ///\param Num The number type of the capacities and the flow values.
2.43 + ///\param CapMap The capacity map type.
2.44 + ///\param FlowMap The flow map type.
2.45 + ///
2.46 + ///\author Jacint Szabo
2.47 + template <typename Graph, typename Num,
2.48 + typename CapMap=typename Graph::template EdgeMap<Num>,
2.49 + typename FlowMap=typename Graph::template EdgeMap<Num> >
2.50 + class Preflow {
2.51 + protected:
2.52 + typedef typename Graph::Node Node;
2.53 + typedef typename Graph::NodeIt NodeIt;
2.54 + typedef typename Graph::EdgeIt EdgeIt;
2.55 + typedef typename Graph::OutEdgeIt OutEdgeIt;
2.56 + typedef typename Graph::InEdgeIt InEdgeIt;
2.57 +
2.58 + typedef typename Graph::template NodeMap<Node> NNMap;
2.59 + typedef typename std::vector<Node> VecNode;
2.60 +
2.61 + const Graph* g;
2.62 + Node s;
2.63 + Node t;
2.64 + const CapMap* capacity;
2.65 + FlowMap* flow;
2.66 + int n; //the number of nodes of G
2.67 +
2.68 + typename Graph::template NodeMap<int> level;
2.69 + typename Graph::template NodeMap<Num> excess;
2.70 +
2.71 + // constants used for heuristics
2.72 + static const int H0=20;
2.73 + static const int H1=1;
2.74 +
2.75 + public:
2.76 +
2.77 + ///Indicates the property of the starting flow map.
2.78 +
2.79 + ///Indicates the property of the starting flow map. The meanings are as follows:
2.80 + ///- \c ZERO_FLOW: constant zero flow
2.81 + ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
2.82 + ///the sum of the out-flows in every node except the \e source and
2.83 + ///the \e target.
2.84 + ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
2.85 + ///least the sum of the out-flows in every node except the \e source.
2.86 + ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be
2.87 + ///set to the constant zero flow in the beginning of the algorithm in this case.
2.88 + ///
2.89 + enum FlowEnum{
2.90 + NO_FLOW,
2.91 + ZERO_FLOW,
2.92 + GEN_FLOW,
2.93 + PRE_FLOW
2.94 + };
2.95 +
2.96 + ///Indicates the state of the preflow algorithm.
2.97 +
2.98 + ///Indicates the state of the preflow algorithm. The meanings are as follows:
2.99 + ///- \c AFTER_NOTHING: before running the algorithm or at an unspecified state.
2.100 + ///- \c AFTER_PREFLOW_PHASE_1: right after running \c phase1
2.101 + ///- \c AFTER_PREFLOW_PHASE_2: after running \ref phase2()
2.102 + ///
2.103 + enum StatusEnum {
2.104 + AFTER_NOTHING,
2.105 + AFTER_PREFLOW_PHASE_1,
2.106 + AFTER_PREFLOW_PHASE_2
2.107 + };
2.108 +
2.109 + protected:
2.110 + FlowEnum flow_prop;
2.111 + StatusEnum status; // Do not needle this flag only if necessary.
2.112 +
2.113 + public:
2.114 + ///The constructor of the class.
2.115 +
2.116 + ///The constructor of the class.
2.117 + ///\param _G The directed graph the algorithm runs on.
2.118 + ///\param _s The source node.
2.119 + ///\param _t The target node.
2.120 + ///\param _capacity The capacity of the edges.
2.121 + ///\param _flow The flow of the edges.
2.122 + ///Except the graph, all of these parameters can be reset by
2.123 + ///calling \ref setSource, \ref setTarget, \ref setCap and \ref
2.124 + ///setFlow, resp.
2.125 + Preflow(const Graph& _G, Node _s, Node _t,
2.126 + const CapMap& _capacity, FlowMap& _flow) :
2.127 + g(&_G), s(_s), t(_t), capacity(&_capacity),
2.128 + flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0),
2.129 + flow_prop(NO_FLOW), status(AFTER_NOTHING) { }
2.130 +
2.131 +
2.132 +
2.133 + ///Runs the preflow algorithm.
2.134 +
2.135 + ///Runs the preflow algorithm.
2.136 + void run() {
2.137 + phase1(flow_prop);
2.138 + phase2();
2.139 + }
2.140 +
2.141 + ///Runs the preflow algorithm.
2.142 +
2.143 + ///Runs the preflow algorithm.
2.144 + ///\pre The starting flow map must be
2.145 + /// - a constant zero flow if \c fp is \c ZERO_FLOW,
2.146 + /// - an arbitrary flow if \c fp is \c GEN_FLOW,
2.147 + /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
2.148 + /// - any map if \c fp is NO_FLOW.
2.149 + ///If the starting flow map is a flow or a preflow then
2.150 + ///the algorithm terminates faster.
2.151 + void run(FlowEnum fp) {
2.152 + flow_prop=fp;
2.153 + run();
2.154 + }
2.155 +
2.156 + ///Runs the first phase of the preflow algorithm.
2.157 +
2.158 + ///The preflow algorithm consists of two phases, this method runs the
2.159 + ///first phase. After the first phase the maximum flow value and a
2.160 + ///minimum value cut can already be computed, though a maximum flow
2.161 + ///is not yet obtained. So after calling this method \ref flowValue
2.162 + ///and \ref minCut gives proper results.
2.163 + ///\warning: \ref minMinCut and \ref maxMinCut do not
2.164 + ///give minimum value cuts unless calling \ref phase2.
2.165 + ///\pre The starting flow must be
2.166 + /// - a constant zero flow if \c fp is \c ZERO_FLOW,
2.167 + /// - an arbitary flow if \c fp is \c GEN_FLOW,
2.168 + /// - an arbitary preflow if \c fp is \c PRE_FLOW,
2.169 + /// - any map if \c fp is NO_FLOW.
2.170 + void phase1(FlowEnum fp)
2.171 + {
2.172 + flow_prop=fp;
2.173 + phase1();
2.174 + }
2.175 +
2.176 +
2.177 + ///Runs the first phase of the preflow algorithm.
2.178 +
2.179 + ///The preflow algorithm consists of two phases, this method runs the
2.180 + ///first phase. After the first phase the maximum flow value and a
2.181 + ///minimum value cut can already be computed, though a maximum flow
2.182 + ///is not yet obtained. So after calling this method \ref flowValue
2.183 + ///and \ref actMinCut gives proper results.
2.184 + ///\warning: \ref minCut, \ref minMinCut and \ref maxMinCut do not
2.185 + ///give minimum value cuts unless calling \ref phase2.
2.186 + void phase1()
2.187 + {
2.188 + int heur0=(int)(H0*n); //time while running 'bound decrease'
2.189 + int heur1=(int)(H1*n); //time while running 'highest label'
2.190 + int heur=heur1; //starting time interval (#of relabels)
2.191 + int numrelabel=0;
2.192 +
2.193 + bool what_heur=1;
2.194 + //It is 0 in case 'bound decrease' and 1 in case 'highest label'
2.195 +
2.196 + bool end=false;
2.197 + //Needed for 'bound decrease', true means no active
2.198 + //nodes are above bound b.
2.199 +
2.200 + int k=n-2; //bound on the highest level under n containing a node
2.201 + int b=k; //bound on the highest level under n of an active node
2.202 +
2.203 + VecNode first(n, INVALID);
2.204 + NNMap next(*g, INVALID);
2.205 +
2.206 + NNMap left(*g, INVALID);
2.207 + NNMap right(*g, INVALID);
2.208 + VecNode level_list(n,INVALID);
2.209 + //List of the nodes in level i<n, set to n.
2.210 +
2.211 + preflowPreproc(first, next, level_list, left, right);
2.212 +
2.213 + //Push/relabel on the highest level active nodes.
2.214 + while ( true ) {
2.215 + if ( b == 0 ) {
2.216 + if ( !what_heur && !end && k > 0 ) {
2.217 + b=k;
2.218 + end=true;
2.219 + } else break;
2.220 + }
2.221 +
2.222 + if ( first[b]==INVALID ) --b;
2.223 + else {
2.224 + end=false;
2.225 + Node w=first[b];
2.226 + first[b]=next[w];
2.227 + int newlevel=push(w, next, first);
2.228 + if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list,
2.229 + left, right, b, k, what_heur);
2.230 +
2.231 + ++numrelabel;
2.232 + if ( numrelabel >= heur ) {
2.233 + numrelabel=0;
2.234 + if ( what_heur ) {
2.235 + what_heur=0;
2.236 + heur=heur0;
2.237 + end=false;
2.238 + } else {
2.239 + what_heur=1;
2.240 + heur=heur1;
2.241 + b=k;
2.242 + }
2.243 + }
2.244 + }
2.245 + }
2.246 + flow_prop=PRE_FLOW;
2.247 + status=AFTER_PREFLOW_PHASE_1;
2.248 + }
2.249 + // Heuristics:
2.250 + // 2 phase
2.251 + // gap
2.252 + // list 'level_list' on the nodes on level i implemented by hand
2.253 + // stack 'active' on the active nodes on level i
2.254 + // runs heuristic 'highest label' for H1*n relabels
2.255 + // runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
2.256 + // Parameters H0 and H1 are initialized to 20 and 1.
2.257 +
2.258 +
2.259 + ///Runs the second phase of the preflow algorithm.
2.260 +
2.261 + ///The preflow algorithm consists of two phases, this method runs
2.262 + ///the second phase. After calling \ref phase1 and then
2.263 + ///\ref phase2 the methods \ref flowValue, \ref minCut,
2.264 + ///\ref minMinCut and \ref maxMinCut give proper results.
2.265 + ///\pre \ref phase1 must be called before.
2.266 + void phase2()
2.267 + {
2.268 +
2.269 + int k=n-2; //bound on the highest level under n containing a node
2.270 + int b=k; //bound on the highest level under n of an active node
2.271 +
2.272 +
2.273 + VecNode first(n, INVALID);
2.274 + NNMap next(*g, INVALID);
2.275 + level.set(s,0);
2.276 + std::queue<Node> bfs_queue;
2.277 + bfs_queue.push(s);
2.278 +
2.279 + while ( !bfs_queue.empty() ) {
2.280 +
2.281 + Node v=bfs_queue.front();
2.282 + bfs_queue.pop();
2.283 + int l=level[v]+1;
2.284 +
2.285 + for(InEdgeIt e(*g,v); e!=INVALID; ++e) {
2.286 + if ( (*capacity)[e] <= (*flow)[e] ) continue;
2.287 + Node u=g->tail(e);
2.288 + if ( level[u] >= n ) {
2.289 + bfs_queue.push(u);
2.290 + level.set(u, l);
2.291 + if ( excess[u] > 0 ) {
2.292 + next.set(u,first[l]);
2.293 + first[l]=u;
2.294 + }
2.295 + }
2.296 + }
2.297 +
2.298 + for(OutEdgeIt e(*g,v); e!=INVALID; ++e) {
2.299 + if ( 0 >= (*flow)[e] ) continue;
2.300 + Node u=g->head(e);
2.301 + if ( level[u] >= n ) {
2.302 + bfs_queue.push(u);
2.303 + level.set(u, l);
2.304 + if ( excess[u] > 0 ) {
2.305 + next.set(u,first[l]);
2.306 + first[l]=u;
2.307 + }
2.308 + }
2.309 + }
2.310 + }
2.311 + b=n-2;
2.312 +
2.313 + while ( true ) {
2.314 +
2.315 + if ( b == 0 ) break;
2.316 + if ( first[b]==INVALID ) --b;
2.317 + else {
2.318 + Node w=first[b];
2.319 + first[b]=next[w];
2.320 + int newlevel=push(w,next, first);
2.321 +
2.322 + //relabel
2.323 + if ( excess[w] > 0 ) {
2.324 + level.set(w,++newlevel);
2.325 + next.set(w,first[newlevel]);
2.326 + first[newlevel]=w;
2.327 + b=newlevel;
2.328 + }
2.329 + }
2.330 + } // while(true)
2.331 + flow_prop=GEN_FLOW;
2.332 + status=AFTER_PREFLOW_PHASE_2;
2.333 + }
2.334 +
2.335 + /// Returns the value of the maximum flow.
2.336 +
2.337 + /// Returns the value of the maximum flow by returning the excess
2.338 + /// of the target node \ref t. This value equals to the value of
2.339 + /// the maximum flow already after running \ref phase1.
2.340 + Num flowValue() const {
2.341 + return excess[t];
2.342 + }
2.343 +
2.344 +
2.345 + ///Returns a minimum value cut.
2.346 +
2.347 + ///Sets \c M to the characteristic vector of a minimum value
2.348 + ///cut. This method can be called both after running \ref
2.349 + ///phase1 and \ref phase2. It is much faster after
2.350 + ///\ref phase1. \pre M should be a node map of bools. \pre
2.351 + ///If \ref mincut is called after \ref phase2 then M should
2.352 + ///be initialized to false.
2.353 + template<typename _CutMap>
2.354 + void minCut(_CutMap& M) const {
2.355 + switch ( status ) {
2.356 + case AFTER_PREFLOW_PHASE_1:
2.357 + for(NodeIt v(*g); v!=INVALID; ++v) {
2.358 + if (level[v] < n) {
2.359 + M.set(v, false);
2.360 + } else {
2.361 + M.set(v, true);
2.362 + }
2.363 + }
2.364 + break;
2.365 + case AFTER_PREFLOW_PHASE_2:
2.366 + minMinCut(M);
2.367 + break;
2.368 + case AFTER_NOTHING:
2.369 + break;
2.370 + }
2.371 + }
2.372 +
2.373 + ///Returns the inclusionwise minimum of the minimum value cuts.
2.374 +
2.375 + ///Sets \c M to the characteristic vector of the minimum value cut
2.376 + ///which is inclusionwise minimum. It is computed by processing a
2.377 + ///bfs from the source node \c s in the residual graph. \pre M
2.378 + ///should be a node map of bools initialized to false. \pre \ref
2.379 + ///phase2 should already be run.
2.380 + template<typename _CutMap>
2.381 + void minMinCut(_CutMap& M) const {
2.382 +
2.383 + std::queue<Node> queue;
2.384 + M.set(s,true);
2.385 + queue.push(s);
2.386 +
2.387 + while (!queue.empty()) {
2.388 + Node w=queue.front();
2.389 + queue.pop();
2.390 +
2.391 + for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
2.392 + Node v=g->head(e);
2.393 + if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
2.394 + queue.push(v);
2.395 + M.set(v, true);
2.396 + }
2.397 + }
2.398 +
2.399 + for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
2.400 + Node v=g->tail(e);
2.401 + if (!M[v] && (*flow)[e] > 0 ) {
2.402 + queue.push(v);
2.403 + M.set(v, true);
2.404 + }
2.405 + }
2.406 + }
2.407 + }
2.408 +
2.409 + ///Returns the inclusionwise maximum of the minimum value cuts.
2.410 +
2.411 + ///Sets \c M to the characteristic vector of the minimum value cut
2.412 + ///which is inclusionwise maximum. It is computed by processing a
2.413 + ///backward bfs from the target node \c t in the residual graph.
2.414 + ///\pre \ref phase2() or preflow() should already be run.
2.415 + template<typename _CutMap>
2.416 + void maxMinCut(_CutMap& M) const {
2.417 +
2.418 + for(NodeIt v(*g) ; v!=INVALID; ++v) M.set(v, true);
2.419 +
2.420 + std::queue<Node> queue;
2.421 +
2.422 + M.set(t,false);
2.423 + queue.push(t);
2.424 +
2.425 + while (!queue.empty()) {
2.426 + Node w=queue.front();
2.427 + queue.pop();
2.428 +
2.429 + for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
2.430 + Node v=g->tail(e);
2.431 + if (M[v] && (*flow)[e] < (*capacity)[e] ) {
2.432 + queue.push(v);
2.433 + M.set(v, false);
2.434 + }
2.435 + }
2.436 +
2.437 + for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
2.438 + Node v=g->head(e);
2.439 + if (M[v] && (*flow)[e] > 0 ) {
2.440 + queue.push(v);
2.441 + M.set(v, false);
2.442 + }
2.443 + }
2.444 + }
2.445 + }
2.446 +
2.447 + ///Sets the source node to \c _s.
2.448 +
2.449 + ///Sets the source node to \c _s.
2.450 + ///
2.451 + void setSource(Node _s) {
2.452 + s=_s;
2.453 + if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
2.454 + status=AFTER_NOTHING;
2.455 + }
2.456 +
2.457 + ///Sets the target node to \c _t.
2.458 +
2.459 + ///Sets the target node to \c _t.
2.460 + ///
2.461 + void setTarget(Node _t) {
2.462 + t=_t;
2.463 + if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
2.464 + status=AFTER_NOTHING;
2.465 + }
2.466 +
2.467 + /// Sets the edge map of the capacities to _cap.
2.468 +
2.469 + /// Sets the edge map of the capacities to _cap.
2.470 + ///
2.471 + void setCap(const CapMap& _cap) {
2.472 + capacity=&_cap;
2.473 + status=AFTER_NOTHING;
2.474 + }
2.475 +
2.476 + /// Sets the edge map of the flows to _flow.
2.477 +
2.478 + /// Sets the edge map of the flows to _flow.
2.479 + ///
2.480 + void setFlow(FlowMap& _flow) {
2.481 + flow=&_flow;
2.482 + flow_prop=NO_FLOW;
2.483 + status=AFTER_NOTHING;
2.484 + }
2.485 +
2.486 +
2.487 + private:
2.488 +
2.489 + int push(Node w, NNMap& next, VecNode& first) {
2.490 +
2.491 + int lev=level[w];
2.492 + Num exc=excess[w];
2.493 + int newlevel=n; //bound on the next level of w
2.494 +
2.495 + for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
2.496 + if ( (*flow)[e] >= (*capacity)[e] ) continue;
2.497 + Node v=g->head(e);
2.498 +
2.499 + if( lev > level[v] ) { //Push is allowed now
2.500 +
2.501 + if ( excess[v]<=0 && v!=t && v!=s ) {
2.502 + next.set(v,first[level[v]]);
2.503 + first[level[v]]=v;
2.504 + }
2.505 +
2.506 + Num cap=(*capacity)[e];
2.507 + Num flo=(*flow)[e];
2.508 + Num remcap=cap-flo;
2.509 +
2.510 + if ( remcap >= exc ) { //A nonsaturating push.
2.511 +
2.512 + flow->set(e, flo+exc);
2.513 + excess.set(v, excess[v]+exc);
2.514 + exc=0;
2.515 + break;
2.516 +
2.517 + } else { //A saturating push.
2.518 + flow->set(e, cap);
2.519 + excess.set(v, excess[v]+remcap);
2.520 + exc-=remcap;
2.521 + }
2.522 + } else if ( newlevel > level[v] ) newlevel = level[v];
2.523 + } //for out edges wv
2.524 +
2.525 + if ( exc > 0 ) {
2.526 + for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
2.527 +
2.528 + if( (*flow)[e] <= 0 ) continue;
2.529 + Node v=g->tail(e);
2.530 +
2.531 + if( lev > level[v] ) { //Push is allowed now
2.532 +
2.533 + if ( excess[v]<=0 && v!=t && v!=s ) {
2.534 + next.set(v,first[level[v]]);
2.535 + first[level[v]]=v;
2.536 + }
2.537 +
2.538 + Num flo=(*flow)[e];
2.539 +
2.540 + if ( flo >= exc ) { //A nonsaturating push.
2.541 +
2.542 + flow->set(e, flo-exc);
2.543 + excess.set(v, excess[v]+exc);
2.544 + exc=0;
2.545 + break;
2.546 + } else { //A saturating push.
2.547 +
2.548 + excess.set(v, excess[v]+flo);
2.549 + exc-=flo;
2.550 + flow->set(e,0);
2.551 + }
2.552 + } else if ( newlevel > level[v] ) newlevel = level[v];
2.553 + } //for in edges vw
2.554 +
2.555 + } // if w still has excess after the out edge for cycle
2.556 +
2.557 + excess.set(w, exc);
2.558 +
2.559 + return newlevel;
2.560 + }
2.561 +
2.562 +
2.563 +
2.564 + void preflowPreproc(VecNode& first, NNMap& next,
2.565 + VecNode& level_list, NNMap& left, NNMap& right)
2.566 + {
2.567 + for(NodeIt v(*g); v!=INVALID; ++v) level.set(v,n);
2.568 + std::queue<Node> bfs_queue;
2.569 +
2.570 + if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
2.571 + //Reverse_bfs from t in the residual graph,
2.572 + //to find the starting level.
2.573 + level.set(t,0);
2.574 + bfs_queue.push(t);
2.575 +
2.576 + while ( !bfs_queue.empty() ) {
2.577 +
2.578 + Node v=bfs_queue.front();
2.579 + bfs_queue.pop();
2.580 + int l=level[v]+1;
2.581 +
2.582 + for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
2.583 + if ( (*capacity)[e] <= (*flow)[e] ) continue;
2.584 + Node w=g->tail(e);
2.585 + if ( level[w] == n && w != s ) {
2.586 + bfs_queue.push(w);
2.587 + Node z=level_list[l];
2.588 + if ( z!=INVALID ) left.set(z,w);
2.589 + right.set(w,z);
2.590 + level_list[l]=w;
2.591 + level.set(w, l);
2.592 + }
2.593 + }
2.594 +
2.595 + for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) {
2.596 + if ( 0 >= (*flow)[e] ) continue;
2.597 + Node w=g->head(e);
2.598 + if ( level[w] == n && w != s ) {
2.599 + bfs_queue.push(w);
2.600 + Node z=level_list[l];
2.601 + if ( z!=INVALID ) 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 + } //while
2.608 + } //if
2.609 +
2.610 +
2.611 + switch (flow_prop) {
2.612 + case NO_FLOW:
2.613 + for(EdgeIt e(*g); e!=INVALID; ++e) flow->set(e,0);
2.614 + case ZERO_FLOW:
2.615 + for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
2.616 +
2.617 + //Reverse_bfs from t, to find the starting level.
2.618 + level.set(t,0);
2.619 + bfs_queue.push(t);
2.620 +
2.621 + while ( !bfs_queue.empty() ) {
2.622 +
2.623 + Node v=bfs_queue.front();
2.624 + bfs_queue.pop();
2.625 + int l=level[v]+1;
2.626 +
2.627 + for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
2.628 + Node w=g->tail(e);
2.629 + if ( level[w] == n && w != s ) {
2.630 + bfs_queue.push(w);
2.631 + Node z=level_list[l];
2.632 + if ( z!=INVALID ) left.set(z,w);
2.633 + right.set(w,z);
2.634 + level_list[l]=w;
2.635 + level.set(w, l);
2.636 + }
2.637 + }
2.638 + }
2.639 +
2.640 + //the starting flow
2.641 + for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
2.642 + Num c=(*capacity)[e];
2.643 + if ( c <= 0 ) continue;
2.644 + Node w=g->head(e);
2.645 + if ( level[w] < n ) {
2.646 + if ( excess[w] <= 0 && w!=t ) { //putting into the stack
2.647 + next.set(w,first[level[w]]);
2.648 + first[level[w]]=w;
2.649 + }
2.650 + flow->set(e, c);
2.651 + excess.set(w, excess[w]+c);
2.652 + }
2.653 + }
2.654 + break;
2.655 +
2.656 + case GEN_FLOW:
2.657 + for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
2.658 + {
2.659 + Num exc=0;
2.660 + for(InEdgeIt e(*g,t) ; e!=INVALID; ++e) exc+=(*flow)[e];
2.661 + for(OutEdgeIt e(*g,t) ; e!=INVALID; ++e) exc-=(*flow)[e];
2.662 + excess.set(t,exc);
2.663 + }
2.664 +
2.665 + //the starting flow
2.666 + for(OutEdgeIt e(*g,s); e!=INVALID; ++e) {
2.667 + Num rem=(*capacity)[e]-(*flow)[e];
2.668 + if ( rem <= 0 ) continue;
2.669 + Node w=g->head(e);
2.670 + if ( level[w] < n ) {
2.671 + if ( excess[w] <= 0 && w!=t ) { //putting into the stack
2.672 + next.set(w,first[level[w]]);
2.673 + first[level[w]]=w;
2.674 + }
2.675 + flow->set(e, (*capacity)[e]);
2.676 + excess.set(w, excess[w]+rem);
2.677 + }
2.678 + }
2.679 +
2.680 + for(InEdgeIt e(*g,s); e!=INVALID; ++e) {
2.681 + if ( (*flow)[e] <= 0 ) continue;
2.682 + Node w=g->tail(e);
2.683 + if ( level[w] < n ) {
2.684 + if ( excess[w] <= 0 && w!=t ) {
2.685 + next.set(w,first[level[w]]);
2.686 + first[level[w]]=w;
2.687 + }
2.688 + excess.set(w, excess[w]+(*flow)[e]);
2.689 + flow->set(e, 0);
2.690 + }
2.691 + }
2.692 + break;
2.693 +
2.694 + case PRE_FLOW:
2.695 + //the starting flow
2.696 + for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
2.697 + Num rem=(*capacity)[e]-(*flow)[e];
2.698 + if ( rem <= 0 ) continue;
2.699 + Node w=g->head(e);
2.700 + if ( level[w] < n ) flow->set(e, (*capacity)[e]);
2.701 + }
2.702 +
2.703 + for(InEdgeIt e(*g,s) ; e!=INVALID; ++e) {
2.704 + if ( (*flow)[e] <= 0 ) continue;
2.705 + Node w=g->tail(e);
2.706 + if ( level[w] < n ) flow->set(e, 0);
2.707 + }
2.708 +
2.709 + //computing the excess
2.710 + for(NodeIt w(*g); w!=INVALID; ++w) {
2.711 + Num exc=0;
2.712 + for(InEdgeIt e(*g,w); e!=INVALID; ++e) exc+=(*flow)[e];
2.713 + for(OutEdgeIt e(*g,w); e!=INVALID; ++e) exc-=(*flow)[e];
2.714 + excess.set(w,exc);
2.715 +
2.716 + //putting the active nodes into the stack
2.717 + int lev=level[w];
2.718 + if ( exc > 0 && lev < n && Node(w) != t ) {
2.719 + next.set(w,first[lev]);
2.720 + first[lev]=w;
2.721 + }
2.722 + }
2.723 + break;
2.724 + } //switch
2.725 + } //preflowPreproc
2.726 +
2.727 +
2.728 + void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
2.729 + VecNode& level_list, NNMap& left,
2.730 + NNMap& right, int& b, int& k, bool what_heur )
2.731 + {
2.732 +
2.733 + int lev=level[w];
2.734 +
2.735 + Node right_n=right[w];
2.736 + Node left_n=left[w];
2.737 +
2.738 + //unlacing starts
2.739 + if ( right_n!=INVALID ) {
2.740 + if ( left_n!=INVALID ) {
2.741 + right.set(left_n, right_n);
2.742 + left.set(right_n, left_n);
2.743 + } else {
2.744 + level_list[lev]=right_n;
2.745 + left.set(right_n, INVALID);
2.746 + }
2.747 + } else {
2.748 + if ( left_n!=INVALID ) {
2.749 + right.set(left_n, INVALID);
2.750 + } else {
2.751 + level_list[lev]=INVALID;
2.752 + }
2.753 + }
2.754 + //unlacing ends
2.755 +
2.756 + if ( level_list[lev]==INVALID ) {
2.757 +
2.758 + //gapping starts
2.759 + for (int i=lev; i!=k ; ) {
2.760 + Node v=level_list[++i];
2.761 + while ( v!=INVALID ) {
2.762 + level.set(v,n);
2.763 + v=right[v];
2.764 + }
2.765 + level_list[i]=INVALID;
2.766 + if ( !what_heur ) first[i]=INVALID;
2.767 + }
2.768 +
2.769 + level.set(w,n);
2.770 + b=lev-1;
2.771 + k=b;
2.772 + //gapping ends
2.773 +
2.774 + } else {
2.775 +
2.776 + if ( newlevel == n ) level.set(w,n);
2.777 + else {
2.778 + level.set(w,++newlevel);
2.779 + next.set(w,first[newlevel]);
2.780 + first[newlevel]=w;
2.781 + if ( what_heur ) b=newlevel;
2.782 + if ( k < newlevel ) ++k; //now k=newlevel
2.783 + Node z=level_list[newlevel];
2.784 + if ( z!=INVALID ) left.set(z,w);
2.785 + right.set(w,z);
2.786 + left.set(w,INVALID);
2.787 + level_list[newlevel]=w;
2.788 + }
2.789 + }
2.790 + } //relabel
2.791 +
2.792 + };
2.793 +} //namespace hugo
2.794 +
2.795 +#endif //HUGO_PREFLOW_H
2.796 +
2.797 +
2.798 +
2.799 +