Redesign interface of MaxMatching and UnionFindEnum
authordeba
Tue, 30 Oct 2007 20:21:10 +0000
changeset 25051bb471764ab8
parent 2504 46a82ce84cc6
child 2506 216c6bd5c18c
Redesign interface of MaxMatching and UnionFindEnum
New class ExtendFindEnum

Faster MaxMatching
lemon/max_matching.h
lemon/unionfind.h
test/max_matching_test.cc
test/unionfind_test.cc
     1.1 --- a/lemon/max_matching.h	Tue Oct 30 10:51:07 2007 +0000
     1.2 +++ b/lemon/max_matching.h	Tue Oct 30 20:21:10 2007 +0000
     1.3 @@ -30,25 +30,21 @@
     1.4  
     1.5  namespace lemon {
     1.6  
     1.7 -  /// \ingroup matching
     1.8 +  ///\ingroup matching
     1.9  
    1.10 -  ///Edmonds' alternating forest maximum matching algorithm.
    1.11 -
    1.12 +  ///\brief Edmonds' alternating forest maximum matching algorithm.
    1.13 +  ///
    1.14    ///This class provides Edmonds' alternating forest matching
    1.15    ///algorithm. The starting matching (if any) can be passed to the
    1.16 -  ///algorithm using read-in functions \ref readNMapNode, \ref
    1.17 -  ///readNMapEdge or \ref readEMapBool depending on the container. The
    1.18 -  ///resulting maximum matching can be attained by write-out functions
    1.19 -  ///\ref writeNMapNode, \ref writeNMapEdge or \ref writeEMapBool
    1.20 -  ///depending on the preferred container. 
    1.21 +  ///algorithm using some of init functions.
    1.22    ///
    1.23    ///The dual side of a matching is a map of the nodes to
    1.24 -  ///MaxMatching::pos_enum, having values D, A and C showing the
    1.25 -  ///Gallai-Edmonds decomposition of the graph. The nodes in D induce
    1.26 -  ///a graph with factor-critical components, the nodes in A form the
    1.27 -  ///barrier, and the nodes in C induce a graph having a perfect
    1.28 -  ///matching. This decomposition can be attained by calling \ref
    1.29 -  ///writePos after running the algorithm. 
    1.30 +  ///MaxMatching::DecompType, having values \c D, \c A and \c C
    1.31 +  ///showing the Gallai-Edmonds decomposition of the graph. The nodes
    1.32 +  ///in \c D induce a graph with factor-critical components, the nodes
    1.33 +  ///in \c A form the barrier, and the nodes in \c C induce a graph
    1.34 +  ///having a perfect matching. This decomposition can be attained by
    1.35 +  ///calling \c decomposition() after running the algorithm.
    1.36    ///
    1.37    ///\param Graph The undirected graph type the algorithm runs on.
    1.38    ///
    1.39 @@ -67,17 +63,21 @@
    1.40  
    1.41      typedef typename Graph::template NodeMap<int> UFECrossRef;
    1.42      typedef UnionFindEnum<UFECrossRef> UFE;
    1.43 +    typedef std::vector<Node> NV;
    1.44 +
    1.45 +    typedef typename Graph::template NodeMap<int> EFECrossRef;
    1.46 +    typedef ExtendFindEnum<EFECrossRef> EFE;
    1.47  
    1.48    public:
    1.49      
    1.50 -    ///Indicates the Gallai-Edmonds decomposition of the graph.
    1.51 -
    1.52 +    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
    1.53 +    ///
    1.54      ///Indicates the Gallai-Edmonds decomposition of the graph, which
    1.55      ///shows an upper bound on the size of a maximum matching. The
    1.56 -    ///nodes with pos_enum \c D induce a graph with factor-critical
    1.57 +    ///nodes with DecompType \c D induce a graph with factor-critical
    1.58      ///components, the nodes in \c A form the canonical barrier, and the
    1.59      ///nodes in \c C induce a graph having a perfect matching. 
    1.60 -    enum pos_enum {
    1.61 +    enum DecompType {
    1.62        D=0,
    1.63        A=1,
    1.64        C=2
    1.65 @@ -88,71 +88,32 @@
    1.66      static const int HEUR_density=2;
    1.67      const Graph& g;
    1.68      typename Graph::template NodeMap<Node> _mate;
    1.69 -    typename Graph::template NodeMap<pos_enum> position;
    1.70 +    typename Graph::template NodeMap<DecompType> position;
    1.71       
    1.72    public:
    1.73      
    1.74 -    MaxMatching(const Graph& _g) : g(_g), _mate(_g,INVALID), position(_g) {}
    1.75 +    MaxMatching(const Graph& _g) 
    1.76 +      : g(_g), _mate(_g), position(_g) {}
    1.77  
    1.78 -    ///Runs Edmonds' algorithm.
    1.79 -
    1.80 -    ///Runs Edmonds' algorithm for sparse graphs (number of edges <
    1.81 -    ///2*number of nodes), and a heuristical Edmonds' algorithm with a
    1.82 -    ///heuristic of postponing shrinks for dense graphs. 
    1.83 -    void run() {
    1.84 -      if ( countUEdges(g) < HEUR_density*countNodes(g) ) {
    1.85 -	greedyMatching();
    1.86 -	runEdmonds(0);
    1.87 -      } else runEdmonds(1);
    1.88 -    }
    1.89 -
    1.90 -
    1.91 -    ///Runs Edmonds' algorithm.
    1.92 -    
    1.93 -    ///If heur=0 it runs Edmonds' algorithm. If heur=1 it runs
    1.94 -    ///Edmonds' algorithm with a heuristic of postponing shrinks,
    1.95 -    ///giving a faster algorithm for dense graphs.  
    1.96 -    void runEdmonds( int heur = 1 ) {
    1.97 -      
    1.98 -      //each vertex is put to C
    1.99 -      for(NodeIt v(g); v!=INVALID; ++v)
   1.100 -	position.set(v,C);      
   1.101 -      
   1.102 -      typename Graph::template NodeMap<Node> ear(g,INVALID); 
   1.103 -      //undefined for the base nodes of the blossoms (i.e. for the
   1.104 -      //representative elements of UFE blossom) and for the nodes in C 
   1.105 -      
   1.106 -      UFECrossRef blossom_base(g);
   1.107 -      UFE blossom(blossom_base);
   1.108 -
   1.109 -      UFECrossRef tree_base(g);
   1.110 -      UFE tree(tree_base);
   1.111 -
   1.112 -      //If these UFE's would be members of the class then also
   1.113 -      //blossom_base and tree_base should be a member.
   1.114 -      
   1.115 -      //We build only one tree and the other vertices uncovered by the
   1.116 -      //matching belong to C. (They can be considered as singleton
   1.117 -      //trees.) If this tree can be augmented or no more
   1.118 -      //grow/augmentation/shrink is possible then we return to this
   1.119 -      //"for" cycle.
   1.120 +    ///\brief Sets the actual matching to the empty matching.
   1.121 +    ///
   1.122 +    ///Sets the actual matching to the empty matching.  
   1.123 +    ///
   1.124 +    void init() {
   1.125        for(NodeIt v(g); v!=INVALID; ++v) {
   1.126 -	if ( position[v]==C && _mate[v]==INVALID ) {
   1.127 -	  blossom.insert(v);
   1.128 -	  tree.insert(v); 
   1.129 -	  position.set(v,D);
   1.130 -	  if ( heur == 1 ) lateShrink( v, ear, blossom, tree );
   1.131 -	  else normShrink( v, ear, blossom, tree );
   1.132 -	}
   1.133 +	_mate.set(v,INVALID);      
   1.134 +	position.set(v,C);
   1.135        }
   1.136      }
   1.137  
   1.138 -
   1.139 -    ///Finds a greedy matching starting from the actual matching.
   1.140 -    
   1.141 -    ///Starting form the actual matching stored, it finds a maximal
   1.142 -    ///greedy matching.
   1.143 -    void greedyMatching() {
   1.144 +    ///\brief Finds a greedy matching for initial matching.
   1.145 +    ///
   1.146 +    ///For initial matchig it finds a maximal greedy matching.
   1.147 +    void greedyInit() {
   1.148 +      for(NodeIt v(g); v!=INVALID; ++v) {
   1.149 +	_mate.set(v,INVALID);      
   1.150 +	position.set(v,C);
   1.151 +      }
   1.152        for(NodeIt v(g); v!=INVALID; ++v)
   1.153  	if ( _mate[v]==INVALID ) {
   1.154  	  for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
   1.155 @@ -166,8 +127,153 @@
   1.156  	} 
   1.157      }
   1.158  
   1.159 -    ///Returns the size of the actual matching stored.
   1.160 +    ///\brief Initialize the matching from each nodes' mate.
   1.161 +    ///
   1.162 +    ///Initialize the matching from a \c Node valued \c Node map. This
   1.163 +    ///map must be \e symmetric, i.e. if \c map[u]==v then \c
   1.164 +    ///map[v]==u must hold, and \c uv will be an edge of the initial
   1.165 +    ///matching.
   1.166 +    template <typename MateMap>
   1.167 +    void mateMapInit(MateMap& map) {
   1.168 +      for(NodeIt v(g); v!=INVALID; ++v) {
   1.169 +	_mate.set(v,map[v]);
   1.170 +	position.set(v,C);
   1.171 +      } 
   1.172 +    }
   1.173  
   1.174 +    ///\brief Initialize the matching from a node map with the
   1.175 +    ///incident matching edges.
   1.176 +    ///
   1.177 +    ///Initialize the matching from an \c UEdge valued \c Node map. \c
   1.178 +    ///map[v] must be an \c UEdge incident to \c v. This map must have
   1.179 +    ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
   1.180 +    ///g.oppositeNode(v,map[v])==u holds, and now some edge joining \c
   1.181 +    ///u to \c v will be an edge of the matching.
   1.182 +    template<typename MatchingMap>
   1.183 +    void matchingMapInit(MatchingMap& map) {
   1.184 +      for(NodeIt v(g); v!=INVALID; ++v) {
   1.185 +	position.set(v,C);
   1.186 +	UEdge e=map[v];
   1.187 +	if ( e!=INVALID )
   1.188 +	  _mate.set(v,g.oppositeNode(v,e));
   1.189 +	else 
   1.190 +	  _mate.set(v,INVALID);
   1.191 +      } 
   1.192 +    } 
   1.193 +
   1.194 +    ///\brief Initialize the matching from the map containing the
   1.195 +    ///undirected matching edges.
   1.196 +    ///
   1.197 +    ///Initialize the matching from a \c bool valued \c UEdge map. This
   1.198 +    ///map must have the property that there are no two incident edges
   1.199 +    ///\c e, \c f with \c map[e]==map[f]==true. The edges \c e with \c
   1.200 +    ///map[e]==true form the matching.
   1.201 +    template <typename MatchingMap>
   1.202 +    void matchingInit(MatchingMap& map) {
   1.203 +      for(NodeIt v(g); v!=INVALID; ++v) {
   1.204 +	_mate.set(v,INVALID);      
   1.205 +	position.set(v,C);
   1.206 +      }
   1.207 +      for(UEdgeIt e(g); e!=INVALID; ++e) {
   1.208 +	if ( map[e] ) {
   1.209 +	  Node u=g.source(e);	  
   1.210 +	  Node v=g.target(e);
   1.211 +	  _mate.set(u,v);
   1.212 +	  _mate.set(v,u);
   1.213 +	} 
   1.214 +      } 
   1.215 +    }
   1.216 +
   1.217 +
   1.218 +    ///\brief Runs Edmonds' algorithm.
   1.219 +    ///
   1.220 +    ///Runs Edmonds' algorithm for sparse graphs (number of edges <
   1.221 +    ///2*number of nodes), and a heuristical Edmonds' algorithm with a
   1.222 +    ///heuristic of postponing shrinks for dense graphs. 
   1.223 +    void run() {
   1.224 +      if (countUEdges(g) < HEUR_density * countNodes(g)) {
   1.225 +	greedyInit();
   1.226 +	startSparse();
   1.227 +      } else {
   1.228 +	init();
   1.229 +	startDense();
   1.230 +      }
   1.231 +    }
   1.232 +
   1.233 +
   1.234 +    ///\brief Starts Edmonds' algorithm.
   1.235 +    /// 
   1.236 +    ///If runs the original Edmonds' algorithm.  
   1.237 +    void startSparse() {
   1.238 +            
   1.239 +      typename Graph::template NodeMap<Node> ear(g,INVALID); 
   1.240 +      //undefined for the base nodes of the blossoms (i.e. for the
   1.241 +      //representative elements of UFE blossom) and for the nodes in C 
   1.242 +      
   1.243 +      UFECrossRef blossom_base(g);
   1.244 +      UFE blossom(blossom_base);
   1.245 +      NV rep(countNodes(g));
   1.246 +
   1.247 +      EFECrossRef tree_base(g);
   1.248 +      EFE tree(tree_base);
   1.249 +
   1.250 +      //If these UFE's would be members of the class then also
   1.251 +      //blossom_base and tree_base should be a member.
   1.252 +      
   1.253 +      //We build only one tree and the other vertices uncovered by the
   1.254 +      //matching belong to C. (They can be considered as singleton
   1.255 +      //trees.) If this tree can be augmented or no more
   1.256 +      //grow/augmentation/shrink is possible then we return to this
   1.257 +      //"for" cycle.
   1.258 +      for(NodeIt v(g); v!=INVALID; ++v) {
   1.259 +	if (position[v]==C && _mate[v]==INVALID) {
   1.260 +	  rep[blossom.insert(v)] = v;
   1.261 +	  tree.insert(v); 
   1.262 +	  position.set(v,D);
   1.263 +	  normShrink(v, ear, blossom, rep, tree);
   1.264 +	}
   1.265 +      }
   1.266 +    }
   1.267 +
   1.268 +    ///\brief Starts Edmonds' algorithm.
   1.269 +    /// 
   1.270 +    ///It runs Edmonds' algorithm with a heuristic of postponing
   1.271 +    ///shrinks, giving a faster algorithm for dense graphs.
   1.272 +    void startDense() {
   1.273 +            
   1.274 +      typename Graph::template NodeMap<Node> ear(g,INVALID); 
   1.275 +      //undefined for the base nodes of the blossoms (i.e. for the
   1.276 +      //representative elements of UFE blossom) and for the nodes in C 
   1.277 +      
   1.278 +      UFECrossRef blossom_base(g);
   1.279 +      UFE blossom(blossom_base);
   1.280 +      NV rep(countNodes(g));
   1.281 +
   1.282 +      EFECrossRef tree_base(g);
   1.283 +      EFE tree(tree_base);
   1.284 +
   1.285 +      //If these UFE's would be members of the class then also
   1.286 +      //blossom_base and tree_base should be a member.
   1.287 +      
   1.288 +      //We build only one tree and the other vertices uncovered by the
   1.289 +      //matching belong to C. (They can be considered as singleton
   1.290 +      //trees.) If this tree can be augmented or no more
   1.291 +      //grow/augmentation/shrink is possible then we return to this
   1.292 +      //"for" cycle.
   1.293 +      for(NodeIt v(g); v!=INVALID; ++v) {
   1.294 +	if ( position[v]==C && _mate[v]==INVALID ) {
   1.295 +	  rep[blossom.insert(v)] = v;
   1.296 +	  tree.insert(v); 
   1.297 +	  position.set(v,D);
   1.298 +	  lateShrink(v, ear, blossom, rep, tree);
   1.299 +	}
   1.300 +      }
   1.301 +    }
   1.302 +
   1.303 +
   1.304 +
   1.305 +    ///\brief Returns the size of the actual matching stored.
   1.306 +    ///
   1.307      ///Returns the size of the actual matching stored. After \ref
   1.308      ///run() it returns the size of a maximum matching in the graph.
   1.309      int size() const {
   1.310 @@ -181,75 +287,70 @@
   1.311      }
   1.312  
   1.313  
   1.314 -    ///Resets the actual matching to the empty matching.
   1.315 -
   1.316 -    ///Resets the actual matching to the empty matching.  
   1.317 +    ///\brief Returns the mate of a node in the actual matching.
   1.318      ///
   1.319 -    void resetMatching() {
   1.320 -      for(NodeIt v(g); v!=INVALID; ++v)
   1.321 -	_mate.set(v,INVALID);      
   1.322 -    }
   1.323 -
   1.324 -    ///Returns the mate of a node in the actual matching.
   1.325 -
   1.326      ///Returns the mate of a \c node in the actual matching. 
   1.327      ///Returns INVALID if the \c node is not covered by the actual matching. 
   1.328 -    Node mate(Node& node) const {
   1.329 +    Node mate(const Node& node) const {
   1.330        return _mate[node];
   1.331      } 
   1.332  
   1.333 -    ///Reads a matching from a \c Node valued \c Node map.
   1.334 +    ///\brief Returns the matching edge incident to the given node.
   1.335 +    ///
   1.336 +    ///Returns the matching edge of a \c node in the actual matching. 
   1.337 +    ///Returns INVALID if the \c node is not covered by the actual matching. 
   1.338 +    UEdge matchingEdge(const Node& node) const {
   1.339 +      if (_mate[node] == INVALID) return INVALID;
   1.340 +      Node n = node < _mate[node] ? node : _mate[node];
   1.341 +      for (IncEdgeIt e(g, n); e != INVALID; ++e) {
   1.342 +	if (g.oppositeNode(n, e) == _mate[n]) {
   1.343 +	  return e;
   1.344 +	}
   1.345 +      }
   1.346 +      return INVALID;
   1.347 +    } 
   1.348  
   1.349 -    ///Reads a matching from a \c Node valued \c Node map. This map
   1.350 -    ///must be \e symmetric, i.e. if \c map[u]==v then \c map[v]==u
   1.351 -    ///must hold, and \c uv will be an edge of the matching.
   1.352 -    template<typename NMapN>
   1.353 -    void readNMapNode(NMapN& map) {
   1.354 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.355 -	_mate.set(v,map[v]);   
   1.356 -      } 
   1.357 -    } 
   1.358 +    /// \brief Returns the class of the node in the Edmonds-Gallai
   1.359 +    /// decomposition.
   1.360 +    ///
   1.361 +    /// Returns the class of the node in the Edmonds-Gallai
   1.362 +    /// decomposition.
   1.363 +    DecompType decomposition(const Node& n) {
   1.364 +      return position[n] == A;
   1.365 +    }
   1.366 +
   1.367 +    /// \brief Returns true when the node is in the barrier.
   1.368 +    ///
   1.369 +    /// Returns true when the node is in the barrier.
   1.370 +    bool barrier(const Node& n) {
   1.371 +      return position[n] == A;
   1.372 +    }
   1.373      
   1.374 -    ///Writes the stored matching to a \c Node valued \c Node map.
   1.375 -
   1.376 +    ///\brief Gives back the matching in a \c Node of mates.
   1.377 +    ///
   1.378      ///Writes the stored matching to a \c Node valued \c Node map. The
   1.379      ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
   1.380      ///map[v]==u will hold, and now \c uv is an edge of the matching.
   1.381 -    template<typename NMapN>
   1.382 -    void writeNMapNode (NMapN& map) const {
   1.383 +    template <typename MateMap>
   1.384 +    void mateMap(MateMap& map) const {
   1.385        for(NodeIt v(g); v!=INVALID; ++v) {
   1.386  	map.set(v,_mate[v]);   
   1.387        } 
   1.388      } 
   1.389 -
   1.390 -    ///Reads a matching from an \c UEdge valued \c Node map.
   1.391 -
   1.392 -    ///Reads a matching from an \c UEdge valued \c Node map. \c
   1.393 -    ///map[v] must be an \c UEdge incident to \c v. This map must
   1.394 -    ///have the property that if \c g.oppositeNode(u,map[u])==v then
   1.395 -    ///\c \c g.oppositeNode(v,map[v])==u holds, and now some edge
   1.396 -    ///joining \c u to \c v will be an edge of the matching.
   1.397 -    template<typename NMapE>
   1.398 -    void readNMapEdge(NMapE& map) {
   1.399 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.400 -	UEdge e=map[v];
   1.401 -	if ( e!=INVALID )
   1.402 -	  _mate.set(v,g.oppositeNode(v,e));
   1.403 -      } 
   1.404 -    } 
   1.405      
   1.406 -    ///Writes the matching stored to an \c UEdge valued \c Node map.
   1.407 -
   1.408 +    ///\brief Gives back the matching in an \c UEdge valued \c Node
   1.409 +    ///map.
   1.410 +    ///
   1.411      ///Writes the stored matching to an \c UEdge valued \c Node
   1.412      ///map. \c map[v] will be an \c UEdge incident to \c v. This
   1.413      ///map will have the property that if \c g.oppositeNode(u,map[u])
   1.414      ///== v then \c map[u]==map[v] holds, and now this edge is an edge
   1.415      ///of the matching.
   1.416 -    template<typename NMapE>
   1.417 -    void writeNMapEdge (NMapE& map)  const {
   1.418 +    template <typename MatchingMap>
   1.419 +    void matchingMap(MatchingMap& map)  const {
   1.420        typename Graph::template NodeMap<bool> todo(g,true); 
   1.421        for(NodeIt v(g); v!=INVALID; ++v) {
   1.422 -	if ( todo[v] && _mate[v]!=INVALID ) {
   1.423 +	if (_mate[v]!=INVALID && v < _mate[v]) {
   1.424  	  Node u=_mate[v];
   1.425  	  for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
   1.426  	    if ( g.runningNode(e) == u ) {
   1.427 @@ -265,33 +366,15 @@
   1.428      }
   1.429  
   1.430  
   1.431 -    ///Reads a matching from a \c bool valued \c Edge map.
   1.432 -    
   1.433 -    ///Reads a matching from a \c bool valued \c Edge map. This map
   1.434 -    ///must have the property that there are no two incident edges \c
   1.435 -    ///e, \c f with \c map[e]==map[f]==true. The edges \c e with \c
   1.436 -    ///map[e]==true form the matching.
   1.437 -    template<typename EMapB>
   1.438 -    void readEMapBool(EMapB& map) {
   1.439 -      for(UEdgeIt e(g); e!=INVALID; ++e) {
   1.440 -	if ( map[e] ) {
   1.441 -	  Node u=g.source(e);	  
   1.442 -	  Node v=g.target(e);
   1.443 -	  _mate.set(u,v);
   1.444 -	  _mate.set(v,u);
   1.445 -	} 
   1.446 -      } 
   1.447 -    }
   1.448 -
   1.449 -
   1.450 -    ///Writes the matching stored to a \c bool valued \c Edge map.
   1.451 -
   1.452 +    ///\brief Gives back the matching in a \c bool valued \c UEdge
   1.453 +    ///map.
   1.454 +    ///
   1.455      ///Writes the matching stored to a \c bool valued \c Edge
   1.456      ///map. This map will have the property that there are no two
   1.457      ///incident edges \c e, \c f with \c map[e]==map[f]==true. The
   1.458      ///edges \c e with \c map[e]==true form the matching.
   1.459 -    template<typename EMapB>
   1.460 -    void writeEMapBool (EMapB& map) const {
   1.461 +    template<typename MatchingMap>
   1.462 +    void matching(MatchingMap& map) const {
   1.463        for(UEdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
   1.464  
   1.465        typename Graph::template NodeMap<bool> todo(g,true); 
   1.466 @@ -311,262 +394,230 @@
   1.467      }
   1.468  
   1.469  
   1.470 -    ///Writes the canonical decomposition of the graph after running
   1.471 +    ///\brief Returns the canonical decomposition of the graph after running
   1.472      ///the algorithm.
   1.473 -
   1.474 +    ///
   1.475      ///After calling any run methods of the class, it writes the
   1.476      ///Gallai-Edmonds canonical decomposition of the graph. \c map
   1.477 -    ///must be a node map of \ref pos_enum 's.
   1.478 -    template<typename NMapEnum>
   1.479 -    void writePos (NMapEnum& map) const {
   1.480 -      for(NodeIt v(g); v!=INVALID; ++v)  map.set(v,position[v]);
   1.481 +    ///must be a node map of \ref DecompType 's.
   1.482 +    template <typename DecompositionMap>
   1.483 +    void decomposition(DecompositionMap& map) const {
   1.484 +      for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
   1.485 +    }
   1.486 +
   1.487 +    ///\brief Returns a barrier on the nodes.
   1.488 +    ///
   1.489 +    ///After calling any run methods of the class, it writes a
   1.490 +    ///canonical barrier on the nodes. The odd component number of the
   1.491 +    ///remaining graph minus the barrier size is a lower bound for the
   1.492 +    ///uncovered nodes in the graph. The \c map must be a node map of
   1.493 +    ///bools.
   1.494 +    template <typename BarrierMap>
   1.495 +    void barrier(BarrierMap& barrier) {
   1.496 +      for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);    
   1.497      }
   1.498  
   1.499    private: 
   1.500  
   1.501   
   1.502      void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,  
   1.503 -		    UFE& blossom, UFE& tree);
   1.504 +		    UFE& blossom, NV& rep, EFE& tree) {
   1.505 +      //We have one tree which we grow, and also shrink but only if it
   1.506 +      //cannot be postponed. If we augment then we return to the "for"
   1.507 +      //cycle of runEdmonds().
   1.508 +
   1.509 +      std::queue<Node> Q;   //queue of the totally unscanned nodes
   1.510 +      Q.push(v);  
   1.511 +      std::queue<Node> R;   
   1.512 +      //queue of the nodes which must be scanned for a possible shrink
   1.513 +      
   1.514 +      while ( !Q.empty() ) {
   1.515 +	Node x=Q.front();
   1.516 +	Q.pop();
   1.517 +	for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
   1.518 +	  Node y=g.runningNode(e);
   1.519 +	  //growOrAugment grows if y is covered by the matching and
   1.520 +	  //augments if not. In this latter case it returns 1.
   1.521 +	  if (position[y]==C && 
   1.522 +	      growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
   1.523 +	}
   1.524 +	R.push(x);
   1.525 +      }
   1.526 +      
   1.527 +      while ( !R.empty() ) {
   1.528 +	Node x=R.front();
   1.529 +	R.pop();
   1.530 +	
   1.531 +	for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
   1.532 +	  Node y=g.runningNode(e);
   1.533 +
   1.534 +	  if ( position[y] == D && blossom.find(x) != blossom.find(y) )  
   1.535 +	    //Recall that we have only one tree.
   1.536 +	    shrink( x, y, ear, blossom, rep, tree, Q);	
   1.537 +	
   1.538 +	  while ( !Q.empty() ) {
   1.539 +	    Node z=Q.front();
   1.540 +	    Q.pop();
   1.541 +	    for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
   1.542 +	      Node w=g.runningNode(f);
   1.543 +	      //growOrAugment grows if y is covered by the matching and
   1.544 +	      //augments if not. In this latter case it returns 1.
   1.545 +	      if (position[w]==C && 
   1.546 +		  growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
   1.547 +	    }
   1.548 +	    R.push(z);
   1.549 +	  }
   1.550 +	} //for e
   1.551 +      } // while ( !R.empty() )
   1.552 +    }
   1.553  
   1.554      void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,  
   1.555 -		    UFE& blossom, UFE& tree);
   1.556 +		    UFE& blossom, NV& rep, EFE& tree) {
   1.557 +      //We have one tree, which we grow and shrink. If we augment then we
   1.558 +      //return to the "for" cycle of runEdmonds().
   1.559 +    
   1.560 +      std::queue<Node> Q;   //queue of the unscanned nodes
   1.561 +      Q.push(v);  
   1.562 +      while ( !Q.empty() ) {
   1.563 +
   1.564 +	Node x=Q.front();
   1.565 +	Q.pop();
   1.566 +	
   1.567 +	for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
   1.568 +	  Node y=g.runningNode(e);
   1.569 +	      
   1.570 +	  switch ( position[y] ) {
   1.571 +	  case D:          //x and y must be in the same tree
   1.572 +	    if ( blossom.find(x) != blossom.find(y))
   1.573 +	      //x and y are in the same tree
   1.574 +	      shrink(x, y, ear, blossom, rep, tree, Q);
   1.575 +	    break;
   1.576 +	  case C:
   1.577 +	    //growOrAugment grows if y is covered by the matching and
   1.578 +	    //augments if not. In this latter case it returns 1.
   1.579 +	    if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
   1.580 +	    break;
   1.581 +	  default: break;
   1.582 +	  }
   1.583 +	}
   1.584 +      }
   1.585 +    }
   1.586  
   1.587      void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,  
   1.588 -		UFE& blossom, UFE& tree,std::queue<Node>& Q);
   1.589 +		UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
   1.590 +      //x and y are the two adjacent vertices in two blossoms.
   1.591 +   
   1.592 +      typename Graph::template NodeMap<bool> path(g,false);
   1.593 +    
   1.594 +      Node b=rep[blossom.find(x)];
   1.595 +      path.set(b,true);
   1.596 +      b=_mate[b];
   1.597 +      while ( b!=INVALID ) { 
   1.598 +	b=rep[blossom.find(ear[b])];
   1.599 +	path.set(b,true);
   1.600 +	b=_mate[b];
   1.601 +      } //we go until the root through bases of blossoms and odd vertices
   1.602 +    
   1.603 +      Node top=y;
   1.604 +      Node middle=rep[blossom.find(top)];
   1.605 +      Node bottom=x;
   1.606 +      while ( !path[middle] )
   1.607 +	shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
   1.608 +      //Until we arrive to a node on the path, we update blossom, tree
   1.609 +      //and the positions of the odd nodes.
   1.610 +    
   1.611 +      Node base=middle;
   1.612 +      top=x;
   1.613 +      middle=rep[blossom.find(top)];
   1.614 +      bottom=y;
   1.615 +      Node blossom_base=rep[blossom.find(base)];
   1.616 +      while ( middle!=blossom_base )
   1.617 +	shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
   1.618 +      //Until we arrive to a node on the path, we update blossom, tree
   1.619 +      //and the positions of the odd nodes.
   1.620 +    
   1.621 +      rep[blossom.find(base)] = base;
   1.622 +    }
   1.623  
   1.624      void shrinkStep(Node& top, Node& middle, Node& bottom,
   1.625  		    typename Graph::template NodeMap<Node>& ear,  
   1.626 -		    UFE& blossom, UFE& tree, std::queue<Node>& Q);
   1.627 +		    UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
   1.628 +      //We traverse a blossom and update everything.
   1.629 +    
   1.630 +      ear.set(top,bottom);
   1.631 +      Node t=top;
   1.632 +      while ( t!=middle ) {
   1.633 +	Node u=_mate[t];
   1.634 +	t=ear[u];
   1.635 +	ear.set(t,u);
   1.636 +      } 
   1.637 +      bottom=_mate[middle];
   1.638 +      position.set(bottom,D);
   1.639 +      Q.push(bottom);
   1.640 +      top=ear[bottom];		
   1.641 +      Node oldmiddle=middle;
   1.642 +      middle=rep[blossom.find(top)];
   1.643 +      tree.erase(bottom);
   1.644 +      tree.erase(oldmiddle);
   1.645 +      blossom.insert(bottom);
   1.646 +      blossom.join(bottom, oldmiddle);
   1.647 +      blossom.join(top, oldmiddle);
   1.648 +    }
   1.649 +
   1.650 +
   1.651  
   1.652      bool growOrAugment(Node& y, Node& x, typename Graph::template 
   1.653 -		       NodeMap<Node>& ear, UFE& blossom, UFE& tree, 
   1.654 -		       std::queue<Node>& Q);
   1.655 +		       NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree, 
   1.656 +		       std::queue<Node>& Q) {
   1.657 +      //x is in a blossom in the tree, y is outside. If y is covered by
   1.658 +      //the matching we grow, otherwise we augment. In this case we
   1.659 +      //return 1.
   1.660 +    
   1.661 +      if ( _mate[y]!=INVALID ) {       //grow
   1.662 +	ear.set(y,x);
   1.663 +	Node w=_mate[y];
   1.664 +	rep[blossom.insert(w)] = w;
   1.665 +	position.set(y,A);
   1.666 +	position.set(w,D);
   1.667 +	int t = tree.find(rep[blossom.find(x)]);
   1.668 +	tree.insert(y,t);  
   1.669 +	tree.insert(w,t);  
   1.670 +	Q.push(w);
   1.671 +      } else {                      //augment 
   1.672 +	augment(x, ear, blossom, rep, tree);
   1.673 +	_mate.set(x,y);
   1.674 +	_mate.set(y,x);
   1.675 +	return true;
   1.676 +      }
   1.677 +      return false;
   1.678 +    }
   1.679  
   1.680      void augment(Node x, typename Graph::template NodeMap<Node>& ear,  
   1.681 -		 UFE& blossom, UFE& tree);
   1.682 +		 UFE& blossom, NV& rep, EFE& tree) {
   1.683 +      Node v=_mate[x];
   1.684 +      while ( v!=INVALID ) {
   1.685 +	
   1.686 +	Node u=ear[v];
   1.687 +	_mate.set(v,u);
   1.688 +	Node tmp=v;
   1.689 +	v=_mate[u];
   1.690 +	_mate.set(u,tmp);
   1.691 +      }
   1.692 +      int y = tree.find(rep[blossom.find(x)]);
   1.693 +      for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {   
   1.694 +	if ( position[tit] == D ) {
   1.695 +	  int b = blossom.find(tit);
   1.696 +	  for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) { 
   1.697 +	    position.set(bit, C);
   1.698 +	  }
   1.699 +	  blossom.eraseClass(b);
   1.700 +	} else position.set(tit, C);
   1.701 +      }
   1.702 +      tree.eraseClass(y);
   1.703 +
   1.704 +    }
   1.705  
   1.706    };
   1.707 -
   1.708 -
   1.709 -  // **********************************************************************
   1.710 -  //  IMPLEMENTATIONS
   1.711 -  // **********************************************************************
   1.712 -
   1.713 -
   1.714 -  template <typename Graph>
   1.715 -  void MaxMatching<Graph>::lateShrink(Node v, typename Graph::template 
   1.716 -				      NodeMap<Node>& ear, UFE& blossom, 
   1.717 -                                      UFE& tree) {
   1.718 -    //We have one tree which we grow, and also shrink but only if it cannot be
   1.719 -    //postponed. If we augment then we return to the "for" cycle of
   1.720 -    //runEdmonds().
   1.721 -
   1.722 -    std::queue<Node> Q;   //queue of the totally unscanned nodes
   1.723 -    Q.push(v);  
   1.724 -    std::queue<Node> R;   
   1.725 -    //queue of the nodes which must be scanned for a possible shrink
   1.726 -      
   1.727 -    while ( !Q.empty() ) {
   1.728 -      Node x=Q.front();
   1.729 -      Q.pop();
   1.730 -      for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
   1.731 -	Node y=g.runningNode(e);
   1.732 -	//growOrAugment grows if y is covered by the matching and
   1.733 -	//augments if not. In this latter case it returns 1.
   1.734 -	if ( position[y]==C && growOrAugment(y, x, ear, blossom, tree, Q) ) 
   1.735 -          return;
   1.736 -      }
   1.737 -      R.push(x);
   1.738 -    }
   1.739 -      
   1.740 -    while ( !R.empty() ) {
   1.741 -      Node x=R.front();
   1.742 -      R.pop();
   1.743 -	
   1.744 -      for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
   1.745 -	Node y=g.runningNode(e);
   1.746 -
   1.747 -	if ( position[y] == D && blossom.find(x) != blossom.find(y) )  
   1.748 -	  //Recall that we have only one tree.
   1.749 -	  shrink( x, y, ear, blossom, tree, Q);	
   1.750 -	
   1.751 -	while ( !Q.empty() ) {
   1.752 -	  Node z=Q.front();
   1.753 -	  Q.pop();
   1.754 -	  for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
   1.755 -	    Node w=g.runningNode(f);
   1.756 -	    //growOrAugment grows if y is covered by the matching and
   1.757 -	    //augments if not. In this latter case it returns 1.
   1.758 -	    if ( position[w]==C && growOrAugment(w, z, ear, blossom, tree, Q) )
   1.759 -              return;
   1.760 -	  }
   1.761 -	  R.push(z);
   1.762 -	}
   1.763 -      } //for e
   1.764 -    } // while ( !R.empty() )
   1.765 -  }
   1.766 -
   1.767 -
   1.768 -  template <typename Graph>
   1.769 -  void MaxMatching<Graph>::normShrink(Node v,
   1.770 -				      typename Graph::template
   1.771 -				      NodeMap<Node>& ear,  
   1.772 -				      UFE& blossom, UFE& tree) {
   1.773 -    //We have one tree, which we grow and shrink. If we augment then we
   1.774 -    //return to the "for" cycle of runEdmonds().
   1.775 -    
   1.776 -    std::queue<Node> Q;   //queue of the unscanned nodes
   1.777 -    Q.push(v);  
   1.778 -    while ( !Q.empty() ) {
   1.779 -
   1.780 -      Node x=Q.front();
   1.781 -      Q.pop();
   1.782 -	
   1.783 -      for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
   1.784 -	Node y=g.runningNode(e);
   1.785 -	      
   1.786 -	switch ( position[y] ) {
   1.787 -	case D:          //x and y must be in the same tree
   1.788 -	  if ( blossom.find(x) != blossom.find(y) )
   1.789 -	    //x and y are in the same tree
   1.790 -	    shrink( x, y, ear, blossom, tree, Q);
   1.791 -	  break;
   1.792 -	case C:
   1.793 -	  //growOrAugment grows if y is covered by the matching and
   1.794 -	  //augments if not. In this latter case it returns 1.
   1.795 -	  if ( growOrAugment(y, x, ear, blossom, tree, Q) ) return;
   1.796 -	  break;
   1.797 -	default: break;
   1.798 -       	}
   1.799 -      }
   1.800 -    }
   1.801 -  }
   1.802 -  
   1.803 -
   1.804 -  template <typename Graph>
   1.805 -    void MaxMatching<Graph>::shrink(Node x,Node y, typename 
   1.806 -				    Graph::template NodeMap<Node>& ear,  
   1.807 -				    UFE& blossom, UFE& tree, std::queue<Node>& Q) {
   1.808 -    //x and y are the two adjacent vertices in two blossoms.
   1.809 -   
   1.810 -    typename Graph::template NodeMap<bool> path(g,false);
   1.811 -    
   1.812 -    Node b=blossom.find(x);
   1.813 -    path.set(b,true);
   1.814 -    b=_mate[b];
   1.815 -    while ( b!=INVALID ) { 
   1.816 -      b=blossom.find(ear[b]);
   1.817 -      path.set(b,true);
   1.818 -      b=_mate[b];
   1.819 -    } //we go until the root through bases of blossoms and odd vertices
   1.820 -    
   1.821 -    Node top=y;
   1.822 -    Node middle=blossom.find(top);
   1.823 -    Node bottom=x;
   1.824 -    while ( !path[middle] )
   1.825 -      shrinkStep(top, middle, bottom, ear, blossom, tree, Q);
   1.826 -    //Until we arrive to a node on the path, we update blossom, tree
   1.827 -    //and the positions of the odd nodes.
   1.828 -    
   1.829 -    Node base=middle;
   1.830 -    top=x;
   1.831 -    middle=blossom.find(top);
   1.832 -    bottom=y;
   1.833 -    Node blossom_base=blossom.find(base);
   1.834 -    while ( middle!=blossom_base )
   1.835 -      shrinkStep(top, middle, bottom, ear, blossom, tree, Q);
   1.836 -    //Until we arrive to a node on the path, we update blossom, tree
   1.837 -    //and the positions of the odd nodes.
   1.838 -    
   1.839 -    blossom.makeRep(base);
   1.840 -  }
   1.841 -
   1.842 -
   1.843 -
   1.844 -  template <typename Graph>
   1.845 -  void MaxMatching<Graph>::shrinkStep(Node& top, Node& middle, Node& bottom,
   1.846 -				      typename Graph::template
   1.847 -				      NodeMap<Node>& ear,  
   1.848 -				      UFE& blossom, UFE& tree,
   1.849 -				      std::queue<Node>& Q) {
   1.850 -    //We traverse a blossom and update everything.
   1.851 -    
   1.852 -    ear.set(top,bottom);
   1.853 -    Node t=top;
   1.854 -    while ( t!=middle ) {
   1.855 -      Node u=_mate[t];
   1.856 -      t=ear[u];
   1.857 -      ear.set(t,u);
   1.858 -    } 
   1.859 -    bottom=_mate[middle];
   1.860 -    position.set(bottom,D);
   1.861 -    Q.push(bottom);
   1.862 -    top=ear[bottom];		
   1.863 -    Node oldmiddle=middle;
   1.864 -    middle=blossom.find(top);
   1.865 -    tree.erase(bottom);
   1.866 -    tree.erase(oldmiddle);
   1.867 -    blossom.insert(bottom);
   1.868 -    blossom.join(bottom, oldmiddle);
   1.869 -    blossom.join(top, oldmiddle);
   1.870 -  }
   1.871 -
   1.872 -
   1.873 -  template <typename Graph>
   1.874 -  bool MaxMatching<Graph>::growOrAugment(Node& y, Node& x, typename Graph::template
   1.875 -					 NodeMap<Node>& ear, UFE& blossom, UFE& tree,
   1.876 -					 std::queue<Node>& Q) {
   1.877 -    //x is in a blossom in the tree, y is outside. If y is covered by
   1.878 -    //the matching we grow, otherwise we augment. In this case we
   1.879 -    //return 1.
   1.880 -    
   1.881 -    if ( _mate[y]!=INVALID ) {       //grow
   1.882 -      ear.set(y,x);
   1.883 -      Node w=_mate[y];
   1.884 -      blossom.insert(w);
   1.885 -      position.set(y,A);
   1.886 -      position.set(w,D);
   1.887 -      tree.insert(y);
   1.888 -      tree.insert(w);
   1.889 -      tree.join(y,blossom.find(x));  
   1.890 -      tree.join(w,y);  
   1.891 -      Q.push(w);
   1.892 -    } else {                      //augment 
   1.893 -      augment(x, ear, blossom, tree);
   1.894 -      _mate.set(x,y);
   1.895 -      _mate.set(y,x);
   1.896 -      return true;
   1.897 -    }
   1.898 -    return false;
   1.899 -  }
   1.900 -  
   1.901 -
   1.902 -  template <typename Graph>
   1.903 -  void MaxMatching<Graph>::augment(Node x,
   1.904 -				   typename Graph::template NodeMap<Node>& ear,  
   1.905 -				   UFE& blossom, UFE& tree) { 
   1.906 -    Node v=_mate[x];
   1.907 -    while ( v!=INVALID ) {
   1.908 -	
   1.909 -      Node u=ear[v];
   1.910 -      _mate.set(v,u);
   1.911 -      Node tmp=v;
   1.912 -      v=_mate[u];
   1.913 -      _mate.set(u,tmp);
   1.914 -    }
   1.915 -    Node y=blossom.find(x);
   1.916 -    for (typename UFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {   
   1.917 -      if ( position[tit] == D ) {
   1.918 -	for (typename UFE::ItemIt bit(blossom, tit); bit != INVALID; ++bit) {  
   1.919 -	  position.set( bit ,C);
   1.920 -	}
   1.921 -	blossom.eraseClass(tit);
   1.922 -      } else position.set( tit ,C);
   1.923 -    }
   1.924 -    tree.eraseClass(y);
   1.925 -
   1.926 -  }
   1.927 -
   1.928    
   1.929  } //END OF NAMESPACE LEMON
   1.930  
     2.1 --- a/lemon/unionfind.h	Tue Oct 30 10:51:07 2007 +0000
     2.2 +++ b/lemon/unionfind.h	Tue Oct 30 20:21:10 2007 +0000
     2.3 @@ -178,26 +178,56 @@
     2.4  
     2.5    private:
     2.6      
     2.7 +    ItemIntMap& index;
     2.8 +
     2.9      // If the parent stores negative value for an item then that item
    2.10 -    // is root item and it has -items[it].parent component size.  Else
    2.11 +    // is root item and it has ~(items[it].parent) component id.  Else
    2.12      // the items[it].parent contains the index of the parent.
    2.13      //
    2.14 -    // The \c nextItem and \c prevItem provides the double-linked
    2.15 -    // cyclic list of one component's items. The \c prevClass and
    2.16 -    // \c nextClass gives the double linked list of the representant
    2.17 -    // items. 
    2.18 +    // The \c next and \c prev provides the double-linked
    2.19 +    // cyclic list of one component's items.
    2.20      struct ItemT {
    2.21        int parent;
    2.22        Item item;
    2.23  
    2.24 -      int nextItem, prevItem;
    2.25 -      int nextClass, prevClass;
    2.26 +      int next, prev;
    2.27      };
    2.28  
    2.29      std::vector<ItemT> items;
    2.30 -    ItemIntMap& index;
    2.31 +    int firstFreeItem;
    2.32  
    2.33 -    int firstClass;
    2.34 +    struct ClassT {
    2.35 +      int size;
    2.36 +      int firstItem;
    2.37 +      int next, prev;
    2.38 +    };
    2.39 +    
    2.40 +    std::vector<ClassT> classes;
    2.41 +    int firstClass, firstFreeClass;
    2.42 +
    2.43 +    int newClass() {
    2.44 +      if (firstFreeClass == -1) {
    2.45 +	int cdx = classes.size();
    2.46 +	classes.push_back(ClassT());
    2.47 +	return cdx;
    2.48 +      } else {
    2.49 +	int cdx = firstFreeClass;
    2.50 +	firstFreeClass = classes[firstFreeClass].next;
    2.51 +	return cdx;
    2.52 +      }
    2.53 +    }
    2.54 +
    2.55 +    int newItem() {
    2.56 +      if (firstFreeItem == -1) {
    2.57 +	int idx = items.size();
    2.58 +	items.push_back(ItemT());
    2.59 +	return idx;
    2.60 +      } else {
    2.61 +	int idx = firstFreeItem;
    2.62 +	firstFreeItem = items[firstFreeItem].next;
    2.63 +	return idx;
    2.64 +      }
    2.65 +    }
    2.66  
    2.67  
    2.68      bool rep(int idx) const {
    2.69 @@ -217,79 +247,110 @@
    2.70        return k;
    2.71      }
    2.72  
    2.73 -    void unlaceClass(int k) {
    2.74 -      if (items[k].prevClass != -1) {
    2.75 -        items[items[k].prevClass].nextClass = items[k].nextClass;
    2.76 -      } else {
    2.77 -        firstClass = items[k].nextClass;
    2.78 +    int classIndex(int idx) const {
    2.79 +      return ~(items[repIndex(idx)].parent);
    2.80 +    }
    2.81 +
    2.82 +    void singletonItem(int idx) {
    2.83 +      items[idx].next = idx;
    2.84 +      items[idx].prev = idx;
    2.85 +    }
    2.86 +
    2.87 +    void laceItem(int idx, int rdx) {
    2.88 +      items[idx].prev = rdx;
    2.89 +      items[idx].next = items[rdx].next;
    2.90 +      items[items[rdx].next].prev = idx;
    2.91 +      items[rdx].next = idx;
    2.92 +    }
    2.93 +
    2.94 +    void unlaceItem(int idx) {
    2.95 +      items[items[idx].prev].next = items[idx].next;
    2.96 +      items[items[idx].next].prev = items[idx].prev;
    2.97 +      
    2.98 +      items[idx].next = firstFreeItem;
    2.99 +      firstFreeItem = idx;
   2.100 +    }
   2.101 +
   2.102 +    void spliceItems(int ak, int bk) {
   2.103 +      items[items[ak].prev].next = bk;
   2.104 +      items[items[bk].prev].next = ak;
   2.105 +      int tmp = items[ak].prev;
   2.106 +      items[ak].prev = items[bk].prev;
   2.107 +      items[bk].prev = tmp;
   2.108 +        
   2.109 +    }
   2.110 +
   2.111 +    void laceClass(int cls) {
   2.112 +      if (firstClass != -1) {
   2.113 +        classes[firstClass].prev = cls;
   2.114        }
   2.115 -      if (items[k].nextClass != -1) {
   2.116 -        items[items[k].nextClass].prevClass = items[k].prevClass;
   2.117 -      }
   2.118 +      classes[cls].next = firstClass;
   2.119 +      classes[cls].prev = -1;
   2.120 +      firstClass = cls;
   2.121      } 
   2.122  
   2.123 -    void spliceItems(int ak, int bk) {
   2.124 -      items[items[ak].prevItem].nextItem = bk;
   2.125 -      items[items[bk].prevItem].nextItem = ak;
   2.126 -      int tmp = items[ak].prevItem;
   2.127 -      items[ak].prevItem = items[bk].prevItem;
   2.128 -      items[bk].prevItem = tmp;
   2.129 -        
   2.130 -    }
   2.131 +    void unlaceClass(int cls) {
   2.132 +      if (classes[cls].prev != -1) {
   2.133 +        classes[classes[cls].prev].next = classes[cls].next;
   2.134 +      } else {
   2.135 +        firstClass = classes[cls].next;
   2.136 +      }
   2.137 +      if (classes[cls].next != -1) {
   2.138 +        classes[classes[cls].next].prev = classes[cls].prev;
   2.139 +      }
   2.140 +      
   2.141 +      classes[cls].next = firstFreeClass;
   2.142 +      firstFreeClass = cls;
   2.143 +    } 
   2.144  
   2.145    public:
   2.146  
   2.147      UnionFindEnum(ItemIntMap& _index) 
   2.148 -      : items(), index(_index), firstClass(-1) {}
   2.149 +      : index(_index), items(), firstFreeItem(-1), 
   2.150 +	firstClass(-1), firstFreeClass(-1) {}
   2.151      
   2.152      /// \brief Inserts the given element into a new component.
   2.153      ///
   2.154      /// This method creates a new component consisting only of the
   2.155      /// given element.
   2.156      ///
   2.157 -    void insert(const Item& item) {
   2.158 -      ItemT t;
   2.159 +    int insert(const Item& item) {
   2.160 +      int idx = newItem();
   2.161  
   2.162 -      int idx = items.size();
   2.163        index.set(item, idx);
   2.164  
   2.165 -      t.nextItem = idx;
   2.166 -      t.prevItem = idx;
   2.167 -      t.item = item;
   2.168 -      t.parent = -1;
   2.169 +      singletonItem(idx);
   2.170 +      items[idx].item = item;
   2.171 +
   2.172 +      int cdx = newClass();
   2.173 +
   2.174 +      items[idx].parent = ~cdx;
   2.175 +
   2.176 +      laceClass(cdx);
   2.177 +      classes[cdx].size = 1;
   2.178 +      classes[cdx].firstItem = idx;
   2.179 +
   2.180 +      firstClass = cdx;
   2.181        
   2.182 -      t.nextClass = firstClass;
   2.183 -      if (firstClass != -1) {
   2.184 -        items[firstClass].prevClass = idx;
   2.185 -      }
   2.186 -      t.prevClass = -1;
   2.187 -      firstClass = idx;
   2.188 -
   2.189 -      items.push_back(t);
   2.190 +      return cdx;
   2.191      }
   2.192  
   2.193      /// \brief Inserts the given element into the component of the others.
   2.194      ///
   2.195      /// This methods inserts the element \e a into the component of the
   2.196      /// element \e comp. 
   2.197 -    void insert(const Item& item, const Item& comp) {
   2.198 -      int k = repIndex(index[comp]);
   2.199 -      ItemT t;
   2.200 +    void insert(const Item& item, int cls) {
   2.201 +      int rdx = classes[cls].firstItem;
   2.202 +      int idx = newItem();
   2.203  
   2.204 -      int idx = items.size();
   2.205        index.set(item, idx);
   2.206  
   2.207 -      t.prevItem = k;
   2.208 -      t.nextItem = items[k].nextItem;
   2.209 -      items[items[k].nextItem].prevItem = idx;
   2.210 -      items[k].nextItem = idx;
   2.211 +      laceItem(idx, rdx);
   2.212  
   2.213 -      t.item = item;
   2.214 -      t.parent = k;
   2.215 +      items[idx].item = item;
   2.216 +      items[idx].parent = rdx;
   2.217  
   2.218 -      --items[k].parent;
   2.219 -
   2.220 -      items.push_back(t);
   2.221 +      ++classes[~(items[rdx].parent)].size;
   2.222      }
   2.223  
   2.224      /// \brief Clears the union-find data structure
   2.225 @@ -298,13 +359,14 @@
   2.226      void clear() {
   2.227        items.clear();
   2.228        firstClass = -1;
   2.229 +      firstFreeItem = -1;
   2.230      }
   2.231  
   2.232 -    /// \brief Finds the leader of the component of the given element.
   2.233 +    /// \brief Finds the component of the given element.
   2.234      ///
   2.235 -    /// The method returns the leader of the component of the given element.
   2.236 -    const Item& find(const Item &item) const {
   2.237 -      return items[repIndex(index[item])].item;
   2.238 +    /// The method returns the component id of the given element.
   2.239 +    int find(const Item &item) const {
   2.240 +      return ~(items[repIndex(index[item])].parent);
   2.241      }
   2.242  
   2.243      /// \brief Joining the component of element \e a and element \e b.
   2.244 @@ -312,90 +374,71 @@
   2.245      /// This is the \e union operation of the Union-Find structure. 
   2.246      /// Joins the component of element \e a and component of
   2.247      /// element \e b. If \e a and \e b are in the same component then
   2.248 -    /// returns false else returns true.
   2.249 -    bool join(const Item& a, const Item& b) {
   2.250 +    /// returns -1 else returns the remaining class.
   2.251 +    int join(const Item& a, const Item& b) {
   2.252  
   2.253        int ak = repIndex(index[a]);
   2.254        int bk = repIndex(index[b]);
   2.255  
   2.256        if (ak == bk) {
   2.257 -	return false;
   2.258 +	return -1;
   2.259        }
   2.260  
   2.261 -      if ( items[ak].parent < items[bk].parent ) {
   2.262 -        unlaceClass(bk);
   2.263 -        items[ak].parent += items[bk].parent;
   2.264 +      int acx = ~(items[ak].parent);
   2.265 +      int bcx = ~(items[bk].parent);
   2.266 +
   2.267 +      int rcx;
   2.268 +
   2.269 +      if (classes[acx].size > classes[bcx].size) {
   2.270 +	classes[acx].size += classes[bcx].size;
   2.271  	items[bk].parent = ak;
   2.272 +        unlaceClass(bcx);
   2.273 +	rcx = acx;
   2.274        } else {
   2.275 -        unlaceClass(ak);
   2.276 -        items[bk].parent += items[ak].parent;
   2.277 +	classes[bcx].size += classes[acx].size;
   2.278  	items[ak].parent = bk;
   2.279 +        unlaceClass(acx);
   2.280 +	rcx = bcx;
   2.281        }
   2.282        spliceItems(ak, bk);
   2.283  
   2.284 -      return true;
   2.285 +      return rcx;
   2.286      }
   2.287  
   2.288 -    /// \brief Returns the size of the component of element \e a.
   2.289 +    /// \brief Returns the size of the class.
   2.290      ///
   2.291 -    /// Returns the size of the component of element \e a.
   2.292 -    int size(const Item &item) const {
   2.293 -      return - items[repIndex(index[item])].parent;
   2.294 +    /// Returns the size of the class.
   2.295 +    int size(int cls) const {
   2.296 +      return classes[cls].size;
   2.297      }
   2.298  
   2.299 -    /// \brief Splits up the component of the element. 
   2.300 +    /// \brief Splits up the component. 
   2.301      ///
   2.302 -    /// Splitting the component of the element into sigleton
   2.303 -    /// components (component of size one).
   2.304 -    void split(const Item &item) {
   2.305 -      int k = repIndex(index[item]);
   2.306 -      int idx = items[k].nextItem;
   2.307 -      while (idx != k) {
   2.308 -        int next = items[idx].nextItem;
   2.309 +    /// Splitting the component into singleton components (component
   2.310 +    /// of size one).
   2.311 +    void split(int cls) {
   2.312 +      int fdx = classes[cls].firstItem;
   2.313 +      int idx = items[fdx].next;
   2.314 +      while (idx != fdx) {
   2.315 +        int next = items[idx].next;
   2.316 +
   2.317 +	singletonItem(idx);
   2.318 +
   2.319 +	int cdx = newClass();        
   2.320 +        items[idx].parent = ~cdx;
   2.321 +
   2.322 +	laceClass(cdx);
   2.323 +	classes[cdx].size = 1;
   2.324 +	classes[cdx].firstItem = idx;
   2.325          
   2.326 -        items[idx].parent = -1;
   2.327 -        items[idx].prevItem = idx;
   2.328 -        items[idx].nextItem = idx;
   2.329 -        
   2.330 -        items[idx].nextClass = firstClass;
   2.331 -        items[firstClass].prevClass = idx;
   2.332 -        firstClass = idx;
   2.333 -
   2.334          idx = next;
   2.335        }
   2.336  
   2.337 -      items[idx].parent = -1;
   2.338 -      items[idx].prevItem = idx;
   2.339 -      items[idx].nextItem = idx;
   2.340 +      items[idx].prev = idx;
   2.341 +      items[idx].next = idx;
   2.342 +
   2.343 +      classes[~(items[idx].parent)].size = 1;
   2.344        
   2.345 -      items[firstClass].prevClass = -1;
   2.346 -    }
   2.347 -
   2.348 -    /// \brief Sets the given element to the leader element of its component.
   2.349 -    ///
   2.350 -    /// Sets the given element to the leader element of its component.
   2.351 -    void makeRep(const Item &item) {
   2.352 -      int nk = index[item];
   2.353 -      int k = repIndex(nk);
   2.354 -      if (nk == k) return;
   2.355 -      
   2.356 -      if (items[k].prevClass != -1) {
   2.357 -        items[items[k].prevClass].nextClass = nk;
   2.358 -      } else {
   2.359 -        firstClass = nk;
   2.360 -      }
   2.361 -      if (items[k].nextClass != -1) {
   2.362 -        items[items[k].nextClass].prevClass = nk;
   2.363 -      }
   2.364 -      
   2.365 -      int idx = items[k].nextItem;
   2.366 -      while (idx != k) {
   2.367 -        items[idx].parent = nk;
   2.368 -        idx = items[idx].nextItem;
   2.369 -      }
   2.370 -      
   2.371 -      items[nk].parent = items[k].parent;
   2.372 -      items[k].parent = nk;
   2.373      }
   2.374  
   2.375      /// \brief Removes the given element from the structure.
   2.376 @@ -405,124 +448,97 @@
   2.377      ///
   2.378      /// \warning It is an error to remove an element which is not in
   2.379      /// the structure.
   2.380 -    void erase(const Item &item) {
   2.381 +    /// \warning This running time of this operation is proportional to the
   2.382 +    /// number of the items in this class.
   2.383 +    void erase(const Item& item) {
   2.384        int idx = index[item];
   2.385 -      if (rep(idx)) {
   2.386 -        int k = idx;
   2.387 -        if (items[k].parent == -1) {
   2.388 -          unlaceClass(idx);
   2.389 -          return;
   2.390 -        } else {
   2.391 -          int nk = items[k].nextItem;
   2.392 -          if (items[k].prevClass != -1) {
   2.393 -            items[items[k].prevClass].nextClass = nk;
   2.394 -          } else {
   2.395 -            firstClass = nk;
   2.396 -          }
   2.397 -          if (items[k].nextClass != -1) {
   2.398 -            items[items[k].nextClass].prevClass = nk;
   2.399 -          }
   2.400 -      
   2.401 -          int l = items[k].nextItem;
   2.402 -          while (l != k) {
   2.403 -            items[l].parent = nk;
   2.404 -            l = items[l].nextItem;
   2.405 -          }
   2.406 +      int fdx = items[idx].next;
   2.407 +
   2.408 +      int cdx = classIndex(idx);
   2.409 +      if (idx == fdx) {
   2.410 +	unlaceClass(cdx);
   2.411 +	items[idx].next = firstFreeItem;
   2.412 +	firstFreeItem = idx;
   2.413 +	return;
   2.414 +      } else {
   2.415 +	classes[cdx].firstItem = fdx;
   2.416 +	--classes[cdx].size;
   2.417 +	items[fdx].parent = ~cdx;
   2.418 +
   2.419 +	unlaceItem(idx);
   2.420 +	idx = items[fdx].next;
   2.421 +	while (idx != fdx) {
   2.422 +	  items[idx].parent = fdx;
   2.423 +	  idx = items[idx].next;
   2.424 +	}
   2.425            
   2.426 -          items[nk].parent = items[k].parent + 1;
   2.427 -        }
   2.428 -      } else {
   2.429 -        
   2.430 -        int k = repIndex(idx);
   2.431 -        idx = items[k].nextItem;
   2.432 -        while (idx != k) {
   2.433 -          items[idx].parent = k;
   2.434 -          idx = items[idx].nextItem;
   2.435 -        }
   2.436 -
   2.437 -        ++items[k].parent;
   2.438        }
   2.439  
   2.440 -      idx = index[item];      
   2.441 -      items[items[idx].prevItem].nextItem = items[idx].nextItem;
   2.442 -      items[items[idx].nextItem].prevItem = items[idx].prevItem;
   2.443 -
   2.444      }    
   2.445  
   2.446 -    /// \brief Moves the given element to another component.
   2.447 -    ///
   2.448 -    /// This method moves the element \e a from its component
   2.449 -    /// to the component of \e comp.
   2.450 -    /// If \e a and \e comp are in the same component then
   2.451 -    /// it returns false otherwise it returns true.
   2.452 -    bool move(const Item &item, const Item &comp) {
   2.453 -      if (repIndex(index[item]) == repIndex(index[comp])) return false;
   2.454 -      erase(item);
   2.455 -      insert(item, comp);
   2.456 -      return true;
   2.457 -    }
   2.458 -
   2.459 -
   2.460      /// \brief Removes the component of the given element from the structure.
   2.461      ///
   2.462      /// Removes the component of the given element from the structure.
   2.463      ///
   2.464      /// \warning It is an error to give an element which is not in the
   2.465      /// structure.
   2.466 -    void eraseClass(const Item &item) {
   2.467 -      unlaceClass(repIndex(index[item]));
   2.468 +    void eraseClass(int cls) {
   2.469 +      int fdx = classes[cls].firstItem;
   2.470 +      unlaceClass(cls);
   2.471 +      items[items[fdx].prev].next = firstFreeItem;
   2.472 +      firstFreeItem = fdx;
   2.473      }
   2.474  
   2.475      /// \brief Lemon style iterator for the representant items.
   2.476      ///
   2.477      /// ClassIt is a lemon style iterator for the components. It iterates
   2.478 -    /// on the representant items of the classes.
   2.479 +    /// on the ids of the classes.
   2.480      class ClassIt {
   2.481      public:
   2.482        /// \brief Constructor of the iterator
   2.483        ///
   2.484        /// Constructor of the iterator
   2.485        ClassIt(const UnionFindEnum& ufe) : unionFind(&ufe) {
   2.486 -        idx = unionFind->firstClass;
   2.487 +        cdx = unionFind->firstClass;
   2.488        }
   2.489  
   2.490        /// \brief Constructor to get invalid iterator
   2.491        ///
   2.492        /// Constructor to get invalid iterator
   2.493 -      ClassIt(Invalid) : unionFind(0), idx(-1) {}
   2.494 +      ClassIt(Invalid) : unionFind(0), cdx(-1) {}
   2.495        
   2.496        /// \brief Increment operator
   2.497        ///
   2.498        /// It steps to the next representant item.
   2.499        ClassIt& operator++() {
   2.500 -        idx = unionFind->items[idx].nextClass;
   2.501 +        cdx = unionFind->classes[cdx].next;
   2.502          return *this;
   2.503        }
   2.504        
   2.505        /// \brief Conversion operator
   2.506        ///
   2.507        /// It converts the iterator to the current representant item.
   2.508 -      operator const Item&() const {
   2.509 -        return unionFind->items[idx].item;
   2.510 +      operator int() const {
   2.511 +        return cdx;
   2.512        }
   2.513  
   2.514        /// \brief Equality operator
   2.515        ///
   2.516        /// Equality operator
   2.517        bool operator==(const ClassIt& i) { 
   2.518 -        return i.idx == idx;
   2.519 +        return i.cdx == cdx;
   2.520        }
   2.521  
   2.522        /// \brief Inequality operator
   2.523        ///
   2.524        /// Inequality operator
   2.525        bool operator!=(const ClassIt& i) { 
   2.526 -        return i.idx != idx;
   2.527 +        return i.cdx != cdx;
   2.528        }
   2.529        
   2.530      private:
   2.531        const UnionFindEnum* unionFind;
   2.532 -      int idx;
   2.533 +      int cdx;
   2.534      };
   2.535  
   2.536      /// \brief Lemon style iterator for the items of a component.
   2.537 @@ -545,8 +561,8 @@
   2.538        ///
   2.539        /// Constructor of the iterator. The iterator iterates
   2.540        /// on the class of the \c item.
   2.541 -      ItemIt(const UnionFindEnum& ufe, const Item& item) : unionFind(&ufe) {
   2.542 -        idx = unionFind->repIndex(unionFind->index[item]);
   2.543 +      ItemIt(const UnionFindEnum& ufe, int cls) : unionFind(&ufe) {
   2.544 +        fdx = idx = unionFind->classes[cls].firstItem;
   2.545        }
   2.546  
   2.547        /// \brief Constructor to get invalid iterator
   2.548 @@ -558,8 +574,8 @@
   2.549        ///
   2.550        /// It steps to the next item in the class.
   2.551        ItemIt& operator++() {
   2.552 -        idx = unionFind->items[idx].nextItem;
   2.553 -        if (unionFind->rep(idx)) idx = -1;
   2.554 +        idx = unionFind->items[idx].next;
   2.555 +        if (idx == fdx) idx = -1;
   2.556          return *this;
   2.557        }
   2.558        
   2.559 @@ -586,11 +602,310 @@
   2.560        
   2.561      private:
   2.562        const UnionFindEnum* unionFind;
   2.563 -      int idx;
   2.564 +      int idx, fdx;
   2.565      };
   2.566  
   2.567    };
   2.568  
   2.569 +  /// \ingroup auxdat
   2.570 +  ///
   2.571 +  /// \brief A \e Extend-Find data structure implementation which
   2.572 +  /// is able to enumerate the components.
   2.573 +  ///
   2.574 +  /// The class implements an \e Extend-Find data structure which is
   2.575 +  /// able to enumerate the components and the items in a
   2.576 +  /// component. The data structure is a simplification of the
   2.577 +  /// Union-Find structure, and it does not allow to merge two components.
   2.578 +  ///
   2.579 +  /// \pre You need to add all the elements by the \ref insert()
   2.580 +  /// method.
   2.581 +  template <typename _ItemIntMap>
   2.582 +  class ExtendFindEnum {
   2.583 +  public:
   2.584 +    
   2.585 +    typedef _ItemIntMap ItemIntMap;
   2.586 +    typedef typename ItemIntMap::Key Item;
   2.587 +
   2.588 +  private:
   2.589 +    
   2.590 +    ItemIntMap& index;
   2.591 +
   2.592 +    struct ItemT {
   2.593 +      int cls;
   2.594 +      Item item;
   2.595 +      int next, prev;
   2.596 +    };
   2.597 +
   2.598 +    std::vector<ItemT> items;
   2.599 +    int firstFreeItem;
   2.600 +
   2.601 +    struct ClassT {
   2.602 +      int firstItem;
   2.603 +      int next, prev;
   2.604 +    };
   2.605 +
   2.606 +    std::vector<ClassT> classes;
   2.607 +
   2.608 +    int firstClass, firstFreeClass;
   2.609 +
   2.610 +    int newClass() {
   2.611 +      if (firstFreeClass != -1) {
   2.612 +	int cdx = firstFreeClass;
   2.613 +	firstFreeClass = classes[cdx].next;
   2.614 +	return cdx;
   2.615 +      } else {
   2.616 +	classes.push_back(ClassT());
   2.617 +	return classes.size() - 1;
   2.618 +      }
   2.619 +    }
   2.620 +
   2.621 +    int newItem() {
   2.622 +      if (firstFreeItem != -1) {
   2.623 +	int idx = firstFreeItem;
   2.624 +	firstFreeItem = items[idx].next;
   2.625 +	return idx;
   2.626 +      } else {
   2.627 +	items.push_back(ItemT());
   2.628 +	return items.size() - 1;
   2.629 +      }
   2.630 +    }
   2.631 +
   2.632 +  public:
   2.633 +
   2.634 +    /// \brief Constructor
   2.635 +    ExtendFindEnum(ItemIntMap& _index) 
   2.636 +      : index(_index), items(), firstFreeItem(-1), 
   2.637 +	classes(), firstClass(-1), firstFreeClass(-1) {}
   2.638 +    
   2.639 +    /// \brief Inserts the given element into a new component.
   2.640 +    ///
   2.641 +    /// This method creates a new component consisting only of the
   2.642 +    /// given element.
   2.643 +    int insert(const Item& item) {
   2.644 +      int cdx = newClass();
   2.645 +      classes[cdx].prev = -1;
   2.646 +      classes[cdx].next = firstClass;
   2.647 +      firstClass = cdx;
   2.648 +      
   2.649 +      int idx = newItem();
   2.650 +      items[idx].item = item;
   2.651 +      items[idx].cls = cdx;
   2.652 +      items[idx].prev = idx;
   2.653 +      items[idx].next = idx;
   2.654 +
   2.655 +      classes[cdx].firstItem = idx;
   2.656 +
   2.657 +      index.set(item, idx);
   2.658 +      
   2.659 +      return cdx;
   2.660 +    }
   2.661 +
   2.662 +    /// \brief Inserts the given element into the given component.
   2.663 +    ///
   2.664 +    /// This methods inserts the element \e item a into the \e cls class.
   2.665 +    void insert(const Item& item, int cls) {
   2.666 +      int idx = newItem();
   2.667 +      int rdx = classes[cls].firstItem;
   2.668 +      items[idx].item = item;
   2.669 +      items[idx].cls = cls;
   2.670 +
   2.671 +      items[idx].prev = rdx;
   2.672 +      items[idx].next = items[rdx].next;
   2.673 +      items[items[rdx].next].prev = idx;
   2.674 +      items[rdx].next = idx;
   2.675 +
   2.676 +      index.set(item, idx);
   2.677 +    }
   2.678 +
   2.679 +    /// \brief Clears the union-find data structure
   2.680 +    ///
   2.681 +    /// Erase each item from the data structure.
   2.682 +    void clear() {
   2.683 +      items.clear();
   2.684 +      classes.clear;
   2.685 +      firstClass = firstFreeClass = firstFreeItem = -1;
   2.686 +    }
   2.687 +
   2.688 +    /// \brief Gives back the class of the \e item.
   2.689 +    ///
   2.690 +    /// Gives back the class of the \e item.
   2.691 +    int find(const Item &item) const {
   2.692 +      return items[index[item]].cls;
   2.693 +    }
   2.694 +    
   2.695 +    /// \brief Removes the given element from the structure.
   2.696 +    ///
   2.697 +    /// Removes the element from its component and if the component becomes
   2.698 +    /// empty then removes that component from the component list.
   2.699 +    ///
   2.700 +    /// \warning It is an error to remove an element which is not in
   2.701 +    /// the structure.
   2.702 +    void erase(const Item &item) {
   2.703 +      int idx = index[item];
   2.704 +      int cdx = items[idx].cls;
   2.705 +      
   2.706 +      if (idx == items[idx].next) {
   2.707 +	if (classes[cdx].prev != -1) {
   2.708 +	  classes[classes[cdx].prev].next = classes[cdx].next; 
   2.709 +	} else {
   2.710 +	  firstClass = classes[cdx].next;
   2.711 +	}
   2.712 +	if (classes[cdx].next != -1) {
   2.713 +	  classes[classes[cdx].next].prev = classes[cdx].prev; 
   2.714 +	}
   2.715 +	classes[cdx].next = firstFreeClass;
   2.716 +	firstFreeClass = cdx;
   2.717 +      } else {
   2.718 +	classes[cdx].firstItem = items[idx].next;
   2.719 +	items[items[idx].next].prev = items[idx].prev;
   2.720 +	items[items[idx].prev].next = items[idx].next;
   2.721 +      }
   2.722 +      items[idx].next = firstFreeItem;
   2.723 +      firstFreeItem = idx;
   2.724 +	
   2.725 +    }    
   2.726 +
   2.727 +    
   2.728 +    /// \brief Removes the component of the given element from the structure.
   2.729 +    ///
   2.730 +    /// Removes the component of the given element from the structure.
   2.731 +    ///
   2.732 +    /// \warning It is an error to give an element which is not in the
   2.733 +    /// structure.
   2.734 +    void eraseClass(int cdx) {
   2.735 +      int idx = classes[cdx].firstItem;
   2.736 +      items[items[idx].prev].next = firstFreeItem;
   2.737 +      firstFreeItem = idx;
   2.738 +
   2.739 +      if (classes[cdx].prev != -1) {
   2.740 +	classes[classes[cdx].prev].next = classes[cdx].next; 
   2.741 +      } else {
   2.742 +	firstClass = classes[cdx].next;
   2.743 +      }
   2.744 +      if (classes[cdx].next != -1) {
   2.745 +	classes[classes[cdx].next].prev = classes[cdx].prev; 
   2.746 +      }
   2.747 +      classes[cdx].next = firstFreeClass;
   2.748 +      firstFreeClass = cdx;
   2.749 +    }
   2.750 +
   2.751 +    /// \brief Lemon style iterator for the classes.
   2.752 +    ///
   2.753 +    /// ClassIt is a lemon style iterator for the components. It iterates
   2.754 +    /// on the ids of classes.
   2.755 +    class ClassIt {
   2.756 +    public:
   2.757 +      /// \brief Constructor of the iterator
   2.758 +      ///
   2.759 +      /// Constructor of the iterator
   2.760 +      ClassIt(const ExtendFindEnum& ufe) : extendFind(&ufe) {
   2.761 +        cdx = extendFind->firstClass;
   2.762 +      }
   2.763 +
   2.764 +      /// \brief Constructor to get invalid iterator
   2.765 +      ///
   2.766 +      /// Constructor to get invalid iterator
   2.767 +      ClassIt(Invalid) : extendFind(0), cdx(-1) {}
   2.768 +      
   2.769 +      /// \brief Increment operator
   2.770 +      ///
   2.771 +      /// It steps to the next representant item.
   2.772 +      ClassIt& operator++() {
   2.773 +        cdx = extendFind->classes[cdx].next;
   2.774 +        return *this;
   2.775 +      }
   2.776 +      
   2.777 +      /// \brief Conversion operator
   2.778 +      ///
   2.779 +      /// It converts the iterator to the current class id.
   2.780 +      operator int() const {
   2.781 +        return cdx;
   2.782 +      }
   2.783 +
   2.784 +      /// \brief Equality operator
   2.785 +      ///
   2.786 +      /// Equality operator
   2.787 +      bool operator==(const ClassIt& i) { 
   2.788 +        return i.cdx == cdx;
   2.789 +      }
   2.790 +
   2.791 +      /// \brief Inequality operator
   2.792 +      ///
   2.793 +      /// Inequality operator
   2.794 +      bool operator!=(const ClassIt& i) { 
   2.795 +        return i.cdx != cdx;
   2.796 +      }
   2.797 +      
   2.798 +    private:
   2.799 +      const ExtendFindEnum* extendFind;
   2.800 +      int cdx;
   2.801 +    };
   2.802 +
   2.803 +    /// \brief Lemon style iterator for the items of a component.
   2.804 +    ///
   2.805 +    /// ClassIt is a lemon style iterator for the components. It iterates
   2.806 +    /// on the items of a class. By example if you want to iterate on
   2.807 +    /// each items of each classes then you may write the next code.
   2.808 +    ///\code
   2.809 +    /// for (ClassIt cit(ufe); cit != INVALID; ++cit) {
   2.810 +    ///   std::cout << "Class: ";
   2.811 +    ///   for (ItemIt iit(ufe, cit); iit != INVALID; ++iit) {
   2.812 +    ///     std::cout << toString(iit) << ' ' << std::endl;
   2.813 +    ///   }
   2.814 +    ///   std::cout << std::endl;
   2.815 +    /// }
   2.816 +    ///\endcode
   2.817 +    class ItemIt {
   2.818 +    public:
   2.819 +      /// \brief Constructor of the iterator
   2.820 +      ///
   2.821 +      /// Constructor of the iterator. The iterator iterates
   2.822 +      /// on the class of the \c item.
   2.823 +      ItemIt(const ExtendFindEnum& ufe, int cls) : extendFind(&ufe) {
   2.824 +        fdx = idx = extendFind->classes[cls].firstItem;
   2.825 +      }
   2.826 +
   2.827 +      /// \brief Constructor to get invalid iterator
   2.828 +      ///
   2.829 +      /// Constructor to get invalid iterator
   2.830 +      ItemIt(Invalid) : extendFind(0), idx(-1) {}
   2.831 +      
   2.832 +      /// \brief Increment operator
   2.833 +      ///
   2.834 +      /// It steps to the next item in the class.
   2.835 +      ItemIt& operator++() {
   2.836 +        idx = extendFind->items[idx].next;
   2.837 +	if (fdx == idx) idx = -1;
   2.838 +        return *this;
   2.839 +      }
   2.840 +      
   2.841 +      /// \brief Conversion operator
   2.842 +      ///
   2.843 +      /// It converts the iterator to the current item.
   2.844 +      operator const Item&() const {
   2.845 +        return extendFind->items[idx].item;
   2.846 +      }
   2.847 +
   2.848 +      /// \brief Equality operator
   2.849 +      ///
   2.850 +      /// Equality operator
   2.851 +      bool operator==(const ItemIt& i) { 
   2.852 +        return i.idx == idx;
   2.853 +      }
   2.854 +
   2.855 +      /// \brief Inequality operator
   2.856 +      ///
   2.857 +      /// Inequality operator
   2.858 +      bool operator!=(const ItemIt& i) { 
   2.859 +        return i.idx != idx;
   2.860 +      }
   2.861 +      
   2.862 +    private:
   2.863 +      const ExtendFindEnum* extendFind;
   2.864 +      int idx, fdx;
   2.865 +    };
   2.866 +
   2.867 +  };
   2.868  
   2.869    //! @}
   2.870  
     3.1 --- a/test/max_matching_test.cc	Tue Oct 30 10:51:07 2007 +0000
     3.2 +++ b/test/max_matching_test.cc	Tue Oct 30 20:21:10 2007 +0000
     3.3 @@ -66,11 +66,12 @@
     3.4    g.addEdge(nodes[6], nodes[12]);
     3.5    
     3.6    MaxMatching<Graph> max_matching(g);
     3.7 -  max_matching.runEdmonds(0);
     3.8 +  max_matching.init();
     3.9 +  max_matching.startDense();
    3.10    
    3.11    int s=0;
    3.12    Graph::NodeMap<Node> mate(g,INVALID);
    3.13 -  max_matching.writeNMapNode(mate);
    3.14 +  max_matching.mateMap(mate);
    3.15    for(NodeIt v(g); v!=INVALID; ++v) {
    3.16      if ( mate[v]!=INVALID ) ++s;
    3.17    }
    3.18 @@ -82,43 +83,42 @@
    3.19  
    3.20    check ( size == max_matching.size(), "mate() returns a different size matching than max_matching.size()" );
    3.21  
    3.22 -  Graph::NodeMap<MaxMatching<Graph>::pos_enum> pos0(g);
    3.23 -  max_matching.writePos(pos0);
    3.24 +  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos0(g);
    3.25 +  max_matching.decomposition(pos0);
    3.26    
    3.27 -  max_matching.resetMatching();
    3.28 -  max_matching.runEdmonds(1);
    3.29 +  max_matching.init();
    3.30 +  max_matching.startSparse();
    3.31    s=0;
    3.32 -  max_matching.writeNMapNode(mate);
    3.33 +  max_matching.mateMap(mate);
    3.34    for(NodeIt v(g); v!=INVALID; ++v) {
    3.35      if ( mate[v]!=INVALID ) ++s;
    3.36    }
    3.37    check ( int(s/2) == size, "The size does not equal!" );
    3.38  
    3.39 -  Graph::NodeMap<MaxMatching<Graph>::pos_enum> pos1(g);
    3.40 -  max_matching.writePos(pos1);
    3.41 +  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos1(g);
    3.42 +  max_matching.decomposition(pos1);
    3.43  
    3.44    max_matching.run();
    3.45    s=0;
    3.46 -  max_matching.writeNMapNode(mate);
    3.47 +  max_matching.mateMap(mate);
    3.48    for(NodeIt v(g); v!=INVALID; ++v) {
    3.49      if ( mate[v]!=INVALID ) ++s;
    3.50    }
    3.51    check ( int(s/2) == size, "The size does not equal!" ); 
    3.52    
    3.53 -  Graph::NodeMap<MaxMatching<Graph>::pos_enum> pos2(g);
    3.54 -  max_matching.writePos(pos2);
    3.55 +  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos2(g);
    3.56 +  max_matching.decomposition(pos2);
    3.57  
    3.58 -  max_matching.resetMatching();
    3.59    max_matching.run();
    3.60    s=0;
    3.61 -  max_matching.writeNMapNode(mate);
    3.62 +  max_matching.mateMap(mate);
    3.63    for(NodeIt v(g); v!=INVALID; ++v) {
    3.64      if ( mate[v]!=INVALID ) ++s;
    3.65    }
    3.66    check ( int(s/2) == size, "The size does not equal!" ); 
    3.67    
    3.68 -  Graph::NodeMap<MaxMatching<Graph>::pos_enum> pos(g);
    3.69 -  max_matching.writePos(pos);
    3.70 +  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos(g);
    3.71 +  max_matching.decomposition(pos);
    3.72     
    3.73    bool ismatching=true;
    3.74    for(NodeIt v(g); v!=INVALID; ++v) {
    3.75 @@ -139,8 +139,10 @@
    3.76  
    3.77    bool noedge=true;
    3.78    for(UEdgeIt e(g); e!=INVALID; ++e) {
    3.79 -   if ( (pos[g.target(e)]==max_matching.C && pos[g.source(e)]==max_matching.D) || 
    3.80 -	 (pos[g.target(e)]==max_matching.D && pos[g.source(e)]==max_matching.C) )
    3.81 +   if ( (pos[g.target(e)]==max_matching.C && 
    3.82 +	 pos[g.source(e)]==max_matching.D) || 
    3.83 +	 (pos[g.target(e)]==max_matching.D && 
    3.84 +	  pos[g.source(e)]==max_matching.C) )
    3.85        noedge=false; 
    3.86    }
    3.87    check ( noedge, "There are edges between D and C!" );
     4.1 --- a/test/unionfind_test.cc	Tue Oct 30 10:51:07 2007 +0000
     4.2 +++ b/test/unionfind_test.cc	Tue Oct 30 20:21:10 2007 +0000
     4.3 @@ -49,7 +49,7 @@
     4.4    U.insert(1);
     4.5    U.insert(2);
     4.6  
     4.7 -  check(U.join(1,2),"Test failed.");
     4.8 +  check(U.join(1,2) != -1,"Test failed.");
     4.9  
    4.10    U.insert(3);
    4.11    U.insert(4);
    4.12 @@ -57,66 +57,58 @@
    4.13    U.insert(6);
    4.14    U.insert(7);
    4.15  
    4.16 -  check(U.join(1,4),"Test failed.");
    4.17 -  check(!U.join(2,4),"Test failed.");
    4.18 -  check(U.join(3,5),"Test failed.");
    4.19  
    4.20 -  U.insert(8,5);
    4.21 +  check(U.join(1,4) != -1,"Test failed.");
    4.22 +  check(U.join(2,4) == -1,"Test failed.");
    4.23 +  check(U.join(3,5) != -1,"Test failed.");
    4.24  
    4.25 -  check(U.size(4) == 3,"Test failed.");
    4.26 -  check(U.size(5) == 3,"Test failed.");
    4.27 -  check(U.size(6) == 1,"Test failed.");
    4.28 -  check(U.size(2) == 3,"Test failed.");
    4.29 +
    4.30 +  U.insert(8,U.find(5));
    4.31 +
    4.32 +
    4.33 +  check(U.size(U.find(4)) == 3,"Test failed.");
    4.34 +  check(U.size(U.find(5)) == 3,"Test failed.");
    4.35 +  check(U.size(U.find(6)) == 1,"Test failed.");
    4.36 +  check(U.size(U.find(2)) == 3,"Test failed.");
    4.37 +
    4.38  
    4.39    U.insert(9);
    4.40 -  U.insert(10,9);
    4.41 +  U.insert(10,U.find(9));
    4.42  
    4.43 -  check(U.join(8,10),"Test failed.");
    4.44  
    4.45 -  check(U.move(9,4),"Test failed.");
    4.46 -  check(!U.move(9,2),"Test failed.");
    4.47 +  check(U.join(8,10) != -1,"Test failed.");
    4.48  
    4.49 -  check(U.size(4) == 4,"Test failed.");
    4.50 -  check(U.size(9) == 4,"Test failed.");
    4.51  
    4.52 -  check(U.move(5,6),"Test failed.");
    4.53 +  check(U.size(U.find(4)) == 3,"Test failed.");
    4.54 +  check(U.size(U.find(9)) == 5,"Test failed.");
    4.55  
    4.56 -  check(U.size(5) == 2,"Test failed.");
    4.57 -  check(U.size(8) == 3,"Test failed.");
    4.58 -
    4.59 -  check(U.move(7,10),"Test failed.");
    4.60 -  check(U.size(7) == 4,"Test failed.");
    4.61 +  check(U.size(U.find(8)) == 5,"Test failed.");
    4.62  
    4.63    U.erase(9);
    4.64    U.erase(1);
    4.65  
    4.66 -  check(U.size(4) == 2,"Test failed.");
    4.67 -  check(U.size(2) == 2,"Test failed.");
    4.68 +  check(U.size(U.find(10)) == 4,"Test failed.");
    4.69 +  check(U.size(U.find(2)) == 2,"Test failed.");
    4.70  
    4.71    U.erase(6);
    4.72 -  U.split(8);
    4.73 +  U.split(U.find(8));
    4.74  
    4.75 -  check(U.size(4) == 2,"Test failed.");
    4.76 -  check(U.size(3) == 1,"Test failed.");
    4.77 -  check(U.size(2) == 2,"Test failed.");
    4.78  
    4.79 -  check(U.join(3,4),"Test failed.");
    4.80 -  check(!U.join(2,4),"Test failed.");
    4.81 +  check(U.size(U.find(4)) == 2,"Test failed.");
    4.82 +  check(U.size(U.find(3)) == 1,"Test failed.");
    4.83 +  check(U.size(U.find(2)) == 2,"Test failed.");
    4.84  
    4.85 -  check(U.size(4) == 3,"Test failed.");
    4.86 -  check(U.size(3) == 3,"Test failed.");
    4.87 -  check(U.size(2) == 3,"Test failed.");
    4.88  
    4.89 -  U.makeRep(4);
    4.90 -  U.makeRep(3);
    4.91 -  U.makeRep(2);
    4.92 +  check(U.join(3,4) != -1,"Test failed.");
    4.93 +  check(U.join(2,4) == -1,"Test failed.");
    4.94  
    4.95 -  check(U.size(4) == 3,"Test failed.");
    4.96 -  check(U.size(3) == 3,"Test failed.");
    4.97 -  check(U.size(2) == 3,"Test failed.");
    4.98  
    4.99 +  check(U.size(U.find(4)) == 3,"Test failed.");
   4.100 +  check(U.size(U.find(3)) == 3,"Test failed.");
   4.101 +  check(U.size(U.find(2)) == 3,"Test failed.");
   4.102  
   4.103 -  U.eraseClass(4);
   4.104 -  U.eraseClass(7);
   4.105 +  U.eraseClass(U.find(4));
   4.106 +  U.eraseClass(U.find(7));
   4.107  
   4.108 +  return 0;
   4.109  }