Merge
authorAlpar Juttner <alpar@cs.elte.hu>
Mon, 23 Feb 2009 12:26:21 +0000
changeset 52888bd39ef7d98
parent 527 c458e02723b1
parent 501 7f8560cb9d65
child 530 1c5d6e47921f
Merge
lemon/Makefile.am
test/CMakeLists.txt
test/Makefile.am
     1.1 --- a/lemon/Makefile.am	Mon Feb 23 11:52:45 2009 +0000
     1.2 +++ b/lemon/Makefile.am	Mon Feb 23 12:26:21 2009 +0000
     1.3 @@ -83,6 +83,7 @@
     1.4  	lemon/maps.h \
     1.5  	lemon/math.h \
     1.6  	lemon/max_matching.h \
     1.7 +	lemon/min_cost_arborescence.h \
     1.8  	lemon/nauty_reader.h \
     1.9  	lemon/path.h \
    1.10  	lemon/preflow.h \
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/lemon/min_cost_arborescence.h	Mon Feb 23 12:26:21 2009 +0000
     2.3 @@ -0,0 +1,796 @@
     2.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2.5 + *
     2.6 + * This file is a part of LEMON, a generic C++ optimization library.
     2.7 + *
     2.8 + * Copyright (C) 2003-2008
     2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    2.11 + *
    2.12 + * Permission to use, modify and distribute this software is granted
    2.13 + * provided that this copyright notice appears in all copies. For
    2.14 + * precise terms see the accompanying LICENSE file.
    2.15 + *
    2.16 + * This software is provided "AS IS" with no warranty of any kind,
    2.17 + * express or implied, and with no claim as to its suitability for any
    2.18 + * purpose.
    2.19 + *
    2.20 + */
    2.21 +
    2.22 +#ifndef LEMON_MIN_COST_ARBORESCENCE_H
    2.23 +#define LEMON_MIN_COST_ARBORESCENCE_H
    2.24 +
    2.25 +///\ingroup spantree
    2.26 +///\file
    2.27 +///\brief Minimum Cost Arborescence algorithm.
    2.28 +
    2.29 +#include <vector>
    2.30 +
    2.31 +#include <lemon/list_graph.h>
    2.32 +#include <lemon/bin_heap.h>
    2.33 +#include <lemon/assert.h>
    2.34 +
    2.35 +namespace lemon {
    2.36 +
    2.37 +
    2.38 +  /// \brief Default traits class for MinCostArborescence class.
    2.39 +  ///
    2.40 +  /// Default traits class for MinCostArborescence class.
    2.41 +  /// \param _Digraph Digraph type.
    2.42 +  /// \param _CostMap Type of cost map.
    2.43 +  template <class _Digraph, class _CostMap>
    2.44 +  struct MinCostArborescenceDefaultTraits{
    2.45 +
    2.46 +    /// \brief The digraph type the algorithm runs on.
    2.47 +    typedef _Digraph Digraph;
    2.48 +
    2.49 +    /// \brief The type of the map that stores the arc costs.
    2.50 +    ///
    2.51 +    /// The type of the map that stores the arc costs.
    2.52 +    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
    2.53 +    typedef _CostMap CostMap;
    2.54 +
    2.55 +    /// \brief The value type of the costs.
    2.56 +    ///
    2.57 +    /// The value type of the costs.
    2.58 +    typedef typename CostMap::Value Value;
    2.59 +
    2.60 +    /// \brief The type of the map that stores which arcs are in the
    2.61 +    /// arborescence.
    2.62 +    ///
    2.63 +    /// The type of the map that stores which arcs are in the
    2.64 +    /// arborescence.  It must meet the \ref concepts::WriteMap
    2.65 +    /// "WriteMap" concept.  Initially it will be set to false on each
    2.66 +    /// arc. After it will set all arborescence arcs once.
    2.67 +    typedef typename Digraph::template ArcMap<bool> ArborescenceMap;
    2.68 +
    2.69 +    /// \brief Instantiates a ArborescenceMap.
    2.70 +    ///
    2.71 +    /// This function instantiates a \ref ArborescenceMap.
    2.72 +    /// \param digraph is the graph, to which we would like to
    2.73 +    /// calculate the ArborescenceMap.
    2.74 +    static ArborescenceMap *createArborescenceMap(const Digraph &digraph){
    2.75 +      return new ArborescenceMap(digraph);
    2.76 +    }
    2.77 +
    2.78 +    /// \brief The type of the PredMap
    2.79 +    ///
    2.80 +    /// The type of the PredMap. It is a node map with an arc value type.
    2.81 +    typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
    2.82 +
    2.83 +    /// \brief Instantiates a PredMap.
    2.84 +    ///
    2.85 +    /// This function instantiates a \ref PredMap.
    2.86 +    /// \param _digraph is the digraph, to which we would like to define the
    2.87 +    /// PredMap.
    2.88 +    static PredMap *createPredMap(const Digraph &digraph){
    2.89 +      return new PredMap(digraph);
    2.90 +    }
    2.91 +
    2.92 +  };
    2.93 +
    2.94 +  /// \ingroup spantree
    2.95 +  ///
    2.96 +  /// \brief %MinCostArborescence algorithm class.
    2.97 +  ///
    2.98 +  /// This class provides an efficient implementation of
    2.99 +  /// %MinCostArborescence algorithm. The arborescence is a tree
   2.100 +  /// which is directed from a given source node of the digraph. One or
   2.101 +  /// more sources should be given for the algorithm and it will calculate
   2.102 +  /// the minimum cost subgraph which are union of arborescences with the
   2.103 +  /// given sources and spans all the nodes which are reachable from the
   2.104 +  /// sources. The time complexity of the algorithm is \f$ O(n^2+e) \f$.
   2.105 +  ///
   2.106 +  /// The algorithm provides also an optimal dual solution, therefore
   2.107 +  /// the optimality of the solution can be checked.
   2.108 +  ///
   2.109 +  /// \param _Digraph The digraph type the algorithm runs on. The default value
   2.110 +  /// is \ref ListDigraph.
   2.111 +  /// \param _CostMap This read-only ArcMap determines the costs of the
   2.112 +  /// arcs. It is read once for each arc, so the map may involve in
   2.113 +  /// relatively time consuming process to compute the arc cost if
   2.114 +  /// it is necessary. The default map type is \ref
   2.115 +  /// concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
   2.116 +  /// \param _Traits Traits class to set various data types used
   2.117 +  /// by the algorithm. The default traits class is
   2.118 +  /// \ref MinCostArborescenceDefaultTraits
   2.119 +  /// "MinCostArborescenceDefaultTraits<_Digraph, _CostMap>".  See \ref
   2.120 +  /// MinCostArborescenceDefaultTraits for the documentation of a
   2.121 +  /// MinCostArborescence traits class.
   2.122 +  ///
   2.123 +  /// \author Balazs Dezso
   2.124 +#ifndef DOXYGEN
   2.125 +  template <typename _Digraph = ListDigraph,
   2.126 +            typename _CostMap = typename _Digraph::template ArcMap<int>,
   2.127 +            typename _Traits =
   2.128 +            MinCostArborescenceDefaultTraits<_Digraph, _CostMap> >
   2.129 +#else
   2.130 +  template <typename _Digraph, typename _CostMap, typedef _Traits>
   2.131 +#endif
   2.132 +  class MinCostArborescence {
   2.133 +  public:
   2.134 +
   2.135 +    /// The traits.
   2.136 +    typedef _Traits Traits;
   2.137 +    /// The type of the underlying digraph.
   2.138 +    typedef typename Traits::Digraph Digraph;
   2.139 +    /// The type of the map that stores the arc costs.
   2.140 +    typedef typename Traits::CostMap CostMap;
   2.141 +    ///The type of the costs of the arcs.
   2.142 +    typedef typename Traits::Value Value;
   2.143 +    ///The type of the predecessor map.
   2.144 +    typedef typename Traits::PredMap PredMap;
   2.145 +    ///The type of the map that stores which arcs are in the arborescence.
   2.146 +    typedef typename Traits::ArborescenceMap ArborescenceMap;
   2.147 +
   2.148 +    typedef MinCostArborescence Create;
   2.149 +
   2.150 +  private:
   2.151 +
   2.152 +    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
   2.153 +
   2.154 +    struct CostArc {
   2.155 +
   2.156 +      Arc arc;
   2.157 +      Value value;
   2.158 +
   2.159 +      CostArc() {}
   2.160 +      CostArc(Arc _arc, Value _value) : arc(_arc), value(_value) {}
   2.161 +
   2.162 +    };
   2.163 +
   2.164 +    const Digraph *_digraph;
   2.165 +    const CostMap *_cost;
   2.166 +
   2.167 +    PredMap *_pred;
   2.168 +    bool local_pred;
   2.169 +
   2.170 +    ArborescenceMap *_arborescence;
   2.171 +    bool local_arborescence;
   2.172 +
   2.173 +    typedef typename Digraph::template ArcMap<int> ArcOrder;
   2.174 +    ArcOrder *_arc_order;
   2.175 +
   2.176 +    typedef typename Digraph::template NodeMap<int> NodeOrder;
   2.177 +    NodeOrder *_node_order;
   2.178 +
   2.179 +    typedef typename Digraph::template NodeMap<CostArc> CostArcMap;
   2.180 +    CostArcMap *_cost_arcs;
   2.181 +
   2.182 +    struct StackLevel {
   2.183 +
   2.184 +      std::vector<CostArc> arcs;
   2.185 +      int node_level;
   2.186 +
   2.187 +    };
   2.188 +
   2.189 +    std::vector<StackLevel> level_stack;
   2.190 +    std::vector<Node> queue;
   2.191 +
   2.192 +    typedef std::vector<typename Digraph::Node> DualNodeList;
   2.193 +
   2.194 +    DualNodeList _dual_node_list;
   2.195 +
   2.196 +    struct DualVariable {
   2.197 +      int begin, end;
   2.198 +      Value value;
   2.199 +
   2.200 +      DualVariable(int _begin, int _end, Value _value)
   2.201 +        : begin(_begin), end(_end), value(_value) {}
   2.202 +
   2.203 +    };
   2.204 +
   2.205 +    typedef std::vector<DualVariable> DualVariables;
   2.206 +
   2.207 +    DualVariables _dual_variables;
   2.208 +
   2.209 +    typedef typename Digraph::template NodeMap<int> HeapCrossRef;
   2.210 +
   2.211 +    HeapCrossRef *_heap_cross_ref;
   2.212 +
   2.213 +    typedef BinHeap<int, HeapCrossRef> Heap;
   2.214 +
   2.215 +    Heap *_heap;
   2.216 +
   2.217 +  protected:
   2.218 +
   2.219 +    MinCostArborescence() {}
   2.220 +
   2.221 +  private:
   2.222 +
   2.223 +    void createStructures() {
   2.224 +      if (!_pred) {
   2.225 +        local_pred = true;
   2.226 +        _pred = Traits::createPredMap(*_digraph);
   2.227 +      }
   2.228 +      if (!_arborescence) {
   2.229 +        local_arborescence = true;
   2.230 +        _arborescence = Traits::createArborescenceMap(*_digraph);
   2.231 +      }
   2.232 +      if (!_arc_order) {
   2.233 +        _arc_order = new ArcOrder(*_digraph);
   2.234 +      }
   2.235 +      if (!_node_order) {
   2.236 +        _node_order = new NodeOrder(*_digraph);
   2.237 +      }
   2.238 +      if (!_cost_arcs) {
   2.239 +        _cost_arcs = new CostArcMap(*_digraph);
   2.240 +      }
   2.241 +      if (!_heap_cross_ref) {
   2.242 +        _heap_cross_ref = new HeapCrossRef(*_digraph, -1);
   2.243 +      }
   2.244 +      if (!_heap) {
   2.245 +        _heap = new Heap(*_heap_cross_ref);
   2.246 +      }
   2.247 +    }
   2.248 +
   2.249 +    void destroyStructures() {
   2.250 +      if (local_arborescence) {
   2.251 +        delete _arborescence;
   2.252 +      }
   2.253 +      if (local_pred) {
   2.254 +        delete _pred;
   2.255 +      }
   2.256 +      if (_arc_order) {
   2.257 +        delete _arc_order;
   2.258 +      }
   2.259 +      if (_node_order) {
   2.260 +        delete _node_order;
   2.261 +      }
   2.262 +      if (_cost_arcs) {
   2.263 +        delete _cost_arcs;
   2.264 +      }
   2.265 +      if (_heap) {
   2.266 +        delete _heap;
   2.267 +      }
   2.268 +      if (_heap_cross_ref) {
   2.269 +        delete _heap_cross_ref;
   2.270 +      }
   2.271 +    }
   2.272 +
   2.273 +    Arc prepare(Node node) {
   2.274 +      std::vector<Node> nodes;
   2.275 +      (*_node_order)[node] = _dual_node_list.size();
   2.276 +      StackLevel level;
   2.277 +      level.node_level = _dual_node_list.size();
   2.278 +      _dual_node_list.push_back(node);
   2.279 +      for (InArcIt it(*_digraph, node); it != INVALID; ++it) {
   2.280 +        Arc arc = it;
   2.281 +        Node source = _digraph->source(arc);
   2.282 +        Value value = (*_cost)[it];
   2.283 +        if (source == node || (*_node_order)[source] == -3) continue;
   2.284 +        if ((*_cost_arcs)[source].arc == INVALID) {
   2.285 +          (*_cost_arcs)[source].arc = arc;
   2.286 +          (*_cost_arcs)[source].value = value;
   2.287 +          nodes.push_back(source);
   2.288 +        } else {
   2.289 +          if ((*_cost_arcs)[source].value > value) {
   2.290 +            (*_cost_arcs)[source].arc = arc;
   2.291 +            (*_cost_arcs)[source].value = value;
   2.292 +          }
   2.293 +        }
   2.294 +      }
   2.295 +      CostArc minimum = (*_cost_arcs)[nodes[0]];
   2.296 +      for (int i = 1; i < int(nodes.size()); ++i) {
   2.297 +        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
   2.298 +          minimum = (*_cost_arcs)[nodes[i]];
   2.299 +        }
   2.300 +      }
   2.301 +      _arc_order->set(minimum.arc, _dual_variables.size());
   2.302 +      DualVariable var(_dual_node_list.size() - 1,
   2.303 +                       _dual_node_list.size(), minimum.value);
   2.304 +      _dual_variables.push_back(var);
   2.305 +      for (int i = 0; i < int(nodes.size()); ++i) {
   2.306 +        (*_cost_arcs)[nodes[i]].value -= minimum.value;
   2.307 +        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
   2.308 +        (*_cost_arcs)[nodes[i]].arc = INVALID;
   2.309 +      }
   2.310 +      level_stack.push_back(level);
   2.311 +      return minimum.arc;
   2.312 +    }
   2.313 +
   2.314 +    Arc contract(Node node) {
   2.315 +      int node_bottom = bottom(node);
   2.316 +      std::vector<Node> nodes;
   2.317 +      while (!level_stack.empty() &&
   2.318 +             level_stack.back().node_level >= node_bottom) {
   2.319 +        for (int i = 0; i < int(level_stack.back().arcs.size()); ++i) {
   2.320 +          Arc arc = level_stack.back().arcs[i].arc;
   2.321 +          Node source = _digraph->source(arc);
   2.322 +          Value value = level_stack.back().arcs[i].value;
   2.323 +          if ((*_node_order)[source] >= node_bottom) continue;
   2.324 +          if ((*_cost_arcs)[source].arc == INVALID) {
   2.325 +            (*_cost_arcs)[source].arc = arc;
   2.326 +            (*_cost_arcs)[source].value = value;
   2.327 +            nodes.push_back(source);
   2.328 +          } else {
   2.329 +            if ((*_cost_arcs)[source].value > value) {
   2.330 +              (*_cost_arcs)[source].arc = arc;
   2.331 +              (*_cost_arcs)[source].value = value;
   2.332 +            }
   2.333 +          }
   2.334 +        }
   2.335 +        level_stack.pop_back();
   2.336 +      }
   2.337 +      CostArc minimum = (*_cost_arcs)[nodes[0]];
   2.338 +      for (int i = 1; i < int(nodes.size()); ++i) {
   2.339 +        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
   2.340 +          minimum = (*_cost_arcs)[nodes[i]];
   2.341 +        }
   2.342 +      }
   2.343 +      _arc_order->set(minimum.arc, _dual_variables.size());
   2.344 +      DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
   2.345 +      _dual_variables.push_back(var);
   2.346 +      StackLevel level;
   2.347 +      level.node_level = node_bottom;
   2.348 +      for (int i = 0; i < int(nodes.size()); ++i) {
   2.349 +        (*_cost_arcs)[nodes[i]].value -= minimum.value;
   2.350 +        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
   2.351 +        (*_cost_arcs)[nodes[i]].arc = INVALID;
   2.352 +      }
   2.353 +      level_stack.push_back(level);
   2.354 +      return minimum.arc;
   2.355 +    }
   2.356 +
   2.357 +    int bottom(Node node) {
   2.358 +      int k = level_stack.size() - 1;
   2.359 +      while (level_stack[k].node_level > (*_node_order)[node]) {
   2.360 +        --k;
   2.361 +      }
   2.362 +      return level_stack[k].node_level;
   2.363 +    }
   2.364 +
   2.365 +    void finalize(Arc arc) {
   2.366 +      Node node = _digraph->target(arc);
   2.367 +      _heap->push(node, (*_arc_order)[arc]);
   2.368 +      _pred->set(node, arc);
   2.369 +      while (!_heap->empty()) {
   2.370 +        Node source = _heap->top();
   2.371 +        _heap->pop();
   2.372 +        _node_order->set(source, -1);
   2.373 +        for (OutArcIt it(*_digraph, source); it != INVALID; ++it) {
   2.374 +          if ((*_arc_order)[it] < 0) continue;
   2.375 +          Node target = _digraph->target(it);
   2.376 +          switch(_heap->state(target)) {
   2.377 +          case Heap::PRE_HEAP:
   2.378 +            _heap->push(target, (*_arc_order)[it]);
   2.379 +            _pred->set(target, it);
   2.380 +            break;
   2.381 +          case Heap::IN_HEAP:
   2.382 +            if ((*_arc_order)[it] < (*_heap)[target]) {
   2.383 +              _heap->decrease(target, (*_arc_order)[it]);
   2.384 +              _pred->set(target, it);
   2.385 +            }
   2.386 +            break;
   2.387 +          case Heap::POST_HEAP:
   2.388 +            break;
   2.389 +          }
   2.390 +        }
   2.391 +        _arborescence->set((*_pred)[source], true);
   2.392 +      }
   2.393 +    }
   2.394 +
   2.395 +
   2.396 +  public:
   2.397 +
   2.398 +    /// \name Named template parameters
   2.399 +
   2.400 +    /// @{
   2.401 +
   2.402 +    template <class T>
   2.403 +    struct DefArborescenceMapTraits : public Traits {
   2.404 +      typedef T ArborescenceMap;
   2.405 +      static ArborescenceMap *createArborescenceMap(const Digraph &)
   2.406 +      {
   2.407 +        LEMON_ASSERT(false, "ArborescenceMap is not initialized");
   2.408 +        return 0; // ignore warnings
   2.409 +      }
   2.410 +    };
   2.411 +
   2.412 +    /// \brief \ref named-templ-param "Named parameter" for
   2.413 +    /// setting ArborescenceMap type
   2.414 +    ///
   2.415 +    /// \ref named-templ-param "Named parameter" for setting
   2.416 +    /// ArborescenceMap type
   2.417 +    template <class T>
   2.418 +    struct DefArborescenceMap
   2.419 +      : public MinCostArborescence<Digraph, CostMap,
   2.420 +                                   DefArborescenceMapTraits<T> > {
   2.421 +    };
   2.422 +
   2.423 +    template <class T>
   2.424 +    struct DefPredMapTraits : public Traits {
   2.425 +      typedef T PredMap;
   2.426 +      static PredMap *createPredMap(const Digraph &)
   2.427 +      {
   2.428 +        LEMON_ASSERT(false, "PredMap is not initialized");
   2.429 +      }
   2.430 +    };
   2.431 +
   2.432 +    /// \brief \ref named-templ-param "Named parameter" for
   2.433 +    /// setting PredMap type
   2.434 +    ///
   2.435 +    /// \ref named-templ-param "Named parameter" for setting
   2.436 +    /// PredMap type
   2.437 +    template <class T>
   2.438 +    struct DefPredMap
   2.439 +      : public MinCostArborescence<Digraph, CostMap, DefPredMapTraits<T> > {
   2.440 +    };
   2.441 +
   2.442 +    /// @}
   2.443 +
   2.444 +    /// \brief Constructor.
   2.445 +    ///
   2.446 +    /// \param _digraph The digraph the algorithm will run on.
   2.447 +    /// \param _cost The cost map used by the algorithm.
   2.448 +    MinCostArborescence(const Digraph& digraph, const CostMap& cost)
   2.449 +      : _digraph(&digraph), _cost(&cost), _pred(0), local_pred(false),
   2.450 +        _arborescence(0), local_arborescence(false),
   2.451 +        _arc_order(0), _node_order(0), _cost_arcs(0),
   2.452 +        _heap_cross_ref(0), _heap(0) {}
   2.453 +
   2.454 +    /// \brief Destructor.
   2.455 +    ~MinCostArborescence() {
   2.456 +      destroyStructures();
   2.457 +    }
   2.458 +
   2.459 +    /// \brief Sets the arborescence map.
   2.460 +    ///
   2.461 +    /// Sets the arborescence map.
   2.462 +    /// \return \c (*this)
   2.463 +    MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
   2.464 +      if (local_arborescence) {
   2.465 +        delete _arborescence;
   2.466 +      }
   2.467 +      local_arborescence = false;
   2.468 +      _arborescence = &m;
   2.469 +      return *this;
   2.470 +    }
   2.471 +
   2.472 +    /// \brief Sets the arborescence map.
   2.473 +    ///
   2.474 +    /// Sets the arborescence map.
   2.475 +    /// \return \c (*this)
   2.476 +    MinCostArborescence& predMap(PredMap& m) {
   2.477 +      if (local_pred) {
   2.478 +        delete _pred;
   2.479 +      }
   2.480 +      local_pred = false;
   2.481 +      _pred = &m;
   2.482 +      return *this;
   2.483 +    }
   2.484 +
   2.485 +    /// \name Query Functions
   2.486 +    /// The result of the %MinCostArborescence algorithm can be obtained
   2.487 +    /// using these functions.\n
   2.488 +    /// Before the use of these functions,
   2.489 +    /// either run() or start() must be called.
   2.490 +
   2.491 +    /// @{
   2.492 +
   2.493 +    /// \brief Returns a reference to the arborescence map.
   2.494 +    ///
   2.495 +    /// Returns a reference to the arborescence map.
   2.496 +    const ArborescenceMap& arborescenceMap() const {
   2.497 +      return *_arborescence;
   2.498 +    }
   2.499 +
   2.500 +    /// \brief Returns true if the arc is in the arborescence.
   2.501 +    ///
   2.502 +    /// Returns true if the arc is in the arborescence.
   2.503 +    /// \param arc The arc of the digraph.
   2.504 +    /// \pre \ref run() must be called before using this function.
   2.505 +    bool arborescence(Arc arc) const {
   2.506 +      return (*_pred)[_digraph->target(arc)] == arc;
   2.507 +    }
   2.508 +
   2.509 +    /// \brief Returns a reference to the pred map.
   2.510 +    ///
   2.511 +    /// Returns a reference to the pred map.
   2.512 +    const PredMap& predMap() const {
   2.513 +      return *_pred;
   2.514 +    }
   2.515 +
   2.516 +    /// \brief Returns the predecessor arc of the given node.
   2.517 +    ///
   2.518 +    /// Returns the predecessor arc of the given node.
   2.519 +    Arc pred(Node node) const {
   2.520 +      return (*_pred)[node];
   2.521 +    }
   2.522 +
   2.523 +    /// \brief Returns the cost of the arborescence.
   2.524 +    ///
   2.525 +    /// Returns the cost of the arborescence.
   2.526 +    Value arborescenceValue() const {
   2.527 +      Value sum = 0;
   2.528 +      for (ArcIt it(*_digraph); it != INVALID; ++it) {
   2.529 +        if (arborescence(it)) {
   2.530 +          sum += (*_cost)[it];
   2.531 +        }
   2.532 +      }
   2.533 +      return sum;
   2.534 +    }
   2.535 +
   2.536 +    /// \brief Indicates that a node is reachable from the sources.
   2.537 +    ///
   2.538 +    /// Indicates that a node is reachable from the sources.
   2.539 +    bool reached(Node node) const {
   2.540 +      return (*_node_order)[node] != -3;
   2.541 +    }
   2.542 +
   2.543 +    /// \brief Indicates that a node is processed.
   2.544 +    ///
   2.545 +    /// Indicates that a node is processed. The arborescence path exists
   2.546 +    /// from the source to the given node.
   2.547 +    bool processed(Node node) const {
   2.548 +      return (*_node_order)[node] == -1;
   2.549 +    }
   2.550 +
   2.551 +    /// \brief Returns the number of the dual variables in basis.
   2.552 +    ///
   2.553 +    /// Returns the number of the dual variables in basis.
   2.554 +    int dualNum() const {
   2.555 +      return _dual_variables.size();
   2.556 +    }
   2.557 +
   2.558 +    /// \brief Returns the value of the dual solution.
   2.559 +    ///
   2.560 +    /// Returns the value of the dual solution. It should be
   2.561 +    /// equal to the arborescence value.
   2.562 +    Value dualValue() const {
   2.563 +      Value sum = 0;
   2.564 +      for (int i = 0; i < int(_dual_variables.size()); ++i) {
   2.565 +        sum += _dual_variables[i].value;
   2.566 +      }
   2.567 +      return sum;
   2.568 +    }
   2.569 +
   2.570 +    /// \brief Returns the number of the nodes in the dual variable.
   2.571 +    ///
   2.572 +    /// Returns the number of the nodes in the dual variable.
   2.573 +    int dualSize(int k) const {
   2.574 +      return _dual_variables[k].end - _dual_variables[k].begin;
   2.575 +    }
   2.576 +
   2.577 +    /// \brief Returns the value of the dual variable.
   2.578 +    ///
   2.579 +    /// Returns the the value of the dual variable.
   2.580 +    const Value& dualValue(int k) const {
   2.581 +      return _dual_variables[k].value;
   2.582 +    }
   2.583 +
   2.584 +    /// \brief Lemon iterator for get a dual variable.
   2.585 +    ///
   2.586 +    /// Lemon iterator for get a dual variable. This class provides
   2.587 +    /// a common style lemon iterator which gives back a subset of
   2.588 +    /// the nodes.
   2.589 +    class DualIt {
   2.590 +    public:
   2.591 +
   2.592 +      /// \brief Constructor.
   2.593 +      ///
   2.594 +      /// Constructor for get the nodeset of the variable.
   2.595 +      DualIt(const MinCostArborescence& algorithm, int variable)
   2.596 +        : _algorithm(&algorithm)
   2.597 +      {
   2.598 +        _index = _algorithm->_dual_variables[variable].begin;
   2.599 +        _last = _algorithm->_dual_variables[variable].end;
   2.600 +      }
   2.601 +
   2.602 +      /// \brief Conversion to node.
   2.603 +      ///
   2.604 +      /// Conversion to node.
   2.605 +      operator Node() const {
   2.606 +        return _algorithm->_dual_node_list[_index];
   2.607 +      }
   2.608 +
   2.609 +      /// \brief Increment operator.
   2.610 +      ///
   2.611 +      /// Increment operator.
   2.612 +      DualIt& operator++() {
   2.613 +        ++_index;
   2.614 +        return *this;
   2.615 +      }
   2.616 +
   2.617 +      /// \brief Validity checking
   2.618 +      ///
   2.619 +      /// Checks whether the iterator is invalid.
   2.620 +      bool operator==(Invalid) const {
   2.621 +        return _index == _last;
   2.622 +      }
   2.623 +
   2.624 +      /// \brief Validity checking
   2.625 +      ///
   2.626 +      /// Checks whether the iterator is valid.
   2.627 +      bool operator!=(Invalid) const {
   2.628 +        return _index != _last;
   2.629 +      }
   2.630 +
   2.631 +    private:
   2.632 +      const MinCostArborescence* _algorithm;
   2.633 +      int _index, _last;
   2.634 +    };
   2.635 +
   2.636 +    /// @}
   2.637 +
   2.638 +    /// \name Execution control
   2.639 +    /// The simplest way to execute the algorithm is to use
   2.640 +    /// one of the member functions called \c run(...). \n
   2.641 +    /// If you need more control on the execution,
   2.642 +    /// first you must call \ref init(), then you can add several
   2.643 +    /// source nodes with \ref addSource().
   2.644 +    /// Finally \ref start() will perform the arborescence
   2.645 +    /// computation.
   2.646 +
   2.647 +    ///@{
   2.648 +
   2.649 +    /// \brief Initializes the internal data structures.
   2.650 +    ///
   2.651 +    /// Initializes the internal data structures.
   2.652 +    ///
   2.653 +    void init() {
   2.654 +      createStructures();
   2.655 +      _heap->clear();
   2.656 +      for (NodeIt it(*_digraph); it != INVALID; ++it) {
   2.657 +        (*_cost_arcs)[it].arc = INVALID;
   2.658 +        _node_order->set(it, -3);
   2.659 +        _heap_cross_ref->set(it, Heap::PRE_HEAP);
   2.660 +        _pred->set(it, INVALID);
   2.661 +      }
   2.662 +      for (ArcIt it(*_digraph); it != INVALID; ++it) {
   2.663 +        _arborescence->set(it, false);
   2.664 +        _arc_order->set(it, -1);
   2.665 +      }
   2.666 +      _dual_node_list.clear();
   2.667 +      _dual_variables.clear();
   2.668 +    }
   2.669 +
   2.670 +    /// \brief Adds a new source node.
   2.671 +    ///
   2.672 +    /// Adds a new source node to the algorithm.
   2.673 +    void addSource(Node source) {
   2.674 +      std::vector<Node> nodes;
   2.675 +      nodes.push_back(source);
   2.676 +      while (!nodes.empty()) {
   2.677 +        Node node = nodes.back();
   2.678 +        nodes.pop_back();
   2.679 +        for (OutArcIt it(*_digraph, node); it != INVALID; ++it) {
   2.680 +          Node target = _digraph->target(it);
   2.681 +          if ((*_node_order)[target] == -3) {
   2.682 +            (*_node_order)[target] = -2;
   2.683 +            nodes.push_back(target);
   2.684 +            queue.push_back(target);
   2.685 +          }
   2.686 +        }
   2.687 +      }
   2.688 +      (*_node_order)[source] = -1;
   2.689 +    }
   2.690 +
   2.691 +    /// \brief Processes the next node in the priority queue.
   2.692 +    ///
   2.693 +    /// Processes the next node in the priority queue.
   2.694 +    ///
   2.695 +    /// \return The processed node.
   2.696 +    ///
   2.697 +    /// \warning The queue must not be empty!
   2.698 +    Node processNextNode() {
   2.699 +      Node node = queue.back();
   2.700 +      queue.pop_back();
   2.701 +      if ((*_node_order)[node] == -2) {
   2.702 +        Arc arc = prepare(node);
   2.703 +        Node source = _digraph->source(arc);
   2.704 +        while ((*_node_order)[source] != -1) {
   2.705 +          if ((*_node_order)[source] >= 0) {
   2.706 +            arc = contract(source);
   2.707 +          } else {
   2.708 +            arc = prepare(source);
   2.709 +          }
   2.710 +          source = _digraph->source(arc);
   2.711 +        }
   2.712 +        finalize(arc);
   2.713 +        level_stack.clear();
   2.714 +      }
   2.715 +      return node;
   2.716 +    }
   2.717 +
   2.718 +    /// \brief Returns the number of the nodes to be processed.
   2.719 +    ///
   2.720 +    /// Returns the number of the nodes to be processed.
   2.721 +    int queueSize() const {
   2.722 +      return queue.size();
   2.723 +    }
   2.724 +
   2.725 +    /// \brief Returns \c false if there are nodes to be processed.
   2.726 +    ///
   2.727 +    /// Returns \c false if there are nodes to be processed.
   2.728 +    bool emptyQueue() const {
   2.729 +      return queue.empty();
   2.730 +    }
   2.731 +
   2.732 +    /// \brief Executes the algorithm.
   2.733 +    ///
   2.734 +    /// Executes the algorithm.
   2.735 +    ///
   2.736 +    /// \pre init() must be called and at least one node should be added
   2.737 +    /// with addSource() before using this function.
   2.738 +    ///
   2.739 +    ///\note mca.start() is just a shortcut of the following code.
   2.740 +    ///\code
   2.741 +    ///while (!mca.emptyQueue()) {
   2.742 +    ///  mca.processNextNode();
   2.743 +    ///}
   2.744 +    ///\endcode
   2.745 +    void start() {
   2.746 +      while (!emptyQueue()) {
   2.747 +        processNextNode();
   2.748 +      }
   2.749 +    }
   2.750 +
   2.751 +    /// \brief Runs %MinCostArborescence algorithm from node \c s.
   2.752 +    ///
   2.753 +    /// This method runs the %MinCostArborescence algorithm from
   2.754 +    /// a root node \c s.
   2.755 +    ///
   2.756 +    /// \note mca.run(s) is just a shortcut of the following code.
   2.757 +    /// \code
   2.758 +    /// mca.init();
   2.759 +    /// mca.addSource(s);
   2.760 +    /// mca.start();
   2.761 +    /// \endcode
   2.762 +    void run(Node node) {
   2.763 +      init();
   2.764 +      addSource(node);
   2.765 +      start();
   2.766 +    }
   2.767 +
   2.768 +    ///@}
   2.769 +
   2.770 +  };
   2.771 +
   2.772 +  /// \ingroup spantree
   2.773 +  ///
   2.774 +  /// \brief Function type interface for MinCostArborescence algorithm.
   2.775 +  ///
   2.776 +  /// Function type interface for MinCostArborescence algorithm.
   2.777 +  /// \param digraph The Digraph that the algorithm runs on.
   2.778 +  /// \param cost The CostMap of the arcs.
   2.779 +  /// \param source The source of the arborescence.
   2.780 +  /// \retval arborescence The bool ArcMap which stores the arborescence.
   2.781 +  /// \return The cost of the arborescence.
   2.782 +  ///
   2.783 +  /// \sa MinCostArborescence
   2.784 +  template <typename Digraph, typename CostMap, typename ArborescenceMap>
   2.785 +  typename CostMap::Value minCostArborescence(const Digraph& digraph,
   2.786 +                                              const CostMap& cost,
   2.787 +                                              typename Digraph::Node source,
   2.788 +                                              ArborescenceMap& arborescence) {
   2.789 +    typename MinCostArborescence<Digraph, CostMap>
   2.790 +      ::template DefArborescenceMap<ArborescenceMap>
   2.791 +      ::Create mca(digraph, cost);
   2.792 +    mca.arborescenceMap(arborescence);
   2.793 +    mca.run(source);
   2.794 +    return mca.arborescenceValue();
   2.795 +  }
   2.796 +
   2.797 +}
   2.798 +
   2.799 +#endif
     3.1 --- a/test/CMakeLists.txt	Mon Feb 23 11:52:45 2009 +0000
     3.2 +++ b/test/CMakeLists.txt	Mon Feb 23 12:26:21 2009 +0000
     3.3 @@ -29,6 +29,7 @@
     3.4    kruskal_test
     3.5    maps_test
     3.6    max_matching_test
     3.7 +  min_cost_arborescence_test
     3.8    path_test
     3.9    preflow_test
    3.10    radix_sort_test
     4.1 --- a/test/Makefile.am	Mon Feb 23 11:52:45 2009 +0000
     4.2 +++ b/test/Makefile.am	Mon Feb 23 12:26:21 2009 +0000
     4.3 @@ -25,6 +25,7 @@
     4.4  	test/kruskal_test \
     4.5  	test/maps_test \
     4.6  	test/max_matching_test \
     4.7 +	test/min_cost_arborescence_test \
     4.8  	test/path_test \
     4.9  	test/preflow_test \
    4.10  	test/radix_sort_test \
    4.11 @@ -66,6 +67,7 @@
    4.12  test_maps_test_SOURCES = test/maps_test.cc
    4.13  test_mip_test_SOURCES = test/mip_test.cc
    4.14  test_max_matching_test_SOURCES = test/max_matching_test.cc
    4.15 +test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
    4.16  test_path_test_SOURCES = test/path_test.cc
    4.17  test_preflow_test_SOURCES = test/preflow_test.cc
    4.18  test_radix_sort_test_SOURCES = test/radix_sort_test.cc
     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/test/min_cost_arborescence_test.cc	Mon Feb 23 12:26:21 2009 +0000
     5.3 @@ -0,0 +1,146 @@
     5.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     5.5 + *
     5.6 + * This file is a part of LEMON, a generic C++ optimization library.
     5.7 + *
     5.8 + * Copyright (C) 2003-2008
     5.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    5.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    5.11 + *
    5.12 + * Permission to use, modify and distribute this software is granted
    5.13 + * provided that this copyright notice appears in all copies. For
    5.14 + * precise terms see the accompanying LICENSE file.
    5.15 + *
    5.16 + * This software is provided "AS IS" with no warranty of any kind,
    5.17 + * express or implied, and with no claim as to its suitability for any
    5.18 + * purpose.
    5.19 + *
    5.20 + */
    5.21 +
    5.22 +#include <iostream>
    5.23 +#include <set>
    5.24 +#include <vector>
    5.25 +#include <iterator>
    5.26 +
    5.27 +#include <lemon/smart_graph.h>
    5.28 +#include <lemon/min_cost_arborescence.h>
    5.29 +#include <lemon/lgf_reader.h>
    5.30 +
    5.31 +#include "test_tools.h"
    5.32 +
    5.33 +using namespace lemon;
    5.34 +using namespace std;
    5.35 +
    5.36 +const char test_lgf[] =
    5.37 +  "@nodes\n"
    5.38 +  "label\n"
    5.39 +  "0\n"
    5.40 +  "1\n"
    5.41 +  "2\n"
    5.42 +  "3\n"
    5.43 +  "4\n"
    5.44 +  "5\n"
    5.45 +  "6\n"
    5.46 +  "7\n"
    5.47 +  "8\n"
    5.48 +  "9\n"
    5.49 +  "@arcs\n"
    5.50 +  "     label  cost\n"
    5.51 +  "1 8  0      107\n"
    5.52 +  "0 3  1      70\n"
    5.53 +  "2 1  2      46\n"
    5.54 +  "4 1  3      28\n"
    5.55 +  "4 4  4      91\n"
    5.56 +  "3 9  5      76\n"
    5.57 +  "9 8  6      61\n"
    5.58 +  "8 1  7      39\n"
    5.59 +  "9 8  8      74\n"
    5.60 +  "8 0  9      39\n"
    5.61 +  "4 3  10     45\n"
    5.62 +  "2 2  11     34\n"
    5.63 +  "0 1  12     100\n"
    5.64 +  "6 3  13     95\n"
    5.65 +  "4 1  14     22\n"
    5.66 +  "1 1  15     31\n"
    5.67 +  "7 2  16     51\n"
    5.68 +  "2 6  17     29\n"
    5.69 +  "8 3  18     115\n"
    5.70 +  "6 9  19     32\n"
    5.71 +  "1 1  20     60\n"
    5.72 +  "0 3  21     40\n"
    5.73 +  "@attributes\n"
    5.74 +  "source 0\n";
    5.75 +
    5.76 +int main() {
    5.77 +  typedef SmartDigraph Digraph;
    5.78 +  DIGRAPH_TYPEDEFS(Digraph);
    5.79 +
    5.80 +  typedef Digraph::ArcMap<double> CostMap;
    5.81 +
    5.82 +  Digraph digraph;
    5.83 +  CostMap cost(digraph);
    5.84 +  Node source;
    5.85 +
    5.86 +  std::istringstream is(test_lgf);
    5.87 +  digraphReader(digraph, is).
    5.88 +    arcMap("cost", cost).
    5.89 +    node("source", source).run();
    5.90 +
    5.91 +  MinCostArborescence<Digraph, CostMap> mca(digraph, cost);
    5.92 +  mca.run(source);
    5.93 +
    5.94 +  vector<pair<double, set<Node> > > dualSolution(mca.dualNum());
    5.95 +
    5.96 +  for (int i = 0; i < mca.dualNum(); ++i) {
    5.97 +    dualSolution[i].first = mca.dualValue(i);
    5.98 +    for (MinCostArborescence<Digraph, CostMap>::DualIt it(mca, i);
    5.99 +         it != INVALID; ++it) {
   5.100 +      dualSolution[i].second.insert(it);
   5.101 +    }
   5.102 +  }
   5.103 +
   5.104 +  for (ArcIt it(digraph); it != INVALID; ++it) {
   5.105 +    if (mca.reached(digraph.source(it))) {
   5.106 +      double sum = 0.0;
   5.107 +      for (int i = 0; i < int(dualSolution.size()); ++i) {
   5.108 +        if (dualSolution[i].second.find(digraph.target(it))
   5.109 +            != dualSolution[i].second.end() &&
   5.110 +            dualSolution[i].second.find(digraph.source(it))
   5.111 +            == dualSolution[i].second.end()) {
   5.112 +          sum += dualSolution[i].first;
   5.113 +        }
   5.114 +      }
   5.115 +      if (mca.arborescence(it)) {
   5.116 +        check(sum == cost[it], "INVALID DUAL");
   5.117 +      }
   5.118 +      check(sum <= cost[it], "INVALID DUAL");
   5.119 +    }
   5.120 +  }
   5.121 +
   5.122 +
   5.123 +  check(mca.dualValue() == mca.arborescenceValue(), "INVALID DUAL");
   5.124 +
   5.125 +  check(mca.reached(source), "INVALID ARBORESCENCE");
   5.126 +  for (ArcIt a(digraph); a != INVALID; ++a) {
   5.127 +    check(!mca.reached(digraph.source(a)) ||
   5.128 +          mca.reached(digraph.target(a)), "INVALID ARBORESCENCE");
   5.129 +  }
   5.130 +
   5.131 +  for (NodeIt n(digraph); n != INVALID; ++n) {
   5.132 +    if (!mca.reached(n)) continue;
   5.133 +    int cnt = 0;
   5.134 +    for (InArcIt a(digraph, n); a != INVALID; ++a) {
   5.135 +      if (mca.arborescence(a)) {
   5.136 +        check(mca.pred(n) == a, "INVALID ARBORESCENCE");
   5.137 +        ++cnt;
   5.138 +      }
   5.139 +    }
   5.140 +    check((n == source ? cnt == 0 : cnt == 1), "INVALID ARBORESCENCE");
   5.141 +  }
   5.142 +
   5.143 +  Digraph::ArcMap<bool> arborescence(digraph);
   5.144 +  check(mca.arborescenceValue() ==
   5.145 +        minCostArborescence(digraph, cost, source, arborescence),
   5.146 +        "WRONG FUNCTION");
   5.147 +
   5.148 +  return 0;
   5.149 +}