lemon/gomory_hu.h
author Alpar Juttner <alpar@cs.elte.hu>
Tue, 06 Aug 2013 12:28:37 +0200
changeset 1079 5958cc5c0a98
parent 786 e20173729589
child 1080 c5cd8960df74
permissions -rw-r--r--
Merge further fixes #470
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2010
     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 
    33 namespace 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     {
   125       checkConcept<concepts::ReadMap<Edge, Value>, Capacity>();
   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
   217     ///See also \ref MinCutNodeIt and \ref MinCutEdgeIt.
   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                     CutMap& cutMap
   300                     ) const {
   301       Node sn = s, tn = t;
   302       bool s_root=false;
   303       Node rn = INVALID;
   304       Value value = std::numeric_limits<Value>::max();
   305 
   306       while (sn != tn) {
   307         if ((*_order)[sn] < (*_order)[tn]) {
   308           if ((*_weight)[tn] <= value) {
   309             rn = tn;
   310             s_root = false;
   311             value = (*_weight)[tn];
   312           }
   313           tn = (*_pred)[tn];
   314         } else {
   315           if ((*_weight)[sn] <= value) {
   316             rn = sn;
   317             s_root = true;
   318             value = (*_weight)[sn];
   319           }
   320           sn = (*_pred)[sn];
   321         }
   322       }
   323 
   324       typename Graph::template NodeMap<bool> reached(_graph, false);
   325       reached[_root] = true;
   326       cutMap.set(_root, !s_root);
   327       reached[rn] = true;
   328       cutMap.set(rn, s_root);
   329 
   330       std::vector<Node> st;
   331       for (NodeIt n(_graph); n != INVALID; ++n) {
   332         st.clear();
   333         Node nn = n;
   334         while (!reached[nn]) {
   335           st.push_back(nn);
   336           nn = (*_pred)[nn];
   337         }
   338         while (!st.empty()) {
   339           cutMap.set(st.back(), cutMap[nn]);
   340           st.pop_back();
   341         }
   342       }
   343 
   344       return value;
   345     }
   346 
   347     ///@}
   348 
   349     friend class MinCutNodeIt;
   350 
   351     /// Iterate on the nodes of a minimum cut
   352 
   353     /// This iterator class lists the nodes of a minimum cut found by
   354     /// GomoryHu. Before using it, you must allocate a GomoryHu class
   355     /// and call its \ref GomoryHu::run() "run()" method.
   356     ///
   357     /// This example counts the nodes in the minimum cut separating \c s from
   358     /// \c t.
   359     /// \code
   360     /// GomoryHu<Graph> gom(g, capacities);
   361     /// gom.run();
   362     /// int cnt=0;
   363     /// for(GomoryHu<Graph>::MinCutNodeIt n(gom,s,t); n!=INVALID; ++n) ++cnt;
   364     /// \endcode
   365     class MinCutNodeIt
   366     {
   367       bool _side;
   368       typename Graph::NodeIt _node_it;
   369       typename Graph::template NodeMap<bool> _cut;
   370     public:
   371       /// Constructor
   372 
   373       /// Constructor.
   374       ///
   375       MinCutNodeIt(GomoryHu const &gomory,
   376                    ///< The GomoryHu class. You must call its
   377                    ///  run() method
   378                    ///  before initializing this iterator.
   379                    const Node& s, ///< The base node.
   380                    const Node& t,
   381                    ///< The node you want to separate from node \c s.
   382                    bool side=true
   383                    ///< If it is \c true (default) then the iterator lists
   384                    ///  the nodes of the component containing \c s,
   385                    ///  otherwise it lists the other component.
   386                    /// \note As the minimum cut is not always unique,
   387                    /// \code
   388                    /// MinCutNodeIt(gomory, s, t, true);
   389                    /// \endcode
   390                    /// and
   391                    /// \code
   392                    /// MinCutNodeIt(gomory, t, s, false);
   393                    /// \endcode
   394                    /// does not necessarily give the same set of nodes.
   395                    /// However, it is ensured that
   396                    /// \code
   397                    /// MinCutNodeIt(gomory, s, t, true);
   398                    /// \endcode
   399                    /// and
   400                    /// \code
   401                    /// MinCutNodeIt(gomory, s, t, false);
   402                    /// \endcode
   403                    /// together list each node exactly once.
   404                    )
   405         : _side(side), _cut(gomory._graph)
   406       {
   407         gomory.minCutMap(s,t,_cut);
   408         for(_node_it=typename Graph::NodeIt(gomory._graph);
   409             _node_it!=INVALID && _cut[_node_it]!=_side;
   410             ++_node_it) {}
   411       }
   412       /// Conversion to \c Node
   413 
   414       /// Conversion to \c Node.
   415       ///
   416       operator typename Graph::Node() const
   417       {
   418         return _node_it;
   419       }
   420       bool operator==(Invalid) { return _node_it==INVALID; }
   421       bool operator!=(Invalid) { return _node_it!=INVALID; }
   422       /// Next node
   423 
   424       /// Next node.
   425       ///
   426       MinCutNodeIt &operator++()
   427       {
   428         for(++_node_it;_node_it!=INVALID&&_cut[_node_it]!=_side;++_node_it) {}
   429         return *this;
   430       }
   431       /// Postfix incrementation
   432 
   433       /// Postfix incrementation.
   434       ///
   435       /// \warning This incrementation
   436       /// returns a \c Node, not a \c MinCutNodeIt, as one may
   437       /// expect.
   438       typename Graph::Node operator++(int)
   439       {
   440         typename Graph::Node n=*this;
   441         ++(*this);
   442         return n;
   443       }
   444     };
   445 
   446     friend class MinCutEdgeIt;
   447 
   448     /// Iterate on the edges of a minimum cut
   449 
   450     /// This iterator class lists the edges of a minimum cut found by
   451     /// GomoryHu. Before using it, you must allocate a GomoryHu class
   452     /// and call its \ref GomoryHu::run() "run()" method.
   453     ///
   454     /// This example computes the value of the minimum cut separating \c s from
   455     /// \c t.
   456     /// \code
   457     /// GomoryHu<Graph> gom(g, capacities);
   458     /// gom.run();
   459     /// int value=0;
   460     /// for(GomoryHu<Graph>::MinCutEdgeIt e(gom,s,t); e!=INVALID; ++e)
   461     ///   value+=capacities[e];
   462     /// \endcode
   463     /// The result will be the same as the value returned by
   464     /// \ref GomoryHu::minCutValue() "gom.minCutValue(s,t)".
   465     class MinCutEdgeIt
   466     {
   467       bool _side;
   468       const Graph &_graph;
   469       typename Graph::NodeIt _node_it;
   470       typename Graph::OutArcIt _arc_it;
   471       typename Graph::template NodeMap<bool> _cut;
   472       void step()
   473       {
   474         ++_arc_it;
   475         while(_node_it!=INVALID && _arc_it==INVALID)
   476           {
   477             for(++_node_it;_node_it!=INVALID&&!_cut[_node_it];++_node_it) {}
   478             if(_node_it!=INVALID)
   479               _arc_it=typename Graph::OutArcIt(_graph,_node_it);
   480           }
   481       }
   482 
   483     public:
   484       /// Constructor
   485 
   486       /// Constructor.
   487       ///
   488       MinCutEdgeIt(GomoryHu const &gomory,
   489                    ///< The GomoryHu class. You must call its
   490                    ///  run() method
   491                    ///  before initializing this iterator.
   492                    const Node& s,  ///< The base node.
   493                    const Node& t,
   494                    ///< The node you want to separate from node \c s.
   495                    bool side=true
   496                    ///< If it is \c true (default) then the listed arcs
   497                    ///  will be oriented from the
   498                    ///  nodes of the component containing \c s,
   499                    ///  otherwise they will be oriented in the opposite
   500                    ///  direction.
   501                    )
   502         : _graph(gomory._graph), _cut(_graph)
   503       {
   504         gomory.minCutMap(s,t,_cut);
   505         if(!side)
   506           for(typename Graph::NodeIt n(_graph);n!=INVALID;++n)
   507             _cut[n]=!_cut[n];
   508 
   509         for(_node_it=typename Graph::NodeIt(_graph);
   510             _node_it!=INVALID && !_cut[_node_it];
   511             ++_node_it) {}
   512         _arc_it = _node_it!=INVALID ?
   513           typename Graph::OutArcIt(_graph,_node_it) : INVALID;
   514         while(_node_it!=INVALID && _arc_it == INVALID)
   515           {
   516             for(++_node_it; _node_it!=INVALID&&!_cut[_node_it]; ++_node_it) {}
   517             if(_node_it!=INVALID)
   518               _arc_it= typename Graph::OutArcIt(_graph,_node_it);
   519           }
   520         while(_arc_it!=INVALID && _cut[_graph.target(_arc_it)]) step();
   521       }
   522       /// Conversion to \c Arc
   523 
   524       /// Conversion to \c Arc.
   525       ///
   526       operator typename Graph::Arc() const
   527       {
   528         return _arc_it;
   529       }
   530       /// Conversion to \c Edge
   531 
   532       /// Conversion to \c Edge.
   533       ///
   534       operator typename Graph::Edge() const
   535       {
   536         return _arc_it;
   537       }
   538       bool operator==(Invalid) { return _node_it==INVALID; }
   539       bool operator!=(Invalid) { return _node_it!=INVALID; }
   540       /// Next edge
   541 
   542       /// Next edge.
   543       ///
   544       MinCutEdgeIt &operator++()
   545       {
   546         step();
   547         while(_arc_it!=INVALID && _cut[_graph.target(_arc_it)]) step();
   548         return *this;
   549       }
   550       /// Postfix incrementation
   551 
   552       /// Postfix incrementation.
   553       ///
   554       /// \warning This incrementation
   555       /// returns an \c Arc, not a \c MinCutEdgeIt, as one may expect.
   556       typename Graph::Arc operator++(int)
   557       {
   558         typename Graph::Arc e=*this;
   559         ++(*this);
   560         return e;
   561       }
   562     };
   563 
   564   };
   565 
   566 }
   567 
   568 #endif