COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/jacint/matching.h @ 1057:4588f97ad91f

Last change on this file since 1057:4588f97ad91f was 1057:4588f97ad91f, checked in by jacint, 15 years ago

undirgrafbug

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