lemon/min_cost_arborescence.h
changeset 490 7f8560cb9d65
child 559 c5fd2d996909
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lemon/min_cost_arborescence.h	Tue Dec 02 23:33:47 2008 +0100
     1.3 @@ -0,0 +1,796 @@
     1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     1.5 + *
     1.6 + * This file is a part of LEMON, a generic C++ optimization library.
     1.7 + *
     1.8 + * Copyright (C) 2003-2008
     1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    1.11 + *
    1.12 + * Permission to use, modify and distribute this software is granted
    1.13 + * provided that this copyright notice appears in all copies. For
    1.14 + * precise terms see the accompanying LICENSE file.
    1.15 + *
    1.16 + * This software is provided "AS IS" with no warranty of any kind,
    1.17 + * express or implied, and with no claim as to its suitability for any
    1.18 + * purpose.
    1.19 + *
    1.20 + */
    1.21 +
    1.22 +#ifndef LEMON_MIN_COST_ARBORESCENCE_H
    1.23 +#define LEMON_MIN_COST_ARBORESCENCE_H
    1.24 +
    1.25 +///\ingroup spantree
    1.26 +///\file
    1.27 +///\brief Minimum Cost Arborescence algorithm.
    1.28 +
    1.29 +#include <vector>
    1.30 +
    1.31 +#include <lemon/list_graph.h>
    1.32 +#include <lemon/bin_heap.h>
    1.33 +#include <lemon/assert.h>
    1.34 +
    1.35 +namespace lemon {
    1.36 +
    1.37 +
    1.38 +  /// \brief Default traits class for MinCostArborescence class.
    1.39 +  ///
    1.40 +  /// Default traits class for MinCostArborescence class.
    1.41 +  /// \param _Digraph Digraph type.
    1.42 +  /// \param _CostMap Type of cost map.
    1.43 +  template <class _Digraph, class _CostMap>
    1.44 +  struct MinCostArborescenceDefaultTraits{
    1.45 +
    1.46 +    /// \brief The digraph type the algorithm runs on.
    1.47 +    typedef _Digraph Digraph;
    1.48 +
    1.49 +    /// \brief The type of the map that stores the arc costs.
    1.50 +    ///
    1.51 +    /// The type of the map that stores the arc costs.
    1.52 +    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
    1.53 +    typedef _CostMap CostMap;
    1.54 +
    1.55 +    /// \brief The value type of the costs.
    1.56 +    ///
    1.57 +    /// The value type of the costs.
    1.58 +    typedef typename CostMap::Value Value;
    1.59 +
    1.60 +    /// \brief The type of the map that stores which arcs are in the
    1.61 +    /// arborescence.
    1.62 +    ///
    1.63 +    /// The type of the map that stores which arcs are in the
    1.64 +    /// arborescence.  It must meet the \ref concepts::WriteMap
    1.65 +    /// "WriteMap" concept.  Initially it will be set to false on each
    1.66 +    /// arc. After it will set all arborescence arcs once.
    1.67 +    typedef typename Digraph::template ArcMap<bool> ArborescenceMap;
    1.68 +
    1.69 +    /// \brief Instantiates a ArborescenceMap.
    1.70 +    ///
    1.71 +    /// This function instantiates a \ref ArborescenceMap.
    1.72 +    /// \param digraph is the graph, to which we would like to
    1.73 +    /// calculate the ArborescenceMap.
    1.74 +    static ArborescenceMap *createArborescenceMap(const Digraph &digraph){
    1.75 +      return new ArborescenceMap(digraph);
    1.76 +    }
    1.77 +
    1.78 +    /// \brief The type of the PredMap
    1.79 +    ///
    1.80 +    /// The type of the PredMap. It is a node map with an arc value type.
    1.81 +    typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
    1.82 +
    1.83 +    /// \brief Instantiates a PredMap.
    1.84 +    ///
    1.85 +    /// This function instantiates a \ref PredMap.
    1.86 +    /// \param _digraph is the digraph, to which we would like to define the
    1.87 +    /// PredMap.
    1.88 +    static PredMap *createPredMap(const Digraph &digraph){
    1.89 +      return new PredMap(digraph);
    1.90 +    }
    1.91 +
    1.92 +  };
    1.93 +
    1.94 +  /// \ingroup spantree
    1.95 +  ///
    1.96 +  /// \brief %MinCostArborescence algorithm class.
    1.97 +  ///
    1.98 +  /// This class provides an efficient implementation of
    1.99 +  /// %MinCostArborescence algorithm. The arborescence is a tree
   1.100 +  /// which is directed from a given source node of the digraph. One or
   1.101 +  /// more sources should be given for the algorithm and it will calculate
   1.102 +  /// the minimum cost subgraph which are union of arborescences with the
   1.103 +  /// given sources and spans all the nodes which are reachable from the
   1.104 +  /// sources. The time complexity of the algorithm is \f$ O(n^2+e) \f$.
   1.105 +  ///
   1.106 +  /// The algorithm provides also an optimal dual solution, therefore
   1.107 +  /// the optimality of the solution can be checked.
   1.108 +  ///
   1.109 +  /// \param _Digraph The digraph type the algorithm runs on. The default value
   1.110 +  /// is \ref ListDigraph.
   1.111 +  /// \param _CostMap This read-only ArcMap determines the costs of the
   1.112 +  /// arcs. It is read once for each arc, so the map may involve in
   1.113 +  /// relatively time consuming process to compute the arc cost if
   1.114 +  /// it is necessary. The default map type is \ref
   1.115 +  /// concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
   1.116 +  /// \param _Traits Traits class to set various data types used
   1.117 +  /// by the algorithm. The default traits class is
   1.118 +  /// \ref MinCostArborescenceDefaultTraits
   1.119 +  /// "MinCostArborescenceDefaultTraits<_Digraph, _CostMap>".  See \ref
   1.120 +  /// MinCostArborescenceDefaultTraits for the documentation of a
   1.121 +  /// MinCostArborescence traits class.
   1.122 +  ///
   1.123 +  /// \author Balazs Dezso
   1.124 +#ifndef DOXYGEN
   1.125 +  template <typename _Digraph = ListDigraph,
   1.126 +            typename _CostMap = typename _Digraph::template ArcMap<int>,
   1.127 +            typename _Traits =
   1.128 +            MinCostArborescenceDefaultTraits<_Digraph, _CostMap> >
   1.129 +#else
   1.130 +  template <typename _Digraph, typename _CostMap, typedef _Traits>
   1.131 +#endif
   1.132 +  class MinCostArborescence {
   1.133 +  public:
   1.134 +
   1.135 +    /// The traits.
   1.136 +    typedef _Traits Traits;
   1.137 +    /// The type of the underlying digraph.
   1.138 +    typedef typename Traits::Digraph Digraph;
   1.139 +    /// The type of the map that stores the arc costs.
   1.140 +    typedef typename Traits::CostMap CostMap;
   1.141 +    ///The type of the costs of the arcs.
   1.142 +    typedef typename Traits::Value Value;
   1.143 +    ///The type of the predecessor map.
   1.144 +    typedef typename Traits::PredMap PredMap;
   1.145 +    ///The type of the map that stores which arcs are in the arborescence.
   1.146 +    typedef typename Traits::ArborescenceMap ArborescenceMap;
   1.147 +
   1.148 +    typedef MinCostArborescence Create;
   1.149 +
   1.150 +  private:
   1.151 +
   1.152 +    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
   1.153 +
   1.154 +    struct CostArc {
   1.155 +
   1.156 +      Arc arc;
   1.157 +      Value value;
   1.158 +
   1.159 +      CostArc() {}
   1.160 +      CostArc(Arc _arc, Value _value) : arc(_arc), value(_value) {}
   1.161 +
   1.162 +    };
   1.163 +
   1.164 +    const Digraph *_digraph;
   1.165 +    const CostMap *_cost;
   1.166 +
   1.167 +    PredMap *_pred;
   1.168 +    bool local_pred;
   1.169 +
   1.170 +    ArborescenceMap *_arborescence;
   1.171 +    bool local_arborescence;
   1.172 +
   1.173 +    typedef typename Digraph::template ArcMap<int> ArcOrder;
   1.174 +    ArcOrder *_arc_order;
   1.175 +
   1.176 +    typedef typename Digraph::template NodeMap<int> NodeOrder;
   1.177 +    NodeOrder *_node_order;
   1.178 +
   1.179 +    typedef typename Digraph::template NodeMap<CostArc> CostArcMap;
   1.180 +    CostArcMap *_cost_arcs;
   1.181 +
   1.182 +    struct StackLevel {
   1.183 +
   1.184 +      std::vector<CostArc> arcs;
   1.185 +      int node_level;
   1.186 +
   1.187 +    };
   1.188 +
   1.189 +    std::vector<StackLevel> level_stack;
   1.190 +    std::vector<Node> queue;
   1.191 +
   1.192 +    typedef std::vector<typename Digraph::Node> DualNodeList;
   1.193 +
   1.194 +    DualNodeList _dual_node_list;
   1.195 +
   1.196 +    struct DualVariable {
   1.197 +      int begin, end;
   1.198 +      Value value;
   1.199 +
   1.200 +      DualVariable(int _begin, int _end, Value _value)
   1.201 +        : begin(_begin), end(_end), value(_value) {}
   1.202 +
   1.203 +    };
   1.204 +
   1.205 +    typedef std::vector<DualVariable> DualVariables;
   1.206 +
   1.207 +    DualVariables _dual_variables;
   1.208 +
   1.209 +    typedef typename Digraph::template NodeMap<int> HeapCrossRef;
   1.210 +
   1.211 +    HeapCrossRef *_heap_cross_ref;
   1.212 +
   1.213 +    typedef BinHeap<int, HeapCrossRef> Heap;
   1.214 +
   1.215 +    Heap *_heap;
   1.216 +
   1.217 +  protected:
   1.218 +
   1.219 +    MinCostArborescence() {}
   1.220 +
   1.221 +  private:
   1.222 +
   1.223 +    void createStructures() {
   1.224 +      if (!_pred) {
   1.225 +        local_pred = true;
   1.226 +        _pred = Traits::createPredMap(*_digraph);
   1.227 +      }
   1.228 +      if (!_arborescence) {
   1.229 +        local_arborescence = true;
   1.230 +        _arborescence = Traits::createArborescenceMap(*_digraph);
   1.231 +      }
   1.232 +      if (!_arc_order) {
   1.233 +        _arc_order = new ArcOrder(*_digraph);
   1.234 +      }
   1.235 +      if (!_node_order) {
   1.236 +        _node_order = new NodeOrder(*_digraph);
   1.237 +      }
   1.238 +      if (!_cost_arcs) {
   1.239 +        _cost_arcs = new CostArcMap(*_digraph);
   1.240 +      }
   1.241 +      if (!_heap_cross_ref) {
   1.242 +        _heap_cross_ref = new HeapCrossRef(*_digraph, -1);
   1.243 +      }
   1.244 +      if (!_heap) {
   1.245 +        _heap = new Heap(*_heap_cross_ref);
   1.246 +      }
   1.247 +    }
   1.248 +
   1.249 +    void destroyStructures() {
   1.250 +      if (local_arborescence) {
   1.251 +        delete _arborescence;
   1.252 +      }
   1.253 +      if (local_pred) {
   1.254 +        delete _pred;
   1.255 +      }
   1.256 +      if (_arc_order) {
   1.257 +        delete _arc_order;
   1.258 +      }
   1.259 +      if (_node_order) {
   1.260 +        delete _node_order;
   1.261 +      }
   1.262 +      if (_cost_arcs) {
   1.263 +        delete _cost_arcs;
   1.264 +      }
   1.265 +      if (_heap) {
   1.266 +        delete _heap;
   1.267 +      }
   1.268 +      if (_heap_cross_ref) {
   1.269 +        delete _heap_cross_ref;
   1.270 +      }
   1.271 +    }
   1.272 +
   1.273 +    Arc prepare(Node node) {
   1.274 +      std::vector<Node> nodes;
   1.275 +      (*_node_order)[node] = _dual_node_list.size();
   1.276 +      StackLevel level;
   1.277 +      level.node_level = _dual_node_list.size();
   1.278 +      _dual_node_list.push_back(node);
   1.279 +      for (InArcIt it(*_digraph, node); it != INVALID; ++it) {
   1.280 +        Arc arc = it;
   1.281 +        Node source = _digraph->source(arc);
   1.282 +        Value value = (*_cost)[it];
   1.283 +        if (source == node || (*_node_order)[source] == -3) continue;
   1.284 +        if ((*_cost_arcs)[source].arc == INVALID) {
   1.285 +          (*_cost_arcs)[source].arc = arc;
   1.286 +          (*_cost_arcs)[source].value = value;
   1.287 +          nodes.push_back(source);
   1.288 +        } else {
   1.289 +          if ((*_cost_arcs)[source].value > value) {
   1.290 +            (*_cost_arcs)[source].arc = arc;
   1.291 +            (*_cost_arcs)[source].value = value;
   1.292 +          }
   1.293 +        }
   1.294 +      }
   1.295 +      CostArc minimum = (*_cost_arcs)[nodes[0]];
   1.296 +      for (int i = 1; i < int(nodes.size()); ++i) {
   1.297 +        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
   1.298 +          minimum = (*_cost_arcs)[nodes[i]];
   1.299 +        }
   1.300 +      }
   1.301 +      _arc_order->set(minimum.arc, _dual_variables.size());
   1.302 +      DualVariable var(_dual_node_list.size() - 1,
   1.303 +                       _dual_node_list.size(), minimum.value);
   1.304 +      _dual_variables.push_back(var);
   1.305 +      for (int i = 0; i < int(nodes.size()); ++i) {
   1.306 +        (*_cost_arcs)[nodes[i]].value -= minimum.value;
   1.307 +        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
   1.308 +        (*_cost_arcs)[nodes[i]].arc = INVALID;
   1.309 +      }
   1.310 +      level_stack.push_back(level);
   1.311 +      return minimum.arc;
   1.312 +    }
   1.313 +
   1.314 +    Arc contract(Node node) {
   1.315 +      int node_bottom = bottom(node);
   1.316 +      std::vector<Node> nodes;
   1.317 +      while (!level_stack.empty() &&
   1.318 +             level_stack.back().node_level >= node_bottom) {
   1.319 +        for (int i = 0; i < int(level_stack.back().arcs.size()); ++i) {
   1.320 +          Arc arc = level_stack.back().arcs[i].arc;
   1.321 +          Node source = _digraph->source(arc);
   1.322 +          Value value = level_stack.back().arcs[i].value;
   1.323 +          if ((*_node_order)[source] >= node_bottom) continue;
   1.324 +          if ((*_cost_arcs)[source].arc == INVALID) {
   1.325 +            (*_cost_arcs)[source].arc = arc;
   1.326 +            (*_cost_arcs)[source].value = value;
   1.327 +            nodes.push_back(source);
   1.328 +          } else {
   1.329 +            if ((*_cost_arcs)[source].value > value) {
   1.330 +              (*_cost_arcs)[source].arc = arc;
   1.331 +              (*_cost_arcs)[source].value = value;
   1.332 +            }
   1.333 +          }
   1.334 +        }
   1.335 +        level_stack.pop_back();
   1.336 +      }
   1.337 +      CostArc minimum = (*_cost_arcs)[nodes[0]];
   1.338 +      for (int i = 1; i < int(nodes.size()); ++i) {
   1.339 +        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
   1.340 +          minimum = (*_cost_arcs)[nodes[i]];
   1.341 +        }
   1.342 +      }
   1.343 +      _arc_order->set(minimum.arc, _dual_variables.size());
   1.344 +      DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
   1.345 +      _dual_variables.push_back(var);
   1.346 +      StackLevel level;
   1.347 +      level.node_level = node_bottom;
   1.348 +      for (int i = 0; i < int(nodes.size()); ++i) {
   1.349 +        (*_cost_arcs)[nodes[i]].value -= minimum.value;
   1.350 +        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
   1.351 +        (*_cost_arcs)[nodes[i]].arc = INVALID;
   1.352 +      }
   1.353 +      level_stack.push_back(level);
   1.354 +      return minimum.arc;
   1.355 +    }
   1.356 +
   1.357 +    int bottom(Node node) {
   1.358 +      int k = level_stack.size() - 1;
   1.359 +      while (level_stack[k].node_level > (*_node_order)[node]) {
   1.360 +        --k;
   1.361 +      }
   1.362 +      return level_stack[k].node_level;
   1.363 +    }
   1.364 +
   1.365 +    void finalize(Arc arc) {
   1.366 +      Node node = _digraph->target(arc);
   1.367 +      _heap->push(node, (*_arc_order)[arc]);
   1.368 +      _pred->set(node, arc);
   1.369 +      while (!_heap->empty()) {
   1.370 +        Node source = _heap->top();
   1.371 +        _heap->pop();
   1.372 +        _node_order->set(source, -1);
   1.373 +        for (OutArcIt it(*_digraph, source); it != INVALID; ++it) {
   1.374 +          if ((*_arc_order)[it] < 0) continue;
   1.375 +          Node target = _digraph->target(it);
   1.376 +          switch(_heap->state(target)) {
   1.377 +          case Heap::PRE_HEAP:
   1.378 +            _heap->push(target, (*_arc_order)[it]);
   1.379 +            _pred->set(target, it);
   1.380 +            break;
   1.381 +          case Heap::IN_HEAP:
   1.382 +            if ((*_arc_order)[it] < (*_heap)[target]) {
   1.383 +              _heap->decrease(target, (*_arc_order)[it]);
   1.384 +              _pred->set(target, it);
   1.385 +            }
   1.386 +            break;
   1.387 +          case Heap::POST_HEAP:
   1.388 +            break;
   1.389 +          }
   1.390 +        }
   1.391 +        _arborescence->set((*_pred)[source], true);
   1.392 +      }
   1.393 +    }
   1.394 +
   1.395 +
   1.396 +  public:
   1.397 +
   1.398 +    /// \name Named template parameters
   1.399 +
   1.400 +    /// @{
   1.401 +
   1.402 +    template <class T>
   1.403 +    struct DefArborescenceMapTraits : public Traits {
   1.404 +      typedef T ArborescenceMap;
   1.405 +      static ArborescenceMap *createArborescenceMap(const Digraph &)
   1.406 +      {
   1.407 +        LEMON_ASSERT(false, "ArborescenceMap is not initialized");
   1.408 +        return 0; // ignore warnings
   1.409 +      }
   1.410 +    };
   1.411 +
   1.412 +    /// \brief \ref named-templ-param "Named parameter" for
   1.413 +    /// setting ArborescenceMap type
   1.414 +    ///
   1.415 +    /// \ref named-templ-param "Named parameter" for setting
   1.416 +    /// ArborescenceMap type
   1.417 +    template <class T>
   1.418 +    struct DefArborescenceMap
   1.419 +      : public MinCostArborescence<Digraph, CostMap,
   1.420 +                                   DefArborescenceMapTraits<T> > {
   1.421 +    };
   1.422 +
   1.423 +    template <class T>
   1.424 +    struct DefPredMapTraits : public Traits {
   1.425 +      typedef T PredMap;
   1.426 +      static PredMap *createPredMap(const Digraph &)
   1.427 +      {
   1.428 +        LEMON_ASSERT(false, "PredMap is not initialized");
   1.429 +      }
   1.430 +    };
   1.431 +
   1.432 +    /// \brief \ref named-templ-param "Named parameter" for
   1.433 +    /// setting PredMap type
   1.434 +    ///
   1.435 +    /// \ref named-templ-param "Named parameter" for setting
   1.436 +    /// PredMap type
   1.437 +    template <class T>
   1.438 +    struct DefPredMap
   1.439 +      : public MinCostArborescence<Digraph, CostMap, DefPredMapTraits<T> > {
   1.440 +    };
   1.441 +
   1.442 +    /// @}
   1.443 +
   1.444 +    /// \brief Constructor.
   1.445 +    ///
   1.446 +    /// \param _digraph The digraph the algorithm will run on.
   1.447 +    /// \param _cost The cost map used by the algorithm.
   1.448 +    MinCostArborescence(const Digraph& digraph, const CostMap& cost)
   1.449 +      : _digraph(&digraph), _cost(&cost), _pred(0), local_pred(false),
   1.450 +        _arborescence(0), local_arborescence(false),
   1.451 +        _arc_order(0), _node_order(0), _cost_arcs(0),
   1.452 +        _heap_cross_ref(0), _heap(0) {}
   1.453 +
   1.454 +    /// \brief Destructor.
   1.455 +    ~MinCostArborescence() {
   1.456 +      destroyStructures();
   1.457 +    }
   1.458 +
   1.459 +    /// \brief Sets the arborescence map.
   1.460 +    ///
   1.461 +    /// Sets the arborescence map.
   1.462 +    /// \return \c (*this)
   1.463 +    MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
   1.464 +      if (local_arborescence) {
   1.465 +        delete _arborescence;
   1.466 +      }
   1.467 +      local_arborescence = false;
   1.468 +      _arborescence = &m;
   1.469 +      return *this;
   1.470 +    }
   1.471 +
   1.472 +    /// \brief Sets the arborescence map.
   1.473 +    ///
   1.474 +    /// Sets the arborescence map.
   1.475 +    /// \return \c (*this)
   1.476 +    MinCostArborescence& predMap(PredMap& m) {
   1.477 +      if (local_pred) {
   1.478 +        delete _pred;
   1.479 +      }
   1.480 +      local_pred = false;
   1.481 +      _pred = &m;
   1.482 +      return *this;
   1.483 +    }
   1.484 +
   1.485 +    /// \name Query Functions
   1.486 +    /// The result of the %MinCostArborescence algorithm can be obtained
   1.487 +    /// using these functions.\n
   1.488 +    /// Before the use of these functions,
   1.489 +    /// either run() or start() must be called.
   1.490 +
   1.491 +    /// @{
   1.492 +
   1.493 +    /// \brief Returns a reference to the arborescence map.
   1.494 +    ///
   1.495 +    /// Returns a reference to the arborescence map.
   1.496 +    const ArborescenceMap& arborescenceMap() const {
   1.497 +      return *_arborescence;
   1.498 +    }
   1.499 +
   1.500 +    /// \brief Returns true if the arc is in the arborescence.
   1.501 +    ///
   1.502 +    /// Returns true if the arc is in the arborescence.
   1.503 +    /// \param arc The arc of the digraph.
   1.504 +    /// \pre \ref run() must be called before using this function.
   1.505 +    bool arborescence(Arc arc) const {
   1.506 +      return (*_pred)[_digraph->target(arc)] == arc;
   1.507 +    }
   1.508 +
   1.509 +    /// \brief Returns a reference to the pred map.
   1.510 +    ///
   1.511 +    /// Returns a reference to the pred map.
   1.512 +    const PredMap& predMap() const {
   1.513 +      return *_pred;
   1.514 +    }
   1.515 +
   1.516 +    /// \brief Returns the predecessor arc of the given node.
   1.517 +    ///
   1.518 +    /// Returns the predecessor arc of the given node.
   1.519 +    Arc pred(Node node) const {
   1.520 +      return (*_pred)[node];
   1.521 +    }
   1.522 +
   1.523 +    /// \brief Returns the cost of the arborescence.
   1.524 +    ///
   1.525 +    /// Returns the cost of the arborescence.
   1.526 +    Value arborescenceValue() const {
   1.527 +      Value sum = 0;
   1.528 +      for (ArcIt it(*_digraph); it != INVALID; ++it) {
   1.529 +        if (arborescence(it)) {
   1.530 +          sum += (*_cost)[it];
   1.531 +        }
   1.532 +      }
   1.533 +      return sum;
   1.534 +    }
   1.535 +
   1.536 +    /// \brief Indicates that a node is reachable from the sources.
   1.537 +    ///
   1.538 +    /// Indicates that a node is reachable from the sources.
   1.539 +    bool reached(Node node) const {
   1.540 +      return (*_node_order)[node] != -3;
   1.541 +    }
   1.542 +
   1.543 +    /// \brief Indicates that a node is processed.
   1.544 +    ///
   1.545 +    /// Indicates that a node is processed. The arborescence path exists
   1.546 +    /// from the source to the given node.
   1.547 +    bool processed(Node node) const {
   1.548 +      return (*_node_order)[node] == -1;
   1.549 +    }
   1.550 +
   1.551 +    /// \brief Returns the number of the dual variables in basis.
   1.552 +    ///
   1.553 +    /// Returns the number of the dual variables in basis.
   1.554 +    int dualNum() const {
   1.555 +      return _dual_variables.size();
   1.556 +    }
   1.557 +
   1.558 +    /// \brief Returns the value of the dual solution.
   1.559 +    ///
   1.560 +    /// Returns the value of the dual solution. It should be
   1.561 +    /// equal to the arborescence value.
   1.562 +    Value dualValue() const {
   1.563 +      Value sum = 0;
   1.564 +      for (int i = 0; i < int(_dual_variables.size()); ++i) {
   1.565 +        sum += _dual_variables[i].value;
   1.566 +      }
   1.567 +      return sum;
   1.568 +    }
   1.569 +
   1.570 +    /// \brief Returns the number of the nodes in the dual variable.
   1.571 +    ///
   1.572 +    /// Returns the number of the nodes in the dual variable.
   1.573 +    int dualSize(int k) const {
   1.574 +      return _dual_variables[k].end - _dual_variables[k].begin;
   1.575 +    }
   1.576 +
   1.577 +    /// \brief Returns the value of the dual variable.
   1.578 +    ///
   1.579 +    /// Returns the the value of the dual variable.
   1.580 +    const Value& dualValue(int k) const {
   1.581 +      return _dual_variables[k].value;
   1.582 +    }
   1.583 +
   1.584 +    /// \brief Lemon iterator for get a dual variable.
   1.585 +    ///
   1.586 +    /// Lemon iterator for get a dual variable. This class provides
   1.587 +    /// a common style lemon iterator which gives back a subset of
   1.588 +    /// the nodes.
   1.589 +    class DualIt {
   1.590 +    public:
   1.591 +
   1.592 +      /// \brief Constructor.
   1.593 +      ///
   1.594 +      /// Constructor for get the nodeset of the variable.
   1.595 +      DualIt(const MinCostArborescence& algorithm, int variable)
   1.596 +        : _algorithm(&algorithm)
   1.597 +      {
   1.598 +        _index = _algorithm->_dual_variables[variable].begin;
   1.599 +        _last = _algorithm->_dual_variables[variable].end;
   1.600 +      }
   1.601 +
   1.602 +      /// \brief Conversion to node.
   1.603 +      ///
   1.604 +      /// Conversion to node.
   1.605 +      operator Node() const {
   1.606 +        return _algorithm->_dual_node_list[_index];
   1.607 +      }
   1.608 +
   1.609 +      /// \brief Increment operator.
   1.610 +      ///
   1.611 +      /// Increment operator.
   1.612 +      DualIt& operator++() {
   1.613 +        ++_index;
   1.614 +        return *this;
   1.615 +      }
   1.616 +
   1.617 +      /// \brief Validity checking
   1.618 +      ///
   1.619 +      /// Checks whether the iterator is invalid.
   1.620 +      bool operator==(Invalid) const {
   1.621 +        return _index == _last;
   1.622 +      }
   1.623 +
   1.624 +      /// \brief Validity checking
   1.625 +      ///
   1.626 +      /// Checks whether the iterator is valid.
   1.627 +      bool operator!=(Invalid) const {
   1.628 +        return _index != _last;
   1.629 +      }
   1.630 +
   1.631 +    private:
   1.632 +      const MinCostArborescence* _algorithm;
   1.633 +      int _index, _last;
   1.634 +    };
   1.635 +
   1.636 +    /// @}
   1.637 +
   1.638 +    /// \name Execution control
   1.639 +    /// The simplest way to execute the algorithm is to use
   1.640 +    /// one of the member functions called \c run(...). \n
   1.641 +    /// If you need more control on the execution,
   1.642 +    /// first you must call \ref init(), then you can add several
   1.643 +    /// source nodes with \ref addSource().
   1.644 +    /// Finally \ref start() will perform the arborescence
   1.645 +    /// computation.
   1.646 +
   1.647 +    ///@{
   1.648 +
   1.649 +    /// \brief Initializes the internal data structures.
   1.650 +    ///
   1.651 +    /// Initializes the internal data structures.
   1.652 +    ///
   1.653 +    void init() {
   1.654 +      createStructures();
   1.655 +      _heap->clear();
   1.656 +      for (NodeIt it(*_digraph); it != INVALID; ++it) {
   1.657 +        (*_cost_arcs)[it].arc = INVALID;
   1.658 +        _node_order->set(it, -3);
   1.659 +        _heap_cross_ref->set(it, Heap::PRE_HEAP);
   1.660 +        _pred->set(it, INVALID);
   1.661 +      }
   1.662 +      for (ArcIt it(*_digraph); it != INVALID; ++it) {
   1.663 +        _arborescence->set(it, false);
   1.664 +        _arc_order->set(it, -1);
   1.665 +      }
   1.666 +      _dual_node_list.clear();
   1.667 +      _dual_variables.clear();
   1.668 +    }
   1.669 +
   1.670 +    /// \brief Adds a new source node.
   1.671 +    ///
   1.672 +    /// Adds a new source node to the algorithm.
   1.673 +    void addSource(Node source) {
   1.674 +      std::vector<Node> nodes;
   1.675 +      nodes.push_back(source);
   1.676 +      while (!nodes.empty()) {
   1.677 +        Node node = nodes.back();
   1.678 +        nodes.pop_back();
   1.679 +        for (OutArcIt it(*_digraph, node); it != INVALID; ++it) {
   1.680 +          Node target = _digraph->target(it);
   1.681 +          if ((*_node_order)[target] == -3) {
   1.682 +            (*_node_order)[target] = -2;
   1.683 +            nodes.push_back(target);
   1.684 +            queue.push_back(target);
   1.685 +          }
   1.686 +        }
   1.687 +      }
   1.688 +      (*_node_order)[source] = -1;
   1.689 +    }
   1.690 +
   1.691 +    /// \brief Processes the next node in the priority queue.
   1.692 +    ///
   1.693 +    /// Processes the next node in the priority queue.
   1.694 +    ///
   1.695 +    /// \return The processed node.
   1.696 +    ///
   1.697 +    /// \warning The queue must not be empty!
   1.698 +    Node processNextNode() {
   1.699 +      Node node = queue.back();
   1.700 +      queue.pop_back();
   1.701 +      if ((*_node_order)[node] == -2) {
   1.702 +        Arc arc = prepare(node);
   1.703 +        Node source = _digraph->source(arc);
   1.704 +        while ((*_node_order)[source] != -1) {
   1.705 +          if ((*_node_order)[source] >= 0) {
   1.706 +            arc = contract(source);
   1.707 +          } else {
   1.708 +            arc = prepare(source);
   1.709 +          }
   1.710 +          source = _digraph->source(arc);
   1.711 +        }
   1.712 +        finalize(arc);
   1.713 +        level_stack.clear();
   1.714 +      }
   1.715 +      return node;
   1.716 +    }
   1.717 +
   1.718 +    /// \brief Returns the number of the nodes to be processed.
   1.719 +    ///
   1.720 +    /// Returns the number of the nodes to be processed.
   1.721 +    int queueSize() const {
   1.722 +      return queue.size();
   1.723 +    }
   1.724 +
   1.725 +    /// \brief Returns \c false if there are nodes to be processed.
   1.726 +    ///
   1.727 +    /// Returns \c false if there are nodes to be processed.
   1.728 +    bool emptyQueue() const {
   1.729 +      return queue.empty();
   1.730 +    }
   1.731 +
   1.732 +    /// \brief Executes the algorithm.
   1.733 +    ///
   1.734 +    /// Executes the algorithm.
   1.735 +    ///
   1.736 +    /// \pre init() must be called and at least one node should be added
   1.737 +    /// with addSource() before using this function.
   1.738 +    ///
   1.739 +    ///\note mca.start() is just a shortcut of the following code.
   1.740 +    ///\code
   1.741 +    ///while (!mca.emptyQueue()) {
   1.742 +    ///  mca.processNextNode();
   1.743 +    ///}
   1.744 +    ///\endcode
   1.745 +    void start() {
   1.746 +      while (!emptyQueue()) {
   1.747 +        processNextNode();
   1.748 +      }
   1.749 +    }
   1.750 +
   1.751 +    /// \brief Runs %MinCostArborescence algorithm from node \c s.
   1.752 +    ///
   1.753 +    /// This method runs the %MinCostArborescence algorithm from
   1.754 +    /// a root node \c s.
   1.755 +    ///
   1.756 +    /// \note mca.run(s) is just a shortcut of the following code.
   1.757 +    /// \code
   1.758 +    /// mca.init();
   1.759 +    /// mca.addSource(s);
   1.760 +    /// mca.start();
   1.761 +    /// \endcode
   1.762 +    void run(Node node) {
   1.763 +      init();
   1.764 +      addSource(node);
   1.765 +      start();
   1.766 +    }
   1.767 +
   1.768 +    ///@}
   1.769 +
   1.770 +  };
   1.771 +
   1.772 +  /// \ingroup spantree
   1.773 +  ///
   1.774 +  /// \brief Function type interface for MinCostArborescence algorithm.
   1.775 +  ///
   1.776 +  /// Function type interface for MinCostArborescence algorithm.
   1.777 +  /// \param digraph The Digraph that the algorithm runs on.
   1.778 +  /// \param cost The CostMap of the arcs.
   1.779 +  /// \param source The source of the arborescence.
   1.780 +  /// \retval arborescence The bool ArcMap which stores the arborescence.
   1.781 +  /// \return The cost of the arborescence.
   1.782 +  ///
   1.783 +  /// \sa MinCostArborescence
   1.784 +  template <typename Digraph, typename CostMap, typename ArborescenceMap>
   1.785 +  typename CostMap::Value minCostArborescence(const Digraph& digraph,
   1.786 +                                              const CostMap& cost,
   1.787 +                                              typename Digraph::Node source,
   1.788 +                                              ArborescenceMap& arborescence) {
   1.789 +    typename MinCostArborescence<Digraph, CostMap>
   1.790 +      ::template DefArborescenceMap<ArborescenceMap>
   1.791 +      ::Create mca(digraph, cost);
   1.792 +    mca.arborescenceMap(arborescence);
   1.793 +    mca.run(source);
   1.794 +    return mca.arborescenceValue();
   1.795 +  }
   1.796 +
   1.797 +}
   1.798 +
   1.799 +#endif