lemon/gomory_hu_tree.h
changeset 543 924887566bf2
child 544 ccd2d3a3001e
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lemon/gomory_hu_tree.h	Fri Feb 20 17:17:17 2009 +0100
     1.3 @@ -0,0 +1,298 @@
     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/preflow.h>
    1.28 +#include <lemon/concept_check.h>
    1.29 +#include <lemon/concepts/maps.h>
    1.30 +
    1.31 +/// \ingroup min_cut
    1.32 +/// \file 
    1.33 +/// \brief Gomory-Hu cut tree in graphs.
    1.34 +
    1.35 +namespace lemon {
    1.36 +
    1.37 +  /// \ingroup min_cut
    1.38 +  ///
    1.39 +  /// \brief Gomory-Hu cut tree algorithm
    1.40 +  ///
    1.41 +  /// The Gomory-Hu tree is a tree on the nodeset of the digraph, but it
    1.42 +  /// may contain arcs which are not in the original digraph. It helps
    1.43 +  /// to calculate the minimum cut between all pairs of nodes, because
    1.44 +  /// the minimum capacity arc on the tree path between two nodes has
    1.45 +  /// the same weight as the minimum cut in the digraph between these
    1.46 +  /// nodes. Moreover this arc separates the nodes to two parts which
    1.47 +  /// determine this minimum cut.
    1.48 +  /// 
    1.49 +  /// The algorithm calculates \e n-1 distinict minimum cuts with
    1.50 +  /// preflow algorithm, therefore the algorithm has
    1.51 +  /// \f$(O(n^3\sqrt{e})\f$ overall time complexity. It calculates a
    1.52 +  /// rooted Gomory-Hu tree, the structure of the tree and the weights
    1.53 +  /// can be obtained with \c predNode() and \c predValue()
    1.54 +  /// functions. The \c minCutValue() and \c minCutMap() calculates
    1.55 +  /// the minimum cut and the minimum cut value between any two node
    1.56 +  /// in the digraph.
    1.57 +  template <typename _Graph, 
    1.58 +	    typename _Capacity = typename _Graph::template EdgeMap<int> >
    1.59 +  class GomoryHuTree {
    1.60 +  public:
    1.61 +
    1.62 +    /// The graph type
    1.63 +    typedef _Graph Graph;
    1.64 +    /// The capacity on edges
    1.65 +    typedef _Capacity Capacity;
    1.66 +    /// The value type of capacities
    1.67 +    typedef typename Capacity::Value Value;
    1.68 +    
    1.69 +  private:
    1.70 +
    1.71 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
    1.72 +
    1.73 +    const Graph& _graph;
    1.74 +    const Capacity& _capacity;
    1.75 +
    1.76 +    Node _root;
    1.77 +    typename Graph::template NodeMap<Node>* _pred;
    1.78 +    typename Graph::template NodeMap<Value>* _weight;
    1.79 +    typename Graph::template NodeMap<int>* _order;
    1.80 +
    1.81 +    void createStructures() {
    1.82 +      if (!_pred) {
    1.83 +	_pred = new typename Graph::template NodeMap<Node>(_graph);
    1.84 +      }
    1.85 +      if (!_weight) {
    1.86 +	_weight = new typename Graph::template NodeMap<Value>(_graph);
    1.87 +      }
    1.88 +      if (!_order) {
    1.89 +	_order = new typename Graph::template NodeMap<int>(_graph);
    1.90 +      }
    1.91 +    }
    1.92 +
    1.93 +    void destroyStructures() {
    1.94 +      if (_pred) {
    1.95 +	delete _pred;
    1.96 +      }
    1.97 +      if (_weight) {
    1.98 +	delete _weight;
    1.99 +      }
   1.100 +      if (_order) {
   1.101 +	delete _order;
   1.102 +      }
   1.103 +    }
   1.104 +  
   1.105 +  public:
   1.106 +
   1.107 +    /// \brief Constructor
   1.108 +    ///
   1.109 +    /// Constructor
   1.110 +    /// \param graph The graph type.
   1.111 +    /// \param capacity The capacity map.
   1.112 +    GomoryHuTree(const Graph& graph, const Capacity& capacity) 
   1.113 +      : _graph(graph), _capacity(capacity),
   1.114 +	_pred(0), _weight(0), _order(0) 
   1.115 +    {
   1.116 +      checkConcept<concepts::ReadMap<Edge, Value>, Capacity>();
   1.117 +    }
   1.118 +
   1.119 +
   1.120 +    /// \brief Destructor
   1.121 +    ///
   1.122 +    /// Destructor
   1.123 +    ~GomoryHuTree() {
   1.124 +      destroyStructures();
   1.125 +    }
   1.126 +
   1.127 +    /// \brief Initializes the internal data structures.
   1.128 +    ///
   1.129 +    /// Initializes the internal data structures.
   1.130 +    ///
   1.131 +    void init() {
   1.132 +      createStructures();
   1.133 +
   1.134 +      _root = NodeIt(_graph);
   1.135 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.136 +	_pred->set(n, _root);
   1.137 +	_order->set(n, -1);
   1.138 +      }
   1.139 +      _pred->set(_root, INVALID);
   1.140 +      _weight->set(_root, std::numeric_limits<Value>::max()); 
   1.141 +    }
   1.142 +
   1.143 +
   1.144 +    /// \brief Starts the algorithm
   1.145 +    ///
   1.146 +    /// Starts the algorithm.
   1.147 +    void start() {
   1.148 +      Preflow<Graph, Capacity> fa(_graph, _capacity, _root, INVALID);
   1.149 +
   1.150 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.151 +	if (n == _root) continue;
   1.152 +
   1.153 +	Node pn = (*_pred)[n];
   1.154 +	fa.source(n);
   1.155 +	fa.target(pn);
   1.156 +
   1.157 +	fa.runMinCut();
   1.158 +
   1.159 +	_weight->set(n, fa.flowValue());
   1.160 +
   1.161 +	for (NodeIt nn(_graph); nn != INVALID; ++nn) {
   1.162 +	  if (nn != n && fa.minCut(nn) && (*_pred)[nn] == pn) {
   1.163 +	    _pred->set(nn, n);
   1.164 +	  }
   1.165 +	}
   1.166 +	if ((*_pred)[pn] != INVALID && fa.minCut((*_pred)[pn])) {
   1.167 +	  _pred->set(n, (*_pred)[pn]);
   1.168 +	  _pred->set(pn, n);
   1.169 +	  _weight->set(n, (*_weight)[pn]);
   1.170 +	  _weight->set(pn, fa.flowValue());	
   1.171 +	}
   1.172 +      }
   1.173 +
   1.174 +      _order->set(_root, 0);
   1.175 +      int index = 1;
   1.176 +
   1.177 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.178 +	std::vector<Node> st;
   1.179 +	Node nn = n;
   1.180 +	while ((*_order)[nn] == -1) {
   1.181 +	  st.push_back(nn);
   1.182 +	  nn = (*_pred)[nn];
   1.183 +	}
   1.184 +	while (!st.empty()) {
   1.185 +	  _order->set(st.back(), index++);
   1.186 +	  st.pop_back();
   1.187 +	}
   1.188 +      }
   1.189 +    }
   1.190 +
   1.191 +    /// \brief Runs the Gomory-Hu algorithm.  
   1.192 +    ///
   1.193 +    /// Runs the Gomory-Hu algorithm.
   1.194 +    /// \note gh.run() is just a shortcut of the following code.
   1.195 +    /// \code
   1.196 +    ///   ght.init();
   1.197 +    ///   ght.start();
   1.198 +    /// \endcode
   1.199 +    void run() {
   1.200 +      init();
   1.201 +      start();
   1.202 +    }
   1.203 +
   1.204 +    /// \brief Returns the predecessor node in the Gomory-Hu tree.
   1.205 +    ///
   1.206 +    /// Returns the predecessor node in the Gomory-Hu tree. If the node is
   1.207 +    /// the root of the Gomory-Hu tree, then it returns \c INVALID.
   1.208 +    Node predNode(const Node& node) {
   1.209 +      return (*_pred)[node];
   1.210 +    }
   1.211 +
   1.212 +    /// \brief Returns the weight of the predecessor arc in the
   1.213 +    /// Gomory-Hu tree.
   1.214 +    ///
   1.215 +    /// Returns the weight of the predecessor arc in the Gomory-Hu
   1.216 +    /// tree.  If the node is the root of the Gomory-Hu tree, the
   1.217 +    /// result is undefined.
   1.218 +    Value predValue(const Node& node) {
   1.219 +      return (*_weight)[node];
   1.220 +    }
   1.221 +
   1.222 +    /// \brief Returns the minimum cut value between two nodes
   1.223 +    ///
   1.224 +    /// Returns the minimum cut value between two nodes. The
   1.225 +    /// algorithm finds the nearest common ancestor in the Gomory-Hu
   1.226 +    /// tree and calculates the minimum weight arc on the paths to
   1.227 +    /// the ancestor.
   1.228 +    Value minCutValue(const Node& s, const Node& t) const {
   1.229 +      Node sn = s, tn = t;
   1.230 +      Value value = std::numeric_limits<Value>::max();
   1.231 +      
   1.232 +      while (sn != tn) {
   1.233 +	if ((*_order)[sn] < (*_order)[tn]) {
   1.234 +	  if ((*_weight)[tn] < value) value = (*_weight)[tn];
   1.235 +	  tn = (*_pred)[tn];
   1.236 +	} else {
   1.237 +	  if ((*_weight)[sn] < value) value = (*_weight)[sn];
   1.238 +	  sn = (*_pred)[sn];
   1.239 +	}
   1.240 +      }
   1.241 +      return value;
   1.242 +    }
   1.243 +
   1.244 +    /// \brief Returns the minimum cut between two nodes
   1.245 +    ///
   1.246 +    /// Returns the minimum cut value between two nodes. The
   1.247 +    /// algorithm finds the nearest common ancestor in the Gomory-Hu
   1.248 +    /// tree and calculates the minimum weight arc on the paths to
   1.249 +    /// the ancestor. Then it sets all nodes to the cut determined by
   1.250 +    /// this arc. The \c cutMap should be \ref concepts::ReadWriteMap
   1.251 +    /// "ReadWriteMap".
   1.252 +    template <typename CutMap>
   1.253 +    Value minCutMap(const Node& s, const Node& t, CutMap& cutMap) const {
   1.254 +      Node sn = s, tn = t;
   1.255 +
   1.256 +      Node rn = INVALID;
   1.257 +      Value value = std::numeric_limits<Value>::max();
   1.258 +      
   1.259 +      while (sn != tn) {
   1.260 +	if ((*_order)[sn] < (*_order)[tn]) {
   1.261 +	  if ((*_weight)[tn] < value) {
   1.262 +	    rn = tn;
   1.263 +	    value = (*_weight)[tn];
   1.264 +	  }
   1.265 +	  tn = (*_pred)[tn];
   1.266 +	} else {
   1.267 +	  if ((*_weight)[sn] < value) {
   1.268 +	    rn = sn;
   1.269 +	    value = (*_weight)[sn];
   1.270 +	  }
   1.271 +	  sn = (*_pred)[sn];
   1.272 +	}
   1.273 +      }
   1.274 +
   1.275 +      typename Graph::template NodeMap<bool> reached(_graph, false);
   1.276 +      reached.set(_root, true);
   1.277 +      cutMap.set(_root, false);
   1.278 +      reached.set(rn, true);
   1.279 +      cutMap.set(rn, true);
   1.280 +
   1.281 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.282 +	std::vector<Node> st;
   1.283 +	Node nn = n;
   1.284 +	while (!reached[nn]) {
   1.285 +	  st.push_back(nn);
   1.286 +	  nn = (*_pred)[nn];
   1.287 +	}
   1.288 +	while (!st.empty()) {
   1.289 +	  cutMap.set(st.back(), cutMap[nn]);
   1.290 +	  st.pop_back();
   1.291 +	}
   1.292 +      }
   1.293 +      
   1.294 +      return value;
   1.295 +    }
   1.296 +
   1.297 +  };
   1.298 +
   1.299 +}
   1.300 +
   1.301 +#endif