diff -r c6acc34f98dc -r 1a7fe3bef514 lemon/min_cost_arborescence.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/min_cost_arborescence.h	Thu Nov 05 15:50:01 2009 +0100
@@ -0,0 +1,807 @@
+/* -*- 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 GR Digraph type.
+  /// \param CM Type of the cost map.
+  template <class GR, class CM>
+  struct MinCostArborescenceDefaultTraits{
+
+    /// \brief The digraph type the algorithm runs on.
+    typedef GR Digraph;
+
+    /// \brief The type of the map that stores the arc costs.
+    ///
+    /// The type of the map that stores the arc costs.
+    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
+    typedef CM 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 conform to the \ref concepts::WriteMap
+    /// "WriteMap" concept, and its value type must be \c bool
+    /// (or convertible). Initially it will be set to \c false on each
+    /// arc, then it will be set on each arborescence arc once.
+    typedef typename Digraph::template ArcMap<bool> ArborescenceMap;
+
+    /// \brief Instantiates a \c ArborescenceMap.
+    ///
+    /// This function instantiates a \c ArborescenceMap.
+    /// \param digraph The digraph to which we would like to calculate
+    /// the \c ArborescenceMap.
+    static ArborescenceMap *createArborescenceMap(const Digraph &digraph){
+      return new ArborescenceMap(digraph);
+    }
+
+    /// \brief The type of the \c PredMap
+    ///
+    /// The type of the \c PredMap. It must confrom to the
+    /// \ref concepts::WriteMap "WriteMap" concept, and its value type
+    /// must be the \c Arc type of the digraph.
+    typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
+
+    /// \brief Instantiates a \c PredMap.
+    ///
+    /// This function instantiates a \c PredMap.
+    /// \param digraph The digraph to which we would like to define the
+    /// \c PredMap.
+    static PredMap *createPredMap(const Digraph &digraph){
+      return new PredMap(digraph);
+    }
+
+  };
+
+  /// \ingroup spantree
+  ///
+  /// \brief Minimum Cost Arborescence algorithm class.
+  ///
+  /// This class provides an efficient implementation of the
+  /// Minimum Cost Arborescence algorithm. The arborescence is a tree
+  /// which is directed from a given source node of the digraph. One or
+  /// more sources should be given to the algorithm and it will calculate
+  /// the minimum cost subgraph that is the 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 O(n<sup>2</sup>+e).
+  ///
+  /// The algorithm also provides an optimal dual solution, therefore
+  /// the optimality of the solution can be checked.
+  ///
+  /// \param GR The digraph type the algorithm runs on.
+  /// \param CM A read-only arc map storing 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 costs if
+  /// it is necessary. The default map type is \ref
+  /// concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
+  /// \param TR Traits class to set various data types used
+  /// by the algorithm. The default traits class is
+  /// \ref MinCostArborescenceDefaultTraits
+  /// "MinCostArborescenceDefaultTraits<GR, CM>".
+#ifndef DOXYGEN
+  template <typename GR,
+            typename CM = typename GR::template ArcMap<int>,
+            typename TR =
+              MinCostArborescenceDefaultTraits<GR, CM> >
+#else
+  template <typename GR, typename CM, typedef TR>
+#endif
+  class MinCostArborescence {
+  public:
+
+    /// \brief The \ref MinCostArborescenceDefaultTraits "traits class" 
+    /// of the algorithm. 
+    typedef TR 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)[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)[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)[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 SetArborescenceMapTraits : 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 \c ArborescenceMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting
+    /// \c ArborescenceMap type.
+    /// It must conform to the \ref concepts::WriteMap "WriteMap" concept,
+    /// and its value type must be \c bool (or convertible).
+    /// Initially it will be set to \c false on each arc,
+    /// then it will be set on each arborescence arc once.
+    template <class T>
+    struct SetArborescenceMap
+      : public MinCostArborescence<Digraph, CostMap,
+                                   SetArborescenceMapTraits<T> > {
+    };
+
+    template <class T>
+    struct SetPredMapTraits : public Traits {
+      typedef T PredMap;
+      static PredMap *createPredMap(const Digraph &)
+      {
+        LEMON_ASSERT(false, "PredMap is not initialized");
+        return 0; // ignore warnings
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for
+    /// setting \c PredMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting
+    /// \c PredMap type.
+    /// It must meet the \ref concepts::WriteMap "WriteMap" concept, 
+    /// and its value type must be the \c Arc type of the digraph.
+    template <class T>
+    struct SetPredMap
+      : public MinCostArborescence<Digraph, CostMap, SetPredMapTraits<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 <tt>(*this)</tt>
+    MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
+      if (local_arborescence) {
+        delete _arborescence;
+      }
+      local_arborescence = false;
+      _arborescence = &m;
+      return *this;
+    }
+
+    /// \brief Sets the predecessor map.
+    ///
+    /// Sets the predecessor map.
+    /// \return <tt>(*this)</tt>
+    MinCostArborescence& predMap(PredMap& m) {
+      if (local_pred) {
+        delete _pred;
+      }
+      local_pred = false;
+      _pred = &m;
+      return *this;
+    }
+
+    /// \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 better control on the execution,
+    /// you have to call \ref init() first, 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)[it] = -3;
+        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
+        _pred->set(it, INVALID);
+      }
+      for (ArcIt it(*_digraph); it != INVALID; ++it) {
+        _arborescence->set(it, false);
+        (*_arc_order)[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 in the priority
+    /// queue.
+    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 s) {
+      init();
+      addSource(s);
+      start();
+    }
+
+    ///@}
+
+    /// \name Query Functions
+    /// The result of the %MinCostArborescence algorithm can be obtained
+    /// using these functions.\n
+    /// Either run() or start() must be called before using them.
+
+    /// @{
+
+    /// \brief Returns the cost of the arborescence.
+    ///
+    /// Returns the cost of the arborescence.
+    Value arborescenceCost() const {
+      Value sum = 0;
+      for (ArcIt it(*_digraph); it != INVALID; ++it) {
+        if (arborescence(it)) {
+          sum += (*_cost)[it];
+        }
+      }
+      return sum;
+    }
+
+    /// \brief Returns \c true if the arc is in the arborescence.
+    ///
+    /// Returns \c true if the given arc is in the arborescence.
+    /// \param arc An 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 const reference to the arborescence map.
+    ///
+    /// Returns a const reference to the arborescence map.
+    /// \pre \ref run() must be called before using this function.
+    const ArborescenceMap& arborescenceMap() const {
+      return *_arborescence;
+    }
+
+    /// \brief Returns the predecessor arc of the given node.
+    ///
+    /// Returns the predecessor arc of the given node.
+    /// \pre \ref run() must be called before using this function.
+    Arc pred(Node node) const {
+      return (*_pred)[node];
+    }
+
+    /// \brief Returns a const reference to the pred map.
+    ///
+    /// Returns a const reference to the pred map.
+    /// \pre \ref run() must be called before using this function.
+    const PredMap& predMap() const {
+      return *_pred;
+    }
+
+    /// \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.
+    Value dualValue(int k) const {
+      return _dual_variables[k].value;
+    }
+
+    /// \brief LEMON iterator for getting a dual variable.
+    ///
+    /// This class provides a common style LEMON iterator for getting a
+    /// dual variable of \ref MinCostArborescence algorithm.
+    /// It iterates over a subset of the nodes.
+    class DualIt {
+    public:
+
+      /// \brief Constructor.
+      ///
+      /// Constructor for getting the nodeset of the dual variable
+      /// of \ref MinCostArborescence algorithm.
+      DualIt(const MinCostArborescence& algorithm, int variable)
+        : _algorithm(&algorithm)
+      {
+        _index = _algorithm->_dual_variables[variable].begin;
+        _last = _algorithm->_dual_variables[variable].end;
+      }
+
+      /// \brief Conversion to \c Node.
+      ///
+      /// Conversion to \c 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;
+    };
+
+    /// @}
+
+  };
+
+  /// \ingroup spantree
+  ///
+  /// \brief Function type interface for MinCostArborescence algorithm.
+  ///
+  /// Function type interface for MinCostArborescence algorithm.
+  /// \param digraph The digraph the algorithm runs on.
+  /// \param cost An arc map storing the costs.
+  /// \param source The source node of the arborescence.
+  /// \retval arborescence An arc map with \c bool (or convertible) value
+  /// type that stores the arborescence.
+  /// \return The total 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 SetArborescenceMap<ArborescenceMap>
+      ::Create mca(digraph, cost);
+    mca.arborescenceMap(arborescence);
+    mca.run(source);
+    return mca.arborescenceCost();
+  }
+
+}
+
+#endif