Interface to the cplex MIP solver: it is little, a bit sour but it is ours.
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/bits/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 ///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.
62 ///\param Tolerance The tolerance type.
64 ///\author Jacint Szabo
65 ///\todo Second template parameter is superfluous
66 template <typename Graph, typename Num,
67 typename CapacityMap=typename Graph::template EdgeMap<Num>,
68 typename FlowMap=typename Graph::template EdgeMap<Num>,
69 typename Tolerance=Tolerance<Num> >
72 typedef typename Graph::Node Node;
73 typedef typename Graph::NodeIt NodeIt;
74 typedef typename Graph::EdgeIt EdgeIt;
75 typedef typename Graph::OutEdgeIt OutEdgeIt;
76 typedef typename Graph::InEdgeIt InEdgeIt;
78 typedef typename Graph::template NodeMap<Node> NNMap;
79 typedef typename std::vector<Node> VecNode;
84 const CapacityMap* _capacity;
89 int _node_num; //the number of nodes of G
91 typename Graph::template NodeMap<int> level;
92 typename Graph::template NodeMap<Num> excess;
94 // constants used for heuristics
95 static const int H0=20;
96 static const int H1=1;
100 ///\ref Exception for the case when s=t.
102 ///\ref Exception for the case when the source equals the target.
103 class InvalidArgument : public lemon::LogicError {
105 virtual const char* what() const throw() {
106 return "lemon::Preflow::InvalidArgument";
111 ///Indicates the property of the starting flow map.
113 ///Indicates the property of the starting flow map.
116 ///indicates an unspecified edge map. \c flow will be
117 ///set to the constant zero flow in the beginning of
118 ///the algorithm in this case.
120 ///constant zero flow
122 ///any flow, i.e. the sum of the in-flows equals to
123 ///the sum of the out-flows in every node except the \c source and
126 ///any preflow, i.e. the sum of the in-flows is at
127 ///least the sum of the out-flows in every node except the \c source.
131 ///Indicates the state of the preflow algorithm.
133 ///Indicates the state of the preflow algorithm.
136 ///before running the algorithm or
137 ///at an unspecified state.
139 ///right after running \ref phase1()
140 AFTER_PREFLOW_PHASE_1,
141 ///after running \ref phase2()
142 AFTER_PREFLOW_PHASE_2
147 StatusEnum status; // Do not needle this flag only if necessary.
150 ///The constructor of the class.
152 ///The constructor of the class.
153 ///\param _gr The directed graph the algorithm runs on.
154 ///\param _s The source node.
155 ///\param _t The target node.
156 ///\param _cap The capacity of the edges.
157 ///\param _f The flow of the edges.
158 ///\param tol Tolerance class.
159 ///Except the graph, all of these parameters can be reset by
160 ///calling \ref source, \ref target, \ref capacityMap and \ref
162 Preflow(const Graph& _gr, Node _s, Node _t,
163 const CapacityMap& _cap, FlowMap& _f,
164 const Tolerance &_sr=Tolerance()) :
165 _g(&_gr), _source(_s), _target(_t), _capacity(&_cap),
166 _flow(&_f), _surely(_sr),
167 _node_num(countNodes(_gr)), level(_gr), excess(_gr,0),
168 flow_prop(NO_FLOW), status(AFTER_NOTHING) {
169 if ( _source==_target )
170 throw InvalidArgument();
173 ///Give a reference to the tolerance handler class
175 ///Give a reference to the tolerance handler class
177 Tolerance &tolerance() { return _surely; }
179 ///Runs the preflow algorithm.
181 ///Runs the preflow algorithm.
188 ///Runs the preflow algorithm.
190 ///Runs the preflow algorithm.
191 ///\pre The starting flow map must be
192 /// - a constant zero flow if \c fp is \c ZERO_FLOW,
193 /// - an arbitrary flow if \c fp is \c GEN_FLOW,
194 /// - an arbitrary preflow if \c fp is \c PRE_FLOW,
195 /// - any map if \c fp is NO_FLOW.
196 ///If the starting flow map is a flow or a preflow then
197 ///the algorithm terminates faster.
198 void run(FlowEnum fp) {
203 ///Runs the first phase of the preflow algorithm.
205 ///The preflow algorithm consists of two phases, this method runs
206 ///the first phase. After the first phase the maximum flow value
207 ///and a minimum value cut can already be computed, although a
208 ///maximum flow is not yet obtained. So after calling this method
209 ///\ref flowValue returns the value of a maximum flow and \ref
210 ///minCut returns a minimum cut.
211 ///\warning \ref minMinCut and \ref maxMinCut do not give minimum
212 ///value cuts unless calling \ref phase2.
213 ///\pre The starting flow must be
214 ///- a constant zero flow if \c fp is \c ZERO_FLOW,
215 ///- an arbitary flow if \c fp is \c GEN_FLOW,
216 ///- an arbitary preflow if \c fp is \c PRE_FLOW,
217 ///- any map if \c fp is NO_FLOW.
218 void phase1(FlowEnum fp)
225 ///Runs the first phase of the preflow algorithm.
227 ///The preflow algorithm consists of two phases, this method runs
228 ///the first phase. After the first phase the maximum flow value
229 ///and a minimum value cut can already be computed, although a
230 ///maximum flow is not yet obtained. So after calling this method
231 ///\ref flowValue returns the value of a maximum flow and \ref
232 ///minCut returns a minimum cut.
233 ///\warning \ref minMinCut() and \ref maxMinCut() do not
234 ///give minimum value cuts unless calling \ref phase2().
237 int heur0=(int)(H0*_node_num); //time while running 'bound decrease'
238 int heur1=(int)(H1*_node_num); //time while running 'highest label'
239 int heur=heur1; //starting time interval (#of relabels)
243 //It is 0 in case 'bound decrease' and 1 in case 'highest label'
246 //Needed for 'bound decrease', true means no active
247 //nodes are above bound b.
249 int k=_node_num-2; //bound on the highest level under n containing a node
250 int b=k; //bound on the highest level under n containing an active node
252 VecNode first(_node_num, INVALID);
253 NNMap next(*_g, INVALID);
255 NNMap left(*_g, INVALID);
256 NNMap right(*_g, INVALID);
257 VecNode level_list(_node_num,INVALID);
258 //List of the nodes in level i<n, set to n.
260 preflowPreproc(first, next, level_list, left, right);
262 //Push/relabel on the highest level active nodes.
265 if ( !what_heur && !end && k > 0 ) {
271 if ( first[b]==INVALID ) --b;
276 int newlevel=push(w, next, first);
277 if ( _surely.positive(excess[w]) ) relabel(w, newlevel, first, next, level_list,
278 left, right, b, k, what_heur);
281 if ( numrelabel >= heur ) {
296 status=AFTER_PREFLOW_PHASE_1;
301 // list 'level_list' on the nodes on level i implemented by hand
302 // stack 'active' on the active nodes on level i
303 // runs heuristic 'highest label' for H1*n relabels
304 // runs heuristic 'bound decrease' for H0*n relabels,
305 // starts with 'highest label'
306 // Parameters H0 and H1 are initialized to 20 and 1.
309 ///Runs the second phase of the preflow algorithm.
311 ///The preflow algorithm consists of two phases, this method runs
312 ///the second phase. After calling \ref phase1() and then
314 /// \ref flowMap() return a maximum flow, \ref flowValue
315 ///returns the value of a maximum flow, \ref minCut returns a
316 ///minimum cut, while the methods \ref minMinCut and \ref
317 ///maxMinCut return the inclusionwise minimum and maximum cuts of
318 ///minimum value, resp. \pre \ref phase1 must be called before.
322 int k=_node_num-2; //bound on the highest level under n containing a node
323 int b=k; //bound on the highest level under n of an active node
326 VecNode first(_node_num, INVALID);
327 NNMap next(*_g, INVALID);
328 level.set(_source,0);
329 std::queue<Node> bfs_queue;
330 bfs_queue.push(_source);
332 while ( !bfs_queue.empty() ) {
334 Node v=bfs_queue.front();
338 for(InEdgeIt e(*_g,v); e!=INVALID; ++e) {
339 if ( !_surely.less((*_flow)[e], (*_capacity)[e]) ) continue;
340 Node u=_g->source(e);
341 if ( level[u] >= _node_num ) {
344 if ( _surely.positive(excess[u]) ) {
345 next.set(u,first[l]);
351 for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) {
352 if ( !_surely.positive((*_flow)[e]) ) continue;
353 Node u=_g->target(e);
354 if ( level[u] >= _node_num ) {
357 if ( _surely.positive(excess[u]) ) {
358 next.set(u,first[l]);
369 if ( first[b]==INVALID ) --b;
373 int newlevel=push(w,next, first);
376 if ( _surely.positive(excess[w]) ) {
377 level.set(w,++newlevel);
378 next.set(w,first[newlevel]);
385 status=AFTER_PREFLOW_PHASE_2;
388 /// Returns the value of the maximum flow.
390 /// Returns the value of the maximum flow by returning the excess
391 /// of the target node \c t. This value equals to the value of
392 /// the maximum flow already after running \ref phase1.
393 Num flowValue() const {
394 return excess[_target];
398 ///Returns a minimum value cut.
400 ///Sets \c M to the characteristic vector of a minimum value
401 ///cut. This method can be called both after running \ref
402 ///phase1 and \ref phase2. It is much faster after
403 ///\ref phase1. \pre M should be a bool-valued node-map. \pre
404 ///If \ref minCut() is called after \ref phase2() then M should
405 ///be initialized to false.
406 template<typename _CutMap>
407 void minCut(_CutMap& M) const {
409 case AFTER_PREFLOW_PHASE_1:
410 for(NodeIt v(*_g); v!=INVALID; ++v) {
411 if (level[v] < _node_num) {
418 case AFTER_PREFLOW_PHASE_2:
426 ///Returns the inclusionwise minimum of the minimum value cuts.
428 ///Sets \c M to the characteristic vector of the minimum value cut
429 ///which is inclusionwise minimum. It is computed by processing a
430 ///bfs from the source node \c s in the residual graph. \pre M
431 ///should be a node map of bools initialized to false. \pre \ref
432 ///phase2 should already be run.
433 template<typename _CutMap>
434 void minMinCut(_CutMap& M) const {
436 std::queue<Node> queue;
440 while (!queue.empty()) {
441 Node w=queue.front();
444 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
445 Node v=_g->target(e);
446 if (!M[v] && _surely.less((*_flow)[e] , (*_capacity)[e]) ) {
452 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
453 Node v=_g->source(e);
454 if (!M[v] && _surely.positive((*_flow)[e]) ) {
462 ///Returns the inclusionwise maximum of the minimum value cuts.
464 ///Sets \c M to the characteristic vector of the minimum value cut
465 ///which is inclusionwise maximum. It is computed by processing a
466 ///backward bfs from the target node \c t in the residual graph.
467 ///\pre \ref phase2() or run() should already be run.
468 template<typename _CutMap>
469 void maxMinCut(_CutMap& M) const {
471 for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true);
473 std::queue<Node> queue;
475 M.set(_target,false);
478 while (!queue.empty()) {
479 Node w=queue.front();
482 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
483 Node v=_g->source(e);
484 if (M[v] && _surely.less((*_flow)[e], (*_capacity)[e]) ) {
490 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
491 Node v=_g->target(e);
492 if (M[v] && _surely.positive((*_flow)[e]) ) {
500 ///Sets the source node to \c _s.
502 ///Sets the source node to \c _s.
504 void source(Node _s) {
506 if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW;
507 status=AFTER_NOTHING;
510 ///Returns the source node.
512 ///Returns the source node.
514 Node source() const {
518 ///Sets the target node to \c _t.
520 ///Sets the target node to \c _t.
522 void target(Node _t) {
524 if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW;
525 status=AFTER_NOTHING;
528 ///Returns the target node.
530 ///Returns the target node.
532 Node target() const {
536 /// Sets the edge map of the capacities to _cap.
538 /// Sets the edge map of the capacities to _cap.
540 void capacityMap(const CapacityMap& _cap) {
542 status=AFTER_NOTHING;
544 /// Returns a reference to capacity map.
546 /// Returns a reference to capacity map.
548 const CapacityMap &capacityMap() const {
552 /// Sets the edge map of the flows to _flow.
554 /// Sets the edge map of the flows to _flow.
556 void flowMap(FlowMap& _f) {
559 status=AFTER_NOTHING;
562 /// Returns a reference to flow map.
564 /// Returns a reference to flow map.
566 const FlowMap &flowMap() const {
572 int push(Node w, NNMap& next, VecNode& first) {
576 int newlevel=_node_num; //bound on the next level of w
578 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
579 if ( !_surely.less((*_flow)[e], (*_capacity)[e]) ) continue;
580 Node v=_g->target(e);
582 if( lev > level[v] ) { //Push is allowed now
584 if ( !_surely.positive(excess[v]) && v!=_target && v!=_source ) {
585 next.set(v,first[level[v]]);
589 Num cap=(*_capacity)[e];
593 if ( ! _surely.less(remcap, exc) ) { //A nonsaturating push.
595 _flow->set(e, flo+exc);
596 excess.set(v, excess[v]+exc);
600 } else { //A saturating push.
602 excess.set(v, excess[v]+remcap);
605 } else if ( newlevel > level[v] ) newlevel = level[v];
608 if ( _surely.positive(exc) ) {
609 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) {
611 if ( !_surely.positive((*_flow)[e]) ) continue;
612 Node v=_g->source(e);
614 if( lev > level[v] ) { //Push is allowed now
616 if ( !_surely.positive(excess[v]) && v!=_target && v!=_source ) {
617 next.set(v,first[level[v]]);
623 if ( !_surely.less(flo, exc) ) { //A nonsaturating push.
625 _flow->set(e, flo-exc);
626 excess.set(v, excess[v]+exc);
629 } else { //A saturating push.
631 excess.set(v, excess[v]+flo);
635 } else if ( newlevel > level[v] ) newlevel = level[v];
638 } // if w still has excess after the out edge for cycle
647 void preflowPreproc(VecNode& first, NNMap& next,
648 VecNode& level_list, NNMap& left, NNMap& right)
650 for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num);
651 std::queue<Node> bfs_queue;
653 if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) {
654 //Reverse_bfs from t in the residual graph,
655 //to find the starting level.
656 level.set(_target,0);
657 bfs_queue.push(_target);
659 while ( !bfs_queue.empty() ) {
661 Node v=bfs_queue.front();
665 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
666 if ( !_surely.less((*_flow)[e],(*_capacity)[e]) ) continue;
667 Node w=_g->source(e);
668 if ( level[w] == _node_num && w != _source ) {
670 Node z=level_list[l];
671 if ( z!=INVALID ) left.set(z,w);
678 for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
679 if ( !_surely.positive((*_flow)[e]) ) continue;
680 Node w=_g->target(e);
681 if ( level[w] == _node_num && w != _source ) {
683 Node z=level_list[l];
684 if ( z!=INVALID ) left.set(z,w);
696 for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0);
698 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
700 //Reverse_bfs from t, to find the starting level.
701 level.set(_target,0);
702 bfs_queue.push(_target);
704 while ( !bfs_queue.empty() ) {
706 Node v=bfs_queue.front();
710 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) {
711 Node w=_g->source(e);
712 if ( level[w] == _node_num && w != _source ) {
714 Node z=level_list[l];
715 if ( z!=INVALID ) left.set(z,w);
724 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
725 Num c=(*_capacity)[e];
726 if ( !_surely.positive(c) ) continue;
727 Node w=_g->target(e);
728 if ( level[w] < _node_num ) {
729 if ( !_surely.positive(excess[w]) && w!=_target ) { //putting into the stack
730 next.set(w,first[level[w]]);
734 excess.set(w, excess[w]+c);
740 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0);
743 for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e];
744 for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e];
745 excess.set(_target,exc);
749 for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) {
750 Num rem=(*_capacity)[e]-(*_flow)[e];
751 if ( !_surely.positive(rem) ) continue;
752 Node w=_g->target(e);
753 if ( level[w] < _node_num ) {
754 if ( !_surely.positive(excess[w]) && w!=_target ) { //putting into the stack
755 next.set(w,first[level[w]]);
758 _flow->set(e, (*_capacity)[e]);
759 excess.set(w, excess[w]+rem);
763 for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) {
764 if ( !_surely.positive((*_flow)[e]) ) continue;
765 Node w=_g->source(e);
766 if ( level[w] < _node_num ) {
767 if ( !_surely.positive(excess[w]) && w!=_target ) {
768 next.set(w,first[level[w]]);
771 excess.set(w, excess[w]+(*_flow)[e]);
779 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
780 Num rem=(*_capacity)[e]-(*_flow)[e];
781 if ( !_surely.positive(rem) ) continue;
782 Node w=_g->target(e);
783 if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]);
786 for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) {
787 if ( !_surely.positive((*_flow)[e]) ) continue;
788 Node w=_g->source(e);
789 if ( level[w] < _node_num ) _flow->set(e, 0);
792 //computing the excess
793 for(NodeIt w(*_g); w!=INVALID; ++w) {
795 for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e];
796 for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e];
799 //putting the active nodes into the stack
801 if ( _surely.positive(exc) && lev < _node_num && Node(w) != _target ) {
802 next.set(w,first[lev]);
811 void relabel(Node w, int newlevel, VecNode& first, NNMap& next,
812 VecNode& level_list, NNMap& left,
813 NNMap& right, int& b, int& k, bool what_heur )
818 Node right_n=right[w];
822 if ( right_n!=INVALID ) {
823 if ( left_n!=INVALID ) {
824 right.set(left_n, right_n);
825 left.set(right_n, left_n);
827 level_list[lev]=right_n;
828 left.set(right_n, INVALID);
831 if ( left_n!=INVALID ) {
832 right.set(left_n, INVALID);
834 level_list[lev]=INVALID;
839 if ( level_list[lev]==INVALID ) {
842 for (int i=lev; i!=k ; ) {
843 Node v=level_list[++i];
844 while ( v!=INVALID ) {
845 level.set(v,_node_num);
848 level_list[i]=INVALID;
849 if ( !what_heur ) first[i]=INVALID;
852 level.set(w,_node_num);
859 if ( newlevel == _node_num ) level.set(w,_node_num);
861 level.set(w,++newlevel);
862 next.set(w,first[newlevel]);
864 if ( what_heur ) b=newlevel;
865 if ( k < newlevel ) ++k; //now k=newlevel
866 Node z=level_list[newlevel];
867 if ( z!=INVALID ) left.set(z,w);
870 level_list[newlevel]=w;
878 ///\brief Function type interface for Preflow algorithm.
880 ///Function type interface for Preflow algorithm.
882 template<class GR, class CM, class FM>
883 Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g,
884 typename GR::Node source,
885 typename GR::Node target,
890 return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow);
895 #endif //LEMON_PREFLOW_H