lemon/pr_bipartite_matching.h
author ladanyi
Sat, 13 Oct 2007 08:48:07 +0000
changeset 2495 e4f8367beb41
parent 2463 19651a04d056
child 2512 371cf309fc3c
permissions -rw-r--r--
Added the function isFinite(), and replaced the calls to finite() with it.
This was necessary because finite() is not a standard function. Neither can
we use its standard counterpart isfinite(), because it was introduced only
in C99, and therefore it is not supplied by all C++ implementations.
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2007
     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_PR_BIPARTITE_MATCHING
    20 #define LEMON_PR_BIPARTITE_MATCHING
    21 
    22 #include <lemon/graph_utils.h>
    23 #include <lemon/iterable_maps.h>
    24 #include <iostream>
    25 #include <queue>
    26 #include <lemon/elevator.h>
    27 
    28 ///\ingroup matching
    29 ///\file
    30 ///\brief Push-prelabel maximum matching algorithms in bipartite graphs.
    31 ///
    32 namespace lemon {
    33 
    34   ///Max cardinality matching algorithm based on push-relabel principle
    35 
    36   ///\ingroup matching
    37   ///Bipartite Max Cardinality Matching algorithm. This class uses the
    38   ///push-relabel principle which in several cases has better runtime
    39   ///performance than the augmenting path solutions.
    40   ///
    41   ///\author Alpar Juttner
    42   template<class Graph>
    43   class PrBipartiteMatching {
    44     typedef typename Graph::Node Node;
    45     typedef typename Graph::ANodeIt ANodeIt;
    46     typedef typename Graph::BNodeIt BNodeIt;
    47     typedef typename Graph::UEdge UEdge;
    48     typedef typename Graph::UEdgeIt UEdgeIt;
    49     typedef typename Graph::IncEdgeIt IncEdgeIt;
    50     
    51     const Graph &_g;
    52     int _node_num;
    53     int _matching_size;
    54     int _empty_level;
    55     
    56     typename Graph::template ANodeMap<typename Graph::UEdge> _matching;
    57     Elevator<Graph,typename Graph::BNode> _levels;
    58     typename Graph::template BNodeMap<int> _cov;
    59 
    60   public:
    61 
    62     /// Constructor
    63 
    64     /// Constructor
    65     ///
    66     PrBipartiteMatching(const Graph &g) :
    67       _g(g),
    68       _node_num(countBNodes(g)),
    69       _matching(g),
    70       _levels(g,_node_num),
    71       _cov(g,0)
    72     {
    73     }
    74     
    75     /// \name Execution control 
    76     /// The simplest way to execute the algorithm is to use one of the
    77     /// member functions called \c run(). \n
    78     /// If you need more control on the execution, first
    79     /// you must call \ref init() and then one variant of the start()
    80     /// member. 
    81 
    82     /// @{
    83 
    84     ///Initialize the data structures
    85 
    86     ///This function constructs a prematching first, which is a
    87     ///regular matching on the A-side of the graph, but on the B-side
    88     ///each node could cover more matching edges. After that, the
    89     ///B-nodes which multiple matched, will be pushed into the lowest
    90     ///level of the Elevator. The remaning B-nodes will be pushed to
    91     ///the consequent levels respect to a Bfs on following graph: the
    92     ///nodes are the B-nodes of the original bipartite graph and two
    93     ///nodes are adjacent if a node can pass over a matching edge to
    94     ///an other node. The source of the Bfs are the lowest level
    95     ///nodes. Last, the reached B-nodes without covered matching edge
    96     ///becomes active.
    97     void init() {
    98       _matching_size=0;
    99       _empty_level=_node_num;
   100       for(ANodeIt n(_g);n!=INVALID;++n)
   101 	if((_matching[n]=IncEdgeIt(_g,n))!=INVALID)
   102 	  ++_cov[_g.bNode(_matching[n])];
   103 
   104       std::queue<Node> q;
   105       _levels.initStart();
   106       for(BNodeIt n(_g);n!=INVALID;++n)
   107 	if(_cov[n]>1) {
   108 	  _levels.initAddItem(n);
   109 	  q.push(n);
   110 	}
   111       int hlev=0;
   112       while(!q.empty()) {
   113 	Node n=q.front();
   114 	q.pop();
   115 	int nlev=_levels[n]+1;
   116 	for(IncEdgeIt e(_g,n);e!=INVALID;++e) {
   117 	  Node m=_g.runningNode(e);
   118 	  if(e==_matching[m]) {
   119 	    for(IncEdgeIt f(_g,m);f!=INVALID;++f) {
   120 	      Node r=_g.runningNode(f);
   121 	      if(_levels[r]>nlev) {
   122 		for(;nlev>hlev;hlev++)
   123 		  _levels.initNewLevel();
   124 		_levels.initAddItem(r);
   125 		q.push(r);
   126 	      }
   127 	    }
   128 	  }
   129 	}
   130       }
   131       _levels.initFinish();
   132       for(BNodeIt n(_g);n!=INVALID;++n)
   133 	if(_cov[n]<1&&_levels[n]<_node_num)
   134 	  _levels.activate(n);
   135     }
   136 
   137     ///Start the main phase of the algorithm
   138     
   139     ///This algorithm calculates the maximum matching with the
   140     ///push-relabel principle. This function should be called just
   141     ///after the init() function which already set the initial
   142     ///prematching, the level function on the B-nodes and the active,
   143     ///ie. unmatched B-nodes.
   144     ///
   145     ///The algorithm always takes highest active B-node, and it try to
   146     ///find a B-node which is eligible to pass over one of it's
   147     ///matching edge. This condition holds when the B-node is one
   148     ///level lower, and the opposite node of it's matching edge is
   149     ///adjacent to the highest active node. In this case the current
   150     ///node steals the matching edge and becomes inactive. If there is
   151     ///not eligible node then the highest active node should be lift
   152     ///to the next proper level.
   153     ///
   154     ///The nodes should not lift higher than the number of the
   155     ///B-nodes, if a node reach this level it remains unmatched. If
   156     ///during the execution one level becomes empty the nodes above it
   157     ///can be deactivated and lift to the highest level.
   158     void start() {
   159       Node act;
   160       Node bact=INVALID;
   161       Node last_activated=INVALID;
   162       while((act=_levels.highestActive())!=INVALID) {
   163 	last_activated=INVALID;
   164 	int actlevel=_levels[act];
   165 	
   166 	UEdge bedge=INVALID;
   167 	int nlevel=_node_num;
   168 	{
   169 	  int nnlevel;
   170 	  for(IncEdgeIt tbedge(_g,act);
   171 	      tbedge!=INVALID && nlevel>=actlevel;
   172 	      ++tbedge)
   173 	    if((nnlevel=_levels[_g.bNode(_matching[_g.runningNode(tbedge)])])<
   174 	       nlevel)
   175 	      {
   176 		nlevel=nnlevel;
   177 		bedge=tbedge;
   178 	      }
   179 	}
   180 	if(nlevel<_node_num) {
   181 	  if(nlevel>=actlevel)
   182 	    _levels.liftHighestActiveTo(nlevel+1);
   183 	  bact=_g.bNode(_matching[_g.aNode(bedge)]);
   184 	  if(--_cov[bact]<1) {
   185 	    _levels.activate(bact);
   186 	    last_activated=bact;
   187 	  }
   188 	  _matching[_g.aNode(bedge)]=bedge;
   189 	  _cov[act]=1;
   190 	  _levels.deactivate(act);
   191 	}
   192 	else {
   193 	  if(_node_num>actlevel) 
   194 	    _levels.liftHighestActiveTo(_node_num);
   195 	  _levels.deactivate(act); 
   196 	}
   197 
   198 	if(_levels.onLevel(actlevel)==0)
   199 	  _levels.liftToTop(actlevel);
   200       }
   201       
   202       for(ANodeIt n(_g);n!=INVALID;++n) {
   203 	if (_matching[n]==INVALID)continue;
   204 	if (_cov[_g.bNode(_matching[n])]>1) {
   205 	  _cov[_g.bNode(_matching[n])]--;
   206 	  _matching[n]=INVALID;
   207 	} else {
   208 	  ++_matching_size;
   209 	}
   210       }
   211     }
   212 
   213     ///Start the algorithm to find a perfect matching
   214 
   215     ///This function is close to identical to the simple start()
   216     ///member function but it calculates just perfect matching.
   217     ///However, the perfect property is only checked on the B-side of
   218     ///the graph
   219     ///
   220     ///The main difference between the two function is the handling of
   221     ///the empty levels. The simple start() function let the nodes
   222     ///above the empty levels unmatched while this variant if it find
   223     ///an empty level immediately terminates and gives back false
   224     ///return value.
   225     bool startPerfect() {
   226       Node act;
   227       Node bact=INVALID;
   228       Node last_activated=INVALID;
   229       while((act=_levels.highestActive())!=INVALID) {
   230 	last_activated=INVALID;
   231 	int actlevel=_levels[act];
   232 	
   233 	UEdge bedge=INVALID;
   234 	int nlevel=_node_num;
   235 	{
   236 	  int nnlevel;
   237 	  for(IncEdgeIt tbedge(_g,act);
   238 	      tbedge!=INVALID && nlevel>=actlevel;
   239 	      ++tbedge)
   240 	    if((nnlevel=_levels[_g.bNode(_matching[_g.runningNode(tbedge)])])<
   241 	       nlevel)
   242 	      {
   243 		nlevel=nnlevel;
   244 		bedge=tbedge;
   245 	      }
   246 	}
   247 	if(nlevel<_node_num) {
   248 	  if(nlevel>=actlevel)
   249 	    _levels.liftHighestActiveTo(nlevel+1);
   250 	  bact=_g.bNode(_matching[_g.aNode(bedge)]);
   251 	  if(--_cov[bact]<1) {
   252 	    _levels.activate(bact);
   253 	    last_activated=bact;
   254 	  }
   255 	  _matching[_g.aNode(bedge)]=bedge;
   256 	  _cov[act]=1;
   257 	  _levels.deactivate(act);
   258 	}
   259 	else {
   260 	  if(_node_num>actlevel) 
   261 	    _levels.liftHighestActiveTo(_node_num);
   262 	  _levels.deactivate(act); 
   263 	}
   264 
   265 	if(_levels.onLevel(actlevel)==0)
   266 	  _empty_level=actlevel;
   267 	  return false;
   268       }
   269       _matching_size = _node_num;
   270       return true;
   271     }
   272   
   273     ///Runs the algorithm
   274     
   275     ///Just a shortcut for the next code:
   276     ///\code
   277     /// init();
   278     /// start();
   279     ///\endcode
   280     void run() {
   281       init();
   282       start();
   283     }
   284     
   285     ///Runs the algorithm to find a perfect matching
   286     
   287     ///Just a shortcut for the next code:
   288     ///\code
   289     /// init();
   290     /// startPerfect();
   291     ///\endcode
   292     ///
   293     ///\note If the two nodesets of the graph have different size then
   294     ///this algorithm checks the perfect property on the B-side.
   295     bool runPerfect() {
   296       init();
   297       return startPerfect();
   298     }
   299 
   300     ///Runs the algorithm to find a perfect matching
   301     
   302     ///Just a shortcut for the next code:
   303     ///\code
   304     /// init();
   305     /// startPerfect();
   306     ///\endcode
   307     ///
   308     ///\note It checks that the size of the two nodesets are equal.
   309     bool checkedRunPerfect() {
   310       if (countANodes(_g) != _node_num) return false;
   311       init();
   312       return startPerfect();
   313     }
   314 
   315     ///@}
   316 
   317     /// \name Query Functions
   318     /// The result of the %Matching algorithm can be obtained using these
   319     /// functions.\n
   320     /// Before the use of these functions,
   321     /// either run() or start() must be called.
   322     ///@{
   323 
   324     ///Set true all matching uedge in the map.
   325 
   326     ///Set true all matching uedge in the map. It does not change the
   327     ///value mapped to the other uedges.
   328     ///\return The number of the matching edges.
   329     template <typename MatchingMap>
   330     int quickMatching(MatchingMap& mm) const {
   331       for (ANodeIt n(_g);n!=INVALID;++n) {
   332         if (_matching[n]!=INVALID) mm.set(_matching[n],true);
   333       }
   334       return _matching_size;
   335     }
   336 
   337     ///Set true all matching uedge in the map and the others to false.
   338 
   339     ///Set true all matching uedge in the map and the others to false.
   340     ///\return The number of the matching edges.
   341     template<class MatchingMap>
   342     int matching(MatchingMap& mm) const {
   343       for (UEdgeIt e(_g);e!=INVALID;++e) {
   344         mm.set(e,e==_matching[_g.aNode(e)]);
   345       }
   346       return _matching_size;
   347     }
   348 
   349     ///Gives back the matching in an ANodeMap.
   350 
   351     ///Gives back the matching in an ANodeMap. The parameter should
   352     ///be a write ANodeMap of UEdge values.
   353     ///\return The number of the matching edges.
   354     template<class MatchingMap>
   355     int aMatching(MatchingMap& mm) const {
   356       for (ANodeIt n(_g);n!=INVALID;++n) {
   357         mm.set(n,_matching[n]);
   358       }
   359       return _matching_size;
   360     }
   361 
   362     ///Gives back the matching in a BNodeMap.
   363 
   364     ///Gives back the matching in a BNodeMap. The parameter should
   365     ///be a write BNodeMap of UEdge values.
   366     ///\return The number of the matching edges.
   367     template<class MatchingMap>
   368     int bMatching(MatchingMap& mm) const {
   369       for (BNodeIt n(_g);n!=INVALID;++n) {
   370         mm.set(n,INVALID);
   371       }
   372       for (ANodeIt n(_g);n!=INVALID;++n) {
   373         if (_matching[n]!=INVALID)
   374 	  mm.set(_g.bNode(_matching[n]),_matching[n]);
   375       }
   376       return _matching_size;
   377     }
   378 
   379 
   380     ///Returns true if the given uedge is in the matching.
   381 
   382     ///It returns true if the given uedge is in the matching.
   383     ///
   384     bool matchingEdge(const UEdge& e) const {
   385       return _matching[_g.aNode(e)]==e;
   386     }
   387 
   388     ///Returns the matching edge from the node.
   389 
   390     ///Returns the matching edge from the node. If there is not such
   391     ///edge it gives back \c INVALID.  
   392     ///\note If the parameter node is a B-node then the running time is
   393     ///propotional to the degree of the node.
   394     UEdge matchingEdge(const Node& n) const {
   395       if (_g.aNode(n)) {
   396         return _matching[n];
   397       } else {
   398 	for (IncEdgeIt e(_g,n);e!=INVALID;++e)
   399 	  if (e==_matching[_g.aNode(e)]) return e;
   400 	return INVALID;
   401       }
   402     }
   403 
   404     ///Gives back the number of the matching edges.
   405 
   406     ///Gives back the number of the matching edges.
   407     int matchingSize() const {
   408       return _matching_size;
   409     }
   410 
   411     ///Gives back a barrier on the A-nodes
   412     
   413     ///The barrier is s subset of the nodes on the same side of the
   414     ///graph. If we tried to find a perfect matching and it failed
   415     ///then the barrier size will be greater than its neighbours. If
   416     ///the maximum matching searched then the barrier size minus its
   417     ///neighbours will be exactly the unmatched nodes on the A-side.
   418     ///\retval bar A WriteMap on the ANodes with bool value.
   419     template<class BarrierMap>
   420     void aBarrier(BarrierMap &bar) const 
   421     {
   422       for(ANodeIt n(_g);n!=INVALID;++n)
   423 	bar.set(n,_matching[n]==INVALID ||
   424 	  _levels[_g.bNode(_matching[n])]<_empty_level);  
   425     }  
   426 
   427     ///Gives back a barrier on the B-nodes
   428     
   429     ///The barrier is s subset of the nodes on the same side of the
   430     ///graph. If we tried to find a perfect matching and it failed
   431     ///then the barrier size will be greater than its neighbours. If
   432     ///the maximum matching searched then the barrier size minus its
   433     ///neighbours will be exactly the unmatched nodes on the B-side.
   434     ///\retval bar A WriteMap on the BNodes with bool value.
   435     template<class BarrierMap>
   436     void bBarrier(BarrierMap &bar) const
   437     {
   438       for(BNodeIt n(_g);n!=INVALID;++n) bar.set(n,_levels[n]>=_empty_level);  
   439     }
   440 
   441     ///Returns a minimum covering of the nodes.
   442 
   443     ///The minimum covering set problem is the dual solution of the
   444     ///maximum bipartite matching. It provides a solution for this
   445     ///problem what is proof of the optimality of the matching.
   446     ///\param covering NodeMap of bool values, the nodes of the cover
   447     ///set will set true while the others false.  
   448     ///\return The size of the cover set.
   449     ///\note This function can be called just after the algorithm have
   450     ///already found a matching. 
   451     template<class CoverMap>
   452     int coverSet(CoverMap& covering) const {
   453       int ret=0;
   454       for(BNodeIt n(_g);n!=INVALID;++n) {
   455 	if (_levels[n]<_empty_level) { covering.set(n,true); ++ret; }
   456 	else covering.set(n,false);
   457       }
   458       for(ANodeIt n(_g);n!=INVALID;++n) {
   459 	if (_matching[n]!=INVALID &&
   460 	    _levels[_g.bNode(_matching[n])]>=_empty_level) 
   461 	  { covering.set(n,true); ++ret; }
   462 	else covering.set(n,false);
   463       }
   464       return ret;
   465     }
   466 
   467 
   468     /// @}
   469     
   470   };
   471   
   472   
   473   ///Maximum cardinality of the matchings in a bipartite graph
   474 
   475   ///\ingroup matching
   476   ///This function finds the maximum cardinality of the matchings
   477   ///in a bipartite graph \c g.
   478   ///\param g An undirected bipartite graph.
   479   ///\return The cardinality of the maximum matching.
   480   ///
   481   ///\note The the implementation is based
   482   ///on the push-relabel principle.
   483   template<class Graph>
   484   int prBipartiteMatching(const Graph &g)
   485   {
   486     PrBipartiteMatching<Graph> bpm(g);
   487     return bpm.matchingSize();
   488   }
   489 
   490   ///Maximum cardinality matching in a bipartite graph
   491 
   492   ///\ingroup matching
   493   ///This function finds a maximum cardinality matching
   494   ///in a bipartite graph \c g.
   495   ///\param g An undirected bipartite graph.
   496   ///\retval matching A write ANodeMap of value type \c UEdge.
   497   /// The found edges will be returned in this map,
   498   /// i.e. for an \c ANode \c n the edge <tt>matching[n]</tt> is the one
   499   /// that covers the node \c n.
   500   ///\return The cardinality of the maximum matching.
   501   ///
   502   ///\note The the implementation is based
   503   ///on the push-relabel principle.
   504   template<class Graph,class MT>
   505   int prBipartiteMatching(const Graph &g,MT &matching) 
   506   {
   507     PrBipartiteMatching<Graph> bpm(g);
   508     bpm.run();
   509     bpm.aMatching(matching);
   510     return bpm.matchingSize();
   511   }
   512 
   513   ///Maximum cardinality matching in a bipartite graph
   514 
   515   ///\ingroup matching
   516   ///This function finds a maximum cardinality matching
   517   ///in a bipartite graph \c g.
   518   ///\param g An undirected bipartite graph.
   519   ///\retval matching A write ANodeMap of value type \c UEdge.
   520   /// The found edges will be returned in this map,
   521   /// i.e. for an \c ANode \c n the edge <tt>matching[n]</tt> is the one
   522   /// that covers the node \c n.
   523   ///\retval barrier A \c bool WriteMap on the BNodes. The map will be set
   524   /// exactly once for each BNode. The nodes with \c true value represent
   525   /// a barrier \e B, i.e. the cardinality of \e B minus the number of its
   526   /// neighbor is equal to the number of the <tt>BNode</tt>s minus the
   527   /// cardinality of the maximum matching.
   528   ///\return The cardinality of the maximum matching.
   529   ///
   530   ///\note The the implementation is based
   531   ///on the push-relabel principle.
   532   template<class Graph,class MT, class GT>
   533   int prBipartiteMatching(const Graph &g,MT &matching,GT &barrier) 
   534   {
   535     PrBipartiteMatching<Graph> bpm(g);
   536     bpm.run();
   537     bpm.aMatching(matching);
   538     bpm.bBarrier(barrier);
   539     return bpm.matchingSize();
   540   }  
   541 
   542   ///Perfect matching in a bipartite graph
   543 
   544   ///\ingroup matching
   545   ///This function checks whether the bipartite graph \c g
   546   ///has a perfect matching.
   547   ///\param g An undirected bipartite graph.
   548   ///\return \c true iff \c g has a perfect matching.
   549   ///
   550   ///\note The the implementation is based
   551   ///on the push-relabel principle.
   552   template<class Graph>
   553   bool prPerfectBipartiteMatching(const Graph &g)
   554   {
   555     PrBipartiteMatching<Graph> bpm(g);
   556     return bpm.runPerfect();
   557   }
   558 
   559   ///Perfect matching in a bipartite graph
   560 
   561   ///\ingroup matching
   562   ///This function finds a perfect matching in a bipartite graph \c g.
   563   ///\param g An undirected bipartite graph.
   564   ///\retval matching A write ANodeMap of value type \c UEdge.
   565   /// The found edges will be returned in this map,
   566   /// i.e. for an \c ANode \c n the edge <tt>matching[n]</tt> is the one
   567   /// that covers the node \c n.
   568   /// The values are unchanged if the graph
   569   /// has no perfect matching.
   570   ///\return \c true iff \c g has a perfect matching.
   571   ///
   572   ///\note The the implementation is based
   573   ///on the push-relabel principle.
   574   template<class Graph,class MT>
   575   bool prPerfectBipartiteMatching(const Graph &g,MT &matching) 
   576   {
   577     PrBipartiteMatching<Graph> bpm(g);
   578     bool ret = bpm.checkedRunPerfect();
   579     if (ret) bpm.aMatching(matching);
   580     return ret;
   581   }
   582 
   583   ///Perfect matching in a bipartite graph
   584 
   585   ///\ingroup matching
   586   ///This function finds a perfect matching in a bipartite graph \c g.
   587   ///\param g An undirected bipartite graph.
   588   ///\retval matching A write ANodeMap of value type \c UEdge.
   589   /// The found edges will be returned in this map,
   590   /// i.e. for an \c ANode \c n the edge <tt>matching[n]</tt> is the one
   591   /// that covers the node \c n.
   592   /// The values are unchanged if the graph
   593   /// has no perfect matching.
   594   ///\retval barrier A \c bool WriteMap on the BNodes. The map will only
   595   /// be set if \c g has no perfect matching. In this case it is set 
   596   /// exactly once for each BNode. The nodes with \c true value represent
   597   /// a barrier, i.e. a subset \e B a of BNodes with the property that
   598   /// the cardinality of \e B is greater than the number of its neighbors.
   599   ///\return \c true iff \c g has a perfect matching.
   600   ///
   601   ///\note The the implementation is based
   602   ///on the push-relabel principle.
   603   template<class Graph,class MT, class GT>
   604   bool prPerfectBipartiteMatching(const Graph &g,MT &matching,GT &barrier) 
   605   {
   606     PrBipartiteMatching<Graph> bpm(g);
   607     bool ret=bpm.checkedRunPerfect();
   608     if(ret)
   609       bpm.aMatching(matching);
   610     else
   611       bpm.bBarrier(barrier);
   612     return ret;
   613   }  
   614 }
   615 
   616 #endif