COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/max_matching.h @ 2193:8057a4245685

Last change on this file since 2193:8057a4245685 was 2042:bdc953f2a449, checked in by Balazs Dezso, 18 years ago

New Algorithm group for matchings

LaTeX formulas
Bug fix => /\f$ will cause parsing error in doxygen

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