1.1 --- a/lemon/Makefile.am Mon Mar 27 08:01:10 2006 +0000
1.2 +++ b/lemon/Makefile.am Mon Mar 27 08:12:01 2006 +0000
1.3 @@ -59,6 +59,7 @@
1.4 matrix_maps.h \
1.5 map_iterator.h \
1.6 max_matching.h \
1.7 + min_cost_arborescence.h \
1.8 min_cost_flow.h \
1.9 min_cut.h \
1.10 suurballe.h \
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/lemon/min_cost_arborescence.h Mon Mar 27 08:12:01 2006 +0000
2.3 @@ -0,0 +1,561 @@
2.4 +/* -*- C++ -*-
2.5 + *
2.6 + * This file is a part of LEMON, a generic C++ optimization library
2.7 + *
2.8 + * Copyright (C) 2003-2006
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 +
2.33 +namespace lemon {
2.34 +
2.35 +
2.36 + /// \brief Default traits class of MinCostArborescence class.
2.37 + ///
2.38 + /// Default traits class of MinCostArborescence class.
2.39 + /// \param _Graph Graph type.
2.40 + /// \param _CostMap Type of cost map.
2.41 + template <class _Graph, class _CostMap>
2.42 + struct MinCostArborescenceDefaultTraits{
2.43 +
2.44 + /// \brief The graph type the algorithm runs on.
2.45 + typedef _Graph Graph;
2.46 +
2.47 + /// \brief The type of the map that stores the edge costs.
2.48 + ///
2.49 + /// The type of the map that stores the edge costs.
2.50 + /// It must meet the \ref concept::ReadMap "ReadMap" concept.
2.51 + typedef _CostMap CostMap;
2.52 +
2.53 + /// \brief The value type of the costs.
2.54 + ///
2.55 + /// The value type of the costs.
2.56 + typedef typename CostMap::Value Value;
2.57 +
2.58 + /// \brief The type of the map that stores which edges are
2.59 + /// in the arborescence.
2.60 + ///
2.61 + /// The type of the map that stores which edges are in the arborescence.
2.62 + /// It must meet the \ref concept::ReadWriteMap "ReadWriteMap" concept.
2.63 + /// Initially it will be setted to false on each edge. The algorithm
2.64 + /// may set each value one time to true and maybe after it to false again.
2.65 + /// Therefore you cannot use maps like BackInserteBoolMap with this
2.66 + /// algorithm.
2.67 + typedef typename Graph::template EdgeMap<bool> ArborescenceMap;
2.68 +
2.69 + /// \brief Instantiates a ArborescenceMap.
2.70 + ///
2.71 + /// This function instantiates a \ref ArborescenceMap.
2.72 + /// \param _graph is the graph, to which we would like to define the
2.73 + /// ArborescenceMap.
2.74 + static ArborescenceMap *createArborescenceMap(const Graph &_graph){
2.75 + return new ArborescenceMap(_graph);
2.76 + }
2.77 +
2.78 + };
2.79 +
2.80 + /// \ingroup spantree
2.81 + ///
2.82 + /// \brief %MinCostArborescence algorithm class.
2.83 + ///
2.84 + /// This class provides an efficient implementation of
2.85 + /// %MinCostArborescence algorithm. The arborescence is a tree
2.86 + /// which is directed from a given source node of the graph. One or
2.87 + /// more sources should be given for the algorithm and it will calculate
2.88 + /// the minimum cost subgraph which are union of arborescences with the
2.89 + /// given sources and spans all the nodes which are reachable from the
2.90 + /// sources. The time complexity of the algorithm is O(n^2 + e).
2.91 + ///
2.92 + /// \param _Graph The graph type the algorithm runs on. The default value
2.93 + /// is \ref ListGraph. The value of _Graph is not used directly by
2.94 + /// MinCostArborescence, it is only passed to
2.95 + /// \ref MinCostArborescenceDefaultTraits.
2.96 + /// \param _CostMap This read-only EdgeMap determines the costs of the
2.97 + /// edges. It is read once for each edge, so the map may involve in
2.98 + /// relatively time consuming process to compute the edge cost if
2.99 + /// it is necessary. The default map type is \ref
2.100 + /// concept::StaticGraph::EdgeMap "Graph::EdgeMap<int>". The value
2.101 + /// of _CostMap is not used directly by MinCostArborescence,
2.102 + /// it is only passed to \ref MinCostArborescenceDefaultTraits.
2.103 + /// \param _Traits Traits class to set various data types used
2.104 + /// by the algorithm. The default traits class is
2.105 + /// \ref MinCostArborescenceDefaultTraits
2.106 + /// "MinCostArborescenceDefaultTraits<_Graph,_CostMap>". See \ref
2.107 + /// MinCostArborescenceDefaultTraits for the documentation of a
2.108 + /// MinCostArborescence traits class.
2.109 + ///
2.110 + /// \author Balazs Dezso
2.111 +#ifndef DOXYGEN
2.112 + template <typename _Graph = ListGraph,
2.113 + typename _CostMap = typename _Graph::template EdgeMap<int>,
2.114 + typename _Traits =
2.115 + MinCostArborescenceDefaultTraits<_Graph, _CostMap> >
2.116 +#else
2.117 + template <typename _Graph, typename _CostMap, typedef _Traits>
2.118 +#endif
2.119 + class MinCostArborescence {
2.120 + public:
2.121 +
2.122 + /// \brief \ref Exception for uninitialized parameters.
2.123 + ///
2.124 + /// This error represents problems in the initialization
2.125 + /// of the parameters of the algorithms.
2.126 + class UninitializedParameter : public lemon::UninitializedParameter {
2.127 + public:
2.128 + virtual const char* exceptionName() const {
2.129 + return "lemon::MinCostArborescence::UninitializedParameter";
2.130 + }
2.131 + };
2.132 +
2.133 + /// The traits.
2.134 + typedef _Traits Traits;
2.135 + /// The type of the underlying graph.
2.136 + typedef typename Traits::Graph Graph;
2.137 + /// The type of the map that stores the edge costs.
2.138 + typedef typename Traits::CostMap CostMap;
2.139 + ///The type of the costs of the edges.
2.140 + typedef typename Traits::Value Value;
2.141 + ///The type of the map that stores which edges are in the arborescence.
2.142 + typedef typename Traits::ArborescenceMap ArborescenceMap;
2.143 +
2.144 + protected:
2.145 +
2.146 + typedef typename Graph::Node Node;
2.147 + typedef typename Graph::Edge Edge;
2.148 + typedef typename Graph::NodeIt NodeIt;
2.149 + typedef typename Graph::EdgeIt EdgeIt;
2.150 + typedef typename Graph::InEdgeIt InEdgeIt;
2.151 + typedef typename Graph::OutEdgeIt OutEdgeIt;
2.152 +
2.153 + struct CostEdge {
2.154 +
2.155 + Edge edge;
2.156 + Value value;
2.157 +
2.158 + CostEdge() {}
2.159 + CostEdge(Edge _edge, Value _value) : edge(_edge), value(_value) {}
2.160 +
2.161 + };
2.162 +
2.163 + const Graph* graph;
2.164 + const CostMap* cost;
2.165 +
2.166 + ArborescenceMap* _arborescence_map;
2.167 + bool local_arborescence_map;
2.168 +
2.169 + typedef typename Graph::template NodeMap<int> LevelMap;
2.170 + LevelMap *_level;
2.171 +
2.172 + typedef typename Graph::template NodeMap<CostEdge> CostEdgeMap;
2.173 + CostEdgeMap *_cost_edges;
2.174 +
2.175 + struct StackLevel {
2.176 +
2.177 + std::vector<CostEdge> edges;
2.178 + int node_level;
2.179 +
2.180 + };
2.181 +
2.182 + std::vector<StackLevel> level_stack;
2.183 + std::vector<Node> queue;
2.184 +
2.185 + int node_counter;
2.186 +
2.187 + public:
2.188 +
2.189 + /// \name Named template parameters
2.190 +
2.191 + /// @{
2.192 +
2.193 + template <class T>
2.194 + struct DefArborescenceMapTraits : public Traits {
2.195 + typedef T ArborescenceMap;
2.196 + static ArborescenceMap *createArborescenceMap(const Graph &)
2.197 + {
2.198 + throw UninitializedParameter();
2.199 + }
2.200 + };
2.201 +
2.202 + /// \brief \ref named-templ-param "Named parameter" for
2.203 + /// setting ArborescenceMap type
2.204 + ///
2.205 + /// \ref named-templ-param "Named parameter" for setting
2.206 + /// ArborescenceMap type
2.207 + template <class T>
2.208 + struct DefArborescenceMap
2.209 + : public MinCostArborescence<Graph, CostMap,
2.210 + DefArborescenceMapTraits<T> > {
2.211 + typedef MinCostArborescence<Graph, CostMap,
2.212 + DefArborescenceMapTraits<T> > Create;
2.213 + };
2.214 +
2.215 + /// @}
2.216 +
2.217 + /// \brief Constructor.
2.218 + ///
2.219 + /// \param _graph The graph the algorithm will run on.
2.220 + /// \param _cost The cost map used by the algorithm.
2.221 + MinCostArborescence(const Graph& _graph, const CostMap& _cost)
2.222 + : graph(&_graph), cost(&_cost),
2.223 + _arborescence_map(0), local_arborescence_map(false),
2.224 + _level(0), _cost_edges(0) {}
2.225 +
2.226 + /// \brief Destructor.
2.227 + ~MinCostArborescence() {
2.228 + destroyStructures();
2.229 + }
2.230 +
2.231 + /// \brief Sets the arborescence map.
2.232 + ///
2.233 + /// Sets the arborescence map.
2.234 + /// \return \c (*this)
2.235 + MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
2.236 + _arborescence_map = &m;
2.237 + return *this;
2.238 + }
2.239 +
2.240 + /// \name Query Functions
2.241 + /// The result of the %MinCostArborescence algorithm can be obtained
2.242 + /// using these functions.\n
2.243 + /// Before the use of these functions,
2.244 + /// either run() or start() must be called.
2.245 +
2.246 + /// @{
2.247 +
2.248 + /// \brief Returns a reference to the arborescence map.
2.249 + ///
2.250 + /// Returns a reference to the arborescence map.
2.251 + const ArborescenceMap& arborescenceMap() const {
2.252 + return *_arborescence_map;
2.253 + }
2.254 +
2.255 + /// \brief Returns true if the edge is in the arborescence.
2.256 + ///
2.257 + /// Returns true if the edge is in the arborescence.
2.258 + /// \param edge The edge of the graph.
2.259 + /// \pre \ref run() must be called before using this function.
2.260 + bool arborescenceEdge(Edge edge) const {
2.261 + return (*_arborescence_map)[edge];
2.262 + }
2.263 +
2.264 + /// \brief Returns the cost of the arborescence.
2.265 + ///
2.266 + /// Returns the cost of the arborescence.
2.267 + Value arborescenceCost() const {
2.268 + Value sum = 0;
2.269 + for (EdgeIt it(*graph); it != INVALID; ++it) {
2.270 + if (arborescenceEdge(it)) {
2.271 + sum += (*cost)[it];
2.272 + }
2.273 + }
2.274 + return sum;
2.275 + }
2.276 +
2.277 + /// @}
2.278 +
2.279 + /// \name Execution control
2.280 + /// The simplest way to execute the algorithm is to use
2.281 + /// one of the member functions called \c run(...). \n
2.282 + /// If you need more control on the execution,
2.283 + /// first you must call \ref init(), then you can add several
2.284 + /// source nodes with \ref addSource().
2.285 + /// Finally \ref start() will perform the actual path
2.286 + /// computation.
2.287 +
2.288 + ///@{
2.289 +
2.290 + /// \brief Initializes the internal data structures.
2.291 + ///
2.292 + /// Initializes the internal data structures.
2.293 + ///
2.294 + void init() {
2.295 + initStructures();
2.296 + for (NodeIt it(*graph); it != INVALID; ++it) {
2.297 + (*_cost_edges)[it].edge = INVALID;
2.298 + (*_level)[it] = -3;
2.299 + }
2.300 + for (EdgeIt it(*graph); it != INVALID; ++it) {
2.301 + _arborescence_map->set(it, false);
2.302 + }
2.303 + }
2.304 +
2.305 + /// \brief Adds a new source node.
2.306 + ///
2.307 + /// Adds a new source node to the algorithm.
2.308 + void addSource(Node source) {
2.309 + std::vector<Node> nodes;
2.310 + nodes.push_back(source);
2.311 + while (!nodes.empty()) {
2.312 + Node node = nodes.back();
2.313 + nodes.pop_back();
2.314 + for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
2.315 + if ((*_level)[graph->target(it)] == -3) {
2.316 + (*_level)[graph->target(it)] = -2;
2.317 + nodes.push_back(graph->target(it));
2.318 + queue.push_back(graph->target(it));
2.319 + }
2.320 + }
2.321 + }
2.322 + (*_level)[source] = -1;
2.323 + }
2.324 +
2.325 + /// \brief Processes the next node in the priority queue.
2.326 + ///
2.327 + /// Processes the next node in the priority queue.
2.328 + ///
2.329 + /// \return The processed node.
2.330 + ///
2.331 + /// \warning The queue must not be empty!
2.332 + Node processNextNode() {
2.333 + node_counter = 0;
2.334 + Node node = queue.back();
2.335 + queue.pop_back();
2.336 + if ((*_level)[node] == -2) {
2.337 + Edge edge = prepare(node);
2.338 + while ((*_level)[graph->source(edge)] != -1) {
2.339 + if ((*_level)[graph->source(edge)] >= 0) {
2.340 + edge = contract(bottom((*_level)[graph->source(edge)]));
2.341 + } else {
2.342 + edge = prepare(graph->source(edge));
2.343 + }
2.344 + }
2.345 + finalize(graph->target(edge));
2.346 + level_stack.clear();
2.347 + }
2.348 + return node;
2.349 + }
2.350 +
2.351 + /// \brief Returns the number of the nodes to be processed.
2.352 + ///
2.353 + /// Returns the number of the nodes to be processed.
2.354 + int queueSize() const {
2.355 + return queue.size();
2.356 + }
2.357 +
2.358 + /// \brief Returns \c false if there are nodes to be processed.
2.359 + ///
2.360 + /// Returns \c false if there are nodes to be processed.
2.361 + bool emptyQueue() const {
2.362 + return queue.empty();
2.363 + }
2.364 +
2.365 + /// \brief Executes the algorithm.
2.366 + ///
2.367 + /// Executes the algorithm.
2.368 + ///
2.369 + /// \pre init() must be called and at least one node should be added
2.370 + /// with addSource() before using this function.
2.371 + ///
2.372 + ///\note mca.start() is just a shortcut of the following code.
2.373 + ///\code
2.374 + ///while (!mca.emptyQueue()) {
2.375 + /// mca.processNextNode();
2.376 + ///}
2.377 + ///\endcode
2.378 + void start() {
2.379 + while (!emptyQueue()) {
2.380 + processNextNode();
2.381 + }
2.382 + }
2.383 +
2.384 + /// \brief Runs %MinCostArborescence algorithm from node \c s.
2.385 + ///
2.386 + /// This method runs the %MinCostArborescence algorithm from
2.387 + /// a root node \c s.
2.388 + ///
2.389 + ///\note mca.run(s) is just a shortcut of the following code.
2.390 + ///\code
2.391 + ///mca.init();
2.392 + ///mca.addSource(s);
2.393 + ///mca.start();
2.394 + ///\endcode
2.395 + void run(Node node) {
2.396 + init();
2.397 + addSource(node);
2.398 + start();
2.399 + }
2.400 +
2.401 + ///@}
2.402 +
2.403 + protected:
2.404 +
2.405 + void initStructures() {
2.406 + if (!_arborescence_map) {
2.407 + local_arborescence_map = true;
2.408 + _arborescence_map = Traits::createArborescenceMap(*graph);
2.409 + }
2.410 + if (!_level) {
2.411 + _level = new LevelMap(*graph);
2.412 + }
2.413 + if (!_cost_edges) {
2.414 + _cost_edges = new CostEdgeMap(*graph);
2.415 + }
2.416 + }
2.417 +
2.418 + void destroyStructures() {
2.419 + if (_level) {
2.420 + delete _level;
2.421 + }
2.422 + if (!_cost_edges) {
2.423 + delete _cost_edges;
2.424 + }
2.425 + if (local_arborescence_map) {
2.426 + delete _arborescence_map;
2.427 + }
2.428 + }
2.429 +
2.430 + Edge prepare(Node node) {
2.431 + std::vector<Node> nodes;
2.432 + (*_level)[node] = node_counter;
2.433 + for (InEdgeIt it(*graph, node); it != INVALID; ++it) {
2.434 + Edge edge = it;
2.435 + Value value = (*cost)[it];
2.436 + if (graph->source(edge) == node ||
2.437 + (*_level)[graph->source(edge)] == -3) continue;
2.438 + if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {
2.439 + (*_cost_edges)[graph->source(edge)].edge = edge;
2.440 + (*_cost_edges)[graph->source(edge)].value = value;
2.441 + nodes.push_back(graph->source(edge));
2.442 + } else {
2.443 + if ((*_cost_edges)[graph->source(edge)].value > value) {
2.444 + (*_cost_edges)[graph->source(edge)].edge = edge;
2.445 + (*_cost_edges)[graph->source(edge)].value = value;
2.446 + }
2.447 + }
2.448 + }
2.449 + CostEdge minimum = (*_cost_edges)[nodes[0]];
2.450 + for (int i = 1; i < (int)nodes.size(); ++i) {
2.451 + if ((*_cost_edges)[nodes[i]].value < minimum.value) {
2.452 + minimum = (*_cost_edges)[nodes[i]];
2.453 + }
2.454 + }
2.455 + StackLevel level;
2.456 + level.node_level = node_counter;
2.457 + for (int i = 0; i < (int)nodes.size(); ++i) {
2.458 + (*_cost_edges)[nodes[i]].value -= minimum.value;
2.459 + level.edges.push_back((*_cost_edges)[nodes[i]]);
2.460 + (*_cost_edges)[nodes[i]].edge = INVALID;
2.461 + }
2.462 + level_stack.push_back(level);
2.463 + ++node_counter;
2.464 + _arborescence_map->set(minimum.edge, true);
2.465 + return minimum.edge;
2.466 + }
2.467 +
2.468 + Edge contract(int node_bottom) {
2.469 + std::vector<Node> nodes;
2.470 + while (!level_stack.empty() &&
2.471 + level_stack.back().node_level >= node_bottom) {
2.472 + for (int i = 0; i < (int)level_stack.back().edges.size(); ++i) {
2.473 + Edge edge = level_stack.back().edges[i].edge;
2.474 + Value value = level_stack.back().edges[i].value;
2.475 + if ((*_level)[graph->source(edge)] >= node_bottom) continue;
2.476 + if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {
2.477 + (*_cost_edges)[graph->source(edge)].edge = edge;
2.478 + (*_cost_edges)[graph->source(edge)].value = value;
2.479 + nodes.push_back(graph->source(edge));
2.480 + } else {
2.481 + if ((*_cost_edges)[graph->source(edge)].value > value) {
2.482 + (*_cost_edges)[graph->source(edge)].edge = edge;
2.483 + (*_cost_edges)[graph->source(edge)].value = value;
2.484 + }
2.485 + }
2.486 + }
2.487 + level_stack.pop_back();
2.488 + }
2.489 + CostEdge minimum = (*_cost_edges)[nodes[0]];
2.490 + for (int i = 1; i < (int)nodes.size(); ++i) {
2.491 + if ((*_cost_edges)[nodes[i]].value < minimum.value) {
2.492 + minimum = (*_cost_edges)[nodes[i]];
2.493 + }
2.494 + }
2.495 + StackLevel level;
2.496 + level.node_level = node_bottom;
2.497 + for (int i = 0; i < (int)nodes.size(); ++i) {
2.498 + (*_cost_edges)[nodes[i]].value -= minimum.value;
2.499 + level.edges.push_back((*_cost_edges)[nodes[i]]);
2.500 + (*_cost_edges)[nodes[i]].edge = INVALID;
2.501 + }
2.502 + level_stack.push_back(level);
2.503 + _arborescence_map->set(minimum.edge, true);
2.504 + return minimum.edge;
2.505 + }
2.506 +
2.507 + int bottom(int level) {
2.508 + int k = level_stack.size() - 1;
2.509 + while (level_stack[k].node_level > level) {
2.510 + --k;
2.511 + }
2.512 + return level_stack[k].node_level;
2.513 + }
2.514 +
2.515 + void finalize(Node source) {
2.516 + std::vector<Node> nodes;
2.517 + nodes.push_back(source);
2.518 + while (!nodes.empty()) {
2.519 + Node node = nodes.back();
2.520 + nodes.pop_back();
2.521 + for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
2.522 + if ((*_level)[graph->target(it)] >= 0 && (*_arborescence_map)[it]) {
2.523 + (*_level)[graph->target(it)] = -1;
2.524 + nodes.push_back(graph->target(it));
2.525 + } else {
2.526 + _arborescence_map->set(it, false);
2.527 + }
2.528 + }
2.529 + }
2.530 + (*_level)[source] = -1;
2.531 + }
2.532 +
2.533 + };
2.534 +
2.535 + /// \ingroup spantree
2.536 + ///
2.537 + /// \brief Function type interface for MinCostArborescence algorithm.
2.538 + ///
2.539 + /// Function type interface for MinCostArborescence algorithm.
2.540 + /// \param graph The Graph that the algorithm runs on.
2.541 + /// \param cost The CostMap of the edges.
2.542 + /// \param source The source of the arborescence.
2.543 + /// \retval arborescence The bool EdgeMap which stores the arborescence.
2.544 + /// \return The cost of the arborescence.
2.545 + ///
2.546 + /// \sa MinCostArborescence
2.547 + template <typename Graph, typename CostMap, typename ArborescenceMap>
2.548 + typename CostMap::Value minCostArborescence(const Graph& graph,
2.549 + const CostMap& cost,
2.550 + typename Graph::Node source,
2.551 + ArborescenceMap& arborescence) {
2.552 + typename MinCostArborescence<Graph, CostMap>
2.553 + ::template DefArborescenceMap<ArborescenceMap>
2.554 + ::Create mca(graph, cost);
2.555 + mca.arborescenceMap(arborescence);
2.556 + mca.run(source);
2.557 + return mca.arborescenceCost();
2.558 + }
2.559 +
2.560 +}
2.561 +
2.562 +#endif
2.563 +
2.564 +// Hilbert - Huang