1.1 --- a/lemon/Makefile.am Tue Dec 02 15:33:22 2008 +0000
1.2 +++ b/lemon/Makefile.am Tue Dec 02 23:33:47 2008 +0100
1.3 @@ -44,6 +44,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 Tue Dec 02 23:33:47 2008 +0100
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 Tue Dec 02 15:33:22 2008 +0000
3.2 +++ b/test/CMakeLists.txt Tue Dec 02 23:33:47 2008 +0100
3.3 @@ -18,6 +18,7 @@
3.4 kruskal_test
3.5 maps_test
3.6 max_matching_test
3.7 + min_cost_arborescence_test
3.8 random_test
3.9 path_test
3.10 time_measure_test
4.1 --- a/test/Makefile.am Tue Dec 02 15:33:22 2008 +0000
4.2 +++ b/test/Makefile.am Tue Dec 02 23:33:47 2008 +0100
4.3 @@ -25,6 +25,7 @@
4.4 test/hao_orlin_test \
4.5 test/maps_test \
4.6 test/max_matching_test \
4.7 + test/min_cost_arborescence_test \
4.8 test/random_test \
4.9 test/path_test \
4.10 test/preflow_test \
4.11 @@ -54,6 +55,7 @@
4.12 test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
4.13 test_maps_test_SOURCES = test/maps_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_suurballe_test_SOURCES = test/suurballe_test.cc
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/test/min_cost_arborescence_test.cc Tue Dec 02 23:33:47 2008 +0100
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 +}