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