Redesigned CapacityScaling algorithm with almost the same interface.
The new version does not use the ResidualGraphAdaptor for performance reasons.
Scaling can be enabled and disabled with a parameter of the run() function.
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_DINITZ_SLEATOR_TARJAN_H
20 #define LEMON_DINITZ_SLEATOR_TARJAN_H
24 /// \brief Implementation the dynamic tree data structure of Sleator
27 #include <lemon/graph_utils.h>
28 #include <lemon/tolerance.h>
29 #include <lemon/dynamic_tree.h>
38 /// \brief Default traits class of DinitzSleatorTarjan class.
40 /// Default traits class of DinitzSleatorTarjan class.
41 /// \param _Graph Graph type.
42 /// \param _CapacityMap Type of capacity map.
43 template <typename _Graph, typename _CapacityMap>
44 struct DinitzSleatorTarjanDefaultTraits {
46 /// \brief The graph type the algorithm runs on.
49 /// \brief The type of the map that stores the edge capacities.
51 /// The type of the map that stores the edge capacities.
52 /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
53 typedef _CapacityMap CapacityMap;
55 /// \brief The type of the length of the edges.
56 typedef typename CapacityMap::Value Value;
58 /// \brief The map type that stores the flow values.
60 /// The map type that stores the flow values.
61 /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
62 typedef typename Graph::template EdgeMap<Value> FlowMap;
64 /// \brief Instantiates a FlowMap.
66 /// This function instantiates a \ref FlowMap.
67 /// \param graph The graph, to which we would like to define the flow map.
68 static FlowMap* createFlowMap(const Graph& graph) {
69 return new FlowMap(graph);
72 /// \brief The tolerance used by the algorithm
74 /// The tolerance used by the algorithm to handle inexact computation.
75 typedef Tolerance<Value> Tolerance;
81 /// \brief Dinitz-Sleator-Tarjan algorithms class.
83 /// This class provides an implementation of the \e
84 /// Dinitz-Sleator-Tarjan \e algorithm producing a flow of maximum
85 /// value in a directed graph. The DinitzSleatorTarjan algorithm is
86 /// the fastest known max flow algorithms wich using blocking
87 /// flow. It is an improvement of the Dinitz algorithm by using the
88 /// \ref DynamicTree "dynamic tree" data structure of Sleator and
91 /// This blocking flow algorithms builds a layered graph according
92 /// to \ref Bfs "breadth-first search" distance from the target node
93 /// in the reversed residual graph. The layered graph contains each
94 /// residual edge which steps one level down. After that the
95 /// algorithm constructs a blocking flow on the layered graph and
96 /// augments the overall flow with it. The number of the levels in
97 /// the layered graph is strictly increasing in each augmenting
98 /// phase therefore the number of the augmentings is at most
99 /// \f$n-1\f$. The length of each phase is at most
100 /// \f$O(m\log(n))\f$, that the overall time complexity is
101 /// \f$O(nm\log(n))\f$.
103 /// \param _Graph The directed graph type the algorithm runs on.
104 /// \param _CapacityMap The capacity map type.
105 /// \param _Traits Traits class to set various data types used by
106 /// the algorithm. The default traits class is \ref
107 /// DinitzSleatorTarjanDefaultTraits. See \ref
108 /// DinitzSleatorTarjanDefaultTraits for the documentation of a
109 /// Dinitz-Sleator-Tarjan traits class.
111 /// \author Tamas Hamori and Balazs Dezso
113 template <typename _Graph, typename _CapacityMap, typename _Traits>
115 template <typename _Graph,
116 typename _CapacityMap = typename _Graph::template EdgeMap<int>,
118 DinitzSleatorTarjanDefaultTraits<_Graph, _CapacityMap> >
120 class DinitzSleatorTarjan {
123 typedef _Traits Traits;
124 typedef typename Traits::Graph Graph;
125 typedef typename Traits::CapacityMap CapacityMap;
126 typedef typename Traits::Value Value;
128 typedef typename Traits::FlowMap FlowMap;
129 typedef typename Traits::Tolerance Tolerance;
134 GRAPH_TYPEDEFS(typename Graph);
137 typedef typename Graph::template NodeMap<int> LevelMap;
138 typedef typename Graph::template NodeMap<int> IntNodeMap;
139 typedef typename Graph::template NodeMap<Edge> EdgeNodeMap;
140 typedef DynamicTree<Value, IntNodeMap, Tolerance, false> DynTree;
145 const CapacityMap* _capacity;
147 Node _source, _target;
153 EdgeNodeMap* _dt_edges;
155 IntNodeMap* _dt_index;
158 std::vector<Node> _queue;
160 Tolerance _tolerance;
168 typedef DinitzSleatorTarjan Create;
170 ///\name Named template parameters
174 template <typename _FlowMap>
175 struct DefFlowMapTraits : public Traits {
176 typedef _FlowMap FlowMap;
177 static FlowMap *createFlowMap(const Graph&) {
178 throw UninitializedParameter();
182 /// \brief \ref named-templ-param "Named parameter" for setting
185 /// \ref named-templ-param "Named parameter" for setting FlowMap
187 template <typename _FlowMap>
189 : public DinitzSleatorTarjan<Graph, CapacityMap,
190 DefFlowMapTraits<_FlowMap> > {
191 typedef DinitzSleatorTarjan<Graph, CapacityMap,
192 DefFlowMapTraits<_FlowMap> > Create;
195 template <typename _Elevator>
196 struct DefElevatorTraits : public Traits {
197 typedef _Elevator Elevator;
198 static Elevator *createElevator(const Graph&, int) {
199 throw UninitializedParameter();
205 /// \brief \ref Exception for the case when the source equals the target.
207 /// \ref Exception for the case when the source equals the target.
209 class InvalidArgument : public lemon::LogicError {
211 virtual const char* what() const throw() {
212 return "lemon::DinitzSleatorTarjan::InvalidArgument";
218 DinitzSleatorTarjan() {}
222 /// \brief The constructor of the class.
224 /// The constructor of the class.
225 /// \param graph The directed graph the algorithm runs on.
226 /// \param capacity The capacity of the edges.
227 /// \param source The source node.
228 /// \param target The target node.
229 DinitzSleatorTarjan(const Graph& graph, const CapacityMap& capacity,
230 Node source, Node target)
231 : _graph(graph), _capacity(&capacity),
232 _source(source), _target(target),
233 _flow(0), _local_flow(false),
234 _level(0), _dt_edges(0),
235 _dt_index(0), _dt(0), _queue(),
236 _tolerance(), _flow_value(), _max_value()
238 if (_source == _target) {
239 throw InvalidArgument();
243 /// \brief Destrcutor.
246 ~DinitzSleatorTarjan() {
250 /// \brief Sets the capacity map.
252 /// Sets the capacity map.
253 /// \return \c (*this)
254 DinitzSleatorTarjan& capacityMap(const CapacityMap& map) {
259 /// \brief Sets the flow map.
261 /// Sets the flow map.
262 /// \return \c (*this)
263 DinitzSleatorTarjan& flowMap(FlowMap& map) {
272 /// \brief Returns the flow map.
274 /// \return The flow map.
275 const FlowMap& flowMap() {
279 /// \brief Sets the source node.
281 /// Sets the source node.
282 /// \return \c (*this)
283 DinitzSleatorTarjan& source(const Node& node) {
288 /// \brief Sets the target node.
290 /// Sets the target node.
291 /// \return \c (*this)
292 DinitzSleatorTarjan& target(const Node& node) {
297 /// \brief Sets the tolerance used by algorithm.
299 /// Sets the tolerance used by algorithm.
300 DinitzSleatorTarjan& tolerance(const Tolerance& tolerance) const {
301 _tolerance = tolerance;
303 _dt.tolerance(_tolerance);
308 /// \brief Returns the tolerance used by algorithm.
310 /// Returns the tolerance used by algorithm.
311 const Tolerance& tolerance() const {
317 void createStructures() {
319 _flow = Traits::createFlowMap(_graph);
323 _level = new LevelMap(_graph);
325 if (!_dt_index && !_dt) {
326 _dt_index = new IntNodeMap(_graph);
327 _dt = new DynTree(*_dt_index, _tolerance);
330 _dt_edges = new EdgeNodeMap(_graph);
332 _queue.resize(countNodes(_graph));
333 _max_value = _dt->maxValue();
336 void destroyStructures() {
354 bool createLayeredGraph() {
356 for (NodeIt n(_graph); n != INVALID; ++n) {
363 _level->set(_target, level);
365 int first = 0, last = 1, limit = 0;
367 while (first != last && (*_level)[_source] == -2) {
368 if (first == limit) {
373 Node n = _queue[first++];
375 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
376 Node v = _graph.target(e);
377 if ((*_level)[v] != -2) continue;
378 Value rem = (*_flow)[e];
379 if (!_tolerance.positive(rem)) continue;
380 _level->set(v, level);
384 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
385 Node v = _graph.source(e);
386 if ((*_level)[v] != -2) continue;
387 Value rem = (*_capacity)[e] - (*_flow)[e];
388 if (!_tolerance.positive(rem)) continue;
389 _level->set(v, level);
393 return (*_level)[_source] != -2;
397 for (NodeIt n(_graph); n != INVALID; ++n) {
398 _graph.firstOut((*_dt_edges)[n], n);
405 Node n = _dt->findCost(_source, rem);
407 _dt->addCost(_source, - rem);
410 _dt->addCost(n, _max_value);
412 Edge e = (*_dt_edges)[n];
413 if (_graph.source(e) == n) {
414 _flow->set(e, (*_capacity)[e]);
418 _graph.firstIn(e, n);
424 _dt_edges->set(n, e);
428 bool advance(Node n) {
429 Edge e = (*_dt_edges)[n];
430 if (e == INVALID) return false;
434 if (_graph.source(e) == n) {
435 u = _graph.target(e);
436 while ((*_level)[n] != (*_level)[u] + 1 ||
437 !_tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
439 if (e == INVALID) break;
440 u = _graph.target(e);
443 rem = (*_capacity)[e] - (*_flow)[e];
445 _graph.firstIn(e, n);
447 _dt_edges->set(n, INVALID);
450 u = _graph.source(e);
451 while ((*_level)[n] != (*_level)[u] + 1 ||
452 !_tolerance.positive((*_flow)[e])) {
455 _dt_edges->set(n, INVALID);
458 u = _graph.source(e);
463 u = _graph.source(e);
464 while ((*_level)[n] != (*_level)[u] + 1 ||
465 !_tolerance.positive((*_flow)[e])) {
468 _dt_edges->set(n, INVALID);
471 u = _graph.source(e);
476 _dt->addCost(n, - std::numeric_limits<Value>::max());
477 _dt->addCost(n, rem);
479 _dt_edges->set(n, e);
483 void retreat(Node n) {
486 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
487 Node u = _graph.target(e);
488 if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
490 _dt->findCost(u, rem);
493 _dt->addCost(u, - rem);
494 _dt->addCost(u, _max_value);
497 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
498 Node u = _graph.source(e);
499 if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
501 _dt->findCost(u, rem);
502 _flow->set(e, (*_capacity)[e] - rem);
504 _dt->addCost(u, - rem);
505 _dt->addCost(u, _max_value);
510 void extractTrees() {
511 for (NodeIt n(_graph); n != INVALID; ++n) {
513 Node w = _dt->findRoot(n);
518 Node u = _dt->findCost(n, rem);
521 _dt->addCost(u, - rem);
522 _dt->addCost(u, _max_value);
524 Edge e = (*_dt_edges)[u];
525 _dt_edges->set(u, INVALID);
527 if (u == _graph.source(e)) {
528 _flow->set(e, (*_capacity)[e] - rem);
533 w = _dt->findRoot(n);
541 /// \name Execution control The simplest way to execute the
542 /// algorithm is to use the \c run() member functions.
544 /// If you need more control on initial solution or
545 /// execution then you have to call one \ref init() function and then
546 /// the start() or multiple times the \c augment() member function.
550 /// \brief Initializes the algorithm
552 /// It sets the flow to empty flow.
557 for (NodeIt n(_graph); n != INVALID; ++n) {
559 _dt->addCost(n, _max_value);
562 for (EdgeIt it(_graph); it != INVALID; ++it) {
568 /// \brief Initializes the algorithm
570 /// Initializes the flow to the \c flowMap. The \c flowMap should
571 /// contain a feasible flow, ie. in each node excluding the source
572 /// and the target the incoming flow should be equal to the
574 template <typename FlowMap>
575 void flowInit(const FlowMap& flowMap) {
579 for (NodeIt n(_graph); n != INVALID; ++n) {
581 _dt->addCost(n, _max_value);
584 for (EdgeIt e(_graph); e != INVALID; ++e) {
585 _flow->set(e, flowMap[e]);
588 for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
589 _flow_value += (*_flow)[jt];
591 for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
592 _flow_value -= (*_flow)[jt];
596 /// \brief Initializes the algorithm
598 /// Initializes the flow to the \c flowMap. The \c flowMap should
599 /// contain a feasible flow, ie. in each node excluding the source
600 /// and the target the incoming flow should be equal to the
602 /// \return %False when the given flowMap does not contain
604 template <typename FlowMap>
605 bool checkedFlowInit(const FlowMap& flowMap) {
609 for (NodeIt n(_graph); n != INVALID; ++n) {
611 _dt->addCost(n, _max_value);
614 for (EdgeIt e(_graph); e != INVALID; ++e) {
615 _flow->set(e, flowMap[e]);
617 for (NodeIt it(_graph); it != INVALID; ++it) {
618 if (it == _source || it == _target) continue;
620 for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
621 outFlow += (*_flow)[jt];
624 for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
625 inFlow += (*_flow)[jt];
627 if (_tolerance.different(outFlow, inFlow)) {
631 for (EdgeIt it(_graph); it != INVALID; ++it) {
632 if (_tolerance.less((*_flow)[it], 0)) return false;
633 if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
636 for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
637 _flow_value += (*_flow)[jt];
639 for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
640 _flow_value -= (*_flow)[jt];
645 /// \brief Executes the algorithm
647 /// It runs augmenting phases by adding blocking flow until the
648 /// optimal solution is reached.
653 /// \brief Augments the flow with a blocking flow on a layered
656 /// This function builds a layered graph and then find a blocking
657 /// flow on this graph. The number of the levels in the layered
658 /// graph is strictly increasing in each augmenting phase
659 /// therefore the number of the augmentings is at most \f$ n-1
660 /// \f$. The length of each phase is at most \f$ O(m \log(n))
661 /// \f$, that the overall time complexity is \f$ O(nm \log(n)) \f$.
662 /// \return %False when there is not residual path between the
663 /// source and the target so the current flow is a feasible and
664 /// optimal solution.
668 if (createLayeredGraph()) {
673 n = _dt->findRoot(_source);
678 } else if (!advance(n)) {
685 n = _dt->findRoot(_source);
695 /// \brief runs the algorithm.
697 /// It is just a shorthand for:
710 /// \name Query Functions
711 /// The result of the Dinitz-Sleator-Tarjan algorithm can be
712 /// obtained using these functions.
714 /// Before the use of these functions,
715 /// either run() or start() must be called.
719 /// \brief Returns the value of the maximum flow.
721 /// Returns the value of the maximum flow by returning the excess
722 /// of the target node \c t. This value equals to the value of
723 /// the maximum flow already after the first phase.
724 Value flowValue() const {
729 /// \brief Returns the flow on the edge.
731 /// Sets the \c flowMap to the flow on the edges. This method can
732 /// be called after the second phase of algorithm.
733 Value flow(const Edge& edge) const {
734 return (*_flow)[edge];
737 /// \brief Returns true when the node is on the source side of minimum cut.
740 /// Returns true when the node is on the source side of minimum
741 /// cut. This method can be called both after running \ref
742 /// startFirstPhase() and \ref startSecondPhase().
743 bool minCut(const Node& node) const {
744 return (*_level)[node] == -2;
747 /// \brief Returns a minimum value cut.
749 /// Sets \c cut to the characteristic vector of a minimum value cut
750 /// It simply calls the minMinCut member.
751 /// \retval cut Write node bool map.
752 template <typename CutMap>
753 void minCutMap(CutMap& cutMap) const {
754 for (NodeIt n(_graph); n != INVALID; ++n) {
755 cutMap.set(n, (*_level)[n] == -2);
757 cutMap.set(_source, true);