lemon/min_cost_arborescence.h
author Balazs Dezso <deba@inf.elte.hu>
Thu, 04 Mar 2010 15:20:59 +0100
changeset 951 41d7ac528c3a
parent 631 33c6b6e755cd
child 760 4ac30454f1c1
child 1081 f1398882a928
permissions -rw-r--r--
Uniforming primal scale to 2 (#314)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2008
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     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.
    12  *
    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
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_MIN_COST_ARBORESCENCE_H
    20 #define LEMON_MIN_COST_ARBORESCENCE_H
    21 
    22 ///\ingroup spantree
    23 ///\file
    24 ///\brief Minimum Cost Arborescence algorithm.
    25 
    26 #include <vector>
    27 
    28 #include <lemon/list_graph.h>
    29 #include <lemon/bin_heap.h>
    30 #include <lemon/assert.h>
    31 
    32 namespace lemon {
    33 
    34 
    35   /// \brief Default traits class for MinCostArborescence class.
    36   ///
    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{
    42 
    43     /// \brief The digraph type the algorithm runs on.
    44     typedef GR Digraph;
    45 
    46     /// \brief The type of the map that stores the arc costs.
    47     ///
    48     /// The type of the map that stores the arc costs.
    49     /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
    50     typedef CM CostMap;
    51 
    52     /// \brief The value type of the costs.
    53     ///
    54     /// The value type of the costs.
    55     typedef typename CostMap::Value Value;
    56 
    57     /// \brief The type of the map that stores which arcs are in the
    58     /// arborescence.
    59     ///
    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;
    66 
    67     /// \brief Instantiates a \c ArborescenceMap.
    68     ///
    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);
    74     }
    75 
    76     /// \brief The type of the \c PredMap
    77     ///
    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;
    82 
    83     /// \brief Instantiates a \c PredMap.
    84     ///
    85     /// This function instantiates a \c PredMap.
    86     /// \param digraph The digraph to which we would like to define the
    87     /// \c PredMap.
    88     static PredMap *createPredMap(const Digraph &digraph){
    89       return new PredMap(digraph);
    90     }
    91 
    92   };
    93 
    94   /// \ingroup spantree
    95   ///
    96   /// \brief Minimum Cost Arborescence algorithm class.
    97   ///
    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).
   105   ///
   106   /// The algorithm also provides an optimal dual solution, therefore
   107   /// the optimality of the solution can be checked.
   108   ///
   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>".
   119 #ifndef DOXYGEN
   120   template <typename GR,
   121             typename CM = typename GR::template ArcMap<int>,
   122             typename TR =
   123               MinCostArborescenceDefaultTraits<GR, CM> >
   124 #else
   125   template <typename GR, typename CM, typedef TR>
   126 #endif
   127   class MinCostArborescence {
   128   public:
   129 
   130     /// \brief The \ref MinCostArborescenceDefaultTraits "traits class" 
   131     /// of the algorithm. 
   132     typedef TR Traits;
   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;
   143 
   144     typedef MinCostArborescence Create;
   145 
   146   private:
   147 
   148     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
   149 
   150     struct CostArc {
   151 
   152       Arc arc;
   153       Value value;
   154 
   155       CostArc() {}
   156       CostArc(Arc _arc, Value _value) : arc(_arc), value(_value) {}
   157 
   158     };
   159 
   160     const Digraph *_digraph;
   161     const CostMap *_cost;
   162 
   163     PredMap *_pred;
   164     bool local_pred;
   165 
   166     ArborescenceMap *_arborescence;
   167     bool local_arborescence;
   168 
   169     typedef typename Digraph::template ArcMap<int> ArcOrder;
   170     ArcOrder *_arc_order;
   171 
   172     typedef typename Digraph::template NodeMap<int> NodeOrder;
   173     NodeOrder *_node_order;
   174 
   175     typedef typename Digraph::template NodeMap<CostArc> CostArcMap;
   176     CostArcMap *_cost_arcs;
   177 
   178     struct StackLevel {
   179 
   180       std::vector<CostArc> arcs;
   181       int node_level;
   182 
   183     };
   184 
   185     std::vector<StackLevel> level_stack;
   186     std::vector<Node> queue;
   187 
   188     typedef std::vector<typename Digraph::Node> DualNodeList;
   189 
   190     DualNodeList _dual_node_list;
   191 
   192     struct DualVariable {
   193       int begin, end;
   194       Value value;
   195 
   196       DualVariable(int _begin, int _end, Value _value)
   197         : begin(_begin), end(_end), value(_value) {}
   198 
   199     };
   200 
   201     typedef std::vector<DualVariable> DualVariables;
   202 
   203     DualVariables _dual_variables;
   204 
   205     typedef typename Digraph::template NodeMap<int> HeapCrossRef;
   206 
   207     HeapCrossRef *_heap_cross_ref;
   208 
   209     typedef BinHeap<int, HeapCrossRef> Heap;
   210 
   211     Heap *_heap;
   212 
   213   protected:
   214 
   215     MinCostArborescence() {}
   216 
   217   private:
   218 
   219     void createStructures() {
   220       if (!_pred) {
   221         local_pred = true;
   222         _pred = Traits::createPredMap(*_digraph);
   223       }
   224       if (!_arborescence) {
   225         local_arborescence = true;
   226         _arborescence = Traits::createArborescenceMap(*_digraph);
   227       }
   228       if (!_arc_order) {
   229         _arc_order = new ArcOrder(*_digraph);
   230       }
   231       if (!_node_order) {
   232         _node_order = new NodeOrder(*_digraph);
   233       }
   234       if (!_cost_arcs) {
   235         _cost_arcs = new CostArcMap(*_digraph);
   236       }
   237       if (!_heap_cross_ref) {
   238         _heap_cross_ref = new HeapCrossRef(*_digraph, -1);
   239       }
   240       if (!_heap) {
   241         _heap = new Heap(*_heap_cross_ref);
   242       }
   243     }
   244 
   245     void destroyStructures() {
   246       if (local_arborescence) {
   247         delete _arborescence;
   248       }
   249       if (local_pred) {
   250         delete _pred;
   251       }
   252       if (_arc_order) {
   253         delete _arc_order;
   254       }
   255       if (_node_order) {
   256         delete _node_order;
   257       }
   258       if (_cost_arcs) {
   259         delete _cost_arcs;
   260       }
   261       if (_heap) {
   262         delete _heap;
   263       }
   264       if (_heap_cross_ref) {
   265         delete _heap_cross_ref;
   266       }
   267     }
   268 
   269     Arc prepare(Node node) {
   270       std::vector<Node> nodes;
   271       (*_node_order)[node] = _dual_node_list.size();
   272       StackLevel level;
   273       level.node_level = _dual_node_list.size();
   274       _dual_node_list.push_back(node);
   275       for (InArcIt it(*_digraph, node); it != INVALID; ++it) {
   276         Arc arc = 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);
   284         } else {
   285           if ((*_cost_arcs)[source].value > value) {
   286             (*_cost_arcs)[source].arc = arc;
   287             (*_cost_arcs)[source].value = value;
   288           }
   289         }
   290       }
   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]];
   295         }
   296       }
   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;
   305       }
   306       level_stack.push_back(level);
   307       return minimum.arc;
   308     }
   309 
   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);
   324           } else {
   325             if ((*_cost_arcs)[source].value > value) {
   326               (*_cost_arcs)[source].arc = arc;
   327               (*_cost_arcs)[source].value = value;
   328             }
   329           }
   330         }
   331         level_stack.pop_back();
   332       }
   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]];
   337         }
   338       }
   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);
   342       StackLevel level;
   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;
   348       }
   349       level_stack.push_back(level);
   350       return minimum.arc;
   351     }
   352 
   353     int bottom(Node node) {
   354       int k = level_stack.size() - 1;
   355       while (level_stack[k].node_level > (*_node_order)[node]) {
   356         --k;
   357       }
   358       return level_stack[k].node_level;
   359     }
   360 
   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();
   367         _heap->pop();
   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)) {
   373           case Heap::PRE_HEAP:
   374             _heap->push(target, (*_arc_order)[it]);
   375             _pred->set(target, it);
   376             break;
   377           case Heap::IN_HEAP:
   378             if ((*_arc_order)[it] < (*_heap)[target]) {
   379               _heap->decrease(target, (*_arc_order)[it]);
   380               _pred->set(target, it);
   381             }
   382             break;
   383           case Heap::POST_HEAP:
   384             break;
   385           }
   386         }
   387         _arborescence->set((*_pred)[source], true);
   388       }
   389     }
   390 
   391 
   392   public:
   393 
   394     /// \name Named Template Parameters
   395 
   396     /// @{
   397 
   398     template <class T>
   399     struct SetArborescenceMapTraits : public Traits {
   400       typedef T ArborescenceMap;
   401       static ArborescenceMap *createArborescenceMap(const Digraph &)
   402       {
   403         LEMON_ASSERT(false, "ArborescenceMap is not initialized");
   404         return 0; // ignore warnings
   405       }
   406     };
   407 
   408     /// \brief \ref named-templ-param "Named parameter" for
   409     /// setting \c ArborescenceMap type
   410     ///
   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.
   417     template <class T>
   418     struct SetArborescenceMap
   419       : public MinCostArborescence<Digraph, CostMap,
   420                                    SetArborescenceMapTraits<T> > {
   421     };
   422 
   423     template <class T>
   424     struct SetPredMapTraits : public Traits {
   425       typedef T PredMap;
   426       static PredMap *createPredMap(const Digraph &)
   427       {
   428         LEMON_ASSERT(false, "PredMap is not initialized");
   429         return 0; // ignore warnings
   430       }
   431     };
   432 
   433     /// \brief \ref named-templ-param "Named parameter" for
   434     /// setting \c PredMap type
   435     ///
   436     /// \ref named-templ-param "Named parameter" for setting
   437     /// \c PredMap type.
   438     /// It must meet the \ref concepts::WriteMap "WriteMap" concept, 
   439     /// and its value type must be the \c Arc type of the digraph.
   440     template <class T>
   441     struct SetPredMap
   442       : public MinCostArborescence<Digraph, CostMap, SetPredMapTraits<T> > {
   443     };
   444 
   445     /// @}
   446 
   447     /// \brief Constructor.
   448     ///
   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) {}
   456 
   457     /// \brief Destructor.
   458     ~MinCostArborescence() {
   459       destroyStructures();
   460     }
   461 
   462     /// \brief Sets the arborescence map.
   463     ///
   464     /// Sets the arborescence map.
   465     /// \return <tt>(*this)</tt>
   466     MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
   467       if (local_arborescence) {
   468         delete _arborescence;
   469       }
   470       local_arborescence = false;
   471       _arborescence = &m;
   472       return *this;
   473     }
   474 
   475     /// \brief Sets the predecessor map.
   476     ///
   477     /// Sets the predecessor map.
   478     /// \return <tt>(*this)</tt>
   479     MinCostArborescence& predMap(PredMap& m) {
   480       if (local_pred) {
   481         delete _pred;
   482       }
   483       local_pred = false;
   484       _pred = &m;
   485       return *this;
   486     }
   487 
   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
   495     /// computation.
   496 
   497     ///@{
   498 
   499     /// \brief Initializes the internal data structures.
   500     ///
   501     /// Initializes the internal data structures.
   502     ///
   503     void init() {
   504       createStructures();
   505       _heap->clear();
   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);
   511       }
   512       for (ArcIt it(*_digraph); it != INVALID; ++it) {
   513         _arborescence->set(it, false);
   514         (*_arc_order)[it] = -1;
   515       }
   516       _dual_node_list.clear();
   517       _dual_variables.clear();
   518     }
   519 
   520     /// \brief Adds a new source node.
   521     ///
   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();
   528         nodes.pop_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);
   535           }
   536         }
   537       }
   538       (*_node_order)[source] = -1;
   539     }
   540 
   541     /// \brief Processes the next node in the priority queue.
   542     ///
   543     /// Processes the next node in the priority queue.
   544     ///
   545     /// \return The processed node.
   546     ///
   547     /// \warning The queue must not be empty.
   548     Node processNextNode() {
   549       Node node = queue.back();
   550       queue.pop_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);
   557           } else {
   558             arc = prepare(source);
   559           }
   560           source = _digraph->source(arc);
   561         }
   562         finalize(arc);
   563         level_stack.clear();
   564       }
   565       return node;
   566     }
   567 
   568     /// \brief Returns the number of the nodes to be processed.
   569     ///
   570     /// Returns the number of the nodes to be processed in the priority
   571     /// queue.
   572     int queueSize() const {
   573       return queue.size();
   574     }
   575 
   576     /// \brief Returns \c false if there are nodes to be processed.
   577     ///
   578     /// Returns \c false if there are nodes to be processed.
   579     bool emptyQueue() const {
   580       return queue.empty();
   581     }
   582 
   583     /// \brief Executes the algorithm.
   584     ///
   585     /// Executes the algorithm.
   586     ///
   587     /// \pre init() must be called and at least one node should be added
   588     /// with addSource() before using this function.
   589     ///
   590     ///\note mca.start() is just a shortcut of the following code.
   591     ///\code
   592     ///while (!mca.emptyQueue()) {
   593     ///  mca.processNextNode();
   594     ///}
   595     ///\endcode
   596     void start() {
   597       while (!emptyQueue()) {
   598         processNextNode();
   599       }
   600     }
   601 
   602     /// \brief Runs %MinCostArborescence algorithm from node \c s.
   603     ///
   604     /// This method runs the %MinCostArborescence algorithm from
   605     /// a root node \c s.
   606     ///
   607     /// \note mca.run(s) is just a shortcut of the following code.
   608     /// \code
   609     /// mca.init();
   610     /// mca.addSource(s);
   611     /// mca.start();
   612     /// \endcode
   613     void run(Node s) {
   614       init();
   615       addSource(s);
   616       start();
   617     }
   618 
   619     ///@}
   620 
   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.
   625 
   626     /// @{
   627 
   628     /// \brief Returns the cost of the arborescence.
   629     ///
   630     /// Returns the cost of the arborescence.
   631     Value arborescenceCost() const {
   632       Value sum = 0;
   633       for (ArcIt it(*_digraph); it != INVALID; ++it) {
   634         if (arborescence(it)) {
   635           sum += (*_cost)[it];
   636         }
   637       }
   638       return sum;
   639     }
   640 
   641     /// \brief Returns \c true if the arc is in the arborescence.
   642     ///
   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;
   648     }
   649 
   650     /// \brief Returns a const reference to the arborescence map.
   651     ///
   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;
   656     }
   657 
   658     /// \brief Returns the predecessor arc of the given node.
   659     ///
   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];
   664     }
   665 
   666     /// \brief Returns a const reference to the pred map.
   667     ///
   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 {
   671       return *_pred;
   672     }
   673 
   674     /// \brief Indicates that a node is reachable from the sources.
   675     ///
   676     /// Indicates that a node is reachable from the sources.
   677     bool reached(Node node) const {
   678       return (*_node_order)[node] != -3;
   679     }
   680 
   681     /// \brief Indicates that a node is processed.
   682     ///
   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;
   687     }
   688 
   689     /// \brief Returns the number of the dual variables in basis.
   690     ///
   691     /// Returns the number of the dual variables in basis.
   692     int dualNum() const {
   693       return _dual_variables.size();
   694     }
   695 
   696     /// \brief Returns the value of the dual solution.
   697     ///
   698     /// Returns the value of the dual solution. It should be
   699     /// equal to the arborescence value.
   700     Value dualValue() const {
   701       Value sum = 0;
   702       for (int i = 0; i < int(_dual_variables.size()); ++i) {
   703         sum += _dual_variables[i].value;
   704       }
   705       return sum;
   706     }
   707 
   708     /// \brief Returns the number of the nodes in the dual variable.
   709     ///
   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;
   713     }
   714 
   715     /// \brief Returns the value of the dual variable.
   716     ///
   717     /// Returns the the value of the dual variable.
   718     Value dualValue(int k) const {
   719       return _dual_variables[k].value;
   720     }
   721 
   722     /// \brief LEMON iterator for getting a dual variable.
   723     ///
   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.
   727     class DualIt {
   728     public:
   729 
   730       /// \brief Constructor.
   731       ///
   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)
   736       {
   737         _index = _algorithm->_dual_variables[variable].begin;
   738         _last = _algorithm->_dual_variables[variable].end;
   739       }
   740 
   741       /// \brief Conversion to \c Node.
   742       ///
   743       /// Conversion to \c Node.
   744       operator Node() const {
   745         return _algorithm->_dual_node_list[_index];
   746       }
   747 
   748       /// \brief Increment operator.
   749       ///
   750       /// Increment operator.
   751       DualIt& operator++() {
   752         ++_index;
   753         return *this;
   754       }
   755 
   756       /// \brief Validity checking
   757       ///
   758       /// Checks whether the iterator is invalid.
   759       bool operator==(Invalid) const {
   760         return _index == _last;
   761       }
   762 
   763       /// \brief Validity checking
   764       ///
   765       /// Checks whether the iterator is valid.
   766       bool operator!=(Invalid) const {
   767         return _index != _last;
   768       }
   769 
   770     private:
   771       const MinCostArborescence* _algorithm;
   772       int _index, _last;
   773     };
   774 
   775     /// @}
   776 
   777   };
   778 
   779   /// \ingroup spantree
   780   ///
   781   /// \brief Function type interface for MinCostArborescence algorithm.
   782   ///
   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.
   790   ///
   791   /// \sa MinCostArborescence
   792   template <typename Digraph, typename CostMap, typename ArborescenceMap>
   793   typename CostMap::Value minCostArborescence(const Digraph& digraph,
   794                                               const CostMap& cost,
   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);
   801     mca.run(source);
   802     return mca.arborescenceCost();
   803   }
   804 
   805 }
   806 
   807 #endif