deba@512: /* -*- mode: C++; indent-tabs-mode: nil; -*-
deba@512:  *
deba@512:  * This file is a part of LEMON, a generic C++ optimization library.
deba@512:  *
deba@512:  * Copyright (C) 2003-2008
deba@512:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@512:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@512:  *
deba@512:  * Permission to use, modify and distribute this software is granted
deba@512:  * provided that this copyright notice appears in all copies. For
deba@512:  * precise terms see the accompanying LICENSE file.
deba@512:  *
deba@512:  * This software is provided "AS IS" with no warranty of any kind,
deba@512:  * express or implied, and with no claim as to its suitability for any
deba@512:  * purpose.
deba@512:  *
deba@512:  */
deba@512: 
deba@512: #ifndef LEMON_MIN_COST_ARBORESCENCE_H
deba@512: #define LEMON_MIN_COST_ARBORESCENCE_H
deba@512: 
deba@512: ///\ingroup spantree
deba@512: ///\file
deba@512: ///\brief Minimum Cost Arborescence algorithm.
deba@512: 
deba@512: #include <vector>
deba@512: 
deba@512: #include <lemon/list_graph.h>
deba@512: #include <lemon/bin_heap.h>
deba@512: #include <lemon/assert.h>
deba@512: 
deba@512: namespace lemon {
deba@512: 
deba@512: 
deba@512:   /// \brief Default traits class for MinCostArborescence class.
deba@512:   ///
deba@512:   /// Default traits class for MinCostArborescence class.
kpeter@550:   /// \param GR Digraph type.
kpeter@617:   /// \param CM Type of the cost map.
kpeter@550:   template <class GR, class CM>
deba@512:   struct MinCostArborescenceDefaultTraits{
deba@512: 
deba@512:     /// \brief The digraph type the algorithm runs on.
kpeter@550:     typedef GR Digraph;
deba@512: 
deba@512:     /// \brief The type of the map that stores the arc costs.
deba@512:     ///
deba@512:     /// The type of the map that stores the arc costs.
kpeter@617:     /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
kpeter@550:     typedef CM CostMap;
deba@512: 
deba@512:     /// \brief The value type of the costs.
deba@512:     ///
deba@512:     /// The value type of the costs.
deba@512:     typedef typename CostMap::Value Value;
deba@512: 
deba@512:     /// \brief The type of the map that stores which arcs are in the
deba@512:     /// arborescence.
deba@512:     ///
deba@512:     /// The type of the map that stores which arcs are in the
kpeter@617:     /// arborescence.  It must conform to the \ref concepts::WriteMap
kpeter@617:     /// "WriteMap" concept, and its value type must be \c bool
kpeter@617:     /// (or convertible). Initially it will be set to \c false on each
kpeter@617:     /// arc, then it will be set on each arborescence arc once.
deba@512:     typedef typename Digraph::template ArcMap<bool> ArborescenceMap;
deba@512: 
kpeter@550:     /// \brief Instantiates a \c ArborescenceMap.
deba@512:     ///
kpeter@550:     /// This function instantiates a \c ArborescenceMap.
kpeter@617:     /// \param digraph The digraph to which we would like to calculate
kpeter@617:     /// the \c ArborescenceMap.
deba@512:     static ArborescenceMap *createArborescenceMap(const Digraph &digraph){
deba@512:       return new ArborescenceMap(digraph);
deba@512:     }
deba@512: 
kpeter@550:     /// \brief The type of the \c PredMap
deba@512:     ///
kpeter@617:     /// The type of the \c PredMap. It must confrom to the
kpeter@617:     /// \ref concepts::WriteMap "WriteMap" concept, and its value type
kpeter@617:     /// must be the \c Arc type of the digraph.
deba@512:     typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
deba@512: 
kpeter@550:     /// \brief Instantiates a \c PredMap.
deba@512:     ///
kpeter@550:     /// This function instantiates a \c PredMap.
kpeter@550:     /// \param digraph The digraph to which we would like to define the
kpeter@550:     /// \c PredMap.
deba@512:     static PredMap *createPredMap(const Digraph &digraph){
deba@512:       return new PredMap(digraph);
deba@512:     }
deba@512: 
deba@512:   };
deba@512: 
deba@512:   /// \ingroup spantree
deba@512:   ///
kpeter@576:   /// \brief Minimum Cost Arborescence algorithm class.
deba@512:   ///
kpeter@617:   /// This class provides an efficient implementation of the
kpeter@576:   /// Minimum Cost Arborescence algorithm. The arborescence is a tree
deba@512:   /// which is directed from a given source node of the digraph. One or
kpeter@617:   /// more sources should be given to the algorithm and it will calculate
kpeter@617:   /// the minimum cost subgraph that is the union of arborescences with the
deba@512:   /// given sources and spans all the nodes which are reachable from the
kpeter@550:   /// sources. The time complexity of the algorithm is O(n<sup>2</sup>+e).
deba@512:   ///
kpeter@617:   /// The algorithm also provides an optimal dual solution, therefore
deba@512:   /// the optimality of the solution can be checked.
deba@512:   ///
kpeter@617:   /// \param GR The digraph type the algorithm runs on.
kpeter@617:   /// \param CM A read-only arc map storing the costs of the
deba@512:   /// arcs. It is read once for each arc, so the map may involve in
kpeter@617:   /// relatively time consuming process to compute the arc costs if
deba@512:   /// it is necessary. The default map type is \ref
deba@512:   /// concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
kpeter@550:   /// \param TR Traits class to set various data types used
deba@512:   /// by the algorithm. The default traits class is
deba@512:   /// \ref MinCostArborescenceDefaultTraits
kpeter@617:   /// "MinCostArborescenceDefaultTraits<GR, CM>".
deba@512: #ifndef DOXYGEN
kpeter@617:   template <typename GR,
kpeter@550:             typename CM = typename GR::template ArcMap<int>,
kpeter@550:             typename TR =
kpeter@550:               MinCostArborescenceDefaultTraits<GR, CM> >
deba@512: #else
kpeter@550:   template <typename GR, typename CM, typedef TR>
deba@512: #endif
deba@512:   class MinCostArborescence {
deba@512:   public:
deba@512: 
kpeter@617:     /// \brief The \ref MinCostArborescenceDefaultTraits "traits class" 
kpeter@617:     /// of the algorithm. 
kpeter@550:     typedef TR Traits;
deba@512:     /// The type of the underlying digraph.
deba@512:     typedef typename Traits::Digraph Digraph;
deba@512:     /// The type of the map that stores the arc costs.
deba@512:     typedef typename Traits::CostMap CostMap;
deba@512:     ///The type of the costs of the arcs.
deba@512:     typedef typename Traits::Value Value;
deba@512:     ///The type of the predecessor map.
deba@512:     typedef typename Traits::PredMap PredMap;
deba@512:     ///The type of the map that stores which arcs are in the arborescence.
deba@512:     typedef typename Traits::ArborescenceMap ArborescenceMap;
deba@512: 
deba@512:     typedef MinCostArborescence Create;
deba@512: 
deba@512:   private:
deba@512: 
deba@512:     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
deba@512: 
deba@512:     struct CostArc {
deba@512: 
deba@512:       Arc arc;
deba@512:       Value value;
deba@512: 
deba@512:       CostArc() {}
deba@512:       CostArc(Arc _arc, Value _value) : arc(_arc), value(_value) {}
deba@512: 
deba@512:     };
deba@512: 
deba@512:     const Digraph *_digraph;
deba@512:     const CostMap *_cost;
deba@512: 
deba@512:     PredMap *_pred;
deba@512:     bool local_pred;
deba@512: 
deba@512:     ArborescenceMap *_arborescence;
deba@512:     bool local_arborescence;
deba@512: 
deba@512:     typedef typename Digraph::template ArcMap<int> ArcOrder;
deba@512:     ArcOrder *_arc_order;
deba@512: 
deba@512:     typedef typename Digraph::template NodeMap<int> NodeOrder;
deba@512:     NodeOrder *_node_order;
deba@512: 
deba@512:     typedef typename Digraph::template NodeMap<CostArc> CostArcMap;
deba@512:     CostArcMap *_cost_arcs;
deba@512: 
deba@512:     struct StackLevel {
deba@512: 
deba@512:       std::vector<CostArc> arcs;
deba@512:       int node_level;
deba@512: 
deba@512:     };
deba@512: 
deba@512:     std::vector<StackLevel> level_stack;
deba@512:     std::vector<Node> queue;
deba@512: 
deba@512:     typedef std::vector<typename Digraph::Node> DualNodeList;
deba@512: 
deba@512:     DualNodeList _dual_node_list;
deba@512: 
deba@512:     struct DualVariable {
deba@512:       int begin, end;
deba@512:       Value value;
deba@512: 
deba@512:       DualVariable(int _begin, int _end, Value _value)
deba@512:         : begin(_begin), end(_end), value(_value) {}
deba@512: 
deba@512:     };
deba@512: 
deba@512:     typedef std::vector<DualVariable> DualVariables;
deba@512: 
deba@512:     DualVariables _dual_variables;
deba@512: 
deba@512:     typedef typename Digraph::template NodeMap<int> HeapCrossRef;
deba@512: 
deba@512:     HeapCrossRef *_heap_cross_ref;
deba@512: 
deba@512:     typedef BinHeap<int, HeapCrossRef> Heap;
deba@512: 
deba@512:     Heap *_heap;
deba@512: 
deba@512:   protected:
deba@512: 
deba@512:     MinCostArborescence() {}
deba@512: 
deba@512:   private:
deba@512: 
deba@512:     void createStructures() {
deba@512:       if (!_pred) {
deba@512:         local_pred = true;
deba@512:         _pred = Traits::createPredMap(*_digraph);
deba@512:       }
deba@512:       if (!_arborescence) {
deba@512:         local_arborescence = true;
deba@512:         _arborescence = Traits::createArborescenceMap(*_digraph);
deba@512:       }
deba@512:       if (!_arc_order) {
deba@512:         _arc_order = new ArcOrder(*_digraph);
deba@512:       }
deba@512:       if (!_node_order) {
deba@512:         _node_order = new NodeOrder(*_digraph);
deba@512:       }
deba@512:       if (!_cost_arcs) {
deba@512:         _cost_arcs = new CostArcMap(*_digraph);
deba@512:       }
deba@512:       if (!_heap_cross_ref) {
deba@512:         _heap_cross_ref = new HeapCrossRef(*_digraph, -1);
deba@512:       }
deba@512:       if (!_heap) {
deba@512:         _heap = new Heap(*_heap_cross_ref);
deba@512:       }
deba@512:     }
deba@512: 
deba@512:     void destroyStructures() {
deba@512:       if (local_arborescence) {
deba@512:         delete _arborescence;
deba@512:       }
deba@512:       if (local_pred) {
deba@512:         delete _pred;
deba@512:       }
deba@512:       if (_arc_order) {
deba@512:         delete _arc_order;
deba@512:       }
deba@512:       if (_node_order) {
deba@512:         delete _node_order;
deba@512:       }
deba@512:       if (_cost_arcs) {
deba@512:         delete _cost_arcs;
deba@512:       }
deba@512:       if (_heap) {
deba@512:         delete _heap;
deba@512:       }
deba@512:       if (_heap_cross_ref) {
deba@512:         delete _heap_cross_ref;
deba@512:       }
deba@512:     }
deba@512: 
deba@512:     Arc prepare(Node node) {
deba@512:       std::vector<Node> nodes;
deba@512:       (*_node_order)[node] = _dual_node_list.size();
deba@512:       StackLevel level;
deba@512:       level.node_level = _dual_node_list.size();
deba@512:       _dual_node_list.push_back(node);
deba@512:       for (InArcIt it(*_digraph, node); it != INVALID; ++it) {
deba@512:         Arc arc = it;
deba@512:         Node source = _digraph->source(arc);
deba@512:         Value value = (*_cost)[it];
deba@512:         if (source == node || (*_node_order)[source] == -3) continue;
deba@512:         if ((*_cost_arcs)[source].arc == INVALID) {
deba@512:           (*_cost_arcs)[source].arc = arc;
deba@512:           (*_cost_arcs)[source].value = value;
deba@512:           nodes.push_back(source);
deba@512:         } else {
deba@512:           if ((*_cost_arcs)[source].value > value) {
deba@512:             (*_cost_arcs)[source].arc = arc;
deba@512:             (*_cost_arcs)[source].value = value;
deba@512:           }
deba@512:         }
deba@512:       }
deba@512:       CostArc minimum = (*_cost_arcs)[nodes[0]];
deba@512:       for (int i = 1; i < int(nodes.size()); ++i) {
deba@512:         if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
deba@512:           minimum = (*_cost_arcs)[nodes[i]];
deba@512:         }
deba@512:       }
kpeter@573:       (*_arc_order)[minimum.arc] = _dual_variables.size();
deba@512:       DualVariable var(_dual_node_list.size() - 1,
deba@512:                        _dual_node_list.size(), minimum.value);
deba@512:       _dual_variables.push_back(var);
deba@512:       for (int i = 0; i < int(nodes.size()); ++i) {
deba@512:         (*_cost_arcs)[nodes[i]].value -= minimum.value;
deba@512:         level.arcs.push_back((*_cost_arcs)[nodes[i]]);
deba@512:         (*_cost_arcs)[nodes[i]].arc = INVALID;
deba@512:       }
deba@512:       level_stack.push_back(level);
deba@512:       return minimum.arc;
deba@512:     }
deba@512: 
deba@512:     Arc contract(Node node) {
deba@512:       int node_bottom = bottom(node);
deba@512:       std::vector<Node> nodes;
deba@512:       while (!level_stack.empty() &&
deba@512:              level_stack.back().node_level >= node_bottom) {
deba@512:         for (int i = 0; i < int(level_stack.back().arcs.size()); ++i) {
deba@512:           Arc arc = level_stack.back().arcs[i].arc;
deba@512:           Node source = _digraph->source(arc);
deba@512:           Value value = level_stack.back().arcs[i].value;
deba@512:           if ((*_node_order)[source] >= node_bottom) continue;
deba@512:           if ((*_cost_arcs)[source].arc == INVALID) {
deba@512:             (*_cost_arcs)[source].arc = arc;
deba@512:             (*_cost_arcs)[source].value = value;
deba@512:             nodes.push_back(source);
deba@512:           } else {
deba@512:             if ((*_cost_arcs)[source].value > value) {
deba@512:               (*_cost_arcs)[source].arc = arc;
deba@512:               (*_cost_arcs)[source].value = value;
deba@512:             }
deba@512:           }
deba@512:         }
deba@512:         level_stack.pop_back();
deba@512:       }
deba@512:       CostArc minimum = (*_cost_arcs)[nodes[0]];
deba@512:       for (int i = 1; i < int(nodes.size()); ++i) {
deba@512:         if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
deba@512:           minimum = (*_cost_arcs)[nodes[i]];
deba@512:         }
deba@512:       }
kpeter@573:       (*_arc_order)[minimum.arc] = _dual_variables.size();
deba@512:       DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
deba@512:       _dual_variables.push_back(var);
deba@512:       StackLevel level;
deba@512:       level.node_level = node_bottom;
deba@512:       for (int i = 0; i < int(nodes.size()); ++i) {
deba@512:         (*_cost_arcs)[nodes[i]].value -= minimum.value;
deba@512:         level.arcs.push_back((*_cost_arcs)[nodes[i]]);
deba@512:         (*_cost_arcs)[nodes[i]].arc = INVALID;
deba@512:       }
deba@512:       level_stack.push_back(level);
deba@512:       return minimum.arc;
deba@512:     }
deba@512: 
deba@512:     int bottom(Node node) {
deba@512:       int k = level_stack.size() - 1;
deba@512:       while (level_stack[k].node_level > (*_node_order)[node]) {
deba@512:         --k;
deba@512:       }
deba@512:       return level_stack[k].node_level;
deba@512:     }
deba@512: 
deba@512:     void finalize(Arc arc) {
deba@512:       Node node = _digraph->target(arc);
deba@512:       _heap->push(node, (*_arc_order)[arc]);
deba@512:       _pred->set(node, arc);
deba@512:       while (!_heap->empty()) {
deba@512:         Node source = _heap->top();
deba@512:         _heap->pop();
kpeter@573:         (*_node_order)[source] = -1;
deba@512:         for (OutArcIt it(*_digraph, source); it != INVALID; ++it) {
deba@512:           if ((*_arc_order)[it] < 0) continue;
deba@512:           Node target = _digraph->target(it);
deba@512:           switch(_heap->state(target)) {
deba@512:           case Heap::PRE_HEAP:
deba@512:             _heap->push(target, (*_arc_order)[it]);
deba@512:             _pred->set(target, it);
deba@512:             break;
deba@512:           case Heap::IN_HEAP:
deba@512:             if ((*_arc_order)[it] < (*_heap)[target]) {
deba@512:               _heap->decrease(target, (*_arc_order)[it]);
deba@512:               _pred->set(target, it);
deba@512:             }
deba@512:             break;
deba@512:           case Heap::POST_HEAP:
deba@512:             break;
deba@512:           }
deba@512:         }
deba@512:         _arborescence->set((*_pred)[source], true);
deba@512:       }
deba@512:     }
deba@512: 
deba@512: 
deba@512:   public:
deba@512: 
kpeter@576:     /// \name Named Template Parameters
deba@512: 
deba@512:     /// @{
deba@512: 
deba@512:     template <class T>
kpeter@617:     struct SetArborescenceMapTraits : public Traits {
deba@512:       typedef T ArborescenceMap;
deba@512:       static ArborescenceMap *createArborescenceMap(const Digraph &)
deba@512:       {
deba@512:         LEMON_ASSERT(false, "ArborescenceMap is not initialized");
deba@512:         return 0; // ignore warnings
deba@512:       }
deba@512:     };
deba@512: 
deba@512:     /// \brief \ref named-templ-param "Named parameter" for
kpeter@617:     /// setting \c ArborescenceMap type
deba@512:     ///
deba@512:     /// \ref named-templ-param "Named parameter" for setting
kpeter@617:     /// \c ArborescenceMap type.
kpeter@617:     /// It must conform to the \ref concepts::WriteMap "WriteMap" concept,
kpeter@617:     /// and its value type must be \c bool (or convertible).
kpeter@617:     /// Initially it will be set to \c false on each arc,
kpeter@617:     /// then it will be set on each arborescence arc once.
deba@512:     template <class T>
kpeter@617:     struct SetArborescenceMap
deba@512:       : public MinCostArborescence<Digraph, CostMap,
kpeter@617:                                    SetArborescenceMapTraits<T> > {
deba@512:     };
deba@512: 
deba@512:     template <class T>
kpeter@617:     struct SetPredMapTraits : public Traits {
deba@512:       typedef T PredMap;
deba@512:       static PredMap *createPredMap(const Digraph &)
deba@512:       {
deba@512:         LEMON_ASSERT(false, "PredMap is not initialized");
kpeter@617:         return 0; // ignore warnings
deba@512:       }
deba@512:     };
deba@512: 
deba@512:     /// \brief \ref named-templ-param "Named parameter" for
kpeter@617:     /// setting \c PredMap type
deba@512:     ///
deba@512:     /// \ref named-templ-param "Named parameter" for setting
kpeter@617:     /// \c PredMap type.
kpeter@617:     /// It must meet the \ref concepts::WriteMap "WriteMap" concept, 
kpeter@617:     /// and its value type must be the \c Arc type of the digraph.
deba@512:     template <class T>
kpeter@617:     struct SetPredMap
kpeter@617:       : public MinCostArborescence<Digraph, CostMap, SetPredMapTraits<T> > {
deba@512:     };
deba@512: 
deba@512:     /// @}
deba@512: 
deba@512:     /// \brief Constructor.
deba@512:     ///
kpeter@550:     /// \param digraph The digraph the algorithm will run on.
kpeter@550:     /// \param cost The cost map used by the algorithm.
deba@512:     MinCostArborescence(const Digraph& digraph, const CostMap& cost)
deba@512:       : _digraph(&digraph), _cost(&cost), _pred(0), local_pred(false),
deba@512:         _arborescence(0), local_arborescence(false),
deba@512:         _arc_order(0), _node_order(0), _cost_arcs(0),
deba@512:         _heap_cross_ref(0), _heap(0) {}
deba@512: 
deba@512:     /// \brief Destructor.
deba@512:     ~MinCostArborescence() {
deba@512:       destroyStructures();
deba@512:     }
deba@512: 
deba@512:     /// \brief Sets the arborescence map.
deba@512:     ///
deba@512:     /// Sets the arborescence map.
kpeter@550:     /// \return <tt>(*this)</tt>
deba@512:     MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
deba@512:       if (local_arborescence) {
deba@512:         delete _arborescence;
deba@512:       }
deba@512:       local_arborescence = false;
deba@512:       _arborescence = &m;
deba@512:       return *this;
deba@512:     }
deba@512: 
kpeter@617:     /// \brief Sets the predecessor map.
deba@512:     ///
kpeter@617:     /// Sets the predecessor map.
kpeter@550:     /// \return <tt>(*this)</tt>
deba@512:     MinCostArborescence& predMap(PredMap& m) {
deba@512:       if (local_pred) {
deba@512:         delete _pred;
deba@512:       }
deba@512:       local_pred = false;
deba@512:       _pred = &m;
deba@512:       return *this;
deba@512:     }
deba@512: 
kpeter@576:     /// \name Execution Control
deba@512:     /// The simplest way to execute the algorithm is to use
deba@512:     /// one of the member functions called \c run(...). \n
deba@512:     /// If you need more control on the execution,
deba@512:     /// first you must call \ref init(), then you can add several
deba@512:     /// source nodes with \ref addSource().
deba@512:     /// Finally \ref start() will perform the arborescence
deba@512:     /// computation.
deba@512: 
deba@512:     ///@{
deba@512: 
deba@512:     /// \brief Initializes the internal data structures.
deba@512:     ///
deba@512:     /// Initializes the internal data structures.
deba@512:     ///
deba@512:     void init() {
deba@512:       createStructures();
deba@512:       _heap->clear();
deba@512:       for (NodeIt it(*_digraph); it != INVALID; ++it) {
deba@512:         (*_cost_arcs)[it].arc = INVALID;
kpeter@573:         (*_node_order)[it] = -3;
kpeter@573:         (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
deba@512:         _pred->set(it, INVALID);
deba@512:       }
deba@512:       for (ArcIt it(*_digraph); it != INVALID; ++it) {
deba@512:         _arborescence->set(it, false);
kpeter@573:         (*_arc_order)[it] = -1;
deba@512:       }
deba@512:       _dual_node_list.clear();
deba@512:       _dual_variables.clear();
deba@512:     }
deba@512: 
deba@512:     /// \brief Adds a new source node.
deba@512:     ///
deba@512:     /// Adds a new source node to the algorithm.
deba@512:     void addSource(Node source) {
deba@512:       std::vector<Node> nodes;
deba@512:       nodes.push_back(source);
deba@512:       while (!nodes.empty()) {
deba@512:         Node node = nodes.back();
deba@512:         nodes.pop_back();
deba@512:         for (OutArcIt it(*_digraph, node); it != INVALID; ++it) {
deba@512:           Node target = _digraph->target(it);
deba@512:           if ((*_node_order)[target] == -3) {
deba@512:             (*_node_order)[target] = -2;
deba@512:             nodes.push_back(target);
deba@512:             queue.push_back(target);
deba@512:           }
deba@512:         }
deba@512:       }
deba@512:       (*_node_order)[source] = -1;
deba@512:     }
deba@512: 
deba@512:     /// \brief Processes the next node in the priority queue.
deba@512:     ///
deba@512:     /// Processes the next node in the priority queue.
deba@512:     ///
deba@512:     /// \return The processed node.
deba@512:     ///
kpeter@617:     /// \warning The queue must not be empty.
deba@512:     Node processNextNode() {
deba@512:       Node node = queue.back();
deba@512:       queue.pop_back();
deba@512:       if ((*_node_order)[node] == -2) {
deba@512:         Arc arc = prepare(node);
deba@512:         Node source = _digraph->source(arc);
deba@512:         while ((*_node_order)[source] != -1) {
deba@512:           if ((*_node_order)[source] >= 0) {
deba@512:             arc = contract(source);
deba@512:           } else {
deba@512:             arc = prepare(source);
deba@512:           }
deba@512:           source = _digraph->source(arc);
deba@512:         }
deba@512:         finalize(arc);
deba@512:         level_stack.clear();
deba@512:       }
deba@512:       return node;
deba@512:     }
deba@512: 
deba@512:     /// \brief Returns the number of the nodes to be processed.
deba@512:     ///
kpeter@617:     /// Returns the number of the nodes to be processed in the priority
kpeter@617:     /// queue.
deba@512:     int queueSize() const {
deba@512:       return queue.size();
deba@512:     }
deba@512: 
deba@512:     /// \brief Returns \c false if there are nodes to be processed.
deba@512:     ///
deba@512:     /// Returns \c false if there are nodes to be processed.
deba@512:     bool emptyQueue() const {
deba@512:       return queue.empty();
deba@512:     }
deba@512: 
deba@512:     /// \brief Executes the algorithm.
deba@512:     ///
deba@512:     /// Executes the algorithm.
deba@512:     ///
deba@512:     /// \pre init() must be called and at least one node should be added
deba@512:     /// with addSource() before using this function.
deba@512:     ///
deba@512:     ///\note mca.start() is just a shortcut of the following code.
deba@512:     ///\code
deba@512:     ///while (!mca.emptyQueue()) {
deba@512:     ///  mca.processNextNode();
deba@512:     ///}
deba@512:     ///\endcode
deba@512:     void start() {
deba@512:       while (!emptyQueue()) {
deba@512:         processNextNode();
deba@512:       }
deba@512:     }
deba@512: 
deba@512:     /// \brief Runs %MinCostArborescence algorithm from node \c s.
deba@512:     ///
deba@512:     /// This method runs the %MinCostArborescence algorithm from
deba@512:     /// a root node \c s.
deba@512:     ///
deba@512:     /// \note mca.run(s) is just a shortcut of the following code.
deba@512:     /// \code
deba@512:     /// mca.init();
deba@512:     /// mca.addSource(s);
deba@512:     /// mca.start();
deba@512:     /// \endcode
kpeter@617:     void run(Node s) {
deba@512:       init();
kpeter@617:       addSource(s);
deba@512:       start();
deba@512:     }
deba@512: 
deba@512:     ///@}
deba@512: 
kpeter@617:     /// \name Query Functions
kpeter@617:     /// The result of the %MinCostArborescence algorithm can be obtained
kpeter@617:     /// using these functions.\n
kpeter@617:     /// Either run() or start() must be called before using them.
kpeter@617: 
kpeter@617:     /// @{
kpeter@617: 
kpeter@617:     /// \brief Returns the cost of the arborescence.
kpeter@617:     ///
kpeter@617:     /// Returns the cost of the arborescence.
kpeter@617:     Value arborescenceCost() const {
kpeter@617:       Value sum = 0;
kpeter@617:       for (ArcIt it(*_digraph); it != INVALID; ++it) {
kpeter@617:         if (arborescence(it)) {
kpeter@617:           sum += (*_cost)[it];
kpeter@617:         }
kpeter@617:       }
kpeter@617:       return sum;
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief Returns \c true if the arc is in the arborescence.
kpeter@617:     ///
kpeter@617:     /// Returns \c true if the given arc is in the arborescence.
kpeter@617:     /// \param arc An arc of the digraph.
kpeter@617:     /// \pre \ref run() must be called before using this function.
kpeter@617:     bool arborescence(Arc arc) const {
kpeter@617:       return (*_pred)[_digraph->target(arc)] == arc;
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief Returns a const reference to the arborescence map.
kpeter@617:     ///
kpeter@617:     /// Returns a const reference to the arborescence map.
kpeter@617:     /// \pre \ref run() must be called before using this function.
kpeter@617:     const ArborescenceMap& arborescenceMap() const {
kpeter@617:       return *_arborescence;
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief Returns the predecessor arc of the given node.
kpeter@617:     ///
kpeter@617:     /// Returns the predecessor arc of the given node.
kpeter@617:     /// \pre \ref run() must be called before using this function.
kpeter@617:     Arc pred(Node node) const {
kpeter@617:       return (*_pred)[node];
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief Returns a const reference to the pred map.
kpeter@617:     ///
kpeter@617:     /// Returns a const reference to the pred map.
kpeter@617:     /// \pre \ref run() must be called before using this function.
kpeter@617:     const PredMap& predMap() const {
kpeter@617:       return *_pred;
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief Indicates that a node is reachable from the sources.
kpeter@617:     ///
kpeter@617:     /// Indicates that a node is reachable from the sources.
kpeter@617:     bool reached(Node node) const {
kpeter@617:       return (*_node_order)[node] != -3;
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief Indicates that a node is processed.
kpeter@617:     ///
kpeter@617:     /// Indicates that a node is processed. The arborescence path exists
kpeter@617:     /// from the source to the given node.
kpeter@617:     bool processed(Node node) const {
kpeter@617:       return (*_node_order)[node] == -1;
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief Returns the number of the dual variables in basis.
kpeter@617:     ///
kpeter@617:     /// Returns the number of the dual variables in basis.
kpeter@617:     int dualNum() const {
kpeter@617:       return _dual_variables.size();
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief Returns the value of the dual solution.
kpeter@617:     ///
kpeter@617:     /// Returns the value of the dual solution. It should be
kpeter@617:     /// equal to the arborescence value.
kpeter@617:     Value dualValue() const {
kpeter@617:       Value sum = 0;
kpeter@617:       for (int i = 0; i < int(_dual_variables.size()); ++i) {
kpeter@617:         sum += _dual_variables[i].value;
kpeter@617:       }
kpeter@617:       return sum;
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief Returns the number of the nodes in the dual variable.
kpeter@617:     ///
kpeter@617:     /// Returns the number of the nodes in the dual variable.
kpeter@617:     int dualSize(int k) const {
kpeter@617:       return _dual_variables[k].end - _dual_variables[k].begin;
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief Returns the value of the dual variable.
kpeter@617:     ///
kpeter@617:     /// Returns the the value of the dual variable.
kpeter@617:     Value dualValue(int k) const {
kpeter@617:       return _dual_variables[k].value;
kpeter@617:     }
kpeter@617: 
kpeter@617:     /// \brief LEMON iterator for getting a dual variable.
kpeter@617:     ///
kpeter@617:     /// This class provides a common style LEMON iterator for getting a
kpeter@617:     /// dual variable of \ref MinCostArborescence algorithm.
kpeter@617:     /// It iterates over a subset of the nodes.
kpeter@617:     class DualIt {
kpeter@617:     public:
kpeter@617: 
kpeter@617:       /// \brief Constructor.
kpeter@617:       ///
kpeter@617:       /// Constructor for getting the nodeset of the dual variable
kpeter@617:       /// of \ref MinCostArborescence algorithm.
kpeter@617:       DualIt(const MinCostArborescence& algorithm, int variable)
kpeter@617:         : _algorithm(&algorithm)
kpeter@617:       {
kpeter@617:         _index = _algorithm->_dual_variables[variable].begin;
kpeter@617:         _last = _algorithm->_dual_variables[variable].end;
kpeter@617:       }
kpeter@617: 
kpeter@617:       /// \brief Conversion to \c Node.
kpeter@617:       ///
kpeter@617:       /// Conversion to \c Node.
kpeter@617:       operator Node() const {
kpeter@617:         return _algorithm->_dual_node_list[_index];
kpeter@617:       }
kpeter@617: 
kpeter@617:       /// \brief Increment operator.
kpeter@617:       ///
kpeter@617:       /// Increment operator.
kpeter@617:       DualIt& operator++() {
kpeter@617:         ++_index;
kpeter@617:         return *this;
kpeter@617:       }
kpeter@617: 
kpeter@617:       /// \brief Validity checking
kpeter@617:       ///
kpeter@617:       /// Checks whether the iterator is invalid.
kpeter@617:       bool operator==(Invalid) const {
kpeter@617:         return _index == _last;
kpeter@617:       }
kpeter@617: 
kpeter@617:       /// \brief Validity checking
kpeter@617:       ///
kpeter@617:       /// Checks whether the iterator is valid.
kpeter@617:       bool operator!=(Invalid) const {
kpeter@617:         return _index != _last;
kpeter@617:       }
kpeter@617: 
kpeter@617:     private:
kpeter@617:       const MinCostArborescence* _algorithm;
kpeter@617:       int _index, _last;
kpeter@617:     };
kpeter@617: 
kpeter@617:     /// @}
kpeter@617: 
deba@512:   };
deba@512: 
deba@512:   /// \ingroup spantree
deba@512:   ///
deba@512:   /// \brief Function type interface for MinCostArborescence algorithm.
deba@512:   ///
deba@512:   /// Function type interface for MinCostArborescence algorithm.
kpeter@617:   /// \param digraph The digraph the algorithm runs on.
kpeter@617:   /// \param cost An arc map storing the costs.
kpeter@617:   /// \param source The source node of the arborescence.
kpeter@617:   /// \retval arborescence An arc map with \c bool (or convertible) value
kpeter@617:   /// type that stores the arborescence.
kpeter@617:   /// \return The total cost of the arborescence.
deba@512:   ///
deba@512:   /// \sa MinCostArborescence
deba@512:   template <typename Digraph, typename CostMap, typename ArborescenceMap>
deba@512:   typename CostMap::Value minCostArborescence(const Digraph& digraph,
deba@512:                                               const CostMap& cost,
deba@512:                                               typename Digraph::Node source,
deba@512:                                               ArborescenceMap& arborescence) {
deba@512:     typename MinCostArborescence<Digraph, CostMap>
kpeter@617:       ::template SetArborescenceMap<ArborescenceMap>
deba@512:       ::Create mca(digraph, cost);
deba@512:     mca.arborescenceMap(arborescence);
deba@512:     mca.run(source);
kpeter@617:     return mca.arborescenceCost();
deba@512:   }
deba@512: 
deba@512: }
deba@512: 
deba@512: #endif