Last change on this file since 596:293551ad254f was 596:293551ad254f, checked in by Peter Kovacs <kpeter@…>, 15 years ago

Improvements and fixes for the minimum cut algorithms (#264)

File size: 16.5 KB
Line
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_GOMORY_HU_TREE_H
20#define LEMON_GOMORY_HU_TREE_H
21
22#include <limits>
23
24#include <lemon/core.h>
25#include <lemon/preflow.h>
26#include <lemon/concept_check.h>
27#include <lemon/concepts/maps.h>
28
29/// \ingroup min_cut
30/// \file
31/// \brief Gomory-Hu cut tree in graphs.
32
33namespace lemon {
34
35  /// \ingroup min_cut
36  ///
37  /// \brief Gomory-Hu cut tree algorithm
38  ///
39  /// The Gomory-Hu tree is a tree on the node set of a given graph, but it
40  /// may contain edges which are not in the original graph. It has the
41  /// property that the minimum capacity edge of the path between two nodes
42  /// in this tree has the same weight as the minimum cut in the graph
43  /// between these nodes. Moreover the components obtained by removing
44  /// this edge from the tree determine the corresponding minimum cut.
45  /// Therefore once this tree is computed, the minimum cut between any pair
46  /// of nodes can easily be obtained.
47  ///
48  /// The algorithm calculates \e n-1 distinct minimum cuts (currently with
49  /// the \ref Preflow algorithm), thus it has \f$O(n^3\sqrt{e})\f$ overall
50  /// time complexity. It calculates a rooted Gomory-Hu tree.
51  /// The structure of the tree and the edge weights can be
52  /// obtained using \c predNode(), \c predValue() and \c rootDist().
53  /// The functions \c minCutMap() and \c minCutValue() calculate
54  /// the minimum cut and the minimum cut value between any two nodes
55  /// in the graph. You can also list (iterate on) the nodes and the
56  /// edges of the cuts using \c MinCutNodeIt and \c MinCutEdgeIt.
57  ///
58  /// \tparam GR The type of the undirected graph the algorithm runs on.
59  /// \tparam CAP The type of the edge map containing the capacities.
60  /// The default map type is \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
61#ifdef DOXYGEN
62  template <typename GR,
63            typename CAP>
64#else
65  template <typename GR,
66            typename CAP = typename GR::template EdgeMap<int> >
67#endif
68  class GomoryHu {
69  public:
70
71    /// The graph type of the algorithm
72    typedef GR Graph;
73    /// The capacity map type of the algorithm
74    typedef CAP Capacity;
75    /// The value type of capacities
76    typedef typename Capacity::Value Value;
77
78  private:
79
80    TEMPLATE_GRAPH_TYPEDEFS(Graph);
81
82    const Graph& _graph;
83    const Capacity& _capacity;
84
85    Node _root;
86    typename Graph::template NodeMap<Node>* _pred;
87    typename Graph::template NodeMap<Value>* _weight;
88    typename Graph::template NodeMap<int>* _order;
89
90    void createStructures() {
91      if (!_pred) {
92        _pred = new typename Graph::template NodeMap<Node>(_graph);
93      }
94      if (!_weight) {
95        _weight = new typename Graph::template NodeMap<Value>(_graph);
96      }
97      if (!_order) {
98        _order = new typename Graph::template NodeMap<int>(_graph);
99      }
100    }
101
102    void destroyStructures() {
103      if (_pred) {
104        delete _pred;
105      }
106      if (_weight) {
107        delete _weight;
108      }
109      if (_order) {
110        delete _order;
111      }
112    }
113
114  public:
115
116    /// \brief Constructor
117    ///
118    /// Constructor.
119    /// \param graph The undirected graph the algorithm runs on.
120    /// \param capacity The edge capacity map.
121    GomoryHu(const Graph& graph, const Capacity& capacity)
122      : _graph(graph), _capacity(capacity),
123        _pred(0), _weight(0), _order(0)
124    {
126    }
127
128
129    /// \brief Destructor
130    ///
131    /// Destructor.
132    ~GomoryHu() {
133      destroyStructures();
134    }
135
136  private:
137
138    // Initialize the internal data structures
139    void init() {
140      createStructures();
141
142      _root = NodeIt(_graph);
143      for (NodeIt n(_graph); n != INVALID; ++n) {
144        (*_pred)[n] = _root;
145        (*_order)[n] = -1;
146      }
147      (*_pred)[_root] = INVALID;
148      (*_weight)[_root] = std::numeric_limits<Value>::max();
149    }
150
151
152    // Start the algorithm
153    void start() {
154      Preflow<Graph, Capacity> fa(_graph, _capacity, _root, INVALID);
155
156      for (NodeIt n(_graph); n != INVALID; ++n) {
157        if (n == _root) continue;
158
159        Node pn = (*_pred)[n];
160        fa.source(n);
161        fa.target(pn);
162
163        fa.runMinCut();
164
165        (*_weight)[n] = fa.flowValue();
166
167        for (NodeIt nn(_graph); nn != INVALID; ++nn) {
168          if (nn != n && fa.minCut(nn) && (*_pred)[nn] == pn) {
169            (*_pred)[nn] = n;
170          }
171        }
172        if ((*_pred)[pn] != INVALID && fa.minCut((*_pred)[pn])) {
173          (*_pred)[n] = (*_pred)[pn];
174          (*_pred)[pn] = n;
175          (*_weight)[n] = (*_weight)[pn];
176          (*_weight)[pn] = fa.flowValue();
177        }
178      }
179
180      (*_order)[_root] = 0;
181      int index = 1;
182
183      for (NodeIt n(_graph); n != INVALID; ++n) {
184        std::vector<Node> st;
185        Node nn = n;
186        while ((*_order)[nn] == -1) {
187          st.push_back(nn);
188          nn = (*_pred)[nn];
189        }
190        while (!st.empty()) {
191          (*_order)[st.back()] = index++;
192          st.pop_back();
193        }
194      }
195    }
196
197  public:
198
199    ///\name Execution Control
200
201    ///@{
202
203    /// \brief Run the Gomory-Hu algorithm.
204    ///
205    /// This function runs the Gomory-Hu algorithm.
206    void run() {
207      init();
208      start();
209    }
210
211    /// @}
212
213    ///\name Query Functions
214    ///The results of the algorithm can be obtained using these
215    ///functions.\n
216    ///\ref run() should be called before using them.\n
218
219    ///@{
220
221    /// \brief Return the predecessor node in the Gomory-Hu tree.
222    ///
223    /// This function returns the predecessor node of the given node
224    /// in the Gomory-Hu tree.
225    /// If \c node is the root of the tree, then it returns \c INVALID.
226    ///
227    /// \pre \ref run() must be called before using this function.
228    Node predNode(const Node& node) const {
229      return (*_pred)[node];
230    }
231
232    /// \brief Return the weight of the predecessor edge in the
233    /// Gomory-Hu tree.
234    ///
235    /// This function returns the weight of the predecessor edge of the
236    /// given node in the Gomory-Hu tree.
237    /// If \c node is the root of the tree, the result is undefined.
238    ///
239    /// \pre \ref run() must be called before using this function.
240    Value predValue(const Node& node) const {
241      return (*_weight)[node];
242    }
243
244    /// \brief Return the distance from the root node in the Gomory-Hu tree.
245    ///
246    /// This function returns the distance of the given node from the root
247    /// node in the Gomory-Hu tree.
248    ///
249    /// \pre \ref run() must be called before using this function.
250    int rootDist(const Node& node) const {
251      return (*_order)[node];
252    }
253
254    /// \brief Return the minimum cut value between two nodes
255    ///
256    /// This function returns the minimum cut value between the nodes
257    /// \c s and \c t.
258    /// It finds the nearest common ancestor of the given nodes in the
259    /// Gomory-Hu tree and calculates the minimum weight edge on the
260    /// paths to the ancestor.
261    ///
262    /// \pre \ref run() must be called before using this function.
263    Value minCutValue(const Node& s, const Node& t) const {
264      Node sn = s, tn = t;
265      Value value = std::numeric_limits<Value>::max();
266
267      while (sn != tn) {
268        if ((*_order)[sn] < (*_order)[tn]) {
269          if ((*_weight)[tn] <= value) value = (*_weight)[tn];
270          tn = (*_pred)[tn];
271        } else {
272          if ((*_weight)[sn] <= value) value = (*_weight)[sn];
273          sn = (*_pred)[sn];
274        }
275      }
276      return value;
277    }
278
279    /// \brief Return the minimum cut between two nodes
280    ///
281    /// This function returns the minimum cut between the nodes \c s and \c t
282    /// in the \c cutMap parameter by setting the nodes in the component of
283    /// \c s to \c true and the other nodes to \c false.
284    ///
285    /// For higher level interfaces see MinCutNodeIt and MinCutEdgeIt.
286    ///
287    /// \param s The base node.
288    /// \param t The node you want to separate from node \c s.
289    /// \param cutMap The cut will be returned in this map.
290    /// It must be a \c bool (or convertible) \ref concepts::ReadWriteMap
291    /// "ReadWriteMap" on the graph nodes.
292    ///
293    /// \return The value of the minimum cut between \c s and \c t.
294    ///
295    /// \pre \ref run() must be called before using this function.
296    template <typename CutMap>
297    Value minCutMap(const Node& s, ///<
298                    const Node& t,
299                    ///<
300                    CutMap& cutMap
301                    ///<
302                    ) const {
303      Node sn = s, tn = t;
304      bool s_root=false;
305      Node rn = INVALID;
306      Value value = std::numeric_limits<Value>::max();
307
308      while (sn != tn) {
309        if ((*_order)[sn] < (*_order)[tn]) {
310          if ((*_weight)[tn] <= value) {
311            rn = tn;
312            s_root = false;
313            value = (*_weight)[tn];
314          }
315          tn = (*_pred)[tn];
316        } else {
317          if ((*_weight)[sn] <= value) {
318            rn = sn;
319            s_root = true;
320            value = (*_weight)[sn];
321          }
322          sn = (*_pred)[sn];
323        }
324      }
325
326      typename Graph::template NodeMap<bool> reached(_graph, false);
327      reached[_root] = true;
328      cutMap.set(_root, !s_root);
329      reached[rn] = true;
330      cutMap.set(rn, s_root);
331
332      std::vector<Node> st;
333      for (NodeIt n(_graph); n != INVALID; ++n) {
334        st.clear();
335        Node nn = n;
336        while (!reached[nn]) {
337          st.push_back(nn);
338          nn = (*_pred)[nn];
339        }
340        while (!st.empty()) {
341          cutMap.set(st.back(), cutMap[nn]);
342          st.pop_back();
343        }
344      }
345
346      return value;
347    }
348
349    ///@}
350
351    friend class MinCutNodeIt;
352
353    /// Iterate on the nodes of a minimum cut
354
355    /// This iterator class lists the nodes of a minimum cut found by
356    /// GomoryHu. Before using it, you must allocate a GomoryHu class
357    /// and call its \ref GomoryHu::run() "run()" method.
358    ///
359    /// This example counts the nodes in the minimum cut separating \c s from
360    /// \c t.
361    /// \code
362    /// GomoruHu<Graph> gom(g, capacities);
363    /// gom.run();
364    /// int cnt=0;
365    /// for(GomoruHu<Graph>::MinCutNodeIt n(gom,s,t); n!=INVALID; ++n) ++cnt;
366    /// \endcode
367    class MinCutNodeIt
368    {
369      bool _side;
370      typename Graph::NodeIt _node_it;
371      typename Graph::template NodeMap<bool> _cut;
372    public:
373      /// Constructor
374
375      /// Constructor.
376      ///
377      MinCutNodeIt(GomoryHu const &gomory,
378                   ///< The GomoryHu class. You must call its
379                   ///  run() method
380                   ///  before initializing this iterator.
381                   const Node& s, ///< The base node.
382                   const Node& t,
383                   ///< The node you want to separate from node \c s.
384                   bool side=true
385                   ///< If it is \c true (default) then the iterator lists
386                   ///  the nodes of the component containing \c s,
387                   ///  otherwise it lists the other component.
388                   /// \note As the minimum cut is not always unique,
389                   /// \code
390                   /// MinCutNodeIt(gomory, s, t, true);
391                   /// \endcode
392                   /// and
393                   /// \code
394                   /// MinCutNodeIt(gomory, t, s, false);
395                   /// \endcode
396                   /// does not necessarily give the same set of nodes.
397                   /// However it is ensured that
398                   /// \code
399                   /// MinCutNodeIt(gomory, s, t, true);
400                   /// \endcode
401                   /// and
402                   /// \code
403                   /// MinCutNodeIt(gomory, s, t, false);
404                   /// \endcode
405                   /// together list each node exactly once.
406                   )
407        : _side(side), _cut(gomory._graph)
408      {
409        gomory.minCutMap(s,t,_cut);
410        for(_node_it=typename Graph::NodeIt(gomory._graph);
411            _node_it!=INVALID && _cut[_node_it]!=_side;
412            ++_node_it) {}
413      }
414      /// Conversion to \c Node
415
416      /// Conversion to \c Node.
417      ///
418      operator typename Graph::Node() const
419      {
420        return _node_it;
421      }
422      bool operator==(Invalid) { return _node_it==INVALID; }
423      bool operator!=(Invalid) { return _node_it!=INVALID; }
424      /// Next node
425
426      /// Next node.
427      ///
428      MinCutNodeIt &operator++()
429      {
430        for(++_node_it;_node_it!=INVALID&&_cut[_node_it]!=_side;++_node_it) {}
431        return *this;
432      }
433      /// Postfix incrementation
434
435      /// Postfix incrementation.
436      ///
437      /// \warning This incrementation
438      /// returns a \c Node, not a \c MinCutNodeIt, as one may
439      /// expect.
440      typename Graph::Node operator++(int)
441      {
442        typename Graph::Node n=*this;
443        ++(*this);
444        return n;
445      }
446    };
447
448    friend class MinCutEdgeIt;
449
450    /// Iterate on the edges of a minimum cut
451
452    /// This iterator class lists the edges of a minimum cut found by
453    /// GomoryHu. Before using it, you must allocate a GomoryHu class
454    /// and call its \ref GomoryHu::run() "run()" method.
455    ///
456    /// This example computes the value of the minimum cut separating \c s from
457    /// \c t.
458    /// \code
459    /// GomoruHu<Graph> gom(g, capacities);
460    /// gom.run();
461    /// int value=0;
462    /// for(GomoruHu<Graph>::MinCutEdgeIt e(gom,s,t); e!=INVALID; ++e)
463    ///   value+=capacities[e];
464    /// \endcode
465    /// The result will be the same as the value returned by
466    /// \ref GomoryHu::minCutValue() "gom.minCutValue(s,t)".
467    class MinCutEdgeIt
468    {
469      bool _side;
470      const Graph &_graph;
471      typename Graph::NodeIt _node_it;
472      typename Graph::OutArcIt _arc_it;
473      typename Graph::template NodeMap<bool> _cut;
474      void step()
475      {
476        ++_arc_it;
477        while(_node_it!=INVALID && _arc_it==INVALID)
478          {
479            for(++_node_it;_node_it!=INVALID&&!_cut[_node_it];++_node_it) {}
480            if(_node_it!=INVALID)
481              _arc_it=typename Graph::OutArcIt(_graph,_node_it);
482          }
483      }
484
485    public:
486      /// Constructor
487
488      /// Constructor.
489      ///
490      MinCutEdgeIt(GomoryHu const &gomory,
491                   ///< The GomoryHu class. You must call its
492                   ///  run() method
493                   ///  before initializing this iterator.
494                   const Node& s,  ///< The base node.
495                   const Node& t,
496                   ///< The node you want to separate from node \c s.
497                   bool side=true
498                   ///< If it is \c true (default) then the listed arcs
499                   ///  will be oriented from the
500                   ///  nodes of the component containing \c s,
501                   ///  otherwise they will be oriented in the opposite
502                   ///  direction.
503                   )
504        : _graph(gomory._graph), _cut(_graph)
505      {
506        gomory.minCutMap(s,t,_cut);
507        if(!side)
508          for(typename Graph::NodeIt n(_graph);n!=INVALID;++n)
509            _cut[n]=!_cut[n];
510
511        for(_node_it=typename Graph::NodeIt(_graph);
512            _node_it!=INVALID && !_cut[_node_it];
513            ++_node_it) {}
514        _arc_it = _node_it!=INVALID ?
515          typename Graph::OutArcIt(_graph,_node_it) : INVALID;
516        while(_node_it!=INVALID && _arc_it == INVALID)
517          {
518            for(++_node_it; _node_it!=INVALID&&!_cut[_node_it]; ++_node_it) {}
519            if(_node_it!=INVALID)
520              _arc_it= typename Graph::OutArcIt(_graph,_node_it);
521          }
522        while(_arc_it!=INVALID && _cut[_graph.target(_arc_it)]) step();
523      }
524      /// Conversion to \c Arc
525
526      /// Conversion to \c Arc.
527      ///
528      operator typename Graph::Arc() const
529      {
530        return _arc_it;
531      }
532      /// Conversion to \c Edge
533
534      /// Conversion to \c Edge.
535      ///
536      operator typename Graph::Edge() const
537      {
538        return _arc_it;
539      }
540      bool operator==(Invalid) { return _node_it==INVALID; }
541      bool operator!=(Invalid) { return _node_it!=INVALID; }
542      /// Next edge
543
544      /// Next edge.
545      ///
546      MinCutEdgeIt &operator++()
547      {
548        step();
549        while(_arc_it!=INVALID && _cut[_graph.target(_arc_it)]) step();
550        return *this;
551      }
552      /// Postfix incrementation
553
554      /// Postfix incrementation.
555      ///
556      /// \warning This incrementation
557      /// returns an \c Arc, not a \c MinCutEdgeIt, as one may expect.
558      typename Graph::Arc operator++(int)
559      {
560        typename Graph::Arc e=*this;
561        ++(*this);
562        return e;
563      }
564    };
565
566  };
567
568}
569
570#endif
Note: See TracBrowser for help on using the repository browser.