lemon/max_matching.h
author hegyi
Fri, 06 Jan 2006 13:58:49 +0000
changeset 1881 f40cdc2057c2
parent 1587 8f1c317ebeb4
child 1909 2d806130e700
permissions -rwxr-xr-x
Result of KruskalGUIAlgo is refreshed if displayed, but no more setin a forced way.
     1 /* -*- C++ -*-
     2  * lemon/max_matching.h - Part of LEMON, a generic C++ optimization library
     3  *
     4  * Copyright (C) 2006 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     5  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     6  *
     7  * Permission to use, modify and distribute this software is granted
     8  * provided that this copyright notice appears in all copies. For
     9  * precise terms see the accompanying LICENSE file.
    10  *
    11  * This software is provided "AS IS" with no warranty of any kind,
    12  * express or implied, and with no claim as to its suitability for any
    13  * purpose.
    14  *
    15  */
    16 
    17 #ifndef LEMON_MAX_MATCHING_H
    18 #define LEMON_MAX_MATCHING_H
    19 
    20 #include <queue>
    21 #include <lemon/invalid.h>
    22 #include <lemon/unionfind.h>
    23 #include <lemon/graph_utils.h>
    24 
    25 ///\ingroup galgs
    26 ///\file
    27 ///\brief Maximum matching algorithm.
    28 
    29 namespace lemon {
    30 
    31   /// \addtogroup galgs
    32   /// @{
    33 
    34   ///Edmonds' alternating forest maximum matching algorithm.
    35 
    36   ///This class provides Edmonds' alternating forest matching
    37   ///algorithm. The starting matching (if any) can be passed to the
    38   ///algorithm using read-in functions \ref readNMapNode, \ref
    39   ///readNMapEdge or \ref readEMapBool depending on the container. The
    40   ///resulting maximum matching can be attained by write-out functions
    41   ///\ref writeNMapNode, \ref writeNMapEdge or \ref writeEMapBool
    42   ///depending on the preferred container. 
    43   ///
    44   ///The dual side of a matching is a map of the nodes to
    45   ///MaxMatching::pos_enum, having values D, A and C showing the
    46   ///Gallai-Edmonds decomposition of the graph. The nodes in D induce
    47   ///a graph with factor-critical components, the nodes in A form the
    48   ///barrier, and the nodes in C induce a graph having a perfect
    49   ///matching. This decomposition can be attained by calling \ref
    50   ///writePos after running the algorithm. 
    51   ///
    52   ///\param Graph The undirected graph type the algorithm runs on.
    53   ///
    54   ///\author Jacint Szabo  
    55   template <typename Graph>
    56   class MaxMatching {
    57 
    58   protected:
    59 
    60     typedef typename Graph::Node Node;
    61     typedef typename Graph::Edge Edge;
    62     typedef typename Graph::UndirEdge UndirEdge;
    63     typedef typename Graph::UndirEdgeIt UndirEdgeIt;
    64     typedef typename Graph::NodeIt NodeIt;
    65     typedef typename Graph::IncEdgeIt IncEdgeIt;
    66 
    67     typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
    68 
    69   public:
    70     
    71     ///Indicates the Gallai-Edmonds decomposition of the graph.
    72 
    73     ///Indicates the Gallai-Edmonds decomposition of the graph, which
    74     ///shows an upper bound on the size of a maximum matching. The
    75     ///nodes with pos_enum \c D induce a graph with factor-critical
    76     ///components, the nodes in \c A form the canonical barrier, and the
    77     ///nodes in \c C induce a graph having a perfect matching. 
    78     enum pos_enum {
    79       D=0,
    80       A=1,
    81       C=2
    82     }; 
    83 
    84   protected:
    85 
    86     static const int HEUR_density=2;
    87     const Graph& g;
    88     typename Graph::template NodeMap<Node> _mate;
    89     typename Graph::template NodeMap<pos_enum> position;
    90      
    91   public:
    92     
    93     MaxMatching(const Graph& _g) : g(_g), _mate(_g,INVALID), position(_g) {}
    94 
    95     ///Runs Edmonds' algorithm.
    96 
    97     ///Runs Edmonds' algorithm for sparse graphs (number of edges <
    98     ///2*number of nodes), and a heuristical Edmonds' algorithm with a
    99     ///heuristic of postponing shrinks for dense graphs. 
   100     void run() {
   101       if ( countUndirEdges(g) < HEUR_density*countNodes(g) ) {
   102 	greedyMatching();
   103 	runEdmonds(0);
   104       } else runEdmonds(1);
   105     }
   106 
   107 
   108     ///Runs Edmonds' algorithm.
   109     
   110     ///If heur=0 it runs Edmonds' algorithm. If heur=1 it runs
   111     ///Edmonds' algorithm with a heuristic of postponing shrinks,
   112     ///giving a faster algorithm for dense graphs.  
   113     void runEdmonds( int heur = 1 ) {
   114       
   115       for(NodeIt v(g); v!=INVALID; ++v)
   116 	position.set(v,C);      
   117       
   118       typename Graph::template NodeMap<Node> ear(g,INVALID); 
   119       //undefined for the base nodes of the blossoms (i.e. for the
   120       //representative elements of UFE blossom) and for the nodes in C 
   121       
   122       typename UFE::MapType blossom_base(g);
   123       UFE blossom(blossom_base);
   124       typename UFE::MapType tree_base(g);
   125       UFE tree(tree_base);
   126       //If these UFE's would be members of the class then also
   127       //blossom_base and tree_base should be a member.
   128       
   129       for(NodeIt v(g); v!=INVALID; ++v) {
   130 	if ( position[v]==C && _mate[v]==INVALID ) {
   131 	  blossom.insert(v);
   132 	  tree.insert(v); 
   133 	  position.set(v,D);
   134 	  if ( heur == 1 ) lateShrink( v, ear, blossom, tree );
   135 	  else normShrink( v, ear, blossom, tree );
   136 	}
   137       }
   138     }
   139 
   140 
   141     ///Finds a greedy matching starting from the actual matching.
   142     
   143     ///Starting form the actual matching stored, it finds a maximal
   144     ///greedy matching.
   145     void greedyMatching() {
   146       for(NodeIt v(g); v!=INVALID; ++v)
   147 	if ( _mate[v]==INVALID ) {
   148 	  for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
   149 	    Node y=g.runningNode(e);
   150 	    if ( _mate[y]==INVALID && y!=v ) {
   151 	      _mate.set(v,y);
   152 	      _mate.set(y,v);
   153 	      break;
   154 	    }
   155 	  }
   156 	} 
   157     }
   158 
   159     ///Returns the size of the actual matching stored.
   160 
   161     ///Returns the size of the actual matching stored. After \ref
   162     ///run() it returns the size of a maximum matching in the graph.
   163     int size() const {
   164       int s=0;
   165       for(NodeIt v(g); v!=INVALID; ++v) {
   166 	if ( _mate[v]!=INVALID ) {
   167 	  ++s;
   168 	}
   169       }
   170       return s/2;
   171     }
   172 
   173 
   174     ///Resets the actual matching to the empty matching.
   175 
   176     ///Resets the actual matching to the empty matching.  
   177     ///
   178     void resetMatching() {
   179       for(NodeIt v(g); v!=INVALID; ++v)
   180 	_mate.set(v,INVALID);      
   181     }
   182 
   183     ///Returns the mate of a node in the actual matching.
   184 
   185     ///Returns the mate of a \c node in the actual matching. 
   186     ///Returns INVALID if the \c node is not covered by the actual matching. 
   187     Node mate(Node& node) const {
   188       return _mate[node];
   189     } 
   190 
   191     ///Reads a matching from a \c Node valued \c Node map.
   192 
   193     ///Reads a matching from a \c Node valued \c Node map. This map
   194     ///must be \e symmetric, i.e. if \c map[u]==v then \c map[v]==u
   195     ///must hold, and \c uv will be an edge of the matching.
   196     template<typename NMapN>
   197     void readNMapNode(NMapN& map) {
   198       for(NodeIt v(g); v!=INVALID; ++v) {
   199 	_mate.set(v,map[v]);   
   200       } 
   201     } 
   202     
   203     ///Writes the stored matching to a \c Node valued \c Node map.
   204 
   205     ///Writes the stored matching to a \c Node valued \c Node map. The
   206     ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
   207     ///map[v]==u will hold, and now \c uv is an edge of the matching.
   208     template<typename NMapN>
   209     void writeNMapNode (NMapN& map) const {
   210       for(NodeIt v(g); v!=INVALID; ++v) {
   211 	map.set(v,_mate[v]);   
   212       } 
   213     } 
   214 
   215     ///Reads a matching from an \c UndirEdge valued \c Node map.
   216 
   217     ///Reads a matching from an \c UndirEdge valued \c Node map. \c
   218     ///map[v] must be an \c UndirEdge incident to \c v. This map must
   219     ///have the property that if \c g.oppositeNode(u,map[u])==v then
   220     ///\c \c g.oppositeNode(v,map[v])==u holds, and now some edge
   221     ///joining \c u to \c v will be an edge of the matching.
   222     template<typename NMapE>
   223     void readNMapEdge(NMapE& map) {
   224      for(NodeIt v(g); v!=INVALID; ++v) {
   225        UndirEdge e=map[v];
   226 	if ( e!=INVALID )
   227 	  _mate.set(v,g.oppositeNode(v,e));
   228       } 
   229     } 
   230     
   231     ///Writes the matching stored to an \c UndirEdge valued \c Node map.
   232 
   233     ///Writes the stored matching to an \c UndirEdge valued \c Node
   234     ///map. \c map[v] will be an \c UndirEdge incident to \c v. This
   235     ///map will have the property that if \c g.oppositeNode(u,map[u])
   236     ///== v then \c map[u]==map[v] holds, and now this edge is an edge
   237     ///of the matching.
   238     template<typename NMapE>
   239     void writeNMapEdge (NMapE& map)  const {
   240       typename Graph::template NodeMap<bool> todo(g,true); 
   241       for(NodeIt v(g); v!=INVALID; ++v) {
   242 	if ( todo[v] && _mate[v]!=INVALID ) {
   243 	  Node u=_mate[v];
   244 	  for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
   245 	    if ( g.runningNode(e) == u ) {
   246 	      map.set(u,e);
   247 	      map.set(v,e);
   248 	      todo.set(u,false);
   249 	      todo.set(v,false);
   250 	      break;
   251 	    }
   252 	  }
   253 	}
   254       } 
   255     }
   256 
   257 
   258     ///Reads a matching from a \c bool valued \c Edge map.
   259     
   260     ///Reads a matching from a \c bool valued \c Edge map. This map
   261     ///must have the property that there are no two incident edges \c
   262     ///e, \c f with \c map[e]==map[f]==true. The edges \c e with \c
   263     ///map[e]==true form the matching.
   264     template<typename EMapB>
   265     void readEMapBool(EMapB& map) {
   266       for(UndirEdgeIt e(g); e!=INVALID; ++e) {
   267 	if ( map[e] ) {
   268 	  Node u=g.source(e);	  
   269 	  Node v=g.target(e);
   270 	  _mate.set(u,v);
   271 	  _mate.set(v,u);
   272 	} 
   273       } 
   274     }
   275 
   276 
   277     ///Writes the matching stored to a \c bool valued \c Edge map.
   278 
   279     ///Writes the matching stored to a \c bool valued \c Edge
   280     ///map. This map will have the property that there are no two
   281     ///incident edges \c e, \c f with \c map[e]==map[f]==true. The
   282     ///edges \c e with \c map[e]==true form the matching.
   283     template<typename EMapB>
   284     void writeEMapBool (EMapB& map) const {
   285       for(UndirEdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
   286 
   287       typename Graph::template NodeMap<bool> todo(g,true); 
   288       for(NodeIt v(g); v!=INVALID; ++v) {
   289 	if ( todo[v] && _mate[v]!=INVALID ) {
   290 	  Node u=_mate[v];
   291 	  for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
   292 	    if ( g.runningNode(e) == u ) {
   293 	      map.set(e,true);
   294 	      todo.set(u,false);
   295 	      todo.set(v,false);
   296 	      break;
   297 	    }
   298 	  }
   299 	}
   300       } 
   301     }
   302 
   303 
   304     ///Writes the canonical decomposition of the graph after running
   305     ///the algorithm.
   306 
   307     ///After calling any run methods of the class, it writes the
   308     ///Gallai-Edmonds canonical decomposition of the graph. \c map
   309     ///must be a node map of \ref pos_enum 's.
   310     template<typename NMapEnum>
   311     void writePos (NMapEnum& map) const {
   312       for(NodeIt v(g); v!=INVALID; ++v)  map.set(v,position[v]);
   313     }
   314 
   315   private: 
   316 
   317  
   318     void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,  
   319 		    UFE& blossom, UFE& tree);
   320 
   321     void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,  
   322 		    UFE& blossom, UFE& tree);
   323 
   324     bool noShrinkStep(Node x, typename Graph::template NodeMap<Node>& ear,  
   325 		      UFE& blossom, UFE& tree, std::queue<Node>& Q);
   326 
   327     void shrinkStep(Node& top, Node& middle, Node& bottom,
   328 		    typename Graph::template NodeMap<Node>& ear,  
   329 		    UFE& blossom, UFE& tree, std::queue<Node>& Q);
   330 
   331     void augment(Node x, typename Graph::template NodeMap<Node>& ear,  
   332 		 UFE& blossom, UFE& tree);
   333 
   334   };
   335 
   336 
   337   // **********************************************************************
   338   //  IMPLEMENTATIONS
   339   // **********************************************************************
   340 
   341 
   342   template <typename Graph>
   343   void MaxMatching<Graph>::lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,  
   344 				      UFE& blossom, UFE& tree) {
   345 
   346     std::queue<Node> Q;   //queue of the totally unscanned nodes
   347     Q.push(v);  
   348     std::queue<Node> R;   
   349     //queue of the nodes which must be scanned for a possible shrink
   350       
   351     while ( !Q.empty() ) {
   352       Node x=Q.front();
   353       Q.pop();
   354       if ( noShrinkStep( x, ear, blossom, tree, Q ) ) return;
   355       else R.push(x);
   356     }
   357       
   358     while ( !R.empty() ) {
   359       Node x=R.front();
   360       R.pop();
   361 	
   362       for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
   363 	Node y=g.runningNode(e);
   364 
   365 	if ( position[y] == D && blossom.find(x) != blossom.find(y) ) { 
   366 	  //x and y must be in the same tree
   367 	
   368 	  typename Graph::template NodeMap<bool> path(g,false);
   369 
   370 	  Node b=blossom.find(x);
   371 	  path.set(b,true);
   372 	  b=_mate[b];
   373 	  while ( b!=INVALID ) { 
   374 	    b=blossom.find(ear[b]);
   375 	    path.set(b,true);
   376 	    b=_mate[b];
   377 	  } //going till the root
   378 	
   379 	  Node top=y;
   380 	  Node middle=blossom.find(top);
   381 	  Node bottom=x;
   382 	  while ( !path[middle] )
   383 	    shrinkStep(top, middle, bottom, ear, blossom, tree, Q);
   384 		  
   385 	  Node base=middle;
   386 	  top=x;
   387 	  middle=blossom.find(top);
   388 	  bottom=y;
   389 	  Node blossom_base=blossom.find(base);
   390 	  while ( middle!=blossom_base )
   391 	    shrinkStep(top, middle, bottom, ear, blossom, tree, Q);
   392 		  
   393 	  blossom.makeRep(base);
   394 	} // if shrink is needed
   395 
   396 	while ( !Q.empty() ) {
   397 	  Node x=Q.front();
   398 	  Q.pop();
   399 	  if ( noShrinkStep(x, ear, blossom, tree, Q) ) return;
   400 	  else R.push(x);
   401 	}
   402       } //for e
   403     } // while ( !R.empty() )
   404   }
   405 
   406 
   407   template <typename Graph>
   408   void MaxMatching<Graph>::normShrink(Node v,
   409 				      typename Graph::template
   410 				      NodeMap<Node>& ear,  
   411 				      UFE& blossom, UFE& tree) {
   412     std::queue<Node> Q;   //queue of the unscanned nodes
   413     Q.push(v);  
   414     while ( !Q.empty() ) {
   415 
   416       Node x=Q.front();
   417       Q.pop();
   418 	
   419       for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
   420 	Node y=g.runningNode(e);
   421 	      
   422 	switch ( position[y] ) {
   423 	case D:          //x and y must be in the same tree
   424 
   425 	  if ( blossom.find(x) != blossom.find(y) ) { //shrink
   426 	    typename Graph::template NodeMap<bool> path(g,false);
   427 	      
   428 	    Node b=blossom.find(x);
   429 	    path.set(b,true);
   430 	    b=_mate[b];
   431 	    while ( b!=INVALID ) { 
   432 	      b=blossom.find(ear[b]);
   433 	      path.set(b,true);
   434 	      b=_mate[b];
   435 	    } //going till the root
   436 	
   437 	    Node top=y;
   438 	    Node middle=blossom.find(top);
   439 	    Node bottom=x;
   440 	    while ( !path[middle] )
   441 	      shrinkStep(top, middle, bottom, ear, blossom, tree, Q);
   442 		
   443 	    Node base=middle;
   444 	    top=x;
   445 	    middle=blossom.find(top);
   446 	    bottom=y;
   447 	    Node blossom_base=blossom.find(base);
   448 	    while ( middle!=blossom_base )
   449 	      shrinkStep(top, middle, bottom, ear, blossom, tree, Q);
   450 		
   451 	    blossom.makeRep(base);
   452 	  }
   453 	  break;
   454 	case C:
   455 	  if ( _mate[y]!=INVALID ) {   //grow
   456 
   457 	    ear.set(y,x);
   458 	    Node w=_mate[y];
   459 	    blossom.insert(w);
   460 	    position.set(y,A); 
   461 	    position.set(w,D); 
   462 	    tree.insert(y);
   463 	    tree.insert(w);
   464 	    tree.join(y,blossom.find(x));  
   465 	    tree.join(w,y);  
   466 	    Q.push(w);
   467 	  } else {                 //augment  
   468 	    augment(x, ear, blossom, tree);
   469 	    _mate.set(x,y);
   470 	    _mate.set(y,x);
   471 	    return;
   472 	  } //if 
   473 	  break;
   474 	default: break;
   475 	}
   476       }
   477     }
   478   }
   479 
   480   template <typename Graph>
   481   bool MaxMatching<Graph>::noShrinkStep(Node x,
   482 					typename Graph::template 
   483 					NodeMap<Node>& ear,  
   484 					UFE& blossom, UFE& tree,
   485 					std::queue<Node>& Q) {
   486     for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
   487       Node y=g.runningNode(e);
   488 	
   489       if ( position[y]==C ) {
   490 	if ( _mate[y]!=INVALID ) {       //grow
   491 	  ear.set(y,x);
   492 	  Node w=_mate[y];
   493 	  blossom.insert(w);
   494 	  position.set(y,A);
   495 	  position.set(w,D);
   496 	  tree.insert(y);
   497 	  tree.insert(w);
   498 	  tree.join(y,blossom.find(x));  
   499 	  tree.join(w,y);  
   500 	  Q.push(w);
   501 	} else {                      //augment 
   502 	  augment(x, ear, blossom, tree);
   503 	  _mate.set(x,y);
   504 	  _mate.set(y,x);
   505 	  return true;
   506 	}
   507       }
   508     }
   509     return false;
   510   }
   511 
   512   template <typename Graph>
   513   void MaxMatching<Graph>::shrinkStep(Node& top, Node& middle, Node& bottom,
   514 				      typename Graph::template
   515 				      NodeMap<Node>& ear,  
   516 				      UFE& blossom, UFE& tree,
   517 				      std::queue<Node>& Q) {
   518     ear.set(top,bottom);
   519     Node t=top;
   520     while ( t!=middle ) {
   521       Node u=_mate[t];
   522       t=ear[u];
   523       ear.set(t,u);
   524     } 
   525     bottom=_mate[middle];
   526     position.set(bottom,D);
   527     Q.push(bottom);
   528     top=ear[bottom];		
   529     Node oldmiddle=middle;
   530     middle=blossom.find(top);
   531     tree.erase(bottom);
   532     tree.erase(oldmiddle);
   533     blossom.insert(bottom);
   534     blossom.join(bottom, oldmiddle);
   535     blossom.join(top, oldmiddle);
   536   }
   537 
   538   template <typename Graph>
   539   void MaxMatching<Graph>::augment(Node x,
   540 				   typename Graph::template NodeMap<Node>& ear,  
   541 				   UFE& blossom, UFE& tree) { 
   542     Node v=_mate[x];
   543     while ( v!=INVALID ) {
   544 	
   545       Node u=ear[v];
   546       _mate.set(v,u);
   547       Node tmp=v;
   548       v=_mate[u];
   549       _mate.set(u,tmp);
   550     }
   551     typename UFE::ItemIt it;
   552     for (tree.first(it,blossom.find(x)); tree.valid(it); tree.next(it)) {   
   553       if ( position[it] == D ) {
   554 	typename UFE::ItemIt b_it;
   555 	for (blossom.first(b_it,it); blossom.valid(b_it); blossom.next(b_it)) {  
   556 	  position.set( b_it ,C);
   557 	}
   558 	blossom.eraseClass(it);
   559       } else position.set( it ,C);
   560     }
   561     tree.eraseClass(x);
   562 
   563   }
   564 
   565   /// @}
   566   
   567 } //END OF NAMESPACE LEMON
   568 
   569 #endif //LEMON_MAX_MATCHING_H