1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
3 * This file is a part of LEMON, a generic C++ optimization library.
5 * Copyright (C) 2003-2011
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>
30 #include <lemon/assert.h>
35 /// \brief Default traits class for MinCostArborescence class.
37 /// Default traits class for MinCostArborescence class.
38 /// \param GR Digraph type.
39 /// \param CM Type of the cost map.
40 template <class GR, class CM>
41 struct MinCostArborescenceDefaultTraits{
43 /// \brief The digraph type the algorithm runs on.
46 /// \brief The type of the map that stores the arc costs.
48 /// The type of the map that stores the arc costs.
49 /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
52 /// \brief The value type of the costs.
54 /// The value type of the costs.
55 typedef typename CostMap::Value Value;
57 /// \brief The type of the map that stores which arcs are in the
60 /// The type of the map that stores which arcs are in the
61 /// arborescence. It must conform to the \ref concepts::WriteMap
62 /// "WriteMap" concept, and its value type must be \c bool
63 /// (or convertible). Initially it will be set to \c false on each
64 /// arc, then it will be set on each arborescence arc once.
65 typedef typename Digraph::template ArcMap<bool> ArborescenceMap;
67 /// \brief Instantiates a \c ArborescenceMap.
69 /// This function instantiates a \c ArborescenceMap.
70 /// \param digraph The digraph to which we would like to calculate
71 /// the \c ArborescenceMap.
72 static ArborescenceMap *createArborescenceMap(const Digraph &digraph){
73 return new ArborescenceMap(digraph);
76 /// \brief The type of the \c PredMap
78 /// The type of the \c PredMap. It must confrom to the
79 /// \ref concepts::WriteMap "WriteMap" concept, and its value type
80 /// must be the \c Arc type of the digraph.
81 typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
83 /// \brief Instantiates a \c PredMap.
85 /// This function instantiates a \c PredMap.
86 /// \param digraph The digraph to which we would like to define the
88 static PredMap *createPredMap(const Digraph &digraph){
89 return new PredMap(digraph);
96 /// \brief Minimum Cost Arborescence algorithm class.
98 /// This class provides an efficient implementation of the
99 /// Minimum Cost Arborescence algorithm. The arborescence is a tree
100 /// which is directed from a given source node of the digraph. One or
101 /// more sources should be given to the algorithm and it will calculate
102 /// the minimum cost subgraph that is the union of arborescences with the
103 /// given sources and spans all the nodes which are reachable from the
104 /// sources. The time complexity of the algorithm is O(n<sup>2</sup>+e).
106 /// The algorithm also provides an optimal dual solution, therefore
107 /// the optimality of the solution can be checked.
109 /// \param GR The digraph type the algorithm runs on.
110 /// \param CM A read-only arc map storing the costs of the
111 /// arcs. It is read once for each arc, so the map may involve in
112 /// relatively time consuming process to compute the arc costs if
113 /// it is necessary. The default map type is \ref
114 /// concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
115 /// \param TR Traits class to set various data types used
116 /// by the algorithm. The default traits class is
117 /// \ref MinCostArborescenceDefaultTraits
118 /// "MinCostArborescenceDefaultTraits<GR, CM>".
120 template <typename GR,
121 typename CM = typename GR::template ArcMap<int>,
123 MinCostArborescenceDefaultTraits<GR, CM> >
125 template <typename GR, typename CM, typedef TR>
127 class MinCostArborescence {
130 /// \brief The \ref MinCostArborescenceDefaultTraits "traits class"
131 /// of the algorithm.
133 /// The type of the underlying digraph.
134 typedef typename Traits::Digraph Digraph;
135 /// The type of the map that stores the arc costs.
136 typedef typename Traits::CostMap CostMap;
137 ///The type of the costs of the arcs.
138 typedef typename Traits::Value Value;
139 ///The type of the predecessor map.
140 typedef typename Traits::PredMap PredMap;
141 ///The type of the map that stores which arcs are in the arborescence.
142 typedef typename Traits::ArborescenceMap ArborescenceMap;
144 typedef MinCostArborescence Create;
148 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
156 CostArc(Arc _arc, Value _value) : arc(_arc), value(_value) {}
160 const Digraph *_digraph;
161 const CostMap *_cost;
166 ArborescenceMap *_arborescence;
167 bool local_arborescence;
169 typedef typename Digraph::template ArcMap<int> ArcOrder;
170 ArcOrder *_arc_order;
172 typedef typename Digraph::template NodeMap<int> NodeOrder;
173 NodeOrder *_node_order;
175 typedef typename Digraph::template NodeMap<CostArc> CostArcMap;
176 CostArcMap *_cost_arcs;
180 std::vector<CostArc> arcs;
185 std::vector<StackLevel> level_stack;
186 std::vector<Node> queue;
188 typedef std::vector<typename Digraph::Node> DualNodeList;
190 DualNodeList _dual_node_list;
192 struct DualVariable {
196 DualVariable(int _begin, int _end, Value _value)
197 : begin(_begin), end(_end), value(_value) {}
201 typedef std::vector<DualVariable> DualVariables;
203 DualVariables _dual_variables;
205 typedef typename Digraph::template NodeMap<int> HeapCrossRef;
207 HeapCrossRef *_heap_cross_ref;
209 typedef BinHeap<int, HeapCrossRef> Heap;
215 MinCostArborescence() {}
219 void createStructures() {
222 _pred = Traits::createPredMap(*_digraph);
224 if (!_arborescence) {
225 local_arborescence = true;
226 _arborescence = Traits::createArborescenceMap(*_digraph);
229 _arc_order = new ArcOrder(*_digraph);
232 _node_order = new NodeOrder(*_digraph);
235 _cost_arcs = new CostArcMap(*_digraph);
237 if (!_heap_cross_ref) {
238 _heap_cross_ref = new HeapCrossRef(*_digraph, -1);
241 _heap = new Heap(*_heap_cross_ref);
245 void destroyStructures() {
246 if (local_arborescence) {
247 delete _arborescence;
264 if (_heap_cross_ref) {
265 delete _heap_cross_ref;
269 Arc prepare(Node node) {
270 std::vector<Node> nodes;
271 (*_node_order)[node] = _dual_node_list.size();
273 level.node_level = _dual_node_list.size();
274 _dual_node_list.push_back(node);
275 for (InArcIt it(*_digraph, node); it != INVALID; ++it) {
277 Node source = _digraph->source(arc);
278 Value value = (*_cost)[it];
279 if (source == node || (*_node_order)[source] == -3) continue;
280 if ((*_cost_arcs)[source].arc == INVALID) {
281 (*_cost_arcs)[source].arc = arc;
282 (*_cost_arcs)[source].value = value;
283 nodes.push_back(source);
285 if ((*_cost_arcs)[source].value > value) {
286 (*_cost_arcs)[source].arc = arc;
287 (*_cost_arcs)[source].value = value;
291 CostArc minimum = (*_cost_arcs)[nodes[0]];
292 for (int i = 1; i < int(nodes.size()); ++i) {
293 if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
294 minimum = (*_cost_arcs)[nodes[i]];
297 (*_arc_order)[minimum.arc] = _dual_variables.size();
298 DualVariable var(_dual_node_list.size() - 1,
299 _dual_node_list.size(), minimum.value);
300 _dual_variables.push_back(var);
301 for (int i = 0; i < int(nodes.size()); ++i) {
302 (*_cost_arcs)[nodes[i]].value -= minimum.value;
303 level.arcs.push_back((*_cost_arcs)[nodes[i]]);
304 (*_cost_arcs)[nodes[i]].arc = INVALID;
306 level_stack.push_back(level);
310 Arc contract(Node node) {
311 int node_bottom = bottom(node);
312 std::vector<Node> nodes;
313 while (!level_stack.empty() &&
314 level_stack.back().node_level >= node_bottom) {
315 for (int i = 0; i < int(level_stack.back().arcs.size()); ++i) {
316 Arc arc = level_stack.back().arcs[i].arc;
317 Node source = _digraph->source(arc);
318 Value value = level_stack.back().arcs[i].value;
319 if ((*_node_order)[source] >= node_bottom) continue;
320 if ((*_cost_arcs)[source].arc == INVALID) {
321 (*_cost_arcs)[source].arc = arc;
322 (*_cost_arcs)[source].value = value;
323 nodes.push_back(source);
325 if ((*_cost_arcs)[source].value > value) {
326 (*_cost_arcs)[source].arc = arc;
327 (*_cost_arcs)[source].value = value;
331 level_stack.pop_back();
333 CostArc minimum = (*_cost_arcs)[nodes[0]];
334 for (int i = 1; i < int(nodes.size()); ++i) {
335 if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
336 minimum = (*_cost_arcs)[nodes[i]];
339 (*_arc_order)[minimum.arc] = _dual_variables.size();
340 DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
341 _dual_variables.push_back(var);
343 level.node_level = node_bottom;
344 for (int i = 0; i < int(nodes.size()); ++i) {
345 (*_cost_arcs)[nodes[i]].value -= minimum.value;
346 level.arcs.push_back((*_cost_arcs)[nodes[i]]);
347 (*_cost_arcs)[nodes[i]].arc = INVALID;
349 level_stack.push_back(level);
353 int bottom(Node node) {
354 int k = level_stack.size() - 1;
355 while (level_stack[k].node_level > (*_node_order)[node]) {
358 return level_stack[k].node_level;
361 void finalize(Arc arc) {
362 Node node = _digraph->target(arc);
363 _heap->push(node, (*_arc_order)[arc]);
364 _pred->set(node, arc);
365 while (!_heap->empty()) {
366 Node source = _heap->top();
368 (*_node_order)[source] = -1;
369 for (OutArcIt it(*_digraph, source); it != INVALID; ++it) {
370 if ((*_arc_order)[it] < 0) continue;
371 Node target = _digraph->target(it);
372 switch(_heap->state(target)) {
374 _heap->push(target, (*_arc_order)[it]);
375 _pred->set(target, it);
378 if ((*_arc_order)[it] < (*_heap)[target]) {
379 _heap->decrease(target, (*_arc_order)[it]);
380 _pred->set(target, it);
383 case Heap::POST_HEAP:
387 _arborescence->set((*_pred)[source], true);
394 /// \name Named Template Parameters
399 struct SetArborescenceMapTraits : public Traits {
400 typedef T ArborescenceMap;
401 static ArborescenceMap *createArborescenceMap(const Digraph &)
403 LEMON_ASSERT(false, "ArborescenceMap is not initialized");
404 return 0; // ignore warnings
408 /// \brief \ref named-templ-param "Named parameter" for
409 /// setting \c ArborescenceMap type
411 /// \ref named-templ-param "Named parameter" for setting
412 /// \c ArborescenceMap type.
413 /// It must conform to the \ref concepts::WriteMap "WriteMap" concept,
414 /// and its value type must be \c bool (or convertible).
415 /// Initially it will be set to \c false on each arc,
416 /// then it will be set on each arborescence arc once.
418 struct SetArborescenceMap
419 : public MinCostArborescence<Digraph, CostMap,
420 SetArborescenceMapTraits<T> > {
424 struct SetPredMapTraits : public Traits {
426 static PredMap *createPredMap(const Digraph &)
428 LEMON_ASSERT(false, "PredMap is not initialized");
429 return 0; // ignore warnings
433 /// \brief \ref named-templ-param "Named parameter" for
434 /// setting \c PredMap type
436 /// \ref named-templ-param "Named parameter" for setting
438 /// It must meet the \ref concepts::WriteMap "WriteMap" concept,
439 /// and its value type must be the \c Arc type of the digraph.
442 : public MinCostArborescence<Digraph, CostMap, SetPredMapTraits<T> > {
447 /// \brief Constructor.
449 /// \param digraph The digraph the algorithm will run on.
450 /// \param cost The cost map used by the algorithm.
451 MinCostArborescence(const Digraph& digraph, const CostMap& cost)
452 : _digraph(&digraph), _cost(&cost), _pred(0), local_pred(false),
453 _arborescence(0), local_arborescence(false),
454 _arc_order(0), _node_order(0), _cost_arcs(0),
455 _heap_cross_ref(0), _heap(0) {}
457 /// \brief Destructor.
458 ~MinCostArborescence() {
462 /// \brief Sets the arborescence map.
464 /// Sets the arborescence map.
465 /// \return <tt>(*this)</tt>
466 MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
467 if (local_arborescence) {
468 delete _arborescence;
470 local_arborescence = false;
475 /// \brief Sets the predecessor map.
477 /// Sets the predecessor map.
478 /// \return <tt>(*this)</tt>
479 MinCostArborescence& predMap(PredMap& m) {
488 /// \name Execution Control
489 /// The simplest way to execute the algorithm is to use
490 /// one of the member functions called \c run(...). \n
491 /// If you need more control on the execution,
492 /// first you must call \ref init(), then you can add several
493 /// source nodes with \ref addSource().
494 /// Finally \ref start() will perform the arborescence
499 /// \brief Initializes the internal data structures.
501 /// Initializes the internal data structures.
506 for (NodeIt it(*_digraph); it != INVALID; ++it) {
507 (*_cost_arcs)[it].arc = INVALID;
508 (*_node_order)[it] = -3;
509 (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
510 _pred->set(it, INVALID);
512 for (ArcIt it(*_digraph); it != INVALID; ++it) {
513 _arborescence->set(it, false);
514 (*_arc_order)[it] = -1;
516 _dual_node_list.clear();
517 _dual_variables.clear();
520 /// \brief Adds a new source node.
522 /// Adds a new source node to the algorithm.
523 void addSource(Node source) {
524 std::vector<Node> nodes;
525 nodes.push_back(source);
526 while (!nodes.empty()) {
527 Node node = nodes.back();
529 for (OutArcIt it(*_digraph, node); it != INVALID; ++it) {
530 Node target = _digraph->target(it);
531 if ((*_node_order)[target] == -3) {
532 (*_node_order)[target] = -2;
533 nodes.push_back(target);
534 queue.push_back(target);
538 (*_node_order)[source] = -1;
541 /// \brief Processes the next node in the priority queue.
543 /// Processes the next node in the priority queue.
545 /// \return The processed node.
547 /// \warning The queue must not be empty.
548 Node processNextNode() {
549 Node node = queue.back();
551 if ((*_node_order)[node] == -2) {
552 Arc arc = prepare(node);
553 Node source = _digraph->source(arc);
554 while ((*_node_order)[source] != -1) {
555 if ((*_node_order)[source] >= 0) {
556 arc = contract(source);
558 arc = prepare(source);
560 source = _digraph->source(arc);
568 /// \brief Returns the number of the nodes to be processed.
570 /// Returns the number of the nodes to be processed in the priority
572 int queueSize() const {
576 /// \brief Returns \c false if there are nodes to be processed.
578 /// Returns \c false if there are nodes to be processed.
579 bool emptyQueue() const {
580 return queue.empty();
583 /// \brief Executes the algorithm.
585 /// Executes the algorithm.
587 /// \pre init() must be called and at least one node should be added
588 /// with addSource() before using this function.
590 ///\note mca.start() is just a shortcut of the following code.
592 ///while (!mca.emptyQueue()) {
593 /// mca.processNextNode();
597 while (!emptyQueue()) {
602 /// \brief Runs %MinCostArborescence algorithm from node \c s.
604 /// This method runs the %MinCostArborescence algorithm from
605 /// a root node \c s.
607 /// \note mca.run(s) is just a shortcut of the following code.
610 /// mca.addSource(s);
621 /// \name Query Functions
622 /// The result of the %MinCostArborescence algorithm can be obtained
623 /// using these functions.\n
624 /// Either run() or start() must be called before using them.
628 /// \brief Returns the cost of the arborescence.
630 /// Returns the cost of the arborescence.
631 Value arborescenceCost() const {
633 for (ArcIt it(*_digraph); it != INVALID; ++it) {
634 if (arborescence(it)) {
641 /// \brief Returns \c true if the arc is in the arborescence.
643 /// Returns \c true if the given arc is in the arborescence.
644 /// \param arc An arc of the digraph.
645 /// \pre \ref run() must be called before using this function.
646 bool arborescence(Arc arc) const {
647 return (*_pred)[_digraph->target(arc)] == arc;
650 /// \brief Returns a const reference to the arborescence map.
652 /// Returns a const reference to the arborescence map.
653 /// \pre \ref run() must be called before using this function.
654 const ArborescenceMap& arborescenceMap() const {
655 return *_arborescence;
658 /// \brief Returns the predecessor arc of the given node.
660 /// Returns the predecessor arc of the given node.
661 /// \pre \ref run() must be called before using this function.
662 Arc pred(Node node) const {
663 return (*_pred)[node];
666 /// \brief Returns a const reference to the pred map.
668 /// Returns a const reference to the pred map.
669 /// \pre \ref run() must be called before using this function.
670 const PredMap& predMap() const {
674 /// \brief Indicates that a node is reachable from the sources.
676 /// Indicates that a node is reachable from the sources.
677 bool reached(Node node) const {
678 return (*_node_order)[node] != -3;
681 /// \brief Indicates that a node is processed.
683 /// Indicates that a node is processed. The arborescence path exists
684 /// from the source to the given node.
685 bool processed(Node node) const {
686 return (*_node_order)[node] == -1;
689 /// \brief Returns the number of the dual variables in basis.
691 /// Returns the number of the dual variables in basis.
692 int dualNum() const {
693 return _dual_variables.size();
696 /// \brief Returns the value of the dual solution.
698 /// Returns the value of the dual solution. It should be
699 /// equal to the arborescence value.
700 Value dualValue() const {
702 for (int i = 0; i < int(_dual_variables.size()); ++i) {
703 sum += _dual_variables[i].value;
708 /// \brief Returns the number of the nodes in the dual variable.
710 /// Returns the number of the nodes in the dual variable.
711 int dualSize(int k) const {
712 return _dual_variables[k].end - _dual_variables[k].begin;
715 /// \brief Returns the value of the dual variable.
717 /// Returns the the value of the dual variable.
718 Value dualValue(int k) const {
719 return _dual_variables[k].value;
722 /// \brief LEMON iterator for getting a dual variable.
724 /// This class provides a common style LEMON iterator for getting a
725 /// dual variable of \ref MinCostArborescence algorithm.
726 /// It iterates over a subset of the nodes.
730 /// \brief Constructor.
732 /// Constructor for getting the nodeset of the dual variable
733 /// of \ref MinCostArborescence algorithm.
734 DualIt(const MinCostArborescence& algorithm, int variable)
735 : _algorithm(&algorithm)
737 _index = _algorithm->_dual_variables[variable].begin;
738 _last = _algorithm->_dual_variables[variable].end;
741 /// \brief Conversion to \c Node.
743 /// Conversion to \c Node.
744 operator Node() const {
745 return _algorithm->_dual_node_list[_index];
748 /// \brief Increment operator.
750 /// Increment operator.
751 DualIt& operator++() {
756 /// \brief Validity checking
758 /// Checks whether the iterator is invalid.
759 bool operator==(Invalid) const {
760 return _index == _last;
763 /// \brief Validity checking
765 /// Checks whether the iterator is valid.
766 bool operator!=(Invalid) const {
767 return _index != _last;
771 const MinCostArborescence* _algorithm;
779 /// \ingroup spantree
781 /// \brief Function type interface for MinCostArborescence algorithm.
783 /// Function type interface for MinCostArborescence algorithm.
784 /// \param digraph The digraph the algorithm runs on.
785 /// \param cost An arc map storing the costs.
786 /// \param source The source node of the arborescence.
787 /// \retval arborescence An arc map with \c bool (or convertible) value
788 /// type that stores the arborescence.
789 /// \return The total cost of the arborescence.
791 /// \sa MinCostArborescence
792 template <typename Digraph, typename CostMap, typename ArborescenceMap>
793 typename CostMap::Value minCostArborescence(const Digraph& digraph,
795 typename Digraph::Node source,
796 ArborescenceMap& arborescence) {
797 typename MinCostArborescence<Digraph, CostMap>
798 ::template SetArborescenceMap<ArborescenceMap>
799 ::Create mca(digraph, cost);
800 mca.arborescenceMap(arborescence);
802 return mca.arborescenceCost();