lemon/min_cost_arborescence.h
author alpar
Mon, 30 Oct 2006 16:19:18 +0000
changeset 2274 432d0469a87e
parent 2260 4274224f8a7d
child 2385 096d83158d41
permissions -rw-r--r--
Turn off 32 bit only tests, cont'd.
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2006
     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 
    31 namespace lemon {
    32 
    33 
    34   /// \brief Default traits class of MinCostArborescence class.
    35   /// 
    36   /// Default traits class of MinCostArborescence class.
    37   /// \param _Graph Graph type.
    38   /// \param _CostMap Type of cost map.
    39   template <class _Graph, class _CostMap>
    40   struct MinCostArborescenceDefaultTraits{
    41 
    42     /// \brief The graph type the algorithm runs on. 
    43     typedef _Graph Graph;
    44 
    45     /// \brief The type of the map that stores the edge costs.
    46     ///
    47     /// The type of the map that stores the edge costs.
    48     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
    49     typedef _CostMap CostMap;
    50 
    51     /// \brief The value type of the costs.
    52     ///
    53     /// The value type of the costs.
    54     typedef typename CostMap::Value Value;
    55 
    56     /// \brief The type of the map that stores which edges are 
    57     /// in the arborescence.
    58     ///
    59     /// The type of the map that stores which edges are in the arborescence.
    60     /// It must meet the \ref concepts::WriteMap "WriteMap" concept.
    61     /// Initially it will be set to false on each edge. After it
    62     /// will set all arborescence edges once.
    63     typedef typename Graph::template EdgeMap<bool> ArborescenceMap; 
    64 
    65     /// \brief Instantiates a ArborescenceMap.
    66     ///
    67     /// This function instantiates a \ref ArborescenceMap. 
    68     /// \param _graph is the graph, to which we would like to define the 
    69     /// ArborescenceMap.
    70     static ArborescenceMap *createArborescenceMap(const Graph &_graph){
    71       return new ArborescenceMap(_graph);
    72     }
    73 
    74     /// \brief The type of the PredMap
    75     ///
    76     /// The type of the PredMap. It is a node map with an edge value type.
    77     typedef typename Graph::template NodeMap<typename Graph::Edge> PredMap;
    78 
    79     /// \brief Instantiates a PredMap.
    80     ///
    81     /// This function instantiates a \ref PredMap. 
    82     /// \param _graph is the graph, to which we would like to define the 
    83     /// PredMap.
    84     static PredMap *createPredMap(const Graph &_graph){
    85       return new PredMap(_graph);
    86     }
    87     
    88   };
    89 
    90   /// \ingroup spantree
    91   ///
    92   /// \brief %MinCostArborescence algorithm class.
    93   ///
    94   /// This class provides an efficient implementation of 
    95   /// %MinCostArborescence algorithm. The arborescence is a tree 
    96   /// which is directed from a given source node of the graph. One or
    97   /// more sources should be given for the algorithm and it will calculate
    98   /// the minimum cost subgraph which are union of arborescences with the
    99   /// given sources and spans all the nodes which are reachable from the
   100   /// sources. The time complexity of the algorithm is \f$ O(n^2+e) \f$.
   101   ///
   102   /// The algorithm provides also an optimal dual solution to arborescence
   103   /// that way the optimality of the solution can be proofed easily.
   104   ///
   105   /// \param _Graph The graph type the algorithm runs on. The default value
   106   /// is \ref ListGraph. The value of _Graph is not used directly by
   107   /// MinCostArborescence, it is only passed to 
   108   /// \ref MinCostArborescenceDefaultTraits.
   109   /// \param _CostMap This read-only EdgeMap determines the costs of the
   110   /// edges. It is read once for each edge, so the map may involve in
   111   /// relatively time consuming process to compute the edge cost if
   112   /// it is necessary. The default map type is \ref
   113   /// concepts::Graph::EdgeMap "Graph::EdgeMap<int>".  The value
   114   /// of _CostMap is not used directly by MinCostArborescence, 
   115   /// it is only passed to \ref MinCostArborescenceDefaultTraits.  
   116   /// \param _Traits Traits class to set various data types used 
   117   /// by the algorithm.  The default traits class is 
   118   /// \ref MinCostArborescenceDefaultTraits
   119   /// "MinCostArborescenceDefaultTraits<_Graph,_CostMap>".  See \ref
   120   /// MinCostArborescenceDefaultTraits for the documentation of a 
   121   /// MinCostArborescence traits class.
   122   ///
   123   /// \author Balazs Dezso
   124 #ifndef DOXYGEN
   125   template <typename _Graph = ListGraph, 
   126             typename _CostMap = typename _Graph::template EdgeMap<int>,
   127             typename _Traits = 
   128             MinCostArborescenceDefaultTraits<_Graph, _CostMap> >
   129 #else 
   130   template <typename _Graph, typename _CostMap, typedef _Traits>
   131 #endif
   132   class MinCostArborescence {
   133   public:
   134 
   135     /// \brief \ref Exception for uninitialized parameters.
   136     ///
   137     /// This error represents problems in the initialization
   138     /// of the parameters of the algorithms.    
   139     class UninitializedParameter : public lemon::UninitializedParameter {
   140     public:
   141       virtual const char* what() const throw() {
   142 	return "lemon::MinCostArborescence::UninitializedParameter";
   143       }
   144     };
   145 
   146     /// The traits.
   147     typedef _Traits Traits;
   148     /// The type of the underlying graph.
   149     typedef typename Traits::Graph Graph;
   150     /// The type of the map that stores the edge costs.
   151     typedef typename Traits::CostMap CostMap;
   152     ///The type of the costs of the edges.
   153     typedef typename Traits::Value Value;
   154     ///The type of the predecessor map.
   155     typedef typename Traits::PredMap PredMap;
   156     ///The type of the map that stores which edges are in the arborescence.
   157     typedef typename Traits::ArborescenceMap ArborescenceMap;
   158 
   159   protected:
   160 
   161     typedef typename Graph::Node Node;
   162     typedef typename Graph::Edge Edge;
   163     typedef typename Graph::NodeIt NodeIt;
   164     typedef typename Graph::EdgeIt EdgeIt;
   165     typedef typename Graph::InEdgeIt InEdgeIt;
   166     typedef typename Graph::OutEdgeIt OutEdgeIt;
   167 
   168     struct CostEdge {
   169 
   170       Edge edge;
   171       Value value;
   172 
   173       CostEdge() {}
   174       CostEdge(Edge _edge, Value _value) : edge(_edge), value(_value) {}
   175 
   176     };
   177 
   178     const Graph *graph;
   179     const CostMap *cost;
   180 
   181     PredMap *_pred;
   182     bool local_pred;
   183 
   184     ArborescenceMap *_arborescence;
   185     bool local_arborescence;
   186 
   187     typedef typename Graph::template EdgeMap<int> EdgeOrder;
   188     EdgeOrder *_edge_order;
   189     
   190     typedef typename Graph::template NodeMap<int> NodeOrder;
   191     NodeOrder *_node_order;
   192 
   193     typedef typename Graph::template NodeMap<CostEdge> CostEdgeMap;
   194     CostEdgeMap *_cost_edges; 
   195 
   196     struct StackLevel {
   197 
   198       std::vector<CostEdge> edges;
   199       int node_level;
   200 
   201     };
   202 
   203     std::vector<StackLevel> level_stack;    
   204     std::vector<Node> queue;
   205 
   206     typedef std::vector<typename Graph::Node> DualNodeList;
   207 
   208     DualNodeList _dual_node_list;
   209 
   210     struct DualVariable {
   211       int begin, end;
   212       Value value;
   213       
   214       DualVariable(int _begin, int _end, Value _value)
   215         : begin(_begin), end(_end), value(_value) {}
   216 
   217     };
   218 
   219     typedef std::vector<DualVariable> DualVariables;
   220 
   221     DualVariables _dual_variables;
   222 
   223     typedef typename Graph::template NodeMap<int> HeapCrossRef;
   224     
   225     HeapCrossRef *_heap_cross_ref;
   226 
   227     typedef BinHeap<int, HeapCrossRef> Heap;
   228 
   229     Heap *_heap;
   230 
   231   public:
   232 
   233     /// \name Named template parameters
   234 
   235     /// @{
   236 
   237     template <class T>
   238     struct DefArborescenceMapTraits : public Traits {
   239       typedef T ArborescenceMap;
   240       static ArborescenceMap *createArborescenceMap(const Graph &)
   241       {
   242 	throw UninitializedParameter();
   243       }
   244     };
   245 
   246     /// \brief \ref named-templ-param "Named parameter" for 
   247     /// setting ArborescenceMap type
   248     ///
   249     /// \ref named-templ-param "Named parameter" for setting 
   250     /// ArborescenceMap type
   251     template <class T>
   252     struct DefArborescenceMap 
   253       : public MinCostArborescence<Graph, CostMap,
   254                                    DefArborescenceMapTraits<T> > {
   255       typedef MinCostArborescence<Graph, CostMap, 
   256                                    DefArborescenceMapTraits<T> > Create;
   257     };
   258 
   259     template <class T>
   260     struct DefPredMapTraits : public Traits {
   261       typedef T PredMap;
   262       static PredMap *createPredMap(const Graph &)
   263       {
   264 	throw UninitializedParameter();
   265       }
   266     };
   267 
   268     /// \brief \ref named-templ-param "Named parameter" for 
   269     /// setting PredMap type
   270     ///
   271     /// \ref named-templ-param "Named parameter" for setting 
   272     /// PredMap type
   273     template <class T>
   274     struct DefPredMap 
   275       : public MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > {
   276       typedef MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > Create;
   277     };
   278     
   279     /// @}
   280 
   281     /// \brief Constructor.
   282     ///
   283     /// \param _graph The graph the algorithm will run on.
   284     /// \param _cost The cost map used by the algorithm.
   285     MinCostArborescence(const Graph& _graph, const CostMap& _cost) 
   286       : graph(&_graph), cost(&_cost), _pred(0), local_pred(false),
   287         _arborescence(0), local_arborescence(false), 
   288         _edge_order(0), _node_order(0), _cost_edges(0), 
   289         _heap_cross_ref(0), _heap(0) {}
   290 
   291     /// \brief Destructor.
   292     ~MinCostArborescence() {
   293       destroyStructures();
   294     }
   295 
   296     /// \brief Sets the arborescence map.
   297     /// 
   298     /// Sets the arborescence map.
   299     /// \return \c (*this)
   300     MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
   301       if (local_arborescence) {
   302         delete _arborescence;
   303       }
   304       local_arborescence = false;
   305       _arborescence = &m;
   306       return *this;
   307     }
   308 
   309     /// \brief Sets the arborescence map.
   310     /// 
   311     /// Sets the arborescence map.
   312     /// \return \c (*this)
   313     MinCostArborescence& predMap(PredMap& m) {
   314       if (local_pred) {
   315         delete _pred;
   316       }
   317       local_pred = false;
   318       _pred = &m;
   319       return *this;
   320     }
   321 
   322     /// \name Query Functions
   323     /// The result of the %MinCostArborescence algorithm can be obtained 
   324     /// using these functions.\n
   325     /// Before the use of these functions,
   326     /// either run() or start() must be called.
   327     
   328     /// @{
   329 
   330     /// \brief Returns a reference to the arborescence map.
   331     ///
   332     /// Returns a reference to the arborescence map.
   333     const ArborescenceMap& arborescenceMap() const {
   334       return *_arborescence;
   335     }
   336 
   337     /// \brief Returns true if the edge is in the arborescence.
   338     ///
   339     /// Returns true if the edge is in the arborescence.
   340     /// \param edge The edge of the graph.
   341     /// \pre \ref run() must be called before using this function.
   342     bool arborescence(Edge edge) const {
   343       return (*_pred)[graph->target(edge)] == edge;
   344     }
   345 
   346     /// \brief Returns a reference to the pred map.
   347     ///
   348     /// Returns a reference to the pred map.
   349     const PredMap& predMap() const {
   350       return *_pred;
   351     }
   352 
   353     /// \brief Returns the predecessor edge of the given node.
   354     ///
   355     /// Returns the predecessor edge of the given node.
   356     bool pred(Node node) const {
   357       return (*_pred)[node];
   358     }
   359  
   360     /// \brief Returns the cost of the arborescence.
   361     ///
   362     /// Returns the cost of the arborescence.
   363     Value arborescenceValue() const {
   364       Value sum = 0;
   365       for (EdgeIt it(*graph); it != INVALID; ++it) {
   366         if (arborescence(it)) {
   367           sum += (*cost)[it];
   368         }
   369       }
   370       return sum;
   371     }
   372 
   373     /// \brief Indicates that a node is reachable from the sources.
   374     ///
   375     /// Indicates that a node is reachable from the sources.
   376     bool reached(Node node) const {
   377       return (*_node_order)[node] != -3;
   378     }
   379 
   380     /// \brief Indicates that a node is processed.
   381     ///
   382     /// Indicates that a node is processed. The arborescence path exists 
   383     /// from the source to the given node.
   384     bool processed(Node node) const {
   385       return (*_node_order)[node] == -1;
   386     }
   387 
   388     /// \brief Returns the number of the dual variables in basis.
   389     ///
   390     /// Returns the number of the dual variables in basis.
   391     int dualSize() const {
   392       return _dual_variables.size();
   393     }
   394 
   395     /// \brief Returns the value of the dual solution.
   396     ///
   397     /// Returns the value of the dual solution. It should be
   398     /// equal to the arborescence value.
   399     Value dualValue() const {
   400       Value sum = 0;
   401       for (int i = 0; i < (int)_dual_variables.size(); ++i) {
   402         sum += _dual_variables[i].value;
   403       }
   404       return sum;
   405     }
   406 
   407     /// \brief Returns the number of the nodes in the dual variable.
   408     ///
   409     /// Returns the number of the nodes in the dual variable.
   410     int dualSize(int k) const {
   411       return _dual_variables[k].end - _dual_variables[k].begin;
   412     }
   413 
   414     /// \brief Returns the value of the dual variable.
   415     ///
   416     /// Returns the the value of the dual variable.
   417     const Value& dualValue(int k) const {
   418       return _dual_variables[k].value;
   419     }
   420 
   421     /// \brief Lemon iterator for get a dual variable.
   422     ///
   423     /// Lemon iterator for get a dual variable. This class provides
   424     /// a common style lemon iterator which gives back a subset of
   425     /// the nodes.
   426     class DualIt {
   427     public:
   428 
   429       /// \brief Constructor.
   430       ///
   431       /// Constructor for get the nodeset of the variable. 
   432       DualIt(const MinCostArborescence& algorithm, int variable) 
   433         : _algorithm(&algorithm), _variable(variable) 
   434       {
   435         _index = _algorithm->_dual_variables[_variable].begin;
   436       }
   437 
   438       /// \brief Invalid constructor.
   439       ///
   440       /// Invalid constructor.
   441       DualIt(Invalid) : _algorithm(0) {}
   442 
   443       /// \brief Conversion to node.
   444       ///
   445       /// Conversion to node.
   446       operator Node() const { 
   447         return _algorithm ? _algorithm->_dual_node_list[_index] : INVALID;
   448       }
   449 
   450       /// \brief Increment operator.
   451       ///
   452       /// Increment operator.
   453       DualIt& operator++() {
   454         ++_index;
   455         if (_algorithm->_dual_variables[_variable].end == _index) {
   456           _algorithm = 0;
   457         }
   458         return *this; 
   459       }
   460 
   461       bool operator==(const DualIt& it) const { 
   462         return (Node)(*this) == (Node)it; 
   463       }
   464       bool operator!=(const DualIt& it) const { 
   465         return (Node)(*this) != (Node)it; 
   466       }
   467       bool operator<(const DualIt& it) const { 
   468         return (Node)(*this) < (Node)it; 
   469       }
   470       
   471     private:
   472       const MinCostArborescence* _algorithm;
   473       int _variable;
   474       int _index;
   475     };
   476 
   477     /// @}
   478     
   479     /// \name Execution control
   480     /// The simplest way to execute the algorithm is to use
   481     /// one of the member functions called \c run(...). \n
   482     /// If you need more control on the execution,
   483     /// first you must call \ref init(), then you can add several 
   484     /// source nodes with \ref addSource().
   485     /// Finally \ref start() will perform the arborescence
   486     /// computation.
   487 
   488     ///@{
   489 
   490     /// \brief Initializes the internal data structures.
   491     ///
   492     /// Initializes the internal data structures.
   493     ///
   494     void init() {
   495       initStructures();
   496       _heap->clear();
   497       for (NodeIt it(*graph); it != INVALID; ++it) {
   498         (*_cost_edges)[it].edge = INVALID;
   499         _node_order->set(it, -3); 
   500         _heap_cross_ref->set(it, Heap::PRE_HEAP);
   501       }
   502       for (EdgeIt it(*graph); it != INVALID; ++it) {
   503         _arborescence->set(it, false);
   504         _edge_order->set(it, -1);
   505       }
   506       _dual_node_list.clear();
   507       _dual_variables.clear();
   508     }
   509 
   510     /// \brief Adds a new source node.
   511     ///
   512     /// Adds a new source node to the algorithm.
   513     void addSource(Node source) {
   514       std::vector<Node> nodes;
   515       nodes.push_back(source);
   516       while (!nodes.empty()) {
   517         Node node = nodes.back();
   518         nodes.pop_back();
   519         for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
   520           Node target = graph->target(it);
   521           if ((*_node_order)[target] == -3) {
   522             (*_node_order)[target] = -2;
   523             nodes.push_back(target);
   524             queue.push_back(target);
   525           }
   526         }
   527       }
   528       (*_node_order)[source] = -1;
   529     }
   530 
   531     /// \brief Processes the next node in the priority queue.
   532     ///
   533     /// Processes the next node in the priority queue.
   534     ///
   535     /// \return The processed node.
   536     ///
   537     /// \warning The queue must not be empty!
   538     Node processNextNode() {
   539       Node node = queue.back();
   540       queue.pop_back();
   541       if ((*_node_order)[node] == -2) {
   542         Edge edge = prepare(node);
   543         Node source = graph->source(edge);
   544         while ((*_node_order)[source] != -1) {
   545           if ((*_node_order)[source] >= 0) {
   546             edge = contract(source);
   547           } else {
   548             edge = prepare(source);
   549           }
   550           source = graph->source(edge);
   551         }
   552         finalize(edge);
   553         level_stack.clear();        
   554       }
   555       return node;
   556     }
   557 
   558     /// \brief Returns the number of the nodes to be processed.
   559     ///
   560     /// Returns the number of the nodes to be processed.
   561     int queueSize() const {
   562       return queue.size();
   563     }
   564 
   565     /// \brief Returns \c false if there are nodes to be processed.
   566     ///
   567     /// Returns \c false if there are nodes to be processed.
   568     bool emptyQueue() const {
   569       return queue.empty();
   570     }
   571 
   572     /// \brief Executes the algorithm.
   573     ///
   574     /// Executes the algorithm.
   575     ///
   576     /// \pre init() must be called and at least one node should be added
   577     /// with addSource() before using this function.
   578     ///
   579     ///\note mca.start() is just a shortcut of the following code.
   580     ///\code
   581     ///while (!mca.emptyQueue()) {
   582     ///  mca.processNextNode();
   583     ///}
   584     ///\endcode
   585     void start() {
   586       while (!emptyQueue()) {
   587         processNextNode();
   588       }
   589     }
   590 
   591     /// \brief Runs %MinCostArborescence algorithm from node \c s.
   592     /// 
   593     /// This method runs the %MinCostArborescence algorithm from 
   594     /// a root node \c s.
   595     ///
   596     ///\note mca.run(s) is just a shortcut of the following code.
   597     ///\code
   598     ///mca.init();
   599     ///mca.addSource(s);
   600     ///mca.start();
   601     ///\endcode
   602     void run(Node node) {
   603       init();
   604       addSource(node);
   605       start();
   606     }
   607 
   608     ///@}
   609 
   610   protected:
   611 
   612     void initStructures() {
   613       if (!_pred) {
   614         local_pred = true;
   615         _pred = Traits::createPredMap(*graph);
   616       }
   617       if (!_arborescence) {
   618         local_arborescence = true;
   619         _arborescence = Traits::createArborescenceMap(*graph);
   620       }
   621       if (!_edge_order) {
   622         _edge_order = new EdgeOrder(*graph);
   623       }
   624       if (!_node_order) {
   625         _node_order = new NodeOrder(*graph);
   626       }
   627       if (!_cost_edges) {
   628         _cost_edges = new CostEdgeMap(*graph);
   629       }
   630       if (!_heap_cross_ref) {
   631         _heap_cross_ref = new HeapCrossRef(*graph, -1);
   632       }
   633       if (!_heap) {
   634         _heap = new Heap(*_heap_cross_ref);
   635       }
   636     }
   637 
   638     void destroyStructures() {
   639       if (local_arborescence) {
   640         delete _arborescence;
   641       }
   642       if (local_pred) {
   643         delete _pred;
   644       }
   645       if (!_edge_order) {
   646         delete _edge_order;
   647       }
   648       if (_node_order) {
   649         delete _node_order;
   650       }
   651       if (!_cost_edges) {
   652         delete _cost_edges;
   653       }
   654       if (!_heap) {
   655         delete _heap;
   656       }
   657       if (!_heap_cross_ref) {
   658         delete _heap_cross_ref;
   659       }
   660     }
   661 
   662     Edge prepare(Node node) {
   663       std::vector<Node> nodes;
   664       (*_node_order)[node] = _dual_node_list.size();
   665       StackLevel level;
   666       level.node_level = _dual_node_list.size();
   667       _dual_node_list.push_back(node);
   668       for (InEdgeIt it(*graph, node); it != INVALID; ++it) {
   669         Edge edge = it;
   670         Node source = graph->source(edge);
   671         Value value = (*cost)[it];
   672         if (source == node || (*_node_order)[source] == -3) continue;
   673         if ((*_cost_edges)[source].edge == INVALID) {
   674           (*_cost_edges)[source].edge = edge;
   675           (*_cost_edges)[source].value = value;
   676           nodes.push_back(source);
   677         } else {
   678           if ((*_cost_edges)[source].value > value) {
   679             (*_cost_edges)[source].edge = edge;
   680             (*_cost_edges)[source].value = value;
   681           }
   682         }      
   683       }
   684       CostEdge minimum = (*_cost_edges)[nodes[0]]; 
   685       for (int i = 1; i < (int)nodes.size(); ++i) {
   686         if ((*_cost_edges)[nodes[i]].value < minimum.value) {
   687           minimum = (*_cost_edges)[nodes[i]];
   688         }
   689       }
   690       _edge_order->set(minimum.edge, _dual_variables.size());
   691       DualVariable var(_dual_node_list.size() - 1, 
   692                        _dual_node_list.size(), minimum.value);
   693       _dual_variables.push_back(var);
   694       for (int i = 0; i < (int)nodes.size(); ++i) {
   695         (*_cost_edges)[nodes[i]].value -= minimum.value;
   696         level.edges.push_back((*_cost_edges)[nodes[i]]);
   697         (*_cost_edges)[nodes[i]].edge = INVALID;
   698       }
   699       level_stack.push_back(level);
   700       return minimum.edge;
   701     }
   702   
   703     Edge contract(Node node) {
   704       int node_bottom = bottom(node);
   705       std::vector<Node> nodes;
   706       while (!level_stack.empty() && 
   707              level_stack.back().node_level >= node_bottom) {
   708         for (int i = 0; i < (int)level_stack.back().edges.size(); ++i) {
   709           Edge edge = level_stack.back().edges[i].edge;
   710           Node source = graph->source(edge);
   711           Value value = level_stack.back().edges[i].value;
   712           if ((*_node_order)[source] >= node_bottom) continue;
   713           if ((*_cost_edges)[source].edge == INVALID) {
   714             (*_cost_edges)[source].edge = edge;
   715             (*_cost_edges)[source].value = value;
   716             nodes.push_back(source);
   717           } else {
   718             if ((*_cost_edges)[source].value > value) {
   719               (*_cost_edges)[source].edge = edge;
   720               (*_cost_edges)[source].value = value;
   721             }
   722           }
   723         }
   724         level_stack.pop_back();
   725       }
   726       CostEdge minimum = (*_cost_edges)[nodes[0]]; 
   727       for (int i = 1; i < (int)nodes.size(); ++i) {
   728         if ((*_cost_edges)[nodes[i]].value < minimum.value) {
   729           minimum = (*_cost_edges)[nodes[i]];
   730         }
   731       }
   732       _edge_order->set(minimum.edge, _dual_variables.size());
   733       DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
   734       _dual_variables.push_back(var);
   735       StackLevel level;
   736       level.node_level = node_bottom;
   737       for (int i = 0; i < (int)nodes.size(); ++i) {
   738         (*_cost_edges)[nodes[i]].value -= minimum.value;
   739         level.edges.push_back((*_cost_edges)[nodes[i]]);
   740         (*_cost_edges)[nodes[i]].edge = INVALID;
   741       }
   742       level_stack.push_back(level);
   743       return minimum.edge;
   744     }
   745 
   746     int bottom(Node node) {
   747       int k = level_stack.size() - 1;
   748       while (level_stack[k].node_level > (*_node_order)[node]) {
   749         --k;
   750       }
   751       return level_stack[k].node_level;
   752     }
   753 
   754     void finalize(Edge edge) {
   755       Node node = graph->target(edge);
   756       _heap->push(node, (*_edge_order)[edge]);
   757       _pred->set(node, edge);
   758       while (!_heap->empty()) {
   759         Node source = _heap->top();
   760         _heap->pop();
   761         _node_order->set(source, -1);
   762         for (OutEdgeIt it(*graph, source); it != INVALID; ++it) {
   763           if ((*_edge_order)[it] < 0) continue;
   764           Node target = graph->target(it);
   765           switch(_heap->state(target)) {
   766           case Heap::PRE_HEAP:
   767             _heap->push(target, (*_edge_order)[it]); 
   768             _pred->set(target, it);
   769             break;
   770           case Heap::IN_HEAP:
   771             if ((*_edge_order)[it] < (*_heap)[target]) {
   772               _heap->decrease(target, (*_edge_order)[it]); 
   773               _pred->set(target, it);
   774             }
   775             break;
   776           case Heap::POST_HEAP:
   777             break;
   778           }
   779         }
   780         _arborescence->set((*_pred)[source], true);
   781       }
   782     }
   783 
   784   };
   785 
   786   /// \ingroup spantree
   787   ///
   788   /// \brief Function type interface for MinCostArborescence algorithm.
   789   ///
   790   /// Function type interface for MinCostArborescence algorithm.
   791   /// \param graph The Graph that the algorithm runs on.
   792   /// \param cost The CostMap of the edges.
   793   /// \param source The source of the arborescence.
   794   /// \retval arborescence The bool EdgeMap which stores the arborescence.
   795   /// \return The cost of the arborescence. 
   796   ///
   797   /// \sa MinCostArborescence
   798   template <typename Graph, typename CostMap, typename ArborescenceMap>
   799   typename CostMap::Value minCostArborescence(const Graph& graph, 
   800                                               const CostMap& cost,
   801                                               typename Graph::Node source,
   802                                               ArborescenceMap& arborescence) {
   803     typename MinCostArborescence<Graph, CostMap>
   804       ::template DefArborescenceMap<ArborescenceMap>
   805       ::Create mca(graph, cost);
   806     mca.arborescenceMap(arborescence);
   807     mca.run(source);
   808     return mca.arborescenceValue();
   809   }
   810 
   811 }
   812 
   813 #endif