1.1 --- a/lemon/min_cost_arborescence.h Thu Mar 30 15:34:56 2006 +0000
1.2 +++ b/lemon/min_cost_arborescence.h Fri Mar 31 11:10:44 2006 +0000
1.3 @@ -26,6 +26,7 @@
1.4 #include <vector>
1.5
1.6 #include <lemon/list_graph.h>
1.7 +#include <lemon/bin_heap.h>
1.8
1.9 namespace lemon {
1.10
1.11 @@ -56,11 +57,9 @@
1.12 /// in the arborescence.
1.13 ///
1.14 /// The type of the map that stores which edges are in the arborescence.
1.15 - /// It must meet the \ref concept::ReadWriteMap "ReadWriteMap" concept.
1.16 - /// Initially it will be setted to false on each edge. The algorithm
1.17 - /// may set each value one time to true and maybe after it to false again.
1.18 - /// Therefore you cannot use maps like BackInserteBoolMap with this
1.19 - /// algorithm.
1.20 + /// It must meet the \ref concept::WritedMap "WriteMap" concept.
1.21 + /// Initially it will be setted to false on each edge. After it
1.22 + /// will set all arborescence edges once.
1.23 typedef typename Graph::template EdgeMap<bool> ArborescenceMap;
1.24
1.25 /// \brief Instantiates a ArborescenceMap.
1.26 @@ -72,6 +71,20 @@
1.27 return new ArborescenceMap(_graph);
1.28 }
1.29
1.30 + /// \brief The type of the PredMap
1.31 + ///
1.32 + /// The type of the PredMap. It is a node map with an edge value type.
1.33 + typedef typename Graph::template NodeMap<typename Graph::Edge> PredMap;
1.34 +
1.35 + /// \brief Instantiates a PredMap.
1.36 + ///
1.37 + /// This function instantiates a \ref PredMap.
1.38 + /// \param _graph is the graph, to which we would like to define the
1.39 + /// PredMap.
1.40 + static PredMap *createPredMap(const Graph &_graph){
1.41 + return new PredMap(_graph);
1.42 + }
1.43 +
1.44 };
1.45
1.46 /// \ingroup spantree
1.47 @@ -86,6 +99,9 @@
1.48 /// given sources and spans all the nodes which are reachable from the
1.49 /// sources. The time complexity of the algorithm is O(n^2 + e).
1.50 ///
1.51 + /// The algorithm provides also an optimal dual solution to arborescence
1.52 + /// that way the optimality of the algorithm can be proofed easily.
1.53 + ///
1.54 /// \param _Graph The graph type the algorithm runs on. The default value
1.55 /// is \ref ListGraph. The value of _Graph is not used directly by
1.56 /// MinCostArborescence, it is only passed to
1.57 @@ -135,6 +151,8 @@
1.58 typedef typename Traits::CostMap CostMap;
1.59 ///The type of the costs of the edges.
1.60 typedef typename Traits::Value Value;
1.61 + ///The type of the predecessor map.
1.62 + typedef typename Traits::PredMap PredMap;
1.63 ///The type of the map that stores which edges are in the arborescence.
1.64 typedef typename Traits::ArborescenceMap ArborescenceMap;
1.65
1.66 @@ -157,14 +175,20 @@
1.67
1.68 };
1.69
1.70 - const Graph* graph;
1.71 - const CostMap* cost;
1.72 + const Graph *graph;
1.73 + const CostMap *cost;
1.74
1.75 - ArborescenceMap* _arborescence_map;
1.76 - bool local_arborescence_map;
1.77 + PredMap *_pred;
1.78 + bool local_pred;
1.79
1.80 - typedef typename Graph::template NodeMap<int> LevelMap;
1.81 - LevelMap *_level;
1.82 + ArborescenceMap *_arborescence;
1.83 + bool local_arborescence;
1.84 +
1.85 + typedef typename Graph::template EdgeMap<int> EdgeOrder;
1.86 + EdgeOrder *_edge_order;
1.87 +
1.88 + typedef typename Graph::template NodeMap<int> NodeOrder;
1.89 + NodeOrder *_node_order;
1.90
1.91 typedef typename Graph::template NodeMap<CostEdge> CostEdgeMap;
1.92 CostEdgeMap *_cost_edges;
1.93 @@ -179,7 +203,30 @@
1.94 std::vector<StackLevel> level_stack;
1.95 std::vector<Node> queue;
1.96
1.97 - int node_counter;
1.98 + typedef std::vector<typename Graph::Node> DualNodeList;
1.99 +
1.100 + DualNodeList _dual_node_list;
1.101 +
1.102 + struct DualVariable {
1.103 + int begin, end;
1.104 + Value value;
1.105 +
1.106 + DualVariable(int _begin, int _end, Value _value)
1.107 + : begin(_begin), end(_end), value(_value) {}
1.108 +
1.109 + };
1.110 +
1.111 + typedef std::vector<DualVariable> DualVariables;
1.112 +
1.113 + DualVariables _dual_variables;
1.114 +
1.115 + typedef typename Graph::template NodeMap<int> HeapCrossRef;
1.116 +
1.117 + HeapCrossRef *_heap_cross_ref;
1.118 +
1.119 + typedef BinHeap<Node, int, HeapCrossRef> Heap;
1.120 +
1.121 + Heap *_heap;
1.122
1.123 public:
1.124
1.125 @@ -208,6 +255,26 @@
1.126 typedef MinCostArborescence<Graph, CostMap,
1.127 DefArborescenceMapTraits<T> > Create;
1.128 };
1.129 +
1.130 + template <class T>
1.131 + struct DefPredMapTraits : public Traits {
1.132 + typedef T PredMap;
1.133 + static PredMap *createPredMap(const Graph &)
1.134 + {
1.135 + throw UninitializedParameter();
1.136 + }
1.137 + };
1.138 +
1.139 + /// \brief \ref named-templ-param "Named parameter" for
1.140 + /// setting PredMap type
1.141 + ///
1.142 + /// \ref named-templ-param "Named parameter" for setting
1.143 + /// PredMap type
1.144 + template <class T>
1.145 + struct DefPredMap
1.146 + : public MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > {
1.147 + typedef MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > Create;
1.148 + };
1.149
1.150 /// @}
1.151
1.152 @@ -216,9 +283,10 @@
1.153 /// \param _graph The graph the algorithm will run on.
1.154 /// \param _cost The cost map used by the algorithm.
1.155 MinCostArborescence(const Graph& _graph, const CostMap& _cost)
1.156 - : graph(&_graph), cost(&_cost),
1.157 - _arborescence_map(0), local_arborescence_map(false),
1.158 - _level(0), _cost_edges(0) {}
1.159 + : graph(&_graph), cost(&_cost), _pred(0), local_pred(false),
1.160 + _arborescence(0), local_arborescence(false),
1.161 + _edge_order(0), _node_order(0), _cost_edges(0),
1.162 + _heap_cross_ref(0), _heap(0) {}
1.163
1.164 /// \brief Destructor.
1.165 ~MinCostArborescence() {
1.166 @@ -230,7 +298,24 @@
1.167 /// Sets the arborescence map.
1.168 /// \return \c (*this)
1.169 MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
1.170 - _arborescence_map = &m;
1.171 + if (local_arborescence) {
1.172 + delete _arborescence;
1.173 + }
1.174 + local_arborescence = false;
1.175 + _arborescence = &m;
1.176 + return *this;
1.177 + }
1.178 +
1.179 + /// \brief Sets the arborescence map.
1.180 + ///
1.181 + /// Sets the arborescence map.
1.182 + /// \return \c (*this)
1.183 + MinCostArborescence& predMap(PredMap& m) {
1.184 + if (local_pred) {
1.185 + delete _pred;
1.186 + }
1.187 + local_pred = false;
1.188 + _pred = &m;
1.189 return *this;
1.190 }
1.191
1.192 @@ -246,7 +331,7 @@
1.193 ///
1.194 /// Returns a reference to the arborescence map.
1.195 const ArborescenceMap& arborescenceMap() const {
1.196 - return *_arborescence_map;
1.197 + return *_arborescence;
1.198 }
1.199
1.200 /// \brief Returns true if the edge is in the arborescence.
1.201 @@ -254,23 +339,141 @@
1.202 /// Returns true if the edge is in the arborescence.
1.203 /// \param edge The edge of the graph.
1.204 /// \pre \ref run() must be called before using this function.
1.205 - bool arborescenceEdge(Edge edge) const {
1.206 - return (*_arborescence_map)[edge];
1.207 + bool arborescence(Edge edge) const {
1.208 + return (*_pred)[graph->target(edge)] == edge;
1.209 + }
1.210 +
1.211 + /// \brief Returns a reference to the pred map.
1.212 + ///
1.213 + /// Returns a reference to the pred map.
1.214 + const PredMap& predMap() const {
1.215 + return *_pred;
1.216 + }
1.217 +
1.218 + /// \brief Returns the predecessor edge of the given node.
1.219 + ///
1.220 + /// Returns the predecessor edge of the given node.
1.221 + bool pred(Node node) const {
1.222 + return (*_pred)[node];
1.223 }
1.224
1.225 /// \brief Returns the cost of the arborescence.
1.226 ///
1.227 /// Returns the cost of the arborescence.
1.228 - Value arborescenceCost() const {
1.229 + Value arborescenceValue() const {
1.230 Value sum = 0;
1.231 for (EdgeIt it(*graph); it != INVALID; ++it) {
1.232 - if (arborescenceEdge(it)) {
1.233 + if (arborescence(it)) {
1.234 sum += (*cost)[it];
1.235 }
1.236 }
1.237 return sum;
1.238 }
1.239
1.240 + /// \brief Indicates that a node is reachable from the sources.
1.241 + ///
1.242 + /// Indicates that a node is reachable from the sources.
1.243 + bool reached(Node node) const {
1.244 + return (*_node_order)[node] != -3;
1.245 + }
1.246 +
1.247 + /// \brief Indicates that a node is processed.
1.248 + ///
1.249 + /// Indicates that a node is processed. The arborescence path exists
1.250 + /// from the source to the given node.
1.251 + bool processed(Node node) const {
1.252 + return (*_node_order)[node] == -1;
1.253 + }
1.254 +
1.255 + /// \brief Returns the number of the dual variables in basis.
1.256 + ///
1.257 + /// Returns the number of the dual variables in basis.
1.258 + int dualSize() const {
1.259 + return _dual_variables.size();
1.260 + }
1.261 +
1.262 + /// \brief Returns the value of the dual solution.
1.263 + ///
1.264 + /// Returns the value of the dual solution. It should be
1.265 + /// equal to the arborescence value.
1.266 + Value dualValue() const {
1.267 + Value sum = 0;
1.268 + for (int i = 0; i < (int)_dual_variables.size(); ++i) {
1.269 + sum += _dual_variables[i].value;
1.270 + }
1.271 + return sum;
1.272 + }
1.273 +
1.274 + /// \brief Returns the number of the nodes in the dual variable.
1.275 + ///
1.276 + /// Returns the number of the nodes in the dual variable.
1.277 + int dualSize(int k) const {
1.278 + return _dual_variables[k].end - _dual_variables[k].begin;
1.279 + }
1.280 +
1.281 + /// \brief Returns the value of the dual variable.
1.282 + ///
1.283 + /// Returns the the value of the dual variable.
1.284 + const Value& dualValue(int k) const {
1.285 + return _dual_variables[k].value;
1.286 + }
1.287 +
1.288 + /// \brief Lemon iterator for get a dual variable.
1.289 + ///
1.290 + /// Lemon iterator for get a dual variable. This class provides
1.291 + /// a common style lemon iterator which gives back a subset of
1.292 + /// the nodes.
1.293 + class DualIt {
1.294 + public:
1.295 +
1.296 + /// \brief Constructor.
1.297 + ///
1.298 + /// Constructor for get the nodeset of the variable.
1.299 + DualIt(const MinCostArborescence& algorithm, int variable)
1.300 + : _algorithm(&algorithm), _variable(variable)
1.301 + {
1.302 + _index = _algorithm->_dual_variables[_variable].begin;
1.303 + }
1.304 +
1.305 + /// \brief Invalid constructor.
1.306 + ///
1.307 + /// Invalid constructor.
1.308 + DualIt(Invalid) : _algorithm(0) {}
1.309 +
1.310 + /// \brief Conversion to node.
1.311 + ///
1.312 + /// Conversion to node.
1.313 + operator Node() const {
1.314 + return _algorithm ? _algorithm->_dual_node_list[_index] : INVALID;
1.315 + }
1.316 +
1.317 + /// \brief Increment operator.
1.318 + ///
1.319 + /// Increment operator.
1.320 + DualIt& operator++() {
1.321 + ++_index;
1.322 + if (_algorithm->_dual_variables[_variable].end == _index) {
1.323 + _algorithm = 0;
1.324 + }
1.325 + return *this;
1.326 + }
1.327 +
1.328 + bool operator==(const DualIt& it) const {
1.329 + return (Node)(*this) == (Node)it;
1.330 + }
1.331 + bool operator!=(const DualIt& it) const {
1.332 + return (Node)(*this) != (Node)it;
1.333 + }
1.334 + bool operator<(const DualIt& it) const {
1.335 + return (Node)(*this) < (Node)it;
1.336 + }
1.337 +
1.338 + private:
1.339 + const MinCostArborescence* _algorithm;
1.340 + int _variable;
1.341 + int _index;
1.342 + };
1.343 +
1.344 /// @}
1.345
1.346 /// \name Execution control
1.347 @@ -279,7 +482,7 @@
1.348 /// If you need more control on the execution,
1.349 /// first you must call \ref init(), then you can add several
1.350 /// source nodes with \ref addSource().
1.351 - /// Finally \ref start() will perform the actual path
1.352 + /// Finally \ref start() will perform the arborescence
1.353 /// computation.
1.354
1.355 ///@{
1.356 @@ -290,13 +493,18 @@
1.357 ///
1.358 void init() {
1.359 initStructures();
1.360 + _heap->clear();
1.361 for (NodeIt it(*graph); it != INVALID; ++it) {
1.362 (*_cost_edges)[it].edge = INVALID;
1.363 - (*_level)[it] = -3;
1.364 + _node_order->set(it, -3);
1.365 + _heap_cross_ref->set(it, Heap::PRE_HEAP);
1.366 }
1.367 for (EdgeIt it(*graph); it != INVALID; ++it) {
1.368 - _arborescence_map->set(it, false);
1.369 + _arborescence->set(it, false);
1.370 + _edge_order->set(it, -1);
1.371 }
1.372 + _dual_node_list.clear();
1.373 + _dual_variables.clear();
1.374 }
1.375
1.376 /// \brief Adds a new source node.
1.377 @@ -309,14 +517,15 @@
1.378 Node node = nodes.back();
1.379 nodes.pop_back();
1.380 for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
1.381 - if ((*_level)[graph->target(it)] == -3) {
1.382 - (*_level)[graph->target(it)] = -2;
1.383 - nodes.push_back(graph->target(it));
1.384 - queue.push_back(graph->target(it));
1.385 + Node target = graph->target(it);
1.386 + if ((*_node_order)[target] == -3) {
1.387 + (*_node_order)[target] = -2;
1.388 + nodes.push_back(target);
1.389 + queue.push_back(target);
1.390 }
1.391 }
1.392 }
1.393 - (*_level)[source] = -1;
1.394 + (*_node_order)[source] = -1;
1.395 }
1.396
1.397 /// \brief Processes the next node in the priority queue.
1.398 @@ -327,19 +536,20 @@
1.399 ///
1.400 /// \warning The queue must not be empty!
1.401 Node processNextNode() {
1.402 - node_counter = 0;
1.403 Node node = queue.back();
1.404 queue.pop_back();
1.405 - if ((*_level)[node] == -2) {
1.406 + if ((*_node_order)[node] == -2) {
1.407 Edge edge = prepare(node);
1.408 - while ((*_level)[graph->source(edge)] != -1) {
1.409 - if ((*_level)[graph->source(edge)] >= 0) {
1.410 - edge = contract(bottom((*_level)[graph->source(edge)]));
1.411 + Node source = graph->source(edge);
1.412 + while ((*_node_order)[source] != -1) {
1.413 + if ((*_node_order)[source] >= 0) {
1.414 + edge = contract(source);
1.415 } else {
1.416 - edge = prepare(graph->source(edge));
1.417 + edge = prepare(source);
1.418 }
1.419 + source = graph->source(edge);
1.420 }
1.421 - finalize(graph->target(edge));
1.422 + finalize(edge);
1.423 level_stack.clear();
1.424 }
1.425 return node;
1.426 @@ -400,46 +610,74 @@
1.427 protected:
1.428
1.429 void initStructures() {
1.430 - if (!_arborescence_map) {
1.431 - local_arborescence_map = true;
1.432 - _arborescence_map = Traits::createArborescenceMap(*graph);
1.433 + if (!_pred) {
1.434 + local_pred = true;
1.435 + _pred = Traits::createPredMap(*graph);
1.436 }
1.437 - if (!_level) {
1.438 - _level = new LevelMap(*graph);
1.439 + if (!_arborescence) {
1.440 + local_arborescence = true;
1.441 + _arborescence = Traits::createArborescenceMap(*graph);
1.442 + }
1.443 + if (!_edge_order) {
1.444 + _edge_order = new EdgeOrder(*graph);
1.445 + }
1.446 + if (!_node_order) {
1.447 + _node_order = new NodeOrder(*graph);
1.448 }
1.449 if (!_cost_edges) {
1.450 _cost_edges = new CostEdgeMap(*graph);
1.451 }
1.452 + if (!_heap_cross_ref) {
1.453 + _heap_cross_ref = new HeapCrossRef(*graph, -1);
1.454 + }
1.455 + if (!_heap) {
1.456 + _heap = new Heap(*_heap_cross_ref);
1.457 + }
1.458 }
1.459
1.460 void destroyStructures() {
1.461 - if (_level) {
1.462 - delete _level;
1.463 + if (local_arborescence) {
1.464 + delete _arborescence;
1.465 + }
1.466 + if (local_pred) {
1.467 + delete _pred;
1.468 + }
1.469 + if (!_edge_order) {
1.470 + delete _edge_order;
1.471 + }
1.472 + if (_node_order) {
1.473 + delete _node_order;
1.474 }
1.475 if (!_cost_edges) {
1.476 delete _cost_edges;
1.477 }
1.478 - if (local_arborescence_map) {
1.479 - delete _arborescence_map;
1.480 + if (!_heap) {
1.481 + delete _heap;
1.482 + }
1.483 + if (!_heap_cross_ref) {
1.484 + delete _heap_cross_ref;
1.485 }
1.486 }
1.487
1.488 Edge prepare(Node node) {
1.489 std::vector<Node> nodes;
1.490 - (*_level)[node] = node_counter;
1.491 + (*_node_order)[node] = _dual_node_list.size();
1.492 + StackLevel level;
1.493 + level.node_level = _dual_node_list.size();
1.494 + _dual_node_list.push_back(node);
1.495 for (InEdgeIt it(*graph, node); it != INVALID; ++it) {
1.496 Edge edge = it;
1.497 + Node source = graph->source(edge);
1.498 Value value = (*cost)[it];
1.499 - if (graph->source(edge) == node ||
1.500 - (*_level)[graph->source(edge)] == -3) continue;
1.501 - if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {
1.502 - (*_cost_edges)[graph->source(edge)].edge = edge;
1.503 - (*_cost_edges)[graph->source(edge)].value = value;
1.504 - nodes.push_back(graph->source(edge));
1.505 + if (source == node || (*_node_order)[source] == -3) continue;
1.506 + if ((*_cost_edges)[source].edge == INVALID) {
1.507 + (*_cost_edges)[source].edge = edge;
1.508 + (*_cost_edges)[source].value = value;
1.509 + nodes.push_back(source);
1.510 } else {
1.511 - if ((*_cost_edges)[graph->source(edge)].value > value) {
1.512 - (*_cost_edges)[graph->source(edge)].edge = edge;
1.513 - (*_cost_edges)[graph->source(edge)].value = value;
1.514 + if ((*_cost_edges)[source].value > value) {
1.515 + (*_cost_edges)[source].edge = edge;
1.516 + (*_cost_edges)[source].value = value;
1.517 }
1.518 }
1.519 }
1.520 @@ -449,35 +687,37 @@
1.521 minimum = (*_cost_edges)[nodes[i]];
1.522 }
1.523 }
1.524 - StackLevel level;
1.525 - level.node_level = node_counter;
1.526 + _edge_order->set(minimum.edge, _dual_variables.size());
1.527 + DualVariable var(_dual_node_list.size() - 1,
1.528 + _dual_node_list.size(), minimum.value);
1.529 + _dual_variables.push_back(var);
1.530 for (int i = 0; i < (int)nodes.size(); ++i) {
1.531 (*_cost_edges)[nodes[i]].value -= minimum.value;
1.532 level.edges.push_back((*_cost_edges)[nodes[i]]);
1.533 (*_cost_edges)[nodes[i]].edge = INVALID;
1.534 }
1.535 level_stack.push_back(level);
1.536 - ++node_counter;
1.537 - _arborescence_map->set(minimum.edge, true);
1.538 return minimum.edge;
1.539 }
1.540
1.541 - Edge contract(int node_bottom) {
1.542 + Edge contract(Node node) {
1.543 + int node_bottom = bottom(node);
1.544 std::vector<Node> nodes;
1.545 while (!level_stack.empty() &&
1.546 level_stack.back().node_level >= node_bottom) {
1.547 for (int i = 0; i < (int)level_stack.back().edges.size(); ++i) {
1.548 Edge edge = level_stack.back().edges[i].edge;
1.549 + Node source = graph->source(edge);
1.550 Value value = level_stack.back().edges[i].value;
1.551 - if ((*_level)[graph->source(edge)] >= node_bottom) continue;
1.552 - if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {
1.553 - (*_cost_edges)[graph->source(edge)].edge = edge;
1.554 - (*_cost_edges)[graph->source(edge)].value = value;
1.555 - nodes.push_back(graph->source(edge));
1.556 + if ((*_node_order)[source] >= node_bottom) continue;
1.557 + if ((*_cost_edges)[source].edge == INVALID) {
1.558 + (*_cost_edges)[source].edge = edge;
1.559 + (*_cost_edges)[source].value = value;
1.560 + nodes.push_back(source);
1.561 } else {
1.562 - if ((*_cost_edges)[graph->source(edge)].value > value) {
1.563 - (*_cost_edges)[graph->source(edge)].edge = edge;
1.564 - (*_cost_edges)[graph->source(edge)].value = value;
1.565 + if ((*_cost_edges)[source].value > value) {
1.566 + (*_cost_edges)[source].edge = edge;
1.567 + (*_cost_edges)[source].value = value;
1.568 }
1.569 }
1.570 }
1.571 @@ -489,6 +729,9 @@
1.572 minimum = (*_cost_edges)[nodes[i]];
1.573 }
1.574 }
1.575 + _edge_order->set(minimum.edge, _dual_variables.size());
1.576 + DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
1.577 + _dual_variables.push_back(var);
1.578 StackLevel level;
1.579 level.node_level = node_bottom;
1.580 for (int i = 0; i < (int)nodes.size(); ++i) {
1.581 @@ -497,34 +740,45 @@
1.582 (*_cost_edges)[nodes[i]].edge = INVALID;
1.583 }
1.584 level_stack.push_back(level);
1.585 - _arborescence_map->set(minimum.edge, true);
1.586 return minimum.edge;
1.587 }
1.588
1.589 - int bottom(int level) {
1.590 + int bottom(Node node) {
1.591 int k = level_stack.size() - 1;
1.592 - while (level_stack[k].node_level > level) {
1.593 + while (level_stack[k].node_level > (*_node_order)[node]) {
1.594 --k;
1.595 }
1.596 return level_stack[k].node_level;
1.597 }
1.598
1.599 - void finalize(Node source) {
1.600 - std::vector<Node> nodes;
1.601 - nodes.push_back(source);
1.602 - while (!nodes.empty()) {
1.603 - Node node = nodes.back();
1.604 - nodes.pop_back();
1.605 - for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
1.606 - if ((*_level)[graph->target(it)] >= 0 && (*_arborescence_map)[it]) {
1.607 - (*_level)[graph->target(it)] = -1;
1.608 - nodes.push_back(graph->target(it));
1.609 - } else {
1.610 - _arborescence_map->set(it, false);
1.611 + void finalize(Edge edge) {
1.612 + Node node = graph->target(edge);
1.613 + _heap->push(node, (*_edge_order)[edge]);
1.614 + _pred->set(node, edge);
1.615 + while (!_heap->empty()) {
1.616 + Node source = _heap->top();
1.617 + _heap->pop();
1.618 + _node_order->set(source, -1);
1.619 + for (OutEdgeIt it(*graph, source); it != INVALID; ++it) {
1.620 + if ((*_edge_order)[it] < 0) continue;
1.621 + Node target = graph->target(it);
1.622 + switch(_heap->state(target)) {
1.623 + case Heap::PRE_HEAP:
1.624 + _heap->push(target, (*_edge_order)[it]);
1.625 + _pred->set(target, it);
1.626 + break;
1.627 + case Heap::IN_HEAP:
1.628 + if ((*_edge_order)[it] < (*_heap)[target]) {
1.629 + _heap->decrease(target, (*_edge_order)[it]);
1.630 + _pred->set(target, it);
1.631 + }
1.632 + break;
1.633 + case Heap::POST_HEAP:
1.634 + break;
1.635 }
1.636 }
1.637 + _arborescence->set((*_pred)[source], true);
1.638 }
1.639 - (*_level)[source] = -1;
1.640 }
1.641
1.642 };
1.643 @@ -551,11 +805,9 @@
1.644 ::Create mca(graph, cost);
1.645 mca.arborescenceMap(arborescence);
1.646 mca.run(source);
1.647 - return mca.arborescenceCost();
1.648 + return mca.arborescenceValue();
1.649 }
1.650
1.651 }
1.652
1.653 #endif
1.654 -
1.655 -// Hilbert - Huang
2.1 --- a/test/Makefile.am Thu Mar 30 15:34:56 2006 +0000
2.2 +++ b/test/Makefile.am Fri Mar 31 11:10:44 2006 +0000
2.3 @@ -13,6 +13,7 @@
2.4
2.5 check_PROGRAMS = \
2.6 all_pairs_shortest_path_test \
2.7 + arborescence_test \
2.8 bfs_test \
2.9 counter_test \
2.10 dfs_test \
2.11 @@ -52,6 +53,7 @@
2.12 XFAIL_TESTS = test_tools_fail$(EXEEXT)
2.13
2.14 all_pairs_shortest_path_test_SOURCES = all_pairs_shortest_path_test.cc
2.15 +arborescence_test_SOURCES = arborescence_test.cc
2.16 bfs_test_SOURCES = bfs_test.cc
2.17 counter_test_SOURCES = counter_test.cc
2.18 dfs_test_SOURCES = dfs_test.cc
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
3.2 +++ b/test/arborescence_test.cc Fri Mar 31 11:10:44 2006 +0000
3.3 @@ -0,0 +1,104 @@
3.4 +#include <iostream>
3.5 +#include <set>
3.6 +#include <vector>
3.7 +#include <iterator>
3.8 +
3.9 +#include <cmath>
3.10 +#include <cstdlib>
3.11 +
3.12 +#include <lemon/smart_graph.h>
3.13 +#include <lemon/min_cost_arborescence.h>
3.14 +
3.15 +#include <lemon/graph_utils.h>
3.16 +#include <lemon/time_measure.h>
3.17 +
3.18 +#include <lemon/tolerance.h>
3.19 +
3.20 +using namespace lemon;
3.21 +using namespace std;
3.22 +
3.23 +int main(int argc, const char *argv[]) {
3.24 + srand(time(0));
3.25 + typedef SmartGraph Graph;
3.26 + GRAPH_TYPEDEFS(Graph);
3.27 +
3.28 + typedef Graph::EdgeMap<double> CostMap;
3.29 +
3.30 + const int n = argc > 1 ? atoi(argv[1]) : 100;
3.31 + const int e = argc > 2 ? atoi(argv[2]) : (int)(n * log(n));
3.32 +
3.33 + Graph graph;
3.34 + CostMap cost(graph);
3.35 + vector<Node> nodes;
3.36 +
3.37 + for (int i = 0; i < n; ++i) {
3.38 + nodes.push_back(graph.addNode());
3.39 + }
3.40 +
3.41 + for (int i = 0; i < e; ++i) {
3.42 + int s = (int)(n * (double)rand() / (RAND_MAX + 1.0));
3.43 + int t = (int)(n * (double)rand() / (RAND_MAX + 1.0));
3.44 + double c = rand() / (1.0 + RAND_MAX) * 100.0 + 20.0;
3.45 + Edge edge = graph.addEdge(nodes[s], nodes[t]);
3.46 + cost[edge] = c;
3.47 + }
3.48 +
3.49 + Node source = nodes[(int)(n * (double)rand() / (RAND_MAX + 1.0))];
3.50 +
3.51 + MinCostArborescence<Graph, CostMap> mca(graph, cost);
3.52 + mca.run(source);
3.53 +
3.54 + vector<pair<double, set<Node> > > dualSolution(mca.dualSize());
3.55 +
3.56 + for (int i = 0; i < mca.dualSize(); ++i) {
3.57 + dualSolution[i].first = mca.dualValue(i);
3.58 + for (MinCostArborescence<Graph, CostMap>::DualIt it(mca, i);
3.59 + it != INVALID; ++it) {
3.60 + dualSolution[i].second.insert(it);
3.61 + }
3.62 + }
3.63 +
3.64 + Tolerance<double> tolerance;
3.65 +
3.66 + for (EdgeIt it(graph); it != INVALID; ++it) {
3.67 + if (mca.reached(graph.source(it))) {
3.68 + double sum = 0.0;
3.69 + for (int i = 0; i < (int)dualSolution.size(); ++i) {
3.70 + if (dualSolution[i].second.find(graph.target(it))
3.71 + != dualSolution[i].second.end() &&
3.72 + dualSolution[i].second.find(graph.source(it))
3.73 + == dualSolution[i].second.end()) {
3.74 + sum += dualSolution[i].first;
3.75 + }
3.76 + }
3.77 + if (mca.arborescence(it)) {
3.78 + LEMON_ASSERT(!tolerance.less(sum, cost[it]), "INVALID DUAL");
3.79 + }
3.80 + LEMON_ASSERT(!tolerance.less(cost[it], sum), "INVALID DUAL");
3.81 + }
3.82 + }
3.83 +
3.84 +
3.85 + LEMON_ASSERT(!tolerance.different(mca.dualValue(), mca.arborescenceValue()),
3.86 + "INVALID DUAL");
3.87 +
3.88 +
3.89 + LEMON_ASSERT(mca.reached(source), "INVALID ARBORESCENCE");
3.90 + for (EdgeIt it(graph); it != INVALID; ++it) {
3.91 + LEMON_ASSERT(!mca.reached(graph.source(it)) ||
3.92 + mca.reached(graph.target(it)), "INVALID ARBORESCENCE");
3.93 + }
3.94 +
3.95 + for (NodeIt it(graph); it != INVALID; ++it) {
3.96 + if (!mca.reached(it)) continue;
3.97 + int cnt = 0;
3.98 + for (InEdgeIt jt(graph, it); jt != INVALID; ++jt) {
3.99 + if (mca.arborescence(jt)) {
3.100 + ++cnt;
3.101 + }
3.102 + }
3.103 + LEMON_ASSERT((it == source ? cnt == 0 : cnt == 1), "INVALID ARBORESCENCE");
3.104 + }
3.105 +
3.106 + return 0;
3.107 +}