# HG changeset patch
# User deba
# Date 1143803444 0
# Node ID 93fcadf94ab05bc833f39079a132e5259cce6320
# Parent  4ab8a25def3cbaeb14a8a509e8fd0c67e0ed4b73
Bugfix in the minimum cost arborescence algorithm
Dual solution computation and interface for algorithm
Optimality test on random graph

diff -r 4ab8a25def3c -r 93fcadf94ab0 lemon/min_cost_arborescence.h
--- a/lemon/min_cost_arborescence.h	Thu Mar 30 15:34:56 2006 +0000
+++ b/lemon/min_cost_arborescence.h	Fri Mar 31 11:10:44 2006 +0000
@@ -26,6 +26,7 @@
 #include <vector>
 
 #include <lemon/list_graph.h>
+#include <lemon/bin_heap.h>
 
 namespace lemon {
 
@@ -56,11 +57,9 @@
     /// 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.   
+    /// It must meet the \ref concept::WritedMap "WriteMap" concept.
+    /// Initially it will be setted to false on each edge. After it
+    /// will set all arborescence edges once.
     typedef typename Graph::template EdgeMap<bool> ArborescenceMap; 
 
     /// \brief Instantiates a ArborescenceMap.
@@ -72,6 +71,20 @@
       return new ArborescenceMap(_graph);
     }
 
+    /// \brief The type of the PredMap
+    ///
+    /// The type of the PredMap. It is a node map with an edge value type.
+    typedef typename Graph::template NodeMap<typename Graph::Edge> PredMap;
+
+    /// \brief Instantiates a PredMap.
+    ///
+    /// This function instantiates a \ref PredMap. 
+    /// \param _graph is the graph, to which we would like to define the 
+    /// PredMap.
+    static PredMap *createPredMap(const Graph &_graph){
+      return new PredMap(_graph);
+    }
+    
   };
 
   /// \ingroup spantree
@@ -86,6 +99,9 @@
   /// given sources and spans all the nodes which are reachable from the
   /// sources. The time complexity of the algorithm is O(n^2 + e).
   ///
+  /// The algorithm provides also an optimal dual solution to arborescence
+  /// that way the optimality of the algorithm can be proofed easily.
+  ///
   /// \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 
@@ -135,6 +151,8 @@
     typedef typename Traits::CostMap CostMap;
     ///The type of the costs of the edges.
     typedef typename Traits::Value Value;
+    ///The type of the predecessor map.
+    typedef typename Traits::PredMap PredMap;
     ///The type of the map that stores which edges are in the arborescence.
     typedef typename Traits::ArborescenceMap ArborescenceMap;
 
@@ -157,14 +175,20 @@
 
     };
 
-    const Graph* graph;
-    const CostMap* cost;
+    const Graph *graph;
+    const CostMap *cost;
 
-    ArborescenceMap* _arborescence_map;
-    bool local_arborescence_map;
+    PredMap *_pred;
+    bool local_pred;
 
-    typedef typename Graph::template NodeMap<int> LevelMap;
-    LevelMap *_level;
+    ArborescenceMap *_arborescence;
+    bool local_arborescence;
+
+    typedef typename Graph::template EdgeMap<int> EdgeOrder;
+    EdgeOrder *_edge_order;
+    
+    typedef typename Graph::template NodeMap<int> NodeOrder;
+    NodeOrder *_node_order;
 
     typedef typename Graph::template NodeMap<CostEdge> CostEdgeMap;
     CostEdgeMap *_cost_edges; 
@@ -179,7 +203,30 @@
     std::vector<StackLevel> level_stack;    
     std::vector<Node> queue;
 
-    int node_counter;
+    typedef std::vector<typename Graph::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 Graph::template NodeMap<int> HeapCrossRef;
+    
+    HeapCrossRef *_heap_cross_ref;
+
+    typedef BinHeap<Node, int, HeapCrossRef> Heap;
+
+    Heap *_heap;
 
   public:
 
