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