3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2007
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
22 #include <lemon/error.h>
23 #include <lemon/bits/invalid.h>
24 #include <lemon/tolerance.h>
25 #include <lemon/maps.h>
26 #include <lemon/graph_utils.h>
27 #include <lemon/elevator.h>
31 /// \brief Implementation of the preflow algorithm.
35 /// \brief Default traits class of Preflow class.
37 /// Default traits class of Preflow class.
38 /// \param _Graph Graph type.
39 /// \param _CapacityMap Type of capacity map.
40 template <typename _Graph, typename _CapacityMap>
41 struct PreflowDefaultTraits {
43 /// \brief The graph type the algorithm runs on.
46 /// \brief The type of the map that stores the edge capacities.
48 /// The type of the map that stores the edge capacities.
49 /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
50 typedef _CapacityMap CapacityMap;
52 /// \brief The type of the length of the edges.
53 typedef typename CapacityMap::Value Value;
55 /// \brief The map type that stores the flow values.
57 /// The map type that stores the flow values.
58 /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
59 typedef typename Graph::template EdgeMap<Value> FlowMap;
61 /// \brief Instantiates a FlowMap.
63 /// This function instantiates a \ref FlowMap.
64 /// \param graph The graph, to which we would like to define the flow map.
65 static FlowMap* createFlowMap(const Graph& graph) {
66 return new FlowMap(graph);
69 /// \brief The eleavator type used by Preflow algorithm.
71 /// The elevator type used by Preflow algorithm.
74 /// \sa LinkedElevator
75 typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
77 /// \brief Instantiates an Elevator.
79 /// This function instantiates a \ref Elevator.
80 /// \param graph The graph, to which we would like to define the elevator.
81 /// \param max_level The maximum level of the elevator.
82 static Elevator* createElevator(const Graph& graph, int max_level) {
83 return new Elevator(graph, max_level);
86 /// \brief The tolerance used by the algorithm
88 /// The tolerance used by the algorithm to handle inexact computation.
89 typedef Tolerance<Value> Tolerance;
96 /// \brief %Preflow algorithms class.
98 /// This class provides an implementation of the Goldberg's \e
99 /// preflow \e algorithm producing a flow of maximum value in a
100 /// directed graph. The preflow algorithms are the fastest known max
101 /// flow algorithms. The current implementation use a mixture of the
102 /// \e "highest label" and the \e "bound decrease" heuristics.
103 /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
105 /// The algorithm consists from two phases. After the first phase
106 /// the maximal flow value and the minimum cut can be obtained. The
107 /// second phase constructs the feasible maximum flow on each edge.
109 /// \param _Graph The directed graph type the algorithm runs on.
110 /// \param _CapacityMap The flow map type.
111 /// \param _Traits Traits class to set various data types used by
112 /// the algorithm. The default traits class is \ref
113 /// PreflowDefaultTraits. See \ref PreflowDefaultTraits for the
114 /// documentation of a %Preflow traits class.
116 ///\author Jacint Szabo and Balazs Dezso
118 template <typename _Graph, typename _CapacityMap, typename _Traits>
120 template <typename _Graph,
121 typename _CapacityMap = typename _Graph::template EdgeMap<int>,
122 typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
127 typedef _Traits Traits;
128 typedef typename Traits::Graph Graph;
129 typedef typename Traits::CapacityMap CapacityMap;
130 typedef typename Traits::Value Value;
132 typedef typename Traits::FlowMap FlowMap;
133 typedef typename Traits::Elevator Elevator;
134 typedef typename Traits::Tolerance Tolerance;
136 /// \brief \ref Exception for uninitialized parameters.
138 /// This error represents problems in the initialization
139 /// of the parameters of the algorithms.
140 class UninitializedParameter : public lemon::UninitializedParameter {
142 virtual const char* what() const throw() {
143 return "lemon::Preflow::UninitializedParameter";
149 GRAPH_TYPEDEFS(typename Graph);
152 const CapacityMap* _capacity;
156 Node _source, _target;
164 typedef typename Graph::template NodeMap<Value> ExcessMap;
167 Tolerance _tolerance;
171 void createStructures() {
172 _node_num = countNodes(_graph);
175 _flow = Traits::createFlowMap(_graph);
179 _level = Traits::createElevator(_graph, _node_num);
183 _excess = new ExcessMap(_graph);
187 void destroyStructures() {
201 typedef Preflow Create;
203 ///\name Named template parameters
207 template <typename _FlowMap>
208 struct DefFlowMapTraits : public Traits {
209 typedef _FlowMap FlowMap;
210 static FlowMap *createFlowMap(const Graph&) {
211 throw UninitializedParameter();
215 /// \brief \ref named-templ-param "Named parameter" for setting
218 /// \ref named-templ-param "Named parameter" for setting FlowMap
220 template <typename _FlowMap>
222 : public Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
223 typedef Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > Create;
226 template <typename _Elevator>
227 struct DefElevatorTraits : public Traits {
228 typedef _Elevator Elevator;
229 static Elevator *createElevator(const Graph&, int) {
230 throw UninitializedParameter();
234 /// \brief \ref named-templ-param "Named parameter" for setting
237 /// \ref named-templ-param "Named parameter" for setting Elevator
239 template <typename _Elevator>
241 : public Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > {
242 typedef Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > Create;
245 template <typename _Elevator>
246 struct DefStandardElevatorTraits : public Traits {
247 typedef _Elevator Elevator;
248 static Elevator *createElevator(const Graph& graph, int max_level) {
249 return new Elevator(graph, max_level);
253 /// \brief \ref named-templ-param "Named parameter" for setting
256 /// \ref named-templ-param "Named parameter" for setting Elevator
257 /// type. The Elevator should be standard constructor interface, ie.
258 /// the graph and the maximum level should be passed to it.
259 template <typename _Elevator>
260 struct DefStandardElevator
261 : public Preflow<Graph, CapacityMap,
262 DefStandardElevatorTraits<_Elevator> > {
263 typedef Preflow<Graph, CapacityMap,
264 DefStandardElevatorTraits<_Elevator> > Create;
269 /// \brief The constructor of the class.
271 /// The constructor of the class.
272 /// \param graph The directed graph the algorithm runs on.
273 /// \param capacity The capacity of the edges.
274 /// \param source The source node.
275 /// \param target The target node.
276 Preflow(const Graph& graph, const CapacityMap& capacity,
277 Node source, Node target)
278 : _graph(graph), _capacity(&capacity),
279 _node_num(0), _source(source), _target(target),
280 _flow(0), _local_flow(false),
281 _level(0), _local_level(false),
282 _excess(0), _tolerance(), _phase() {}
284 /// \brief Destrcutor.
291 /// \brief Sets the capacity map.
293 /// Sets the capacity map.
294 /// \return \c (*this)
295 Preflow& capacityMap(const CapacityMap& map) {
300 /// \brief Sets the flow map.
302 /// Sets the flow map.
303 /// \return \c (*this)
304 Preflow& flowMap(FlowMap& map) {
313 /// \brief Returns the flow map.
315 /// \return The flow map.
316 const FlowMap& flowMap() {
320 /// \brief Sets the elevator.
322 /// Sets the elevator.
323 /// \return \c (*this)
324 Preflow& elevator(Elevator& elevator) {
327 _local_level = false;
333 /// \brief Returns the elevator.
335 /// \return The elevator.
336 const Elevator& elevator() {
340 /// \brief Sets the source node.
342 /// Sets the source node.
343 /// \return \c (*this)
344 Preflow& source(const Node& node) {
349 /// \brief Sets the target node.
351 /// Sets the target node.
352 /// \return \c (*this)
353 Preflow& target(const Node& node) {
358 /// \brief Sets the tolerance used by algorithm.
360 /// Sets the tolerance used by algorithm.
361 Preflow& tolerance(const Tolerance& tolerance) const {
362 _tolerance = tolerance;
366 /// \brief Returns the tolerance used by algorithm.
368 /// Returns the tolerance used by algorithm.
369 const Tolerance& tolerance() const {
373 /// \name Execution control The simplest way to execute the
374 /// algorithm is to use one of the member functions called \c
377 /// If you need more control on initial solution or
378 /// execution then you have to call one \ref init() function and then
379 /// the startFirstPhase() and if you need the startSecondPhase().
383 /// \brief Initializes the internal data structures.
385 /// Initializes the internal data structures.
391 for (NodeIt n(_graph); n != INVALID; ++n) {
395 for (EdgeIt e(_graph); e != INVALID; ++e) {
399 typename Graph::template NodeMap<bool> reached(_graph, false);
402 _level->initAddItem(_target);
404 std::vector<Node> queue;
405 reached.set(_source, true);
407 queue.push_back(_target);
408 reached.set(_target, true);
409 while (!queue.empty()) {
410 _level->initNewLevel();
411 std::vector<Node> nqueue;
412 for (int i = 0; i < int(queue.size()); ++i) {
414 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
415 Node u = _graph.source(e);
416 if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
417 reached.set(u, true);
418 _level->initAddItem(u);
425 _level->initFinish();
427 for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
428 if (_tolerance.positive((*_capacity)[e])) {
429 Node u = _graph.target(e);
430 if ((*_level)[u] == _level->maxLevel()) continue;
431 _flow->set(e, (*_capacity)[e]);
432 _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
433 if (u != _target && !_level->active(u)) {
440 /// \brief Initializes the internal data structures.
442 /// Initializes the internal data structures and sets the initial
443 /// flow to the given \c flowMap. The \c flowMap should contain a
444 /// flow or at least a preflow, ie. in each node excluding the
445 /// target the incoming flow should greater or equal to the
447 /// \return %False when the given \c flowMap is not a preflow.
448 template <typename FlowMap>
449 bool flowInit(const FlowMap& flowMap) {
452 for (EdgeIt e(_graph); e != INVALID; ++e) {
453 _flow->set(e, flowMap[e]);
456 for (NodeIt n(_graph); n != INVALID; ++n) {
458 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
459 excess += (*_flow)[e];
461 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
462 excess -= (*_flow)[e];
464 if (excess < 0 && n != _source) return false;
465 _excess->set(n, excess);
468 typename Graph::template NodeMap<bool> reached(_graph, false);
471 _level->initAddItem(_target);
473 std::vector<Node> queue;
474 reached.set(_source, true);
476 queue.push_back(_target);
477 reached.set(_target, true);
478 while (!queue.empty()) {
479 _level->initNewLevel();
480 std::vector<Node> nqueue;
481 for (int i = 0; i < int(queue.size()); ++i) {
483 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
484 Node u = _graph.source(e);
486 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
487 reached.set(u, true);
488 _level->initAddItem(u);
492 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
493 Node v = _graph.target(e);
494 if (!reached[v] && _tolerance.positive((*_flow)[e])) {
495 reached.set(v, true);
496 _level->initAddItem(v);
503 _level->initFinish();
505 for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
506 Value rem = (*_capacity)[e] - (*_flow)[e];
507 if (_tolerance.positive(rem)) {
508 Node u = _graph.target(e);
509 if ((*_level)[u] == _level->maxLevel()) continue;
510 _flow->set(e, (*_capacity)[e]);
511 _excess->set(u, (*_excess)[u] + rem);
512 if (u != _target && !_level->active(u)) {
517 for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
518 Value rem = (*_flow)[e];
519 if (_tolerance.positive(rem)) {
520 Node v = _graph.source(e);
521 if ((*_level)[v] == _level->maxLevel()) continue;
523 _excess->set(v, (*_excess)[v] + rem);
524 if (v != _target && !_level->active(v)) {
532 /// \brief Starts the first phase of the preflow algorithm.
534 /// The preflow algorithm consists of two phases, this method runs
535 /// the first phase. After the first phase the maximum flow value
536 /// and a minimum value cut can already be computed, although a
537 /// maximum flow is not yet obtained. So after calling this method
538 /// \ref flowValue() returns the value of a maximum flow and \ref
539 /// minCut() returns a minimum cut.
540 /// \pre One of the \ref init() functions should be called.
541 void startFirstPhase() {
544 Node n = _level->highestActive();
545 int level = _level->highestActiveLevel();
546 while (n != INVALID) {
549 while (num > 0 && n != INVALID) {
550 Value excess = (*_excess)[n];
551 int new_level = _level->maxLevel();
553 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
554 Value rem = (*_capacity)[e] - (*_flow)[e];
555 if (!_tolerance.positive(rem)) continue;
556 Node v = _graph.target(e);
557 if ((*_level)[v] < level) {
558 if (!_level->active(v) && v != _target) {
561 if (!_tolerance.less(rem, excess)) {
562 _flow->set(e, (*_flow)[e] + excess);
563 _excess->set(v, (*_excess)[v] + excess);
568 _excess->set(v, (*_excess)[v] + rem);
569 _flow->set(e, (*_capacity)[e]);
571 } else if (new_level > (*_level)[v]) {
572 new_level = (*_level)[v];
576 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
577 Value rem = (*_flow)[e];
578 if (!_tolerance.positive(rem)) continue;
579 Node v = _graph.source(e);
580 if ((*_level)[v] < level) {
581 if (!_level->active(v) && v != _target) {
584 if (!_tolerance.less(rem, excess)) {
585 _flow->set(e, (*_flow)[e] - excess);
586 _excess->set(v, (*_excess)[v] + excess);
591 _excess->set(v, (*_excess)[v] + rem);
594 } else if (new_level > (*_level)[v]) {
595 new_level = (*_level)[v];
601 _excess->set(n, excess);
604 if (new_level + 1 < _level->maxLevel()) {
605 _level->liftHighestActive(new_level + 1);
607 _level->liftHighestActiveToTop();
609 if (_level->emptyLevel(level)) {
610 _level->liftToTop(level);
613 _level->deactivate(n);
616 n = _level->highestActive();
617 level = _level->highestActiveLevel();
621 num = _node_num * 20;
622 while (num > 0 && n != INVALID) {
623 Value excess = (*_excess)[n];
624 int new_level = _level->maxLevel();
626 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
627 Value rem = (*_capacity)[e] - (*_flow)[e];
628 if (!_tolerance.positive(rem)) continue;
629 Node v = _graph.target(e);
630 if ((*_level)[v] < level) {
631 if (!_level->active(v) && v != _target) {
634 if (!_tolerance.less(rem, excess)) {
635 _flow->set(e, (*_flow)[e] + excess);
636 _excess->set(v, (*_excess)[v] + excess);
641 _excess->set(v, (*_excess)[v] + rem);
642 _flow->set(e, (*_capacity)[e]);
644 } else if (new_level > (*_level)[v]) {
645 new_level = (*_level)[v];
649 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
650 Value rem = (*_flow)[e];
651 if (!_tolerance.positive(rem)) continue;
652 Node v = _graph.source(e);
653 if ((*_level)[v] < level) {
654 if (!_level->active(v) && v != _target) {
657 if (!_tolerance.less(rem, excess)) {
658 _flow->set(e, (*_flow)[e] - excess);
659 _excess->set(v, (*_excess)[v] + excess);
664 _excess->set(v, (*_excess)[v] + rem);
667 } else if (new_level > (*_level)[v]) {
668 new_level = (*_level)[v];
674 _excess->set(n, excess);
677 if (new_level + 1 < _level->maxLevel()) {
678 _level->liftActiveOn(level, new_level + 1);
680 _level->liftActiveToTop(level);
682 if (_level->emptyLevel(level)) {
683 _level->liftToTop(level);
686 _level->deactivate(n);
689 while (level >= 0 && _level->activeFree(level)) {
693 n = _level->highestActive();
694 level = _level->highestActiveLevel();
696 n = _level->activeOn(level);
703 /// \brief Starts the second phase of the preflow algorithm.
705 /// The preflow algorithm consists of two phases, this method runs
706 /// the second phase. After calling \ref init() and \ref
707 /// startFirstPhase() and then \ref startSecondPhase(), \ref
708 /// flowMap() return a maximum flow, \ref flowValue() returns the
709 /// value of a maximum flow, \ref minCut() returns a minimum cut
710 /// \pre The \ref init() and startFirstPhase() functions should be
712 void startSecondPhase() {
715 typename Graph::template NodeMap<bool> reached(_graph);
716 for (NodeIt n(_graph); n != INVALID; ++n) {
717 reached.set(n, (*_level)[n] < _level->maxLevel());
721 _level->initAddItem(_source);
723 std::vector<Node> queue;
724 queue.push_back(_source);
725 reached.set(_source, true);
727 while (!queue.empty()) {
728 _level->initNewLevel();
729 std::vector<Node> nqueue;
730 for (int i = 0; i < int(queue.size()); ++i) {
732 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
733 Node v = _graph.target(e);
734 if (!reached[v] && _tolerance.positive((*_flow)[e])) {
735 reached.set(v, true);
736 _level->initAddItem(v);
740 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
741 Node u = _graph.source(e);
743 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
744 reached.set(u, true);
745 _level->initAddItem(u);
752 _level->initFinish();
754 for (NodeIt n(_graph); n != INVALID; ++n) {
756 _level->markToBottom(n);
757 } else if ((*_excess)[n] > 0 && _target != n) {
763 while ((n = _level->highestActive()) != INVALID) {
764 Value excess = (*_excess)[n];
765 int level = _level->highestActiveLevel();
766 int new_level = _level->maxLevel();
768 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
769 Value rem = (*_capacity)[e] - (*_flow)[e];
770 if (!_tolerance.positive(rem)) continue;
771 Node v = _graph.target(e);
772 if ((*_level)[v] < level) {
773 if (!_level->active(v) && v != _source) {
776 if (!_tolerance.less(rem, excess)) {
777 _flow->set(e, (*_flow)[e] + excess);
778 _excess->set(v, (*_excess)[v] + excess);
783 _excess->set(v, (*_excess)[v] + rem);
784 _flow->set(e, (*_capacity)[e]);
786 } else if (new_level > (*_level)[v]) {
787 new_level = (*_level)[v];
791 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
792 Value rem = (*_flow)[e];
793 if (!_tolerance.positive(rem)) continue;
794 Node v = _graph.source(e);
795 if ((*_level)[v] < level) {
796 if (!_level->active(v) && v != _source) {
799 if (!_tolerance.less(rem, excess)) {
800 _flow->set(e, (*_flow)[e] - excess);
801 _excess->set(v, (*_excess)[v] + excess);
806 _excess->set(v, (*_excess)[v] + rem);
809 } else if (new_level > (*_level)[v]) {
810 new_level = (*_level)[v];
816 _excess->set(n, excess);
819 if (new_level + 1 < _level->maxLevel()) {
820 _level->liftHighestActive(new_level + 1);
823 _level->liftHighestActiveToTop();
825 if (_level->emptyLevel(level)) {
827 _level->liftToTop(level);
830 _level->deactivate(n);
836 /// \brief Runs the preflow algorithm.
838 /// Runs the preflow algorithm.
839 /// \note pf.run() is just a shortcut of the following code.
842 /// pf.startFirstPhase();
843 /// pf.startSecondPhase();
851 /// \brief Runs the preflow algorithm to compute the minimum cut.
853 /// Runs the preflow algorithm to compute the minimum cut.
854 /// \note pf.runMinCut() is just a shortcut of the following code.
857 /// pf.startFirstPhase();
866 /// \name Query Functions
867 /// The result of the %Preflow algorithm can be obtained using these
869 /// Before the use of these functions,
870 /// either run() or start() must be called.
874 /// \brief Returns the value of the maximum flow.
876 /// Returns the value of the maximum flow by returning the excess
877 /// of the target node \c t. This value equals to the value of
878 /// the maximum flow already after the first phase.
879 Value flowValue() const {
880 return (*_excess)[_target];
883 /// \brief Returns true when the node is on the source side of minimum cut.
885 /// Returns true when the node is on the source side of minimum
886 /// cut. This method can be called both after running \ref
887 /// startFirstPhase() and \ref startSecondPhase().
888 bool minCut(const Node& node) const {
889 return ((*_level)[node] == _level->maxLevel()) == _phase;
892 /// \brief Returns a minimum value cut.
894 /// Sets the \c cutMap to the characteristic vector of a minimum value
895 /// cut. This method can be called both after running \ref
896 /// startFirstPhase() and \ref startSecondPhase(). The result after second
897 /// phase could be changed slightly if inexact computation is used.
898 /// \pre The \c cutMap should be a bool-valued node-map.
899 template <typename CutMap>
900 void minCutMap(CutMap& cutMap) const {
901 for (NodeIt n(_graph); n != INVALID; ++n) {
902 cutMap.set(n, minCut(n));
906 /// \brief Returns the flow on the edge.
908 /// Sets the \c flowMap to the flow on the edges. This method can
909 /// be called after the second phase of algorithm.
910 Value flow(const Edge& edge) const {
911 return (*_flow)[edge];