@@ -208,6 +255,26 @@
       typedef MinCostArborescence<Graph, CostMap, 
                                    DefArborescenceMapTraits<T> > Create;
     };
+
+    template <class T>
+    struct DefPredMapTraits : public Traits {
+      typedef T PredMap;
+      static PredMap *createPredMap(const Graph &)
+      {
+	throw UninitializedParameter();
+      }
+    };
+
+    /// \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<Graph, CostMap, DefPredMapTraits<T> > {
+      typedef MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > Create;
+    };
     
     /// @}
 
@@ -216,9 +283,10 @@
     /// \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) {}
+      : graph(&_graph), cost(&_cost), _pred(0), local_pred(false),
+        _arborescence(0), local_arborescence(false), 
+        _edge_order(0), _node_order(0), _cost_edges(0), 
+        _heap_cross_ref(0), _heap(0) {}
 
     /// \brief Destructor.
     ~MinCostArborescence() {
@@ -230,7 +298,24 @@
     /// Sets the arborescence map.
     /// \return \c (*this)
     MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
-      _arborescence_map = &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;
     }
 
@@ -246,7 +331,7 @@
     ///
     /// Returns a reference to the arborescence map.
     const ArborescenceMap& arborescenceMap() const {
-      return *_arborescence_map;
+      return *_arborescence;
     }
 
     /// \brief Returns true if the edge is in the arborescence.
@@ -254,23 +339,141 @@
     /// 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];
+    bool arborescence(Edge edge) const {
+      return (*_pred)[graph->target(edge)] == edge;
+    }
+
+    /// \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 edge of the given node.
+    ///
+    /// Returns the predecessor edge of the given node.
+    bool pred(Node node) const {
+      return (*_pred)[node];
     }
  
     /// \brief Returns the cost of the arborescence.
     ///
     /// Returns the cost of the arborescence.
-   Value arborescenceCost() const {
+    Value arborescenceValue() const {
       Value sum = 0;
       for (EdgeIt it(*graph); it != INVALID; ++it) {
-        if (arborescenceEdge(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 dualSize() 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), _variable(variable) 
+      {
+        _index = _algorithm->_dual_variables[_variable].begin;
+      }
+
+      /// \brief Invalid constructor.
+      ///
+      /// Invalid constructor.
+      DualIt(Invalid) : _algorithm(0) {}
+
+      /// \brief Conversion to node.
+      ///
+      /// Conversion to node.
+      operator Node() const { 
+        return _algorithm ? _algorithm->_dual_node_list[_index] : INVALID;
+      }
+
+      /// \brief Increment operator.
+      ///
+      /// Increment operator.
+      DualIt& operator++() {
+        ++_index;
+        if (_algorithm->_dual_variables[_variable].end == _index) {
+          _algorithm = 0;
+        }
+        return *this; 
+      }
+
+      bool operator==(const DualIt& it) const { 
+        return (Node)(*this) == (Node)it; 
+      }
+      bool operator!=(const DualIt& it) const { 
+        return (Node)(*this) != (Node)it; 
+      }
+      bool operator<(const DualIt& it) const { 
+        return (Node)(*this) < (Node)it; 
+      }
+      
+    private:
+      const MinCostArborescence* _algorithm;
+      int _variable;
+      int _index;
+    };
+
     /// @}
     
     /// \name Execution control
@@ -279,7 +482,7 @@
     /// 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
+    /// Finally \ref start() will perform the arborescence
     /// computation.
 
     ///@{
@@ -290,13 +493,18 @@
     ///
     void init() {
       initStructures();
+      _heap->clear();
       for (NodeIt it(*graph); it != INVALID; ++it) {
         (*_cost_edges)[it].edge = INVALID;
-        (*_level)[it] = -3; 
+        _node_order->set(it, -3); 
+        _heap_cross_ref->set(it, Heap::PRE_HEAP);
       }
       for (EdgeIt it(*graph); it != INVALID; ++it) {
-        _arborescence_map->set(it, false);
+        _arborescence->set(it, false);
+        _edge_order->set(it, -1);
       }
+      _dual_node_list.clear();
+      _dual_variables.clear();
     }
 
     /// \brief Adds a new source node.
@@ -309,14 +517,15 @@
         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));
+          Node target = graph->target(it);
+          if ((*_node_order)[target] == -3) {
+            (*_node_order)[target] = -2;
+            nodes.push_back(target);
+            queue.push_back(target);
           }
         }
       }
