/* -*- C++ -*- * lemon/preflow.h - Part of LEMON, a generic C++ optimization library * * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport * (Egervary Research Group on Combinatorial Optimization, EGRES). * * Permission to use, modify and distribute this software is granted * provided that this copyright notice appears in all copies. For * precise terms see the accompanying LICENSE file. * * This software is provided "AS IS" with no warranty of any kind, * express or implied, and with no claim as to its suitability for any * purpose. * */ #ifndef LEMON_PREFLOW_H #define LEMON_PREFLOW_H #include #include #include #include #include #include /// \file /// \ingroup flowalgs /// \brief Implementation of the preflow algorithm. namespace lemon { ///\ingroup flowalgs ///\brief %Preflow algorithms class. /// ///This class provides an implementation of the \e preflow \e ///algorithm producing a flow of maximum value in a directed ///graph. The preflow algorithms are the fastest known max flow algorithms ///up to now. The \e source node, the \e target node, the \e ///capacity of the edges and the \e starting \e flow value of the ///edges should be passed to the algorithm through the ///constructor. It is possible to change these quantities using the ///functions \ref source, \ref target, \ref capacityMap and \ref ///flowMap. /// ///After running \ref lemon::Preflow::phase1() "phase1()" ///or \ref lemon::Preflow::run() "run()", the maximal flow ///value can be obtained by calling \ref flowValue(). The minimum ///value cut can be written into a bool node map by ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes ///the inclusionwise minimum and maximum of the minimum value cuts, ///resp.) /// ///\param Graph The directed graph type the algorithm runs on. ///\param Num The number type of the capacities and the flow values. ///\param CapacityMap The capacity map type. ///\param FlowMap The flow map type. /// ///\author Jacint Szabo ///\todo Second template parameter is superfluous template , typename FlowMap=typename Graph::template EdgeMap > class Preflow { protected: typedef typename Graph::Node Node; typedef typename Graph::NodeIt NodeIt; typedef typename Graph::EdgeIt EdgeIt; typedef typename Graph::OutEdgeIt OutEdgeIt; typedef typename Graph::InEdgeIt InEdgeIt; typedef typename Graph::template NodeMap NNMap; typedef typename std::vector VecNode; const Graph* _g; Node _source; Node _target; const CapacityMap* _capacity; FlowMap* _flow; int _node_num; //the number of nodes of G typename Graph::template NodeMap level; typename Graph::template NodeMap excess; // constants used for heuristics static const int H0=20; static const int H1=1; public: ///\ref Exception for the case when s=t. ///\ref Exception for the case when the source equals the target. class InvalidArgument : public lemon::LogicError { public: virtual const char* exceptionName() const { return "lemon::Preflow::InvalidArgument"; } }; ///Indicates the property of the starting flow map. ///Indicates the property of the starting flow map. ///The meanings are as follows: ///- \c ZERO_FLOW: constant zero flow ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to ///the sum of the out-flows in every node except the \e source and ///the \e target. ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at ///least the sum of the out-flows in every node except the \e source. ///- \c NO_FLOW: indicates an unspecified edge map. \c flow will be ///set to the constant zero flow in the beginning of ///the algorithm in this case. /// enum FlowEnum{ NO_FLOW, ZERO_FLOW, GEN_FLOW, PRE_FLOW }; ///Indicates the state of the preflow algorithm. ///Indicates the state of the preflow algorithm. ///The meanings are as follows: ///- \c AFTER_NOTHING: before running the algorithm or /// at an unspecified state. ///- \c AFTER_PREFLOW_PHASE_1: right after running \c phase1 ///- \c AFTER_PREFLOW_PHASE_2: after running \ref phase2() /// enum StatusEnum { AFTER_NOTHING, AFTER_PREFLOW_PHASE_1, AFTER_PREFLOW_PHASE_2 }; protected: FlowEnum flow_prop; StatusEnum status; // Do not needle this flag only if necessary. public: ///The constructor of the class. ///The constructor of the class. ///\param _gr The directed graph the algorithm runs on. ///\param _s The source node. ///\param _t The target node. ///\param _cap The capacity of the edges. ///\param _f The flow of the edges. ///Except the graph, all of these parameters can be reset by ///calling \ref source, \ref target, \ref capacityMap and \ref ///flowMap, resp. Preflow(const Graph& _gr, Node _s, Node _t, const CapacityMap& _cap, FlowMap& _f) : _g(&_gr), _source(_s), _target(_t), _capacity(&_cap), _flow(&_f), _node_num(countNodes(_gr)), level(_gr), excess(_gr,0), flow_prop(NO_FLOW), status(AFTER_NOTHING) { if ( _source==_target ) throw InvalidArgument(); } ///Runs the preflow algorithm. ///Runs the preflow algorithm. /// void run() { phase1(flow_prop); phase2(); } ///Runs the preflow algorithm. ///Runs the preflow algorithm. ///\pre The starting flow map must be /// - a constant zero flow if \c fp is \c ZERO_FLOW, /// - an arbitrary flow if \c fp is \c GEN_FLOW, /// - an arbitrary preflow if \c fp is \c PRE_FLOW, /// - any map if \c fp is NO_FLOW. ///If the starting flow map is a flow or a preflow then ///the algorithm terminates faster. void run(FlowEnum fp) { flow_prop=fp; run(); } ///Runs the first phase of the preflow algorithm. ///The preflow algorithm consists of two phases, this method runs ///the first phase. After the first phase the maximum flow value ///and a minimum value cut can already be computed, although a ///maximum flow is not yet obtained. So after calling this method ///\ref flowValue returns the value of a maximum flow and \ref ///minCut returns a minimum cut. ///\warning \ref minMinCut and \ref maxMinCut do not give minimum ///value cuts unless calling \ref phase2. ///\pre The starting flow must be ///- a constant zero flow if \c fp is \c ZERO_FLOW, ///- an arbitary flow if \c fp is \c GEN_FLOW, ///- an arbitary preflow if \c fp is \c PRE_FLOW, ///- any map if \c fp is NO_FLOW. void phase1(FlowEnum fp) { flow_prop=fp; phase1(); } ///Runs the first phase of the preflow algorithm. ///The preflow algorithm consists of two phases, this method runs ///the first phase. After the first phase the maximum flow value ///and a minimum value cut can already be computed, although a ///maximum flow is not yet obtained. So after calling this method ///\ref flowValue returns the value of a maximum flow and \ref ///minCut returns a minimum cut. ///\warning \ref minMinCut() and \ref maxMinCut() do not ///give minimum value cuts unless calling \ref phase2(). void phase1() { int heur0=(int)(H0*_node_num); //time while running 'bound decrease' int heur1=(int)(H1*_node_num); //time while running 'highest label' int heur=heur1; //starting time interval (#of relabels) int numrelabel=0; bool what_heur=1; //It is 0 in case 'bound decrease' and 1 in case 'highest label' bool end=false; //Needed for 'bound decrease', true means no active //nodes are above bound b. int k=_node_num-2; //bound on the highest level under n containing a node int b=k; //bound on the highest level under n of an active node VecNode first(_node_num, INVALID); NNMap next(*_g, INVALID); NNMap left(*_g, INVALID); NNMap right(*_g, INVALID); VecNode level_list(_node_num,INVALID); //List of the nodes in level i 0 ) { b=k; end=true; } else break; } if ( first[b]==INVALID ) --b; else { end=false; Node w=first[b]; first[b]=next[w]; int newlevel=push(w, next, first); if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list, left, right, b, k, what_heur); ++numrelabel; if ( numrelabel >= heur ) { numrelabel=0; if ( what_heur ) { what_heur=0; heur=heur0; end=false; } else { what_heur=1; heur=heur1; b=k; } } } } flow_prop=PRE_FLOW; status=AFTER_PREFLOW_PHASE_1; } // Heuristics: // 2 phase // gap // list 'level_list' on the nodes on level i implemented by hand // stack 'active' on the active nodes on level i // runs heuristic 'highest label' for H1*n relabels // runs heuristic 'bound decrease' for H0*n relabels, // starts with 'highest label' // Parameters H0 and H1 are initialized to 20 and 1. ///Runs the second phase of the preflow algorithm. ///The preflow algorithm consists of two phases, this method runs ///the second phase. After calling \ref phase1() and then ///\ref phase2(), /// \ref flowMap() return a maximum flow, \ref flowValue ///returns the value of a maximum flow, \ref minCut returns a ///minimum cut, while the methods \ref minMinCut and \ref ///maxMinCut return the inclusionwise minimum and maximum cuts of ///minimum value, resp. \pre \ref phase1 must be called before. void phase2() { int k=_node_num-2; //bound on the highest level under n containing a node int b=k; //bound on the highest level under n of an active node VecNode first(_node_num, INVALID); NNMap next(*_g, INVALID); level.set(_source,0); std::queue bfs_queue; bfs_queue.push(_source); while ( !bfs_queue.empty() ) { Node v=bfs_queue.front(); bfs_queue.pop(); int l=level[v]+1; for(InEdgeIt e(*_g,v); e!=INVALID; ++e) { if ( (*_capacity)[e] <= (*_flow)[e] ) continue; Node u=_g->source(e); if ( level[u] >= _node_num ) { bfs_queue.push(u); level.set(u, l); if ( excess[u] > 0 ) { next.set(u,first[l]); first[l]=u; } } } for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) { if ( 0 >= (*_flow)[e] ) continue; Node u=_g->target(e); if ( level[u] >= _node_num ) { bfs_queue.push(u); level.set(u, l); if ( excess[u] > 0 ) { next.set(u,first[l]); first[l]=u; } } } } b=_node_num-2; while ( true ) { if ( b == 0 ) break; if ( first[b]==INVALID ) --b; else { Node w=first[b]; first[b]=next[w]; int newlevel=push(w,next, first); //relabel if ( excess[w] > 0 ) { level.set(w,++newlevel); next.set(w,first[newlevel]); first[newlevel]=w; b=newlevel; } } } // while(true) flow_prop=GEN_FLOW; status=AFTER_PREFLOW_PHASE_2; } /// Returns the value of the maximum flow. /// Returns the value of the maximum flow by returning the excess /// of the target node \c t. This value equals to the value of /// the maximum flow already after running \ref phase1. Num flowValue() const { return excess[_target]; } ///Returns a minimum value cut. ///Sets \c M to the characteristic vector of a minimum value ///cut. This method can be called both after running \ref ///phase1 and \ref phase2. It is much faster after ///\ref phase1. \pre M should be a bool-valued node-map. \pre ///If \ref minCut() is called after \ref phase2() then M should ///be initialized to false. template void minCut(_CutMap& M) const { switch ( status ) { case AFTER_PREFLOW_PHASE_1: for(NodeIt v(*_g); v!=INVALID; ++v) { if (level[v] < _node_num) { M.set(v, false); } else { M.set(v, true); } } break; case AFTER_PREFLOW_PHASE_2: minMinCut(M); break; case AFTER_NOTHING: break; } } ///Returns the inclusionwise minimum of the minimum value cuts. ///Sets \c M to the characteristic vector of the minimum value cut ///which is inclusionwise minimum. It is computed by processing a ///bfs from the source node \c s in the residual graph. \pre M ///should be a node map of bools initialized to false. \pre \ref ///phase2 should already be run. template void minMinCut(_CutMap& M) const { std::queue queue; M.set(_source,true); queue.push(_source); while (!queue.empty()) { Node w=queue.front(); queue.pop(); for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { Node v=_g->target(e); if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) { queue.push(v); M.set(v, true); } } for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { Node v=_g->source(e); if (!M[v] && (*_flow)[e] > 0 ) { queue.push(v); M.set(v, true); } } } } ///Returns the inclusionwise maximum of the minimum value cuts. ///Sets \c M to the characteristic vector of the minimum value cut ///which is inclusionwise maximum. It is computed by processing a ///backward bfs from the target node \c t in the residual graph. ///\pre \ref phase2() or run() should already be run. template void maxMinCut(_CutMap& M) const { for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true); std::queue queue; M.set(_target,false); queue.push(_target); while (!queue.empty()) { Node w=queue.front(); queue.pop(); for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { Node v=_g->source(e); if (M[v] && (*_flow)[e] < (*_capacity)[e] ) { queue.push(v); M.set(v, false); } } for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { Node v=_g->target(e); if (M[v] && (*_flow)[e] > 0 ) { queue.push(v); M.set(v, false); } } } } ///Sets the source node to \c _s. ///Sets the source node to \c _s. /// void source(Node _s) { _source=_s; if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW; status=AFTER_NOTHING; } ///Returns the source node. ///Returns the source node. /// Node source() const { return _source; } ///Sets the target node to \c _t. ///Sets the target node to \c _t. /// void target(Node _t) { _target=_t; if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW; status=AFTER_NOTHING; } ///Returns the target node. ///Returns the target node. /// Node target() const { return _target; } /// Sets the edge map of the capacities to _cap. /// Sets the edge map of the capacities to _cap. /// void capacityMap(const CapacityMap& _cap) { _capacity=&_cap; status=AFTER_NOTHING; } /// Returns a reference to capacity map. /// Returns a reference to capacity map. /// const CapacityMap &capacityMap() const { return *_capacity; } /// Sets the edge map of the flows to _flow. /// Sets the edge map of the flows to _flow. /// void flowMap(FlowMap& _f) { _flow=&_f; flow_prop=NO_FLOW; status=AFTER_NOTHING; } /// Returns a reference to flow map. /// Returns a reference to flow map. /// const FlowMap &flowMap() const { return *_flow; } private: int push(Node w, NNMap& next, VecNode& first) { int lev=level[w]; Num exc=excess[w]; int newlevel=_node_num; //bound on the next level of w for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { if ( (*_flow)[e] >= (*_capacity)[e] ) continue; Node v=_g->target(e); if( lev > level[v] ) { //Push is allowed now if ( excess[v]<=0 && v!=_target && v!=_source ) { next.set(v,first[level[v]]); first[level[v]]=v; } Num cap=(*_capacity)[e]; Num flo=(*_flow)[e]; Num remcap=cap-flo; if ( remcap >= exc ) { //A nonsaturating push. _flow->set(e, flo+exc); excess.set(v, excess[v]+exc); exc=0; break; } else { //A saturating push. _flow->set(e, cap); excess.set(v, excess[v]+remcap); exc-=remcap; } } else if ( newlevel > level[v] ) newlevel = level[v]; } //for out edges wv if ( exc > 0 ) { for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { if( (*_flow)[e] <= 0 ) continue; Node v=_g->source(e); if( lev > level[v] ) { //Push is allowed now if ( excess[v]<=0 && v!=_target && v!=_source ) { next.set(v,first[level[v]]); first[level[v]]=v; } Num flo=(*_flow)[e]; if ( flo >= exc ) { //A nonsaturating push. _flow->set(e, flo-exc); excess.set(v, excess[v]+exc); exc=0; break; } else { //A saturating push. excess.set(v, excess[v]+flo); exc-=flo; _flow->set(e,0); } } else if ( newlevel > level[v] ) newlevel = level[v]; } //for in edges vw } // if w still has excess after the out edge for cycle excess.set(w, exc); return newlevel; } void preflowPreproc(VecNode& first, NNMap& next, VecNode& level_list, NNMap& left, NNMap& right) { for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num); std::queue bfs_queue; if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) { //Reverse_bfs from t in the residual graph, //to find the starting level. level.set(_target,0); bfs_queue.push(_target); while ( !bfs_queue.empty() ) { Node v=bfs_queue.front(); bfs_queue.pop(); int l=level[v]+1; for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) { if ( (*_capacity)[e] <= (*_flow)[e] ) continue; Node w=_g->source(e); if ( level[w] == _node_num && w != _source ) { bfs_queue.push(w); Node z=level_list[l]; if ( z!=INVALID ) left.set(z,w); right.set(w,z); level_list[l]=w; level.set(w, l); } } for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) { if ( 0 >= (*_flow)[e] ) continue; Node w=_g->target(e); if ( level[w] == _node_num && w != _source ) { bfs_queue.push(w); Node z=level_list[l]; if ( z!=INVALID ) left.set(z,w); right.set(w,z); level_list[l]=w; level.set(w, l); } } } //while } //if switch (flow_prop) { case NO_FLOW: for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0); case ZERO_FLOW: for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0); //Reverse_bfs from t, to find the starting level. level.set(_target,0); bfs_queue.push(_target); while ( !bfs_queue.empty() ) { Node v=bfs_queue.front(); bfs_queue.pop(); int l=level[v]+1; for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) { Node w=_g->source(e); if ( level[w] == _node_num && w != _source ) { bfs_queue.push(w); Node z=level_list[l]; if ( z!=INVALID ) left.set(z,w); right.set(w,z); level_list[l]=w; level.set(w, l); } } } //the starting flow for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { Num c=(*_capacity)[e]; if ( c <= 0 ) continue; Node w=_g->target(e); if ( level[w] < _node_num ) { if ( excess[w] <= 0 && w!=_target ) { //putting into the stack next.set(w,first[level[w]]); first[level[w]]=w; } _flow->set(e, c); excess.set(w, excess[w]+c); } } break; case GEN_FLOW: for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0); { Num exc=0; for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e]; for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e]; excess.set(_target,exc); } //the starting flow for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) { Num rem=(*_capacity)[e]-(*_flow)[e]; if ( rem <= 0 ) continue; Node w=_g->target(e); if ( level[w] < _node_num ) { if ( excess[w] <= 0 && w!=_target ) { //putting into the stack next.set(w,first[level[w]]); first[level[w]]=w; } _flow->set(e, (*_capacity)[e]); excess.set(w, excess[w]+rem); } } for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) { if ( (*_flow)[e] <= 0 ) continue; Node w=_g->source(e); if ( level[w] < _node_num ) { if ( excess[w] <= 0 && w!=_target ) { next.set(w,first[level[w]]); first[level[w]]=w; } excess.set(w, excess[w]+(*_flow)[e]); _flow->set(e, 0); } } break; case PRE_FLOW: //the starting flow for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { Num rem=(*_capacity)[e]-(*_flow)[e]; if ( rem <= 0 ) continue; Node w=_g->target(e); if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]); } for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { if ( (*_flow)[e] <= 0 ) continue; Node w=_g->source(e); if ( level[w] < _node_num ) _flow->set(e, 0); } //computing the excess for(NodeIt w(*_g); w!=INVALID; ++w) { Num exc=0; for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e]; for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e]; excess.set(w,exc); //putting the active nodes into the stack int lev=level[w]; if ( exc > 0 && lev < _node_num && Node(w) != _target ) { next.set(w,first[lev]); first[lev]=w; } } break; } //switch } //preflowPreproc void relabel(Node w, int newlevel, VecNode& first, NNMap& next, VecNode& level_list, NNMap& left, NNMap& right, int& b, int& k, bool what_heur ) { int lev=level[w]; Node right_n=right[w]; Node left_n=left[w]; //unlacing starts if ( right_n!=INVALID ) { if ( left_n!=INVALID ) { right.set(left_n, right_n); left.set(right_n, left_n); } else { level_list[lev]=right_n; left.set(right_n, INVALID); } } else { if ( left_n!=INVALID ) { right.set(left_n, INVALID); } else { level_list[lev]=INVALID; } } //unlacing ends if ( level_list[lev]==INVALID ) { //gapping starts for (int i=lev; i!=k ; ) { Node v=level_list[++i]; while ( v!=INVALID ) { level.set(v,_node_num); v=right[v]; } level_list[i]=INVALID; if ( !what_heur ) first[i]=INVALID; } level.set(w,_node_num); b=lev-1; k=b; //gapping ends } else { if ( newlevel == _node_num ) level.set(w,_node_num); else { level.set(w,++newlevel); next.set(w,first[newlevel]); first[newlevel]=w; if ( what_heur ) b=newlevel; if ( k < newlevel ) ++k; //now k=newlevel Node z=level_list[newlevel]; if ( z!=INVALID ) left.set(z,w); right.set(w,z); left.set(w,INVALID); level_list[newlevel]=w; } } } //relabel }; ///\ingroup flowalgs ///\brief Function type interface for Preflow algorithm. /// ///Function type interface for Preflow algorithm. ///\sa Preflow template Preflow preflow(const GR &g, typename GR::Node source, typename GR::Node target, const CM &cap, FM &flow ) { return Preflow(g,source,target,cap,flow); } } //namespace lemon #endif //LEMON_PREFLOW_H