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_MIN_COST_ARBORESCENCE_H
20 #define LEMON_MIN_COST_ARBORESCENCE_H
24 ///\brief Minimum Cost Arborescence algorithm.
28 #include <lemon/list_graph.h>
29 #include <lemon/bin_heap.h>
34 /// \brief Default traits class of MinCostArborescence class.
36 /// Default traits class of MinCostArborescence class.
37 /// \param _Graph Graph type.
38 /// \param _CostMap Type of cost map.
39 template <class _Graph, class _CostMap>
40 struct MinCostArborescenceDefaultTraits{
42 /// \brief The graph type the algorithm runs on.
45 /// \brief The type of the map that stores the edge costs.
47 /// The type of the map that stores the edge costs.
48 /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
49 typedef _CostMap CostMap;
51 /// \brief The value type of the costs.
53 /// The value type of the costs.
54 typedef typename CostMap::Value Value;
56 /// \brief The type of the map that stores which edges are
57 /// in the arborescence.
59 /// The type of the map that stores which edges are in the arborescence.
60 /// It must meet the \ref concepts::WriteMap "WriteMap" concept.
61 /// Initially it will be set to false on each edge. After it
62 /// will set all arborescence edges once.
63 typedef typename Graph::template EdgeMap<bool> ArborescenceMap;
65 /// \brief Instantiates a ArborescenceMap.
67 /// This function instantiates a \ref ArborescenceMap.
68 /// \param _graph is the graph, to which we would like to define the
70 static ArborescenceMap *createArborescenceMap(const Graph &_graph){
71 return new ArborescenceMap(_graph);
74 /// \brief The type of the PredMap
76 /// The type of the PredMap. It is a node map with an edge value type.
77 typedef typename Graph::template NodeMap<typename Graph::Edge> PredMap;
79 /// \brief Instantiates a PredMap.
81 /// This function instantiates a \ref PredMap.
82 /// \param _graph is the graph, to which we would like to define the
84 static PredMap *createPredMap(const Graph &_graph){
85 return new PredMap(_graph);
92 /// \brief %MinCostArborescence algorithm class.
94 /// This class provides an efficient implementation of
95 /// %MinCostArborescence algorithm. The arborescence is a tree
96 /// which is directed from a given source node of the graph. One or
97 /// more sources should be given for the algorithm and it will calculate
98 /// the minimum cost subgraph which are union of arborescences with the
99 /// given sources and spans all the nodes which are reachable from the
100 /// sources. The time complexity of the algorithm is \f$ O(n^2+e) \f$.
102 /// The algorithm provides also an optimal dual solution to arborescence
103 /// that way the optimality of the solution can be proofed easily.
105 /// \param _Graph The graph type the algorithm runs on. The default value
106 /// is \ref ListGraph. The value of _Graph is not used directly by
107 /// MinCostArborescence, it is only passed to
108 /// \ref MinCostArborescenceDefaultTraits.
109 /// \param _CostMap This read-only EdgeMap determines the costs of the
110 /// edges. It is read once for each edge, so the map may involve in
111 /// relatively time consuming process to compute the edge cost if
112 /// it is necessary. The default map type is \ref
113 /// concepts::Graph::EdgeMap "Graph::EdgeMap<int>". The value
114 /// of _CostMap is not used directly by MinCostArborescence,
115 /// it is only passed to \ref MinCostArborescenceDefaultTraits.
116 /// \param _Traits Traits class to set various data types used
117 /// by the algorithm. The default traits class is
118 /// \ref MinCostArborescenceDefaultTraits
119 /// "MinCostArborescenceDefaultTraits<_Graph,_CostMap>". See \ref
120 /// MinCostArborescenceDefaultTraits for the documentation of a
121 /// MinCostArborescence traits class.
123 /// \author Balazs Dezso
125 template <typename _Graph = ListGraph,
126 typename _CostMap = typename _Graph::template EdgeMap<int>,
128 MinCostArborescenceDefaultTraits<_Graph, _CostMap> >
130 template <typename _Graph, typename _CostMap, typedef _Traits>
132 class MinCostArborescence {
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::UninitializedParameter {
141 virtual const char* what() const throw() {
142 return "lemon::MinCostArborescence::UninitializedParameter";
147 typedef _Traits Traits;
148 /// The type of the underlying graph.
149 typedef typename Traits::Graph Graph;
150 /// The type of the map that stores the edge costs.
151 typedef typename Traits::CostMap CostMap;
152 ///The type of the costs of the edges.
153 typedef typename Traits::Value Value;
154 ///The type of the predecessor map.
155 typedef typename Traits::PredMap PredMap;
156 ///The type of the map that stores which edges are in the arborescence.
157 typedef typename Traits::ArborescenceMap ArborescenceMap;
161 typedef typename Graph::Node Node;
162 typedef typename Graph::Edge Edge;
163 typedef typename Graph::NodeIt NodeIt;
164 typedef typename Graph::EdgeIt EdgeIt;
165 typedef typename Graph::InEdgeIt InEdgeIt;
166 typedef typename Graph::OutEdgeIt OutEdgeIt;
174 CostEdge(Edge _edge, Value _value) : edge(_edge), value(_value) {}
184 ArborescenceMap *_arborescence;
185 bool local_arborescence;
187 typedef typename Graph::template EdgeMap<int> EdgeOrder;
188 EdgeOrder *_edge_order;
190 typedef typename Graph::template NodeMap<int> NodeOrder;
191 NodeOrder *_node_order;
193 typedef typename Graph::template NodeMap<CostEdge> CostEdgeMap;
194 CostEdgeMap *_cost_edges;
198 std::vector<CostEdge> edges;
203 std::vector<StackLevel> level_stack;
204 std::vector<Node> queue;
206 typedef std::vector<typename Graph::Node> DualNodeList;
208 DualNodeList _dual_node_list;
210 struct DualVariable {
214 DualVariable(int _begin, int _end, Value _value)
215 : begin(_begin), end(_end), value(_value) {}
219 typedef std::vector<DualVariable> DualVariables;
221 DualVariables _dual_variables;
223 typedef typename Graph::template NodeMap<int> HeapCrossRef;
225 HeapCrossRef *_heap_cross_ref;
227 typedef BinHeap<int, HeapCrossRef> Heap;
233 /// \name Named template parameters
238 struct DefArborescenceMapTraits : public Traits {
239 typedef T ArborescenceMap;
240 static ArborescenceMap *createArborescenceMap(const Graph &)
242 throw UninitializedParameter();
246 /// \brief \ref named-templ-param "Named parameter" for
247 /// setting ArborescenceMap type
249 /// \ref named-templ-param "Named parameter" for setting
250 /// ArborescenceMap type
252 struct DefArborescenceMap
253 : public MinCostArborescence<Graph, CostMap,
254 DefArborescenceMapTraits<T> > {
255 typedef MinCostArborescence<Graph, CostMap,
256 DefArborescenceMapTraits<T> > Create;
260 struct DefPredMapTraits : public Traits {
262 static PredMap *createPredMap(const Graph &)
264 throw UninitializedParameter();
268 /// \brief \ref named-templ-param "Named parameter" for
269 /// setting PredMap type
271 /// \ref named-templ-param "Named parameter" for setting
275 : public MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > {
276 typedef MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > Create;
281 /// \brief Constructor.
283 /// \param _graph The graph the algorithm will run on.
284 /// \param _cost The cost map used by the algorithm.
285 MinCostArborescence(const Graph& _graph, const CostMap& _cost)
286 : graph(&_graph), cost(&_cost), _pred(0), local_pred(false),
287 _arborescence(0), local_arborescence(false),
288 _edge_order(0), _node_order(0), _cost_edges(0),
289 _heap_cross_ref(0), _heap(0) {}
291 /// \brief Destructor.
292 ~MinCostArborescence() {
296 /// \brief Sets the arborescence map.
298 /// Sets the arborescence map.
299 /// \return \c (*this)
300 MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
301 if (local_arborescence) {
302 delete _arborescence;
304 local_arborescence = false;
309 /// \brief Sets the arborescence map.
311 /// Sets the arborescence map.
312 /// \return \c (*this)
313 MinCostArborescence& predMap(PredMap& m) {
322 /// \name Query Functions
323 /// The result of the %MinCostArborescence algorithm can be obtained
324 /// using these functions.\n
325 /// Before the use of these functions,
326 /// either run() or start() must be called.
330 /// \brief Returns a reference to the arborescence map.
332 /// Returns a reference to the arborescence map.
333 const ArborescenceMap& arborescenceMap() const {
334 return *_arborescence;
337 /// \brief Returns true if the edge is in the arborescence.
339 /// Returns true if the edge is in the arborescence.
340 /// \param edge The edge of the graph.
341 /// \pre \ref run() must be called before using this function.
342 bool arborescence(Edge edge) const {
343 return (*_pred)[graph->target(edge)] == edge;
346 /// \brief Returns a reference to the pred map.
348 /// Returns a reference to the pred map.
349 const PredMap& predMap() const {
353 /// \brief Returns the predecessor edge of the given node.
355 /// Returns the predecessor edge of the given node.
356 bool pred(Node node) const {
357 return (*_pred)[node];
360 /// \brief Returns the cost of the arborescence.
362 /// Returns the cost of the arborescence.
363 Value arborescenceValue() const {
365 for (EdgeIt it(*graph); it != INVALID; ++it) {
366 if (arborescence(it)) {
373 /// \brief Indicates that a node is reachable from the sources.
375 /// Indicates that a node is reachable from the sources.
376 bool reached(Node node) const {
377 return (*_node_order)[node] != -3;
380 /// \brief Indicates that a node is processed.
382 /// Indicates that a node is processed. The arborescence path exists
383 /// from the source to the given node.
384 bool processed(Node node) const {
385 return (*_node_order)[node] == -1;
388 /// \brief Returns the number of the dual variables in basis.
390 /// Returns the number of the dual variables in basis.
391 int dualSize() const {
392 return _dual_variables.size();
395 /// \brief Returns the value of the dual solution.
397 /// Returns the value of the dual solution. It should be
398 /// equal to the arborescence value.
399 Value dualValue() const {
401 for (int i = 0; i < int(_dual_variables.size()); ++i) {
402 sum += _dual_variables[i].value;
407 /// \brief Returns the number of the nodes in the dual variable.
409 /// Returns the number of the nodes in the dual variable.
410 int dualSize(int k) const {
411 return _dual_variables[k].end - _dual_variables[k].begin;
414 /// \brief Returns the value of the dual variable.
416 /// Returns the the value of the dual variable.
417 const Value& dualValue(int k) const {
418 return _dual_variables[k].value;
421 /// \brief Lemon iterator for get a dual variable.
423 /// Lemon iterator for get a dual variable. This class provides
424 /// a common style lemon iterator which gives back a subset of
429 /// \brief Constructor.
431 /// Constructor for get the nodeset of the variable.
432 DualIt(const MinCostArborescence& algorithm, int variable)
433 : _algorithm(&algorithm), _variable(variable)
435 _index = _algorithm->_dual_variables[_variable].begin;
438 /// \brief Invalid constructor.
440 /// Invalid constructor.
441 DualIt(Invalid) : _algorithm(0) {}
443 /// \brief Conversion to node.
445 /// Conversion to node.
446 operator Node() const {
447 return _algorithm ? _algorithm->_dual_node_list[_index] : INVALID;
450 /// \brief Increment operator.
452 /// Increment operator.
453 DualIt& operator++() {
455 if (_algorithm->_dual_variables[_variable].end == _index) {
461 bool operator==(const DualIt& it) const {
462 return static_cast<Node>(*this) == static_cast<Node>(it);
464 bool operator!=(const DualIt& it) const {
465 return static_cast<Node>(*this) != static_cast<Node>(it);
467 bool operator<(const DualIt& it) const {
468 return static_cast<Node>(*this) < static_cast<Node>(it);
472 const MinCostArborescence* _algorithm;
479 /// \name Execution control
480 /// The simplest way to execute the algorithm is to use
481 /// one of the member functions called \c run(...). \n
482 /// If you need more control on the execution,
483 /// first you must call \ref init(), then you can add several
484 /// source nodes with \ref addSource().
485 /// Finally \ref start() will perform the arborescence
490 /// \brief Initializes the internal data structures.
492 /// Initializes the internal data structures.
497 for (NodeIt it(*graph); it != INVALID; ++it) {
498 (*_cost_edges)[it].edge = INVALID;
499 _node_order->set(it, -3);
500 _heap_cross_ref->set(it, Heap::PRE_HEAP);
501 _pred->set(it, INVALID);
503 for (EdgeIt it(*graph); it != INVALID; ++it) {
504 _arborescence->set(it, false);
505 _edge_order->set(it, -1);
507 _dual_node_list.clear();
508 _dual_variables.clear();
511 /// \brief Adds a new source node.
513 /// Adds a new source node to the algorithm.
514 void addSource(Node source) {
515 std::vector<Node> nodes;
516 nodes.push_back(source);
517 while (!nodes.empty()) {
518 Node node = nodes.back();
520 for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
521 Node target = graph->target(it);
522 if ((*_node_order)[target] == -3) {
523 (*_node_order)[target] = -2;
524 nodes.push_back(target);
525 queue.push_back(target);
529 (*_node_order)[source] = -1;
532 /// \brief Processes the next node in the priority queue.
534 /// Processes the next node in the priority queue.
536 /// \return The processed node.
538 /// \warning The queue must not be empty!
539 Node processNextNode() {
540 Node node = queue.back();
542 if ((*_node_order)[node] == -2) {
543 Edge edge = prepare(node);
544 Node source = graph->source(edge);
545 while ((*_node_order)[source] != -1) {
546 if ((*_node_order)[source] >= 0) {
547 edge = contract(source);
549 edge = prepare(source);
551 source = graph->source(edge);
559 /// \brief Returns the number of the nodes to be processed.
561 /// Returns the number of the nodes to be processed.
562 int queueSize() const {
566 /// \brief Returns \c false if there are nodes to be processed.
568 /// Returns \c false if there are nodes to be processed.
569 bool emptyQueue() const {
570 return queue.empty();
573 /// \brief Executes the algorithm.
575 /// Executes the algorithm.
577 /// \pre init() must be called and at least one node should be added
578 /// with addSource() before using this function.
580 ///\note mca.start() is just a shortcut of the following code.
582 ///while (!mca.emptyQueue()) {
583 /// mca.processNextNode();
587 while (!emptyQueue()) {
592 /// \brief Runs %MinCostArborescence algorithm from node \c s.
594 /// This method runs the %MinCostArborescence algorithm from
595 /// a root node \c s.
597 ///\note mca.run(s) is just a shortcut of the following code.
603 void run(Node node) {
613 void initStructures() {
616 _pred = Traits::createPredMap(*graph);
618 if (!_arborescence) {
619 local_arborescence = true;
620 _arborescence = Traits::createArborescenceMap(*graph);
623 _edge_order = new EdgeOrder(*graph);
626 _node_order = new NodeOrder(*graph);
629 _cost_edges = new CostEdgeMap(*graph);
631 if (!_heap_cross_ref) {
632 _heap_cross_ref = new HeapCrossRef(*graph, -1);
635 _heap = new Heap(*_heap_cross_ref);
639 void destroyStructures() {
640 if (local_arborescence) {
641 delete _arborescence;
658 if (!_heap_cross_ref) {
659 delete _heap_cross_ref;
663 Edge prepare(Node node) {
664 std::vector<Node> nodes;
665 (*_node_order)[node] = _dual_node_list.size();
667 level.node_level = _dual_node_list.size();
668 _dual_node_list.push_back(node);
669 for (InEdgeIt it(*graph, node); it != INVALID; ++it) {
671 Node source = graph->source(edge);
672 Value value = (*cost)[it];
673 if (source == node || (*_node_order)[source] == -3) continue;
674 if ((*_cost_edges)[source].edge == INVALID) {
675 (*_cost_edges)[source].edge = edge;
676 (*_cost_edges)[source].value = value;
677 nodes.push_back(source);
679 if ((*_cost_edges)[source].value > value) {
680 (*_cost_edges)[source].edge = edge;
681 (*_cost_edges)[source].value = value;
685 CostEdge minimum = (*_cost_edges)[nodes[0]];
686 for (int i = 1; i < int(nodes.size()); ++i) {
687 if ((*_cost_edges)[nodes[i]].value < minimum.value) {
688 minimum = (*_cost_edges)[nodes[i]];
691 _edge_order->set(minimum.edge, _dual_variables.size());
692 DualVariable var(_dual_node_list.size() - 1,
693 _dual_node_list.size(), minimum.value);
694 _dual_variables.push_back(var);
695 for (int i = 0; i < int(nodes.size()); ++i) {
696 (*_cost_edges)[nodes[i]].value -= minimum.value;
697 level.edges.push_back((*_cost_edges)[nodes[i]]);
698 (*_cost_edges)[nodes[i]].edge = INVALID;
700 level_stack.push_back(level);
704 Edge contract(Node node) {
705 int node_bottom = bottom(node);
706 std::vector<Node> nodes;
707 while (!level_stack.empty() &&
708 level_stack.back().node_level >= node_bottom) {
709 for (int i = 0; i < int(level_stack.back().edges.size()); ++i) {
710 Edge edge = level_stack.back().edges[i].edge;
711 Node source = graph->source(edge);
712 Value value = level_stack.back().edges[i].value;
713 if ((*_node_order)[source] >= node_bottom) continue;
714 if ((*_cost_edges)[source].edge == INVALID) {
715 (*_cost_edges)[source].edge = edge;
716 (*_cost_edges)[source].value = value;
717 nodes.push_back(source);
719 if ((*_cost_edges)[source].value > value) {
720 (*_cost_edges)[source].edge = edge;
721 (*_cost_edges)[source].value = value;
725 level_stack.pop_back();
727 CostEdge minimum = (*_cost_edges)[nodes[0]];
728 for (int i = 1; i < int(nodes.size()); ++i) {
729 if ((*_cost_edges)[nodes[i]].value < minimum.value) {
730 minimum = (*_cost_edges)[nodes[i]];
733 _edge_order->set(minimum.edge, _dual_variables.size());
734 DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
735 _dual_variables.push_back(var);
737 level.node_level = node_bottom;
738 for (int i = 0; i < int(nodes.size()); ++i) {
739 (*_cost_edges)[nodes[i]].value -= minimum.value;
740 level.edges.push_back((*_cost_edges)[nodes[i]]);
741 (*_cost_edges)[nodes[i]].edge = INVALID;
743 level_stack.push_back(level);
747 int bottom(Node node) {
748 int k = level_stack.size() - 1;
749 while (level_stack[k].node_level > (*_node_order)[node]) {
752 return level_stack[k].node_level;
755 void finalize(Edge edge) {
756 Node node = graph->target(edge);
757 _heap->push(node, (*_edge_order)[edge]);
758 _pred->set(node, edge);
759 while (!_heap->empty()) {
760 Node source = _heap->top();
762 _node_order->set(source, -1);
763 for (OutEdgeIt it(*graph, source); it != INVALID; ++it) {
764 if ((*_edge_order)[it] < 0) continue;
765 Node target = graph->target(it);
766 switch(_heap->state(target)) {
768 _heap->push(target, (*_edge_order)[it]);
769 _pred->set(target, it);
772 if ((*_edge_order)[it] < (*_heap)[target]) {
773 _heap->decrease(target, (*_edge_order)[it]);
774 _pred->set(target, it);
777 case Heap::POST_HEAP:
781 _arborescence->set((*_pred)[source], true);
787 /// \ingroup spantree
789 /// \brief Function type interface for MinCostArborescence algorithm.
791 /// Function type interface for MinCostArborescence algorithm.
792 /// \param graph The Graph that the algorithm runs on.
793 /// \param cost The CostMap of the edges.
794 /// \param source The source of the arborescence.
795 /// \retval arborescence The bool EdgeMap which stores the arborescence.
796 /// \return The cost of the arborescence.
798 /// \sa MinCostArborescence
799 template <typename Graph, typename CostMap, typename ArborescenceMap>
800 typename CostMap::Value minCostArborescence(const Graph& graph,
802 typename Graph::Node source,
803 ArborescenceMap& arborescence) {
804 typename MinCostArborescence<Graph, CostMap>
805 ::template DefArborescenceMap<ArborescenceMap>
806 ::Create mca(graph, cost);
807 mca.arborescenceMap(arborescence);
809 return mca.arborescenceValue();