-      (*_level)[source] = -1;
+      (*_node_order)[source] = -1;
     }
 
     /// \brief Processes the next node in the priority queue.
@@ -327,19 +536,20 @@
     ///
     /// \warning The queue must not be empty!
     Node processNextNode() {
-      node_counter = 0;
       Node node = queue.back();
       queue.pop_back();
-      if ((*_level)[node] == -2) {
+      if ((*_node_order)[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)]));
+        Node source = graph->source(edge);
+        while ((*_node_order)[source] != -1) {
+          if ((*_node_order)[source] >= 0) {
+            edge = contract(source);
           } else {
-            edge = prepare(graph->source(edge));
+            edge = prepare(source);
           }
+          source = graph->source(edge);
         }
-        finalize(graph->target(edge));
+        finalize(edge);
         level_stack.clear();        
       }
       return node;
@@ -400,46 +610,74 @@
   protected:
 
     void initStructures() {
-      if (!_arborescence_map) {
-        local_arborescence_map = true;
-        _arborescence_map = Traits::createArborescenceMap(*graph);
+      if (!_pred) {
+        local_pred = true;
+        _pred = Traits::createPredMap(*graph);
       }
-      if (!_level) {
-        _level = new LevelMap(*graph);
+      if (!_arborescence) {
+        local_arborescence = true;
+        _arborescence = Traits::createArborescenceMap(*graph);
+      }
+      if (!_edge_order) {
+        _edge_order = new EdgeOrder(*graph);
+      }
+      if (!_node_order) {
+        _node_order = new NodeOrder(*graph);
       }
       if (!_cost_edges) {
         _cost_edges = new CostEdgeMap(*graph);
       }
+      if (!_heap_cross_ref) {
+        _heap_cross_ref = new HeapCrossRef(*graph, -1);
+      }
+      if (!_heap) {
+        _heap = new Heap(*_heap_cross_ref);
+      }
     }
 
     void destroyStructures() {
-      if (_level) {
-        delete _level;
+      if (local_arborescence) {
+        delete _arborescence;
+      }
+      if (local_pred) {
+        delete _pred;
+      }
+      if (!_edge_order) {
+        delete _edge_order;
+      }
+      if (_node_order) {
+        delete _node_order;
       }
       if (!_cost_edges) {
         delete _cost_edges;
       }
-      if (local_arborescence_map) {
-        delete _arborescence_map;
+      if (!_heap) {
+        delete _heap;
+      }
+      if (!_heap_cross_ref) {
+        delete _heap_cross_ref;
       }
     }
 
     Edge prepare(Node node) {
       std::vector<Node> nodes;
-      (*_level)[node] = node_counter;
+      (*_node_order)[node] = _dual_node_list.size();
+      StackLevel level;
+      level.node_level = _dual_node_list.size();
+      _dual_node_list.push_back(node);
       for (InEdgeIt it(*graph, node); it != INVALID; ++it) {
         Edge edge = it;
+        Node source = graph->source(edge);
         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));
+        if (source == node || (*_node_order)[source] == -3) continue;
+        if ((*_cost_edges)[source].edge == INVALID) {
+          (*_cost_edges)[source].edge = edge;
+          (*_cost_edges)[source].value = value;
+          nodes.push_back(source);
         } else {
-          if ((*_cost_edges)[graph->source(edge)].value > value) {
-            (*_cost_edges)[graph->source(edge)].edge = edge;
-            (*_cost_edges)[graph->source(edge)].value = value;
+          if ((*_cost_edges)[source].value > value) {
+            (*_cost_edges)[source].edge = edge;
+            (*_cost_edges)[source].value = value;
           }
         }      
       }
@@ -449,35 +687,37 @@
           minimum = (*_cost_edges)[nodes[i]];
         }
       }
-      StackLevel level;
-      level.node_level = node_counter;
+      _edge_order->set(minimum.edge, _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_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) {
+    Edge 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().edges.size(); ++i) {
           Edge edge = level_stack.back().edges[i].edge;
+          Node source = graph->source(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));
+          if ((*_node_order)[source] >= node_bottom) continue;
+          if ((*_cost_edges)[source].edge == INVALID) {
+            (*_cost_edges)[source].edge = edge;
+            (*_cost_edges)[source].value = value;
+            nodes.push_back(source);
           } else {
-            if ((*_cost_edges)[graph->source(edge)].value > value) {
-              (*_cost_edges)[graph->source(edge)].edge = edge;
-              (*_cost_edges)[graph->source(edge)].value = value;
+            if ((*_cost_edges)[source].value > value) {
+              (*_cost_edges)[source].edge = edge;
+              (*_cost_edges)[source].value = value;
             }
           }
         }
@@ -489,6 +729,9 @@
           minimum = (*_cost_edges)[nodes[i]];
         }
       }
