The graph adadptors can be alteration observed.
In most cases it uses the adapted graph alteration notifiers.
Only special case is now the UndirGraphAdaptor, where
we have to proxy the signals from the graph.
The SubBidirGraphAdaptor is removed, because it doest not
gives more feature than the EdgeSubGraphAdaptor<UndirGraphAdaptor<Graph>>.
The ResGraphAdaptor is based on this composition.
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2006
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
19 #ifndef LEMON_PREFLOW_H
20 #define LEMON_PREFLOW_H
25 #include <lemon/error.h>
26 #include <lemon/invalid.h>
27 #include <lemon/tolerance.h>
28 #include <lemon/maps.h>
29 #include <lemon/graph_utils.h>
33 /// \brief Implementation of the preflow algorithm.
38 ///\brief %Preflow algorithms class.
40 ///This class provides an implementation of the \e preflow \e
41 ///algorithm producing a flow of maximum value in a directed
42 ///graph. The preflow algorithms are the fastest known max flow algorithms
43 ///up to now. The \e source node, the \e target node, the \e
44 ///capacity of the edges and the \e starting \e flow value of the
45 ///edges should be passed to the algorithm through the
46 ///constructor. It is possible to change these quantities using the
47 ///functions \ref source, \ref target, \ref capacityMap and \ref
50 ///After running \ref lemon::Preflow::phase1() "phase1()"
51 ///or \ref lemon::Preflow::run() "run()", the maximal flow
52 ///value can be obtained by calling \ref flowValue(). The minimum
53 ///value cut can be written into a <tt>bool</tt> node map by
54 ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes
55 ///the inclusionwise minimum and maximum of the minimum value cuts,
58 ///\param Graph The directed graph type the algorithm runs on.
59 ///\param Num The number type of the capacities and the flow values.
60 ///\param CapacityMap The capacity map type.
61 ///\param FlowMap The flow map type.
63 ///\author Jacint Szabo
64 ///\todo Second template parameter is superfluous
65 template <typename Graph, typename Num,
66 typename CapacityMap=typename Graph::template EdgeMap<Num>,
67 typename FlowMap=typename Graph::template EdgeMap<Num>,
68 typename TOL=Tolerance<Num> >
71 typedef typename Graph::Node Node;
72 typedef typename Graph::NodeIt NodeIt;
73 typedef typename Graph::EdgeIt EdgeIt;
74 typedef typename Graph::OutEdgeIt OutEdgeIt;
75 typedef typename Graph::InEdgeIt InEdgeIt;
77 typedef typename Graph::template NodeMap<Node> NNMap;
78 typedef typename std::vector<Node> VecNode;
83 const CapacityMap* _capacity;
88 int _node_num; //the number of nodes of G
90 typename Graph::template NodeMap<int> level;
91 typename Graph::template NodeMap<Num> excess;
93 // constants used for heuristics
94 static const int H0=20;
95 static const int H1=1;
99 ///\ref Exception for the case when s=t.
101 ///\ref Exception for the case when the source equals the target.
102 class InvalidArgument : public lemon::LogicError {
104 virtual const char* exceptionName() const {
105 return "lemon::Preflow::InvalidArgument";
110 ///Indicates the property of the starting flow map.
112 ///Indicates the property of the starting flow map.
115 ///indicates an unspecified edge map. \c flow will be
116 ///set to the constant zero flow in the beginning of
117 ///the algorithm in this case.
119 ///constant zero flow
121 ///any flow, i.e. the sum of the in-flows equals to
122 ///the sum of the out-flows in every node except the \c source and
125 ///any preflow, i.e. the sum of the in-flows is at
126 ///least the sum of the out-flows in every node except the \c source.
130 ///Indicates the state of the preflow algorithm.
132 ///Indicates the state of the preflow algorithm.
135 ///before running the algorithm or
136 ///at an unspecified state.
138 ///right after running \ref phase1()
139 AFTER_PREFLOW_PHASE_1,
140 ///after running \ref phase2()
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 ///\param tol Tolerance class.
158 ///Except the graph, all of these parameters can be reset by
159 ///calling \ref source, \ref target, \ref capacityMap and \ref
161 Preflow(const Graph& _gr, Node _s, Node _t,
162 const CapacityMap& _cap, FlowMap& _f,
163 const TOL &tol=TOL()) :
164 _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
165 _flow(&_f), surely(tol),
166 _node_num(countNodes(_gr)), level(_gr), excess(_gr,0),
167 flow_prop(NO_FLOW), status(AFTER_NOTHING) {
168 if ( _source==_target )
169 throw InvalidArgument();
172 ///Give a reference to the tolerance handler class
174 ///Give a reference to the tolerance handler class
176 TOL &tolerance() { return surely; }
178 ///Runs the preflow algorithm.
180 ///Runs the preflow algorithm.
187 ///Runs the preflow algorithm.
189 ///Runs the preflow algorithm.
190 ///\pre The starting flow map must be
191 /// - a constant zero flow if \c fp is \c ZERO_FLOW,
192 /// - an arbitrary flow if \c fp is \c GEN_FLOW,
193 /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
194 /// - any map if \c fp is NO_FLOW.
195 ///If the starting flow map is a flow or a preflow then
196 ///the algorithm terminates faster.
197 void run(FlowEnum fp) {
202 ///Runs the first phase of the preflow algorithm.
204 ///The preflow algorithm consists of two phases, this method runs
205 ///the first phase. After the first phase the maximum flow value
206 ///and a minimum value cut can already be computed, although a
207 ///maximum flow is not yet obtained. So after calling this method
208 ///\ref flowValue returns the value of a maximum flow and \ref
209 ///minCut returns a minimum cut.
210 ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
211 ///value cuts unless calling \ref phase2.
212 ///\pre The starting flow must be
213 ///- a constant zero flow if \c fp is \c ZERO_FLOW,
214 ///- an arbitary flow if \c fp is \c GEN_FLOW,
215 ///- an arbitary preflow if \c fp is \c PRE_FLOW,
216 ///- any map if \c fp is NO_FLOW.
217 void phase1(FlowEnum fp)
224 ///Runs the first phase of the preflow algorithm.
226 ///The preflow algorithm consists of two phases, this method runs
227 ///the first phase. After the first phase the maximum flow value
228 ///and a minimum value cut can already be computed, although a
229 ///maximum flow is not yet obtained. So after calling this method
230 ///\ref flowValue returns the value of a maximum flow and \ref
231 ///minCut returns a minimum cut.
232 ///\warning \ref minMinCut() and \ref maxMinCut() do not
233 ///give minimum value cuts unless calling \ref phase2().
236 int heur0=(int)(H0*_node_num); //time while running 'bound decrease'
237 int heur1=(int)(H1*_node_num); //time while running 'highest label'
238 int heur=heur1; //starting time interval (#of relabels)
242 //It is 0 in case 'bound decrease' and 1 in case 'highest label'
245 //Needed for 'bound decrease', true means no active
246 //nodes are above bound b.
248 int k=_node_num-2; //bound on the highest level under n containing a node
249 int b=k; //bound on the highest level under n of an active node
251 VecNode first(_node_num, INVALID);
252 NNMap next(*_g, INVALID);
254 NNMap left(*_g, INVALID);
255 NNMap right(*_g, INVALID);
256 VecNode level_list(_node_num,INVALID);
257 //List of the nodes in level i<n, set to n.
259 preflowPreproc(first, next, level_list, left, right);
261 //Push/relabel on the highest level active nodes.
264 if ( !what_heur && !end && k > 0 ) {
270 if ( first[b]==INVALID ) --b;
275 int newlevel=push(w, next, first);
276 if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list,
277 left, right, b, k, what_heur);
280 if ( numrelabel >= heur ) {
295 status=AFTER_PREFLOW_PHASE_1;
300 // list 'level_list' on the nodes on level i implemented by hand
301 // stack 'active' on the active nodes on level i
302 // runs heuristic 'highest label' for H1*n relabels
303 // runs heuristic 'bound decrease' for H0*n relabels,
304 // starts with 'highest label'
305 // Parameters H0 and H1 are initialized to 20 and 1.
308 ///Runs the second phase of the preflow algorithm.
310 ///The preflow algorithm consists of two phases, this method runs
311 ///the second phase. After calling \ref phase1() and then
313 /// \ref flowMap() return a maximum flow, \ref flowValue
314 ///returns the value of a maximum flow, \ref minCut returns a
315 ///minimum cut, while the methods \ref minMinCut and \ref
316 ///maxMinCut return the inclusionwise minimum and maximum cuts of
317 ///minimum value, resp. \pre \ref phase1 must be called before.
321 int k=_node_num-2; //bound on the highest level under n containing a node
322 int b=k; //bound on the highest level under n of an active node
325 VecNode first(_node_num, INVALID);
326 NNMap next(*_g, INVALID);
327 level.set(_source,0);
328 std::queue<Node> bfs_queue;
329 bfs_queue.push(_source);
331 while ( !bfs_queue.empty() ) {
333 Node v=bfs_queue.front();
337 for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
338 if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
339 Node u=_g->source(e);
340 if ( level[u] >= _node_num ) {
343 if ( excess[u] > 0 ) {
344 next.set(u,first[l]);
350 for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
351 if ( 0 >= (*_flow)[e] ) continue;
352 Node u=_g->target(e);
353 if ( level[u] >= _node_num ) {
356 if ( excess[u] > 0 ) {
357 next.set(u,first[l]);
368 if ( first[b]==INVALID ) --b;
372 int newlevel=push(w,next, first);
375 if ( excess[w] > 0 ) {
376 level.set(w,++newlevel);
377 next.set(w,first[newlevel]);
384 status=AFTER_PREFLOW_PHASE_2;
387 /// Returns the value of the maximum flow.
389 /// Returns the value of the maximum flow by returning the excess
390 /// of the target node \c t. This value equals to the value of
391 /// the maximum flow already after running \ref phase1.
392 Num flowValue() const {
393 return excess[_target];
397 ///Returns a minimum value cut.
399 ///Sets \c M to the characteristic vector of a minimum value
400 ///cut. This method can be called both after running \ref
401 ///phase1 and \ref phase2. It is much faster after
402 ///\ref phase1. \pre M should be a bool-valued node-map. \pre
403 ///If \ref minCut() is called after \ref phase2() then M should
404 ///be initialized to false.
405 template<typename _CutMap>
406 void minCut(_CutMap& M) const {
408 case AFTER_PREFLOW_PHASE_1:
409 for(NodeIt v(*_g); v!=INVALID; ++v) {
410 if (level[v] < _node_num) {
417 case AFTER_PREFLOW_PHASE_2:
425 ///Returns the inclusionwise minimum of the minimum value cuts.
427 ///Sets \c M to the characteristic vector of the minimum value cut
428 ///which is inclusionwise minimum. It is computed by processing a
429 ///bfs from the source node \c s in the residual graph. \pre M
430 ///should be a node map of bools initialized to false. \pre \ref
431 ///phase2 should already be run.
432 template<typename _CutMap>
433 void minMinCut(_CutMap& M) const {
435 std::queue<Node> queue;
439 while (!queue.empty()) {
440 Node w=queue.front();
443 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
444 Node v=_g->target(e);
445 if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) {
451 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
452 Node v=_g->source(e);
453 if (!M[v] && (*_flow)[e] > 0 ) {
461 ///Returns the inclusionwise maximum of the minimum value cuts.
463 ///Sets \c M to the characteristic vector of the minimum value cut
464 ///which is inclusionwise maximum. It is computed by processing a
465 ///backward bfs from the target node \c t in the residual graph.
466 ///\pre \ref phase2() or run() should already be run.
467 template<typename _CutMap>
468 void maxMinCut(_CutMap& M) const {
470 for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
472 std::queue<Node> queue;
474 M.set(_target,false);
477 while (!queue.empty()) {
478 Node w=queue.front();
481 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
482 Node v=_g->source(e);
483 if (M[v] && (*_flow)[e] < (*_capacity)[e] ) {
489 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
490 Node v=_g->target(e);
491 if (M[v] && (*_flow)[e] > 0 ) {
499 ///Sets the source node to \c _s.
501 ///Sets the source node to \c _s.
503 void source(Node _s) {
505 if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
506 status=AFTER_NOTHING;
509 ///Returns the source node.
511 ///Returns the source node.
513 Node source() const {
517 ///Sets the target node to \c _t.
519 ///Sets the target node to \c _t.
521 void target(Node _t) {
523 if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
524 status=AFTER_NOTHING;
527 ///Returns the target node.
529 ///Returns the target node.
531 Node target() const {
535 /// Sets the edge map of the capacities to _cap.
537 /// Sets the edge map of the capacities to _cap.
539 void capacityMap(const CapacityMap& _cap) {
541 status=AFTER_NOTHING;
543 /// Returns a reference to capacity map.
545 /// Returns a reference to capacity map.
547 const CapacityMap &capacityMap() const {
551 /// Sets the edge map of the flows to _flow.
553 /// Sets the edge map of the flows to _flow.
555 void flowMap(FlowMap& _f) {
558 status=AFTER_NOTHING;
561 /// Returns a reference to flow map.
563 /// Returns a reference to flow map.
565 const FlowMap &flowMap() const {
571 int push(Node w, NNMap& next, VecNode& first) {
575 int newlevel=_node_num; //bound on the next level of w
577 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
578 if ( (*_flow)[e] >= (*_capacity)[e] ) continue;
579 Node v=_g->target(e);
581 if( lev > level[v] ) { //Push is allowed now
583 if ( excess[v]<=0 && v!=_target && v!=_source ) {
584 next.set(v,first[level[v]]);
588 Num cap=(*_capacity)[e];
592 if ( remcap >= exc ) { //A nonsaturating push.
594 _flow->set(e, flo+exc);
595 excess.set(v, excess[v]+exc);
599 } else { //A saturating push.
601 excess.set(v, excess[v]+remcap);
604 } else if ( newlevel > level[v] ) newlevel = level[v];
608 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
610 if( (*_flow)[e] <= 0 ) continue;
611 Node v=_g->source(e);
613 if( lev > level[v] ) { //Push is allowed now
615 if ( excess[v]<=0 && v!=_target && v!=_source ) {
616 next.set(v,first[level[v]]);
622 if ( flo >= exc ) { //A nonsaturating push.
624 _flow->set(e, flo-exc);
625 excess.set(v, excess[v]+exc);
628 } else { //A saturating push.
630 excess.set(v, excess[v]+flo);
634 } else if ( newlevel > level[v] ) newlevel = level[v];
637 } // if w still has excess after the out edge for cycle
646 void preflowPreproc(VecNode& first, NNMap& next,
647 VecNode& level_list, NNMap& left, NNMap& right)
649 for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
650 std::queue<Node> bfs_queue;
652 if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
653 //Reverse_bfs from t in the residual graph,
654 //to find the starting level.
655 level.set(_target,0);
656 bfs_queue.push(_target);
658 while ( !bfs_queue.empty() ) {
660 Node v=bfs_queue.front();
664 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
665 if ( (*_capacity)[e] <= (*_flow)[e] ) continue;
666 Node w=_g->source(e);
667 if ( level[w] == _node_num && w != _source ) {
669 Node z=level_list[l];
670 if ( z!=INVALID ) left.set(z,w);
677 for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
678 if ( 0 >= (*_flow)[e] ) continue;
679 Node w=_g->target(e);
680 if ( level[w] == _node_num && w != _source ) {
682 Node z=level_list[l];
683 if ( z!=INVALID ) left.set(z,w);
695 for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
697 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
699 //Reverse_bfs from t, to find the starting level.
700 level.set(_target,0);
701 bfs_queue.push(_target);
703 while ( !bfs_queue.empty() ) {
705 Node v=bfs_queue.front();
709 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
710 Node w=_g->source(e);
711 if ( level[w] == _node_num && w != _source ) {
713 Node z=level_list[l];
714 if ( z!=INVALID ) left.set(z,w);
723 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
724 Num c=(*_capacity)[e];
725 if ( c <= 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]]);
733 excess.set(w, excess[w]+c);
739 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
742 for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
743 for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
744 excess.set(_target,exc);
748 for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) {
749 Num rem=(*_capacity)[e]-(*_flow)[e];
750 if ( rem <= 0 ) continue;
751 Node w=_g->target(e);
752 if ( level[w] < _node_num ) {
753 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack
754 next.set(w,first[level[w]]);
757 _flow->set(e, (*_capacity)[e]);
758 excess.set(w, excess[w]+rem);
762 for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
763 if ( (*_flow)[e] <= 0 ) continue;
764 Node w=_g->source(e);
765 if ( level[w] < _node_num ) {
766 if ( excess[w] <= 0 && w!=_target ) {
767 next.set(w,first[level[w]]);
770 excess.set(w, excess[w]+(*_flow)[e]);
778 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
779 Num rem=(*_capacity)[e]-(*_flow)[e];
780 if ( rem <= 0 ) continue;
781 Node w=_g->target(e);
782 if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
785 for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
786 if ( (*_flow)[e] <= 0 ) continue;
787 Node w=_g->source(e);
788 if ( level[w] < _node_num ) _flow->set(e, 0);
791 //computing the excess
792 for(NodeIt w(*_g); w!=INVALID; ++w) {
794 for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
795 for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
798 //putting the active nodes into the stack
800 if ( exc > 0 && lev < _node_num && Node(w) != _target ) {
801 next.set(w,first[lev]);
810 void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
811 VecNode& level_list, NNMap& left,
812 NNMap& right, int& b, int& k, bool what_heur )
817 Node right_n=right[w];
821 if ( right_n!=INVALID ) {
822 if ( left_n!=INVALID ) {
823 right.set(left_n, right_n);
824 left.set(right_n, left_n);
826 level_list[lev]=right_n;
827 left.set(right_n, INVALID);
830 if ( left_n!=INVALID ) {
831 right.set(left_n, INVALID);
833 level_list[lev]=INVALID;
838 if ( level_list[lev]==INVALID ) {
841 for (int i=lev; i!=k ; ) {
842 Node v=level_list[++i];
843 while ( v!=INVALID ) {
844 level.set(v,_node_num);
847 level_list[i]=INVALID;
848 if ( !what_heur ) first[i]=INVALID;
851 level.set(w,_node_num);
858 if ( newlevel == _node_num ) level.set(w,_node_num);
860 level.set(w,++newlevel);
861 next.set(w,first[newlevel]);
863 if ( what_heur ) b=newlevel;
864 if ( k < newlevel ) ++k; //now k=newlevel
865 Node z=level_list[newlevel];
866 if ( z!=INVALID ) left.set(z,w);
869 level_list[newlevel]=w;
877 ///\brief Function type interface for Preflow algorithm.
879 ///Function type interface for Preflow algorithm.
881 template<class GR, class CM, class FM>
882 Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
883 typename GR::Node source,
884 typename GR::Node target,
889 return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
894 #endif //LEMON_PREFLOW_H