lemon/max_matching.h
changeset 2538 7bdd328de87a
parent 2391 14a343be7a5a
child 2548 a3ba22ebccc6
equal deleted inserted replaced
11:e0f9529fb650 12:5fe37616f660
    28 ///\file
    28 ///\file
    29 ///\brief Maximum matching algorithm in undirected graph.
    29 ///\brief Maximum matching algorithm in undirected graph.
    30 
    30 
    31 namespace lemon {
    31 namespace lemon {
    32 
    32 
    33   /// \ingroup matching
    33   ///\ingroup matching
    34 
    34 
    35   ///Edmonds' alternating forest maximum matching algorithm.
    35   ///\brief Edmonds' alternating forest maximum matching algorithm.
    36 
    36   ///
    37   ///This class provides Edmonds' alternating forest matching
    37   ///This class provides Edmonds' alternating forest matching
    38   ///algorithm. The starting matching (if any) can be passed to the
    38   ///algorithm. The starting matching (if any) can be passed to the
    39   ///algorithm using read-in functions \ref readNMapNode, \ref
    39   ///algorithm using some of init functions.
    40   ///readNMapEdge or \ref readEMapBool depending on the container. The
       
    41   ///resulting maximum matching can be attained by write-out functions
       
    42   ///\ref writeNMapNode, \ref writeNMapEdge or \ref writeEMapBool
       
    43   ///depending on the preferred container. 
       
    44   ///
    40   ///
    45   ///The dual side of a matching is a map of the nodes to
    41   ///The dual side of a matching is a map of the nodes to
    46   ///MaxMatching::pos_enum, having values D, A and C showing the
    42   ///MaxMatching::DecompType, having values \c D, \c A and \c C
    47   ///Gallai-Edmonds decomposition of the graph. The nodes in D induce
    43   ///showing the Gallai-Edmonds decomposition of the graph. The nodes
    48   ///a graph with factor-critical components, the nodes in A form the
    44   ///in \c D induce a graph with factor-critical components, the nodes
    49   ///barrier, and the nodes in C induce a graph having a perfect
    45   ///in \c A form the barrier, and the nodes in \c C induce a graph
    50   ///matching. This decomposition can be attained by calling \ref
    46   ///having a perfect matching. This decomposition can be attained by
    51   ///writePos after running the algorithm. 
    47   ///calling \c decomposition() after running the algorithm.
    52   ///
    48   ///
    53   ///\param Graph The undirected graph type the algorithm runs on.
    49   ///\param Graph The undirected graph type the algorithm runs on.
    54   ///
    50   ///
    55   ///\author Jacint Szabo  
    51   ///\author Jacint Szabo  
    56   template <typename Graph>
    52   template <typename Graph>
    65     typedef typename Graph::NodeIt NodeIt;
    61     typedef typename Graph::NodeIt NodeIt;
    66     typedef typename Graph::IncEdgeIt IncEdgeIt;
    62     typedef typename Graph::IncEdgeIt IncEdgeIt;
    67 
    63 
    68     typedef typename Graph::template NodeMap<int> UFECrossRef;
    64     typedef typename Graph::template NodeMap<int> UFECrossRef;
    69     typedef UnionFindEnum<UFECrossRef> UFE;
    65     typedef UnionFindEnum<UFECrossRef> UFE;
       
    66     typedef std::vector<Node> NV;
       
    67 
       
    68     typedef typename Graph::template NodeMap<int> EFECrossRef;
       
    69     typedef ExtendFindEnum<EFECrossRef> EFE;
    70 
    70 
    71   public:
    71   public:
    72     
    72     
    73     ///Indicates the Gallai-Edmonds decomposition of the graph.
    73     ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
    74 
    74     ///
    75     ///Indicates the Gallai-Edmonds decomposition of the graph, which
    75     ///Indicates the Gallai-Edmonds decomposition of the graph, which
    76     ///shows an upper bound on the size of a maximum matching. The
    76     ///shows an upper bound on the size of a maximum matching. The
    77     ///nodes with pos_enum \c D induce a graph with factor-critical
    77     ///nodes with DecompType \c D induce a graph with factor-critical
    78     ///components, the nodes in \c A form the canonical barrier, and the
    78     ///components, the nodes in \c A form the canonical barrier, and the
    79     ///nodes in \c C induce a graph having a perfect matching. 
    79     ///nodes in \c C induce a graph having a perfect matching. 
    80     enum pos_enum {
    80     enum DecompType {
    81       D=0,
    81       D=0,
    82       A=1,
    82       A=1,
    83       C=2
    83       C=2
    84     }; 
    84     }; 
    85 
    85 
    86   protected:
    86   protected:
    87 
    87 
    88     static const int HEUR_density=2;
    88     static const int HEUR_density=2;
    89     const Graph& g;
    89     const Graph& g;
    90     typename Graph::template NodeMap<Node> _mate;
    90     typename Graph::template NodeMap<Node> _mate;
    91     typename Graph::template NodeMap<pos_enum> position;
    91     typename Graph::template NodeMap<DecompType> position;
    92      
    92      
    93   public:
    93   public:
    94     
    94     
    95     MaxMatching(const Graph& _g) : g(_g), _mate(_g,INVALID), position(_g) {}
    95     MaxMatching(const Graph& _g) 
    96 
    96       : g(_g), _mate(_g), position(_g) {}
    97     ///Runs Edmonds' algorithm.
    97 
    98 
    98     ///\brief Sets the actual matching to the empty matching.
    99     ///Runs Edmonds' algorithm for sparse graphs (number of edges <
    99     ///
   100     ///2*number of nodes), and a heuristical Edmonds' algorithm with a
   100     ///Sets the actual matching to the empty matching.  
   101     ///heuristic of postponing shrinks for dense graphs. 
   101     ///
   102     void run() {
   102     void init() {
   103       if ( countUEdges(g) < HEUR_density*countNodes(g) ) {
   103       for(NodeIt v(g); v!=INVALID; ++v) {
   104 	greedyMatching();
   104 	_mate.set(v,INVALID);      
   105 	runEdmonds(0);
   105 	position.set(v,C);
   106       } else runEdmonds(1);
   106       }
   107     }
   107     }
   108 
   108 
   109 
   109     ///\brief Finds a greedy matching for initial matching.
   110     ///Runs Edmonds' algorithm.
   110     ///
   111     
   111     ///For initial matchig it finds a maximal greedy matching.
   112     ///If heur=0 it runs Edmonds' algorithm. If heur=1 it runs
   112     void greedyInit() {
   113     ///Edmonds' algorithm with a heuristic of postponing shrinks,
   113       for(NodeIt v(g); v!=INVALID; ++v) {
   114     ///giving a faster algorithm for dense graphs.  
   114 	_mate.set(v,INVALID);      
   115     void runEdmonds( int heur = 1 ) {
   115 	position.set(v,C);
   116       
   116       }
   117       //each vertex is put to C
       
   118       for(NodeIt v(g); v!=INVALID; ++v)
       
   119 	position.set(v,C);      
       
   120       
       
   121       typename Graph::template NodeMap<Node> ear(g,INVALID); 
       
   122       //undefined for the base nodes of the blossoms (i.e. for the
       
   123       //representative elements of UFE blossom) and for the nodes in C 
       
   124       
       
   125       UFECrossRef blossom_base(g);
       
   126       UFE blossom(blossom_base);
       
   127 
       
   128       UFECrossRef tree_base(g);
       
   129       UFE tree(tree_base);
       
   130 
       
   131       //If these UFE's would be members of the class then also
       
   132       //blossom_base and tree_base should be a member.
       
   133       
       
   134       //We build only one tree and the other vertices uncovered by the
       
   135       //matching belong to C. (They can be considered as singleton
       
   136       //trees.) If this tree can be augmented or no more
       
   137       //grow/augmentation/shrink is possible then we return to this
       
   138       //"for" cycle.
       
   139       for(NodeIt v(g); v!=INVALID; ++v) {
       
   140 	if ( position[v]==C && _mate[v]==INVALID ) {
       
   141 	  blossom.insert(v);
       
   142 	  tree.insert(v); 
       
   143 	  position.set(v,D);
       
   144 	  if ( heur == 1 ) lateShrink( v, ear, blossom, tree );
       
   145 	  else normShrink( v, ear, blossom, tree );
       
   146 	}
       
   147       }
       
   148     }
       
   149 
       
   150 
       
   151     ///Finds a greedy matching starting from the actual matching.
       
   152     
       
   153     ///Starting form the actual matching stored, it finds a maximal
       
   154     ///greedy matching.
       
   155     void greedyMatching() {
       
   156       for(NodeIt v(g); v!=INVALID; ++v)
   117       for(NodeIt v(g); v!=INVALID; ++v)
   157 	if ( _mate[v]==INVALID ) {
   118 	if ( _mate[v]==INVALID ) {
   158 	  for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
   119 	  for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
   159 	    Node y=g.runningNode(e);
   120 	    Node y=g.runningNode(e);
   160 	    if ( _mate[y]==INVALID && y!=v ) {
   121 	    if ( _mate[y]==INVALID && y!=v ) {
   164 	    }
   125 	    }
   165 	  }
   126 	  }
   166 	} 
   127 	} 
   167     }
   128     }
   168 
   129 
   169     ///Returns the size of the actual matching stored.
   130     ///\brief Initialize the matching from each nodes' mate.
   170 
   131     ///
       
   132     ///Initialize the matching from a \c Node valued \c Node map. This
       
   133     ///map must be \e symmetric, i.e. if \c map[u]==v then \c
       
   134     ///map[v]==u must hold, and \c uv will be an edge of the initial
       
   135     ///matching.
       
   136     template <typename MateMap>
       
   137     void mateMapInit(MateMap& map) {
       
   138       for(NodeIt v(g); v!=INVALID; ++v) {
       
   139 	_mate.set(v,map[v]);
       
   140 	position.set(v,C);
       
   141       } 
       
   142     }
       
   143 
       
   144     ///\brief Initialize the matching from a node map with the
       
   145     ///incident matching edges.
       
   146     ///
       
   147     ///Initialize the matching from an \c UEdge valued \c Node map. \c
       
   148     ///map[v] must be an \c UEdge incident to \c v. This map must have
       
   149     ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
       
   150     ///g.oppositeNode(v,map[v])==u holds, and now some edge joining \c
       
   151     ///u to \c v will be an edge of the matching.
       
   152     template<typename MatchingMap>
       
   153     void matchingMapInit(MatchingMap& map) {
       
   154       for(NodeIt v(g); v!=INVALID; ++v) {
       
   155 	position.set(v,C);
       
   156 	UEdge e=map[v];
       
   157 	if ( e!=INVALID )
       
   158 	  _mate.set(v,g.oppositeNode(v,e));
       
   159 	else 
       
   160 	  _mate.set(v,INVALID);
       
   161       } 
       
   162     } 
       
   163 
       
   164     ///\brief Initialize the matching from the map containing the
       
   165     ///undirected matching edges.
       
   166     ///
       
   167     ///Initialize the matching from a \c bool valued \c UEdge map. This
       
   168     ///map must have the property that there are no two incident edges
       
   169     ///\c e, \c f with \c map[e]==map[f]==true. The edges \c e with \c
       
   170     ///map[e]==true form the matching.
       
   171     template <typename MatchingMap>
       
   172     void matchingInit(MatchingMap& map) {
       
   173       for(NodeIt v(g); v!=INVALID; ++v) {
       
   174 	_mate.set(v,INVALID);      
       
   175 	position.set(v,C);
       
   176       }
       
   177       for(UEdgeIt e(g); e!=INVALID; ++e) {
       
   178 	if ( map[e] ) {
       
   179 	  Node u=g.source(e);	  
       
   180 	  Node v=g.target(e);
       
   181 	  _mate.set(u,v);
       
   182 	  _mate.set(v,u);
       
   183 	} 
       
   184       } 
       
   185     }
       
   186 
       
   187 
       
   188     ///\brief Runs Edmonds' algorithm.
       
   189     ///
       
   190     ///Runs Edmonds' algorithm for sparse graphs (number of edges <
       
   191     ///2*number of nodes), and a heuristical Edmonds' algorithm with a
       
   192     ///heuristic of postponing shrinks for dense graphs. 
       
   193     void run() {
       
   194       if (countUEdges(g) < HEUR_density * countNodes(g)) {
       
   195 	greedyInit();
       
   196 	startSparse();
       
   197       } else {
       
   198 	init();
       
   199 	startDense();
       
   200       }
       
   201     }
       
   202 
       
   203 
       
   204     ///\brief Starts Edmonds' algorithm.
       
   205     /// 
       
   206     ///If runs the original Edmonds' algorithm.  
       
   207     void startSparse() {
       
   208             
       
   209       typename Graph::template NodeMap<Node> ear(g,INVALID); 
       
   210       //undefined for the base nodes of the blossoms (i.e. for the
       
   211       //representative elements of UFE blossom) and for the nodes in C 
       
   212       
       
   213       UFECrossRef blossom_base(g);
       
   214       UFE blossom(blossom_base);
       
   215       NV rep(countNodes(g));
       
   216 
       
   217       EFECrossRef tree_base(g);
       
   218       EFE tree(tree_base);
       
   219 
       
   220       //If these UFE's would be members of the class then also
       
   221       //blossom_base and tree_base should be a member.
       
   222       
       
   223       //We build only one tree and the other vertices uncovered by the
       
   224       //matching belong to C. (They can be considered as singleton
       
   225       //trees.) If this tree can be augmented or no more
       
   226       //grow/augmentation/shrink is possible then we return to this
       
   227       //"for" cycle.
       
   228       for(NodeIt v(g); v!=INVALID; ++v) {
       
   229 	if (position[v]==C && _mate[v]==INVALID) {
       
   230 	  rep[blossom.insert(v)] = v;
       
   231 	  tree.insert(v); 
       
   232 	  position.set(v,D);
       
   233 	  normShrink(v, ear, blossom, rep, tree);
       
   234 	}
       
   235       }
       
   236     }
       
   237 
       
   238     ///\brief Starts Edmonds' algorithm.
       
   239     /// 
       
   240     ///It runs Edmonds' algorithm with a heuristic of postponing
       
   241     ///shrinks, giving a faster algorithm for dense graphs.
       
   242     void startDense() {
       
   243             
       
   244       typename Graph::template NodeMap<Node> ear(g,INVALID); 
       
   245       //undefined for the base nodes of the blossoms (i.e. for the
       
   246       //representative elements of UFE blossom) and for the nodes in C 
       
   247       
       
   248       UFECrossRef blossom_base(g);
       
   249       UFE blossom(blossom_base);
       
   250       NV rep(countNodes(g));
       
   251 
       
   252       EFECrossRef tree_base(g);
       
   253       EFE tree(tree_base);
       
   254 
       
   255       //If these UFE's would be members of the class then also
       
   256       //blossom_base and tree_base should be a member.
       
   257       
       
   258       //We build only one tree and the other vertices uncovered by the
       
   259       //matching belong to C. (They can be considered as singleton
       
   260       //trees.) If this tree can be augmented or no more
       
   261       //grow/augmentation/shrink is possible then we return to this
       
   262       //"for" cycle.
       
   263       for(NodeIt v(g); v!=INVALID; ++v) {
       
   264 	if ( position[v]==C && _mate[v]==INVALID ) {
       
   265 	  rep[blossom.insert(v)] = v;
       
   266 	  tree.insert(v); 
       
   267 	  position.set(v,D);
       
   268 	  lateShrink(v, ear, blossom, rep, tree);
       
   269 	}
       
   270       }
       
   271     }
       
   272 
       
   273 
       
   274 
       
   275     ///\brief Returns the size of the actual matching stored.
       
   276     ///
   171     ///Returns the size of the actual matching stored. After \ref
   277     ///Returns the size of the actual matching stored. After \ref
   172     ///run() it returns the size of a maximum matching in the graph.
   278     ///run() it returns the size of a maximum matching in the graph.
   173     int size() const {
   279     int size() const {
   174       int s=0;
   280       int s=0;
   175       for(NodeIt v(g); v!=INVALID; ++v) {
   281       for(NodeIt v(g); v!=INVALID; ++v) {
   179       }
   285       }
   180       return s/2;
   286       return s/2;
   181     }
   287     }
   182 
   288 
   183 
   289 
   184     ///Resets the actual matching to the empty matching.
   290     ///\brief Returns the mate of a node in the actual matching.
   185 
   291     ///
   186     ///Resets the actual matching to the empty matching.  
       
   187     ///
       
   188     void resetMatching() {
       
   189       for(NodeIt v(g); v!=INVALID; ++v)
       
   190 	_mate.set(v,INVALID);      
       
   191     }
       
   192 
       
   193     ///Returns the mate of a node in the actual matching.
       
   194 
       
   195     ///Returns the mate of a \c node in the actual matching. 
   292     ///Returns the mate of a \c node in the actual matching. 
   196     ///Returns INVALID if the \c node is not covered by the actual matching. 
   293     ///Returns INVALID if the \c node is not covered by the actual matching. 
   197     Node mate(Node& node) const {
   294     Node mate(const Node& node) const {
   198       return _mate[node];
   295       return _mate[node];
   199     } 
   296     } 
   200 
   297 
   201     ///Reads a matching from a \c Node valued \c Node map.
   298     ///\brief Returns the matching edge incident to the given node.
   202 
   299     ///
   203     ///Reads a matching from a \c Node valued \c Node map. This map
   300     ///Returns the matching edge of a \c node in the actual matching. 
   204     ///must be \e symmetric, i.e. if \c map[u]==v then \c map[v]==u
   301     ///Returns INVALID if the \c node is not covered by the actual matching. 
   205     ///must hold, and \c uv will be an edge of the matching.
   302     UEdge matchingEdge(const Node& node) const {
   206     template<typename NMapN>
   303       if (_mate[node] == INVALID) return INVALID;
   207     void readNMapNode(NMapN& map) {
   304       Node n = node < _mate[node] ? node : _mate[node];
   208       for(NodeIt v(g); v!=INVALID; ++v) {
   305       for (IncEdgeIt e(g, n); e != INVALID; ++e) {
   209 	_mate.set(v,map[v]);   
   306 	if (g.oppositeNode(n, e) == _mate[n]) {
   210       } 
   307 	  return e;
       
   308 	}
       
   309       }
       
   310       return INVALID;
   211     } 
   311     } 
   212     
   312 
   213     ///Writes the stored matching to a \c Node valued \c Node map.
   313     /// \brief Returns the class of the node in the Edmonds-Gallai
   214 
   314     /// decomposition.
       
   315     ///
       
   316     /// Returns the class of the node in the Edmonds-Gallai
       
   317     /// decomposition.
       
   318     DecompType decomposition(const Node& n) {
       
   319       return position[n] == A;
       
   320     }
       
   321 
       
   322     /// \brief Returns true when the node is in the barrier.
       
   323     ///
       
   324     /// Returns true when the node is in the barrier.
       
   325     bool barrier(const Node& n) {
       
   326       return position[n] == A;
       
   327     }
       
   328     
       
   329     ///\brief Gives back the matching in a \c Node of mates.
       
   330     ///
   215     ///Writes the stored matching to a \c Node valued \c Node map. The
   331     ///Writes the stored matching to a \c Node valued \c Node map. The
   216     ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
   332     ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
   217     ///map[v]==u will hold, and now \c uv is an edge of the matching.
   333     ///map[v]==u will hold, and now \c uv is an edge of the matching.
   218     template<typename NMapN>
   334     template <typename MateMap>
   219     void writeNMapNode (NMapN& map) const {
   335     void mateMap(MateMap& map) const {
   220       for(NodeIt v(g); v!=INVALID; ++v) {
   336       for(NodeIt v(g); v!=INVALID; ++v) {
   221 	map.set(v,_mate[v]);   
   337 	map.set(v,_mate[v]);   
   222       } 
   338       } 
   223     } 
   339     } 
   224 
   340     
   225     ///Reads a matching from an \c UEdge valued \c Node map.
   341     ///\brief Gives back the matching in an \c UEdge valued \c Node
   226 
   342     ///map.
   227     ///Reads a matching from an \c UEdge valued \c Node map. \c
   343     ///
   228     ///map[v] must be an \c UEdge incident to \c v. This map must
       
   229     ///have the property that if \c g.oppositeNode(u,map[u])==v then
       
   230     ///\c \c g.oppositeNode(v,map[v])==u holds, and now some edge
       
   231     ///joining \c u to \c v will be an edge of the matching.
       
   232     template<typename NMapE>
       
   233     void readNMapEdge(NMapE& map) {
       
   234       for(NodeIt v(g); v!=INVALID; ++v) {
       
   235 	UEdge e=map[v];
       
   236 	if ( e!=INVALID )
       
   237 	  _mate.set(v,g.oppositeNode(v,e));
       
   238       } 
       
   239     } 
       
   240     
       
   241     ///Writes the matching stored to an \c UEdge valued \c Node map.
       
   242 
       
   243     ///Writes the stored matching to an \c UEdge valued \c Node
   344     ///Writes the stored matching to an \c UEdge valued \c Node
   244     ///map. \c map[v] will be an \c UEdge incident to \c v. This
   345     ///map. \c map[v] will be an \c UEdge incident to \c v. This
   245     ///map will have the property that if \c g.oppositeNode(u,map[u])
   346     ///map will have the property that if \c g.oppositeNode(u,map[u])
   246     ///== v then \c map[u]==map[v] holds, and now this edge is an edge
   347     ///== v then \c map[u]==map[v] holds, and now this edge is an edge
   247     ///of the matching.
   348     ///of the matching.
   248     template<typename NMapE>
   349     template <typename MatchingMap>
   249     void writeNMapEdge (NMapE& map)  const {
   350     void matchingMap(MatchingMap& map)  const {
   250       typename Graph::template NodeMap<bool> todo(g,true); 
   351       typename Graph::template NodeMap<bool> todo(g,true); 
   251       for(NodeIt v(g); v!=INVALID; ++v) {
   352       for(NodeIt v(g); v!=INVALID; ++v) {
   252 	if ( todo[v] && _mate[v]!=INVALID ) {
   353 	if (_mate[v]!=INVALID && v < _mate[v]) {
   253 	  Node u=_mate[v];
   354 	  Node u=_mate[v];
   254 	  for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
   355 	  for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
   255 	    if ( g.runningNode(e) == u ) {
   356 	    if ( g.runningNode(e) == u ) {
   256 	      map.set(u,e);
   357 	      map.set(u,e);
   257 	      map.set(v,e);
   358 	      map.set(v,e);
   263 	}
   364 	}
   264       } 
   365       } 
   265     }
   366     }
   266 
   367 
   267 
   368 
   268     ///Reads a matching from a \c bool valued \c Edge map.
   369     ///\brief Gives back the matching in a \c bool valued \c UEdge
   269     
   370     ///map.
   270     ///Reads a matching from a \c bool valued \c Edge map. This map
   371     ///
   271     ///must have the property that there are no two incident edges \c
       
   272     ///e, \c f with \c map[e]==map[f]==true. The edges \c e with \c
       
   273     ///map[e]==true form the matching.
       
   274     template<typename EMapB>
       
   275     void readEMapBool(EMapB& map) {
       
   276       for(UEdgeIt e(g); e!=INVALID; ++e) {
       
   277 	if ( map[e] ) {
       
   278 	  Node u=g.source(e);	  
       
   279 	  Node v=g.target(e);
       
   280 	  _mate.set(u,v);
       
   281 	  _mate.set(v,u);
       
   282 	} 
       
   283       } 
       
   284     }
       
   285 
       
   286 
       
   287     ///Writes the matching stored to a \c bool valued \c Edge map.
       
   288 
       
   289     ///Writes the matching stored to a \c bool valued \c Edge
   372     ///Writes the matching stored to a \c bool valued \c Edge
   290     ///map. This map will have the property that there are no two
   373     ///map. This map will have the property that there are no two
   291     ///incident edges \c e, \c f with \c map[e]==map[f]==true. The
   374     ///incident edges \c e, \c f with \c map[e]==map[f]==true. The
   292     ///edges \c e with \c map[e]==true form the matching.
   375     ///edges \c e with \c map[e]==true form the matching.
   293     template<typename EMapB>
   376     template<typename MatchingMap>
   294     void writeEMapBool (EMapB& map) const {
   377     void matching(MatchingMap& map) const {
   295       for(UEdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
   378       for(UEdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
   296 
   379 
   297       typename Graph::template NodeMap<bool> todo(g,true); 
   380       typename Graph::template NodeMap<bool> todo(g,true); 
   298       for(NodeIt v(g); v!=INVALID; ++v) {
   381       for(NodeIt v(g); v!=INVALID; ++v) {
   299 	if ( todo[v] && _mate[v]!=INVALID ) {
   382 	if ( todo[v] && _mate[v]!=INVALID ) {
   309 	}
   392 	}
   310       } 
   393       } 
   311     }
   394     }
   312 
   395 
   313 
   396 
   314     ///Writes the canonical decomposition of the graph after running
   397     ///\brief Returns the canonical decomposition of the graph after running
   315     ///the algorithm.
   398     ///the algorithm.
   316 
   399     ///
   317     ///After calling any run methods of the class, it writes the
   400     ///After calling any run methods of the class, it writes the
   318     ///Gallai-Edmonds canonical decomposition of the graph. \c map
   401     ///Gallai-Edmonds canonical decomposition of the graph. \c map
   319     ///must be a node map of \ref pos_enum 's.
   402     ///must be a node map of \ref DecompType 's.
   320     template<typename NMapEnum>
   403     template <typename DecompositionMap>
   321     void writePos (NMapEnum& map) const {
   404     void decomposition(DecompositionMap& map) const {
   322       for(NodeIt v(g); v!=INVALID; ++v)  map.set(v,position[v]);
   405       for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
       
   406     }
       
   407 
       
   408     ///\brief Returns a barrier on the nodes.
       
   409     ///
       
   410     ///After calling any run methods of the class, it writes a
       
   411     ///canonical barrier on the nodes. The odd component number of the
       
   412     ///remaining graph minus the barrier size is a lower bound for the
       
   413     ///uncovered nodes in the graph. The \c map must be a node map of
       
   414     ///bools.
       
   415     template <typename BarrierMap>
       
   416     void barrier(BarrierMap& barrier) {
       
   417       for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);    
   323     }
   418     }
   324 
   419 
   325   private: 
   420   private: 
   326 
   421 
   327  
   422  
   328     void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,  
   423     void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,  
   329 		    UFE& blossom, UFE& tree);
   424 		    UFE& blossom, NV& rep, EFE& tree) {
       
   425       //We have one tree which we grow, and also shrink but only if it
       
   426       //cannot be postponed. If we augment then we return to the "for"
       
   427       //cycle of runEdmonds().
       
   428 
       
   429       std::queue<Node> Q;   //queue of the totally unscanned nodes
       
   430       Q.push(v);  
       
   431       std::queue<Node> R;   
       
   432       //queue of the nodes which must be scanned for a possible shrink
       
   433       
       
   434       while ( !Q.empty() ) {
       
   435 	Node x=Q.front();
       
   436 	Q.pop();
       
   437 	for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
       
   438 	  Node y=g.runningNode(e);
       
   439 	  //growOrAugment grows if y is covered by the matching and
       
   440 	  //augments if not. In this latter case it returns 1.
       
   441 	  if (position[y]==C && 
       
   442 	      growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
       
   443 	}
       
   444 	R.push(x);
       
   445       }
       
   446       
       
   447       while ( !R.empty() ) {
       
   448 	Node x=R.front();
       
   449 	R.pop();
       
   450 	
       
   451 	for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
       
   452 	  Node y=g.runningNode(e);
       
   453 
       
   454 	  if ( position[y] == D && blossom.find(x) != blossom.find(y) )  
       
   455 	    //Recall that we have only one tree.
       
   456 	    shrink( x, y, ear, blossom, rep, tree, Q);	
       
   457 	
       
   458 	  while ( !Q.empty() ) {
       
   459 	    Node z=Q.front();
       
   460 	    Q.pop();
       
   461 	    for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
       
   462 	      Node w=g.runningNode(f);
       
   463 	      //growOrAugment grows if y is covered by the matching and
       
   464 	      //augments if not. In this latter case it returns 1.
       
   465 	      if (position[w]==C && 
       
   466 		  growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
       
   467 	    }
       
   468 	    R.push(z);
       
   469 	  }
       
   470 	} //for e
       
   471       } // while ( !R.empty() )
       
   472     }
   330 
   473 
   331     void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,  
   474     void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,  
   332 		    UFE& blossom, UFE& tree);
   475 		    UFE& blossom, NV& rep, EFE& tree) {
       
   476       //We have one tree, which we grow and shrink. If we augment then we
       
   477       //return to the "for" cycle of runEdmonds().
       
   478     
       
   479       std::queue<Node> Q;   //queue of the unscanned nodes
       
   480       Q.push(v);  
       
   481       while ( !Q.empty() ) {
       
   482 
       
   483 	Node x=Q.front();
       
   484 	Q.pop();
       
   485 	
       
   486 	for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
       
   487 	  Node y=g.runningNode(e);
       
   488 	      
       
   489 	  switch ( position[y] ) {
       
   490 	  case D:          //x and y must be in the same tree
       
   491 	    if ( blossom.find(x) != blossom.find(y))
       
   492 	      //x and y are in the same tree
       
   493 	      shrink(x, y, ear, blossom, rep, tree, Q);
       
   494 	    break;
       
   495 	  case C:
       
   496 	    //growOrAugment grows if y is covered by the matching and
       
   497 	    //augments if not. In this latter case it returns 1.
       
   498 	    if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
       
   499 	    break;
       
   500 	  default: break;
       
   501 	  }
       
   502 	}
       
   503       }
       
   504     }
   333 
   505 
   334     void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,  
   506     void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,  
   335 		UFE& blossom, UFE& tree,std::queue<Node>& Q);
   507 		UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
       
   508       //x and y are the two adjacent vertices in two blossoms.
       
   509    
       
   510       typename Graph::template NodeMap<bool> path(g,false);
       
   511     
       
   512       Node b=rep[blossom.find(x)];
       
   513       path.set(b,true);
       
   514       b=_mate[b];
       
   515       while ( b!=INVALID ) { 
       
   516 	b=rep[blossom.find(ear[b])];
       
   517 	path.set(b,true);
       
   518 	b=_mate[b];
       
   519       } //we go until the root through bases of blossoms and odd vertices
       
   520     
       
   521       Node top=y;
       
   522       Node middle=rep[blossom.find(top)];
       
   523       Node bottom=x;
       
   524       while ( !path[middle] )
       
   525 	shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
       
   526       //Until we arrive to a node on the path, we update blossom, tree
       
   527       //and the positions of the odd nodes.
       
   528     
       
   529       Node base=middle;
       
   530       top=x;
       
   531       middle=rep[blossom.find(top)];
       
   532       bottom=y;
       
   533       Node blossom_base=rep[blossom.find(base)];
       
   534       while ( middle!=blossom_base )
       
   535 	shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
       
   536       //Until we arrive to a node on the path, we update blossom, tree
       
   537       //and the positions of the odd nodes.
       
   538     
       
   539       rep[blossom.find(base)] = base;
       
   540     }
   336 
   541 
   337     void shrinkStep(Node& top, Node& middle, Node& bottom,
   542     void shrinkStep(Node& top, Node& middle, Node& bottom,
   338 		    typename Graph::template NodeMap<Node>& ear,  
   543 		    typename Graph::template NodeMap<Node>& ear,  
   339 		    UFE& blossom, UFE& tree, std::queue<Node>& Q);
   544 		    UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
       
   545       //We traverse a blossom and update everything.
       
   546     
       
   547       ear.set(top,bottom);
       
   548       Node t=top;
       
   549       while ( t!=middle ) {
       
   550 	Node u=_mate[t];
       
   551 	t=ear[u];
       
   552 	ear.set(t,u);
       
   553       } 
       
   554       bottom=_mate[middle];
       
   555       position.set(bottom,D);
       
   556       Q.push(bottom);
       
   557       top=ear[bottom];		
       
   558       Node oldmiddle=middle;
       
   559       middle=rep[blossom.find(top)];
       
   560       tree.erase(bottom);
       
   561       tree.erase(oldmiddle);
       
   562       blossom.insert(bottom);
       
   563       blossom.join(bottom, oldmiddle);
       
   564       blossom.join(top, oldmiddle);
       
   565     }
       
   566 
       
   567 
   340 
   568 
   341     bool growOrAugment(Node& y, Node& x, typename Graph::template 
   569     bool growOrAugment(Node& y, Node& x, typename Graph::template 
   342 		       NodeMap<Node>& ear, UFE& blossom, UFE& tree, 
   570 		       NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree, 
   343 		       std::queue<Node>& Q);
   571 		       std::queue<Node>& Q) {
       
   572       //x is in a blossom in the tree, y is outside. If y is covered by
       
   573       //the matching we grow, otherwise we augment. In this case we
       
   574       //return 1.
       
   575     
       
   576       if ( _mate[y]!=INVALID ) {       //grow
       
   577 	ear.set(y,x);
       
   578 	Node w=_mate[y];
       
   579 	rep[blossom.insert(w)] = w;
       
   580 	position.set(y,A);
       
   581 	position.set(w,D);
       
   582 	int t = tree.find(rep[blossom.find(x)]);
       
   583 	tree.insert(y,t);  
       
   584 	tree.insert(w,t);  
       
   585 	Q.push(w);
       
   586       } else {                      //augment 
       
   587 	augment(x, ear, blossom, rep, tree);
       
   588 	_mate.set(x,y);
       
   589 	_mate.set(y,x);
       
   590 	return true;
       
   591       }
       
   592       return false;
       
   593     }
   344 
   594 
   345     void augment(Node x, typename Graph::template NodeMap<Node>& ear,  
   595     void augment(Node x, typename Graph::template NodeMap<Node>& ear,  
   346 		 UFE& blossom, UFE& tree);
   596 		 UFE& blossom, NV& rep, EFE& tree) {
       
   597       Node v=_mate[x];
       
   598       while ( v!=INVALID ) {
       
   599 	
       
   600 	Node u=ear[v];
       
   601 	_mate.set(v,u);
       
   602 	Node tmp=v;
       
   603 	v=_mate[u];
       
   604 	_mate.set(u,tmp);
       
   605       }
       
   606       int y = tree.find(rep[blossom.find(x)]);
       
   607       for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {   
       
   608 	if ( position[tit] == D ) {
       
   609 	  int b = blossom.find(tit);
       
   610 	  for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) { 
       
   611 	    position.set(bit, C);
       
   612 	  }
       
   613 	  blossom.eraseClass(b);
       
   614 	} else position.set(tit, C);
       
   615       }
       
   616       tree.eraseClass(y);
       
   617 
       
   618     }
   347 
   619 
   348   };
   620   };
   349 
       
   350 
       
   351   // **********************************************************************
       
   352   //  IMPLEMENTATIONS
       
   353   // **********************************************************************
       
   354 
       
   355 
       
   356   template <typename Graph>
       
   357   void MaxMatching<Graph>::lateShrink(Node v, typename Graph::template 
       
   358 				      NodeMap<Node>& ear, UFE& blossom, 
       
   359                                       UFE& tree) {
       
   360     //We have one tree which we grow, and also shrink but only if it cannot be
       
   361     //postponed. If we augment then we return to the "for" cycle of
       
   362     //runEdmonds().
       
   363 
       
   364     std::queue<Node> Q;   //queue of the totally unscanned nodes
       
   365     Q.push(v);  
       
   366     std::queue<Node> R;   
       
   367     //queue of the nodes which must be scanned for a possible shrink
       
   368       
       
   369     while ( !Q.empty() ) {
       
   370       Node x=Q.front();
       
   371       Q.pop();
       
   372       for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
       
   373 	Node y=g.runningNode(e);
       
   374 	//growOrAugment grows if y is covered by the matching and
       
   375 	//augments if not. In this latter case it returns 1.
       
   376 	if ( position[y]==C && growOrAugment(y, x, ear, blossom, tree, Q) ) 
       
   377           return;
       
   378       }
       
   379       R.push(x);
       
   380     }
       
   381       
       
   382     while ( !R.empty() ) {
       
   383       Node x=R.front();
       
   384       R.pop();
       
   385 	
       
   386       for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
       
   387 	Node y=g.runningNode(e);
       
   388 
       
   389 	if ( position[y] == D && blossom.find(x) != blossom.find(y) )  
       
   390 	  //Recall that we have only one tree.
       
   391 	  shrink( x, y, ear, blossom, tree, Q);	
       
   392 	
       
   393 	while ( !Q.empty() ) {
       
   394 	  Node z=Q.front();
       
   395 	  Q.pop();
       
   396 	  for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
       
   397 	    Node w=g.runningNode(f);
       
   398 	    //growOrAugment grows if y is covered by the matching and
       
   399 	    //augments if not. In this latter case it returns 1.
       
   400 	    if ( position[w]==C && growOrAugment(w, z, ear, blossom, tree, Q) )
       
   401               return;
       
   402 	  }
       
   403 	  R.push(z);
       
   404 	}
       
   405       } //for e
       
   406     } // while ( !R.empty() )
       
   407   }
       
   408 
       
   409 
       
   410   template <typename Graph>
       
   411   void MaxMatching<Graph>::normShrink(Node v,
       
   412 				      typename Graph::template
       
   413 				      NodeMap<Node>& ear,  
       
   414 				      UFE& blossom, UFE& tree) {
       
   415     //We have one tree, which we grow and shrink. If we augment then we
       
   416     //return to the "for" cycle of runEdmonds().
       
   417     
       
   418     std::queue<Node> Q;   //queue of the unscanned nodes
       
   419     Q.push(v);  
       
   420     while ( !Q.empty() ) {
       
   421 
       
   422       Node x=Q.front();
       
   423       Q.pop();
       
   424 	
       
   425       for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
       
   426 	Node y=g.runningNode(e);
       
   427 	      
       
   428 	switch ( position[y] ) {
       
   429 	case D:          //x and y must be in the same tree
       
   430 	  if ( blossom.find(x) != blossom.find(y) )
       
   431 	    //x and y are in the same tree
       
   432 	    shrink( x, y, ear, blossom, tree, Q);
       
   433 	  break;
       
   434 	case C:
       
   435 	  //growOrAugment grows if y is covered by the matching and
       
   436 	  //augments if not. In this latter case it returns 1.
       
   437 	  if ( growOrAugment(y, x, ear, blossom, tree, Q) ) return;
       
   438 	  break;
       
   439 	default: break;
       
   440        	}
       
   441       }
       
   442     }
       
   443   }
       
   444   
       
   445 
       
   446   template <typename Graph>
       
   447     void MaxMatching<Graph>::shrink(Node x,Node y, typename 
       
   448 				    Graph::template NodeMap<Node>& ear,  
       
   449 				    UFE& blossom, UFE& tree, std::queue<Node>& Q) {
       
   450     //x and y are the two adjacent vertices in two blossoms.
       
   451    
       
   452     typename Graph::template NodeMap<bool> path(g,false);
       
   453     
       
   454     Node b=blossom.find(x);
       
   455     path.set(b,true);
       
   456     b=_mate[b];
       
   457     while ( b!=INVALID ) { 
       
   458       b=blossom.find(ear[b]);
       
   459       path.set(b,true);
       
   460       b=_mate[b];
       
   461     } //we go until the root through bases of blossoms and odd vertices
       
   462     
       
   463     Node top=y;
       
   464     Node middle=blossom.find(top);
       
   465     Node bottom=x;
       
   466     while ( !path[middle] )
       
   467       shrinkStep(top, middle, bottom, ear, blossom, tree, Q);
       
   468     //Until we arrive to a node on the path, we update blossom, tree
       
   469     //and the positions of the odd nodes.
       
   470     
       
   471     Node base=middle;
       
   472     top=x;
       
   473     middle=blossom.find(top);
       
   474     bottom=y;
       
   475     Node blossom_base=blossom.find(base);
       
   476     while ( middle!=blossom_base )
       
   477       shrinkStep(top, middle, bottom, ear, blossom, tree, Q);
       
   478     //Until we arrive to a node on the path, we update blossom, tree
       
   479     //and the positions of the odd nodes.
       
   480     
       
   481     blossom.makeRep(base);
       
   482   }
       
   483 
       
   484 
       
   485 
       
   486   template <typename Graph>
       
   487   void MaxMatching<Graph>::shrinkStep(Node& top, Node& middle, Node& bottom,
       
   488 				      typename Graph::template
       
   489 				      NodeMap<Node>& ear,  
       
   490 				      UFE& blossom, UFE& tree,
       
   491 				      std::queue<Node>& Q) {
       
   492     //We traverse a blossom and update everything.
       
   493     
       
   494     ear.set(top,bottom);
       
   495     Node t=top;
       
   496     while ( t!=middle ) {
       
   497       Node u=_mate[t];
       
   498       t=ear[u];
       
   499       ear.set(t,u);
       
   500     } 
       
   501     bottom=_mate[middle];
       
   502     position.set(bottom,D);
       
   503     Q.push(bottom);
       
   504     top=ear[bottom];		
       
   505     Node oldmiddle=middle;
       
   506     middle=blossom.find(top);
       
   507     tree.erase(bottom);
       
   508     tree.erase(oldmiddle);
       
   509     blossom.insert(bottom);
       
   510     blossom.join(bottom, oldmiddle);
       
   511     blossom.join(top, oldmiddle);
       
   512   }
       
   513 
       
   514 
       
   515   template <typename Graph>
       
   516   bool MaxMatching<Graph>::growOrAugment(Node& y, Node& x, typename Graph::template
       
   517 					 NodeMap<Node>& ear, UFE& blossom, UFE& tree,
       
   518 					 std::queue<Node>& Q) {
       
   519     //x is in a blossom in the tree, y is outside. If y is covered by
       
   520     //the matching we grow, otherwise we augment. In this case we
       
   521     //return 1.
       
   522     
       
   523     if ( _mate[y]!=INVALID ) {       //grow
       
   524       ear.set(y,x);
       
   525       Node w=_mate[y];
       
   526       blossom.insert(w);
       
   527       position.set(y,A);
       
   528       position.set(w,D);
       
   529       tree.insert(y);
       
   530       tree.insert(w);
       
   531       tree.join(y,blossom.find(x));  
       
   532       tree.join(w,y);  
       
   533       Q.push(w);
       
   534     } else {                      //augment 
       
   535       augment(x, ear, blossom, tree);
       
   536       _mate.set(x,y);
       
   537       _mate.set(y,x);
       
   538       return true;
       
   539     }
       
   540     return false;
       
   541   }
       
   542   
       
   543 
       
   544   template <typename Graph>
       
   545   void MaxMatching<Graph>::augment(Node x,
       
   546 				   typename Graph::template NodeMap<Node>& ear,  
       
   547 				   UFE& blossom, UFE& tree) { 
       
   548     Node v=_mate[x];
       
   549     while ( v!=INVALID ) {
       
   550 	
       
   551       Node u=ear[v];
       
   552       _mate.set(v,u);
       
   553       Node tmp=v;
       
   554       v=_mate[u];
       
   555       _mate.set(u,tmp);
       
   556     }
       
   557     Node y=blossom.find(x);
       
   558     for (typename UFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {   
       
   559       if ( position[tit] == D ) {
       
   560 	for (typename UFE::ItemIt bit(blossom, tit); bit != INVALID; ++bit) {  
       
   561 	  position.set( bit ,C);
       
   562 	}
       
   563 	blossom.eraseClass(tit);
       
   564       } else position.set( tit ,C);
       
   565     }
       
   566     tree.eraseClass(y);
       
   567 
       
   568   }
       
   569 
       
   570   
   621   
   571 } //END OF NAMESPACE LEMON
   622 } //END OF NAMESPACE LEMON
   572 
   623 
   573 #endif //LEMON_MAX_MATCHING_H
   624 #endif //LEMON_MAX_MATCHING_H