Revert to long long int since currently I don't know a better solution.
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/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;
172 void createStructures() {
173 _node_num = countNodes(_graph);
176 _flow = Traits::createFlowMap(_graph);
180 _level = Traits::createElevator(_graph, _node_num);
184 _excess = new ExcessMap(_graph);
188 void destroyStructures() {
202 typedef Preflow Create;
204 ///\name Named template parameters
208 template <typename _FlowMap>
209 struct DefFlowMapTraits : public Traits {
210 typedef _FlowMap FlowMap;
211 static FlowMap *createFlowMap(const Graph&) {
212 throw UninitializedParameter();
216 /// \brief \ref named-templ-param "Named parameter" for setting
219 /// \ref named-templ-param "Named parameter" for setting FlowMap
221 template <typename _FlowMap>
223 : public Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
224 typedef Preflow<Graph, CapacityMap, DefFlowMapTraits<_FlowMap> > Create;
227 template <typename _Elevator>
228 struct DefElevatorTraits : public Traits {
229 typedef _Elevator Elevator;
230 static Elevator *createElevator(const Graph&, 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<Graph, CapacityMap, DefElevatorTraits<_Elevator> > {
243 typedef Preflow<Graph, CapacityMap, DefElevatorTraits<_Elevator> > Create;
246 template <typename _Elevator>
247 struct DefStandardElevatorTraits : public Traits {
248 typedef _Elevator Elevator;
249 static Elevator *createElevator(const Graph& graph, int max_level) {
250 return new Elevator(graph, max_level);
254 /// \brief \ref named-templ-param "Named parameter" for setting
257 /// \ref named-templ-param "Named parameter" for setting Elevator
258 /// type. The Elevator should be standard constructor interface, ie.
259 /// the graph and the maximum level should be passed to it.
260 template <typename _Elevator>
261 struct DefStandardElevator
262 : public Preflow<Graph, CapacityMap,
263 DefStandardElevatorTraits<_Elevator> > {
264 typedef Preflow<Graph, CapacityMap,
265 DefStandardElevatorTraits<_Elevator> > Create;
277 /// \brief The constructor of the class.
279 /// The constructor of the class.
280 /// \param graph The directed graph the algorithm runs on.
281 /// \param capacity The capacity of the edges.
282 /// \param source The source node.
283 /// \param target The target node.
284 Preflow(const Graph& graph, const CapacityMap& capacity,
285 Node source, Node target)
286 : _graph(graph), _capacity(&capacity),
287 _node_num(0), _source(source), _target(target),
288 _flow(0), _local_flow(false),
289 _level(0), _local_level(false),
290 _excess(0), _tolerance(), _phase() {}
292 /// \brief Destrcutor.
299 /// \brief Sets the capacity map.
301 /// Sets the capacity map.
302 /// \return \c (*this)
303 Preflow& capacityMap(const CapacityMap& map) {
308 /// \brief Sets the flow map.
310 /// Sets the flow map.
311 /// \return \c (*this)
312 Preflow& flowMap(FlowMap& map) {
321 /// \brief Returns the flow map.
323 /// \return The flow map.
324 const FlowMap& flowMap() {
328 /// \brief Sets the elevator.
330 /// Sets the elevator.
331 /// \return \c (*this)
332 Preflow& elevator(Elevator& elevator) {
335 _local_level = false;
341 /// \brief Returns the elevator.
343 /// \return The elevator.
344 const Elevator& elevator() {
348 /// \brief Sets the source node.
350 /// Sets the source node.
351 /// \return \c (*this)
352 Preflow& source(const Node& node) {
357 /// \brief Sets the target node.
359 /// Sets the target node.
360 /// \return \c (*this)
361 Preflow& target(const Node& node) {
366 /// \brief Sets the tolerance used by algorithm.
368 /// Sets the tolerance used by algorithm.
369 Preflow& tolerance(const Tolerance& tolerance) const {
370 _tolerance = tolerance;
374 /// \brief Returns the tolerance used by algorithm.
376 /// Returns the tolerance used by algorithm.
377 const Tolerance& tolerance() const {
381 /// \name Execution control The simplest way to execute the
382 /// algorithm is to use one of the member functions called \c
385 /// If you need more control on initial solution or
386 /// execution then you have to call one \ref init() function and then
387 /// the startFirstPhase() and if you need the startSecondPhase().
391 /// \brief Initializes the internal data structures.
393 /// Initializes the internal data structures.
399 for (NodeIt n(_graph); n != INVALID; ++n) {
403 for (EdgeIt e(_graph); e != INVALID; ++e) {
407 typename Graph::template NodeMap<bool> reached(_graph, false);
410 _level->initAddItem(_target);
412 std::vector<Node> queue;
413 reached.set(_source, true);
415 queue.push_back(_target);
416 reached.set(_target, true);
417 while (!queue.empty()) {
418 _level->initNewLevel();
419 std::vector<Node> nqueue;
420 for (int i = 0; i < int(queue.size()); ++i) {
422 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
423 Node u = _graph.source(e);
424 if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
425 reached.set(u, true);
426 _level->initAddItem(u);
433 _level->initFinish();
435 for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
436 if (_tolerance.positive((*_capacity)[e])) {
437 Node u = _graph.target(e);
438 if ((*_level)[u] == _level->maxLevel()) continue;
439 _flow->set(e, (*_capacity)[e]);
440 _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
441 if (u != _target && !_level->active(u)) {
448 /// \brief Initializes the internal data structures.
450 /// Initializes the internal data structures and sets the initial
451 /// flow to the given \c flowMap. The \c flowMap should contain a
452 /// flow or at least a preflow, ie. in each node excluding the
453 /// target the incoming flow should greater or equal to the
455 /// \return %False when the given \c flowMap is not a preflow.
456 template <typename FlowMap>
457 bool flowInit(const FlowMap& flowMap) {
460 for (EdgeIt e(_graph); e != INVALID; ++e) {
461 _flow->set(e, flowMap[e]);
464 for (NodeIt n(_graph); n != INVALID; ++n) {
466 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
467 excess += (*_flow)[e];
469 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
470 excess -= (*_flow)[e];
472 if (excess < 0 && n != _source) return false;
473 _excess->set(n, excess);
476 typename Graph::template NodeMap<bool> reached(_graph, false);
479 _level->initAddItem(_target);
481 std::vector<Node> queue;
482 reached.set(_source, true);
484 queue.push_back(_target);
485 reached.set(_target, true);
486 while (!queue.empty()) {
487 _level->initNewLevel();
488 std::vector<Node> nqueue;
489 for (int i = 0; i < int(queue.size()); ++i) {
491 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
492 Node u = _graph.source(e);
494 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
495 reached.set(u, true);
496 _level->initAddItem(u);
500 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
501 Node v = _graph.target(e);
502 if (!reached[v] && _tolerance.positive((*_flow)[e])) {
503 reached.set(v, true);
504 _level->initAddItem(v);
511 _level->initFinish();
513 for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
514 Value rem = (*_capacity)[e] - (*_flow)[e];
515 if (_tolerance.positive(rem)) {
516 Node u = _graph.target(e);
517 if ((*_level)[u] == _level->maxLevel()) continue;
518 _flow->set(e, (*_capacity)[e]);
519 _excess->set(u, (*_excess)[u] + rem);
520 if (u != _target && !_level->active(u)) {
525 for (InEdgeIt e(_graph, _source); e != INVALID; ++e) {
526 Value rem = (*_flow)[e];
527 if (_tolerance.positive(rem)) {
528 Node v = _graph.source(e);
529 if ((*_level)[v] == _level->maxLevel()) continue;
531 _excess->set(v, (*_excess)[v] + rem);
532 if (v != _target && !_level->active(v)) {
540 /// \brief Starts the first phase of the preflow algorithm.
542 /// The preflow algorithm consists of two phases, this method runs
543 /// the first phase. After the first phase the maximum flow value
544 /// and a minimum value cut can already be computed, although a
545 /// maximum flow is not yet obtained. So after calling this method
546 /// \ref flowValue() returns the value of a maximum flow and \ref
547 /// minCut() returns a minimum cut.
548 /// \pre One of the \ref init() functions should be called.
549 void startFirstPhase() {
552 Node n = _level->highestActive();
553 int level = _level->highestActiveLevel();
554 while (n != INVALID) {
557 while (num > 0 && n != INVALID) {
558 Value excess = (*_excess)[n];
559 int new_level = _level->maxLevel();
561 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
562 Value rem = (*_capacity)[e] - (*_flow)[e];
563 if (!_tolerance.positive(rem)) continue;
564 Node v = _graph.target(e);
565 if ((*_level)[v] < level) {
566 if (!_level->active(v) && v != _target) {
569 if (!_tolerance.less(rem, excess)) {
570 _flow->set(e, (*_flow)[e] + excess);
571 _excess->set(v, (*_excess)[v] + excess);
576 _excess->set(v, (*_excess)[v] + rem);
577 _flow->set(e, (*_capacity)[e]);
579 } else if (new_level > (*_level)[v]) {
580 new_level = (*_level)[v];
584 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
585 Value rem = (*_flow)[e];
586 if (!_tolerance.positive(rem)) continue;
587 Node v = _graph.source(e);
588 if ((*_level)[v] < level) {
589 if (!_level->active(v) && v != _target) {
592 if (!_tolerance.less(rem, excess)) {
593 _flow->set(e, (*_flow)[e] - excess);
594 _excess->set(v, (*_excess)[v] + excess);
599 _excess->set(v, (*_excess)[v] + rem);
602 } else if (new_level > (*_level)[v]) {
603 new_level = (*_level)[v];
609 _excess->set(n, excess);
612 if (new_level + 1 < _level->maxLevel()) {
613 _level->liftHighestActive(new_level + 1);
615 _level->liftHighestActiveToTop();
617 if (_level->emptyLevel(level)) {
618 _level->liftToTop(level);
621 _level->deactivate(n);
624 n = _level->highestActive();
625 level = _level->highestActiveLevel();
629 num = _node_num * 20;
630 while (num > 0 && n != INVALID) {
631 Value excess = (*_excess)[n];
632 int new_level = _level->maxLevel();
634 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
635 Value rem = (*_capacity)[e] - (*_flow)[e];
636 if (!_tolerance.positive(rem)) continue;
637 Node v = _graph.target(e);
638 if ((*_level)[v] < level) {
639 if (!_level->active(v) && v != _target) {
642 if (!_tolerance.less(rem, excess)) {
643 _flow->set(e, (*_flow)[e] + excess);
644 _excess->set(v, (*_excess)[v] + excess);
649 _excess->set(v, (*_excess)[v] + rem);
650 _flow->set(e, (*_capacity)[e]);
652 } else if (new_level > (*_level)[v]) {
653 new_level = (*_level)[v];
657 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
658 Value rem = (*_flow)[e];
659 if (!_tolerance.positive(rem)) continue;
660 Node v = _graph.source(e);
661 if ((*_level)[v] < level) {
662 if (!_level->active(v) && v != _target) {
665 if (!_tolerance.less(rem, excess)) {
666 _flow->set(e, (*_flow)[e] - excess);
667 _excess->set(v, (*_excess)[v] + excess);
672 _excess->set(v, (*_excess)[v] + rem);
675 } else if (new_level > (*_level)[v]) {
676 new_level = (*_level)[v];
682 _excess->set(n, excess);
685 if (new_level + 1 < _level->maxLevel()) {
686 _level->liftActiveOn(level, new_level + 1);
688 _level->liftActiveToTop(level);
690 if (_level->emptyLevel(level)) {
691 _level->liftToTop(level);
694 _level->deactivate(n);
697 while (level >= 0 && _level->activeFree(level)) {
701 n = _level->highestActive();
702 level = _level->highestActiveLevel();
704 n = _level->activeOn(level);
711 /// \brief Starts the second phase of the preflow algorithm.
713 /// The preflow algorithm consists of two phases, this method runs
714 /// the second phase. After calling \ref init() and \ref
715 /// startFirstPhase() and then \ref startSecondPhase(), \ref
716 /// flowMap() return a maximum flow, \ref flowValue() returns the
717 /// value of a maximum flow, \ref minCut() returns a minimum cut
718 /// \pre The \ref init() and startFirstPhase() functions should be
720 void startSecondPhase() {
723 typename Graph::template NodeMap<bool> reached(_graph);
724 for (NodeIt n(_graph); n != INVALID; ++n) {
725 reached.set(n, (*_level)[n] < _level->maxLevel());
729 _level->initAddItem(_source);
731 std::vector<Node> queue;
732 queue.push_back(_source);
733 reached.set(_source, true);
735 while (!queue.empty()) {
736 _level->initNewLevel();
737 std::vector<Node> nqueue;
738 for (int i = 0; i < int(queue.size()); ++i) {
740 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
741 Node v = _graph.target(e);
742 if (!reached[v] && _tolerance.positive((*_flow)[e])) {
743 reached.set(v, true);
744 _level->initAddItem(v);
748 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
749 Node u = _graph.source(e);
751 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
752 reached.set(u, true);
753 _level->initAddItem(u);
760 _level->initFinish();
762 for (NodeIt n(_graph); n != INVALID; ++n) {
764 _level->markToBottom(n);
765 } else if ((*_excess)[n] > 0 && _target != n) {
771 while ((n = _level->highestActive()) != INVALID) {
772 Value excess = (*_excess)[n];
773 int level = _level->highestActiveLevel();
774 int new_level = _level->maxLevel();
776 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
777 Value rem = (*_capacity)[e] - (*_flow)[e];
778 if (!_tolerance.positive(rem)) continue;
779 Node v = _graph.target(e);
780 if ((*_level)[v] < level) {
781 if (!_level->active(v) && v != _source) {
784 if (!_tolerance.less(rem, excess)) {
785 _flow->set(e, (*_flow)[e] + excess);
786 _excess->set(v, (*_excess)[v] + excess);
791 _excess->set(v, (*_excess)[v] + rem);
792 _flow->set(e, (*_capacity)[e]);
794 } else if (new_level > (*_level)[v]) {
795 new_level = (*_level)[v];
799 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
800 Value rem = (*_flow)[e];
801 if (!_tolerance.positive(rem)) continue;
802 Node v = _graph.source(e);
803 if ((*_level)[v] < level) {
804 if (!_level->active(v) && v != _source) {
807 if (!_tolerance.less(rem, excess)) {
808 _flow->set(e, (*_flow)[e] - excess);
809 _excess->set(v, (*_excess)[v] + excess);
814 _excess->set(v, (*_excess)[v] + rem);
817 } else if (new_level > (*_level)[v]) {
818 new_level = (*_level)[v];
824 _excess->set(n, excess);
827 if (new_level + 1 < _level->maxLevel()) {
828 _level->liftHighestActive(new_level + 1);
831 _level->liftHighestActiveToTop();
833 if (_level->emptyLevel(level)) {
835 _level->liftToTop(level);
838 _level->deactivate(n);
844 /// \brief Runs the preflow algorithm.
846 /// Runs the preflow algorithm.
847 /// \note pf.run() is just a shortcut of the following code.
850 /// pf.startFirstPhase();
851 /// pf.startSecondPhase();
859 /// \brief Runs the preflow algorithm to compute the minimum cut.
861 /// Runs the preflow algorithm to compute the minimum cut.
862 /// \note pf.runMinCut() is just a shortcut of the following code.
865 /// pf.startFirstPhase();
874 /// \name Query Functions
875 /// The result of the %Preflow algorithm can be obtained using these
877 /// Before the use of these functions,
878 /// either run() or start() must be called.
882 /// \brief Returns the value of the maximum flow.
884 /// Returns the value of the maximum flow by returning the excess
885 /// of the target node \c t. This value equals to the value of
886 /// the maximum flow already after the first phase.
887 Value flowValue() const {
888 return (*_excess)[_target];
891 /// \brief Returns true when the node is on the source side of minimum cut.
893 /// Returns true when the node is on the source side of minimum
894 /// cut. This method can be called both after running \ref
895 /// startFirstPhase() and \ref startSecondPhase().
896 bool minCut(const Node& node) const {
897 return ((*_level)[node] == _level->maxLevel()) == _phase;
900 /// \brief Returns a minimum value cut.
902 /// Sets the \c cutMap to the characteristic vector of a minimum value
903 /// cut. This method can be called both after running \ref
904 /// startFirstPhase() and \ref startSecondPhase(). The result after second
905 /// phase could be changed slightly if inexact computation is used.
906 /// \pre The \c cutMap should be a bool-valued node-map.
907 template <typename CutMap>
908 void minCutMap(CutMap& cutMap) const {
909 for (NodeIt n(_graph); n != INVALID; ++n) {
910 cutMap.set(n, minCut(n));
914 /// \brief Returns the flow on the edge.
916 /// Sets the \c flowMap to the flow on the edges. This method can
917 /// be called after the second phase of algorithm.
918 Value flow(const Edge& edge) const {
919 return (*_flow)[edge];