lemon/gomory_hu.h
author Peter Kovacs <kpeter@inf.elte.hu>
Thu, 06 Aug 2009 20:12:43 +0200
changeset 760 83ce7ce39f21
parent 581 aa1804409f29
child 713 4ac30454f1c1
permissions -rw-r--r--
Rework and fix the implementation of MinMeanCycle (#179)

- Fix the handling of the cycle means.
- Many implementation improvements:
- More efficient data storage for the strongly connected
components.
- Better handling of BFS queues.
- Merge consecutive BFS searches (perform two BFS searches
instead of three).

This version is about two times faster on average and an order of
magnitude faster if there are a lot of strongly connected components.
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2008
     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                     ///< 
   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