COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/max_matching.h @ 2543:a0443c411220

Last change on this file since 2543:a0443c411220 was 2505:1bb471764ab8, checked in by Balazs Dezso, 12 years ago

Redesign interface of MaxMatching? and UnionFindEnum?
New class ExtendFindEnum?

Faster MaxMatching?

  • Property exe set to *
File size: 18.6 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  ///\brief 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 some of init functions.
40  ///
41  ///The dual side of a matching is a map of the nodes to
42  ///MaxMatching::DecompType, having values \c D, \c A and \c C
43  ///showing the Gallai-Edmonds decomposition of the graph. The nodes
44  ///in \c D induce a graph with factor-critical components, the nodes
45  ///in \c A form the barrier, and the nodes in \c C induce a graph
46  ///having a perfect matching. This decomposition can be attained by
47  ///calling \c decomposition() after running the algorithm.
48  ///
49  ///\param Graph The undirected graph type the algorithm runs on.
50  ///
51  ///\author Jacint Szabo 
52  template <typename Graph>
53  class MaxMatching {
54
55  protected:
56
57    typedef typename Graph::Node Node;
58    typedef typename Graph::Edge Edge;
59    typedef typename Graph::UEdge UEdge;
60    typedef typename Graph::UEdgeIt UEdgeIt;
61    typedef typename Graph::NodeIt NodeIt;
62    typedef typename Graph::IncEdgeIt IncEdgeIt;
63
64    typedef typename Graph::template NodeMap<int> UFECrossRef;
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
71  public:
72   
73    ///\brief 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 DecompType \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 DecompType {
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<DecompType> position;
92     
93  public:
94   
95    MaxMatching(const Graph& _g)
96      : g(_g), _mate(_g), position(_g) {}
97
98    ///\brief Sets the actual matching to the empty matching.
99    ///
100    ///Sets the actual matching to the empty matching. 
101    ///
102    void init() {
103      for(NodeIt v(g); v!=INVALID; ++v) {
104        _mate.set(v,INVALID);     
105        position.set(v,C);
106      }
107    }
108
109    ///\brief Finds a greedy matching for initial matching.
110    ///
111    ///For initial matchig it finds a maximal greedy matching.
112    void greedyInit() {
113      for(NodeIt v(g); v!=INVALID; ++v) {
114        _mate.set(v,INVALID);     
115        position.set(v,C);
116      }
117      for(NodeIt v(g); v!=INVALID; ++v)
118        if ( _mate[v]==INVALID ) {
119          for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
120            Node y=g.runningNode(e);
121            if ( _mate[y]==INVALID && y!=v ) {
122              _mate.set(v,y);
123              _mate.set(y,v);
124              break;
125            }
126          }
127        }
128    }
129
130    ///\brief Initialize the matching from each nodes' mate.
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    ///
277    ///Returns the size of the actual matching stored. After \ref
278    ///run() it returns the size of a maximum matching in the graph.
279    int size() const {
280      int s=0;
281      for(NodeIt v(g); v!=INVALID; ++v) {
282        if ( _mate[v]!=INVALID ) {
283          ++s;
284        }
285      }
286      return s/2;
287    }
288
289
290    ///\brief Returns the mate of a node in the actual matching.
291    ///
292    ///Returns the mate of a \c node in the actual matching.
293    ///Returns INVALID if the \c node is not covered by the actual matching.
294    Node mate(const Node& node) const {
295      return _mate[node];
296    }
297
298    ///\brief Returns the matching edge incident to the given node.
299    ///
300    ///Returns the matching edge of a \c node in the actual matching.
301    ///Returns INVALID if the \c node is not covered by the actual matching.
302    UEdge matchingEdge(const Node& node) const {
303      if (_mate[node] == INVALID) return INVALID;
304      Node n = node < _mate[node] ? node : _mate[node];
305      for (IncEdgeIt e(g, n); e != INVALID; ++e) {
306        if (g.oppositeNode(n, e) == _mate[n]) {
307          return e;
308        }
309      }
310      return INVALID;
311    }
312
313    /// \brief Returns the class of the node in the Edmonds-Gallai
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    ///
331    ///Writes the stored matching to a \c Node valued \c Node map. The
332    ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
333    ///map[v]==u will hold, and now \c uv is an edge of the matching.
334    template <typename MateMap>
335    void mateMap(MateMap& map) const {
336      for(NodeIt v(g); v!=INVALID; ++v) {
337        map.set(v,_mate[v]);   
338      }
339    }
340   
341    ///\brief Gives back the matching in an \c UEdge valued \c Node
342    ///map.
343    ///
344    ///Writes the stored matching to an \c UEdge valued \c Node
345    ///map. \c map[v] will be an \c UEdge incident to \c v. This
346    ///map will have the property that if \c g.oppositeNode(u,map[u])
347    ///== v then \c map[u]==map[v] holds, and now this edge is an edge
348    ///of the matching.
349    template <typename MatchingMap>
350    void matchingMap(MatchingMap& map)  const {
351      typename Graph::template NodeMap<bool> todo(g,true);
352      for(NodeIt v(g); v!=INVALID; ++v) {
353        if (_mate[v]!=INVALID && v < _mate[v]) {
354          Node u=_mate[v];
355          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
356            if ( g.runningNode(e) == u ) {
357              map.set(u,e);
358              map.set(v,e);
359              todo.set(u,false);
360              todo.set(v,false);
361              break;
362            }
363          }
364        }
365      }
366    }
367
368
369    ///\brief Gives back the matching in a \c bool valued \c UEdge
370    ///map.
371    ///
372    ///Writes the matching stored to a \c bool valued \c Edge
373    ///map. This map will have the property that there are no two
374    ///incident edges \c e, \c f with \c map[e]==map[f]==true. The
375    ///edges \c e with \c map[e]==true form the matching.
376    template<typename MatchingMap>
377    void matching(MatchingMap& map) const {
378      for(UEdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
379
380      typename Graph::template NodeMap<bool> todo(g,true);
381      for(NodeIt v(g); v!=INVALID; ++v) {
382        if ( todo[v] && _mate[v]!=INVALID ) {
383          Node u=_mate[v];
384          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
385            if ( g.runningNode(e) == u ) {
386              map.set(e,true);
387              todo.set(u,false);
388              todo.set(v,false);
389              break;
390            }
391          }
392        }
393      }
394    }
395
396
397    ///\brief Returns the canonical decomposition of the graph after running
398    ///the algorithm.
399    ///
400    ///After calling any run methods of the class, it writes the
401    ///Gallai-Edmonds canonical decomposition of the graph. \c map
402    ///must be a node map of \ref DecompType 's.
403    template <typename DecompositionMap>
404    void decomposition(DecompositionMap& map) const {
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);   
418    }
419
420  private:
421
422 
423    void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear, 
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    }
473
474    void normShrink(Node v, typename Graph::template NodeMap<Node>& ear, 
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    }
505
506    void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear, 
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    }
541
542    void shrinkStep(Node& top, Node& middle, Node& bottom,
543                    typename Graph::template NodeMap<Node>& ear, 
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
568
569    bool growOrAugment(Node& y, Node& x, typename Graph::template
570                       NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree,
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    }
594
595    void augment(Node x, typename Graph::template NodeMap<Node>& ear, 
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    }
619
620  };
621 
622} //END OF NAMESPACE LEMON
623
624#endif //LEMON_MAX_MATCHING_H
Note: See TracBrowser for help on using the repository browser.