Communication with algorithm window is developed.
2 * lemon/preflow.h - Part of LEMON, a generic C++ optimization library
4 * Copyright (C) 2006 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/error.h>
24 #include <lemon/invalid.h>
25 #include <lemon/tolerance.h>
26 #include <lemon/maps.h>
27 #include <lemon/graph_utils.h>
31 /// \brief Implementation of the preflow algorithm.
36 ///\brief %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>,
66 typename TOL=Tolerance<Num> >
69 typedef typename Graph::Node Node;
70 typedef typename Graph::NodeIt NodeIt;
71 typedef typename Graph::EdgeIt EdgeIt;
72 typedef typename Graph::OutEdgeIt OutEdgeIt;
73 typedef typename Graph::InEdgeIt InEdgeIt;
75 typedef typename Graph::template NodeMap<Node> NNMap;
76 typedef typename std::vector<Node> VecNode;
81 const CapacityMap* _capacity;
86 int _node_num; //the number of nodes of G
88 typename Graph::template NodeMap<int> level;
89 typename Graph::template NodeMap<Num> excess;
91 // constants used for heuristics
92 static const int H0=20;
93 static const int H1=1;
97 ///\ref Exception for the case when s=t.
99 ///\ref Exception for the case when the source equals the target.
100 class InvalidArgument : public lemon::LogicError {
102 virtual const char* exceptionName() const {
103 return "lemon::Preflow::InvalidArgument";
108 ///Indicates the property of the starting flow map.
110 ///Indicates the property of the starting flow map.
111 ///The meanings are as follows:
112 ///- \c ZERO_FLOW: constant zero flow
113 ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
114 ///the sum of the out-flows in every node except the \e source and
116 ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
117 ///least the sum of the out-flows in every node except the \e source.
118 ///- \c NO_FLOW: indicates an unspecified edge map. \c flow will be
119 ///set to the constant zero flow in the beginning of
120 ///the algorithm in this case.
129 ///Indicates the state of the preflow algorithm.
131 ///Indicates the state of the preflow algorithm.
132 ///The meanings are as follows:
133 ///- \c AFTER_NOTHING: before running the algorithm or
134 /// at an unspecified state.
135 ///- \c AFTER_PREFLOW_PHASE_1: right after running \c phase1
136 ///- \c AFTER_PREFLOW_PHASE_2: after running \ref phase2()
140 AFTER_PREFLOW_PHASE_1,
141 AFTER_PREFLOW_PHASE_2
146 StatusEnum status; // Do not needle this flag only if necessary.
149 ///The constructor of the class.
151 ///The constructor of the class.
152 ///\param _gr The directed graph the algorithm runs on.
153 ///\param _s The source node.
154 ///\param _t The target node.
155 ///\param _cap The capacity of the edges.
156 ///\param _f The flow of the edges.
157 ///Except the graph, all of these parameters can be reset by
158 ///calling \ref source, \ref target, \ref capacityMap and \ref
160 Preflow(const Graph& _gr, Node _s, Node _t,
161 const CapacityMap& _cap, FlowMap& _f,
162 const TOL &tol=TOL()) :
163 _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
164 _flow(&_f), surely(tol),
165 _node_num(countNodes(_gr)), level(_gr), excess(_gr,0),
166 flow_prop(NO_FLOW), status(AFTER_NOTHING) {
167 if ( _source==_target )
168 throw InvalidArgument();
173 ///Runs the preflow algorithm.
175 ///Runs the preflow algorithm.
182 ///Runs the preflow algorithm.
184 ///Runs the preflow algorithm.
185 ///\pre The starting flow map must be
186 /// - a constant zero flow if \c fp is \c ZERO_FLOW,
187 /// - an arbitrary flow if \c fp is \c GEN_FLOW,
188 /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
189 /// - any map if \c fp is NO_FLOW.
190 ///If the starting flow map is a flow or a preflow then
191 ///the algorithm terminates faster.
192 void run(FlowEnum fp) {
197 ///Runs the first phase of the preflow algorithm.
199 ///The preflow algorithm consists of two phases, this method runs
200 ///the first phase. After the first phase the maximum flow value
201 ///and a minimum value cut can already be computed, although a
202 ///maximum flow is not yet obtained. So after calling this method
203 ///\ref flowValue returns the value of a maximum flow and \ref
204 ///minCut returns a minimum cut.
205 ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
206 ///value cuts unless calling \ref phase2.
207 ///\pre The starting flow must be
208 ///- a constant zero flow if \c fp is \c ZERO_FLOW,
209 ///- an arbitary flow if \c fp is \c GEN_FLOW,
210 ///- an arbitary preflow if \c fp is \c PRE_FLOW,
211 ///- any map if \c fp is NO_FLOW.
212 void phase1(FlowEnum fp)
219 ///Runs the first phase of the preflow algorithm.
221 ///The preflow algorithm consists of two phases, this method runs
222 ///the first phase. After the first phase the maximum flow value
223 ///and a minimum value cut can already be computed, although a
224 ///maximum flow is not yet obtained. So after calling this method
225 ///\ref flowValue returns the value of a maximum flow and \ref
226 ///minCut returns a minimum cut.
227 ///\warning \ref minMinCut() and \ref maxMinCut() do not
228 ///give minimum value cuts unless calling \ref phase2().
231 int heur0=(int)(H0*_node_num); //time while running 'bound decrease'
232 int heur1=(int)(H1*_node_num); //time while running 'highest label'
233 int heur=heur1; //starting time interval (#of relabels)
237 //It is 0 in case 'bound decrease' and 1 in case 'highest label'
240 //Needed for 'bound decrease', true means no active
241 //nodes are above bound b.
243 int k=_node_num-2; //bound on the highest level under n containing a node
244 int b=k; //bound on the highest level under n of an active node
246 VecNode first(_node_num, INVALID);
247 NNMap next(*_g, INVALID);
249 NNMap left(*_g, INVALID);
250 NNMap right(*_g, INVALID);
251 VecNode level_list(_node_num,INVALID);
252 //List of the nodes in level i<n, set to n.
254 preflowPreproc(first, next, level_list, left, right);
256 //Push/relabel on the highest level active nodes.
259 if ( !what_heur && !end && k > 0 ) {
265 if ( first[b]==INVALID ) --b;
270 int newlevel=push(w, next, first);
271 if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list,
272 left, right, b, k, what_heur);
275 if ( numrelabel >= heur ) {
290 status=AFTER_PREFLOW_PHASE_1;
295 // list 'level_list' on the nodes on level i implemented by hand
296 // stack 'active' on the active nodes on level i
297 // runs heuristic 'highest label' for H1*n relabels
298 // runs heuristic 'bound decrease' for H0*n relabels,
299 // starts with 'highest label'
300 // Parameters H0 and H1 are initialized to 20 and 1.
303 ///Runs the second phase of the preflow algorithm.
305 ///The preflow algorithm consists of two phases, this method runs
306 ///the second phase. After calling \ref phase1() and then
308 /// \ref flowMap() return a maximum flow, \ref flowValue
309 ///returns the value of a maximum flow, \ref minCut returns a
310 ///minimum cut, while the methods \ref minMinCut and \ref
311 ///maxMinCut return the inclusionwise minimum and maximum cuts of
312 ///minimum value, resp. \pre \ref phase1 must be called before.
316 int k=_node_num-2; //bound on the highest level under n containing a node
317 int b=k; //bound on the highest level under n of an active node
320 VecNode first(_node_num, INVALID);
321 NNMap next(*_g, INVALID);
322 level.set(_source,0);
323 std::queue<Node> bfs_queue;
324 bfs_queue.push(_source);
326 while ( !bfs_queue.empty() ) {
328 Node v=bfs_queue.front();
332 for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
333 if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
334 Node u=_g->source(e);
335 if ( level[u] >= _node_num ) {
338 if ( excess[u] > 0 ) {
339 next.set(u,first[l]);
345 for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
346 if ( 0 >= (*_flow)[e] ) continue;
347 Node u=_g->target(e);
348 if ( level[u] >= _node_num ) {
351 if ( excess[u] > 0 ) {
352 next.set(u,first[l]);
363 if ( first[b]==INVALID ) --b;
367 int newlevel=push(w,next, first);
370 if ( excess[w] > 0 ) {
371 level.set(w,++newlevel);
372 next.set(w,first[newlevel]);
379 status=AFTER_PREFLOW_PHASE_2;
382 /// Returns the value of the maximum flow.
384 /// Returns the value of the maximum flow by returning the excess
385 /// of the target node \c t. This value equals to the value of
386 /// the maximum flow already after running \ref phase1.
387 Num flowValue() const {
388 return excess[_target];
392 ///Returns a minimum value cut.
394 ///Sets \c M to the characteristic vector of a minimum value
395 ///cut. This method can be called both after running \ref
396 ///phase1 and \ref phase2. It is much faster after
397 ///\ref phase1. \pre M should be a bool-valued node-map. \pre
398 ///If \ref minCut() is called after \ref phase2() then M should
399 ///be initialized to false.
400 template<typename _CutMap>
401 void minCut(_CutMap& M) const {
403 case AFTER_PREFLOW_PHASE_1:
404 for(NodeIt v(*_g); v!=INVALID; ++v) {
405 if (level[v] < _node_num) {
412 case AFTER_PREFLOW_PHASE_2:
420 ///Returns the inclusionwise minimum of the minimum value cuts.
422 ///Sets \c M to the characteristic vector of the minimum value cut
423 ///which is inclusionwise minimum. It is computed by processing a
424 ///bfs from the source node \c s in the residual graph. \pre M
425 ///should be a node map of bools initialized to false. \pre \ref
426 ///phase2 should already be run.
427 template<typename _CutMap>
428 void minMinCut(_CutMap& M) const {
430 std::queue<Node> queue;
434 while (!queue.empty()) {
435 Node w=queue.front();
438 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
439 Node v=_g->target(e);
440 if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) {
446 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
447 Node v=_g->source(e);
448 if (!M[v] && (*_flow)[e] > 0 ) {
456 ///Returns the inclusionwise maximum of the minimum value cuts.
458 ///Sets \c M to the characteristic vector of the minimum value cut
459 ///which is inclusionwise maximum. It is computed by processing a
460 ///backward bfs from the target node \c t in the residual graph.
461 ///\pre \ref phase2() or run() should already be run.
462 template<typename _CutMap>
463 void maxMinCut(_CutMap& M) const {
465 for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
467 std::queue<Node> queue;
469 M.set(_target,false);
472 while (!queue.empty()) {
473 Node w=queue.front();
476 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
477 Node v=_g->source(e);
478 if (M[v] && (*_flow)[e] < (*_capacity)[e] ) {
484 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
485 Node v=_g->target(e);
486 if (M[v] && (*_flow)[e] > 0 ) {
494 ///Sets the source node to \c _s.
496 ///Sets the source node to \c _s.
498 void source(Node _s) {
500 if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
501 status=AFTER_NOTHING;
504 ///Returns the source node.
506 ///Returns the source node.
508 Node source() const {
512 ///Sets the target node to \c _t.
514 ///Sets the target node to \c _t.
516 void target(Node _t) {
518 if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
519 status=AFTER_NOTHING;
522 ///Returns the target node.
524 ///Returns the target node.
526 Node target() const {
530 /// Sets the edge map of the capacities to _cap.
532 /// Sets the edge map of the capacities to _cap.
534 void capacityMap(const CapacityMap& _cap) {
536 status=AFTER_NOTHING;
538 /// Returns a reference to capacity map.
540 /// Returns a reference to capacity map.
542 const CapacityMap &capacityMap() const {
546 /// Sets the edge map of the flows to _flow.
548 /// Sets the edge map of the flows to _flow.
550 void flowMap(FlowMap& _f) {
553 status=AFTER_NOTHING;
556 /// Returns a reference to flow map.
558 /// Returns a reference to flow map.
560 const FlowMap &flowMap() const {
566 int push(Node w, NNMap& next, VecNode& first) {
570 int newlevel=_node_num; //bound on the next level of w
572 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
573 if ( (*_flow)[e] >= (*_capacity)[e] ) continue;
574 Node v=_g->target(e);
576 if( lev > level[v] ) { //Push is allowed now
578 if ( excess[v]<=0 && v!=_target && v!=_source ) {
579 next.set(v,first[level[v]]);
583 Num cap=(*_capacity)[e];
587 if ( remcap >= exc ) { //A nonsaturating push.
589 _flow->set(e, flo+exc);
590 excess.set(v, excess[v]+exc);
594 } else { //A saturating push.
596 excess.set(v, excess[v]+remcap);
599 } else if ( newlevel > level[v] ) newlevel = level[v];
603 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
605 if( (*_flow)[e] <= 0 ) continue;
606 Node v=_g->source(e);
608 if( lev > level[v] ) { //Push is allowed now
610 if ( excess[v]<=0 && v!=_target && v!=_source ) {
611 next.set(v,first[level[v]]);
617 if ( flo >= exc ) { //A nonsaturating push.
619 _flow->set(e, flo-exc);
620 excess.set(v, excess[v]+exc);
623 } else { //A saturating push.
625 excess.set(v, excess[v]+flo);
629 } else if ( newlevel > level[v] ) newlevel = level[v];
632 } // if w still has excess after the out edge for cycle
641 void preflowPreproc(VecNode& first, NNMap& next,
642 VecNode& level_list, NNMap& left, NNMap& right)
644 for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
645 std::queue<Node> bfs_queue;
647 if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
648 //Reverse_bfs from t in the residual graph,
649 //to find the starting level.
650 level.set(_target,0);
651 bfs_queue.push(_target);
653 while ( !bfs_queue.empty() ) {
655 Node v=bfs_queue.front();
659 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
660 if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
661 Node w=_g->source(e);
662 if ( level[w] == _node_num && w != _source ) {
664 Node z=level_list[l];
665 if ( z!=INVALID ) left.set(z,w);
672 for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
673 if ( 0 >= (*_flow)[e] ) continue;
674 Node w=_g->target(e);
675 if ( level[w] == _node_num && w != _source ) {
677 Node z=level_list[l];
678 if ( z!=INVALID ) left.set(z,w);
690 for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
692 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
694 //Reverse_bfs from t, to find the starting level.
695 level.set(_target,0);
696 bfs_queue.push(_target);
698 while ( !bfs_queue.empty() ) {
700 Node v=bfs_queue.front();
704 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
705 Node w=_g->source(e);
706 if ( level[w] == _node_num && w != _source ) {
708 Node z=level_list[l];
709 if ( z!=INVALID ) left.set(z,w);
718 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
719 Num c=(*_capacity)[e];
720 if ( c <= 0 ) continue;
721 Node w=_g->target(e);
722 if ( level[w] < _node_num ) {
723 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
724 next.set(w,first[level[w]]);
728 excess.set(w, excess[w]+c);
734 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
737 for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
738 for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
739 excess.set(_target,exc);
743 for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) {
744 Num rem=(*_capacity)[e]-(*_flow)[e];
745 if ( rem <= 0 ) continue;
746 Node w=_g->target(e);
747 if ( level[w] < _node_num ) {
748 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
749 next.set(w,first[level[w]]);
752 _flow->set(e, (*_capacity)[e]);
753 excess.set(w, excess[w]+rem);
757 for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
758 if ( (*_flow)[e] <= 0 ) continue;
759 Node w=_g->source(e);
760 if ( level[w] < _node_num ) {
761 if ( excess[w] <= 0 && w!=_target ) {
762 next.set(w,first[level[w]]);
765 excess.set(w, excess[w]+(*_flow)[e]);
773 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
774 Num rem=(*_capacity)[e]-(*_flow)[e];
775 if ( rem <= 0 ) continue;
776 Node w=_g->target(e);
777 if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
780 for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
781 if ( (*_flow)[e] <= 0 ) continue;
782 Node w=_g->source(e);
783 if ( level[w] < _node_num ) _flow->set(e, 0);
786 //computing the excess
787 for(NodeIt w(*_g); w!=INVALID; ++w) {
789 for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
790 for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
793 //putting the active nodes into the stack
795 if ( exc > 0 && lev < _node_num && Node(w) != _target ) {
796 next.set(w,first[lev]);
805 void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
806 VecNode& level_list, NNMap& left,
807 NNMap& right, int& b, int& k, bool what_heur )
812 Node right_n=right[w];
816 if ( right_n!=INVALID ) {
817 if ( left_n!=INVALID ) {
818 right.set(left_n, right_n);
819 left.set(right_n, left_n);
821 level_list[lev]=right_n;
822 left.set(right_n, INVALID);
825 if ( left_n!=INVALID ) {
826 right.set(left_n, INVALID);
828 level_list[lev]=INVALID;
833 if ( level_list[lev]==INVALID ) {
836 for (int i=lev; i!=k ; ) {
837 Node v=level_list[++i];
838 while ( v!=INVALID ) {
839 level.set(v,_node_num);
842 level_list[i]=INVALID;
843 if ( !what_heur ) first[i]=INVALID;
846 level.set(w,_node_num);
853 if ( newlevel == _node_num ) level.set(w,_node_num);
855 level.set(w,++newlevel);
856 next.set(w,first[newlevel]);
858 if ( what_heur ) b=newlevel;
859 if ( k < newlevel ) ++k; //now k=newlevel
860 Node z=level_list[newlevel];
861 if ( z!=INVALID ) left.set(z,w);
864 level_list[newlevel]=w;
872 ///\brief Function type interface for Preflow algorithm.
874 ///Function type interface for Preflow algorithm.
876 template<class GR, class CM, class FM>
877 Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
878 typename GR::Node source,
879 typename GR::Node target,
884 return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
889 #endif //LEMON_PREFLOW_H