1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/gomory_hu.h Thu Nov 05 15:50:01 2009 +0100
1.3 @@ -0,0 +1,570 @@
1.4 +/* -*- C++ -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library
1.7 + *
1.8 + * Copyright (C) 2003-2008
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_GOMORY_HU_TREE_H
1.23 +#define LEMON_GOMORY_HU_TREE_H
1.24 +
1.25 +#include <limits>
1.26 +
1.27 +#include <lemon/core.h>
1.28 +#include <lemon/preflow.h>
1.29 +#include <lemon/concept_check.h>
1.30 +#include <lemon/concepts/maps.h>
1.31 +
1.32 +/// \ingroup min_cut
1.33 +/// \file
1.34 +/// \brief Gomory-Hu cut tree in graphs.
1.35 +
1.36 +namespace lemon {
1.37 +
1.38 + /// \ingroup min_cut
1.39 + ///
1.40 + /// \brief Gomory-Hu cut tree algorithm
1.41 + ///
1.42 + /// The Gomory-Hu tree is a tree on the node set of a given graph, but it
1.43 + /// may contain edges which are not in the original graph. It has the
1.44 + /// property that the minimum capacity edge of the path between two nodes
1.45 + /// in this tree has the same weight as the minimum cut in the graph
1.46 + /// between these nodes. Moreover the components obtained by removing
1.47 + /// this edge from the tree determine the corresponding minimum cut.
1.48 + /// Therefore once this tree is computed, the minimum cut between any pair
1.49 + /// of nodes can easily be obtained.
1.50 + ///
1.51 + /// The algorithm calculates \e n-1 distinct minimum cuts (currently with
1.52 + /// the \ref Preflow algorithm), thus it has \f$O(n^3\sqrt{e})\f$ overall
1.53 + /// time complexity. It calculates a rooted Gomory-Hu tree.
1.54 + /// The structure of the tree and the edge weights can be
1.55 + /// obtained using \c predNode(), \c predValue() and \c rootDist().
1.56 + /// The functions \c minCutMap() and \c minCutValue() calculate
1.57 + /// the minimum cut and the minimum cut value between any two nodes
1.58 + /// in the graph. You can also list (iterate on) the nodes and the
1.59 + /// edges of the cuts using \c MinCutNodeIt and \c MinCutEdgeIt.
1.60 + ///
1.61 + /// \tparam GR The type of the undirected graph the algorithm runs on.
1.62 + /// \tparam CAP The type of the edge map containing the capacities.
1.63 + /// The default map type is \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
1.64 +#ifdef DOXYGEN
1.65 + template <typename GR,
1.66 + typename CAP>
1.67 +#else
1.68 + template <typename GR,
1.69 + typename CAP = typename GR::template EdgeMap<int> >
1.70 +#endif
1.71 + class GomoryHu {
1.72 + public:
1.73 +
1.74 + /// The graph type of the algorithm
1.75 + typedef GR Graph;
1.76 + /// The capacity map type of the algorithm
1.77 + typedef CAP Capacity;
1.78 + /// The value type of capacities
1.79 + typedef typename Capacity::Value Value;
1.80 +
1.81 + private:
1.82 +
1.83 + TEMPLATE_GRAPH_TYPEDEFS(Graph);
1.84 +
1.85 + const Graph& _graph;
1.86 + const Capacity& _capacity;
1.87 +
1.88 + Node _root;
1.89 + typename Graph::template NodeMap<Node>* _pred;
1.90 + typename Graph::template NodeMap<Value>* _weight;
1.91 + typename Graph::template NodeMap<int>* _order;
1.92 +
1.93 + void createStructures() {
1.94 + if (!_pred) {
1.95 + _pred = new typename Graph::template NodeMap<Node>(_graph);
1.96 + }
1.97 + if (!_weight) {
1.98 + _weight = new typename Graph::template NodeMap<Value>(_graph);
1.99 + }
1.100 + if (!_order) {
1.101 + _order = new typename Graph::template NodeMap<int>(_graph);
1.102 + }
1.103 + }
1.104 +
1.105 + void destroyStructures() {
1.106 + if (_pred) {
1.107 + delete _pred;
1.108 + }
1.109 + if (_weight) {
1.110 + delete _weight;
1.111 + }
1.112 + if (_order) {
1.113 + delete _order;
1.114 + }
1.115 + }
1.116 +
1.117 + public:
1.118 +
1.119 + /// \brief Constructor
1.120 + ///
1.121 + /// Constructor.
1.122 + /// \param graph The undirected graph the algorithm runs on.
1.123 + /// \param capacity The edge capacity map.
1.124 + GomoryHu(const Graph& graph, const Capacity& capacity)
1.125 + : _graph(graph), _capacity(capacity),
1.126 + _pred(0), _weight(0), _order(0)
1.127 + {
1.128 + checkConcept<concepts::ReadMap<Edge, Value>, Capacity>();
1.129 + }
1.130 +
1.131 +
1.132 + /// \brief Destructor
1.133 + ///
1.134 + /// Destructor.
1.135 + ~GomoryHu() {
1.136 + destroyStructures();
1.137 + }
1.138 +
1.139 + private:
1.140 +
1.141 + // Initialize the internal data structures
1.142 + void init() {
1.143 + createStructures();
1.144 +
1.145 + _root = NodeIt(_graph);
1.146 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.147 + (*_pred)[n] = _root;
1.148 + (*_order)[n] = -1;
1.149 + }
1.150 + (*_pred)[_root] = INVALID;
1.151 + (*_weight)[_root] = std::numeric_limits<Value>::max();
1.152 + }
1.153 +
1.154 +
1.155 + // Start the algorithm
1.156 + void start() {
1.157 + Preflow<Graph, Capacity> fa(_graph, _capacity, _root, INVALID);
1.158 +
1.159 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.160 + if (n == _root) continue;
1.161 +
1.162 + Node pn = (*_pred)[n];
1.163 + fa.source(n);
1.164 + fa.target(pn);
1.165 +
1.166 + fa.runMinCut();
1.167 +
1.168 + (*_weight)[n] = fa.flowValue();
1.169 +
1.170 + for (NodeIt nn(_graph); nn != INVALID; ++nn) {
1.171 + if (nn != n && fa.minCut(nn) && (*_pred)[nn] == pn) {
1.172 + (*_pred)[nn] = n;
1.173 + }
1.174 + }
1.175 + if ((*_pred)[pn] != INVALID && fa.minCut((*_pred)[pn])) {
1.176 + (*_pred)[n] = (*_pred)[pn];
1.177 + (*_pred)[pn] = n;
1.178 + (*_weight)[n] = (*_weight)[pn];
1.179 + (*_weight)[pn] = fa.flowValue();
1.180 + }
1.181 + }
1.182 +
1.183 + (*_order)[_root] = 0;
1.184 + int index = 1;
1.185 +
1.186 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.187 + std::vector<Node> st;
1.188 + Node nn = n;
1.189 + while ((*_order)[nn] == -1) {
1.190 + st.push_back(nn);
1.191 + nn = (*_pred)[nn];
1.192 + }
1.193 + while (!st.empty()) {
1.194 + (*_order)[st.back()] = index++;
1.195 + st.pop_back();
1.196 + }
1.197 + }
1.198 + }
1.199 +
1.200 + public:
1.201 +
1.202 + ///\name Execution Control
1.203 +
1.204 + ///@{
1.205 +
1.206 + /// \brief Run the Gomory-Hu algorithm.
1.207 + ///
1.208 + /// This function runs the Gomory-Hu algorithm.
1.209 + void run() {
1.210 + init();
1.211 + start();
1.212 + }
1.213 +
1.214 + /// @}
1.215 +
1.216 + ///\name Query Functions
1.217 + ///The results of the algorithm can be obtained using these
1.218 + ///functions.\n
1.219 + ///\ref run() should be called before using them.\n
1.220 + ///See also \ref MinCutNodeIt and \ref MinCutEdgeIt.
1.221 +
1.222 + ///@{
1.223 +
1.224 + /// \brief Return the predecessor node in the Gomory-Hu tree.
1.225 + ///
1.226 + /// This function returns the predecessor node of the given node
1.227 + /// in the Gomory-Hu tree.
1.228 + /// If \c node is the root of the tree, then it returns \c INVALID.
1.229 + ///
1.230 + /// \pre \ref run() must be called before using this function.
1.231 + Node predNode(const Node& node) const {
1.232 + return (*_pred)[node];
1.233 + }
1.234 +
1.235 + /// \brief Return the weight of the predecessor edge in the
1.236 + /// Gomory-Hu tree.
1.237 + ///
1.238 + /// This function returns the weight of the predecessor edge of the
1.239 + /// given node in the Gomory-Hu tree.
1.240 + /// If \c node is the root of the tree, the result is undefined.
1.241 + ///
1.242 + /// \pre \ref run() must be called before using this function.
1.243 + Value predValue(const Node& node) const {
1.244 + return (*_weight)[node];
1.245 + }
1.246 +
1.247 + /// \brief Return the distance from the root node in the Gomory-Hu tree.
1.248 + ///
1.249 + /// This function returns the distance of the given node from the root
1.250 + /// node in the Gomory-Hu tree.
1.251 + ///
1.252 + /// \pre \ref run() must be called before using this function.
1.253 + int rootDist(const Node& node) const {
1.254 + return (*_order)[node];
1.255 + }
1.256 +
1.257 + /// \brief Return the minimum cut value between two nodes
1.258 + ///
1.259 + /// This function returns the minimum cut value between the nodes
1.260 + /// \c s and \c t.
1.261 + /// It finds the nearest common ancestor of the given nodes in the
1.262 + /// Gomory-Hu tree and calculates the minimum weight edge on the
1.263 + /// paths to the ancestor.
1.264 + ///
1.265 + /// \pre \ref run() must be called before using this function.
1.266 + Value minCutValue(const Node& s, const Node& t) const {
1.267 + Node sn = s, tn = t;
1.268 + Value value = std::numeric_limits<Value>::max();
1.269 +
1.270 + while (sn != tn) {
1.271 + if ((*_order)[sn] < (*_order)[tn]) {
1.272 + if ((*_weight)[tn] <= value) value = (*_weight)[tn];
1.273 + tn = (*_pred)[tn];
1.274 + } else {
1.275 + if ((*_weight)[sn] <= value) value = (*_weight)[sn];
1.276 + sn = (*_pred)[sn];
1.277 + }
1.278 + }
1.279 + return value;
1.280 + }
1.281 +
1.282 + /// \brief Return the minimum cut between two nodes
1.283 + ///
1.284 + /// This function returns the minimum cut between the nodes \c s and \c t
1.285 + /// in the \c cutMap parameter by setting the nodes in the component of
1.286 + /// \c s to \c true and the other nodes to \c false.
1.287 + ///
1.288 + /// For higher level interfaces see MinCutNodeIt and MinCutEdgeIt.
1.289 + ///
1.290 + /// \param s The base node.
1.291 + /// \param t The node you want to separate from node \c s.
1.292 + /// \param cutMap The cut will be returned in this map.
1.293 + /// It must be a \c bool (or convertible) \ref concepts::ReadWriteMap
1.294 + /// "ReadWriteMap" on the graph nodes.
1.295 + ///
1.296 + /// \return The value of the minimum cut between \c s and \c t.
1.297 + ///
1.298 + /// \pre \ref run() must be called before using this function.
1.299 + template <typename CutMap>
1.300 + Value minCutMap(const Node& s, ///<
1.301 + const Node& t,
1.302 + ///<
1.303 + CutMap& cutMap
1.304 + ///<
1.305 + ) const {
1.306 + Node sn = s, tn = t;
1.307 + bool s_root=false;
1.308 + Node rn = INVALID;
1.309 + Value value = std::numeric_limits<Value>::max();
1.310 +
1.311 + while (sn != tn) {
1.312 + if ((*_order)[sn] < (*_order)[tn]) {
1.313 + if ((*_weight)[tn] <= value) {
1.314 + rn = tn;
1.315 + s_root = false;
1.316 + value = (*_weight)[tn];
1.317 + }
1.318 + tn = (*_pred)[tn];
1.319 + } else {
1.320 + if ((*_weight)[sn] <= value) {
1.321 + rn = sn;
1.322 + s_root = true;
1.323 + value = (*_weight)[sn];
1.324 + }
1.325 + sn = (*_pred)[sn];
1.326 + }
1.327 + }
1.328 +
1.329 + typename Graph::template NodeMap<bool> reached(_graph, false);
1.330 + reached[_root] = true;
1.331 + cutMap.set(_root, !s_root);
1.332 + reached[rn] = true;
1.333 + cutMap.set(rn, s_root);
1.334 +
1.335 + std::vector<Node> st;
1.336 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.337 + st.clear();
1.338 + Node nn = n;
1.339 + while (!reached[nn]) {
1.340 + st.push_back(nn);
1.341 + nn = (*_pred)[nn];
1.342 + }
1.343 + while (!st.empty()) {
1.344 + cutMap.set(st.back(), cutMap[nn]);
1.345 + st.pop_back();
1.346 + }
1.347 + }
1.348 +
1.349 + return value;
1.350 + }
1.351 +
1.352 + ///@}
1.353 +
1.354 + friend class MinCutNodeIt;
1.355 +
1.356 + /// Iterate on the nodes of a minimum cut
1.357 +
1.358 + /// This iterator class lists the nodes of a minimum cut found by
1.359 + /// GomoryHu. Before using it, you must allocate a GomoryHu class
1.360 + /// and call its \ref GomoryHu::run() "run()" method.
1.361 + ///
1.362 + /// This example counts the nodes in the minimum cut separating \c s from
1.363 + /// \c t.
1.364 + /// \code
1.365 + /// GomoryHu<Graph> gom(g, capacities);
1.366 + /// gom.run();
1.367 + /// int cnt=0;
1.368 + /// for(GomoryHu<Graph>::MinCutNodeIt n(gom,s,t); n!=INVALID; ++n) ++cnt;
1.369 + /// \endcode
1.370 + class MinCutNodeIt
1.371 + {
1.372 + bool _side;
1.373 + typename Graph::NodeIt _node_it;
1.374 + typename Graph::template NodeMap<bool> _cut;
1.375 + public:
1.376 + /// Constructor
1.377 +
1.378 + /// Constructor.
1.379 + ///
1.380 + MinCutNodeIt(GomoryHu const &gomory,
1.381 + ///< The GomoryHu class. You must call its
1.382 + /// run() method
1.383 + /// before initializing this iterator.
1.384 + const Node& s, ///< The base node.
1.385 + const Node& t,
1.386 + ///< The node you want to separate from node \c s.
1.387 + bool side=true
1.388 + ///< If it is \c true (default) then the iterator lists
1.389 + /// the nodes of the component containing \c s,
1.390 + /// otherwise it lists the other component.
1.391 + /// \note As the minimum cut is not always unique,
1.392 + /// \code
1.393 + /// MinCutNodeIt(gomory, s, t, true);
1.394 + /// \endcode
1.395 + /// and
1.396 + /// \code
1.397 + /// MinCutNodeIt(gomory, t, s, false);
1.398 + /// \endcode
1.399 + /// does not necessarily give the same set of nodes.
1.400 + /// However it is ensured that
1.401 + /// \code
1.402 + /// MinCutNodeIt(gomory, s, t, true);
1.403 + /// \endcode
1.404 + /// and
1.405 + /// \code
1.406 + /// MinCutNodeIt(gomory, s, t, false);
1.407 + /// \endcode
1.408 + /// together list each node exactly once.
1.409 + )
1.410 + : _side(side), _cut(gomory._graph)
1.411 + {
1.412 + gomory.minCutMap(s,t,_cut);
1.413 + for(_node_it=typename Graph::NodeIt(gomory._graph);
1.414 + _node_it!=INVALID && _cut[_node_it]!=_side;
1.415 + ++_node_it) {}
1.416 + }
1.417 + /// Conversion to \c Node
1.418 +
1.419 + /// Conversion to \c Node.
1.420 + ///
1.421 + operator typename Graph::Node() const
1.422 + {
1.423 + return _node_it;
1.424 + }
1.425 + bool operator==(Invalid) { return _node_it==INVALID; }
1.426 + bool operator!=(Invalid) { return _node_it!=INVALID; }
1.427 + /// Next node
1.428 +
1.429 + /// Next node.
1.430 + ///
1.431 + MinCutNodeIt &operator++()
1.432 + {
1.433 + for(++_node_it;_node_it!=INVALID&&_cut[_node_it]!=_side;++_node_it) {}
1.434 + return *this;
1.435 + }
1.436 + /// Postfix incrementation
1.437 +
1.438 + /// Postfix incrementation.
1.439 + ///
1.440 + /// \warning This incrementation
1.441 + /// returns a \c Node, not a \c MinCutNodeIt, as one may
1.442 + /// expect.
1.443 + typename Graph::Node operator++(int)
1.444 + {
1.445 + typename Graph::Node n=*this;
1.446 + ++(*this);
1.447 + return n;
1.448 + }
1.449 + };
1.450 +
1.451 + friend class MinCutEdgeIt;
1.452 +
1.453 + /// Iterate on the edges of a minimum cut
1.454 +
1.455 + /// This iterator class lists the edges of a minimum cut found by
1.456 + /// GomoryHu. Before using it, you must allocate a GomoryHu class
1.457 + /// and call its \ref GomoryHu::run() "run()" method.
1.458 + ///
1.459 + /// This example computes the value of the minimum cut separating \c s from
1.460 + /// \c t.
1.461 + /// \code
1.462 + /// GomoryHu<Graph> gom(g, capacities);
1.463 + /// gom.run();
1.464 + /// int value=0;
1.465 + /// for(GomoryHu<Graph>::MinCutEdgeIt e(gom,s,t); e!=INVALID; ++e)
1.466 + /// value+=capacities[e];
1.467 + /// \endcode
1.468 + /// The result will be the same as the value returned by
1.469 + /// \ref GomoryHu::minCutValue() "gom.minCutValue(s,t)".
1.470 + class MinCutEdgeIt
1.471 + {
1.472 + bool _side;
1.473 + const Graph &_graph;
1.474 + typename Graph::NodeIt _node_it;
1.475 + typename Graph::OutArcIt _arc_it;
1.476 + typename Graph::template NodeMap<bool> _cut;
1.477 + void step()
1.478 + {
1.479 + ++_arc_it;
1.480 + while(_node_it!=INVALID && _arc_it==INVALID)
1.481 + {
1.482 + for(++_node_it;_node_it!=INVALID&&!_cut[_node_it];++_node_it) {}
1.483 + if(_node_it!=INVALID)
1.484 + _arc_it=typename Graph::OutArcIt(_graph,_node_it);
1.485 + }
1.486 + }
1.487 +
1.488 + public:
1.489 + /// Constructor
1.490 +
1.491 + /// Constructor.
1.492 + ///
1.493 + MinCutEdgeIt(GomoryHu const &gomory,
1.494 + ///< The GomoryHu class. You must call its
1.495 + /// run() method
1.496 + /// before initializing this iterator.
1.497 + const Node& s, ///< The base node.
1.498 + const Node& t,
1.499 + ///< The node you want to separate from node \c s.
1.500 + bool side=true
1.501 + ///< If it is \c true (default) then the listed arcs
1.502 + /// will be oriented from the
1.503 + /// nodes of the component containing \c s,
1.504 + /// otherwise they will be oriented in the opposite
1.505 + /// direction.
1.506 + )
1.507 + : _graph(gomory._graph), _cut(_graph)
1.508 + {
1.509 + gomory.minCutMap(s,t,_cut);
1.510 + if(!side)
1.511 + for(typename Graph::NodeIt n(_graph);n!=INVALID;++n)
1.512 + _cut[n]=!_cut[n];
1.513 +
1.514 + for(_node_it=typename Graph::NodeIt(_graph);
1.515 + _node_it!=INVALID && !_cut[_node_it];
1.516 + ++_node_it) {}
1.517 + _arc_it = _node_it!=INVALID ?
1.518 + typename Graph::OutArcIt(_graph,_node_it) : INVALID;
1.519 + while(_node_it!=INVALID && _arc_it == INVALID)
1.520 + {
1.521 + for(++_node_it; _node_it!=INVALID&&!_cut[_node_it]; ++_node_it) {}
1.522 + if(_node_it!=INVALID)
1.523 + _arc_it= typename Graph::OutArcIt(_graph,_node_it);
1.524 + }
1.525 + while(_arc_it!=INVALID && _cut[_graph.target(_arc_it)]) step();
1.526 + }
1.527 + /// Conversion to \c Arc
1.528 +
1.529 + /// Conversion to \c Arc.
1.530 + ///
1.531 + operator typename Graph::Arc() const
1.532 + {
1.533 + return _arc_it;
1.534 + }
1.535 + /// Conversion to \c Edge
1.536 +
1.537 + /// Conversion to \c Edge.
1.538 + ///
1.539 + operator typename Graph::Edge() const
1.540 + {
1.541 + return _arc_it;
1.542 + }
1.543 + bool operator==(Invalid) { return _node_it==INVALID; }
1.544 + bool operator!=(Invalid) { return _node_it!=INVALID; }
1.545 + /// Next edge
1.546 +
1.547 + /// Next edge.
1.548 + ///
1.549 + MinCutEdgeIt &operator++()
1.550 + {
1.551 + step();
1.552 + while(_arc_it!=INVALID && _cut[_graph.target(_arc_it)]) step();
1.553 + return *this;
1.554 + }
1.555 + /// Postfix incrementation
1.556 +
1.557 + /// Postfix incrementation.
1.558 + ///
1.559 + /// \warning This incrementation
1.560 + /// returns an \c Arc, not a \c MinCutEdgeIt, as one may expect.
1.561 + typename Graph::Arc operator++(int)
1.562 + {
1.563 + typename Graph::Arc e=*this;
1.564 + ++(*this);
1.565 + return e;
1.566 + }
1.567 + };
1.568 +
1.569 + };
1.570 +
1.571 +}
1.572 +
1.573 +#endif