1.1 --- a/lemon/max_matching.h Tue Oct 30 10:51:07 2007 +0000
1.2 +++ b/lemon/max_matching.h Tue Oct 30 20:21:10 2007 +0000
1.3 @@ -30,25 +30,21 @@
1.4
1.5 namespace lemon {
1.6
1.7 - /// \ingroup matching
1.8 + ///\ingroup matching
1.9
1.10 - ///Edmonds' alternating forest maximum matching algorithm.
1.11 -
1.12 + ///\brief Edmonds' alternating forest maximum matching algorithm.
1.13 + ///
1.14 ///This class provides Edmonds' alternating forest matching
1.15 ///algorithm. The starting matching (if any) can be passed to the
1.16 - ///algorithm using read-in functions \ref readNMapNode, \ref
1.17 - ///readNMapEdge or \ref readEMapBool depending on the container. The
1.18 - ///resulting maximum matching can be attained by write-out functions
1.19 - ///\ref writeNMapNode, \ref writeNMapEdge or \ref writeEMapBool
1.20 - ///depending on the preferred container.
1.21 + ///algorithm using some of init functions.
1.22 ///
1.23 ///The dual side of a matching is a map of the nodes to
1.24 - ///MaxMatching::pos_enum, having values D, A and C showing the
1.25 - ///Gallai-Edmonds decomposition of the graph. The nodes in D induce
1.26 - ///a graph with factor-critical components, the nodes in A form the
1.27 - ///barrier, and the nodes in C induce a graph having a perfect
1.28 - ///matching. This decomposition can be attained by calling \ref
1.29 - ///writePos after running the algorithm.
1.30 + ///MaxMatching::DecompType, having values \c D, \c A and \c C
1.31 + ///showing the Gallai-Edmonds decomposition of the graph. The nodes
1.32 + ///in \c D induce a graph with factor-critical components, the nodes
1.33 + ///in \c A form the barrier, and the nodes in \c C induce a graph
1.34 + ///having a perfect matching. This decomposition can be attained by
1.35 + ///calling \c decomposition() after running the algorithm.
1.36 ///
1.37 ///\param Graph The undirected graph type the algorithm runs on.
1.38 ///
1.39 @@ -67,17 +63,21 @@
1.40
1.41 typedef typename Graph::template NodeMap<int> UFECrossRef;
1.42 typedef UnionFindEnum<UFECrossRef> UFE;
1.43 + typedef std::vector<Node> NV;
1.44 +
1.45 + typedef typename Graph::template NodeMap<int> EFECrossRef;
1.46 + typedef ExtendFindEnum<EFECrossRef> EFE;
1.47
1.48 public:
1.49
1.50 - ///Indicates the Gallai-Edmonds decomposition of the graph.
1.51 -
1.52 + ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
1.53 + ///
1.54 ///Indicates the Gallai-Edmonds decomposition of the graph, which
1.55 ///shows an upper bound on the size of a maximum matching. The
1.56 - ///nodes with pos_enum \c D induce a graph with factor-critical
1.57 + ///nodes with DecompType \c D induce a graph with factor-critical
1.58 ///components, the nodes in \c A form the canonical barrier, and the
1.59 ///nodes in \c C induce a graph having a perfect matching.
1.60 - enum pos_enum {
1.61 + enum DecompType {
1.62 D=0,
1.63 A=1,
1.64 C=2
1.65 @@ -88,71 +88,32 @@
1.66 static const int HEUR_density=2;
1.67 const Graph& g;
1.68 typename Graph::template NodeMap<Node> _mate;
1.69 - typename Graph::template NodeMap<pos_enum> position;
1.70 + typename Graph::template NodeMap<DecompType> position;
1.71
1.72 public:
1.73
1.74 - MaxMatching(const Graph& _g) : g(_g), _mate(_g,INVALID), position(_g) {}
1.75 + MaxMatching(const Graph& _g)
1.76 + : g(_g), _mate(_g), position(_g) {}
1.77
1.78 - ///Runs Edmonds' algorithm.
1.79 -
1.80 - ///Runs Edmonds' algorithm for sparse graphs (number of edges <
1.81 - ///2*number of nodes), and a heuristical Edmonds' algorithm with a
1.82 - ///heuristic of postponing shrinks for dense graphs.
1.83 - void run() {
1.84 - if ( countUEdges(g) < HEUR_density*countNodes(g) ) {
1.85 - greedyMatching();
1.86 - runEdmonds(0);
1.87 - } else runEdmonds(1);
1.88 - }
1.89 -
1.90 -
1.91 - ///Runs Edmonds' algorithm.
1.92 -
1.93 - ///If heur=0 it runs Edmonds' algorithm. If heur=1 it runs
1.94 - ///Edmonds' algorithm with a heuristic of postponing shrinks,
1.95 - ///giving a faster algorithm for dense graphs.
1.96 - void runEdmonds( int heur = 1 ) {
1.97 -
1.98 - //each vertex is put to C
1.99 - for(NodeIt v(g); v!=INVALID; ++v)
1.100 - position.set(v,C);
1.101 -
1.102 - typename Graph::template NodeMap<Node> ear(g,INVALID);
1.103 - //undefined for the base nodes of the blossoms (i.e. for the
1.104 - //representative elements of UFE blossom) and for the nodes in C
1.105 -
1.106 - UFECrossRef blossom_base(g);
1.107 - UFE blossom(blossom_base);
1.108 -
1.109 - UFECrossRef tree_base(g);
1.110 - UFE tree(tree_base);
1.111 -
1.112 - //If these UFE's would be members of the class then also
1.113 - //blossom_base and tree_base should be a member.
1.114 -
1.115 - //We build only one tree and the other vertices uncovered by the
1.116 - //matching belong to C. (They can be considered as singleton
1.117 - //trees.) If this tree can be augmented or no more
1.118 - //grow/augmentation/shrink is possible then we return to this
1.119 - //"for" cycle.
1.120 + ///\brief Sets the actual matching to the empty matching.
1.121 + ///
1.122 + ///Sets the actual matching to the empty matching.
1.123 + ///
1.124 + void init() {
1.125 for(NodeIt v(g); v!=INVALID; ++v) {
1.126 - if ( position[v]==C && _mate[v]==INVALID ) {
1.127 - blossom.insert(v);
1.128 - tree.insert(v);
1.129 - position.set(v,D);
1.130 - if ( heur == 1 ) lateShrink( v, ear, blossom, tree );
1.131 - else normShrink( v, ear, blossom, tree );
1.132 - }
1.133 + _mate.set(v,INVALID);
1.134 + position.set(v,C);
1.135 }
1.136 }
1.137
1.138 -
1.139 - ///Finds a greedy matching starting from the actual matching.
1.140 -
1.141 - ///Starting form the actual matching stored, it finds a maximal
1.142 - ///greedy matching.
1.143 - void greedyMatching() {
1.144 + ///\brief Finds a greedy matching for initial matching.
1.145 + ///
1.146 + ///For initial matchig it finds a maximal greedy matching.
1.147 + void greedyInit() {
1.148 + for(NodeIt v(g); v!=INVALID; ++v) {
1.149 + _mate.set(v,INVALID);
1.150 + position.set(v,C);
1.151 + }
1.152 for(NodeIt v(g); v!=INVALID; ++v)
1.153 if ( _mate[v]==INVALID ) {
1.154 for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
1.155 @@ -166,8 +127,153 @@
1.156 }
1.157 }
1.158
1.159 - ///Returns the size of the actual matching stored.
1.160 + ///\brief Initialize the matching from each nodes' mate.
1.161 + ///
1.162 + ///Initialize the matching from a \c Node valued \c Node map. This
1.163 + ///map must be \e symmetric, i.e. if \c map[u]==v then \c
1.164 + ///map[v]==u must hold, and \c uv will be an edge of the initial
1.165 + ///matching.
1.166 + template <typename MateMap>
1.167 + void mateMapInit(MateMap& map) {
1.168 + for(NodeIt v(g); v!=INVALID; ++v) {
1.169 + _mate.set(v,map[v]);
1.170 + position.set(v,C);
1.171 + }
1.172 + }
1.173
1.174 + ///\brief Initialize the matching from a node map with the
1.175 + ///incident matching edges.
1.176 + ///
1.177 + ///Initialize the matching from an \c UEdge valued \c Node map. \c
1.178 + ///map[v] must be an \c UEdge incident to \c v. This map must have
1.179 + ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
1.180 + ///g.oppositeNode(v,map[v])==u holds, and now some edge joining \c
1.181 + ///u to \c v will be an edge of the matching.
1.182 + template<typename MatchingMap>
1.183 + void matchingMapInit(MatchingMap& map) {
1.184 + for(NodeIt v(g); v!=INVALID; ++v) {
1.185 + position.set(v,C);
1.186 + UEdge e=map[v];
1.187 + if ( e!=INVALID )
1.188 + _mate.set(v,g.oppositeNode(v,e));
1.189 + else
1.190 + _mate.set(v,INVALID);
1.191 + }
1.192 + }
1.193 +
1.194 + ///\brief Initialize the matching from the map containing the
1.195 + ///undirected matching edges.
1.196 + ///
1.197 + ///Initialize the matching from a \c bool valued \c UEdge map. This
1.198 + ///map must have the property that there are no two incident edges
1.199 + ///\c e, \c f with \c map[e]==map[f]==true. The edges \c e with \c
1.200 + ///map[e]==true form the matching.
1.201 + template <typename MatchingMap>
1.202 + void matchingInit(MatchingMap& map) {
1.203 + for(NodeIt v(g); v!=INVALID; ++v) {
1.204 + _mate.set(v,INVALID);
1.205 + position.set(v,C);
1.206 + }
1.207 + for(UEdgeIt e(g); e!=INVALID; ++e) {
1.208 + if ( map[e] ) {
1.209 + Node u=g.source(e);
1.210 + Node v=g.target(e);
1.211 + _mate.set(u,v);
1.212 + _mate.set(v,u);
1.213 + }
1.214 + }
1.215 + }
1.216 +
1.217 +
1.218 + ///\brief Runs Edmonds' algorithm.
1.219 + ///
1.220 + ///Runs Edmonds' algorithm for sparse graphs (number of edges <
1.221 + ///2*number of nodes), and a heuristical Edmonds' algorithm with a
1.222 + ///heuristic of postponing shrinks for dense graphs.
1.223 + void run() {
1.224 + if (countUEdges(g) < HEUR_density * countNodes(g)) {
1.225 + greedyInit();
1.226 + startSparse();
1.227 + } else {
1.228 + init();
1.229 + startDense();
1.230 + }
1.231 + }
1.232 +
1.233 +
1.234 + ///\brief Starts Edmonds' algorithm.
1.235 + ///
1.236 + ///If runs the original Edmonds' algorithm.
1.237 + void startSparse() {
1.238 +
1.239 + typename Graph::template NodeMap<Node> ear(g,INVALID);
1.240 + //undefined for the base nodes of the blossoms (i.e. for the
1.241 + //representative elements of UFE blossom) and for the nodes in C
1.242 +
1.243 + UFECrossRef blossom_base(g);
1.244 + UFE blossom(blossom_base);
1.245 + NV rep(countNodes(g));
1.246 +
1.247 + EFECrossRef tree_base(g);
1.248 + EFE tree(tree_base);
1.249 +
1.250 + //If these UFE's would be members of the class then also
1.251 + //blossom_base and tree_base should be a member.
1.252 +
1.253 + //We build only one tree and the other vertices uncovered by the
1.254 + //matching belong to C. (They can be considered as singleton
1.255 + //trees.) If this tree can be augmented or no more
1.256 + //grow/augmentation/shrink is possible then we return to this
1.257 + //"for" cycle.
1.258 + for(NodeIt v(g); v!=INVALID; ++v) {
1.259 + if (position[v]==C && _mate[v]==INVALID) {
1.260 + rep[blossom.insert(v)] = v;
1.261 + tree.insert(v);
1.262 + position.set(v,D);
1.263 + normShrink(v, ear, blossom, rep, tree);
1.264 + }
1.265 + }
1.266 + }
1.267 +
1.268 + ///\brief Starts Edmonds' algorithm.
1.269 + ///
1.270 + ///It runs Edmonds' algorithm with a heuristic of postponing
1.271 + ///shrinks, giving a faster algorithm for dense graphs.
1.272 + void startDense() {
1.273 +
1.274 + typename Graph::template NodeMap<Node> ear(g,INVALID);
1.275 + //undefined for the base nodes of the blossoms (i.e. for the
1.276 + //representative elements of UFE blossom) and for the nodes in C
1.277 +
1.278 + UFECrossRef blossom_base(g);
1.279 + UFE blossom(blossom_base);
1.280 + NV rep(countNodes(g));
1.281 +
1.282 + EFECrossRef tree_base(g);
1.283 + EFE tree(tree_base);
1.284 +
1.285 + //If these UFE's would be members of the class then also
1.286 + //blossom_base and tree_base should be a member.
1.287 +
1.288 + //We build only one tree and the other vertices uncovered by the
1.289 + //matching belong to C. (They can be considered as singleton
1.290 + //trees.) If this tree can be augmented or no more
1.291 + //grow/augmentation/shrink is possible then we return to this
1.292 + //"for" cycle.
1.293 + for(NodeIt v(g); v!=INVALID; ++v) {
1.294 + if ( position[v]==C && _mate[v]==INVALID ) {
1.295 + rep[blossom.insert(v)] = v;
1.296 + tree.insert(v);
1.297 + position.set(v,D);
1.298 + lateShrink(v, ear, blossom, rep, tree);
1.299 + }
1.300 + }
1.301 + }
1.302 +
1.303 +
1.304 +
1.305 + ///\brief Returns the size of the actual matching stored.
1.306 + ///
1.307 ///Returns the size of the actual matching stored. After \ref
1.308 ///run() it returns the size of a maximum matching in the graph.
1.309 int size() const {
1.310 @@ -181,75 +287,70 @@
1.311 }
1.312
1.313
1.314 - ///Resets the actual matching to the empty matching.
1.315 -
1.316 - ///Resets the actual matching to the empty matching.
1.317 + ///\brief Returns the mate of a node in the actual matching.
1.318 ///
1.319 - void resetMatching() {
1.320 - for(NodeIt v(g); v!=INVALID; ++v)
1.321 - _mate.set(v,INVALID);
1.322 - }
1.323 -
1.324 - ///Returns the mate of a node in the actual matching.
1.325 -
1.326 ///Returns the mate of a \c node in the actual matching.
1.327 ///Returns INVALID if the \c node is not covered by the actual matching.
1.328 - Node mate(Node& node) const {
1.329 + Node mate(const Node& node) const {
1.330 return _mate[node];
1.331 }
1.332
1.333 - ///Reads a matching from a \c Node valued \c Node map.
1.334 + ///\brief Returns the matching edge incident to the given node.
1.335 + ///
1.336 + ///Returns the matching edge of a \c node in the actual matching.
1.337 + ///Returns INVALID if the \c node is not covered by the actual matching.
1.338 + UEdge matchingEdge(const Node& node) const {
1.339 + if (_mate[node] == INVALID) return INVALID;
1.340 + Node n = node < _mate[node] ? node : _mate[node];
1.341 + for (IncEdgeIt e(g, n); e != INVALID; ++e) {
1.342 + if (g.oppositeNode(n, e) == _mate[n]) {
1.343 + return e;
1.344 + }
1.345 + }
1.346 + return INVALID;
1.347 + }
1.348
1.349 - ///Reads a matching from a \c Node valued \c Node map. This map
1.350 - ///must be \e symmetric, i.e. if \c map[u]==v then \c map[v]==u
1.351 - ///must hold, and \c uv will be an edge of the matching.
1.352 - template<typename NMapN>
1.353 - void readNMapNode(NMapN& map) {
1.354 - for(NodeIt v(g); v!=INVALID; ++v) {
1.355 - _mate.set(v,map[v]);
1.356 - }
1.357 - }
1.358 + /// \brief Returns the class of the node in the Edmonds-Gallai
1.359 + /// decomposition.
1.360 + ///
1.361 + /// Returns the class of the node in the Edmonds-Gallai
1.362 + /// decomposition.
1.363 + DecompType decomposition(const Node& n) {
1.364 + return position[n] == A;
1.365 + }
1.366 +
1.367 + /// \brief Returns true when the node is in the barrier.
1.368 + ///
1.369 + /// Returns true when the node is in the barrier.
1.370 + bool barrier(const Node& n) {
1.371 + return position[n] == A;
1.372 + }
1.373
1.374 - ///Writes the stored matching to a \c Node valued \c Node map.
1.375 -
1.376 + ///\brief Gives back the matching in a \c Node of mates.
1.377 + ///
1.378 ///Writes the stored matching to a \c Node valued \c Node map. The
1.379 ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
1.380 ///map[v]==u will hold, and now \c uv is an edge of the matching.
1.381 - template<typename NMapN>
1.382 - void writeNMapNode (NMapN& map) const {
1.383 + template <typename MateMap>
1.384 + void mateMap(MateMap& map) const {
1.385 for(NodeIt v(g); v!=INVALID; ++v) {
1.386 map.set(v,_mate[v]);
1.387 }
1.388 }
1.389 -
1.390 - ///Reads a matching from an \c UEdge valued \c Node map.
1.391 -
1.392 - ///Reads a matching from an \c UEdge valued \c Node map. \c
1.393 - ///map[v] must be an \c UEdge incident to \c v. This map must
1.394 - ///have the property that if \c g.oppositeNode(u,map[u])==v then
1.395 - ///\c \c g.oppositeNode(v,map[v])==u holds, and now some edge
1.396 - ///joining \c u to \c v will be an edge of the matching.
1.397 - template<typename NMapE>
1.398 - void readNMapEdge(NMapE& map) {
1.399 - for(NodeIt v(g); v!=INVALID; ++v) {
1.400 - UEdge e=map[v];
1.401 - if ( e!=INVALID )
1.402 - _mate.set(v,g.oppositeNode(v,e));
1.403 - }
1.404 - }
1.405
1.406 - ///Writes the matching stored to an \c UEdge valued \c Node map.
1.407 -
1.408 + ///\brief Gives back the matching in an \c UEdge valued \c Node
1.409 + ///map.
1.410 + ///
1.411 ///Writes the stored matching to an \c UEdge valued \c Node
1.412 ///map. \c map[v] will be an \c UEdge incident to \c v. This
1.413 ///map will have the property that if \c g.oppositeNode(u,map[u])
1.414 ///== v then \c map[u]==map[v] holds, and now this edge is an edge
1.415 ///of the matching.
1.416 - template<typename NMapE>
1.417 - void writeNMapEdge (NMapE& map) const {
1.418 + template <typename MatchingMap>
1.419 + void matchingMap(MatchingMap& map) const {
1.420 typename Graph::template NodeMap<bool> todo(g,true);
1.421 for(NodeIt v(g); v!=INVALID; ++v) {
1.422 - if ( todo[v] && _mate[v]!=INVALID ) {
1.423 + if (_mate[v]!=INVALID && v < _mate[v]) {
1.424 Node u=_mate[v];
1.425 for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
1.426 if ( g.runningNode(e) == u ) {
1.427 @@ -265,33 +366,15 @@
1.428 }
1.429
1.430
1.431 - ///Reads a matching from a \c bool valued \c Edge map.
1.432 -
1.433 - ///Reads a matching from a \c bool valued \c Edge map. This map
1.434 - ///must have the property that there are no two incident edges \c
1.435 - ///e, \c f with \c map[e]==map[f]==true. The edges \c e with \c
1.436 - ///map[e]==true form the matching.
1.437 - template<typename EMapB>
1.438 - void readEMapBool(EMapB& map) {
1.439 - for(UEdgeIt e(g); e!=INVALID; ++e) {
1.440 - if ( map[e] ) {
1.441 - Node u=g.source(e);
1.442 - Node v=g.target(e);
1.443 - _mate.set(u,v);
1.444 - _mate.set(v,u);
1.445 - }
1.446 - }
1.447 - }
1.448 -
1.449 -
1.450 - ///Writes the matching stored to a \c bool valued \c Edge map.
1.451 -
1.452 + ///\brief Gives back the matching in a \c bool valued \c UEdge
1.453 + ///map.
1.454 + ///
1.455 ///Writes the matching stored to a \c bool valued \c Edge
1.456 ///map. This map will have the property that there are no two
1.457 ///incident edges \c e, \c f with \c map[e]==map[f]==true. The
1.458 ///edges \c e with \c map[e]==true form the matching.
1.459 - template<typename EMapB>
1.460 - void writeEMapBool (EMapB& map) const {
1.461 + template<typename MatchingMap>
1.462 + void matching(MatchingMap& map) const {
1.463 for(UEdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
1.464
1.465 typename Graph::template NodeMap<bool> todo(g,true);
1.466 @@ -311,262 +394,230 @@
1.467 }
1.468
1.469
1.470 - ///Writes the canonical decomposition of the graph after running
1.471 + ///\brief Returns the canonical decomposition of the graph after running
1.472 ///the algorithm.
1.473 -
1.474 + ///
1.475 ///After calling any run methods of the class, it writes the
1.476 ///Gallai-Edmonds canonical decomposition of the graph. \c map
1.477 - ///must be a node map of \ref pos_enum 's.
1.478 - template<typename NMapEnum>
1.479 - void writePos (NMapEnum& map) const {
1.480 - for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
1.481 + ///must be a node map of \ref DecompType 's.
1.482 + template <typename DecompositionMap>
1.483 + void decomposition(DecompositionMap& map) const {
1.484 + for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
1.485 + }
1.486 +
1.487 + ///\brief Returns a barrier on the nodes.
1.488 + ///
1.489 + ///After calling any run methods of the class, it writes a
1.490 + ///canonical barrier on the nodes. The odd component number of the
1.491 + ///remaining graph minus the barrier size is a lower bound for the
1.492 + ///uncovered nodes in the graph. The \c map must be a node map of
1.493 + ///bools.
1.494 + template <typename BarrierMap>
1.495 + void barrier(BarrierMap& barrier) {
1.496 + for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);
1.497 }
1.498
1.499 private:
1.500
1.501
1.502 void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,
1.503 - UFE& blossom, UFE& tree);
1.504 + UFE& blossom, NV& rep, EFE& tree) {
1.505 + //We have one tree which we grow, and also shrink but only if it
1.506 + //cannot be postponed. If we augment then we return to the "for"
1.507 + //cycle of runEdmonds().
1.508 +
1.509 + std::queue<Node> Q; //queue of the totally unscanned nodes
1.510 + Q.push(v);
1.511 + std::queue<Node> R;
1.512 + //queue of the nodes which must be scanned for a possible shrink
1.513 +
1.514 + while ( !Q.empty() ) {
1.515 + Node x=Q.front();
1.516 + Q.pop();
1.517 + for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
1.518 + Node y=g.runningNode(e);
1.519 + //growOrAugment grows if y is covered by the matching and
1.520 + //augments if not. In this latter case it returns 1.
1.521 + if (position[y]==C &&
1.522 + growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
1.523 + }
1.524 + R.push(x);
1.525 + }
1.526 +
1.527 + while ( !R.empty() ) {
1.528 + Node x=R.front();
1.529 + R.pop();
1.530 +
1.531 + for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
1.532 + Node y=g.runningNode(e);
1.533 +
1.534 + if ( position[y] == D && blossom.find(x) != blossom.find(y) )
1.535 + //Recall that we have only one tree.
1.536 + shrink( x, y, ear, blossom, rep, tree, Q);
1.537 +
1.538 + while ( !Q.empty() ) {
1.539 + Node z=Q.front();
1.540 + Q.pop();
1.541 + for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
1.542 + Node w=g.runningNode(f);
1.543 + //growOrAugment grows if y is covered by the matching and
1.544 + //augments if not. In this latter case it returns 1.
1.545 + if (position[w]==C &&
1.546 + growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
1.547 + }
1.548 + R.push(z);
1.549 + }
1.550 + } //for e
1.551 + } // while ( !R.empty() )
1.552 + }
1.553
1.554 void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,
1.555 - UFE& blossom, UFE& tree);
1.556 + UFE& blossom, NV& rep, EFE& tree) {
1.557 + //We have one tree, which we grow and shrink. If we augment then we
1.558 + //return to the "for" cycle of runEdmonds().
1.559 +
1.560 + std::queue<Node> Q; //queue of the unscanned nodes
1.561 + Q.push(v);
1.562 + while ( !Q.empty() ) {
1.563 +
1.564 + Node x=Q.front();
1.565 + Q.pop();
1.566 +
1.567 + for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
1.568 + Node y=g.runningNode(e);
1.569 +
1.570 + switch ( position[y] ) {
1.571 + case D: //x and y must be in the same tree
1.572 + if ( blossom.find(x) != blossom.find(y))
1.573 + //x and y are in the same tree
1.574 + shrink(x, y, ear, blossom, rep, tree, Q);
1.575 + break;
1.576 + case C:
1.577 + //growOrAugment grows if y is covered by the matching and
1.578 + //augments if not. In this latter case it returns 1.
1.579 + if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
1.580 + break;
1.581 + default: break;
1.582 + }
1.583 + }
1.584 + }
1.585 + }
1.586
1.587 void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,
1.588 - UFE& blossom, UFE& tree,std::queue<Node>& Q);
1.589 + UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
1.590 + //x and y are the two adjacent vertices in two blossoms.
1.591 +
1.592 + typename Graph::template NodeMap<bool> path(g,false);
1.593 +
1.594 + Node b=rep[blossom.find(x)];
1.595 + path.set(b,true);
1.596 + b=_mate[b];
1.597 + while ( b!=INVALID ) {
1.598 + b=rep[blossom.find(ear[b])];
1.599 + path.set(b,true);
1.600 + b=_mate[b];
1.601 + } //we go until the root through bases of blossoms and odd vertices
1.602 +
1.603 + Node top=y;
1.604 + Node middle=rep[blossom.find(top)];
1.605 + Node bottom=x;
1.606 + while ( !path[middle] )
1.607 + shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
1.608 + //Until we arrive to a node on the path, we update blossom, tree
1.609 + //and the positions of the odd nodes.
1.610 +
1.611 + Node base=middle;
1.612 + top=x;
1.613 + middle=rep[blossom.find(top)];
1.614 + bottom=y;
1.615 + Node blossom_base=rep[blossom.find(base)];
1.616 + while ( middle!=blossom_base )
1.617 + shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
1.618 + //Until we arrive to a node on the path, we update blossom, tree
1.619 + //and the positions of the odd nodes.
1.620 +
1.621 + rep[blossom.find(base)] = base;
1.622 + }
1.623
1.624 void shrinkStep(Node& top, Node& middle, Node& bottom,
1.625 typename Graph::template NodeMap<Node>& ear,
1.626 - UFE& blossom, UFE& tree, std::queue<Node>& Q);
1.627 + UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
1.628 + //We traverse a blossom and update everything.
1.629 +
1.630 + ear.set(top,bottom);
1.631 + Node t=top;
1.632 + while ( t!=middle ) {
1.633 + Node u=_mate[t];
1.634 + t=ear[u];
1.635 + ear.set(t,u);
1.636 + }
1.637 + bottom=_mate[middle];
1.638 + position.set(bottom,D);
1.639 + Q.push(bottom);
1.640 + top=ear[bottom];
1.641 + Node oldmiddle=middle;
1.642 + middle=rep[blossom.find(top)];
1.643 + tree.erase(bottom);
1.644 + tree.erase(oldmiddle);
1.645 + blossom.insert(bottom);
1.646 + blossom.join(bottom, oldmiddle);
1.647 + blossom.join(top, oldmiddle);
1.648 + }
1.649 +
1.650 +
1.651
1.652 bool growOrAugment(Node& y, Node& x, typename Graph::template
1.653 - NodeMap<Node>& ear, UFE& blossom, UFE& tree,
1.654 - std::queue<Node>& Q);
1.655 + NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree,
1.656 + std::queue<Node>& Q) {
1.657 + //x is in a blossom in the tree, y is outside. If y is covered by
1.658 + //the matching we grow, otherwise we augment. In this case we
1.659 + //return 1.
1.660 +
1.661 + if ( _mate[y]!=INVALID ) { //grow
1.662 + ear.set(y,x);
1.663 + Node w=_mate[y];
1.664 + rep[blossom.insert(w)] = w;
1.665 + position.set(y,A);
1.666 + position.set(w,D);
1.667 + int t = tree.find(rep[blossom.find(x)]);
1.668 + tree.insert(y,t);
1.669 + tree.insert(w,t);
1.670 + Q.push(w);
1.671 + } else { //augment
1.672 + augment(x, ear, blossom, rep, tree);
1.673 + _mate.set(x,y);
1.674 + _mate.set(y,x);
1.675 + return true;
1.676 + }
1.677 + return false;
1.678 + }
1.679
1.680 void augment(Node x, typename Graph::template NodeMap<Node>& ear,
1.681 - UFE& blossom, UFE& tree);
1.682 + UFE& blossom, NV& rep, EFE& tree) {
1.683 + Node v=_mate[x];
1.684 + while ( v!=INVALID ) {
1.685 +
1.686 + Node u=ear[v];
1.687 + _mate.set(v,u);
1.688 + Node tmp=v;
1.689 + v=_mate[u];
1.690 + _mate.set(u,tmp);
1.691 + }
1.692 + int y = tree.find(rep[blossom.find(x)]);
1.693 + for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {
1.694 + if ( position[tit] == D ) {
1.695 + int b = blossom.find(tit);
1.696 + for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) {
1.697 + position.set(bit, C);
1.698 + }
1.699 + blossom.eraseClass(b);
1.700 + } else position.set(tit, C);
1.701 + }
1.702 + tree.eraseClass(y);
1.703 +
1.704 + }
1.705
1.706 };
1.707 -
1.708 -
1.709 - // **********************************************************************
1.710 - // IMPLEMENTATIONS
1.711 - // **********************************************************************
1.712 -
1.713 -
1.714 - template <typename Graph>
1.715 - void MaxMatching<Graph>::lateShrink(Node v, typename Graph::template
1.716 - NodeMap<Node>& ear, UFE& blossom,
1.717 - UFE& tree) {
1.718 - //We have one tree which we grow, and also shrink but only if it cannot be
1.719 - //postponed. If we augment then we return to the "for" cycle of
1.720 - //runEdmonds().
1.721 -
1.722 - std::queue<Node> Q; //queue of the totally unscanned nodes
1.723 - Q.push(v);
1.724 - std::queue<Node> R;
1.725 - //queue of the nodes which must be scanned for a possible shrink
1.726 -
1.727 - while ( !Q.empty() ) {
1.728 - Node x=Q.front();
1.729 - Q.pop();
1.730 - for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
1.731 - Node y=g.runningNode(e);
1.732 - //growOrAugment grows if y is covered by the matching and
1.733 - //augments if not. In this latter case it returns 1.
1.734 - if ( position[y]==C && growOrAugment(y, x, ear, blossom, tree, Q) )
1.735 - return;
1.736 - }
1.737 - R.push(x);
1.738 - }
1.739 -
1.740 - while ( !R.empty() ) {
1.741 - Node x=R.front();
1.742 - R.pop();
1.743 -
1.744 - for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
1.745 - Node y=g.runningNode(e);
1.746 -
1.747 - if ( position[y] == D && blossom.find(x) != blossom.find(y) )
1.748 - //Recall that we have only one tree.
1.749 - shrink( x, y, ear, blossom, tree, Q);
1.750 -
1.751 - while ( !Q.empty() ) {
1.752 - Node z=Q.front();
1.753 - Q.pop();
1.754 - for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
1.755 - Node w=g.runningNode(f);
1.756 - //growOrAugment grows if y is covered by the matching and
1.757 - //augments if not. In this latter case it returns 1.
1.758 - if ( position[w]==C && growOrAugment(w, z, ear, blossom, tree, Q) )
1.759 - return;
1.760 - }
1.761 - R.push(z);
1.762 - }
1.763 - } //for e
1.764 - } // while ( !R.empty() )
1.765 - }
1.766 -
1.767 -
1.768 - template <typename Graph>
1.769 - void MaxMatching<Graph>::normShrink(Node v,
1.770 - typename Graph::template
1.771 - NodeMap<Node>& ear,
1.772 - UFE& blossom, UFE& tree) {
1.773 - //We have one tree, which we grow and shrink. If we augment then we
1.774 - //return to the "for" cycle of runEdmonds().
1.775 -
1.776 - std::queue<Node> Q; //queue of the unscanned nodes
1.777 - Q.push(v);
1.778 - while ( !Q.empty() ) {
1.779 -
1.780 - Node x=Q.front();
1.781 - Q.pop();
1.782 -
1.783 - for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
1.784 - Node y=g.runningNode(e);
1.785 -
1.786 - switch ( position[y] ) {
1.787 - case D: //x and y must be in the same tree
1.788 - if ( blossom.find(x) != blossom.find(y) )
1.789 - //x and y are in the same tree
1.790 - shrink( x, y, ear, blossom, tree, Q);
1.791 - break;
1.792 - case C:
1.793 - //growOrAugment grows if y is covered by the matching and
1.794 - //augments if not. In this latter case it returns 1.
1.795 - if ( growOrAugment(y, x, ear, blossom, tree, Q) ) return;
1.796 - break;
1.797 - default: break;
1.798 - }
1.799 - }
1.800 - }
1.801 - }
1.802 -
1.803 -
1.804 - template <typename Graph>
1.805 - void MaxMatching<Graph>::shrink(Node x,Node y, typename
1.806 - Graph::template NodeMap<Node>& ear,
1.807 - UFE& blossom, UFE& tree, std::queue<Node>& Q) {
1.808 - //x and y are the two adjacent vertices in two blossoms.
1.809 -
1.810 - typename Graph::template NodeMap<bool> path(g,false);
1.811 -
1.812 - Node b=blossom.find(x);
1.813 - path.set(b,true);
1.814 - b=_mate[b];
1.815 - while ( b!=INVALID ) {
1.816 - b=blossom.find(ear[b]);
1.817 - path.set(b,true);
1.818 - b=_mate[b];
1.819 - } //we go until the root through bases of blossoms and odd vertices
1.820 -
1.821 - Node top=y;
1.822 - Node middle=blossom.find(top);
1.823 - Node bottom=x;
1.824 - while ( !path[middle] )
1.825 - shrinkStep(top, middle, bottom, ear, blossom, tree, Q);
1.826 - //Until we arrive to a node on the path, we update blossom, tree
1.827 - //and the positions of the odd nodes.
1.828 -
1.829 - Node base=middle;
1.830 - top=x;
1.831 - middle=blossom.find(top);
1.832 - bottom=y;
1.833 - Node blossom_base=blossom.find(base);
1.834 - while ( middle!=blossom_base )
1.835 - shrinkStep(top, middle, bottom, ear, blossom, tree, Q);
1.836 - //Until we arrive to a node on the path, we update blossom, tree
1.837 - //and the positions of the odd nodes.
1.838 -
1.839 - blossom.makeRep(base);
1.840 - }
1.841 -
1.842 -
1.843 -
1.844 - template <typename Graph>
1.845 - void MaxMatching<Graph>::shrinkStep(Node& top, Node& middle, Node& bottom,
1.846 - typename Graph::template
1.847 - NodeMap<Node>& ear,
1.848 - UFE& blossom, UFE& tree,
1.849 - std::queue<Node>& Q) {
1.850 - //We traverse a blossom and update everything.
1.851 -
1.852 - ear.set(top,bottom);
1.853 - Node t=top;
1.854 - while ( t!=middle ) {
1.855 - Node u=_mate[t];
1.856 - t=ear[u];
1.857 - ear.set(t,u);
1.858 - }
1.859 - bottom=_mate[middle];
1.860 - position.set(bottom,D);
1.861 - Q.push(bottom);
1.862 - top=ear[bottom];
1.863 - Node oldmiddle=middle;
1.864 - middle=blossom.find(top);
1.865 - tree.erase(bottom);
1.866 - tree.erase(oldmiddle);
1.867 - blossom.insert(bottom);
1.868 - blossom.join(bottom, oldmiddle);
1.869 - blossom.join(top, oldmiddle);
1.870 - }
1.871 -
1.872 -
1.873 - template <typename Graph>
1.874 - bool MaxMatching<Graph>::growOrAugment(Node& y, Node& x, typename Graph::template
1.875 - NodeMap<Node>& ear, UFE& blossom, UFE& tree,
1.876 - std::queue<Node>& Q) {
1.877 - //x is in a blossom in the tree, y is outside. If y is covered by
1.878 - //the matching we grow, otherwise we augment. In this case we
1.879 - //return 1.
1.880 -
1.881 - if ( _mate[y]!=INVALID ) { //grow
1.882 - ear.set(y,x);
1.883 - Node w=_mate[y];
1.884 - blossom.insert(w);
1.885 - position.set(y,A);
1.886 - position.set(w,D);
1.887 - tree.insert(y);
1.888 - tree.insert(w);
1.889 - tree.join(y,blossom.find(x));
1.890 - tree.join(w,y);
1.891 - Q.push(w);
1.892 - } else { //augment
1.893 - augment(x, ear, blossom, tree);
1.894 - _mate.set(x,y);
1.895 - _mate.set(y,x);
1.896 - return true;
1.897 - }
1.898 - return false;
1.899 - }
1.900 -
1.901 -
1.902 - template <typename Graph>
1.903 - void MaxMatching<Graph>::augment(Node x,
1.904 - typename Graph::template NodeMap<Node>& ear,
1.905 - UFE& blossom, UFE& tree) {
1.906 - Node v=_mate[x];
1.907 - while ( v!=INVALID ) {
1.908 -
1.909 - Node u=ear[v];
1.910 - _mate.set(v,u);
1.911 - Node tmp=v;
1.912 - v=_mate[u];
1.913 - _mate.set(u,tmp);
1.914 - }
1.915 - Node y=blossom.find(x);
1.916 - for (typename UFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {
1.917 - if ( position[tit] == D ) {
1.918 - for (typename UFE::ItemIt bit(blossom, tit); bit != INVALID; ++bit) {
1.919 - position.set( bit ,C);
1.920 - }
1.921 - blossom.eraseClass(tit);
1.922 - } else position.set( tit ,C);
1.923 - }
1.924 - tree.eraseClass(y);
1.925 -
1.926 - }
1.927 -
1.928
1.929 } //END OF NAMESPACE LEMON
1.930