COIN-OR::LEMON - Graph Library

Changeset 2505:1bb471764ab8 in lemon-0.x


Ignore:
Timestamp:
10/30/07 21:21:10 (12 years ago)
Author:
Balazs Dezso
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3351
Message:

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

Faster MaxMatching?

Files:
4 edited

Legend:

Unmodified
Added
Removed
  • lemon/max_matching.h

    r2391 r2505  
    3131namespace lemon {
    3232
    33   /// \ingroup matching
    34 
    35   ///Edmonds' alternating forest maximum matching algorithm.
    36 
     33  ///\ingroup matching
     34
     35  ///\brief Edmonds' alternating forest maximum matching algorithm.
     36  ///
    3737  ///This class provides Edmonds' alternating forest matching
    3838  ///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.
     39  ///algorithm using some of init functions.
    4440  ///
    4541  ///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.
     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.
    5248  ///
    5349  ///\param Graph The undirected graph type the algorithm runs on.
     
    6864    typedef typename Graph::template NodeMap<int> UFECrossRef;
    6965    typedef UnionFindEnum<UFECrossRef> UFE;
     66    typedef std::vector<Node> NV;
     67
     68    typedef typename Graph::template NodeMap<int> EFECrossRef;
     69    typedef ExtendFindEnum<EFECrossRef> EFE;
    7070
    7171  public:
    7272   
    73     ///Indicates the Gallai-Edmonds decomposition of the graph.
    74 
     73    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
     74    ///
    7575    ///Indicates the Gallai-Edmonds decomposition of the graph, which
    7676    ///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
     77    ///nodes with DecompType \c D induce a graph with factor-critical
    7878    ///components, the nodes in \c A form the canonical barrier, and the
    7979    ///nodes in \c C induce a graph having a perfect matching.
    80     enum pos_enum {
     80    enum DecompType {
    8181      D=0,
    8282      A=1,
     
    8989    const Graph& g;
    9090    typename Graph::template NodeMap<Node> _mate;
    91     typename Graph::template NodeMap<pos_enum> position;
     91    typename Graph::template NodeMap<DecompType> position;
    9292     
    9393  public:
    9494   
    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() {
     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      }
    156117      for(NodeIt v(g); v!=INVALID; ++v)
    157118        if ( _mate[v]==INVALID ) {
     
    167128    }
    168129
    169     ///Returns the size of the actual matching stored.
    170 
     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    ///
    171277    ///Returns the size of the actual matching stored. After \ref
    172278    ///run() it returns the size of a maximum matching in the graph.
     
    182288
    183289
    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 
     290    ///\brief Returns the mate of a node in the actual matching.
     291    ///
    195292    ///Returns the mate of a \c node in the actual matching.
    196293    ///Returns INVALID if the \c node is not covered by the actual matching.
    197     Node mate(Node& node) const {
     294    Node mate(const Node& node) const {
    198295      return _mate[node];
    199296    }
    200297
    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       }
     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;
    211311    }
    212    
    213     ///Writes the stored matching to a \c Node valued \c Node map.
    214 
     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    ///
    215331    ///Writes the stored matching to a \c Node valued \c Node map. The
    216332    ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
    217333    ///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 {
     334    template <typename MateMap>
     335    void mateMap(MateMap& map) const {
    220336      for(NodeIt v(g); v!=INVALID; ++v) {
    221337        map.set(v,_mate[v]);   
    222338      }
    223339    }
    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 
     340   
     341    ///\brief Gives back the matching in an \c UEdge valued \c Node
     342    ///map.
     343    ///
    243344    ///Writes the stored matching to an \c UEdge valued \c Node
    244345    ///map. \c map[v] will be an \c UEdge incident to \c v. This
     
    246347    ///== v then \c map[u]==map[v] holds, and now this edge is an edge
    247348    ///of the matching.
    248     template<typename NMapE>
    249     void writeNMapEdge (NMapE& map)  const {
     349    template <typename MatchingMap>
     350    void matchingMap(MatchingMap& map)  const {
    250351      typename Graph::template NodeMap<bool> todo(g,true);
    251352      for(NodeIt v(g); v!=INVALID; ++v) {
    252         if ( todo[v] && _mate[v]!=INVALID ) {
     353        if (_mate[v]!=INVALID && v < _mate[v]) {
    253354          Node u=_mate[v];
    254355          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
     
    266367
    267368
    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 
     369    ///\brief Gives back the matching in a \c bool valued \c UEdge
     370    ///map.
     371    ///
    289372    ///Writes the matching stored to a \c bool valued \c Edge
    290373    ///map. This map will have the property that there are no two
    291374    ///incident edges \c e, \c f with \c map[e]==map[f]==true. The
    292375    ///edges \c e with \c map[e]==true form the matching.
    293     template<typename EMapB>
    294     void writeEMapBool (EMapB& map) const {
     376    template<typename MatchingMap>
     377    void matching(MatchingMap& map) const {
    295378      for(UEdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
    296379
     
    312395
    313396
    314     ///Writes the canonical decomposition of the graph after running
     397    ///\brief Returns the canonical decomposition of the graph after running
    315398    ///the algorithm.
    316 
     399    ///
    317400    ///After calling any run methods of the class, it writes the
    318401    ///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]);
     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);   
    323418    }
    324419
     
    327422 
    328423    void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear, 
    329                     UFE& blossom, UFE& tree);
     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    }
    330473
    331474    void normShrink(Node v, typename Graph::template NodeMap<Node>& ear, 
    332                     UFE& blossom, UFE& tree);
     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    }
    333505
    334506    void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear, 
    335                 UFE& blossom, UFE& tree,std::queue<Node>& Q);
     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    }
    336541
    337542    void shrinkStep(Node& top, Node& middle, Node& bottom,
    338543                    typename Graph::template NodeMap<Node>& ear, 
    339                     UFE& blossom, UFE& tree, std::queue<Node>& Q);
     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
    340568
    341569    bool growOrAugment(Node& y, Node& x, typename Graph::template
    342                        NodeMap<Node>& ear, UFE& blossom, UFE& tree,
    343                        std::queue<Node>& Q);
     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    }
    344594
    345595    void augment(Node x, typename Graph::template NodeMap<Node>& ear, 
    346                  UFE& blossom, UFE& tree);
     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    }
    347619
    348620  };
    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 
    570621 
    571622} //END OF NAMESPACE LEMON
  • lemon/unionfind.h

    r2427 r2505  
    179179  private:
    180180   
     181    ItemIntMap& index;
     182
    181183    // If the parent stores negative value for an item then that item
    182     // is root item and it has -items[it].parent component size.  Else
     184    // is root item and it has ~(items[it].parent) component id.  Else
    183185    // the items[it].parent contains the index of the parent.
    184186    //
    185     // The \c nextItem and \c prevItem provides the double-linked
    186     // cyclic list of one component's items. The \c prevClass and
    187     // \c nextClass gives the double linked list of the representant
    188     // items.
     187    // The \c next and \c prev provides the double-linked
     188    // cyclic list of one component's items.
    189189    struct ItemT {
    190190      int parent;
    191191      Item item;
    192192
    193       int nextItem, prevItem;
    194       int nextClass, prevClass;
     193      int next, prev;
    195194    };
    196195
    197196    std::vector<ItemT> items;
    198     ItemIntMap& index;
    199 
    200     int firstClass;
     197    int firstFreeItem;
     198
     199    struct ClassT {
     200      int size;
     201      int firstItem;
     202      int next, prev;
     203    };
     204   
     205    std::vector<ClassT> classes;
     206    int firstClass, firstFreeClass;
     207
     208    int newClass() {
     209      if (firstFreeClass == -1) {
     210        int cdx = classes.size();
     211        classes.push_back(ClassT());
     212        return cdx;
     213      } else {
     214        int cdx = firstFreeClass;
     215        firstFreeClass = classes[firstFreeClass].next;
     216        return cdx;
     217      }
     218    }
     219
     220    int newItem() {
     221      if (firstFreeItem == -1) {
     222        int idx = items.size();
     223        items.push_back(ItemT());
     224        return idx;
     225      } else {
     226        int idx = firstFreeItem;
     227        firstFreeItem = items[firstFreeItem].next;
     228        return idx;
     229      }
     230    }
    201231
    202232
     
    218248    }
    219249
    220     void unlaceClass(int k) {
    221       if (items[k].prevClass != -1) {
    222         items[items[k].prevClass].nextClass = items[k].nextClass;
     250    int classIndex(int idx) const {
     251      return ~(items[repIndex(idx)].parent);
     252    }
     253
     254    void singletonItem(int idx) {
     255      items[idx].next = idx;
     256      items[idx].prev = idx;
     257    }
     258
     259    void laceItem(int idx, int rdx) {
     260      items[idx].prev = rdx;
     261      items[idx].next = items[rdx].next;
     262      items[items[rdx].next].prev = idx;
     263      items[rdx].next = idx;
     264    }
     265
     266    void unlaceItem(int idx) {
     267      items[items[idx].prev].next = items[idx].next;
     268      items[items[idx].next].prev = items[idx].prev;
     269     
     270      items[idx].next = firstFreeItem;
     271      firstFreeItem = idx;
     272    }
     273
     274    void spliceItems(int ak, int bk) {
     275      items[items[ak].prev].next = bk;
     276      items[items[bk].prev].next = ak;
     277      int tmp = items[ak].prev;
     278      items[ak].prev = items[bk].prev;
     279      items[bk].prev = tmp;
     280       
     281    }
     282
     283    void laceClass(int cls) {
     284      if (firstClass != -1) {
     285        classes[firstClass].prev = cls;
     286      }
     287      classes[cls].next = firstClass;
     288      classes[cls].prev = -1;
     289      firstClass = cls;
     290    }
     291
     292    void unlaceClass(int cls) {
     293      if (classes[cls].prev != -1) {
     294        classes[classes[cls].prev].next = classes[cls].next;
    223295      } else {
    224         firstClass = items[k].nextClass;
    225       }
    226       if (items[k].nextClass != -1) {
    227         items[items[k].nextClass].prevClass = items[k].prevClass;
    228       }
     296        firstClass = classes[cls].next;
     297      }
     298      if (classes[cls].next != -1) {
     299        classes[classes[cls].next].prev = classes[cls].prev;
     300      }
     301     
     302      classes[cls].next = firstFreeClass;
     303      firstFreeClass = cls;
    229304    }
    230305
    231     void spliceItems(int ak, int bk) {
    232       items[items[ak].prevItem].nextItem = bk;
    233       items[items[bk].prevItem].nextItem = ak;
    234       int tmp = items[ak].prevItem;
    235       items[ak].prevItem = items[bk].prevItem;
    236       items[bk].prevItem = tmp;
    237        
    238     }
    239 
    240306  public:
    241307
    242308    UnionFindEnum(ItemIntMap& _index)
    243       : items(), index(_index), firstClass(-1) {}
     309      : index(_index), items(), firstFreeItem(-1),
     310        firstClass(-1), firstFreeClass(-1) {}
    244311   
    245312    /// \brief Inserts the given element into a new component.
     
    248315    /// given element.
    249316    ///
    250     void insert(const Item& item) {
    251       ItemT t;
    252 
    253       int idx = items.size();
     317    int insert(const Item& item) {
     318      int idx = newItem();
     319
    254320      index.set(item, idx);
    255321
    256       t.nextItem = idx;
    257       t.prevItem = idx;
    258       t.item = item;
    259       t.parent = -1;
    260      
    261       t.nextClass = firstClass;
    262       if (firstClass != -1) {
    263         items[firstClass].prevClass = idx;
    264       }
    265       t.prevClass = -1;
    266       firstClass = idx;
    267 
    268       items.push_back(t);
     322      singletonItem(idx);
     323      items[idx].item = item;
     324
     325      int cdx = newClass();
     326
     327      items[idx].parent = ~cdx;
     328
     329      laceClass(cdx);
     330      classes[cdx].size = 1;
     331      classes[cdx].firstItem = idx;
     332
     333      firstClass = cdx;
     334     
     335      return cdx;
    269336    }
    270337
     
    273340    /// This methods inserts the element \e a into the component of the
    274341    /// element \e comp.
    275     void insert(const Item& item, const Item& comp) {
    276       int k = repIndex(index[comp]);
    277       ItemT t;
    278 
    279       int idx = items.size();
     342    void insert(const Item& item, int cls) {
     343      int rdx = classes[cls].firstItem;
     344      int idx = newItem();
     345
    280346      index.set(item, idx);
    281347
    282       t.prevItem = k;
    283       t.nextItem = items[k].nextItem;
    284       items[items[k].nextItem].prevItem = idx;
    285       items[k].nextItem = idx;
    286 
    287       t.item = item;
    288       t.parent = k;
    289 
    290       --items[k].parent;
    291 
    292       items.push_back(t);
     348      laceItem(idx, rdx);
     349
     350      items[idx].item = item;
     351      items[idx].parent = rdx;
     352
     353      ++classes[~(items[rdx].parent)].size;
    293354    }
    294355
     
    299360      items.clear();
    300361      firstClass = -1;
    301     }
    302 
    303     /// \brief Finds the leader of the component of the given element.
    304     ///
    305     /// The method returns the leader of the component of the given element.
    306     const Item& find(const Item &item) const {
    307       return items[repIndex(index[item])].item;
     362      firstFreeItem = -1;
     363    }
     364
     365    /// \brief Finds the component of the given element.
     366    ///
     367    /// The method returns the component id of the given element.
     368    int find(const Item &item) const {
     369      return ~(items[repIndex(index[item])].parent);
    308370    }
    309371
     
    313375    /// Joins the component of element \e a and component of
    314376    /// element \e b. If \e a and \e b are in the same component then
    315     /// returns false else returns true.
    316     bool join(const Item& a, const Item& b) {
     377    /// returns -1 else returns the remaining class.
     378    int join(const Item& a, const Item& b) {
    317379
    318380      int ak = repIndex(index[a]);
     
    320382
    321383      if (ak == bk) {
    322         return false;
    323       }
    324 
    325       if ( items[ak].parent < items[bk].parent ) {
    326         unlaceClass(bk);
    327         items[ak].parent += items[bk].parent;
     384        return -1;
     385      }
     386
     387      int acx = ~(items[ak].parent);
     388      int bcx = ~(items[bk].parent);
     389
     390      int rcx;
     391
     392      if (classes[acx].size > classes[bcx].size) {
     393        classes[acx].size += classes[bcx].size;
    328394        items[bk].parent = ak;
     395        unlaceClass(bcx);
     396        rcx = acx;
    329397      } else {
    330         unlaceClass(ak);
    331         items[bk].parent += items[ak].parent;
     398        classes[bcx].size += classes[acx].size;
    332399        items[ak].parent = bk;
     400        unlaceClass(acx);
     401        rcx = bcx;
    333402      }
    334403      spliceItems(ak, bk);
    335404
    336       return true;
    337     }
    338 
    339     /// \brief Returns the size of the component of element \e a.
    340     ///
    341     /// Returns the size of the component of element \e a.
    342     int size(const Item &item) const {
    343       return - items[repIndex(index[item])].parent;
    344     }
    345 
    346     /// \brief Splits up the component of the element.
    347     ///
    348     /// Splitting the component of the element into sigleton
    349     /// components (component of size one).
    350     void split(const Item &item) {
    351       int k = repIndex(index[item]);
    352       int idx = items[k].nextItem;
    353       while (idx != k) {
    354         int next = items[idx].nextItem;
     405      return rcx;
     406    }
     407
     408    /// \brief Returns the size of the class.
     409    ///
     410    /// Returns the size of the class.
     411    int size(int cls) const {
     412      return classes[cls].size;
     413    }
     414
     415    /// \brief Splits up the component.
     416    ///
     417    /// Splitting the component into singleton components (component
     418    /// of size one).
     419    void split(int cls) {
     420      int fdx = classes[cls].firstItem;
     421      int idx = items[fdx].next;
     422      while (idx != fdx) {
     423        int next = items[idx].next;
     424
     425        singletonItem(idx);
     426
     427        int cdx = newClass();       
     428        items[idx].parent = ~cdx;
     429
     430        laceClass(cdx);
     431        classes[cdx].size = 1;
     432        classes[cdx].firstItem = idx;
    355433       
    356         items[idx].parent = -1;
    357         items[idx].prevItem = idx;
    358         items[idx].nextItem = idx;
    359        
    360         items[idx].nextClass = firstClass;
    361         items[firstClass].prevClass = idx;
    362         firstClass = idx;
    363 
    364434        idx = next;
    365435      }
    366436
    367       items[idx].parent = -1;
    368       items[idx].prevItem = idx;
    369       items[idx].nextItem = idx;
    370      
    371       items[firstClass].prevClass = -1;
    372     }
    373 
    374     /// \brief Sets the given element to the leader element of its component.
    375     ///
    376     /// Sets the given element to the leader element of its component.
    377     void makeRep(const Item &item) {
    378       int nk = index[item];
    379       int k = repIndex(nk);
    380       if (nk == k) return;
    381      
    382       if (items[k].prevClass != -1) {
    383         items[items[k].prevClass].nextClass = nk;
    384       } else {
    385         firstClass = nk;
    386       }
    387       if (items[k].nextClass != -1) {
    388         items[items[k].nextClass].prevClass = nk;
    389       }
    390      
    391       int idx = items[k].nextItem;
    392       while (idx != k) {
    393         items[idx].parent = nk;
    394         idx = items[idx].nextItem;
    395       }
    396      
    397       items[nk].parent = items[k].parent;
    398       items[k].parent = nk;
     437      items[idx].prev = idx;
     438      items[idx].next = idx;
     439
     440      classes[~(items[idx].parent)].size = 1;
     441     
    399442    }
    400443
     
    406449    /// \warning It is an error to remove an element which is not in
    407450    /// the structure.
    408     void erase(const Item &item) {
     451    /// \warning This running time of this operation is proportional to the
     452    /// number of the items in this class.
     453    void erase(const Item& item) {
    409454      int idx = index[item];
    410       if (rep(idx)) {
    411         int k = idx;
    412         if (items[k].parent == -1) {
    413           unlaceClass(idx);
    414           return;
    415         } else {
    416           int nk = items[k].nextItem;
    417           if (items[k].prevClass != -1) {
    418             items[items[k].prevClass].nextClass = nk;
    419           } else {
    420             firstClass = nk;
    421           }
    422           if (items[k].nextClass != -1) {
    423             items[items[k].nextClass].prevClass = nk;
    424           }
    425      
    426           int l = items[k].nextItem;
    427           while (l != k) {
    428             items[l].parent = nk;
    429             l = items[l].nextItem;
    430           }
     455      int fdx = items[idx].next;
     456
     457      int cdx = classIndex(idx);
     458      if (idx == fdx) {
     459        unlaceClass(cdx);
     460        items[idx].next = firstFreeItem;
     461        firstFreeItem = idx;
     462        return;
     463      } else {
     464        classes[cdx].firstItem = fdx;
     465        --classes[cdx].size;
     466        items[fdx].parent = ~cdx;
     467
     468        unlaceItem(idx);
     469        idx = items[fdx].next;
     470        while (idx != fdx) {
     471          items[idx].parent = fdx;
     472          idx = items[idx].next;
     473        }
    431474         
    432           items[nk].parent = items[k].parent + 1;
    433         }
    434       } else {
    435        
    436         int k = repIndex(idx);
    437         idx = items[k].nextItem;
    438         while (idx != k) {
    439           items[idx].parent = k;
    440           idx = items[idx].nextItem;
    441         }
    442 
    443         ++items[k].parent;
    444       }
    445 
    446       idx = index[item];     
    447       items[items[idx].prevItem].nextItem = items[idx].nextItem;
    448       items[items[idx].nextItem].prevItem = items[idx].prevItem;
     475      }
    449476
    450477    }   
    451 
    452     /// \brief Moves the given element to another component.
    453     ///
    454     /// This method moves the element \e a from its component
    455     /// to the component of \e comp.
    456     /// If \e a and \e comp are in the same component then
    457     /// it returns false otherwise it returns true.
    458     bool move(const Item &item, const Item &comp) {
    459       if (repIndex(index[item]) == repIndex(index[comp])) return false;
    460       erase(item);
    461       insert(item, comp);
    462       return true;
    463     }
    464 
    465478
    466479    /// \brief Removes the component of the given element from the structure.
     
    470483    /// \warning It is an error to give an element which is not in the
    471484    /// structure.
    472     void eraseClass(const Item &item) {
    473       unlaceClass(repIndex(index[item]));
     485    void eraseClass(int cls) {
     486      int fdx = classes[cls].firstItem;
     487      unlaceClass(cls);
     488      items[items[fdx].prev].next = firstFreeItem;
     489      firstFreeItem = fdx;
    474490    }
    475491
     
    477493    ///
    478494    /// ClassIt is a lemon style iterator for the components. It iterates
    479     /// on the representant items of the classes.
     495    /// on the ids of the classes.
    480496    class ClassIt {
    481497    public:
     
    484500      /// Constructor of the iterator
    485501      ClassIt(const UnionFindEnum& ufe) : unionFind(&ufe) {
    486         idx = unionFind->firstClass;
     502        cdx = unionFind->firstClass;
    487503      }
    488504
     
    490506      ///
    491507      /// Constructor to get invalid iterator
    492       ClassIt(Invalid) : unionFind(0), idx(-1) {}
     508      ClassIt(Invalid) : unionFind(0), cdx(-1) {}
    493509     
    494510      /// \brief Increment operator
     
    496512      /// It steps to the next representant item.
    497513      ClassIt& operator++() {
    498         idx = unionFind->items[idx].nextClass;
     514        cdx = unionFind->classes[cdx].next;
    499515        return *this;
    500516      }
     
    503519      ///
    504520      /// It converts the iterator to the current representant item.
    505       operator const Item&() const {
    506         return unionFind->items[idx].item;
     521      operator int() const {
     522        return cdx;
    507523      }
    508524
     
    511527      /// Equality operator
    512528      bool operator==(const ClassIt& i) {
    513         return i.idx == idx;
     529        return i.cdx == cdx;
    514530      }
    515531
     
    518534      /// Inequality operator
    519535      bool operator!=(const ClassIt& i) {
    520         return i.idx != idx;
     536        return i.cdx != cdx;
    521537      }
    522538     
    523539    private:
    524540      const UnionFindEnum* unionFind;
    525       int idx;
     541      int cdx;
    526542    };
    527543
     
    546562      /// Constructor of the iterator. The iterator iterates
    547563      /// on the class of the \c item.
    548       ItemIt(const UnionFindEnum& ufe, const Item& item) : unionFind(&ufe) {
    549         idx = unionFind->repIndex(unionFind->index[item]);
     564      ItemIt(const UnionFindEnum& ufe, int cls) : unionFind(&ufe) {
     565        fdx = idx = unionFind->classes[cls].firstItem;
    550566      }
    551567
     
    559575      /// It steps to the next item in the class.
    560576      ItemIt& operator++() {
    561         idx = unionFind->items[idx].nextItem;
    562         if (unionFind->rep(idx)) idx = -1;
     577        idx = unionFind->items[idx].next;
     578        if (idx == fdx) idx = -1;
    563579        return *this;
    564580      }
     
    587603    private:
    588604      const UnionFindEnum* unionFind;
    589       int idx;
     605      int idx, fdx;
    590606    };
    591607
    592608  };
    593609
     610  /// \ingroup auxdat
     611  ///
     612  /// \brief A \e Extend-Find data structure implementation which
     613  /// is able to enumerate the components.
     614  ///
     615  /// The class implements an \e Extend-Find data structure which is
     616  /// able to enumerate the components and the items in a
     617  /// component. The data structure is a simplification of the
     618  /// Union-Find structure, and it does not allow to merge two components.
     619  ///
     620  /// \pre You need to add all the elements by the \ref insert()
     621  /// method.
     622  template <typename _ItemIntMap>
     623  class ExtendFindEnum {
     624  public:
     625   
     626    typedef _ItemIntMap ItemIntMap;
     627    typedef typename ItemIntMap::Key Item;
     628
     629  private:
     630   
     631    ItemIntMap& index;
     632
     633    struct ItemT {
     634      int cls;
     635      Item item;
     636      int next, prev;
     637    };
     638
     639    std::vector<ItemT> items;
     640    int firstFreeItem;
     641
     642    struct ClassT {
     643      int firstItem;
     644      int next, prev;
     645    };
     646
     647    std::vector<ClassT> classes;
     648
     649    int firstClass, firstFreeClass;
     650
     651    int newClass() {
     652      if (firstFreeClass != -1) {
     653        int cdx = firstFreeClass;
     654        firstFreeClass = classes[cdx].next;
     655        return cdx;
     656      } else {
     657        classes.push_back(ClassT());
     658        return classes.size() - 1;
     659      }
     660    }
     661
     662    int newItem() {
     663      if (firstFreeItem != -1) {
     664        int idx = firstFreeItem;
     665        firstFreeItem = items[idx].next;
     666        return idx;
     667      } else {
     668        items.push_back(ItemT());
     669        return items.size() - 1;
     670      }
     671    }
     672
     673  public:
     674
     675    /// \brief Constructor
     676    ExtendFindEnum(ItemIntMap& _index)
     677      : index(_index), items(), firstFreeItem(-1),
     678        classes(), firstClass(-1), firstFreeClass(-1) {}
     679   
     680    /// \brief Inserts the given element into a new component.
     681    ///
     682    /// This method creates a new component consisting only of the
     683    /// given element.
     684    int insert(const Item& item) {
     685      int cdx = newClass();
     686      classes[cdx].prev = -1;
     687      classes[cdx].next = firstClass;
     688      firstClass = cdx;
     689     
     690      int idx = newItem();
     691      items[idx].item = item;
     692      items[idx].cls = cdx;
     693      items[idx].prev = idx;
     694      items[idx].next = idx;
     695
     696      classes[cdx].firstItem = idx;
     697
     698      index.set(item, idx);
     699     
     700      return cdx;
     701    }
     702
     703    /// \brief Inserts the given element into the given component.
     704    ///
     705    /// This methods inserts the element \e item a into the \e cls class.
     706    void insert(const Item& item, int cls) {
     707      int idx = newItem();
     708      int rdx = classes[cls].firstItem;
     709      items[idx].item = item;
     710      items[idx].cls = cls;
     711
     712      items[idx].prev = rdx;
     713      items[idx].next = items[rdx].next;
     714      items[items[rdx].next].prev = idx;
     715      items[rdx].next = idx;
     716
     717      index.set(item, idx);
     718    }
     719
     720    /// \brief Clears the union-find data structure
     721    ///
     722    /// Erase each item from the data structure.
     723    void clear() {
     724      items.clear();
     725      classes.clear;
     726      firstClass = firstFreeClass = firstFreeItem = -1;
     727    }
     728
     729    /// \brief Gives back the class of the \e item.
     730    ///
     731    /// Gives back the class of the \e item.
     732    int find(const Item &item) const {
     733      return items[index[item]].cls;
     734    }
     735   
     736    /// \brief Removes the given element from the structure.
     737    ///
     738    /// Removes the element from its component and if the component becomes
     739    /// empty then removes that component from the component list.
     740    ///
     741    /// \warning It is an error to remove an element which is not in
     742    /// the structure.
     743    void erase(const Item &item) {
     744      int idx = index[item];
     745      int cdx = items[idx].cls;
     746     
     747      if (idx == items[idx].next) {
     748        if (classes[cdx].prev != -1) {
     749          classes[classes[cdx].prev].next = classes[cdx].next;
     750        } else {
     751          firstClass = classes[cdx].next;
     752        }
     753        if (classes[cdx].next != -1) {
     754          classes[classes[cdx].next].prev = classes[cdx].prev;
     755        }
     756        classes[cdx].next = firstFreeClass;
     757        firstFreeClass = cdx;
     758      } else {
     759        classes[cdx].firstItem = items[idx].next;
     760        items[items[idx].next].prev = items[idx].prev;
     761        items[items[idx].prev].next = items[idx].next;
     762      }
     763      items[idx].next = firstFreeItem;
     764      firstFreeItem = idx;
     765       
     766    }   
     767
     768   
     769    /// \brief Removes the component of the given element from the structure.
     770    ///
     771    /// Removes the component of the given element from the structure.
     772    ///
     773    /// \warning It is an error to give an element which is not in the
     774    /// structure.
     775    void eraseClass(int cdx) {
     776      int idx = classes[cdx].firstItem;
     777      items[items[idx].prev].next = firstFreeItem;
     778      firstFreeItem = idx;
     779
     780      if (classes[cdx].prev != -1) {
     781        classes[classes[cdx].prev].next = classes[cdx].next;
     782      } else {
     783        firstClass = classes[cdx].next;
     784      }
     785      if (classes[cdx].next != -1) {
     786        classes[classes[cdx].next].prev = classes[cdx].prev;
     787      }
     788      classes[cdx].next = firstFreeClass;
     789      firstFreeClass = cdx;
     790    }
     791
     792    /// \brief Lemon style iterator for the classes.
     793    ///
     794    /// ClassIt is a lemon style iterator for the components. It iterates
     795    /// on the ids of classes.
     796    class ClassIt {
     797    public:
     798      /// \brief Constructor of the iterator
     799      ///
     800      /// Constructor of the iterator
     801      ClassIt(const ExtendFindEnum& ufe) : extendFind(&ufe) {
     802        cdx = extendFind->firstClass;
     803      }
     804
     805      /// \brief Constructor to get invalid iterator
     806      ///
     807      /// Constructor to get invalid iterator
     808      ClassIt(Invalid) : extendFind(0), cdx(-1) {}
     809     
     810      /// \brief Increment operator
     811      ///
     812      /// It steps to the next representant item.
     813      ClassIt& operator++() {
     814        cdx = extendFind->classes[cdx].next;
     815        return *this;
     816      }
     817     
     818      /// \brief Conversion operator
     819      ///
     820      /// It converts the iterator to the current class id.
     821      operator int() const {
     822        return cdx;
     823      }
     824
     825      /// \brief Equality operator
     826      ///
     827      /// Equality operator
     828      bool operator==(const ClassIt& i) {
     829        return i.cdx == cdx;
     830      }
     831
     832      /// \brief Inequality operator
     833      ///
     834      /// Inequality operator
     835      bool operator!=(const ClassIt& i) {
     836        return i.cdx != cdx;
     837      }
     838     
     839    private:
     840      const ExtendFindEnum* extendFind;
     841      int cdx;
     842    };
     843
     844    /// \brief Lemon style iterator for the items of a component.
     845    ///
     846    /// ClassIt is a lemon style iterator for the components. It iterates
     847    /// on the items of a class. By example if you want to iterate on
     848    /// each items of each classes then you may write the next code.
     849    ///\code
     850    /// for (ClassIt cit(ufe); cit != INVALID; ++cit) {
     851    ///   std::cout << "Class: ";
     852    ///   for (ItemIt iit(ufe, cit); iit != INVALID; ++iit) {
     853    ///     std::cout << toString(iit) << ' ' << std::endl;
     854    ///   }
     855    ///   std::cout << std::endl;
     856    /// }
     857    ///\endcode
     858    class ItemIt {
     859    public:
     860      /// \brief Constructor of the iterator
     861      ///
     862      /// Constructor of the iterator. The iterator iterates
     863      /// on the class of the \c item.
     864      ItemIt(const ExtendFindEnum& ufe, int cls) : extendFind(&ufe) {
     865        fdx = idx = extendFind->classes[cls].firstItem;
     866      }
     867
     868      /// \brief Constructor to get invalid iterator
     869      ///
     870      /// Constructor to get invalid iterator
     871      ItemIt(Invalid) : extendFind(0), idx(-1) {}
     872     
     873      /// \brief Increment operator
     874      ///
     875      /// It steps to the next item in the class.
     876      ItemIt& operator++() {
     877        idx = extendFind->items[idx].next;
     878        if (fdx == idx) idx = -1;
     879        return *this;
     880      }
     881     
     882      /// \brief Conversion operator
     883      ///
     884      /// It converts the iterator to the current item.
     885      operator const Item&() const {
     886        return extendFind->items[idx].item;
     887      }
     888
     889      /// \brief Equality operator
     890      ///
     891      /// Equality operator
     892      bool operator==(const ItemIt& i) {
     893        return i.idx == idx;
     894      }
     895
     896      /// \brief Inequality operator
     897      ///
     898      /// Inequality operator
     899      bool operator!=(const ItemIt& i) {
     900        return i.idx != idx;
     901      }
     902     
     903    private:
     904      const ExtendFindEnum* extendFind;
     905      int idx, fdx;
     906    };
     907
     908  };
    594909
    595910  //! @}
  • test/max_matching_test.cc

    r2391 r2505  
    6767 
    6868  MaxMatching<Graph> max_matching(g);
    69   max_matching.runEdmonds(0);
     69  max_matching.init();
     70  max_matching.startDense();
    7071 
    7172  int s=0;
    7273  Graph::NodeMap<Node> mate(g,INVALID);
    73   max_matching.writeNMapNode(mate);
     74  max_matching.mateMap(mate);
    7475  for(NodeIt v(g); v!=INVALID; ++v) {
    7576    if ( mate[v]!=INVALID ) ++s;
     
    8384  check ( size == max_matching.size(), "mate() returns a different size matching than max_matching.size()" );
    8485
    85   Graph::NodeMap<MaxMatching<Graph>::pos_enum> pos0(g);
    86   max_matching.writePos(pos0);
     86  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos0(g);
     87  max_matching.decomposition(pos0);
    8788 
    88   max_matching.resetMatching();
    89   max_matching.runEdmonds(1);
     89  max_matching.init();
     90  max_matching.startSparse();
    9091  s=0;
    91   max_matching.writeNMapNode(mate);
     92  max_matching.mateMap(mate);
    9293  for(NodeIt v(g); v!=INVALID; ++v) {
    9394    if ( mate[v]!=INVALID ) ++s;
     
    9596  check ( int(s/2) == size, "The size does not equal!" );
    9697
    97   Graph::NodeMap<MaxMatching<Graph>::pos_enum> pos1(g);
    98   max_matching.writePos(pos1);
     98  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos1(g);
     99  max_matching.decomposition(pos1);
    99100
    100101  max_matching.run();
    101102  s=0;
    102   max_matching.writeNMapNode(mate);
     103  max_matching.mateMap(mate);
    103104  for(NodeIt v(g); v!=INVALID; ++v) {
    104105    if ( mate[v]!=INVALID ) ++s;
     
    106107  check ( int(s/2) == size, "The size does not equal!" );
    107108 
    108   Graph::NodeMap<MaxMatching<Graph>::pos_enum> pos2(g);
    109   max_matching.writePos(pos2);
     109  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos2(g);
     110  max_matching.decomposition(pos2);
    110111
    111   max_matching.resetMatching();
    112112  max_matching.run();
    113113  s=0;
    114   max_matching.writeNMapNode(mate);
     114  max_matching.mateMap(mate);
    115115  for(NodeIt v(g); v!=INVALID; ++v) {
    116116    if ( mate[v]!=INVALID ) ++s;
     
    118118  check ( int(s/2) == size, "The size does not equal!" );
    119119 
    120   Graph::NodeMap<MaxMatching<Graph>::pos_enum> pos(g);
    121   max_matching.writePos(pos);
     120  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos(g);
     121  max_matching.decomposition(pos);
    122122   
    123123  bool ismatching=true;
     
    140140  bool noedge=true;
    141141  for(UEdgeIt e(g); e!=INVALID; ++e) {
    142    if ( (pos[g.target(e)]==max_matching.C && pos[g.source(e)]==max_matching.D) ||
    143          (pos[g.target(e)]==max_matching.D && pos[g.source(e)]==max_matching.C) )
     142   if ( (pos[g.target(e)]==max_matching.C &&
     143         pos[g.source(e)]==max_matching.D) ||
     144         (pos[g.target(e)]==max_matching.D &&
     145          pos[g.source(e)]==max_matching.C) )
    144146      noedge=false;
    145147  }
  • test/unionfind_test.cc

    r2391 r2505  
    5050  U.insert(2);
    5151
    52   check(U.join(1,2),"Test failed.");
     52  check(U.join(1,2) != -1,"Test failed.");
    5353
    5454  U.insert(3);
     
    5858  U.insert(7);
    5959
    60   check(U.join(1,4),"Test failed.");
    61   check(!U.join(2,4),"Test failed.");
    62   check(U.join(3,5),"Test failed.");
    6360
    64   U.insert(8,5);
     61  check(U.join(1,4) != -1,"Test failed.");
     62  check(U.join(2,4) == -1,"Test failed.");
     63  check(U.join(3,5) != -1,"Test failed.");
    6564
    66   check(U.size(4) == 3,"Test failed.");
    67   check(U.size(5) == 3,"Test failed.");
    68   check(U.size(6) == 1,"Test failed.");
    69   check(U.size(2) == 3,"Test failed.");
     65
     66  U.insert(8,U.find(5));
     67
     68
     69  check(U.size(U.find(4)) == 3,"Test failed.");
     70  check(U.size(U.find(5)) == 3,"Test failed.");
     71  check(U.size(U.find(6)) == 1,"Test failed.");
     72  check(U.size(U.find(2)) == 3,"Test failed.");
     73
    7074
    7175  U.insert(9);
    72   U.insert(10,9);
     76  U.insert(10,U.find(9));
    7377
    74   check(U.join(8,10),"Test failed.");
    7578
    76   check(U.move(9,4),"Test failed.");
    77   check(!U.move(9,2),"Test failed.");
     79  check(U.join(8,10) != -1,"Test failed.");
    7880
    79   check(U.size(4) == 4,"Test failed.");
    80   check(U.size(9) == 4,"Test failed.");
    8181
    82   check(U.move(5,6),"Test failed.");
     82  check(U.size(U.find(4)) == 3,"Test failed.");
     83  check(U.size(U.find(9)) == 5,"Test failed.");
    8384
    84   check(U.size(5) == 2,"Test failed.");
    85   check(U.size(8) == 3,"Test failed.");
    86 
    87   check(U.move(7,10),"Test failed.");
    88   check(U.size(7) == 4,"Test failed.");
     85  check(U.size(U.find(8)) == 5,"Test failed.");
    8986
    9087  U.erase(9);
    9188  U.erase(1);
    9289
    93   check(U.size(4) == 2,"Test failed.");
    94   check(U.size(2) == 2,"Test failed.");
     90  check(U.size(U.find(10)) == 4,"Test failed.");
     91  check(U.size(U.find(2)) == 2,"Test failed.");
    9592
    9693  U.erase(6);
    97   U.split(8);
    98 
    99   check(U.size(4) == 2,"Test failed.");
    100   check(U.size(3) == 1,"Test failed.");
    101   check(U.size(2) == 2,"Test failed.");
    102 
    103   check(U.join(3,4),"Test failed.");
    104   check(!U.join(2,4),"Test failed.");
    105 
    106   check(U.size(4) == 3,"Test failed.");
    107   check(U.size(3) == 3,"Test failed.");
    108   check(U.size(2) == 3,"Test failed.");
    109 
    110   U.makeRep(4);
    111   U.makeRep(3);
    112   U.makeRep(2);
    113 
    114   check(U.size(4) == 3,"Test failed.");
    115   check(U.size(3) == 3,"Test failed.");
    116   check(U.size(2) == 3,"Test failed.");
     94  U.split(U.find(8));
    11795
    11896
    119   U.eraseClass(4);
    120   U.eraseClass(7);
     97  check(U.size(U.find(4)) == 2,"Test failed.");
     98  check(U.size(U.find(3)) == 1,"Test failed.");
     99  check(U.size(U.find(2)) == 2,"Test failed.");
    121100
     101
     102  check(U.join(3,4) != -1,"Test failed.");
     103  check(U.join(2,4) == -1,"Test failed.");
     104
     105
     106  check(U.size(U.find(4)) == 3,"Test failed.");
     107  check(U.size(U.find(3)) == 3,"Test failed.");
     108  check(U.size(U.find(2)) == 3,"Test failed.");
     109
     110  U.eraseClass(U.find(4));
     111  U.eraseClass(U.find(7));
     112
     113  return 0;
    122114}
Note: See TracChangeset for help on using the changeset viewer.