COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/max_matching.h @ 2476:059dcdda37c5

Last change on this file since 2476:059dcdda37c5 was 2391:14a343be7a5a, checked in by Alpar Juttner, 13 years ago

Happy New Year to all source files!

  • Property exe set to *
File size: 16.9 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2007
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_MAX_MATCHING_H
20#define LEMON_MAX_MATCHING_H
21
22#include <queue>
23#include <lemon/bits/invalid.h>
24#include <lemon/unionfind.h>
25#include <lemon/graph_utils.h>
26
27///\ingroup matching
28///\file
29///\brief Maximum matching algorithm in undirected graph.
30
31namespace lemon {
32
33  /// \ingroup matching
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
51  ///writePos after running the algorithm.
52  ///
53  ///\param Graph The undirected graph type the algorithm runs on.
54  ///
55  ///\author Jacint Szabo 
56  template <typename Graph>
57  class MaxMatching {
58
59  protected:
60
61    typedef typename Graph::Node Node;
62    typedef typename Graph::Edge Edge;
63    typedef typename Graph::UEdge UEdge;
64    typedef typename Graph::UEdgeIt UEdgeIt;
65    typedef typename Graph::NodeIt NodeIt;
66    typedef typename Graph::IncEdgeIt IncEdgeIt;
67
68    typedef typename Graph::template NodeMap<int> UFECrossRef;
69    typedef UnionFindEnum<UFECrossRef> UFE;
70
71  public:
72   
73    ///Indicates the Gallai-Edmonds decomposition of the graph.
74
75    ///Indicates the Gallai-Edmonds decomposition of the graph, which
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
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.
80    enum pos_enum {
81      D=0,
82      A=1,
83      C=2
84    };
85
86  protected:
87
88    static const int HEUR_density=2;
89    const Graph& g;
90    typename Graph::template NodeMap<Node> _mate;
91    typename Graph::template NodeMap<pos_enum> position;
92     
93  public:
94   
95    MaxMatching(const Graph& _g) : g(_g), _mate(_g,INVALID), position(_g) {}
96
97    ///Runs Edmonds' algorithm.
98
99    ///Runs Edmonds' algorithm for sparse graphs (number of edges <
100    ///2*number of nodes), and a heuristical Edmonds' algorithm with a
101    ///heuristic of postponing shrinks for dense graphs.
102    void run() {
103      if ( countUEdges(g) < HEUR_density*countNodes(g) ) {
104        greedyMatching();
105        runEdmonds(0);
106      } else runEdmonds(1);
107    }
108
109
110    ///Runs Edmonds' algorithm.
111   
112    ///If heur=0 it runs Edmonds' algorithm. If heur=1 it runs
113    ///Edmonds' algorithm with a heuristic of postponing shrinks,
114    ///giving a faster algorithm for dense graphs. 
115    void runEdmonds( int heur = 1 ) {
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)
157        if ( _mate[v]==INVALID ) {
158          for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
159            Node y=g.runningNode(e);
160            if ( _mate[y]==INVALID && y!=v ) {
161              _mate.set(v,y);
162              _mate.set(y,v);
163              break;
164            }
165          }
166        }
167    }
168
169    ///Returns the size of the actual matching stored.
170
171    ///Returns the size of the actual matching stored. After \ref
172    ///run() it returns the size of a maximum matching in the graph.
173    int size() const {
174      int s=0;
175      for(NodeIt v(g); v!=INVALID; ++v) {
176        if ( _mate[v]!=INVALID ) {
177          ++s;
178        }
179      }
180      return s/2;
181    }
182
183
184    ///Resets the actual matching to the empty matching.
185
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.
196    ///Returns INVALID if the \c node is not covered by the actual matching.
197    Node mate(Node& node) const {
198      return _mate[node];
199    }
200
201    ///Reads a matching from a \c Node valued \c Node map.
202
203    ///Reads a matching from a \c Node valued \c Node map. This map
204    ///must be \e symmetric, i.e. if \c map[u]==v then \c map[v]==u
205    ///must hold, and \c uv will be an edge of the matching.
206    template<typename NMapN>
207    void readNMapNode(NMapN& map) {
208      for(NodeIt v(g); v!=INVALID; ++v) {
209        _mate.set(v,map[v]);   
210      }
211    }
212   
213    ///Writes the stored matching to a \c Node valued \c Node map.
214
215    ///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
217    ///map[v]==u will hold, and now \c uv is an edge of the matching.
218    template<typename NMapN>
219    void writeNMapNode (NMapN& map) const {
220      for(NodeIt v(g); v!=INVALID; ++v) {
221        map.set(v,_mate[v]);   
222      }
223    }
224
225    ///Reads a matching from an \c UEdge valued \c Node map.
226
227    ///Reads a matching from an \c UEdge valued \c Node map. \c
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
244    ///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])
246    ///== v then \c map[u]==map[v] holds, and now this edge is an edge
247    ///of the matching.
248    template<typename NMapE>
249    void writeNMapEdge (NMapE& map)  const {
250      typename Graph::template NodeMap<bool> todo(g,true);
251      for(NodeIt v(g); v!=INVALID; ++v) {
252        if ( todo[v] && _mate[v]!=INVALID ) {
253          Node u=_mate[v];
254          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
255            if ( g.runningNode(e) == u ) {
256              map.set(u,e);
257              map.set(v,e);
258              todo.set(u,false);
259              todo.set(v,false);
260              break;
261            }
262          }
263        }
264      }
265    }
266
267
268    ///Reads a matching from a \c bool valued \c Edge map.
269   
270    ///Reads a matching from a \c bool valued \c Edge map. This map
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
290    ///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
292    ///edges \c e with \c map[e]==true form the matching.
293    template<typename EMapB>
294    void writeEMapBool (EMapB& map) const {
295      for(UEdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
296
297      typename Graph::template NodeMap<bool> todo(g,true);
298      for(NodeIt v(g); v!=INVALID; ++v) {
299        if ( todo[v] && _mate[v]!=INVALID ) {
300          Node u=_mate[v];
301          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
302            if ( g.runningNode(e) == u ) {
303              map.set(e,true);
304              todo.set(u,false);
305              todo.set(v,false);
306              break;
307            }
308          }
309        }
310      }
311    }
312
313
314    ///Writes the canonical decomposition of the graph after running
315    ///the algorithm.
316
317    ///After calling any run methods of the class, it writes the
318    ///Gallai-Edmonds canonical decomposition of the graph. \c map
319    ///must be a node map of \ref pos_enum 's.
320    template<typename NMapEnum>
321    void writePos (NMapEnum& map) const {
322      for(NodeIt v(g); v!=INVALID; ++v)  map.set(v,position[v]);
323    }
324
325  private:
326
327 
328    void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear, 
329                    UFE& blossom, UFE& tree);
330
331    void normShrink(Node v, typename Graph::template NodeMap<Node>& ear, 
332                    UFE& blossom, UFE& tree);
333
334    void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear, 
335                UFE& blossom, UFE& tree,std::queue<Node>& Q);
336
337    void shrinkStep(Node& top, Node& middle, Node& bottom,
338                    typename Graph::template NodeMap<Node>& ear, 
339                    UFE& blossom, UFE& tree, std::queue<Node>& Q);
340
341    bool growOrAugment(Node& y, Node& x, typename Graph::template
342                       NodeMap<Node>& ear, UFE& blossom, UFE& tree,
343                       std::queue<Node>& Q);
344
345    void augment(Node x, typename Graph::template NodeMap<Node>& ear, 
346                 UFE& blossom, UFE& tree);
347
348  };
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 
571} //END OF NAMESPACE LEMON
572
573#endif //LEMON_MAX_MATCHING_H
Note: See TracBrowser for help on using the repository browser.