Index: lemon/Makefile.am
===================================================================
--- lemon/Makefile.am	(revision 1999)
+++ lemon/Makefile.am	(revision 2017)
@@ -60,4 +60,5 @@
 	map_iterator.h \
 	max_matching.h \
+	min_cost_arborescence.h \
 	min_cost_flow.h \
 	min_cut.h \
Index: lemon/min_cost_arborescence.h
===================================================================
--- lemon/min_cost_arborescence.h	(revision 2017)
+++ lemon/min_cost_arborescence.h	(revision 2017)
@@ -0,0 +1,561 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2006
+ * 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>
+
+namespace lemon {
+
+
+  /// \brief Default traits class of MinCostArborescence class.
+  /// 
+  /// Default traits class of MinCostArborescence class.
+  /// \param _Graph Graph type.
+  /// \param _CostMap Type of cost map.
+  template <class _Graph, class _CostMap>
+  struct MinCostArborescenceDefaultTraits{
+
+    /// \brief The graph type the algorithm runs on. 
+    typedef _Graph Graph;
+
+    /// \brief The type of the map that stores the edge costs.
+    ///
+    /// The type of the map that stores the edge costs.
+    /// It must meet the \ref concept::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 edges are 
+    /// in the arborescence.
+    ///
+    /// The type of the map that stores which edges are in the arborescence.
+    /// It must meet the \ref concept::ReadWriteMap "ReadWriteMap" concept.
+    /// Initially it will be setted to false on each edge. The algorithm
+    /// may set each value one time to true and maybe after it to false again.
+    /// Therefore you cannot use maps like BackInserteBoolMap with this
+    /// algorithm.   
+    typedef typename Graph::template EdgeMap<bool> ArborescenceMap; 
+
+    /// \brief Instantiates a ArborescenceMap.
+    ///
+    /// This function instantiates a \ref ArborescenceMap. 
+    /// \param _graph is the graph, to which we would like to define the 
+    /// ArborescenceMap.
+    static ArborescenceMap *createArborescenceMap(const Graph &_graph){
+      return new ArborescenceMap(_graph);
+    }
+
+  };
+
+  /// \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 graph. 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 O(n^2 + e).
+  ///
+  /// \param _Graph The graph type the algorithm runs on. The default value
+  /// is \ref ListGraph. The value of _Graph is not used directly by
+  /// MinCostArborescence, it is only passed to 
+  /// \ref MinCostArborescenceDefaultTraits.
+  /// \param _CostMap This read-only EdgeMap determines the costs of the
+  /// edges. It is read once for each edge, so the map may involve in
+  /// relatively time consuming process to compute the edge cost if
+  /// it is necessary. The default map type is \ref
+  /// concept::StaticGraph::EdgeMap "Graph::EdgeMap<int>".  The value
+  /// of _CostMap is not used directly by MinCostArborescence, 
+  /// it is only passed to \ref MinCostArborescenceDefaultTraits.  
+  /// \param _Traits Traits class to set various data types used 
+  /// by the algorithm.  The default traits class is 
+  /// \ref MinCostArborescenceDefaultTraits
+  /// "MinCostArborescenceDefaultTraits<_Graph,_CostMap>".  See \ref
+  /// MinCostArborescenceDefaultTraits for the documentation of a 
+  /// MinCostArborescence traits class.
+  ///
+  /// \author Balazs Dezso
+#ifndef DOXYGEN
+  template <typename _Graph = ListGraph, 
+            typename _CostMap = typename _Graph::template EdgeMap<int>,
+            typename _Traits = 
+            MinCostArborescenceDefaultTraits<_Graph, _CostMap> >
+#else 
+  template <typename _Graph, typename _CostMap, typedef _Traits>
+#endif
+  class MinCostArborescence {
+  public:
+
+    /// \brief \ref Exception for uninitialized parameters.
+    ///
+    /// This error represents problems in the initialization
+    /// of the parameters of the algorithms.    
+    class UninitializedParameter : public lemon::UninitializedParameter {
+    public:
+      virtual const char* exceptionName() const {
+	return "lemon::MinCostArborescence::UninitializedParameter";
+      }
+    };
+
+    /// The traits.
+    typedef _Traits Traits;
+    /// The type of the underlying graph.
+    typedef typename Traits::Graph Graph;
+    /// The type of the map that stores the edge costs.
+    typedef typename Traits::CostMap CostMap;
+    ///The type of the costs of the edges.
+    typedef typename Traits::Value Value;
+    ///The type of the map that stores which edges are in the arborescence.
+    typedef typename Traits::ArborescenceMap ArborescenceMap;
+
+  protected:
+
+    typedef typename Graph::Node Node;
+    typedef typename Graph::Edge Edge;
+    typedef typename Graph::NodeIt NodeIt;
+    typedef typename Graph::EdgeIt EdgeIt;
+    typedef typename Graph::InEdgeIt InEdgeIt;
+    typedef typename Graph::OutEdgeIt OutEdgeIt;
+
+    struct CostEdge {
+
+      Edge edge;
+      Value value;
+
+      CostEdge() {}
+      CostEdge(Edge _edge, Value _value) : edge(_edge), value(_value) {}
+
+    };
+
+    const Graph* graph;
+    const CostMap* cost;
+
+    ArborescenceMap* _arborescence_map;
+    bool local_arborescence_map;
+
+    typedef typename Graph::template NodeMap<int> LevelMap;
+    LevelMap *_level;
+
+    typedef typename Graph::template NodeMap<CostEdge> CostEdgeMap;
+    CostEdgeMap *_cost_edges; 
+
+    struct StackLevel {
+
+      std::vector<CostEdge> edges;
+      int node_level;
+
+    };
+
+    std::vector<StackLevel> level_stack;    
+    std::vector<Node> queue;
+
+    int node_counter;
+
+  public:
+
+    /// \name Named template parameters
+
+    /// @{
+
+    template <class T>
+    struct DefArborescenceMapTraits : public Traits {
+      typedef T ArborescenceMap;
+      static ArborescenceMap *createArborescenceMap(const Graph &)
+      {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// \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<Graph, CostMap,
+                                   DefArborescenceMapTraits<T> > {
+      typedef MinCostArborescence<Graph, CostMap, 
+                                   DefArborescenceMapTraits<T> > Create;
+    };
+    
+    /// @}
+
+    /// \brief Constructor.
+    ///
+    /// \param _graph The graph the algorithm will run on.
+    /// \param _cost The cost map used by the algorithm.
+    MinCostArborescence(const Graph& _graph, const CostMap& _cost) 
+      : graph(&_graph), cost(&_cost),
+        _arborescence_map(0), local_arborescence_map(false), 
+        _level(0), _cost_edges(0) {}
+
+    /// \brief Destructor.
+    ~MinCostArborescence() {
+      destroyStructures();
+    }
+
+    /// \brief Sets the arborescence map.
+    /// 
+    /// Sets the arborescence map.
+    /// \return \c (*this)
+    MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
+      _arborescence_map = &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_map;
+    }
+
+    /// \brief Returns true if the edge is in the arborescence.
+    ///
+    /// Returns true if the edge is in the arborescence.
+    /// \param edge The edge of the graph.
+    /// \pre \ref run() must be called before using this function.
+    bool arborescenceEdge(Edge edge) const {
+      return (*_arborescence_map)[edge];
+    }
+ 
+    /// \brief Returns the cost of the arborescence.
+    ///
+    /// Returns the cost of the arborescence.
+   Value arborescenceCost() const {
+      Value sum = 0;
+      for (EdgeIt it(*graph); it != INVALID; ++it) {
+        if (arborescenceEdge(it)) {
+          sum += (*cost)[it];
+        }
+      }
+      return sum;
+    }
+
+    /// @}
+    
+    /// \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 actual path
+    /// computation.
+
+    ///@{
+
+    /// \brief Initializes the internal data structures.
+    ///
+    /// Initializes the internal data structures.
+    ///
+    void init() {
+      initStructures();
+      for (NodeIt it(*graph); it != INVALID; ++it) {
+        (*_cost_edges)[it].edge = INVALID;
+        (*_level)[it] = -3; 
+      }
+      for (EdgeIt it(*graph); it != INVALID; ++it) {
+        _arborescence_map->set(it, false);
+      }
+    }
+
+    /// \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 (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
+          if ((*_level)[graph->target(it)] == -3) {
+            (*_level)[graph->target(it)] = -2;
+            nodes.push_back(graph->target(it));
+            queue.push_back(graph->target(it));
+          }
+        }
+      }
+      (*_level)[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_counter = 0;
+      Node node = queue.back();
+      queue.pop_back();
+      if ((*_level)[node] == -2) {
+        Edge edge = prepare(node);
+        while ((*_level)[graph->source(edge)] != -1) {
+          if ((*_level)[graph->source(edge)] >= 0) {
+            edge = contract(bottom((*_level)[graph->source(edge)]));
+          } else {
+            edge = prepare(graph->source(edge));
+          }
+        }
+        finalize(graph->target(edge));
+        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();
+    }
+
+    ///@}
+
+  protected:
+
+    void initStructures() {
+      if (!_arborescence_map) {
+        local_arborescence_map = true;
+        _arborescence_map = Traits::createArborescenceMap(*graph);
+      }
+      if (!_level) {
+        _level = new LevelMap(*graph);
+      }
+      if (!_cost_edges) {
+        _cost_edges = new CostEdgeMap(*graph);
+      }
+    }
+
+    void destroyStructures() {
+      if (_level) {
+        delete _level;
+      }
+      if (!_cost_edges) {
+        delete _cost_edges;
+      }
+      if (local_arborescence_map) {
+        delete _arborescence_map;
+      }
+    }
+
+    Edge prepare(Node node) {
+      std::vector<Node> nodes;
+      (*_level)[node] = node_counter;
+      for (InEdgeIt it(*graph, node); it != INVALID; ++it) {
+        Edge edge = it;
+        Value value = (*cost)[it];
+        if (graph->source(edge) == node || 
+            (*_level)[graph->source(edge)] == -3) continue;
+        if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {
+          (*_cost_edges)[graph->source(edge)].edge = edge;
+          (*_cost_edges)[graph->source(edge)].value = value;
+          nodes.push_back(graph->source(edge));
+        } else {
+          if ((*_cost_edges)[graph->source(edge)].value > value) {
+            (*_cost_edges)[graph->source(edge)].edge = edge;
+            (*_cost_edges)[graph->source(edge)].value = value;
+          }
+        }      
+      }
+      CostEdge minimum = (*_cost_edges)[nodes[0]]; 
+      for (int i = 1; i < (int)nodes.size(); ++i) {
+        if ((*_cost_edges)[nodes[i]].value < minimum.value) {
+          minimum = (*_cost_edges)[nodes[i]];
+        }
+      }
+      StackLevel level;
+      level.node_level = node_counter;
+      for (int i = 0; i < (int)nodes.size(); ++i) {
+        (*_cost_edges)[nodes[i]].value -= minimum.value;
+        level.edges.push_back((*_cost_edges)[nodes[i]]);
+        (*_cost_edges)[nodes[i]].edge = INVALID;
+      }
+      level_stack.push_back(level);
+      ++node_counter;
+      _arborescence_map->set(minimum.edge, true);
+      return minimum.edge;
+    }
+  
+    Edge contract(int node_bottom) {
+      std::vector<Node> nodes;
+      while (!level_stack.empty() && 
+             level_stack.back().node_level >= node_bottom) {
+        for (int i = 0; i < (int)level_stack.back().edges.size(); ++i) {
+          Edge edge = level_stack.back().edges[i].edge;
+          Value value = level_stack.back().edges[i].value;
+          if ((*_level)[graph->source(edge)] >= node_bottom) continue;
+          if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {
+            (*_cost_edges)[graph->source(edge)].edge = edge;
+            (*_cost_edges)[graph->source(edge)].value = value;
+            nodes.push_back(graph->source(edge));
+          } else {
+            if ((*_cost_edges)[graph->source(edge)].value > value) {
+              (*_cost_edges)[graph->source(edge)].edge = edge;
+              (*_cost_edges)[graph->source(edge)].value = value;
+            }
+          }
+        }
+        level_stack.pop_back();
+      }
+      CostEdge minimum = (*_cost_edges)[nodes[0]]; 
+      for (int i = 1; i < (int)nodes.size(); ++i) {
+        if ((*_cost_edges)[nodes[i]].value < minimum.value) {
+          minimum = (*_cost_edges)[nodes[i]];
+        }
+      }
+      StackLevel level;
+      level.node_level = node_bottom;
+      for (int i = 0; i < (int)nodes.size(); ++i) {
+        (*_cost_edges)[nodes[i]].value -= minimum.value;
+        level.edges.push_back((*_cost_edges)[nodes[i]]);
+        (*_cost_edges)[nodes[i]].edge = INVALID;
+      }
+      level_stack.push_back(level);
+      _arborescence_map->set(minimum.edge, true);
+      return minimum.edge;
+    }
+
+    int bottom(int level) {
+      int k = level_stack.size() - 1;
+      while (level_stack[k].node_level > level) {
+        --k;
+      }
+      return level_stack[k].node_level;
+    }
+
+    void finalize(Node source) {
+      std::vector<Node> nodes;
+      nodes.push_back(source);
+      while (!nodes.empty()) {
+        Node node = nodes.back();
+        nodes.pop_back();
+        for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
+          if ((*_level)[graph->target(it)] >= 0 && (*_arborescence_map)[it]) {
+            (*_level)[graph->target(it)] = -1;
+            nodes.push_back(graph->target(it));
+          } else {
+            _arborescence_map->set(it, false);
+          }
+        }
+      }
+      (*_level)[source] = -1;      
+    }
+
+  };
+
+  /// \ingroup spantree
+  ///
+  /// \brief Function type interface for MinCostArborescence algorithm.
+  ///
+  /// Function type interface for MinCostArborescence algorithm.
+  /// \param graph The Graph that the algorithm runs on.
+  /// \param cost The CostMap of the edges.
+  /// \param source The source of the arborescence.
+  /// \retval arborescence The bool EdgeMap which stores the arborescence.
+  /// \return The cost of the arborescence. 
+  ///
+  /// \sa MinCostArborescence
+  template <typename Graph, typename CostMap, typename ArborescenceMap>
+  typename CostMap::Value minCostArborescence(const Graph& graph, 
+                                              const CostMap& cost,
+                                              typename Graph::Node source,
+                                              ArborescenceMap& arborescence) {
+    typename MinCostArborescence<Graph, CostMap>
+      ::template DefArborescenceMap<ArborescenceMap>
+      ::Create mca(graph, cost);
+    mca.arborescenceMap(arborescence);
+    mca.run(source);
+    return mca.arborescenceCost();
+  }
+
+}
+
+#endif
+
+// Hilbert - Huang
