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