Port preflow push max flow alg. from svn -r3516 (#176)
Namely,
- port the files
- apply the migrate script
- apply the unify script
- break the long lines in lemon/preflow.h
- convert the .dim test file to .lgf
- fix compilation problems
1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2008
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/tolerance.h>
24 #include <lemon/elevator.h>
28 /// \brief Implementation of the preflow algorithm.
32 /// \brief Default traits class of Preflow class.
34 /// Default traits class of Preflow class.
35 /// \param _Graph Digraph type.
36 /// \param _CapacityMap Type of capacity map.
37 template <typename _Graph, typename _CapacityMap>
38 struct PreflowDefaultTraits {
40 /// \brief The digraph type the algorithm runs on.
41 typedef _Graph Digraph;
43 /// \brief The type of the map that stores the arc capacities.
45 /// The type of the map that stores the arc capacities.
46 /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
47 typedef _CapacityMap CapacityMap;
49 /// \brief The type of the length of the arcs.
50 typedef typename CapacityMap::Value Value;
52 /// \brief The map type that stores the flow values.
54 /// The map type that stores the flow values.
55 /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
56 typedef typename Digraph::template ArcMap<Value> FlowMap;
58 /// \brief Instantiates a FlowMap.
60 /// This function instantiates a \ref FlowMap.
61 /// \param digraph The digraph, to which we would like to define
63 static FlowMap* createFlowMap(const Digraph& digraph) {
64 return new FlowMap(digraph);
67 /// \brief The eleavator type used by Preflow algorithm.
69 /// The elevator type used by Preflow algorithm.
72 /// \sa LinkedElevator
73 typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;
75 /// \brief Instantiates an Elevator.
77 /// This function instantiates a \ref Elevator.
78 /// \param digraph The digraph, to which we would like to define
80 /// \param max_level The maximum level of the elevator.
81 static Elevator* createElevator(const Digraph& digraph, int max_level) {
82 return new Elevator(digraph, max_level);
85 /// \brief The tolerance used by the algorithm
87 /// The tolerance used by the algorithm to handle inexact computation.
88 typedef lemon::Tolerance<Value> Tolerance;
95 /// \brief %Preflow algorithms class.
97 /// This class provides an implementation of the Goldberg's \e
98 /// preflow \e algorithm producing a flow of maximum value in a
99 /// digraph. The preflow algorithms are the fastest known max
100 /// flow algorithms. The current implementation use a mixture of the
101 /// \e "highest label" and the \e "bound decrease" heuristics.
102 /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
104 /// The algorithm consists from two phases. After the first phase
105 /// the maximal flow value and the minimum cut can be obtained. The
106 /// second phase constructs the feasible maximum flow on each arc.
108 /// \param _Graph The digraph type the algorithm runs on.
109 /// \param _CapacityMap The flow map type.
110 /// \param _Traits Traits class to set various data types used by
111 /// the algorithm. The default traits class is \ref
112 /// PreflowDefaultTraits. See \ref PreflowDefaultTraits for the
113 /// documentation of a %Preflow traits class.
115 ///\author Jacint Szabo and Balazs Dezso
117 template <typename _Graph, typename _CapacityMap, typename _Traits>
119 template <typename _Graph,
120 typename _CapacityMap = typename _Graph::template ArcMap<int>,
121 typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
126 typedef _Traits Traits;
127 typedef typename Traits::Digraph Digraph;
128 typedef typename Traits::CapacityMap CapacityMap;
129 typedef typename Traits::Value Value;
131 typedef typename Traits::FlowMap FlowMap;
132 typedef typename Traits::Elevator Elevator;
133 typedef typename Traits::Tolerance Tolerance;
135 /// \brief \ref Exception for uninitialized parameters.
137 /// This error represents problems in the initialization
138 /// of the parameters of the algorithms.
139 class UninitializedParameter : public lemon::Exception {
141 virtual const char* what() const throw() {
142 return "lemon::Preflow::UninitializedParameter";
148 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
150 const Digraph& _graph;
151 const CapacityMap* _capacity;
155 Node _source, _target;
163 typedef typename Digraph::template NodeMap<Value> ExcessMap;
166 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 Digraph&) {
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<Digraph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
223 typedef Preflow<Digraph, CapacityMap,
224 DefFlowMapTraits<_FlowMap> > Create;
227 template <typename _Elevator>
228 struct DefElevatorTraits : public Traits {
229 typedef _Elevator Elevator;
230 static Elevator *createElevator(const Digraph&, int) {
231 throw UninitializedParameter();
235 /// \brief \ref named-templ-param "Named parameter" for setting
238 /// \ref named-templ-param "Named parameter" for setting Elevator
240 template <typename _Elevator>
242 : public Preflow<Digraph, CapacityMap, DefElevatorTraits<_Elevator> > {
243 typedef Preflow<Digraph, CapacityMap,
244 DefElevatorTraits<_Elevator> > Create;
247 template <typename _Elevator>
248 struct DefStandardElevatorTraits : public Traits {
249 typedef _Elevator Elevator;
250 static Elevator *createElevator(const Digraph& digraph, int max_level) {
251 return new Elevator(digraph, max_level);
255 /// \brief \ref named-templ-param "Named parameter" for setting
258 /// \ref named-templ-param "Named parameter" for setting Elevator
259 /// type. The Elevator should be standard constructor interface, ie.
260 /// the digraph and the maximum level should be passed to it.
261 template <typename _Elevator>
262 struct DefStandardElevator
263 : public Preflow<Digraph, CapacityMap,
264 DefStandardElevatorTraits<_Elevator> > {
265 typedef Preflow<Digraph, CapacityMap,
266 DefStandardElevatorTraits<_Elevator> > Create;
278 /// \brief The constructor of the class.
280 /// The constructor of the class.
281 /// \param digraph The digraph the algorithm runs on.
282 /// \param capacity The capacity of the arcs.
283 /// \param source The source node.
284 /// \param target The target node.
285 Preflow(const Digraph& digraph, const CapacityMap& capacity,
286 Node source, Node target)
287 : _graph(digraph), _capacity(&capacity),
288 _node_num(0), _source(source), _target(target),
289 _flow(0), _local_flow(false),
290 _level(0), _local_level(false),
291 _excess(0), _tolerance(), _phase() {}
293 /// \brief Destrcutor.
300 /// \brief Sets the capacity map.
302 /// Sets the capacity map.
303 /// \return \c (*this)
304 Preflow& capacityMap(const CapacityMap& map) {
309 /// \brief Sets the flow map.
311 /// Sets the flow map.
312 /// \return \c (*this)
313 Preflow& flowMap(FlowMap& map) {
322 /// \brief Returns the flow map.
324 /// \return The flow map.
325 const FlowMap& flowMap() {
329 /// \brief Sets the elevator.
331 /// Sets the elevator.
332 /// \return \c (*this)
333 Preflow& elevator(Elevator& elevator) {
336 _local_level = false;
342 /// \brief Returns the elevator.
344 /// \return The elevator.
345 const Elevator& elevator() {
349 /// \brief Sets the source node.
351 /// Sets the source node.
352 /// \return \c (*this)
353 Preflow& source(const Node& node) {
358 /// \brief Sets the target node.
360 /// Sets the target node.
361 /// \return \c (*this)
362 Preflow& target(const Node& node) {
367 /// \brief Sets the tolerance used by algorithm.
369 /// Sets the tolerance used by algorithm.
370 Preflow& tolerance(const Tolerance& tolerance) const {
371 _tolerance = tolerance;
375 /// \brief Returns the tolerance used by algorithm.
377 /// Returns the tolerance used by algorithm.
378 const Tolerance& tolerance() const {
382 /// \name Execution control The simplest way to execute the
383 /// algorithm is to use one of the member functions called \c
386 /// If you need more control on initial solution or
387 /// execution then you have to call one \ref init() function and then
388 /// the startFirstPhase() and if you need the startSecondPhase().
392 /// \brief Initializes the internal data structures.
394 /// Initializes the internal data structures.
400 for (NodeIt n(_graph); n != INVALID; ++n) {
404 for (ArcIt e(_graph); e != INVALID; ++e) {
408 typename Digraph::template NodeMap<bool> reached(_graph, false);
411 _level->initAddItem(_target);
413 std::vector<Node> queue;
414 reached.set(_source, true);
416 queue.push_back(_target);
417 reached.set(_target, true);
418 while (!queue.empty()) {
419 _level->initNewLevel();
420 std::vector<Node> nqueue;
421 for (int i = 0; i < int(queue.size()); ++i) {
423 for (InArcIt e(_graph, n); e != INVALID; ++e) {
424 Node u = _graph.source(e);
425 if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
426 reached.set(u, true);
427 _level->initAddItem(u);
434 _level->initFinish();
436 for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
437 if (_tolerance.positive((*_capacity)[e])) {
438 Node u = _graph.target(e);
439 if ((*_level)[u] == _level->maxLevel()) continue;
440 _flow->set(e, (*_capacity)[e]);
441 _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
442 if (u != _target && !_level->active(u)) {
449 /// \brief Initializes the internal data structures.
451 /// Initializes the internal data structures and sets the initial
452 /// flow to the given \c flowMap. The \c flowMap should contain a
453 /// flow or at least a preflow, ie. in each node excluding the
454 /// target the incoming flow should greater or equal to the
456 /// \return %False when the given \c flowMap is not a preflow.
457 template <typename FlowMap>
458 bool flowInit(const FlowMap& flowMap) {
461 for (ArcIt e(_graph); e != INVALID; ++e) {
462 _flow->set(e, flowMap[e]);
465 for (NodeIt n(_graph); n != INVALID; ++n) {
467 for (InArcIt e(_graph, n); e != INVALID; ++e) {
468 excess += (*_flow)[e];
470 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
471 excess -= (*_flow)[e];
473 if (excess < 0 && n != _source) return false;
474 _excess->set(n, excess);
477 typename Digraph::template NodeMap<bool> reached(_graph, false);
480 _level->initAddItem(_target);
482 std::vector<Node> queue;
483 reached.set(_source, true);
485 queue.push_back(_target);
486 reached.set(_target, true);
487 while (!queue.empty()) {
488 _level->initNewLevel();
489 std::vector<Node> nqueue;
490 for (int i = 0; i < int(queue.size()); ++i) {
492 for (InArcIt e(_graph, n); e != INVALID; ++e) {
493 Node u = _graph.source(e);
495 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
496 reached.set(u, true);
497 _level->initAddItem(u);
501 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
502 Node v = _graph.target(e);
503 if (!reached[v] && _tolerance.positive((*_flow)[e])) {
504 reached.set(v, true);
505 _level->initAddItem(v);
512 _level->initFinish();
514 for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
515 Value rem = (*_capacity)[e] - (*_flow)[e];
516 if (_tolerance.positive(rem)) {
517 Node u = _graph.target(e);
518 if ((*_level)[u] == _level->maxLevel()) continue;
519 _flow->set(e, (*_capacity)[e]);
520 _excess->set(u, (*_excess)[u] + rem);
521 if (u != _target && !_level->active(u)) {
526 for (InArcIt e(_graph, _source); e != INVALID; ++e) {
527 Value rem = (*_flow)[e];
528 if (_tolerance.positive(rem)) {
529 Node v = _graph.source(e);
530 if ((*_level)[v] == _level->maxLevel()) continue;
532 _excess->set(v, (*_excess)[v] + rem);
533 if (v != _target && !_level->active(v)) {
541 /// \brief Starts the first phase of the preflow algorithm.
543 /// The preflow algorithm consists of two phases, this method runs
544 /// the first phase. After the first phase the maximum flow value
545 /// and a minimum value cut can already be computed, although a
546 /// maximum flow is not yet obtained. So after calling this method
547 /// \ref flowValue() returns the value of a maximum flow and \ref
548 /// minCut() returns a minimum cut.
549 /// \pre One of the \ref init() functions should be called.
550 void startFirstPhase() {
553 Node n = _level->highestActive();
554 int level = _level->highestActiveLevel();
555 while (n != INVALID) {
558 while (num > 0 && n != INVALID) {
559 Value excess = (*_excess)[n];
560 int new_level = _level->maxLevel();
562 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
563 Value rem = (*_capacity)[e] - (*_flow)[e];
564 if (!_tolerance.positive(rem)) continue;
565 Node v = _graph.target(e);
566 if ((*_level)[v] < level) {
567 if (!_level->active(v) && v != _target) {
570 if (!_tolerance.less(rem, excess)) {
571 _flow->set(e, (*_flow)[e] + excess);
572 _excess->set(v, (*_excess)[v] + excess);
577 _excess->set(v, (*_excess)[v] + rem);
578 _flow->set(e, (*_capacity)[e]);
580 } else if (new_level > (*_level)[v]) {
581 new_level = (*_level)[v];
585 for (InArcIt e(_graph, n); e != INVALID; ++e) {
586 Value rem = (*_flow)[e];
587 if (!_tolerance.positive(rem)) continue;
588 Node v = _graph.source(e);
589 if ((*_level)[v] < level) {
590 if (!_level->active(v) && v != _target) {
593 if (!_tolerance.less(rem, excess)) {
594 _flow->set(e, (*_flow)[e] - excess);
595 _excess->set(v, (*_excess)[v] + excess);
600 _excess->set(v, (*_excess)[v] + rem);
603 } else if (new_level > (*_level)[v]) {
604 new_level = (*_level)[v];
610 _excess->set(n, excess);
613 if (new_level + 1 < _level->maxLevel()) {
614 _level->liftHighestActive(new_level + 1);
616 _level->liftHighestActiveToTop();
618 if (_level->emptyLevel(level)) {
619 _level->liftToTop(level);
622 _level->deactivate(n);
625 n = _level->highestActive();
626 level = _level->highestActiveLevel();
630 num = _node_num * 20;
631 while (num > 0 && n != INVALID) {
632 Value excess = (*_excess)[n];
633 int new_level = _level->maxLevel();
635 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
636 Value rem = (*_capacity)[e] - (*_flow)[e];
637 if (!_tolerance.positive(rem)) continue;
638 Node v = _graph.target(e);
639 if ((*_level)[v] < level) {
640 if (!_level->active(v) && v != _target) {
643 if (!_tolerance.less(rem, excess)) {
644 _flow->set(e, (*_flow)[e] + excess);
645 _excess->set(v, (*_excess)[v] + excess);
650 _excess->set(v, (*_excess)[v] + rem);
651 _flow->set(e, (*_capacity)[e]);
653 } else if (new_level > (*_level)[v]) {
654 new_level = (*_level)[v];
658 for (InArcIt e(_graph, n); e != INVALID; ++e) {
659 Value rem = (*_flow)[e];
660 if (!_tolerance.positive(rem)) continue;
661 Node v = _graph.source(e);
662 if ((*_level)[v] < level) {
663 if (!_level->active(v) && v != _target) {
666 if (!_tolerance.less(rem, excess)) {
667 _flow->set(e, (*_flow)[e] - excess);
668 _excess->set(v, (*_excess)[v] + excess);
673 _excess->set(v, (*_excess)[v] + rem);
676 } else if (new_level > (*_level)[v]) {
677 new_level = (*_level)[v];
683 _excess->set(n, excess);
686 if (new_level + 1 < _level->maxLevel()) {
687 _level->liftActiveOn(level, new_level + 1);
689 _level->liftActiveToTop(level);
691 if (_level->emptyLevel(level)) {
692 _level->liftToTop(level);
695 _level->deactivate(n);
698 while (level >= 0 && _level->activeFree(level)) {
702 n = _level->highestActive();
703 level = _level->highestActiveLevel();
705 n = _level->activeOn(level);
712 /// \brief Starts the second phase of the preflow algorithm.
714 /// The preflow algorithm consists of two phases, this method runs
715 /// the second phase. After calling \ref init() and \ref
716 /// startFirstPhase() and then \ref startSecondPhase(), \ref
717 /// flowMap() return a maximum flow, \ref flowValue() returns the
718 /// value of a maximum flow, \ref minCut() returns a minimum cut
719 /// \pre The \ref init() and startFirstPhase() functions should be
721 void startSecondPhase() {
724 typename Digraph::template NodeMap<bool> reached(_graph);
725 for (NodeIt n(_graph); n != INVALID; ++n) {
726 reached.set(n, (*_level)[n] < _level->maxLevel());
730 _level->initAddItem(_source);
732 std::vector<Node> queue;
733 queue.push_back(_source);
734 reached.set(_source, true);
736 while (!queue.empty()) {
737 _level->initNewLevel();
738 std::vector<Node> nqueue;
739 for (int i = 0; i < int(queue.size()); ++i) {
741 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
742 Node v = _graph.target(e);
743 if (!reached[v] && _tolerance.positive((*_flow)[e])) {
744 reached.set(v, true);
745 _level->initAddItem(v);
749 for (InArcIt e(_graph, n); e != INVALID; ++e) {
750 Node u = _graph.source(e);
752 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
753 reached.set(u, true);
754 _level->initAddItem(u);
761 _level->initFinish();
763 for (NodeIt n(_graph); n != INVALID; ++n) {
765 _level->dirtyTopButOne(n);
766 } else if ((*_excess)[n] > 0 && _target != n) {
772 while ((n = _level->highestActive()) != INVALID) {
773 Value excess = (*_excess)[n];
774 int level = _level->highestActiveLevel();
775 int new_level = _level->maxLevel();
777 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
778 Value rem = (*_capacity)[e] - (*_flow)[e];
779 if (!_tolerance.positive(rem)) continue;
780 Node v = _graph.target(e);
781 if ((*_level)[v] < level) {
782 if (!_level->active(v) && v != _source) {
785 if (!_tolerance.less(rem, excess)) {
786 _flow->set(e, (*_flow)[e] + excess);
787 _excess->set(v, (*_excess)[v] + excess);
792 _excess->set(v, (*_excess)[v] + rem);
793 _flow->set(e, (*_capacity)[e]);
795 } else if (new_level > (*_level)[v]) {
796 new_level = (*_level)[v];
800 for (InArcIt e(_graph, n); e != INVALID; ++e) {
801 Value rem = (*_flow)[e];
802 if (!_tolerance.positive(rem)) continue;
803 Node v = _graph.source(e);
804 if ((*_level)[v] < level) {
805 if (!_level->active(v) && v != _source) {
808 if (!_tolerance.less(rem, excess)) {
809 _flow->set(e, (*_flow)[e] - excess);
810 _excess->set(v, (*_excess)[v] + excess);
815 _excess->set(v, (*_excess)[v] + rem);
818 } else if (new_level > (*_level)[v]) {
819 new_level = (*_level)[v];
825 _excess->set(n, excess);
828 if (new_level + 1 < _level->maxLevel()) {
829 _level->liftHighestActive(new_level + 1);
832 _level->liftHighestActiveToTop();
834 if (_level->emptyLevel(level)) {
836 _level->liftToTop(level);
839 _level->deactivate(n);
845 /// \brief Runs the preflow algorithm.
847 /// Runs the preflow algorithm.
848 /// \note pf.run() is just a shortcut of the following code.
851 /// pf.startFirstPhase();
852 /// pf.startSecondPhase();
860 /// \brief Runs the preflow algorithm to compute the minimum cut.
862 /// Runs the preflow algorithm to compute the minimum cut.
863 /// \note pf.runMinCut() is just a shortcut of the following code.
866 /// pf.startFirstPhase();
875 /// \name Query Functions
876 /// The result of the %Preflow algorithm can be obtained using these
878 /// Before the use of these functions,
879 /// either run() or start() must be called.
883 /// \brief Returns the value of the maximum flow.
885 /// Returns the value of the maximum flow by returning the excess
886 /// of the target node \c t. This value equals to the value of
887 /// the maximum flow already after the first phase.
888 Value flowValue() const {
889 return (*_excess)[_target];
892 /// \brief Returns true when the node is on the source side of minimum cut.
894 /// Returns true when the node is on the source side of minimum
895 /// cut. This method can be called both after running \ref
896 /// startFirstPhase() and \ref startSecondPhase().
897 bool minCut(const Node& node) const {
898 return ((*_level)[node] == _level->maxLevel()) == _phase;
901 /// \brief Returns a minimum value cut.
903 /// Sets the \c cutMap to the characteristic vector of a minimum value
904 /// cut. This method can be called both after running \ref
905 /// startFirstPhase() and \ref startSecondPhase(). The result after second
906 /// phase could be changed slightly if inexact computation is used.
907 /// \pre The \c cutMap should be a bool-valued node-map.
908 template <typename CutMap>
909 void minCutMap(CutMap& cutMap) const {
910 for (NodeIt n(_graph); n != INVALID; ++n) {
911 cutMap.set(n, minCut(n));
915 /// \brief Returns the flow on the arc.
917 /// Sets the \c flowMap to the flow on the arcs. This method can
918 /// be called after the second phase of algorithm.
919 Value flow(const Arc& arc) const {
920 return (*_flow)[arc];