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/tolerance.h>
23 #include <lemon/elevator.h>
27 /// \brief Implementation of the preflow algorithm.
31 /// \brief Default traits class of Preflow class.
33 /// Default traits class of Preflow class.
34 /// \param _Graph Digraph type.
35 /// \param _CapacityMap Type of capacity map.
36 template <typename _Graph, typename _CapacityMap>
37 struct PreflowDefaultTraits {
39 /// \brief The digraph type the algorithm runs on.
40 typedef _Graph Digraph;
42 /// \brief The type of the map that stores the arc capacities.
44 /// The type of the map that stores the arc capacities.
45 /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
46 typedef _CapacityMap CapacityMap;
48 /// \brief The type of the length of the arcs.
49 typedef typename CapacityMap::Value Value;
51 /// \brief The map type that stores the flow values.
53 /// The map type that stores the flow values.
54 /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
55 typedef typename Digraph::template ArcMap<Value> FlowMap;
57 /// \brief Instantiates a FlowMap.
59 /// This function instantiates a \ref FlowMap.
60 /// \param digraph The digraph, to which we would like to define
62 static FlowMap* createFlowMap(const Digraph& digraph) {
63 return new FlowMap(digraph);
66 /// \brief The eleavator type used by Preflow algorithm.
68 /// The elevator type used by Preflow algorithm.
71 /// \sa LinkedElevator
72 typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;
74 /// \brief Instantiates an Elevator.
76 /// This function instantiates a \ref Elevator.
77 /// \param digraph The digraph, to which we would like to define
79 /// \param max_level The maximum level of the elevator.
80 static Elevator* createElevator(const Digraph& digraph, int max_level) {
81 return new Elevator(digraph, max_level);
84 /// \brief The tolerance used by the algorithm
86 /// The tolerance used by the algorithm to handle inexact computation.
87 typedef lemon::Tolerance<Value> Tolerance;
94 /// \brief %Preflow algorithms class.
96 /// This class provides an implementation of the Goldberg's \e
97 /// preflow \e algorithm producing a flow of maximum value in a
98 /// digraph. The preflow algorithms are the fastest known max
99 /// flow algorithms. The current implementation use a mixture of the
100 /// \e "highest label" and the \e "bound decrease" heuristics.
101 /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
103 /// The algorithm consists from two phases. After the first phase
104 /// the maximal flow value and the minimum cut can be obtained. The
105 /// second phase constructs the feasible maximum flow on each arc.
107 /// \param _Graph The digraph type the algorithm runs on.
108 /// \param _CapacityMap The flow map type.
109 /// \param _Traits Traits class to set various data types used by
110 /// the algorithm. The default traits class is \ref
111 /// PreflowDefaultTraits. See \ref PreflowDefaultTraits for the
112 /// documentation of a %Preflow traits class.
114 ///\author Jacint Szabo and Balazs Dezso
116 template <typename _Graph, typename _CapacityMap, typename _Traits>
118 template <typename _Graph,
119 typename _CapacityMap = typename _Graph::template ArcMap<int>,
120 typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
125 typedef _Traits Traits;
126 typedef typename Traits::Digraph Digraph;
127 typedef typename Traits::CapacityMap CapacityMap;
128 typedef typename Traits::Value Value;
130 typedef typename Traits::FlowMap FlowMap;
131 typedef typename Traits::Elevator Elevator;
132 typedef typename Traits::Tolerance Tolerance;
136 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
138 const Digraph& _graph;
139 const CapacityMap* _capacity;
143 Node _source, _target;
151 typedef typename Digraph::template NodeMap<Value> ExcessMap;
154 Tolerance _tolerance;
159 void createStructures() {
160 _node_num = countNodes(_graph);
163 _flow = Traits::createFlowMap(_graph);
167 _level = Traits::createElevator(_graph, _node_num);
171 _excess = new ExcessMap(_graph);
175 void destroyStructures() {
189 typedef Preflow Create;
191 ///\name Named template parameters
195 template <typename _FlowMap>
196 struct SetFlowMapTraits : public Traits {
197 typedef _FlowMap FlowMap;
198 static FlowMap *createFlowMap(const Digraph&) {
199 LEMON_ASSERT(false, "FlowMap is not initialized");
200 return 0; // ignore warnings
204 /// \brief \ref named-templ-param "Named parameter" for setting
207 /// \ref named-templ-param "Named parameter" for setting FlowMap
209 template <typename _FlowMap>
211 : public Preflow<Digraph, CapacityMap, SetFlowMapTraits<_FlowMap> > {
212 typedef Preflow<Digraph, CapacityMap,
213 SetFlowMapTraits<_FlowMap> > Create;
216 template <typename _Elevator>
217 struct SetElevatorTraits : public Traits {
218 typedef _Elevator Elevator;
219 static Elevator *createElevator(const Digraph&, int) {
220 LEMON_ASSERT(false, "Elevator is not initialized");
221 return 0; // ignore warnings
225 /// \brief \ref named-templ-param "Named parameter" for setting
228 /// \ref named-templ-param "Named parameter" for setting Elevator
230 template <typename _Elevator>
232 : public Preflow<Digraph, CapacityMap, SetElevatorTraits<_Elevator> > {
233 typedef Preflow<Digraph, CapacityMap,
234 SetElevatorTraits<_Elevator> > Create;
237 template <typename _Elevator>
238 struct SetStandardElevatorTraits : public Traits {
239 typedef _Elevator Elevator;
240 static Elevator *createElevator(const Digraph& digraph, int max_level) {
241 return new Elevator(digraph, max_level);
245 /// \brief \ref named-templ-param "Named parameter" for setting
248 /// \ref named-templ-param "Named parameter" for setting Elevator
249 /// type. The Elevator should be standard constructor interface, ie.
250 /// the digraph and the maximum level should be passed to it.
251 template <typename _Elevator>
252 struct SetStandardElevator
253 : public Preflow<Digraph, CapacityMap,
254 SetStandardElevatorTraits<_Elevator> > {
255 typedef Preflow<Digraph, CapacityMap,
256 SetStandardElevatorTraits<_Elevator> > Create;
268 /// \brief The constructor of the class.
270 /// The constructor of the class.
271 /// \param digraph The digraph the algorithm runs on.
272 /// \param capacity The capacity of the arcs.
273 /// \param source The source node.
274 /// \param target The target node.
275 Preflow(const Digraph& digraph, const CapacityMap& capacity,
276 Node source, Node target)
277 : _graph(digraph), _capacity(&capacity),
278 _node_num(0), _source(source), _target(target),
279 _flow(0), _local_flow(false),
280 _level(0), _local_level(false),
281 _excess(0), _tolerance(), _phase() {}
283 /// \brief Destrcutor.
290 /// \brief Sets the capacity map.
292 /// Sets the capacity map.
293 /// \return \c (*this)
294 Preflow& capacityMap(const CapacityMap& map) {
299 /// \brief Sets the flow map.
301 /// Sets the flow map.
302 /// \return \c (*this)
303 Preflow& flowMap(FlowMap& map) {
312 /// \brief Returns the flow map.
314 /// \return The flow map.
315 const FlowMap& flowMap() {
319 /// \brief Sets the elevator.
321 /// Sets the elevator.
322 /// \return \c (*this)
323 Preflow& elevator(Elevator& elevator) {
326 _local_level = false;
332 /// \brief Returns the elevator.
334 /// \return The elevator.
335 const Elevator& elevator() {
339 /// \brief Sets the source node.
341 /// Sets the source node.
342 /// \return \c (*this)
343 Preflow& source(const Node& node) {
348 /// \brief Sets the target node.
350 /// Sets the target node.
351 /// \return \c (*this)
352 Preflow& target(const Node& node) {
357 /// \brief Sets the tolerance used by algorithm.
359 /// Sets the tolerance used by algorithm.
360 Preflow& tolerance(const Tolerance& tolerance) const {
361 _tolerance = tolerance;
365 /// \brief Returns the tolerance used by algorithm.
367 /// Returns the tolerance used by algorithm.
368 const Tolerance& tolerance() const {
372 /// \name Execution control The simplest way to execute the
373 /// algorithm is to use one of the member functions called \c
376 /// If you need more control on initial solution or
377 /// execution then you have to call one \ref init() function and then
378 /// the startFirstPhase() and if you need the startSecondPhase().
382 /// \brief Initializes the internal data structures.
384 /// Initializes the internal data structures.
390 for (NodeIt n(_graph); n != INVALID; ++n) {
394 for (ArcIt e(_graph); e != INVALID; ++e) {
398 typename Digraph::template NodeMap<bool> reached(_graph, false);
401 _level->initAddItem(_target);
403 std::vector<Node> queue;
404 reached.set(_source, true);
406 queue.push_back(_target);
407 reached.set(_target, true);
408 while (!queue.empty()) {
409 _level->initNewLevel();
410 std::vector<Node> nqueue;
411 for (int i = 0; i < int(queue.size()); ++i) {
413 for (InArcIt e(_graph, n); e != INVALID; ++e) {
414 Node u = _graph.source(e);
415 if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
416 reached.set(u, true);
417 _level->initAddItem(u);
424 _level->initFinish();
426 for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
427 if (_tolerance.positive((*_capacity)[e])) {
428 Node u = _graph.target(e);
429 if ((*_level)[u] == _level->maxLevel()) continue;
430 _flow->set(e, (*_capacity)[e]);
431 _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
432 if (u != _target && !_level->active(u)) {
439 /// \brief Initializes the internal data structures.
441 /// Initializes the internal data structures and sets the initial
442 /// flow to the given \c flowMap. The \c flowMap should contain a
443 /// flow or at least a preflow, ie. in each node excluding the
444 /// target the incoming flow should greater or equal to the
446 /// \return %False when the given \c flowMap is not a preflow.
447 template <typename FlowMap>
448 bool flowInit(const FlowMap& flowMap) {
451 for (ArcIt e(_graph); e != INVALID; ++e) {
452 _flow->set(e, flowMap[e]);
455 for (NodeIt n(_graph); n != INVALID; ++n) {
457 for (InArcIt e(_graph, n); e != INVALID; ++e) {
458 excess += (*_flow)[e];
460 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
461 excess -= (*_flow)[e];
463 if (excess < 0 && n != _source) return false;
464 _excess->set(n, excess);
467 typename Digraph::template NodeMap<bool> reached(_graph, false);
470 _level->initAddItem(_target);
472 std::vector<Node> queue;
473 reached.set(_source, true);
475 queue.push_back(_target);
476 reached.set(_target, true);
477 while (!queue.empty()) {
478 _level->initNewLevel();
479 std::vector<Node> nqueue;
480 for (int i = 0; i < int(queue.size()); ++i) {
482 for (InArcIt e(_graph, n); e != INVALID; ++e) {
483 Node u = _graph.source(e);
485 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
486 reached.set(u, true);
487 _level->initAddItem(u);
491 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
492 Node v = _graph.target(e);
493 if (!reached[v] && _tolerance.positive((*_flow)[e])) {
494 reached.set(v, true);
495 _level->initAddItem(v);
502 _level->initFinish();
504 for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
505 Value rem = (*_capacity)[e] - (*_flow)[e];
506 if (_tolerance.positive(rem)) {
507 Node u = _graph.target(e);
508 if ((*_level)[u] == _level->maxLevel()) continue;
509 _flow->set(e, (*_capacity)[e]);
510 _excess->set(u, (*_excess)[u] + rem);
511 if (u != _target && !_level->active(u)) {
516 for (InArcIt e(_graph, _source); e != INVALID; ++e) {
517 Value rem = (*_flow)[e];
518 if (_tolerance.positive(rem)) {
519 Node v = _graph.source(e);
520 if ((*_level)[v] == _level->maxLevel()) continue;
522 _excess->set(v, (*_excess)[v] + rem);
523 if (v != _target && !_level->active(v)) {
531 /// \brief Starts the first phase of the preflow algorithm.
533 /// The preflow algorithm consists of two phases, this method runs
534 /// the first phase. After the first phase the maximum flow value
535 /// and a minimum value cut can already be computed, although a
536 /// maximum flow is not yet obtained. So after calling this method
537 /// \ref flowValue() returns the value of a maximum flow and \ref
538 /// minCut() returns a minimum cut.
539 /// \pre One of the \ref init() functions should be called.
540 void startFirstPhase() {
543 Node n = _level->highestActive();
544 int level = _level->highestActiveLevel();
545 while (n != INVALID) {
548 while (num > 0 && n != INVALID) {
549 Value excess = (*_excess)[n];
550 int new_level = _level->maxLevel();
552 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
553 Value rem = (*_capacity)[e] - (*_flow)[e];
554 if (!_tolerance.positive(rem)) continue;
555 Node v = _graph.target(e);
556 if ((*_level)[v] < level) {
557 if (!_level->active(v) && v != _target) {
560 if (!_tolerance.less(rem, excess)) {
561 _flow->set(e, (*_flow)[e] + excess);
562 _excess->set(v, (*_excess)[v] + excess);
567 _excess->set(v, (*_excess)[v] + rem);
568 _flow->set(e, (*_capacity)[e]);
570 } else if (new_level > (*_level)[v]) {
571 new_level = (*_level)[v];
575 for (InArcIt e(_graph, n); e != INVALID; ++e) {
576 Value rem = (*_flow)[e];
577 if (!_tolerance.positive(rem)) continue;
578 Node v = _graph.source(e);
579 if ((*_level)[v] < level) {
580 if (!_level->active(v) && v != _target) {
583 if (!_tolerance.less(rem, excess)) {
584 _flow->set(e, (*_flow)[e] - excess);
585 _excess->set(v, (*_excess)[v] + excess);
590 _excess->set(v, (*_excess)[v] + rem);
593 } else if (new_level > (*_level)[v]) {
594 new_level = (*_level)[v];
600 _excess->set(n, excess);
603 if (new_level + 1 < _level->maxLevel()) {
604 _level->liftHighestActive(new_level + 1);
606 _level->liftHighestActiveToTop();
608 if (_level->emptyLevel(level)) {
609 _level->liftToTop(level);
612 _level->deactivate(n);
615 n = _level->highestActive();
616 level = _level->highestActiveLevel();
620 num = _node_num * 20;
621 while (num > 0 && n != INVALID) {
622 Value excess = (*_excess)[n];
623 int new_level = _level->maxLevel();
625 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
626 Value rem = (*_capacity)[e] - (*_flow)[e];
627 if (!_tolerance.positive(rem)) continue;
628 Node v = _graph.target(e);
629 if ((*_level)[v] < level) {
630 if (!_level->active(v) && v != _target) {
633 if (!_tolerance.less(rem, excess)) {
634 _flow->set(e, (*_flow)[e] + excess);
635 _excess->set(v, (*_excess)[v] + excess);
640 _excess->set(v, (*_excess)[v] + rem);
641 _flow->set(e, (*_capacity)[e]);
643 } else if (new_level > (*_level)[v]) {
644 new_level = (*_level)[v];
648 for (InArcIt e(_graph, n); e != INVALID; ++e) {
649 Value rem = (*_flow)[e];
650 if (!_tolerance.positive(rem)) continue;
651 Node v = _graph.source(e);
652 if ((*_level)[v] < level) {
653 if (!_level->active(v) && v != _target) {
656 if (!_tolerance.less(rem, excess)) {
657 _flow->set(e, (*_flow)[e] - excess);
658 _excess->set(v, (*_excess)[v] + excess);
663 _excess->set(v, (*_excess)[v] + rem);
666 } else if (new_level > (*_level)[v]) {
667 new_level = (*_level)[v];
673 _excess->set(n, excess);
676 if (new_level + 1 < _level->maxLevel()) {
677 _level->liftActiveOn(level, new_level + 1);
679 _level->liftActiveToTop(level);
681 if (_level->emptyLevel(level)) {
682 _level->liftToTop(level);
685 _level->deactivate(n);
688 while (level >= 0 && _level->activeFree(level)) {
692 n = _level->highestActive();
693 level = _level->highestActiveLevel();
695 n = _level->activeOn(level);
702 /// \brief Starts the second phase of the preflow algorithm.
704 /// The preflow algorithm consists of two phases, this method runs
705 /// the second phase. After calling \ref init() and \ref
706 /// startFirstPhase() and then \ref startSecondPhase(), \ref
707 /// flowMap() return a maximum flow, \ref flowValue() returns the
708 /// value of a maximum flow, \ref minCut() returns a minimum cut
709 /// \pre The \ref init() and startFirstPhase() functions should be
711 void startSecondPhase() {
714 typename Digraph::template NodeMap<bool> reached(_graph);
715 for (NodeIt n(_graph); n != INVALID; ++n) {
716 reached.set(n, (*_level)[n] < _level->maxLevel());
720 _level->initAddItem(_source);
722 std::vector<Node> queue;
723 queue.push_back(_source);
724 reached.set(_source, true);
726 while (!queue.empty()) {
727 _level->initNewLevel();
728 std::vector<Node> nqueue;
729 for (int i = 0; i < int(queue.size()); ++i) {
731 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
732 Node v = _graph.target(e);
733 if (!reached[v] && _tolerance.positive((*_flow)[e])) {
734 reached.set(v, true);
735 _level->initAddItem(v);
739 for (InArcIt e(_graph, n); e != INVALID; ++e) {
740 Node u = _graph.source(e);
742 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
743 reached.set(u, true);
744 _level->initAddItem(u);
751 _level->initFinish();
753 for (NodeIt n(_graph); n != INVALID; ++n) {
755 _level->dirtyTopButOne(n);
756 } else if ((*_excess)[n] > 0 && _target != n) {
762 while ((n = _level->highestActive()) != INVALID) {
763 Value excess = (*_excess)[n];
764 int level = _level->highestActiveLevel();
765 int new_level = _level->maxLevel();
767 for (OutArcIt e(_graph, n); e != INVALID; ++e) {
768 Value rem = (*_capacity)[e] - (*_flow)[e];
769 if (!_tolerance.positive(rem)) continue;
770 Node v = _graph.target(e);
771 if ((*_level)[v] < level) {
772 if (!_level->active(v) && v != _source) {
775 if (!_tolerance.less(rem, excess)) {
776 _flow->set(e, (*_flow)[e] + excess);
777 _excess->set(v, (*_excess)[v] + excess);
782 _excess->set(v, (*_excess)[v] + rem);
783 _flow->set(e, (*_capacity)[e]);
785 } else if (new_level > (*_level)[v]) {
786 new_level = (*_level)[v];
790 for (InArcIt e(_graph, n); e != INVALID; ++e) {
791 Value rem = (*_flow)[e];
792 if (!_tolerance.positive(rem)) continue;
793 Node v = _graph.source(e);
794 if ((*_level)[v] < level) {
795 if (!_level->active(v) && v != _source) {
798 if (!_tolerance.less(rem, excess)) {
799 _flow->set(e, (*_flow)[e] - excess);
800 _excess->set(v, (*_excess)[v] + excess);
805 _excess->set(v, (*_excess)[v] + rem);
808 } else if (new_level > (*_level)[v]) {
809 new_level = (*_level)[v];
815 _excess->set(n, excess);
818 if (new_level + 1 < _level->maxLevel()) {
819 _level->liftHighestActive(new_level + 1);
822 _level->liftHighestActiveToTop();
824 if (_level->emptyLevel(level)) {
826 _level->liftToTop(level);
829 _level->deactivate(n);
835 /// \brief Runs the preflow algorithm.
837 /// Runs the preflow algorithm.
838 /// \note pf.run() is just a shortcut of the following code.
841 /// pf.startFirstPhase();
842 /// pf.startSecondPhase();
850 /// \brief Runs the preflow algorithm to compute the minimum cut.
852 /// Runs the preflow algorithm to compute the minimum cut.
853 /// \note pf.runMinCut() is just a shortcut of the following code.
856 /// pf.startFirstPhase();
865 /// \name Query Functions
866 /// The result of the %Preflow algorithm can be obtained using these
868 /// Before the use of these functions,
869 /// either run() or start() must be called.
873 /// \brief Returns the value of the maximum flow.
875 /// Returns the value of the maximum flow by returning the excess
876 /// of the target node \c t. This value equals to the value of
877 /// the maximum flow already after the first phase.
878 Value flowValue() const {
879 return (*_excess)[_target];
882 /// \brief Returns true when the node is on the source side of minimum cut.
884 /// Returns true when the node is on the source side of minimum
885 /// cut. This method can be called both after running \ref
886 /// startFirstPhase() and \ref startSecondPhase().
887 bool minCut(const Node& node) const {
888 return ((*_level)[node] == _level->maxLevel()) == _phase;
891 /// \brief Returns a minimum value cut.
893 /// Sets the \c cutMap to the characteristic vector of a minimum value
894 /// cut. This method can be called both after running \ref
895 /// startFirstPhase() and \ref startSecondPhase(). The result after second
896 /// phase could be changed slightly if inexact computation is used.
897 /// \pre The \c cutMap should be a bool-valued node-map.
898 template <typename CutMap>
899 void minCutMap(CutMap& cutMap) const {
900 for (NodeIt n(_graph); n != INVALID; ++n) {
901 cutMap.set(n, minCut(n));
905 /// \brief Returns the flow on the arc.
907 /// Sets the \c flowMap to the flow on the arcs. This method can
908 /// be called after the second phase of algorithm.
909 Value flow(const Arc& arc) const {
910 return (*_flow)[arc];