+      _edge_order->set(minimum.edge, _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) {
@@ -497,34 +740,45 @@
         (*_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 bottom(Node node) {
       int k = level_stack.size() - 1;
-      while (level_stack[k].node_level > level) {
+      while (level_stack[k].node_level > (*_node_order)[node]) {
         --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);
+    void finalize(Edge edge) {
+      Node node = graph->target(edge);
+      _heap->push(node, (*_edge_order)[edge]);
+      _pred->set(node, edge);
+      while (!_heap->empty()) {
+        Node source = _heap->top();
+        _heap->pop();
+        _node_order->set(source, -1);
+        for (OutEdgeIt it(*graph, source); it != INVALID; ++it) {
+          if ((*_edge_order)[it] < 0) continue;
+          Node target = graph->target(it);
+          switch(_heap->state(target)) {
+          case Heap::PRE_HEAP:
+            _heap->push(target, (*_edge_order)[it]); 
+            _pred->set(target, it);
+            break;
+          case Heap::IN_HEAP:
+            if ((*_edge_order)[it] < (*_heap)[target]) {
+              _heap->decrease(target, (*_edge_order)[it]); 
+              _pred->set(target, it);
+            }
+            break;
+          case Heap::POST_HEAP:
+            break;
           }
         }
+        _arborescence->set((*_pred)[source], true);
       }
-      (*_level)[source] = -1;      
     }
 
   };
@@ -551,11 +805,9 @@
       ::Create mca(graph, cost);
     mca.arborescenceMap(arborescence);
     mca.run(source);
-    return mca.arborescenceCost();
+    return mca.arborescenceValue();
   }
 
 }
 
 #endif
-
-// Hilbert - Huang
diff -r 4ab8a25def3c -r 93fcadf94ab0 test/Makefile.am
--- a/test/Makefile.am	Thu Mar 30 15:34:56 2006 +0000
+++ b/test/Makefile.am	Fri Mar 31 11:10:44 2006 +0000
@@ -13,6 +13,7 @@
 
 check_PROGRAMS = \
 	all_pairs_shortest_path_test \
+	arborescence_test \
 	bfs_test \
 	counter_test \
 	dfs_test \
@@ -52,6 +53,7 @@
 XFAIL_TESTS = test_tools_fail$(EXEEXT)
 
 all_pairs_shortest_path_test_SOURCES = all_pairs_shortest_path_test.cc
+arborescence_test_SOURCES = arborescence_test.cc
 bfs_test_SOURCES = bfs_test.cc
 counter_test_SOURCES = counter_test.cc
 dfs_test_SOURCES = dfs_test.cc
diff -r 4ab8a25def3c -r 93fcadf94ab0 test/arborescence_test.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/arborescence_test.cc	Fri Mar 31 11:10:44 2006 +0000
@@ -0,0 +1,104 @@
+#include <iostream>
+#include <set>
+#include <vector>
+#include <iterator>
+
+#include <cmath>
+#include <cstdlib>
+
+#include <lemon/smart_graph.h>
+#include <lemon/min_cost_arborescence.h>
+
+#include <lemon/graph_utils.h>
+#include <lemon/time_measure.h>
+
+#include <lemon/tolerance.h>
+
+using namespace lemon;
+using namespace std;
+
+int main(int argc, const char *argv[]) {
+  srand(time(0));
+  typedef SmartGraph Graph;
+  GRAPH_TYPEDEFS(Graph);
+
+  typedef Graph::EdgeMap<double> CostMap;
+
+  const int n = argc > 1 ? atoi(argv[1]) : 100;
+  const int e = argc > 2 ? atoi(argv[2]) : (int)(n * log(n));
+
+  Graph graph;
+  CostMap cost(graph);
+  vector<Node> nodes;
+  
+  for (int i = 0; i < n; ++i) {
+    nodes.push_back(graph.addNode());
+  }
+
+  for (int i = 0; i < e; ++i) {
+    int s = (int)(n * (double)rand() / (RAND_MAX + 1.0));
+    int t = (int)(n * (double)rand() / (RAND_MAX + 1.0));
+    double c = rand() / (1.0 + RAND_MAX) * 100.0 + 20.0;
+    Edge edge = graph.addEdge(nodes[s], nodes[t]);
+    cost[edge] = c;
+  }
+
+  Node source = nodes[(int)(n * (double)rand() / (RAND_MAX + 1.0))];
+
+  MinCostArborescence<Graph, CostMap> mca(graph, cost);
+  mca.run(source);
+
+  vector<pair<double, set<Node> > > dualSolution(mca.dualSize());
+
+  for (int i = 0; i < mca.dualSize(); ++i) {
+    dualSolution[i].first = mca.dualValue(i);
+    for (MinCostArborescence<Graph, CostMap>::DualIt it(mca, i); 
+         it != INVALID; ++it) {
+      dualSolution[i].second.insert(it);
+    }
+  }
+
+  Tolerance<double> tolerance;
+
+  for (EdgeIt it(graph); it != INVALID; ++it) {
+    if (mca.reached(graph.source(it))) {
+      double sum = 0.0;
+      for (int i = 0; i < (int)dualSolution.size(); ++i) {
+        if (dualSolution[i].second.find(graph.target(it)) 
+            != dualSolution[i].second.end() &&
+            dualSolution[i].second.find(graph.source(it)) 
+            == dualSolution[i].second.end()) {
+          sum += dualSolution[i].first;
+        }
+      }
+      if (mca.arborescence(it)) {
+        LEMON_ASSERT(!tolerance.less(sum, cost[it]), "INVALID DUAL");
+      }
+      LEMON_ASSERT(!tolerance.less(cost[it], sum), "INVALID DUAL");
+    }
+  }
+
+
+  LEMON_ASSERT(!tolerance.different(mca.dualValue(), mca.arborescenceValue()),
+               "INVALID DUAL");
+
+
+  LEMON_ASSERT(mca.reached(source), "INVALID ARBORESCENCE");
+  for (EdgeIt it(graph); it != INVALID; ++it) {
+    LEMON_ASSERT(!mca.reached(graph.source(it)) || 
+                 mca.reached(graph.target(it)), "INVALID ARBORESCENCE");
+  }
+
+  for (NodeIt it(graph); it != INVALID; ++it) {
+    if (!mca.reached(it)) continue;
+    int cnt = 0;
+    for (InEdgeIt jt(graph, it); jt != INVALID; ++jt) {
+      if (mca.arborescence(jt)) {
+        ++cnt;
+      }
+    }
+    LEMON_ASSERT((it == source ? cnt == 0 : cnt == 1), "INVALID ARBORESCENCE");
+  }
+  
+  return 0;
+}