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