# HG changeset patch
# User Balazs Dezso <deba@inf.elte.hu>
# Date 2008-12-02 23:33:47
# Node ID 7f8560cb9d650ee9445d399eb3a7b3d1c05eb931
# Parent  ad483acf1654bd4d0c093b807b858ed59529f78e

Port MinCostArborescence algorithm from SVN #3509

diff --git a/lemon/Makefile.am b/lemon/Makefile.am
--- a/lemon/Makefile.am
+++ b/lemon/Makefile.am
@@ -44,6 +44,7 @@
 	lemon/maps.h \
 	lemon/math.h \
 	lemon/max_matching.h \
+	lemon/min_cost_arborescence.h \
 	lemon/nauty_reader.h \
 	lemon/path.h \
 	lemon/preflow.h \
diff --git a/lemon/min_cost_arborescence.h b/lemon/min_cost_arborescence.h
new file mode 100644
--- /dev/null
+++ b/lemon/min_cost_arborescence.h
@@ -0,0 +1,796 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_MIN_COST_ARBORESCENCE_H
+#define LEMON_MIN_COST_ARBORESCENCE_H
+
+///\ingroup spantree
+///\file
+///\brief Minimum Cost Arborescence algorithm.
+
+#include <vector>
+
+#include <lemon/list_graph.h>
+#include <lemon/bin_heap.h>
+#include <lemon/assert.h>
+
+namespace lemon {
+
+
+  /// \brief Default traits class for MinCostArborescence class.
+  ///
+  /// Default traits class for MinCostArborescence class.
+  /// \param _Digraph Digraph type.
+  /// \param _CostMap Type of cost map.
+  template <class _Digraph, class _CostMap>
+  struct MinCostArborescenceDefaultTraits{
+
+    /// \brief The digraph type the algorithm runs on.
+    typedef _Digraph Digraph;
+
+    /// \brief The type of the map that stores the arc costs.
+    ///
+    /// The type of the map that stores the arc costs.
+    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
+    typedef _CostMap CostMap;
+
+    /// \brief The value type of the costs.
+    ///
+    /// The value type of the costs.
+    typedef typename CostMap::Value Value;
+
+    /// \brief The type of the map that stores which arcs are in the
+    /// arborescence.
+    ///
+    /// The type of the map that stores which arcs are in the
+    /// arborescence.  It must meet the \ref concepts::WriteMap
+    /// "WriteMap" concept.  Initially it will be set to false on each
+    /// arc. After it will set all arborescence arcs once.
+    typedef typename Digraph::template ArcMap<bool> ArborescenceMap;
+
+    /// \brief Instantiates a ArborescenceMap.
+    ///
+    /// This function instantiates a \ref ArborescenceMap.
+    /// \param digraph is the graph, to which we would like to
+    /// calculate the ArborescenceMap.
+    static ArborescenceMap *createArborescenceMap(const Digraph &digraph){
+      return new ArborescenceMap(digraph);
+    }
+
+    /// \brief The type of the PredMap
+    ///
+    /// The type of the PredMap. It is a node map with an arc value type.
+    typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
+
+    /// \brief Instantiates a PredMap.
+    ///
+    /// This function instantiates a \ref PredMap.
+    /// \param _digraph is the digraph, to which we would like to define the
+    /// PredMap.
+    static PredMap *createPredMap(const Digraph &digraph){
+      return new PredMap(digraph);
+    }
+
+  };
+
+  /// \ingroup spantree
+  ///
+  /// \brief %MinCostArborescence algorithm class.
+  ///
+  /// This class provides an efficient implementation of
+  /// %MinCostArborescence algorithm. The arborescence is a tree
+  /// which is directed from a given source node of the digraph. One or
+  /// more sources should be given for the algorithm and it will calculate
+  /// the minimum cost subgraph which are union of arborescences with the
+  /// given sources and spans all the nodes which are reachable from the
+  /// sources. The time complexity of the algorithm is \f$ O(n^2+e) \f$.
+  ///
+  /// The algorithm provides also an optimal dual solution, therefore
+  /// the optimality of the solution can be checked.
+  ///
+  /// \param _Digraph The digraph type the algorithm runs on. The default value
+  /// is \ref ListDigraph.
+  /// \param _CostMap This read-only ArcMap determines the costs of the
+  /// arcs. It is read once for each arc, so the map may involve in
+  /// relatively time consuming process to compute the arc cost if
+  /// it is necessary. The default map type is \ref
+  /// concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
+  /// \param _Traits Traits class to set various data types used
+  /// by the algorithm. The default traits class is
+  /// \ref MinCostArborescenceDefaultTraits
+  /// "MinCostArborescenceDefaultTraits<_Digraph, _CostMap>".  See \ref
+  /// MinCostArborescenceDefaultTraits for the documentation of a
+  /// MinCostArborescence traits class.
+  ///
+  /// \author Balazs Dezso
+#ifndef DOXYGEN
+  template <typename _Digraph = ListDigraph,
+            typename _CostMap = typename _Digraph::template ArcMap<int>,
+            typename _Traits =
+            MinCostArborescenceDefaultTraits<_Digraph, _CostMap> >
+#else
+  template <typename _Digraph, typename _CostMap, typedef _Traits>
+#endif
+  class MinCostArborescence {
+  public:
+
+    /// The traits.
+    typedef _Traits Traits;
+    /// The type of the underlying digraph.
+    typedef typename Traits::Digraph Digraph;
+    /// The type of the map that stores the arc costs.
+    typedef typename Traits::CostMap CostMap;
+    ///The type of the costs of the arcs.
+    typedef typename Traits::Value Value;
+    ///The type of the predecessor map.
+    typedef typename Traits::PredMap PredMap;
+    ///The type of the map that stores which arcs are in the arborescence.
+    typedef typename Traits::ArborescenceMap ArborescenceMap;
+
+    typedef MinCostArborescence Create;
+
+  private:
+
+    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+
+    struct CostArc {
+
+      Arc arc;
+      Value value;
+
+      CostArc() {}
+      CostArc(Arc _arc, Value _value) : arc(_arc), value(_value) {}
+
+    };
+
+    const Digraph *_digraph;
+    const CostMap *_cost;
+
+    PredMap *_pred;
+    bool local_pred;
+
+    ArborescenceMap *_arborescence;
+    bool local_arborescence;
+
+    typedef typename Digraph::template ArcMap<int> ArcOrder;
+    ArcOrder *_arc_order;
+
+    typedef typename Digraph::template NodeMap<int> NodeOrder;
+    NodeOrder *_node_order;
+
+    typedef typename Digraph::template NodeMap<CostArc> CostArcMap;
+    CostArcMap *_cost_arcs;
+
+    struct StackLevel {
+
+      std::vector<CostArc> arcs;
+      int node_level;
+
+    };
+
+    std::vector<StackLevel> level_stack;
+    std::vector<Node> queue;
+
+    typedef std::vector<typename Digraph::Node> DualNodeList;
+
+    DualNodeList _dual_node_list;
+
+    struct DualVariable {
+      int begin, end;
+      Value value;
+
+      DualVariable(int _begin, int _end, Value _value)
+        : begin(_begin), end(_end), value(_value) {}
+
+    };
+
+    typedef std::vector<DualVariable> DualVariables;
+
+    DualVariables _dual_variables;
+
+    typedef typename Digraph::template NodeMap<int> HeapCrossRef;
+
+    HeapCrossRef *_heap_cross_ref;
+
+    typedef BinHeap<int, HeapCrossRef> Heap;
+
+    Heap *_heap;
+
+  protected:
+
+    MinCostArborescence() {}
+
+  private:
+
+    void createStructures() {
+      if (!_pred) {
+        local_pred = true;
+        _pred = Traits::createPredMap(*_digraph);
+      }
+      if (!_arborescence) {
+        local_arborescence = true;
+        _arborescence = Traits::createArborescenceMap(*_digraph);
+      }
+      if (!_arc_order) {
+        _arc_order = new ArcOrder(*_digraph);
+      }
+      if (!_node_order) {
+        _node_order = new NodeOrder(*_digraph);
+      }
+      if (!_cost_arcs) {
+        _cost_arcs = new CostArcMap(*_digraph);
+      }
+      if (!_heap_cross_ref) {
+        _heap_cross_ref = new HeapCrossRef(*_digraph, -1);
+      }
+      if (!_heap) {
+        _heap = new Heap(*_heap_cross_ref);
+      }
+    }
+
+    void destroyStructures() {
+      if (local_arborescence) {
+        delete _arborescence;
+      }
+      if (local_pred) {
+        delete _pred;
+      }
+      if (_arc_order) {
+        delete _arc_order;
+      }
+      if (_node_order) {
+        delete _node_order;
+      }
+      if (_cost_arcs) {
+        delete _cost_arcs;
+      }
+      if (_heap) {
+        delete _heap;
+      }
+      if (_heap_cross_ref) {
+        delete _heap_cross_ref;
+      }
+    }
+
+    Arc prepare(Node node) {
+      std::vector<Node> nodes;
+      (*_node_order)[node] = _dual_node_list.size();
+      StackLevel level;
+      level.node_level = _dual_node_list.size();
+      _dual_node_list.push_back(node);
+      for (InArcIt it(*_digraph, node); it != INVALID; ++it) {
+        Arc arc = it;
+        Node source = _digraph->source(arc);
+        Value value = (*_cost)[it];
+        if (source == node || (*_node_order)[source] == -3) continue;
+        if ((*_cost_arcs)[source].arc == INVALID) {
+          (*_cost_arcs)[source].arc = arc;
+          (*_cost_arcs)[source].value = value;
+          nodes.push_back(source);
+        } else {
+          if ((*_cost_arcs)[source].value > value) {
+            (*_cost_arcs)[source].arc = arc;
+            (*_cost_arcs)[source].value = value;
+          }
+        }
+      }
+      CostArc minimum = (*_cost_arcs)[nodes[0]];
+      for (int i = 1; i < int(nodes.size()); ++i) {
+        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
+          minimum = (*_cost_arcs)[nodes[i]];
+        }
+      }
+      _arc_order->set(minimum.arc, _dual_variables.size());
+      DualVariable var(_dual_node_list.size() - 1,
+                       _dual_node_list.size(), minimum.value);
+      _dual_variables.push_back(var);
+      for (int i = 0; i < int(nodes.size()); ++i) {
+        (*_cost_arcs)[nodes[i]].value -= minimum.value;
+        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
+        (*_cost_arcs)[nodes[i]].arc = INVALID;
+      }
+      level_stack.push_back(level);
+      return minimum.arc;
+    }
+
+    Arc contract(Node node) {
+      int node_bottom = bottom(node);
+      std::vector<Node> nodes;
+      while (!level_stack.empty() &&
+             level_stack.back().node_level >= node_bottom) {
+        for (int i = 0; i < int(level_stack.back().arcs.size()); ++i) {
+          Arc arc = level_stack.back().arcs[i].arc;
+          Node source = _digraph->source(arc);
+          Value value = level_stack.back().arcs[i].value;
+          if ((*_node_order)[source] >= node_bottom) continue;
+          if ((*_cost_arcs)[source].arc == INVALID) {
+            (*_cost_arcs)[source].arc = arc;
+            (*_cost_arcs)[source].value = value;
+            nodes.push_back(source);
+          } else {
+            if ((*_cost_arcs)[source].value > value) {
+              (*_cost_arcs)[source].arc = arc;
+              (*_cost_arcs)[source].value = value;
+            }
+          }
+        }
+        level_stack.pop_back();
+      }
+      CostArc minimum = (*_cost_arcs)[nodes[0]];
+      for (int i = 1; i < int(nodes.size()); ++i) {
+        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
+          minimum = (*_cost_arcs)[nodes[i]];
+        }
+      }
+      _arc_order->set(minimum.arc, _dual_variables.size());
+      DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
+      _dual_variables.push_back(var);
+      StackLevel level;
+      level.node_level = node_bottom;
+      for (int i = 0; i < int(nodes.size()); ++i) {
+        (*_cost_arcs)[nodes[i]].value -= minimum.value;
+        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
+        (*_cost_arcs)[nodes[i]].arc = INVALID;
+      }
+      level_stack.push_back(level);
+      return minimum.arc;
+    }
+
+    int bottom(Node node) {
+      int k = level_stack.size() - 1;
+      while (level_stack[k].node_level > (*_node_order)[node]) {
+        --k;
+      }
+      return level_stack[k].node_level;
+    }
+
+    void finalize(Arc arc) {
+      Node node = _digraph->target(arc);
+      _heap->push(node, (*_arc_order)[arc]);
+      _pred->set(node, arc);
+      while (!_heap->empty()) {
+        Node source = _heap->top();
+        _heap->pop();
+        _node_order->set(source, -1);
+        for (OutArcIt it(*_digraph, source); it != INVALID; ++it) {
+          if ((*_arc_order)[it] < 0) continue;
+          Node target = _digraph->target(it);
+          switch(_heap->state(target)) {
+          case Heap::PRE_HEAP:
+            _heap->push(target, (*_arc_order)[it]);
+            _pred->set(target, it);
+            break;
+          case Heap::IN_HEAP:
+            if ((*_arc_order)[it] < (*_heap)[target]) {
+              _heap->decrease(target, (*_arc_order)[it]);
+              _pred->set(target, it);
+            }
+            break;
+          case Heap::POST_HEAP:
+            break;
+          }
+        }
+        _arborescence->set((*_pred)[source], true);
+      }
+    }
+
+
+  public:
+
+    /// \name Named template parameters
+
+    /// @{
+
+    template <class T>
+    struct DefArborescenceMapTraits : public Traits {
+      typedef T ArborescenceMap;
+      static ArborescenceMap *createArborescenceMap(const Digraph &)
+      {
+        LEMON_ASSERT(false, "ArborescenceMap is not initialized");
+        return 0; // ignore warnings
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for
+    /// setting ArborescenceMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting
+    /// ArborescenceMap type
+    template <class T>
+    struct DefArborescenceMap
+      : public MinCostArborescence<Digraph, CostMap,
+                                   DefArborescenceMapTraits<T> > {
+    };
+
+    template <class T>
+    struct DefPredMapTraits : public Traits {
+      typedef T PredMap;
+      static PredMap *createPredMap(const Digraph &)
+      {
+        LEMON_ASSERT(false, "PredMap is not initialized");
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for
+    /// setting PredMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting
+    /// PredMap type
+    template <class T>
+    struct DefPredMap
+      : public MinCostArborescence<Digraph, CostMap, DefPredMapTraits<T> > {
+    };
+
+    /// @}
+
+    /// \brief Constructor.
+    ///
+    /// \param _digraph The digraph the algorithm will run on.
+    /// \param _cost The cost map used by the algorithm.
+    MinCostArborescence(const Digraph& digraph, const CostMap& cost)
+      : _digraph(&digraph), _cost(&cost), _pred(0), local_pred(false),
+        _arborescence(0), local_arborescence(false),
+        _arc_order(0), _node_order(0), _cost_arcs(0),
+        _heap_cross_ref(0), _heap(0) {}
+
+    /// \brief Destructor.
+    ~MinCostArborescence() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the arborescence map.
+    ///
+    /// Sets the arborescence map.
+    /// \return \c (*this)
+    MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
+      if (local_arborescence) {
+        delete _arborescence;
+      }
+      local_arborescence = false;
+      _arborescence = &m;
+      return *this;
+    }
+
+    /// \brief Sets the arborescence map.
+    ///
+    /// Sets the arborescence map.
+    /// \return \c (*this)
+    MinCostArborescence& predMap(PredMap& m) {
+      if (local_pred) {
+        delete _pred;
+      }
+      local_pred = false;
+      _pred = &m;
+      return *this;
+    }
+
+    /// \name Query Functions
+    /// The result of the %MinCostArborescence algorithm can be obtained
+    /// using these functions.\n
+    /// Before the use of these functions,
+    /// either run() or start() must be called.
+
+    /// @{
+
+    /// \brief Returns a reference to the arborescence map.
+    ///
+    /// Returns a reference to the arborescence map.
+    const ArborescenceMap& arborescenceMap() const {
+      return *_arborescence;
+    }
+
+    /// \brief Returns true if the arc is in the arborescence.
+    ///
+    /// Returns true if the arc is in the arborescence.
+    /// \param arc The arc of the digraph.
+    /// \pre \ref run() must be called before using this function.
+    bool arborescence(Arc arc) const {
+      return (*_pred)[_digraph->target(arc)] == arc;
+    }
+
+    /// \brief Returns a reference to the pred map.
+    ///
+    /// Returns a reference to the pred map.
+    const PredMap& predMap() const {
+      return *_pred;
+    }
+
+    /// \brief Returns the predecessor arc of the given node.
+    ///
+    /// Returns the predecessor arc of the given node.
+    Arc pred(Node node) const {
+      return (*_pred)[node];
+    }
+
+    /// \brief Returns the cost of the arborescence.
+    ///
+    /// Returns the cost of the arborescence.
+    Value arborescenceValue() const {
+      Value sum = 0;
+      for (ArcIt it(*_digraph); it != INVALID; ++it) {
+        if (arborescence(it)) {
+          sum += (*_cost)[it];
+        }
+      }
+      return sum;
+    }
+
+    /// \brief Indicates that a node is reachable from the sources.
+    ///
+    /// Indicates that a node is reachable from the sources.
+    bool reached(Node node) const {
+      return (*_node_order)[node] != -3;
+    }
+
+    /// \brief Indicates that a node is processed.
+    ///
+    /// Indicates that a node is processed. The arborescence path exists
+    /// from the source to the given node.
+    bool processed(Node node) const {
+      return (*_node_order)[node] == -1;
+    }
+
+    /// \brief Returns the number of the dual variables in basis.
+    ///
+    /// Returns the number of the dual variables in basis.
+    int dualNum() const {
+      return _dual_variables.size();
+    }
+
+    /// \brief Returns the value of the dual solution.
+    ///
+    /// Returns the value of the dual solution. It should be
+    /// equal to the arborescence value.
+    Value dualValue() const {
+      Value sum = 0;
+      for (int i = 0; i < int(_dual_variables.size()); ++i) {
+        sum += _dual_variables[i].value;
+      }
+      return sum;
+    }
+
+    /// \brief Returns the number of the nodes in the dual variable.
+    ///
+    /// Returns the number of the nodes in the dual variable.
+    int dualSize(int k) const {
+      return _dual_variables[k].end - _dual_variables[k].begin;
+    }
+
+    /// \brief Returns the value of the dual variable.
+    ///
+    /// Returns the the value of the dual variable.
+    const Value& dualValue(int k) const {
+      return _dual_variables[k].value;
+    }
+
+    /// \brief Lemon iterator for get a dual variable.
+    ///
+    /// Lemon iterator for get a dual variable. This class provides
+    /// a common style lemon iterator which gives back a subset of
+    /// the nodes.
+    class DualIt {
+    public:
+
+      /// \brief Constructor.
+      ///
+      /// Constructor for get the nodeset of the variable.
+      DualIt(const MinCostArborescence& algorithm, int variable)
+        : _algorithm(&algorithm)
+      {
+        _index = _algorithm->_dual_variables[variable].begin;
+        _last = _algorithm->_dual_variables[variable].end;
+      }
+
+      /// \brief Conversion to node.
+      ///
+      /// Conversion to node.
+      operator Node() const {
+        return _algorithm->_dual_node_list[_index];
+      }
+
+      /// \brief Increment operator.
+      ///
+      /// Increment operator.
+      DualIt& operator++() {
+        ++_index;
+        return *this;
+      }
+
+      /// \brief Validity checking
+      ///
+      /// Checks whether the iterator is invalid.
+      bool operator==(Invalid) const {
+        return _index == _last;
+      }
+
+      /// \brief Validity checking
+      ///
+      /// Checks whether the iterator is valid.
+      bool operator!=(Invalid) const {
+        return _index != _last;
+      }
+
+    private:
+      const MinCostArborescence* _algorithm;
+      int _index, _last;
+    };
+
+    /// @}
+
+    /// \name Execution control
+    /// The simplest way to execute the algorithm is to use
+    /// one of the member functions called \c run(...). \n
+    /// If you need more control on the execution,
+    /// first you must call \ref init(), then you can add several
+    /// source nodes with \ref addSource().
+    /// Finally \ref start() will perform the arborescence
+    /// computation.
+
+    ///@{
+
+    /// \brief Initializes the internal data structures.
+    ///
+    /// Initializes the internal data structures.
+    ///
+    void init() {
+      createStructures();
+      _heap->clear();
+      for (NodeIt it(*_digraph); it != INVALID; ++it) {
+        (*_cost_arcs)[it].arc = INVALID;
+        _node_order->set(it, -3);
+        _heap_cross_ref->set(it, Heap::PRE_HEAP);
+        _pred->set(it, INVALID);
+      }
+      for (ArcIt it(*_digraph); it != INVALID; ++it) {
+        _arborescence->set(it, false);
+        _arc_order->set(it, -1);
+      }
+      _dual_node_list.clear();
+      _dual_variables.clear();
+    }
+
+    /// \brief Adds a new source node.
+    ///
+    /// Adds a new source node to the algorithm.
+    void addSource(Node source) {
+      std::vector<Node> nodes;
+      nodes.push_back(source);
+      while (!nodes.empty()) {
+        Node node = nodes.back();
+        nodes.pop_back();
+        for (OutArcIt it(*_digraph, node); it != INVALID; ++it) {
+          Node target = _digraph->target(it);
+          if ((*_node_order)[target] == -3) {
+            (*_node_order)[target] = -2;
+            nodes.push_back(target);
+            queue.push_back(target);
+          }
+        }
+      }
+      (*_node_order)[source] = -1;
+    }
+
+    /// \brief Processes the next node in the priority queue.
+    ///
+    /// Processes the next node in the priority queue.
+    ///
+    /// \return The processed node.
+    ///
+    /// \warning The queue must not be empty!
+    Node processNextNode() {
+      Node node = queue.back();
+      queue.pop_back();
+      if ((*_node_order)[node] == -2) {
+        Arc arc = prepare(node);
+        Node source = _digraph->source(arc);
+        while ((*_node_order)[source] != -1) {
+          if ((*_node_order)[source] >= 0) {
+            arc = contract(source);
+          } else {
+            arc = prepare(source);
+          }
+          source = _digraph->source(arc);
+        }
+        finalize(arc);
+        level_stack.clear();
+      }
+      return node;
+    }
+
+    /// \brief Returns the number of the nodes to be processed.
+    ///
+    /// Returns the number of the nodes to be processed.
+    int queueSize() const {
+      return queue.size();
+    }
+
+    /// \brief Returns \c false if there are nodes to be processed.
+    ///
+    /// Returns \c false if there are nodes to be processed.
+    bool emptyQueue() const {
+      return queue.empty();
+    }
+
+    /// \brief Executes the algorithm.
+    ///
+    /// Executes the algorithm.
+    ///
+    /// \pre init() must be called and at least one node should be added
+    /// with addSource() before using this function.
+    ///
+    ///\note mca.start() is just a shortcut of the following code.
+    ///\code
+    ///while (!mca.emptyQueue()) {
+    ///  mca.processNextNode();
+    ///}
+    ///\endcode
+    void start() {
+      while (!emptyQueue()) {
+        processNextNode();
+      }
+    }
+
+    /// \brief Runs %MinCostArborescence algorithm from node \c s.
+    ///
+    /// This method runs the %MinCostArborescence algorithm from
+    /// a root node \c s.
+    ///
+    /// \note mca.run(s) is just a shortcut of the following code.
+    /// \code
+    /// mca.init();
+    /// mca.addSource(s);
+    /// mca.start();
+    /// \endcode
+    void run(Node node) {
+      init();
+      addSource(node);
+      start();
+    }
+
+    ///@}
+
+  };
+
+  /// \ingroup spantree
+  ///
+  /// \brief Function type interface for MinCostArborescence algorithm.
+  ///
+  /// Function type interface for MinCostArborescence algorithm.
+  /// \param digraph The Digraph that the algorithm runs on.
+  /// \param cost The CostMap of the arcs.
+  /// \param source The source of the arborescence.
+  /// \retval arborescence The bool ArcMap which stores the arborescence.
+  /// \return The cost of the arborescence.
+  ///
+  /// \sa MinCostArborescence
+  template <typename Digraph, typename CostMap, typename ArborescenceMap>
+  typename CostMap::Value minCostArborescence(const Digraph& digraph,
+                                              const CostMap& cost,
+                                              typename Digraph::Node source,
+                                              ArborescenceMap& arborescence) {
+    typename MinCostArborescence<Digraph, CostMap>
+      ::template DefArborescenceMap<ArborescenceMap>
+      ::Create mca(digraph, cost);
+    mca.arborescenceMap(arborescence);
+    mca.run(source);
+    return mca.arborescenceValue();
+  }
+
+}
+
+#endif
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -18,6 +18,7 @@
   kruskal_test
   maps_test
   max_matching_test
+  min_cost_arborescence_test
   random_test
   path_test
   time_measure_test
diff --git a/test/Makefile.am b/test/Makefile.am
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -25,6 +25,7 @@
 	test/hao_orlin_test \
         test/maps_test \
 	test/max_matching_test \
+	test/min_cost_arborescence_test \
         test/random_test \
         test/path_test \
         test/preflow_test \
@@ -54,6 +55,7 @@
 test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
 test_maps_test_SOURCES = test/maps_test.cc
 test_max_matching_test_SOURCES = test/max_matching_test.cc
+test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
 test_path_test_SOURCES = test/path_test.cc
 test_preflow_test_SOURCES = test/preflow_test.cc
 test_suurballe_test_SOURCES = test/suurballe_test.cc
diff --git a/test/min_cost_arborescence_test.cc b/test/min_cost_arborescence_test.cc
new file mode 100644
--- /dev/null
+++ b/test/min_cost_arborescence_test.cc
@@ -0,0 +1,146 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2008
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include <iostream>
+#include <set>
+#include <vector>
+#include <iterator>
+
+#include <lemon/smart_graph.h>
+#include <lemon/min_cost_arborescence.h>
+#include <lemon/lgf_reader.h>
+
+#include "test_tools.h"
+
+using namespace lemon;
+using namespace std;
+
+const char test_lgf[] =
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "8\n"
+  "9\n"
+  "@arcs\n"
+  "     label  cost\n"
+  "1 8  0      107\n"
+  "0 3  1      70\n"
+  "2 1  2      46\n"
+  "4 1  3      28\n"
+  "4 4  4      91\n"
+  "3 9  5      76\n"
+  "9 8  6      61\n"
+  "8 1  7      39\n"
+  "9 8  8      74\n"
+  "8 0  9      39\n"
+  "4 3  10     45\n"
+  "2 2  11     34\n"
+  "0 1  12     100\n"
+  "6 3  13     95\n"
+  "4 1  14     22\n"
+  "1 1  15     31\n"
+  "7 2  16     51\n"
+  "2 6  17     29\n"
+  "8 3  18     115\n"
+  "6 9  19     32\n"
+  "1 1  20     60\n"
+  "0 3  21     40\n"
+  "@attributes\n"
+  "source 0\n";
+
+int main() {
+  typedef SmartDigraph Digraph;
+  DIGRAPH_TYPEDEFS(Digraph);
+
+  typedef Digraph::ArcMap<double> CostMap;
+
+  Digraph digraph;
+  CostMap cost(digraph);
+  Node source;
+
+  std::istringstream is(test_lgf);
+  digraphReader(digraph, is).
+    arcMap("cost", cost).
+    node("source", source).run();
+
+  MinCostArborescence<Digraph, CostMap> mca(digraph, cost);
+  mca.run(source);
+
+  vector<pair<double, set<Node> > > dualSolution(mca.dualNum());
+
+  for (int i = 0; i < mca.dualNum(); ++i) {
+    dualSolution[i].first = mca.dualValue(i);
+    for (MinCostArborescence<Digraph, CostMap>::DualIt it(mca, i);
+         it != INVALID; ++it) {
+      dualSolution[i].second.insert(it);
+    }
+  }
+
+  for (ArcIt it(digraph); it != INVALID; ++it) {
+    if (mca.reached(digraph.source(it))) {
+      double sum = 0.0;
+      for (int i = 0; i < int(dualSolution.size()); ++i) {
+        if (dualSolution[i].second.find(digraph.target(it))
+            != dualSolution[i].second.end() &&
+            dualSolution[i].second.find(digraph.source(it))
+            == dualSolution[i].second.end()) {
+          sum += dualSolution[i].first;
+        }
+      }
+      if (mca.arborescence(it)) {
+        check(sum == cost[it], "INVALID DUAL");
+      }
+      check(sum <= cost[it], "INVALID DUAL");
+    }
+  }
+
+
+  check(mca.dualValue() == mca.arborescenceValue(), "INVALID DUAL");
+
+  check(mca.reached(source), "INVALID ARBORESCENCE");
+  for (ArcIt a(digraph); a != INVALID; ++a) {
+    check(!mca.reached(digraph.source(a)) ||
+          mca.reached(digraph.target(a)), "INVALID ARBORESCENCE");
+  }
+
+  for (NodeIt n(digraph); n != INVALID; ++n) {
+    if (!mca.reached(n)) continue;
+    int cnt = 0;
+    for (InArcIt a(digraph, n); a != INVALID; ++a) {
+      if (mca.arborescence(a)) {
+        check(mca.pred(n) == a, "INVALID ARBORESCENCE");
+        ++cnt;
+      }
+    }
+    check((n == source ? cnt == 0 : cnt == 1), "INVALID ARBORESCENCE");
+  }
+
+  Digraph::ArcMap<bool> arborescence(digraph);
+  check(mca.arborescenceValue() ==
+        minCostArborescence(digraph, cost, source, arborescence),
+        "WRONG FUNCTION");
+
+  return 0;
+}