Toolbar has custom icons now. Works only when installed.
2 * lemon/preflow.h - Part of LEMON, a generic C++ optimization library
4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
5 * (Egervary Research Group on Combinatorial Optimization, EGRES).
7 * Permission to use, modify and distribute this software is granted
8 * provided that this copyright notice appears in all copies. For
9 * precise terms see the accompanying LICENSE file.
11 * This software is provided "AS IS" with no warranty of any kind,
12 * express or implied, and with no claim as to its suitability for any
17 #ifndef LEMON_PREFLOW_H
18 #define LEMON_PREFLOW_H
23 #include <lemon/invalid.h>
24 #include <lemon/maps.h>
25 #include <lemon/graph_utils.h>
29 /// Implementation of the preflow algorithm.
33 /// \addtogroup flowalgs
36 ///%Preflow algorithms class.
38 ///This class provides an implementation of the \e preflow \e
39 ///algorithm producing a flow of maximum value in a directed
40 ///graph. The preflow algorithms are the fastest known max flow algorithms
41 ///up to now. The \e source node, the \e target node, the \e
42 ///capacity of the edges and the \e starting \e flow value of the
43 ///edges should be passed to the algorithm through the
44 ///constructor. It is possible to change these quantities using the
45 ///functions \ref source, \ref target, \ref capacityMap and \ref
48 ///After running \ref lemon::Preflow::phase1() "phase1()"
49 ///or \ref lemon::Preflow::run() "run()", the maximal flow
50 ///value can be obtained by calling \ref flowValue(). The minimum
51 ///value cut can be written into a <tt>bool</tt> node map by
52 ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes
53 ///the inclusionwise minimum and maximum of the minimum value cuts,
56 ///\param Graph The directed graph type the algorithm runs on.
57 ///\param Num The number type of the capacities and the flow values.
58 ///\param CapacityMap The capacity map type.
59 ///\param FlowMap The flow map type.
61 ///\author Jacint Szabo
62 ///\todo Second template parameter is superfluous
63 template <typename Graph, typename Num,
64 typename CapacityMap=typename Graph::template EdgeMap<Num>,
65 typename FlowMap=typename Graph::template EdgeMap<Num> >
68 typedef typename Graph::Node Node;
69 typedef typename Graph::NodeIt NodeIt;
70 typedef typename Graph::EdgeIt EdgeIt;
71 typedef typename Graph::OutEdgeIt OutEdgeIt;
72 typedef typename Graph::InEdgeIt InEdgeIt;
74 typedef typename Graph::template NodeMap<Node> NNMap;
75 typedef typename std::vector<Node> VecNode;
80 const CapacityMap* _capacity;
82 int _node_num; //the number of nodes of G
84 typename Graph::template NodeMap<int> level;
85 typename Graph::template NodeMap<Num> excess;
87 // constants used for heuristics
88 static const int H0=20;
89 static const int H1=1;
93 ///Indicates the property of the starting flow map.
95 ///Indicates the property of the starting flow map.
96 ///The meanings are as follows:
97 ///- \c ZERO_FLOW: constant zero flow
98 ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
99 ///the sum of the out-flows in every node except the \e source and
101 ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
102 ///least the sum of the out-flows in every node except the \e source.
103 ///- \c NO_FLOW: indicates an unspecified edge map. \c flow will be
104 ///set to the constant zero flow in the beginning of
105 ///the algorithm in this case.
114 ///Indicates the state of the preflow algorithm.
116 ///Indicates the state of the preflow algorithm.
117 ///The meanings are as follows:
118 ///- \c AFTER_NOTHING: before running the algorithm or
119 /// at an unspecified state.
120 ///- \c AFTER_PREFLOW_PHASE_1: right after running \c phase1
121 ///- \c AFTER_PREFLOW_PHASE_2: after running \ref phase2()
125 AFTER_PREFLOW_PHASE_1,
126 AFTER_PREFLOW_PHASE_2
131 StatusEnum status; // Do not needle this flag only if necessary.
134 ///The constructor of the class.
136 ///The constructor of the class.
137 ///\param _gr The directed graph the algorithm runs on.
138 ///\param _s The source node.
139 ///\param _t The target node.
140 ///\param _cap The capacity of the edges.
141 ///\param _f The flow of the edges.
142 ///Except the graph, all of these parameters can be reset by
143 ///calling \ref source, \ref target, \ref capacityMap and \ref
145 Preflow(const Graph& _gr, Node _s, Node _t,
146 const CapacityMap& _cap, FlowMap& _f) :
147 _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
148 _flow(&_f), _node_num(countNodes(_gr)), level(_gr), excess(_gr,0),
149 flow_prop(NO_FLOW), status(AFTER_NOTHING) { }
153 ///Runs the preflow algorithm.
155 ///Runs the preflow algorithm.
162 ///Runs the preflow algorithm.
164 ///Runs the preflow algorithm.
165 ///\pre The starting flow map must be
166 /// - a constant zero flow if \c fp is \c ZERO_FLOW,
167 /// - an arbitrary flow if \c fp is \c GEN_FLOW,
168 /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
169 /// - any map if \c fp is NO_FLOW.
170 ///If the starting flow map is a flow or a preflow then
171 ///the algorithm terminates faster.
172 void run(FlowEnum fp) {
177 ///Runs the first phase of the preflow algorithm.
179 ///The preflow algorithm consists of two phases, this method runs
180 ///the first phase. After the first phase the maximum flow value
181 ///and a minimum value cut can already be computed, although a
182 ///maximum flow is not yet obtained. So after calling this method
183 ///\ref flowValue returns the value of a maximum flow and \ref
184 ///minCut returns a minimum cut.
185 ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
186 ///value cuts unless calling \ref phase2.
187 ///\pre The starting flow must be
188 ///- a constant zero flow if \c fp is \c ZERO_FLOW,
189 ///- an arbitary flow if \c fp is \c GEN_FLOW,
190 ///- an arbitary preflow if \c fp is \c PRE_FLOW,
191 ///- any map if \c fp is NO_FLOW.
192 void phase1(FlowEnum fp)
199 ///Runs the first phase of the preflow algorithm.
201 ///The preflow algorithm consists of two phases, this method runs
202 ///the first phase. After the first phase the maximum flow value
203 ///and a minimum value cut can already be computed, although a
204 ///maximum flow is not yet obtained. So after calling this method
205 ///\ref flowValue returns the value of a maximum flow and \ref
206 ///minCut returns a minimum cut.
207 ///\warning \ref minCut(), \ref minMinCut() and \ref maxMinCut() do not
208 ///give minimum value cuts unless calling \ref phase2().
211 int heur0=(int)(H0*_node_num); //time while running 'bound decrease'
212 int heur1=(int)(H1*_node_num); //time while running 'highest label'
213 int heur=heur1; //starting time interval (#of relabels)
217 //It is 0 in case 'bound decrease' and 1 in case 'highest label'
220 //Needed for 'bound decrease', true means no active
221 //nodes are above bound b.
223 int k=_node_num-2; //bound on the highest level under n containing a node
224 int b=k; //bound on the highest level under n of an active node
226 VecNode first(_node_num, INVALID);
227 NNMap next(*_g, INVALID);
229 NNMap left(*_g, INVALID);
230 NNMap right(*_g, INVALID);
231 VecNode level_list(_node_num,INVALID);
232 //List of the nodes in level i<n, set to n.
234 preflowPreproc(first, next, level_list, left, right);
236 //Push/relabel on the highest level active nodes.
239 if ( !what_heur && !end && k > 0 ) {
245 if ( first[b]==INVALID ) --b;
250 int newlevel=push(w, next, first);
251 if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list,
252 left, right, b, k, what_heur);
255 if ( numrelabel >= heur ) {
270 status=AFTER_PREFLOW_PHASE_1;
275 // list 'level_list' on the nodes on level i implemented by hand
276 // stack 'active' on the active nodes on level i
277 // runs heuristic 'highest label' for H1*n relabels
278 // runs heuristic 'bound decrease' for H0*n relabels,
279 // starts with 'highest label'
280 // Parameters H0 and H1 are initialized to 20 and 1.
283 ///Runs the second phase of the preflow algorithm.
285 ///The preflow algorithm consists of two phases, this method runs
286 ///the second phase. After calling \ref phase1() and then
288 /// \ref flowMap() return a maximum flow, \ref flowValue
289 ///returns the value of a maximum flow, \ref minCut returns a
290 ///minimum cut, while the methods \ref minMinCut and \ref
291 ///maxMinCut return the inclusionwise minimum and maximum cuts of
292 ///minimum value, resp. \pre \ref phase1 must be called before.
296 int k=_node_num-2; //bound on the highest level under n containing a node
297 int b=k; //bound on the highest level under n of an active node
300 VecNode first(_node_num, INVALID);
301 NNMap next(*_g, INVALID);
302 level.set(_source,0);
303 std::queue<Node> bfs_queue;
304 bfs_queue.push(_source);
306 while ( !bfs_queue.empty() ) {
308 Node v=bfs_queue.front();
312 for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
313 if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
314 Node u=_g->source(e);
315 if ( level[u] >= _node_num ) {
318 if ( excess[u] > 0 ) {
319 next.set(u,first[l]);
325 for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
326 if ( 0 >= (*_flow)[e] ) continue;
327 Node u=_g->target(e);
328 if ( level[u] >= _node_num ) {
331 if ( excess[u] > 0 ) {
332 next.set(u,first[l]);
343 if ( first[b]==INVALID ) --b;
347 int newlevel=push(w,next, first);
350 if ( excess[w] > 0 ) {
351 level.set(w,++newlevel);
352 next.set(w,first[newlevel]);
359 status=AFTER_PREFLOW_PHASE_2;
362 /// Returns the value of the maximum flow.
364 /// Returns the value of the maximum flow by returning the excess
365 /// of the target node \c t. This value equals to the value of
366 /// the maximum flow already after running \ref phase1.
367 Num flowValue() const {
368 return excess[_target];
372 ///Returns a minimum value cut.
374 ///Sets \c M to the characteristic vector of a minimum value
375 ///cut. This method can be called both after running \ref
376 ///phase1 and \ref phase2. It is much faster after
377 ///\ref phase1. \pre M should be a bool-valued node-map. \pre
378 ///If \ref minCut() is called after \ref phase2() then M should
379 ///be initialized to false.
380 template<typename _CutMap>
381 void minCut(_CutMap& M) const {
383 case AFTER_PREFLOW_PHASE_1:
384 for(NodeIt v(*_g); v!=INVALID; ++v) {
385 if (level[v] < _node_num) {
392 case AFTER_PREFLOW_PHASE_2:
400 ///Returns the inclusionwise minimum of the minimum value cuts.
402 ///Sets \c M to the characteristic vector of the minimum value cut
403 ///which is inclusionwise minimum. It is computed by processing a
404 ///bfs from the source node \c s in the residual graph. \pre M
405 ///should be a node map of bools initialized to false. \pre \ref
406 ///phase2 should already be run.
407 template<typename _CutMap>
408 void minMinCut(_CutMap& M) const {
410 std::queue<Node> queue;
414 while (!queue.empty()) {
415 Node w=queue.front();
418 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
419 Node v=_g->target(e);
420 if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) {
426 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
427 Node v=_g->source(e);
428 if (!M[v] && (*_flow)[e] > 0 ) {
436 ///Returns the inclusionwise maximum of the minimum value cuts.
438 ///Sets \c M to the characteristic vector of the minimum value cut
439 ///which is inclusionwise maximum. It is computed by processing a
440 ///backward bfs from the target node \c t in the residual graph.
441 ///\pre \ref phase2() or run() should already be run.
442 template<typename _CutMap>
443 void maxMinCut(_CutMap& M) const {
445 for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
447 std::queue<Node> queue;
449 M.set(_target,false);
452 while (!queue.empty()) {
453 Node w=queue.front();
456 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
457 Node v=_g->source(e);
458 if (M[v] && (*_flow)[e] < (*_capacity)[e] ) {
464 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
465 Node v=_g->target(e);
466 if (M[v] && (*_flow)[e] > 0 ) {
474 ///Sets the source node to \c _s.
476 ///Sets the source node to \c _s.
478 void source(Node _s) {
480 if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
481 status=AFTER_NOTHING;
484 ///Returns the source node.
486 ///Returns the source node.
488 Node source() const {
492 ///Sets the target node to \c _t.
494 ///Sets the target node to \c _t.
496 void target(Node _t) {
498 if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
499 status=AFTER_NOTHING;
502 ///Returns the target node.
504 ///Returns the target node.
506 Node target() const {
510 /// Sets the edge map of the capacities to _cap.
512 /// Sets the edge map of the capacities to _cap.
514 void capacityMap(const CapacityMap& _cap) {
516 status=AFTER_NOTHING;
518 /// Returns a reference to capacity map.
520 /// Returns a reference to capacity map.
522 const CapacityMap &capacityMap() const {
526 /// Sets the edge map of the flows to _flow.
528 /// Sets the edge map of the flows to _flow.
530 void flowMap(FlowMap& _f) {
533 status=AFTER_NOTHING;
536 /// Returns a reference to flow map.
538 /// Returns a reference to flow map.
540 const FlowMap &flowMap() const {
546 int push(Node w, NNMap& next, VecNode& first) {
550 int newlevel=_node_num; //bound on the next level of w
552 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
553 if ( (*_flow)[e] >= (*_capacity)[e] ) continue;
554 Node v=_g->target(e);
556 if( lev > level[v] ) { //Push is allowed now
558 if ( excess[v]<=0 && v!=_target && v!=_source ) {
559 next.set(v,first[level[v]]);
563 Num cap=(*_capacity)[e];
567 if ( remcap >= exc ) { //A nonsaturating push.
569 _flow->set(e, flo+exc);
570 excess.set(v, excess[v]+exc);
574 } else { //A saturating push.
576 excess.set(v, excess[v]+remcap);
579 } else if ( newlevel > level[v] ) newlevel = level[v];
583 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
585 if( (*_flow)[e] <= 0 ) continue;
586 Node v=_g->source(e);
588 if( lev > level[v] ) { //Push is allowed now
590 if ( excess[v]<=0 && v!=_target && v!=_source ) {
591 next.set(v,first[level[v]]);
597 if ( flo >= exc ) { //A nonsaturating push.
599 _flow->set(e, flo-exc);
600 excess.set(v, excess[v]+exc);
603 } else { //A saturating push.
605 excess.set(v, excess[v]+flo);
609 } else if ( newlevel > level[v] ) newlevel = level[v];
612 } // if w still has excess after the out edge for cycle
621 void preflowPreproc(VecNode& first, NNMap& next,
622 VecNode& level_list, NNMap& left, NNMap& right)
624 for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
625 std::queue<Node> bfs_queue;
627 if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
628 //Reverse_bfs from t in the residual graph,
629 //to find the starting level.
630 level.set(_target,0);
631 bfs_queue.push(_target);
633 while ( !bfs_queue.empty() ) {
635 Node v=bfs_queue.front();
639 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
640 if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
641 Node w=_g->source(e);
642 if ( level[w] == _node_num && w != _source ) {
644 Node z=level_list[l];
645 if ( z!=INVALID ) left.set(z,w);
652 for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
653 if ( 0 >= (*_flow)[e] ) continue;
654 Node w=_g->target(e);
655 if ( level[w] == _node_num && w != _source ) {
657 Node z=level_list[l];
658 if ( z!=INVALID ) left.set(z,w);
670 for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
672 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
674 //Reverse_bfs from t, to find the starting level.
675 level.set(_target,0);
676 bfs_queue.push(_target);
678 while ( !bfs_queue.empty() ) {
680 Node v=bfs_queue.front();
684 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
685 Node w=_g->source(e);
686 if ( level[w] == _node_num && w != _source ) {
688 Node z=level_list[l];
689 if ( z!=INVALID ) left.set(z,w);
698 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
699 Num c=(*_capacity)[e];
700 if ( c <= 0 ) continue;
701 Node w=_g->target(e);
702 if ( level[w] < _node_num ) {
703 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
704 next.set(w,first[level[w]]);
708 excess.set(w, excess[w]+c);
714 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
717 for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
718 for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
719 excess.set(_target,exc);
723 for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) {
724 Num rem=(*_capacity)[e]-(*_flow)[e];
725 if ( rem <= 0 ) continue;
726 Node w=_g->target(e);
727 if ( level[w] < _node_num ) {
728 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
729 next.set(w,first[level[w]]);
732 _flow->set(e, (*_capacity)[e]);
733 excess.set(w, excess[w]+rem);
737 for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
738 if ( (*_flow)[e] <= 0 ) continue;
739 Node w=_g->source(e);
740 if ( level[w] < _node_num ) {
741 if ( excess[w] <= 0 && w!=_target ) {
742 next.set(w,first[level[w]]);
745 excess.set(w, excess[w]+(*_flow)[e]);
753 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
754 Num rem=(*_capacity)[e]-(*_flow)[e];
755 if ( rem <= 0 ) continue;
756 Node w=_g->target(e);
757 if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
760 for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
761 if ( (*_flow)[e] <= 0 ) continue;
762 Node w=_g->source(e);
763 if ( level[w] < _node_num ) _flow->set(e, 0);
766 //computing the excess
767 for(NodeIt w(*_g); w!=INVALID; ++w) {
769 for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
770 for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
773 //putting the active nodes into the stack
775 if ( exc > 0 && lev < _node_num && Node(w) != _target ) {
776 next.set(w,first[lev]);
785 void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
786 VecNode& level_list, NNMap& left,
787 NNMap& right, int& b, int& k, bool what_heur )
792 Node right_n=right[w];
796 if ( right_n!=INVALID ) {
797 if ( left_n!=INVALID ) {
798 right.set(left_n, right_n);
799 left.set(right_n, left_n);
801 level_list[lev]=right_n;
802 left.set(right_n, INVALID);
805 if ( left_n!=INVALID ) {
806 right.set(left_n, INVALID);
808 level_list[lev]=INVALID;
813 if ( level_list[lev]==INVALID ) {
816 for (int i=lev; i!=k ; ) {
817 Node v=level_list[++i];
818 while ( v!=INVALID ) {
819 level.set(v,_node_num);
822 level_list[i]=INVALID;
823 if ( !what_heur ) first[i]=INVALID;
826 level.set(w,_node_num);
833 if ( newlevel == _node_num ) level.set(w,_node_num);
835 level.set(w,++newlevel);
836 next.set(w,first[newlevel]);
838 if ( what_heur ) b=newlevel;
839 if ( k < newlevel ) ++k; //now k=newlevel
840 Node z=level_list[newlevel];
841 if ( z!=INVALID ) left.set(z,w);
844 level_list[newlevel]=w;
851 ///Function type interface for Preflow algorithm.
853 /// \ingroup flowalgs
854 ///Function type interface for Preflow algorithm.
856 template<class GR, class CM, class FM>
857 Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
858 typename GR::Node source,
859 typename GR::Node target,
864 return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
869 #endif //LEMON_PREFLOW_H