Porting Gomory-Hu algorithm (#66)
authorJanos Tapolcai <tapolcai@tmit.bme.hu>
Fri, 20 Feb 2009 17:17:17 +0100
changeset 531924887566bf2
parent 524 97674155c135
child 532 ccd2d3a3001e
Porting Gomory-Hu algorithm (#66)
lemon/Makefile.am
lemon/gomory_hu_tree.h
test/CMakeLists.txt
test/Makefile.am
test/gomory_hu_test.cc
     1.1 --- a/lemon/Makefile.am	Tue Feb 24 06:52:17 2009 +0000
     1.2 +++ b/lemon/Makefile.am	Fri Feb 20 17:17:17 2009 +0100
     1.3 @@ -68,6 +68,7 @@
     1.4  	lemon/euler.h \
     1.5  	lemon/full_graph.h \
     1.6  	lemon/glpk.h \
     1.7 +	lemon/gomory_hu_tree.h \
     1.8  	lemon/graph_to_eps.h \
     1.9  	lemon/grid_graph.h \
    1.10  	lemon/hypercube_graph.h \
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/lemon/gomory_hu_tree.h	Fri Feb 20 17:17:17 2009 +0100
     2.3 @@ -0,0 +1,298 @@
     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-2008
     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_GOMORY_HU_TREE_H
    2.23 +#define LEMON_GOMORY_HU_TREE_H
    2.24 +
    2.25 +#include <limits>
    2.26 +
    2.27 +#include <lemon/preflow.h>
    2.28 +#include <lemon/concept_check.h>
    2.29 +#include <lemon/concepts/maps.h>
    2.30 +
    2.31 +/// \ingroup min_cut
    2.32 +/// \file 
    2.33 +/// \brief Gomory-Hu cut tree in graphs.
    2.34 +
    2.35 +namespace lemon {
    2.36 +
    2.37 +  /// \ingroup min_cut
    2.38 +  ///
    2.39 +  /// \brief Gomory-Hu cut tree algorithm
    2.40 +  ///
    2.41 +  /// The Gomory-Hu tree is a tree on the nodeset of the digraph, but it
    2.42 +  /// may contain arcs which are not in the original digraph. It helps
    2.43 +  /// to calculate the minimum cut between all pairs of nodes, because
    2.44 +  /// the minimum capacity arc on the tree path between two nodes has
    2.45 +  /// the same weight as the minimum cut in the digraph between these
    2.46 +  /// nodes. Moreover this arc separates the nodes to two parts which
    2.47 +  /// determine this minimum cut.
    2.48 +  /// 
    2.49 +  /// The algorithm calculates \e n-1 distinict minimum cuts with
    2.50 +  /// preflow algorithm, therefore the algorithm has
    2.51 +  /// \f$(O(n^3\sqrt{e})\f$ overall time complexity. It calculates a
    2.52 +  /// rooted Gomory-Hu tree, the structure of the tree and the weights
    2.53 +  /// can be obtained with \c predNode() and \c predValue()
    2.54 +  /// functions. The \c minCutValue() and \c minCutMap() calculates
    2.55 +  /// the minimum cut and the minimum cut value between any two node
    2.56 +  /// in the digraph.
    2.57 +  template <typename _Graph, 
    2.58 +	    typename _Capacity = typename _Graph::template EdgeMap<int> >
    2.59 +  class GomoryHuTree {
    2.60 +  public:
    2.61 +
    2.62 +    /// The graph type
    2.63 +    typedef _Graph Graph;
    2.64 +    /// The capacity on edges
    2.65 +    typedef _Capacity Capacity;
    2.66 +    /// The value type of capacities
    2.67 +    typedef typename Capacity::Value Value;
    2.68 +    
    2.69 +  private:
    2.70 +
    2.71 +    TEMPLATE_GRAPH_TYPEDEFS(Graph);
    2.72 +
    2.73 +    const Graph& _graph;
    2.74 +    const Capacity& _capacity;
    2.75 +
    2.76 +    Node _root;
    2.77 +    typename Graph::template NodeMap<Node>* _pred;
    2.78 +    typename Graph::template NodeMap<Value>* _weight;
    2.79 +    typename Graph::template NodeMap<int>* _order;
    2.80 +
    2.81 +    void createStructures() {
    2.82 +      if (!_pred) {
    2.83 +	_pred = new typename Graph::template NodeMap<Node>(_graph);
    2.84 +      }
    2.85 +      if (!_weight) {
    2.86 +	_weight = new typename Graph::template NodeMap<Value>(_graph);
    2.87 +      }
    2.88 +      if (!_order) {
    2.89 +	_order = new typename Graph::template NodeMap<int>(_graph);
    2.90 +      }
    2.91 +    }
    2.92 +
    2.93 +    void destroyStructures() {
    2.94 +      if (_pred) {
    2.95 +	delete _pred;
    2.96 +      }
    2.97 +      if (_weight) {
    2.98 +	delete _weight;
    2.99 +      }
   2.100 +      if (_order) {
   2.101 +	delete _order;
   2.102 +      }
   2.103 +    }
   2.104 +  
   2.105 +  public:
   2.106 +
   2.107 +    /// \brief Constructor
   2.108 +    ///
   2.109 +    /// Constructor
   2.110 +    /// \param graph The graph type.
   2.111 +    /// \param capacity The capacity map.
   2.112 +    GomoryHuTree(const Graph& graph, const Capacity& capacity) 
   2.113 +      : _graph(graph), _capacity(capacity),
   2.114 +	_pred(0), _weight(0), _order(0) 
   2.115 +    {
   2.116 +      checkConcept<concepts::ReadMap<Edge, Value>, Capacity>();
   2.117 +    }
   2.118 +
   2.119 +
   2.120 +    /// \brief Destructor
   2.121 +    ///
   2.122 +    /// Destructor
   2.123 +    ~GomoryHuTree() {
   2.124 +      destroyStructures();
   2.125 +    }
   2.126 +
   2.127 +    /// \brief Initializes the internal data structures.
   2.128 +    ///
   2.129 +    /// Initializes the internal data structures.
   2.130 +    ///
   2.131 +    void init() {
   2.132 +      createStructures();
   2.133 +
   2.134 +      _root = NodeIt(_graph);
   2.135 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   2.136 +	_pred->set(n, _root);
   2.137 +	_order->set(n, -1);
   2.138 +      }
   2.139 +      _pred->set(_root, INVALID);
   2.140 +      _weight->set(_root, std::numeric_limits<Value>::max()); 
   2.141 +    }
   2.142 +
   2.143 +
   2.144 +    /// \brief Starts the algorithm
   2.145 +    ///
   2.146 +    /// Starts the algorithm.
   2.147 +    void start() {
   2.148 +      Preflow<Graph, Capacity> fa(_graph, _capacity, _root, INVALID);
   2.149 +
   2.150 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   2.151 +	if (n == _root) continue;
   2.152 +
   2.153 +	Node pn = (*_pred)[n];
   2.154 +	fa.source(n);
   2.155 +	fa.target(pn);
   2.156 +
   2.157 +	fa.runMinCut();
   2.158 +
   2.159 +	_weight->set(n, fa.flowValue());
   2.160 +
   2.161 +	for (NodeIt nn(_graph); nn != INVALID; ++nn) {
   2.162 +	  if (nn != n && fa.minCut(nn) && (*_pred)[nn] == pn) {
   2.163 +	    _pred->set(nn, n);
   2.164 +	  }
   2.165 +	}
   2.166 +	if ((*_pred)[pn] != INVALID && fa.minCut((*_pred)[pn])) {
   2.167 +	  _pred->set(n, (*_pred)[pn]);
   2.168 +	  _pred->set(pn, n);
   2.169 +	  _weight->set(n, (*_weight)[pn]);
   2.170 +	  _weight->set(pn, fa.flowValue());	
   2.171 +	}
   2.172 +      }
   2.173 +
   2.174 +      _order->set(_root, 0);
   2.175 +      int index = 1;
   2.176 +
   2.177 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   2.178 +	std::vector<Node> st;
   2.179 +	Node nn = n;
   2.180 +	while ((*_order)[nn] == -1) {
   2.181 +	  st.push_back(nn);
   2.182 +	  nn = (*_pred)[nn];
   2.183 +	}
   2.184 +	while (!st.empty()) {
   2.185 +	  _order->set(st.back(), index++);
   2.186 +	  st.pop_back();
   2.187 +	}
   2.188 +      }
   2.189 +    }
   2.190 +
   2.191 +    /// \brief Runs the Gomory-Hu algorithm.  
   2.192 +    ///
   2.193 +    /// Runs the Gomory-Hu algorithm.
   2.194 +    /// \note gh.run() is just a shortcut of the following code.
   2.195 +    /// \code
   2.196 +    ///   ght.init();
   2.197 +    ///   ght.start();
   2.198 +    /// \endcode
   2.199 +    void run() {
   2.200 +      init();
   2.201 +      start();
   2.202 +    }
   2.203 +
   2.204 +    /// \brief Returns the predecessor node in the Gomory-Hu tree.
   2.205 +    ///
   2.206 +    /// Returns the predecessor node in the Gomory-Hu tree. If the node is
   2.207 +    /// the root of the Gomory-Hu tree, then it returns \c INVALID.
   2.208 +    Node predNode(const Node& node) {
   2.209 +      return (*_pred)[node];
   2.210 +    }
   2.211 +
   2.212 +    /// \brief Returns the weight of the predecessor arc in the
   2.213 +    /// Gomory-Hu tree.
   2.214 +    ///
   2.215 +    /// Returns the weight of the predecessor arc in the Gomory-Hu
   2.216 +    /// tree.  If the node is the root of the Gomory-Hu tree, the
   2.217 +    /// result is undefined.
   2.218 +    Value predValue(const Node& node) {
   2.219 +      return (*_weight)[node];
   2.220 +    }
   2.221 +
   2.222 +    /// \brief Returns the minimum cut value between two nodes
   2.223 +    ///
   2.224 +    /// Returns the minimum cut value between two nodes. The
   2.225 +    /// algorithm finds the nearest common ancestor in the Gomory-Hu
   2.226 +    /// tree and calculates the minimum weight arc on the paths to
   2.227 +    /// the ancestor.
   2.228 +    Value minCutValue(const Node& s, const Node& t) const {
   2.229 +      Node sn = s, tn = t;
   2.230 +      Value value = std::numeric_limits<Value>::max();
   2.231 +      
   2.232 +      while (sn != tn) {
   2.233 +	if ((*_order)[sn] < (*_order)[tn]) {
   2.234 +	  if ((*_weight)[tn] < value) value = (*_weight)[tn];
   2.235 +	  tn = (*_pred)[tn];
   2.236 +	} else {
   2.237 +	  if ((*_weight)[sn] < value) value = (*_weight)[sn];
   2.238 +	  sn = (*_pred)[sn];
   2.239 +	}
   2.240 +      }
   2.241 +      return value;
   2.242 +    }
   2.243 +
   2.244 +    /// \brief Returns the minimum cut between two nodes
   2.245 +    ///
   2.246 +    /// Returns the minimum cut value between two nodes. The
   2.247 +    /// algorithm finds the nearest common ancestor in the Gomory-Hu
   2.248 +    /// tree and calculates the minimum weight arc on the paths to
   2.249 +    /// the ancestor. Then it sets all nodes to the cut determined by
   2.250 +    /// this arc. The \c cutMap should be \ref concepts::ReadWriteMap
   2.251 +    /// "ReadWriteMap".
   2.252 +    template <typename CutMap>
   2.253 +    Value minCutMap(const Node& s, const Node& t, CutMap& cutMap) const {
   2.254 +      Node sn = s, tn = t;
   2.255 +
   2.256 +      Node rn = INVALID;
   2.257 +      Value value = std::numeric_limits<Value>::max();
   2.258 +      
   2.259 +      while (sn != tn) {
   2.260 +	if ((*_order)[sn] < (*_order)[tn]) {
   2.261 +	  if ((*_weight)[tn] < value) {
   2.262 +	    rn = tn;
   2.263 +	    value = (*_weight)[tn];
   2.264 +	  }
   2.265 +	  tn = (*_pred)[tn];
   2.266 +	} else {
   2.267 +	  if ((*_weight)[sn] < value) {
   2.268 +	    rn = sn;
   2.269 +	    value = (*_weight)[sn];
   2.270 +	  }
   2.271 +	  sn = (*_pred)[sn];
   2.272 +	}
   2.273 +      }
   2.274 +
   2.275 +      typename Graph::template NodeMap<bool> reached(_graph, false);
   2.276 +      reached.set(_root, true);
   2.277 +      cutMap.set(_root, false);
   2.278 +      reached.set(rn, true);
   2.279 +      cutMap.set(rn, true);
   2.280 +
   2.281 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   2.282 +	std::vector<Node> st;
   2.283 +	Node nn = n;
   2.284 +	while (!reached[nn]) {
   2.285 +	  st.push_back(nn);
   2.286 +	  nn = (*_pred)[nn];
   2.287 +	}
   2.288 +	while (!st.empty()) {
   2.289 +	  cutMap.set(st.back(), cutMap[nn]);
   2.290 +	  st.pop_back();
   2.291 +	}
   2.292 +      }
   2.293 +      
   2.294 +      return value;
   2.295 +    }
   2.296 +
   2.297 +  };
   2.298 +
   2.299 +}
   2.300 +
   2.301 +#endif
     3.1 --- a/test/CMakeLists.txt	Tue Feb 24 06:52:17 2009 +0000
     3.2 +++ b/test/CMakeLists.txt	Fri Feb 20 17:17:17 2009 +0100
     3.3 @@ -21,6 +21,7 @@
     3.4    edge_set_test
     3.5    error_test
     3.6    euler_test
     3.7 +  gomory_hu_test
     3.8    graph_copy_test
     3.9    graph_test
    3.10    graph_utils_test
     4.1 --- a/test/Makefile.am	Tue Feb 24 06:52:17 2009 +0000
     4.2 +++ b/test/Makefile.am	Fri Feb 20 17:17:17 2009 +0100
     4.3 @@ -17,6 +17,7 @@
     4.4  	test/edge_set_test \
     4.5  	test/error_test \
     4.6  	test/euler_test \
     4.7 +	test/gomory_hu_test \
     4.8  	test/graph_copy_test \
     4.9  	test/graph_test \
    4.10  	test/graph_utils_test \
    4.11 @@ -57,6 +58,7 @@
    4.12  test_edge_set_test_SOURCES = test/edge_set_test.cc
    4.13  test_error_test_SOURCES = test/error_test.cc
    4.14  test_euler_test_SOURCES = test/euler_test.cc
    4.15 +test_gomory_hu_test_SOURCES = test/gomory_hu_test.cc
    4.16  test_graph_copy_test_SOURCES = test/graph_copy_test.cc
    4.17  test_graph_test_SOURCES = test/graph_test.cc
    4.18  test_graph_utils_test_SOURCES = test/graph_utils_test.cc
     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/test/gomory_hu_test.cc	Fri Feb 20 17:17:17 2009 +0100
     5.3 @@ -0,0 +1,85 @@
     5.4 +#include <iostream>
     5.5 +
     5.6 +#include "test_tools.h"
     5.7 +#include <lemon/smart_graph.h>
     5.8 +#include <lemon/adaptors.h>
     5.9 +#include <lemon/lgf_reader.h>
    5.10 +#include <lemon/lgf_writer.h>
    5.11 +#include <lemon/dimacs.h>
    5.12 +#include <lemon/time_measure.h>
    5.13 +#include <lemon/gomory_hu_tree.h>
    5.14 +#include <cstdlib>
    5.15 +
    5.16 +using namespace std;
    5.17 +using namespace lemon;
    5.18 +
    5.19 +typedef SmartGraph Graph;
    5.20 +
    5.21 +char test_lgf[] =
    5.22 +  "@nodes\n"
    5.23 +  "label\n"
    5.24 +  "0\n"
    5.25 +  "1\n"
    5.26 +  "2\n"
    5.27 +  "3\n"
    5.28 +  "4\n"
    5.29 +  "@arcs\n"
    5.30 +  "     label capacity\n"
    5.31 +  "0 1  0     1\n"
    5.32 +  "1 2  1     1\n"
    5.33 +  "2 3  2     1\n"
    5.34 +  "0 3  4     5\n"
    5.35 +  "0 3  5     10\n"
    5.36 +  "0 3  6     7\n"
    5.37 +  "4 2  7     1\n"
    5.38 +  "@attributes\n"
    5.39 +  "source 0\n"
    5.40 +  "target 3\n";
    5.41 +  
    5.42 +GRAPH_TYPEDEFS(Graph);
    5.43 +typedef Graph::EdgeMap<int> IntEdgeMap;
    5.44 +typedef Graph::NodeMap<bool> BoolNodeMap;
    5.45 +
    5.46 +int cutValue(const Graph& graph, const BoolNodeMap& cut,
    5.47 +	     const IntEdgeMap& capacity) {
    5.48 +
    5.49 +  int sum = 0;
    5.50 +  for (EdgeIt e(graph); e != INVALID; ++e) {
    5.51 +    Node s = graph.u(e);
    5.52 +    Node t = graph.v(e);
    5.53 +
    5.54 +    if (cut[s] != cut[t]) {
    5.55 +      sum += capacity[e];
    5.56 +    }
    5.57 +  }
    5.58 +  return sum;
    5.59 +}
    5.60 +
    5.61 +
    5.62 +int main() {
    5.63 +  Graph graph;
    5.64 +  IntEdgeMap capacity(graph);
    5.65 +
    5.66 +  std::istringstream input(test_lgf);
    5.67 +  GraphReader<Graph>(graph, input).
    5.68 +    edgeMap("capacity", capacity).run();
    5.69 +
    5.70 +  GomoryHuTree<Graph> ght(graph, capacity);
    5.71 +  ght.init();
    5.72 +  ght.run();
    5.73 +
    5.74 +  for (NodeIt u(graph); u != INVALID; ++u) {
    5.75 +    for (NodeIt v(graph); v != u; ++v) {
    5.76 +      Preflow<Graph, IntEdgeMap> pf(graph, capacity, u, v);
    5.77 +      pf.runMinCut();
    5.78 +      BoolNodeMap cm(graph);
    5.79 +      ght.minCutMap(u, v, cm);
    5.80 +      check(pf.flowValue() == ght.minCutValue(u, v), "Wrong cut 1");
    5.81 +      check(cm[u] != cm[v], "Wrong cut 3");
    5.82 +      check(pf.flowValue() == cutValue(graph, cm, capacity), "Wrong cut 2");
    5.83 +      
    5.84 +    }
    5.85 +  }
    5.86 +  
    5.87 +  return 0;
    5.88 +}