1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2013
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_CIRCULATION_H
20 #define LEMON_CIRCULATION_H
22 #include <lemon/tolerance.h>
23 #include <lemon/elevator.h>
28 ///\brief Push-relabel algorithm for finding a feasible circulation.
32 /// \brief Default traits class of Circulation class.
34 /// Default traits class of Circulation class.
36 /// \tparam GR Type of the digraph the algorithm runs on.
37 /// \tparam LM The type of the lower bound map.
38 /// \tparam UM The type of the upper bound (capacity) map.
39 /// \tparam SM The type of the supply map.
40 template <typename GR, typename LM,
41 typename UM, typename SM>
42 struct CirculationDefaultTraits {
44 /// \brief The type of the digraph the algorithm runs on.
47 /// \brief The type of the lower bound map.
49 /// The type of the map that stores the lower bounds on the arcs.
50 /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
53 /// \brief The type of the upper bound (capacity) map.
55 /// The type of the map that stores the upper bounds (capacities)
57 /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
60 /// \brief The type of supply map.
62 /// The type of the map that stores the signed supply values of the
64 /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
67 /// \brief The type of the flow and supply values.
68 typedef typename SupplyMap::Value Value;
70 /// \brief The type of the map that stores the flow values.
72 /// The type of the map that stores the flow values.
73 /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
76 typedef GR::ArcMap<Value> FlowMap;
78 typedef typename Digraph::template ArcMap<Value> FlowMap;
81 /// \brief Instantiates a FlowMap.
83 /// This function instantiates a \ref FlowMap.
84 /// \param digraph The digraph for which we would like to define
86 static FlowMap* createFlowMap(const Digraph& digraph) {
87 return new FlowMap(digraph);
90 /// \brief The elevator type used by the algorithm.
92 /// The elevator type used by the algorithm.
94 /// \sa Elevator, LinkedElevator
96 typedef lemon::Elevator<GR, GR::Node> Elevator;
98 typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
101 /// \brief Instantiates an Elevator.
103 /// This function instantiates an \ref Elevator.
104 /// \param digraph The digraph for which we would like to define
106 /// \param max_level The maximum level of the elevator.
107 static Elevator* createElevator(const Digraph& digraph, int max_level) {
108 return new Elevator(digraph, max_level);
111 /// \brief The tolerance used by the algorithm
113 /// The tolerance used by the algorithm to handle inexact computation.
114 typedef lemon::Tolerance<Value> Tolerance;
119 \brief Push-relabel algorithm for the network circulation problem.
122 This class implements a push-relabel algorithm for the \e network
123 \e circulation problem.
124 It is to find a feasible circulation when lower and upper bounds
125 are given for the flow values on the arcs and lower bounds are
126 given for the difference between the outgoing and incoming flow
129 The exact formulation of this problem is the following.
130 Let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$
131 \f$upper: A\rightarrow\mathbf{R}\cup\{\infty\}\f$ denote the lower and
132 upper bounds on the arcs, for which \f$lower(uv) \leq upper(uv)\f$
133 holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$
134 denotes the signed supply values of the nodes.
135 If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
136 supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
137 \f$-sup(u)\f$ demand.
138 A feasible circulation is an \f$f: A\rightarrow\mathbf{R}\f$
139 solution of the following problem.
141 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu)
142 \geq sup(u) \quad \forall u\in V, \f]
143 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f]
145 The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
146 zero or negative in order to have a feasible solution (since the sum
147 of the expressions on the left-hand side of the inequalities is zero).
148 It means that the total demand must be greater or equal to the total
149 supply and all the supplies have to be carried out from the supply nodes,
150 but there could be demands that are not satisfied.
151 If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
152 constraints have to be satisfied with equality, i.e. all demands
153 have to be satisfied and all supplies have to be used.
155 If you need the opposite inequalities in the supply/demand constraints
156 (i.e. the total demand is less than the total supply and all the demands
157 have to be satisfied while there could be supplies that are not used),
158 then you could easily transform the problem to the above form by reversing
159 the direction of the arcs and taking the negative of the supply values
160 (e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
162 This algorithm either calculates a feasible circulation, or provides
163 a \ref barrier() "barrier", which prooves that a feasible soultion
166 Note that this algorithm also provides a feasible solution for the
167 \ref min_cost_flow "minimum cost flow problem".
169 \tparam GR The type of the digraph the algorithm runs on.
170 \tparam LM The type of the lower bound map. The default
171 map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
172 \tparam UM The type of the upper bound (capacity) map.
173 The default map type is \c LM.
174 \tparam SM The type of the supply map. The default map type is
175 \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
176 \tparam TR The traits class that defines various types used by the
177 algorithm. By default, it is \ref CirculationDefaultTraits
178 "CirculationDefaultTraits<GR, LM, UM, SM>".
179 In most cases, this parameter should not be set directly,
180 consider to use the named template parameters instead.
183 template< typename GR,
189 template< typename GR,
190 typename LM = typename GR::template ArcMap<int>,
192 typename SM = typename GR::template NodeMap<typename UM::Value>,
193 typename TR = CirculationDefaultTraits<GR, LM, UM, SM> >
198 /// \brief The \ref lemon::CirculationDefaultTraits "traits class"
199 /// of the algorithm.
201 ///The type of the digraph the algorithm runs on.
202 typedef typename Traits::Digraph Digraph;
203 ///The type of the flow and supply values.
204 typedef typename Traits::Value Value;
206 ///The type of the lower bound map.
207 typedef typename Traits::LowerMap LowerMap;
208 ///The type of the upper bound (capacity) map.
209 typedef typename Traits::UpperMap UpperMap;
210 ///The type of the supply map.
211 typedef typename Traits::SupplyMap SupplyMap;
212 ///The type of the flow map.
213 typedef typename Traits::FlowMap FlowMap;
215 ///The type of the elevator.
216 typedef typename Traits::Elevator Elevator;
217 ///The type of the tolerance.
218 typedef typename Traits::Tolerance Tolerance;
222 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
229 const SupplyMap *_supply;
237 typedef typename Digraph::template NodeMap<Value> ExcessMap;
245 typedef Circulation Create;
247 ///\name Named Template Parameters
251 template <typename T>
252 struct SetFlowMapTraits : public Traits {
254 static FlowMap *createFlowMap(const Digraph&) {
255 LEMON_ASSERT(false, "FlowMap is not initialized");
256 return 0; // ignore warnings
260 /// \brief \ref named-templ-param "Named parameter" for setting
263 /// \ref named-templ-param "Named parameter" for setting FlowMap
265 template <typename T>
267 : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
268 SetFlowMapTraits<T> > {
269 typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
270 SetFlowMapTraits<T> > Create;
273 template <typename T>
274 struct SetElevatorTraits : public Traits {
276 static Elevator *createElevator(const Digraph&, int) {
277 LEMON_ASSERT(false, "Elevator is not initialized");
278 return 0; // ignore warnings
282 /// \brief \ref named-templ-param "Named parameter" for setting
285 /// \ref named-templ-param "Named parameter" for setting Elevator
286 /// type. If this named parameter is used, then an external
287 /// elevator object must be passed to the algorithm using the
288 /// \ref elevator(Elevator&) "elevator()" function before calling
289 /// \ref run() or \ref init().
290 /// \sa SetStandardElevator
291 template <typename T>
293 : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
294 SetElevatorTraits<T> > {
295 typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
296 SetElevatorTraits<T> > Create;
299 template <typename T>
300 struct SetStandardElevatorTraits : public Traits {
302 static Elevator *createElevator(const Digraph& digraph, int max_level) {
303 return new Elevator(digraph, max_level);
307 /// \brief \ref named-templ-param "Named parameter" for setting
308 /// Elevator type with automatic allocation
310 /// \ref named-templ-param "Named parameter" for setting Elevator
311 /// type with automatic allocation.
312 /// The Elevator should have standard constructor interface to be
313 /// able to automatically created by the algorithm (i.e. the
314 /// digraph and the maximum level should be passed to it).
315 /// However, an external elevator object could also be passed to the
316 /// algorithm with the \ref elevator(Elevator&) "elevator()" function
317 /// before calling \ref run() or \ref init().
319 template <typename T>
320 struct SetStandardElevator
321 : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
322 SetStandardElevatorTraits<T> > {
323 typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
324 SetStandardElevatorTraits<T> > Create;
337 /// The constructor of the class.
339 /// \param graph The digraph the algorithm runs on.
340 /// \param lower The lower bounds for the flow values on the arcs.
341 /// \param upper The upper bounds (capacities) for the flow values
343 /// \param supply The signed supply values of the nodes.
344 Circulation(const Digraph &graph, const LowerMap &lower,
345 const UpperMap &upper, const SupplyMap &supply)
346 : _g(graph), _lo(&lower), _up(&upper), _supply(&supply),
347 _flow(NULL), _local_flow(false), _level(NULL), _local_level(false),
358 bool checkBoundMaps() {
359 for (ArcIt e(_g);e!=INVALID;++e) {
360 if (_tol.less((*_up)[e], (*_lo)[e])) return false;
365 void createStructures() {
366 _node_num = _el = countNodes(_g);
369 _flow = Traits::createFlowMap(_g);
373 _level = Traits::createElevator(_g, _node_num);
377 _excess = new ExcessMap(_g);
381 void destroyStructures() {
395 /// Sets the lower bound map.
397 /// Sets the lower bound map.
398 /// \return <tt>(*this)</tt>
399 Circulation& lowerMap(const LowerMap& map) {
404 /// Sets the upper bound (capacity) map.
406 /// Sets the upper bound (capacity) map.
407 /// \return <tt>(*this)</tt>
408 Circulation& upperMap(const UpperMap& map) {
413 /// Sets the supply map.
415 /// Sets the supply map.
416 /// \return <tt>(*this)</tt>
417 Circulation& supplyMap(const SupplyMap& map) {
422 /// \brief Sets the flow map.
424 /// Sets the flow map.
425 /// If you don't use this function before calling \ref run() or
426 /// \ref init(), an instance will be allocated automatically.
427 /// The destructor deallocates this automatically allocated map,
429 /// \return <tt>(*this)</tt>
430 Circulation& flowMap(FlowMap& map) {
439 /// \brief Sets the elevator used by algorithm.
441 /// Sets the elevator used by algorithm.
442 /// If you don't use this function before calling \ref run() or
443 /// \ref init(), an instance will be allocated automatically.
444 /// The destructor deallocates this automatically allocated elevator,
446 /// \return <tt>(*this)</tt>
447 Circulation& elevator(Elevator& elevator) {
450 _local_level = false;
456 /// \brief Returns a const reference to the elevator.
458 /// Returns a const reference to the elevator.
460 /// \pre Either \ref run() or \ref init() must be called before
461 /// using this function.
462 const Elevator& elevator() const {
466 /// \brief Sets the tolerance used by the algorithm.
468 /// Sets the tolerance object used by the algorithm.
469 /// \return <tt>(*this)</tt>
470 Circulation& tolerance(const Tolerance& tolerance) {
475 /// \brief Returns a const reference to the tolerance.
477 /// Returns a const reference to the tolerance object used by
479 const Tolerance& tolerance() const {
483 /// \name Execution Control
484 /// The simplest way to execute the algorithm is to call \ref run().\n
485 /// If you need better control on the initial solution or the execution,
486 /// you have to call one of the \ref init() functions first, then
487 /// the \ref start() function.
491 /// Initializes the internal data structures.
493 /// Initializes the internal data structures and sets all flow values
494 /// to the lower bound.
497 LEMON_DEBUG(checkBoundMaps(),
498 "Upper bounds must be greater or equal to the lower bounds");
502 for(NodeIt n(_g);n!=INVALID;++n) {
503 (*_excess)[n] = (*_supply)[n];
506 for (ArcIt e(_g);e!=INVALID;++e) {
507 _flow->set(e, (*_lo)[e]);
508 (*_excess)[_g.target(e)] += (*_flow)[e];
509 (*_excess)[_g.source(e)] -= (*_flow)[e];
512 // global relabeling tested, but in general case it provides
513 // worse performance for random digraphs
515 for(NodeIt n(_g);n!=INVALID;++n)
516 _level->initAddItem(n);
517 _level->initFinish();
518 for(NodeIt n(_g);n!=INVALID;++n)
519 if(_tol.positive((*_excess)[n]))
523 /// Initializes the internal data structures using a greedy approach.
525 /// Initializes the internal data structures using a greedy approach
526 /// to construct the initial solution.
529 LEMON_DEBUG(checkBoundMaps(),
530 "Upper bounds must be greater or equal to the lower bounds");
534 for(NodeIt n(_g);n!=INVALID;++n) {
535 (*_excess)[n] = (*_supply)[n];
538 for (ArcIt e(_g);e!=INVALID;++e) {
539 if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {
540 _flow->set(e, (*_up)[e]);
541 (*_excess)[_g.target(e)] += (*_up)[e];
542 (*_excess)[_g.source(e)] -= (*_up)[e];
543 } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {
544 _flow->set(e, (*_lo)[e]);
545 (*_excess)[_g.target(e)] += (*_lo)[e];
546 (*_excess)[_g.source(e)] -= (*_lo)[e];
548 Value fc = -(*_excess)[_g.target(e)];
550 (*_excess)[_g.target(e)] = 0;
551 (*_excess)[_g.source(e)] -= fc;
556 for(NodeIt n(_g);n!=INVALID;++n)
557 _level->initAddItem(n);
558 _level->initFinish();
559 for(NodeIt n(_g);n!=INVALID;++n)
560 if(_tol.positive((*_excess)[n]))
564 ///Executes the algorithm
566 ///This function executes the algorithm.
568 ///\return \c true if a feasible circulation is found.
576 while((act=_level->highestActive())!=INVALID) {
577 int actlevel=(*_level)[act];
578 int mlevel=_node_num;
579 Value exc=(*_excess)[act];
581 for(OutArcIt e(_g,act);e!=INVALID; ++e) {
582 Node v = _g.target(e);
583 Value fc=(*_up)[e]-(*_flow)[e];
584 if(!_tol.positive(fc)) continue;
585 if((*_level)[v]<actlevel) {
586 if(!_tol.less(fc, exc)) {
587 _flow->set(e, (*_flow)[e] + exc);
588 (*_excess)[v] += exc;
589 if(!_level->active(v) && _tol.positive((*_excess)[v]))
592 _level->deactivate(act);
596 _flow->set(e, (*_up)[e]);
598 if(!_level->active(v) && _tol.positive((*_excess)[v]))
603 else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
605 for(InArcIt e(_g,act);e!=INVALID; ++e) {
606 Node v = _g.source(e);
607 Value fc=(*_flow)[e]-(*_lo)[e];
608 if(!_tol.positive(fc)) continue;
609 if((*_level)[v]<actlevel) {
610 if(!_tol.less(fc, exc)) {
611 _flow->set(e, (*_flow)[e] - exc);
612 (*_excess)[v] += exc;
613 if(!_level->active(v) && _tol.positive((*_excess)[v]))
616 _level->deactivate(act);
620 _flow->set(e, (*_lo)[e]);
622 if(!_level->active(v) && _tol.positive((*_excess)[v]))
627 else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
630 (*_excess)[act] = exc;
631 if(!_tol.positive(exc)) _level->deactivate(act);
632 else if(mlevel==_node_num) {
633 _level->liftHighestActiveToTop();
638 _level->liftHighestActive(mlevel+1);
639 if(_level->onLevel(actlevel)==0) {
650 /// Runs the algorithm.
652 /// This function runs the algorithm.
654 /// \return \c true if a feasible circulation is found.
656 /// \note Apart from the return value, c.run() is just a shortcut of
657 /// the following code.
669 /// \name Query Functions
670 /// The results of the circulation algorithm can be obtained using
671 /// these functions.\n
672 /// Either \ref run() or \ref start() should be called before
677 /// \brief Returns the flow value on the given arc.
679 /// Returns the flow value on the given arc.
681 /// \pre Either \ref run() or \ref init() must be called before
682 /// using this function.
683 Value flow(const Arc& arc) const {
684 return (*_flow)[arc];
687 /// \brief Returns a const reference to the flow map.
689 /// Returns a const reference to the arc map storing the found flow.
691 /// \pre Either \ref run() or \ref init() must be called before
692 /// using this function.
693 const FlowMap& flowMap() const {
698 \brief Returns \c true if the given node is in a barrier.
700 Barrier is a set \e B of nodes for which
702 \f[ \sum_{uv\in A: u\in B} upper(uv) -
703 \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
705 holds. The existence of a set with this property prooves that a
706 feasible circualtion cannot exist.
708 This function returns \c true if the given node is in the found
709 barrier. If a feasible circulation is found, the function
710 gives back \c false for every node.
712 \pre Either \ref run() or \ref init() must be called before
718 bool barrier(const Node& node) const
720 return (*_level)[node] >= _el;
723 /// \brief Gives back a barrier.
725 /// This function sets \c bar to the characteristic vector of the
726 /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
727 /// node map with \c bool (or convertible) value type.
729 /// If a feasible circulation is found, the function gives back an
730 /// empty set, so \c bar[v] will be \c false for all nodes \c v.
732 /// \note This function calls \ref barrier() for each node,
733 /// so it runs in O(n) time.
735 /// \pre Either \ref run() or \ref init() must be called before
736 /// using this function.
739 /// \sa checkBarrier()
740 template<class BarrierMap>
741 void barrierMap(BarrierMap &bar) const
743 for(NodeIt n(_g);n!=INVALID;++n)
744 bar.set(n, (*_level)[n] >= _el);
749 /// \name Checker Functions
750 /// The feasibility of the results can be checked using
751 /// these functions.\n
752 /// Either \ref run() or \ref start() should be called before
757 ///Check if the found flow is a feasible circulation
759 ///Check if the found flow is a feasible circulation,
761 bool checkFlow() const {
762 for(ArcIt e(_g);e!=INVALID;++e)
763 if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
764 for(NodeIt n(_g);n!=INVALID;++n)
766 Value dif=-(*_supply)[n];
767 for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
768 for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
769 if(_tol.negative(dif)) return false;
774 ///Check whether or not the last execution provides a barrier
776 ///Check whether or not the last execution provides a barrier.
779 bool checkBarrier() const
782 Value inf_cap = std::numeric_limits<Value>::has_infinity ?
783 std::numeric_limits<Value>::infinity() :
784 std::numeric_limits<Value>::max();
785 for(NodeIt n(_g);n!=INVALID;++n)
787 delta-=(*_supply)[n];
788 for(ArcIt e(_g);e!=INVALID;++e)
792 if(barrier(s)&&!barrier(t)) {
793 if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
796 else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
798 return _tol.negative(delta);