COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/lemon/max_matching.h @ 1080:568ce2be7fe0

Last change on this file since 1080:568ce2be7fe0 was 1077:784a8bc07316, checked in by jacint, 16 years ago

Edmonds matching alg

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