lemon/unionfind.h
changeset 207 574b963d0275
child 209 765619b7cbb2
equal deleted inserted replaced
-1:000000000000 0:a2fb3411c97d
       
     1 /* -*- C++ -*-
       
     2  *
       
     3  * This file is a part of LEMON, a generic C++ optimization library
       
     4  *
       
     5  * Copyright (C) 2003-2008
       
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
       
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
       
     8  *
       
     9  * Permission to use, modify and distribute this software is granted
       
    10  * provided that this copyright notice appears in all copies. For
       
    11  * precise terms see the accompanying LICENSE file.
       
    12  *
       
    13  * This software is provided "AS IS" with no warranty of any kind,
       
    14  * express or implied, and with no claim as to its suitability for any
       
    15  * purpose.
       
    16  *
       
    17  */
       
    18 
       
    19 #ifndef LEMON_UNION_FIND_H
       
    20 #define LEMON_UNION_FIND_H
       
    21 
       
    22 //!\ingroup auxdat
       
    23 //!\file
       
    24 //!\brief Union-Find data structures.
       
    25 //!
       
    26 
       
    27 #include <vector>
       
    28 #include <list>
       
    29 #include <utility>
       
    30 #include <algorithm>
       
    31 #include <functional>
       
    32 
       
    33 #include <lemon/bits/invalid.h>
       
    34 
       
    35 namespace lemon {
       
    36 
       
    37   /// \ingroup auxdat
       
    38   ///
       
    39   /// \brief A \e Union-Find data structure implementation
       
    40   ///
       
    41   /// The class implements the \e Union-Find data structure. 
       
    42   /// The union operation uses rank heuristic, while
       
    43   /// the find operation uses path compression.
       
    44   /// This is a very simple but efficient implementation, providing 
       
    45   /// only four methods: join (union), find, insert and size.
       
    46   /// For more features see the \ref UnionFindEnum class.
       
    47   ///
       
    48   /// It is primarily used in Kruskal algorithm for finding minimal
       
    49   /// cost spanning tree in a graph.
       
    50   /// \sa kruskal()
       
    51   ///
       
    52   /// \pre You need to add all the elements by the \ref insert()
       
    53   /// method.  
       
    54   template <typename _ItemIntMap>
       
    55   class UnionFind {
       
    56   public:
       
    57 
       
    58     typedef _ItemIntMap ItemIntMap;
       
    59     typedef typename ItemIntMap::Key Item;
       
    60 
       
    61   private:
       
    62     // If the items vector stores negative value for an item then
       
    63     // that item is root item and it has -items[it] component size.
       
    64     // Else the items[it] contains the index of the parent.
       
    65     std::vector<int> items;
       
    66     ItemIntMap& index;
       
    67 
       
    68     bool rep(int idx) const {
       
    69       return items[idx] < 0;
       
    70     }
       
    71 
       
    72     int repIndex(int idx) const {
       
    73       int k = idx;
       
    74       while (!rep(k)) {
       
    75         k = items[k] ;
       
    76       }
       
    77       while (idx != k) {
       
    78         int next = items[idx];
       
    79         const_cast<int&>(items[idx]) = k;
       
    80         idx = next;
       
    81       }
       
    82       return k;
       
    83     }
       
    84 
       
    85   public:
       
    86 
       
    87     /// \brief Constructor
       
    88     ///
       
    89     /// Constructor of the UnionFind class. You should give an item to
       
    90     /// integer map which will be used from the data structure. If you
       
    91     /// modify directly this map that may cause segmentation fault,
       
    92     /// invalid data structure, or infinite loop when you use again
       
    93     /// the union-find.
       
    94     UnionFind(ItemIntMap& m) : index(m) {}
       
    95 
       
    96     /// \brief Returns the index of the element's component.
       
    97     ///
       
    98     /// The method returns the index of the element's component.
       
    99     /// This is an integer between zero and the number of inserted elements.
       
   100     ///
       
   101     int find(const Item& a) {
       
   102       return repIndex(index[a]);
       
   103     }
       
   104 
       
   105     /// \brief Clears the union-find data structure
       
   106     ///
       
   107     /// Erase each item from the data structure.
       
   108     void clear() {
       
   109       items.clear();
       
   110     }
       
   111 
       
   112     /// \brief Inserts a new element into the structure.
       
   113     ///
       
   114     /// This method inserts a new element into the data structure. 
       
   115     ///
       
   116     /// The method returns the index of the new component.
       
   117     int insert(const Item& a) {
       
   118       int n = items.size();
       
   119       items.push_back(-1);
       
   120       index.set(a,n);
       
   121       return n;
       
   122     }
       
   123 
       
   124     /// \brief Joining the components of element \e a and element \e b.
       
   125     ///
       
   126     /// This is the \e union operation of the Union-Find structure. 
       
   127     /// Joins the component of element \e a and component of
       
   128     /// element \e b. If \e a and \e b are in the same component then
       
   129     /// it returns false otherwise it returns true.
       
   130     bool join(const Item& a, const Item& b) {
       
   131       int ka = repIndex(index[a]);
       
   132       int kb = repIndex(index[b]);
       
   133 
       
   134       if ( ka == kb ) 
       
   135 	return false;
       
   136 
       
   137       if (items[ka] < items[kb]) {
       
   138 	items[ka] += items[kb];
       
   139 	items[kb] = ka;
       
   140       } else {
       
   141 	items[kb] += items[ka];
       
   142 	items[ka] = kb;
       
   143       }
       
   144       return true;
       
   145     }
       
   146 
       
   147     /// \brief Returns the size of the component of element \e a.
       
   148     ///
       
   149     /// Returns the size of the component of element \e a.
       
   150     int size(const Item& a) {
       
   151       int k = repIndex(index[a]);
       
   152       return - items[k];
       
   153     }
       
   154 
       
   155   };
       
   156 
       
   157   /// \ingroup auxdat
       
   158   ///
       
   159   /// \brief A \e Union-Find data structure implementation which
       
   160   /// is able to enumerate the components.
       
   161   ///
       
   162   /// The class implements a \e Union-Find data structure
       
   163   /// which is able to enumerate the components and the items in
       
   164   /// a component. If you don't need this feature then perhaps it's
       
   165   /// better to use the \ref UnionFind class which is more efficient.
       
   166   ///
       
   167   /// The union operation uses rank heuristic, while
       
   168   /// the find operation uses path compression.
       
   169   ///
       
   170   /// \pre You need to add all the elements by the \ref insert()
       
   171   /// method.
       
   172   ///
       
   173   template <typename _ItemIntMap>
       
   174   class UnionFindEnum {
       
   175   public:
       
   176     
       
   177     typedef _ItemIntMap ItemIntMap;
       
   178     typedef typename ItemIntMap::Key Item;
       
   179 
       
   180   private:
       
   181     
       
   182     ItemIntMap& index;
       
   183 
       
   184     // If the parent stores negative value for an item then that item
       
   185     // is root item and it has ~(items[it].parent) component id.  Else
       
   186     // the items[it].parent contains the index of the parent.
       
   187     //
       
   188     // The \c next and \c prev provides the double-linked
       
   189     // cyclic list of one component's items.
       
   190     struct ItemT {
       
   191       int parent;
       
   192       Item item;
       
   193 
       
   194       int next, prev;
       
   195     };
       
   196 
       
   197     std::vector<ItemT> items;
       
   198     int firstFreeItem;
       
   199 
       
   200     struct ClassT {
       
   201       int size;
       
   202       int firstItem;
       
   203       int next, prev;
       
   204     };
       
   205     
       
   206     std::vector<ClassT> classes;
       
   207     int firstClass, firstFreeClass;
       
   208 
       
   209     int newClass() {
       
   210       if (firstFreeClass == -1) {
       
   211 	int cdx = classes.size();
       
   212 	classes.push_back(ClassT());
       
   213 	return cdx;
       
   214       } else {
       
   215 	int cdx = firstFreeClass;
       
   216 	firstFreeClass = classes[firstFreeClass].next;
       
   217 	return cdx;
       
   218       }
       
   219     }
       
   220 
       
   221     int newItem() {
       
   222       if (firstFreeItem == -1) {
       
   223 	int idx = items.size();
       
   224 	items.push_back(ItemT());
       
   225 	return idx;
       
   226       } else {
       
   227 	int idx = firstFreeItem;
       
   228 	firstFreeItem = items[firstFreeItem].next;
       
   229 	return idx;
       
   230       }
       
   231     }
       
   232 
       
   233 
       
   234     bool rep(int idx) const {
       
   235       return items[idx].parent < 0;
       
   236     }
       
   237 
       
   238     int repIndex(int idx) const {
       
   239       int k = idx;
       
   240       while (!rep(k)) {
       
   241         k = items[k].parent;
       
   242       }
       
   243       while (idx != k) {
       
   244         int next = items[idx].parent;
       
   245         const_cast<int&>(items[idx].parent) = k;
       
   246         idx = next;
       
   247       }
       
   248       return k;
       
   249     }
       
   250 
       
   251     int classIndex(int idx) const {
       
   252       return ~(items[repIndex(idx)].parent);
       
   253     }
       
   254 
       
   255     void singletonItem(int idx) {
       
   256       items[idx].next = idx;
       
   257       items[idx].prev = idx;
       
   258     }
       
   259 
       
   260     void laceItem(int idx, int rdx) {
       
   261       items[idx].prev = rdx;
       
   262       items[idx].next = items[rdx].next;
       
   263       items[items[rdx].next].prev = idx;
       
   264       items[rdx].next = idx;
       
   265     }
       
   266 
       
   267     void unlaceItem(int idx) {
       
   268       items[items[idx].prev].next = items[idx].next;
       
   269       items[items[idx].next].prev = items[idx].prev;
       
   270       
       
   271       items[idx].next = firstFreeItem;
       
   272       firstFreeItem = idx;
       
   273     }
       
   274 
       
   275     void spliceItems(int ak, int bk) {
       
   276       items[items[ak].prev].next = bk;
       
   277       items[items[bk].prev].next = ak;
       
   278       int tmp = items[ak].prev;
       
   279       items[ak].prev = items[bk].prev;
       
   280       items[bk].prev = tmp;
       
   281         
       
   282     }
       
   283 
       
   284     void laceClass(int cls) {
       
   285       if (firstClass != -1) {
       
   286         classes[firstClass].prev = cls;
       
   287       }
       
   288       classes[cls].next = firstClass;
       
   289       classes[cls].prev = -1;
       
   290       firstClass = cls;
       
   291     } 
       
   292 
       
   293     void unlaceClass(int cls) {
       
   294       if (classes[cls].prev != -1) {
       
   295         classes[classes[cls].prev].next = classes[cls].next;
       
   296       } else {
       
   297         firstClass = classes[cls].next;
       
   298       }
       
   299       if (classes[cls].next != -1) {
       
   300         classes[classes[cls].next].prev = classes[cls].prev;
       
   301       }
       
   302       
       
   303       classes[cls].next = firstFreeClass;
       
   304       firstFreeClass = cls;
       
   305     } 
       
   306 
       
   307   public:
       
   308 
       
   309     UnionFindEnum(ItemIntMap& _index) 
       
   310       : index(_index), items(), firstFreeItem(-1), 
       
   311 	firstClass(-1), firstFreeClass(-1) {}
       
   312     
       
   313     /// \brief Inserts the given element into a new component.
       
   314     ///
       
   315     /// This method creates a new component consisting only of the
       
   316     /// given element.
       
   317     ///
       
   318     int insert(const Item& item) {
       
   319       int idx = newItem();
       
   320 
       
   321       index.set(item, idx);
       
   322 
       
   323       singletonItem(idx);
       
   324       items[idx].item = item;
       
   325 
       
   326       int cdx = newClass();
       
   327 
       
   328       items[idx].parent = ~cdx;
       
   329 
       
   330       laceClass(cdx);
       
   331       classes[cdx].size = 1;
       
   332       classes[cdx].firstItem = idx;
       
   333 
       
   334       firstClass = cdx;
       
   335       
       
   336       return cdx;
       
   337     }
       
   338 
       
   339     /// \brief Inserts the given element into the component of the others.
       
   340     ///
       
   341     /// This methods inserts the element \e a into the component of the
       
   342     /// element \e comp. 
       
   343     void insert(const Item& item, int cls) {
       
   344       int rdx = classes[cls].firstItem;
       
   345       int idx = newItem();
       
   346 
       
   347       index.set(item, idx);
       
   348 
       
   349       laceItem(idx, rdx);
       
   350 
       
   351       items[idx].item = item;
       
   352       items[idx].parent = rdx;
       
   353 
       
   354       ++classes[~(items[rdx].parent)].size;
       
   355     }
       
   356 
       
   357     /// \brief Clears the union-find data structure
       
   358     ///
       
   359     /// Erase each item from the data structure.
       
   360     void clear() {
       
   361       items.clear();
       
   362       firstClass = -1;
       
   363       firstFreeItem = -1;
       
   364     }
       
   365 
       
   366     /// \brief Finds the component of the given element.
       
   367     ///
       
   368     /// The method returns the component id of the given element.
       
   369     int find(const Item &item) const {
       
   370       return ~(items[repIndex(index[item])].parent);
       
   371     }
       
   372 
       
   373     /// \brief Joining the component of element \e a and element \e b.
       
   374     ///
       
   375     /// This is the \e union operation of the Union-Find structure. 
       
   376     /// Joins the component of element \e a and component of
       
   377     /// element \e b. If \e a and \e b are in the same component then
       
   378     /// returns -1 else returns the remaining class.
       
   379     int join(const Item& a, const Item& b) {
       
   380 
       
   381       int ak = repIndex(index[a]);
       
   382       int bk = repIndex(index[b]);
       
   383 
       
   384       if (ak == bk) {
       
   385 	return -1;
       
   386       }
       
   387 
       
   388       int acx = ~(items[ak].parent);
       
   389       int bcx = ~(items[bk].parent);
       
   390 
       
   391       int rcx;
       
   392 
       
   393       if (classes[acx].size > classes[bcx].size) {
       
   394 	classes[acx].size += classes[bcx].size;
       
   395 	items[bk].parent = ak;
       
   396         unlaceClass(bcx);
       
   397 	rcx = acx;
       
   398       } else {
       
   399 	classes[bcx].size += classes[acx].size;
       
   400 	items[ak].parent = bk;
       
   401         unlaceClass(acx);
       
   402 	rcx = bcx;
       
   403       }
       
   404       spliceItems(ak, bk);
       
   405 
       
   406       return rcx;
       
   407     }
       
   408 
       
   409     /// \brief Returns the size of the class.
       
   410     ///
       
   411     /// Returns the size of the class.
       
   412     int size(int cls) const {
       
   413       return classes[cls].size;
       
   414     }
       
   415 
       
   416     /// \brief Splits up the component. 
       
   417     ///
       
   418     /// Splitting the component into singleton components (component
       
   419     /// of size one).
       
   420     void split(int cls) {
       
   421       int fdx = classes[cls].firstItem;
       
   422       int idx = items[fdx].next;
       
   423       while (idx != fdx) {
       
   424         int next = items[idx].next;
       
   425 
       
   426 	singletonItem(idx);
       
   427 
       
   428 	int cdx = newClass();        
       
   429         items[idx].parent = ~cdx;
       
   430 
       
   431 	laceClass(cdx);
       
   432 	classes[cdx].size = 1;
       
   433 	classes[cdx].firstItem = idx;
       
   434         
       
   435         idx = next;
       
   436       }
       
   437 
       
   438       items[idx].prev = idx;
       
   439       items[idx].next = idx;
       
   440 
       
   441       classes[~(items[idx].parent)].size = 1;
       
   442       
       
   443     }
       
   444 
       
   445     /// \brief Removes the given element from the structure.
       
   446     ///
       
   447     /// Removes the element from its component and if the component becomes
       
   448     /// empty then removes that component from the component list.
       
   449     ///
       
   450     /// \warning It is an error to remove an element which is not in
       
   451     /// the structure.
       
   452     /// \warning This running time of this operation is proportional to the
       
   453     /// number of the items in this class.
       
   454     void erase(const Item& item) {
       
   455       int idx = index[item];
       
   456       int fdx = items[idx].next;
       
   457 
       
   458       int cdx = classIndex(idx);
       
   459       if (idx == fdx) {
       
   460 	unlaceClass(cdx);
       
   461 	items[idx].next = firstFreeItem;
       
   462 	firstFreeItem = idx;
       
   463 	return;
       
   464       } else {
       
   465 	classes[cdx].firstItem = fdx;
       
   466 	--classes[cdx].size;
       
   467 	items[fdx].parent = ~cdx;
       
   468 
       
   469 	unlaceItem(idx);
       
   470 	idx = items[fdx].next;
       
   471 	while (idx != fdx) {
       
   472 	  items[idx].parent = fdx;
       
   473 	  idx = items[idx].next;
       
   474 	}
       
   475           
       
   476       }
       
   477 
       
   478     }
       
   479 
       
   480     /// \brief Gives back a representant item of the component.
       
   481     ///
       
   482     /// Gives back a representant item of the component.
       
   483     Item item(int cls) const {
       
   484       return items[classes[cls].firstItem].item;
       
   485     }
       
   486 
       
   487     /// \brief Removes the component of the given element from the structure.
       
   488     ///
       
   489     /// Removes the component of the given element from the structure.
       
   490     ///
       
   491     /// \warning It is an error to give an element which is not in the
       
   492     /// structure.
       
   493     void eraseClass(int cls) {
       
   494       int fdx = classes[cls].firstItem;
       
   495       unlaceClass(cls);
       
   496       items[items[fdx].prev].next = firstFreeItem;
       
   497       firstFreeItem = fdx;
       
   498     }
       
   499 
       
   500     /// \brief Lemon style iterator for the representant items.
       
   501     ///
       
   502     /// ClassIt is a lemon style iterator for the components. It iterates
       
   503     /// on the ids of the classes.
       
   504     class ClassIt {
       
   505     public:
       
   506       /// \brief Constructor of the iterator
       
   507       ///
       
   508       /// Constructor of the iterator
       
   509       ClassIt(const UnionFindEnum& ufe) : unionFind(&ufe) {
       
   510         cdx = unionFind->firstClass;
       
   511       }
       
   512 
       
   513       /// \brief Constructor to get invalid iterator
       
   514       ///
       
   515       /// Constructor to get invalid iterator
       
   516       ClassIt(Invalid) : unionFind(0), cdx(-1) {}
       
   517       
       
   518       /// \brief Increment operator
       
   519       ///
       
   520       /// It steps to the next representant item.
       
   521       ClassIt& operator++() {
       
   522         cdx = unionFind->classes[cdx].next;
       
   523         return *this;
       
   524       }
       
   525       
       
   526       /// \brief Conversion operator
       
   527       ///
       
   528       /// It converts the iterator to the current representant item.
       
   529       operator int() const {
       
   530         return cdx;
       
   531       }
       
   532 
       
   533       /// \brief Equality operator
       
   534       ///
       
   535       /// Equality operator
       
   536       bool operator==(const ClassIt& i) { 
       
   537         return i.cdx == cdx;
       
   538       }
       
   539 
       
   540       /// \brief Inequality operator
       
   541       ///
       
   542       /// Inequality operator
       
   543       bool operator!=(const ClassIt& i) { 
       
   544         return i.cdx != cdx;
       
   545       }
       
   546       
       
   547     private:
       
   548       const UnionFindEnum* unionFind;
       
   549       int cdx;
       
   550     };
       
   551 
       
   552     /// \brief Lemon style iterator for the items of a component.
       
   553     ///
       
   554     /// ClassIt is a lemon style iterator for the components. It iterates
       
   555     /// on the items of a class. By example if you want to iterate on
       
   556     /// each items of each classes then you may write the next code.
       
   557     ///\code
       
   558     /// for (ClassIt cit(ufe); cit != INVALID; ++cit) {
       
   559     ///   std::cout << "Class: ";
       
   560     ///   for (ItemIt iit(ufe, cit); iit != INVALID; ++iit) {
       
   561     ///     std::cout << toString(iit) << ' ' << std::endl;
       
   562     ///   }
       
   563     ///   std::cout << std::endl;
       
   564     /// }
       
   565     ///\endcode
       
   566     class ItemIt {
       
   567     public:
       
   568       /// \brief Constructor of the iterator
       
   569       ///
       
   570       /// Constructor of the iterator. The iterator iterates
       
   571       /// on the class of the \c item.
       
   572       ItemIt(const UnionFindEnum& ufe, int cls) : unionFind(&ufe) {
       
   573         fdx = idx = unionFind->classes[cls].firstItem;
       
   574       }
       
   575 
       
   576       /// \brief Constructor to get invalid iterator
       
   577       ///
       
   578       /// Constructor to get invalid iterator
       
   579       ItemIt(Invalid) : unionFind(0), idx(-1) {}
       
   580       
       
   581       /// \brief Increment operator
       
   582       ///
       
   583       /// It steps to the next item in the class.
       
   584       ItemIt& operator++() {
       
   585         idx = unionFind->items[idx].next;
       
   586         if (idx == fdx) idx = -1;
       
   587         return *this;
       
   588       }
       
   589       
       
   590       /// \brief Conversion operator
       
   591       ///
       
   592       /// It converts the iterator to the current item.
       
   593       operator const Item&() const {
       
   594         return unionFind->items[idx].item;
       
   595       }
       
   596 
       
   597       /// \brief Equality operator
       
   598       ///
       
   599       /// Equality operator
       
   600       bool operator==(const ItemIt& i) { 
       
   601         return i.idx == idx;
       
   602       }
       
   603 
       
   604       /// \brief Inequality operator
       
   605       ///
       
   606       /// Inequality operator
       
   607       bool operator!=(const ItemIt& i) { 
       
   608         return i.idx != idx;
       
   609       }
       
   610       
       
   611     private:
       
   612       const UnionFindEnum* unionFind;
       
   613       int idx, fdx;
       
   614     };
       
   615 
       
   616   };
       
   617 
       
   618   /// \ingroup auxdat
       
   619   ///
       
   620   /// \brief A \e Extend-Find data structure implementation which
       
   621   /// is able to enumerate the components.
       
   622   ///
       
   623   /// The class implements an \e Extend-Find data structure which is
       
   624   /// able to enumerate the components and the items in a
       
   625   /// component. The data structure is a simplification of the
       
   626   /// Union-Find structure, and it does not allow to merge two components.
       
   627   ///
       
   628   /// \pre You need to add all the elements by the \ref insert()
       
   629   /// method.
       
   630   template <typename _ItemIntMap>
       
   631   class ExtendFindEnum {
       
   632   public:
       
   633     
       
   634     typedef _ItemIntMap ItemIntMap;
       
   635     typedef typename ItemIntMap::Key Item;
       
   636 
       
   637   private:
       
   638     
       
   639     ItemIntMap& index;
       
   640 
       
   641     struct ItemT {
       
   642       int cls;
       
   643       Item item;
       
   644       int next, prev;
       
   645     };
       
   646 
       
   647     std::vector<ItemT> items;
       
   648     int firstFreeItem;
       
   649 
       
   650     struct ClassT {
       
   651       int firstItem;
       
   652       int next, prev;
       
   653     };
       
   654 
       
   655     std::vector<ClassT> classes;
       
   656 
       
   657     int firstClass, firstFreeClass;
       
   658 
       
   659     int newClass() {
       
   660       if (firstFreeClass != -1) {
       
   661 	int cdx = firstFreeClass;
       
   662 	firstFreeClass = classes[cdx].next;
       
   663 	return cdx;
       
   664       } else {
       
   665 	classes.push_back(ClassT());
       
   666 	return classes.size() - 1;
       
   667       }
       
   668     }
       
   669 
       
   670     int newItem() {
       
   671       if (firstFreeItem != -1) {
       
   672 	int idx = firstFreeItem;
       
   673 	firstFreeItem = items[idx].next;
       
   674 	return idx;
       
   675       } else {
       
   676 	items.push_back(ItemT());
       
   677 	return items.size() - 1;
       
   678       }
       
   679     }
       
   680 
       
   681   public:
       
   682 
       
   683     /// \brief Constructor
       
   684     ExtendFindEnum(ItemIntMap& _index) 
       
   685       : index(_index), items(), firstFreeItem(-1), 
       
   686 	classes(), firstClass(-1), firstFreeClass(-1) {}
       
   687     
       
   688     /// \brief Inserts the given element into a new component.
       
   689     ///
       
   690     /// This method creates a new component consisting only of the
       
   691     /// given element.
       
   692     int insert(const Item& item) {
       
   693       int cdx = newClass();
       
   694       classes[cdx].prev = -1;
       
   695       classes[cdx].next = firstClass;
       
   696       if (firstClass != -1) {
       
   697 	classes[firstClass].prev = cdx;
       
   698       }
       
   699       firstClass = cdx;
       
   700       
       
   701       int idx = newItem();
       
   702       items[idx].item = item;
       
   703       items[idx].cls = cdx;
       
   704       items[idx].prev = idx;
       
   705       items[idx].next = idx;
       
   706 
       
   707       classes[cdx].firstItem = idx;
       
   708 
       
   709       index.set(item, idx);
       
   710       
       
   711       return cdx;
       
   712     }
       
   713 
       
   714     /// \brief Inserts the given element into the given component.
       
   715     ///
       
   716     /// This methods inserts the element \e item a into the \e cls class.
       
   717     void insert(const Item& item, int cls) {
       
   718       int idx = newItem();
       
   719       int rdx = classes[cls].firstItem;
       
   720       items[idx].item = item;
       
   721       items[idx].cls = cls;
       
   722 
       
   723       items[idx].prev = rdx;
       
   724       items[idx].next = items[rdx].next;
       
   725       items[items[rdx].next].prev = idx;
       
   726       items[rdx].next = idx;
       
   727 
       
   728       index.set(item, idx);
       
   729     }
       
   730 
       
   731     /// \brief Clears the union-find data structure
       
   732     ///
       
   733     /// Erase each item from the data structure.
       
   734     void clear() {
       
   735       items.clear();
       
   736       classes.clear;
       
   737       firstClass = firstFreeClass = firstFreeItem = -1;
       
   738     }
       
   739 
       
   740     /// \brief Gives back the class of the \e item.
       
   741     ///
       
   742     /// Gives back the class of the \e item.
       
   743     int find(const Item &item) const {
       
   744       return items[index[item]].cls;
       
   745     }
       
   746 
       
   747     /// \brief Gives back a representant item of the component.
       
   748     ///
       
   749     /// Gives back a representant item of the component.
       
   750     Item item(int cls) const {
       
   751       return items[classes[cls].firstItem].item;
       
   752     }
       
   753     
       
   754     /// \brief Removes the given element from the structure.
       
   755     ///
       
   756     /// Removes the element from its component and if the component becomes
       
   757     /// empty then removes that component from the component list.
       
   758     ///
       
   759     /// \warning It is an error to remove an element which is not in
       
   760     /// the structure.
       
   761     void erase(const Item &item) {
       
   762       int idx = index[item];
       
   763       int cdx = items[idx].cls;
       
   764       
       
   765       if (idx == items[idx].next) {
       
   766 	if (classes[cdx].prev != -1) {
       
   767 	  classes[classes[cdx].prev].next = classes[cdx].next; 
       
   768 	} else {
       
   769 	  firstClass = classes[cdx].next;
       
   770 	}
       
   771 	if (classes[cdx].next != -1) {
       
   772 	  classes[classes[cdx].next].prev = classes[cdx].prev; 
       
   773 	}
       
   774 	classes[cdx].next = firstFreeClass;
       
   775 	firstFreeClass = cdx;
       
   776       } else {
       
   777 	classes[cdx].firstItem = items[idx].next;
       
   778 	items[items[idx].next].prev = items[idx].prev;
       
   779 	items[items[idx].prev].next = items[idx].next;
       
   780       }
       
   781       items[idx].next = firstFreeItem;
       
   782       firstFreeItem = idx;
       
   783 	
       
   784     }    
       
   785 
       
   786     
       
   787     /// \brief Removes the component of the given element from the structure.
       
   788     ///
       
   789     /// Removes the component of the given element from the structure.
       
   790     ///
       
   791     /// \warning It is an error to give an element which is not in the
       
   792     /// structure.
       
   793     void eraseClass(int cdx) {
       
   794       int idx = classes[cdx].firstItem;
       
   795       items[items[idx].prev].next = firstFreeItem;
       
   796       firstFreeItem = idx;
       
   797 
       
   798       if (classes[cdx].prev != -1) {
       
   799 	classes[classes[cdx].prev].next = classes[cdx].next; 
       
   800       } else {
       
   801 	firstClass = classes[cdx].next;
       
   802       }
       
   803       if (classes[cdx].next != -1) {
       
   804 	classes[classes[cdx].next].prev = classes[cdx].prev; 
       
   805       }
       
   806       classes[cdx].next = firstFreeClass;
       
   807       firstFreeClass = cdx;
       
   808     }
       
   809 
       
   810     /// \brief Lemon style iterator for the classes.
       
   811     ///
       
   812     /// ClassIt is a lemon style iterator for the components. It iterates
       
   813     /// on the ids of classes.
       
   814     class ClassIt {
       
   815     public:
       
   816       /// \brief Constructor of the iterator
       
   817       ///
       
   818       /// Constructor of the iterator
       
   819       ClassIt(const ExtendFindEnum& ufe) : extendFind(&ufe) {
       
   820         cdx = extendFind->firstClass;
       
   821       }
       
   822 
       
   823       /// \brief Constructor to get invalid iterator
       
   824       ///
       
   825       /// Constructor to get invalid iterator
       
   826       ClassIt(Invalid) : extendFind(0), cdx(-1) {}
       
   827       
       
   828       /// \brief Increment operator
       
   829       ///
       
   830       /// It steps to the next representant item.
       
   831       ClassIt& operator++() {
       
   832         cdx = extendFind->classes[cdx].next;
       
   833         return *this;
       
   834       }
       
   835       
       
   836       /// \brief Conversion operator
       
   837       ///
       
   838       /// It converts the iterator to the current class id.
       
   839       operator int() const {
       
   840         return cdx;
       
   841       }
       
   842 
       
   843       /// \brief Equality operator
       
   844       ///
       
   845       /// Equality operator
       
   846       bool operator==(const ClassIt& i) { 
       
   847         return i.cdx == cdx;
       
   848       }
       
   849 
       
   850       /// \brief Inequality operator
       
   851       ///
       
   852       /// Inequality operator
       
   853       bool operator!=(const ClassIt& i) { 
       
   854         return i.cdx != cdx;
       
   855       }
       
   856       
       
   857     private:
       
   858       const ExtendFindEnum* extendFind;
       
   859       int cdx;
       
   860     };
       
   861 
       
   862     /// \brief Lemon style iterator for the items of a component.
       
   863     ///
       
   864     /// ClassIt is a lemon style iterator for the components. It iterates
       
   865     /// on the items of a class. By example if you want to iterate on
       
   866     /// each items of each classes then you may write the next code.
       
   867     ///\code
       
   868     /// for (ClassIt cit(ufe); cit != INVALID; ++cit) {
       
   869     ///   std::cout << "Class: ";
       
   870     ///   for (ItemIt iit(ufe, cit); iit != INVALID; ++iit) {
       
   871     ///     std::cout << toString(iit) << ' ' << std::endl;
       
   872     ///   }
       
   873     ///   std::cout << std::endl;
       
   874     /// }
       
   875     ///\endcode
       
   876     class ItemIt {
       
   877     public:
       
   878       /// \brief Constructor of the iterator
       
   879       ///
       
   880       /// Constructor of the iterator. The iterator iterates
       
   881       /// on the class of the \c item.
       
   882       ItemIt(const ExtendFindEnum& ufe, int cls) : extendFind(&ufe) {
       
   883         fdx = idx = extendFind->classes[cls].firstItem;
       
   884       }
       
   885 
       
   886       /// \brief Constructor to get invalid iterator
       
   887       ///
       
   888       /// Constructor to get invalid iterator
       
   889       ItemIt(Invalid) : extendFind(0), idx(-1) {}
       
   890       
       
   891       /// \brief Increment operator
       
   892       ///
       
   893       /// It steps to the next item in the class.
       
   894       ItemIt& operator++() {
       
   895         idx = extendFind->items[idx].next;
       
   896 	if (fdx == idx) idx = -1;
       
   897         return *this;
       
   898       }
       
   899       
       
   900       /// \brief Conversion operator
       
   901       ///
       
   902       /// It converts the iterator to the current item.
       
   903       operator const Item&() const {
       
   904         return extendFind->items[idx].item;
       
   905       }
       
   906 
       
   907       /// \brief Equality operator
       
   908       ///
       
   909       /// Equality operator
       
   910       bool operator==(const ItemIt& i) { 
       
   911         return i.idx == idx;
       
   912       }
       
   913 
       
   914       /// \brief Inequality operator
       
   915       ///
       
   916       /// Inequality operator
       
   917       bool operator!=(const ItemIt& i) { 
       
   918         return i.idx != idx;
       
   919       }
       
   920       
       
   921     private:
       
   922       const ExtendFindEnum* extendFind;
       
   923       int idx, fdx;
       
   924     };
       
   925 
       
   926   };
       
   927 
       
   928   /// \ingroup auxdat
       
   929   ///
       
   930   /// \brief A \e Union-Find data structure implementation which
       
   931   /// is able to store a priority for each item and retrieve the minimum of
       
   932   /// each class.
       
   933   ///
       
   934   /// A \e Union-Find data structure implementation which is able to
       
   935   /// store a priority for each item and retrieve the minimum of each
       
   936   /// class. In addition, it supports the joining and splitting the
       
   937   /// components. If you don't need this feature then you makes
       
   938   /// better to use the \ref UnionFind class which is more efficient.
       
   939   ///
       
   940   /// The union-find data strcuture based on a (2, 16)-tree with a
       
   941   /// tournament minimum selection on the internal nodes. The insert
       
   942   /// operation takes O(1), the find, set, decrease and increase takes
       
   943   /// O(log(n)), where n is the number of nodes in the current
       
   944   /// component.  The complexity of join and split is O(log(n)*k),
       
   945   /// where n is the sum of the number of the nodes and k is the
       
   946   /// number of joined components or the number of the components
       
   947   /// after the split.
       
   948   ///
       
   949   /// \pre You need to add all the elements by the \ref insert()
       
   950   /// method.
       
   951   ///
       
   952   template <typename _Value, typename _ItemIntMap, 
       
   953             typename _Comp = std::less<_Value> >
       
   954   class HeapUnionFind {
       
   955   public:
       
   956     
       
   957     typedef _Value Value;
       
   958     typedef typename _ItemIntMap::Key Item;
       
   959 
       
   960     typedef _ItemIntMap ItemIntMap;
       
   961 
       
   962     typedef _Comp Comp;
       
   963 
       
   964   private:
       
   965 
       
   966     static const int cmax = 16;
       
   967 
       
   968     ItemIntMap& index;
       
   969 
       
   970     struct ClassNode {
       
   971       int parent;
       
   972       int depth;
       
   973 
       
   974       int left, right;
       
   975       int next, prev;
       
   976     };
       
   977 
       
   978     int first_class;
       
   979     int first_free_class;
       
   980     std::vector<ClassNode> classes;
       
   981 
       
   982     int newClass() {
       
   983       if (first_free_class < 0) {
       
   984         int id = classes.size();
       
   985         classes.push_back(ClassNode());
       
   986         return id;
       
   987       } else {
       
   988         int id = first_free_class;
       
   989         first_free_class = classes[id].next;
       
   990         return id;
       
   991       }
       
   992     }
       
   993 
       
   994     void deleteClass(int id) {
       
   995       classes[id].next = first_free_class;
       
   996       first_free_class = id;
       
   997     }
       
   998 
       
   999     struct ItemNode {
       
  1000       int parent;
       
  1001       Item item;
       
  1002       Value prio;
       
  1003       int next, prev;
       
  1004       int left, right;
       
  1005       int size;
       
  1006     };
       
  1007 
       
  1008     int first_free_node;
       
  1009     std::vector<ItemNode> nodes;
       
  1010 
       
  1011     int newNode() {
       
  1012       if (first_free_node < 0) {
       
  1013         int id = nodes.size();
       
  1014         nodes.push_back(ItemNode());
       
  1015         return id;
       
  1016       } else {
       
  1017         int id = first_free_node;
       
  1018         first_free_node = nodes[id].next;
       
  1019         return id;
       
  1020       }
       
  1021     }
       
  1022 
       
  1023     void deleteNode(int id) {
       
  1024       nodes[id].next = first_free_node;
       
  1025       first_free_node = id;
       
  1026     }
       
  1027 
       
  1028     Comp comp;
       
  1029 
       
  1030     int findClass(int id) const {
       
  1031       int kd = id;
       
  1032       while (kd >= 0) {
       
  1033         kd = nodes[kd].parent;
       
  1034       }
       
  1035       return ~kd;
       
  1036     }
       
  1037 
       
  1038     int leftNode(int id) const {
       
  1039       int kd = ~(classes[id].parent);
       
  1040       for (int i = 0; i < classes[id].depth; ++i) {
       
  1041         kd = nodes[kd].left;
       
  1042       }
       
  1043       return kd;
       
  1044     }
       
  1045 
       
  1046     int nextNode(int id) const {
       
  1047       int depth = 0;
       
  1048       while (id >= 0 && nodes[id].next == -1) {
       
  1049         id = nodes[id].parent;
       
  1050         ++depth;
       
  1051       }
       
  1052       if (id < 0) {
       
  1053         return -1;
       
  1054       }
       
  1055       id = nodes[id].next;
       
  1056       while (depth--) {
       
  1057         id = nodes[id].left;      
       
  1058       }
       
  1059       return id;
       
  1060     }
       
  1061 
       
  1062 
       
  1063     void setPrio(int id) {
       
  1064       int jd = nodes[id].left;
       
  1065       nodes[id].prio = nodes[jd].prio;
       
  1066       nodes[id].item = nodes[jd].item;
       
  1067       jd = nodes[jd].next;
       
  1068       while (jd != -1) {
       
  1069         if (comp(nodes[jd].prio, nodes[id].prio)) {
       
  1070           nodes[id].prio = nodes[jd].prio;
       
  1071           nodes[id].item = nodes[jd].item;
       
  1072         }
       
  1073         jd = nodes[jd].next;
       
  1074       }
       
  1075     }
       
  1076 
       
  1077     void push(int id, int jd) {
       
  1078       nodes[id].size = 1;
       
  1079       nodes[id].left = nodes[id].right = jd;
       
  1080       nodes[jd].next = nodes[jd].prev = -1;
       
  1081       nodes[jd].parent = id;
       
  1082     }
       
  1083 
       
  1084     void pushAfter(int id, int jd) {
       
  1085       int kd = nodes[id].parent;
       
  1086       if (nodes[id].next != -1) {
       
  1087         nodes[nodes[id].next].prev = jd;
       
  1088         if (kd >= 0) {
       
  1089           nodes[kd].size += 1;
       
  1090         }
       
  1091       } else {
       
  1092         if (kd >= 0) {
       
  1093           nodes[kd].right = jd;
       
  1094           nodes[kd].size += 1;
       
  1095         }
       
  1096       }
       
  1097       nodes[jd].next = nodes[id].next;
       
  1098       nodes[jd].prev = id;
       
  1099       nodes[id].next = jd;
       
  1100       nodes[jd].parent = kd;
       
  1101     }
       
  1102 
       
  1103     void pushRight(int id, int jd) {
       
  1104       nodes[id].size += 1;
       
  1105       nodes[jd].prev = nodes[id].right;
       
  1106       nodes[jd].next = -1;
       
  1107       nodes[nodes[id].right].next = jd;
       
  1108       nodes[id].right = jd;
       
  1109       nodes[jd].parent = id;
       
  1110     }
       
  1111 
       
  1112     void popRight(int id) {
       
  1113       nodes[id].size -= 1;
       
  1114       int jd = nodes[id].right;
       
  1115       nodes[nodes[jd].prev].next = -1;
       
  1116       nodes[id].right = nodes[jd].prev;
       
  1117     }
       
  1118 
       
  1119     void splice(int id, int jd) {
       
  1120       nodes[id].size += nodes[jd].size;
       
  1121       nodes[nodes[id].right].next = nodes[jd].left;
       
  1122       nodes[nodes[jd].left].prev = nodes[id].right;
       
  1123       int kd = nodes[jd].left;
       
  1124       while (kd != -1) {
       
  1125         nodes[kd].parent = id;
       
  1126         kd = nodes[kd].next;
       
  1127       }
       
  1128       nodes[id].right = nodes[jd].right;
       
  1129     }
       
  1130 
       
  1131     void split(int id, int jd) {
       
  1132       int kd = nodes[id].parent;
       
  1133       nodes[kd].right = nodes[id].prev;
       
  1134       nodes[nodes[id].prev].next = -1;
       
  1135       
       
  1136       nodes[jd].left = id;
       
  1137       nodes[id].prev = -1;
       
  1138       int num = 0;
       
  1139       while (id != -1) {
       
  1140         nodes[id].parent = jd;
       
  1141         nodes[jd].right = id;
       
  1142         id = nodes[id].next;
       
  1143         ++num;
       
  1144       }      
       
  1145       nodes[kd].size -= num;
       
  1146       nodes[jd].size = num;
       
  1147     }
       
  1148 
       
  1149     void pushLeft(int id, int jd) {
       
  1150       nodes[id].size += 1;
       
  1151       nodes[jd].next = nodes[id].left;
       
  1152       nodes[jd].prev = -1;
       
  1153       nodes[nodes[id].left].prev = jd;
       
  1154       nodes[id].left = jd;
       
  1155       nodes[jd].parent = id;
       
  1156     }
       
  1157 
       
  1158     void popLeft(int id) {
       
  1159       nodes[id].size -= 1;
       
  1160       int jd = nodes[id].left;
       
  1161       nodes[nodes[jd].next].prev = -1;
       
  1162       nodes[id].left = nodes[jd].next;
       
  1163     }
       
  1164 
       
  1165     void repairLeft(int id) {
       
  1166       int jd = ~(classes[id].parent);
       
  1167       while (nodes[jd].left != -1) {
       
  1168 	int kd = nodes[jd].left;
       
  1169 	if (nodes[jd].size == 1) {
       
  1170 	  if (nodes[jd].parent < 0) {
       
  1171 	    classes[id].parent = ~kd;
       
  1172 	    classes[id].depth -= 1;
       
  1173 	    nodes[kd].parent = ~id;
       
  1174 	    deleteNode(jd);
       
  1175 	    jd = kd;
       
  1176 	  } else {
       
  1177 	    int pd = nodes[jd].parent;
       
  1178 	    if (nodes[nodes[jd].next].size < cmax) {
       
  1179 	      pushLeft(nodes[jd].next, nodes[jd].left);
       
  1180 	      if (less(nodes[jd].left, nodes[jd].next)) {
       
  1181 		nodes[nodes[jd].next].prio = nodes[nodes[jd].left].prio;
       
  1182 		nodes[nodes[jd].next].item = nodes[nodes[jd].left].item;
       
  1183 	      } 
       
  1184 	      popLeft(pd);
       
  1185 	      deleteNode(jd);
       
  1186 	      jd = pd;
       
  1187 	    } else {
       
  1188 	      int ld = nodes[nodes[jd].next].left;
       
  1189 	      popLeft(nodes[jd].next);
       
  1190 	      pushRight(jd, ld);
       
  1191 	      if (less(ld, nodes[jd].left)) {
       
  1192 		nodes[jd].item = nodes[ld].item;
       
  1193 		nodes[jd].prio = nodes[jd].prio;
       
  1194 	      }
       
  1195 	      if (nodes[nodes[jd].next].item == nodes[ld].item) {
       
  1196 		setPrio(nodes[jd].next);
       
  1197 	      }
       
  1198 	      jd = nodes[jd].left;
       
  1199 	    }
       
  1200 	  }
       
  1201 	} else {
       
  1202 	  jd = nodes[jd].left;
       
  1203 	}
       
  1204       }
       
  1205     }    
       
  1206 
       
  1207     void repairRight(int id) {
       
  1208       int jd = ~(classes[id].parent);
       
  1209       while (nodes[jd].right != -1) {
       
  1210 	int kd = nodes[jd].right;
       
  1211 	if (nodes[jd].size == 1) {
       
  1212 	  if (nodes[jd].parent < 0) {
       
  1213 	    classes[id].parent = ~kd;
       
  1214 	    classes[id].depth -= 1;
       
  1215 	    nodes[kd].parent = ~id;
       
  1216 	    deleteNode(jd);
       
  1217 	    jd = kd;
       
  1218 	  } else {
       
  1219 	    int pd = nodes[jd].parent;
       
  1220 	    if (nodes[nodes[jd].prev].size < cmax) {
       
  1221 	      pushRight(nodes[jd].prev, nodes[jd].right);
       
  1222 	      if (less(nodes[jd].right, nodes[jd].prev)) {
       
  1223 		nodes[nodes[jd].prev].prio = nodes[nodes[jd].right].prio;
       
  1224 		nodes[nodes[jd].prev].item = nodes[nodes[jd].right].item;
       
  1225 	      } 
       
  1226 	      popRight(pd);
       
  1227 	      deleteNode(jd);
       
  1228 	      jd = pd;
       
  1229 	    } else {
       
  1230 	      int ld = nodes[nodes[jd].prev].right;
       
  1231 	      popRight(nodes[jd].prev);
       
  1232 	      pushLeft(jd, ld);
       
  1233 	      if (less(ld, nodes[jd].right)) {
       
  1234 		nodes[jd].item = nodes[ld].item;
       
  1235 		nodes[jd].prio = nodes[jd].prio;
       
  1236 	      }
       
  1237 	      if (nodes[nodes[jd].prev].item == nodes[ld].item) {
       
  1238 		setPrio(nodes[jd].prev);
       
  1239 	      }
       
  1240 	      jd = nodes[jd].right;
       
  1241 	    }
       
  1242 	  }
       
  1243 	} else {
       
  1244 	  jd = nodes[jd].right;
       
  1245 	}
       
  1246       }
       
  1247     }
       
  1248 
       
  1249 
       
  1250     bool less(int id, int jd) const {
       
  1251       return comp(nodes[id].prio, nodes[jd].prio);
       
  1252     }
       
  1253 
       
  1254     bool equal(int id, int jd) const {
       
  1255       return !less(id, jd) && !less(jd, id);
       
  1256     }
       
  1257 
       
  1258 
       
  1259   public:
       
  1260 
       
  1261     /// \brief Returns true when the given class is alive.
       
  1262     ///
       
  1263     /// Returns true when the given class is alive, ie. the class is
       
  1264     /// not nested into other class.
       
  1265     bool alive(int cls) const {
       
  1266       return classes[cls].parent < 0;
       
  1267     }
       
  1268 
       
  1269     /// \brief Returns true when the given class is trivial.
       
  1270     ///
       
  1271     /// Returns true when the given class is trivial, ie. the class
       
  1272     /// contains just one item directly.
       
  1273     bool trivial(int cls) const {
       
  1274       return classes[cls].left == -1;
       
  1275     }
       
  1276 
       
  1277     /// \brief Constructs the union-find.
       
  1278     ///
       
  1279     /// Constructs the union-find.  
       
  1280     /// \brief _index The index map of the union-find. The data
       
  1281     /// structure uses internally for store references.
       
  1282     HeapUnionFind(ItemIntMap& _index) 
       
  1283       : index(_index), first_class(-1), 
       
  1284 	first_free_class(-1), first_free_node(-1) {}
       
  1285 
       
  1286     /// \brief Insert a new node into a new component.
       
  1287     ///
       
  1288     /// Insert a new node into a new component.
       
  1289     /// \param item The item of the new node.
       
  1290     /// \param prio The priority of the new node.
       
  1291     /// \return The class id of the one-item-heap.
       
  1292     int insert(const Item& item, const Value& prio) {
       
  1293       int id = newNode();
       
  1294       nodes[id].item = item;
       
  1295       nodes[id].prio = prio;
       
  1296       nodes[id].size = 0;
       
  1297 
       
  1298       nodes[id].prev = -1;
       
  1299       nodes[id].next = -1;
       
  1300 
       
  1301       nodes[id].left = -1;
       
  1302       nodes[id].right = -1;
       
  1303 
       
  1304       nodes[id].item = item;
       
  1305       index[item] = id;
       
  1306       
       
  1307       int class_id = newClass();
       
  1308       classes[class_id].parent = ~id;
       
  1309       classes[class_id].depth = 0;
       
  1310 
       
  1311       classes[class_id].left = -1;
       
  1312       classes[class_id].right = -1;
       
  1313       
       
  1314       if (first_class != -1) {
       
  1315         classes[first_class].prev = class_id;
       
  1316       }
       
  1317       classes[class_id].next = first_class;
       
  1318       classes[class_id].prev = -1;
       
  1319       first_class = class_id;
       
  1320 
       
  1321       nodes[id].parent = ~class_id;
       
  1322       
       
  1323       return class_id;
       
  1324     }
       
  1325 
       
  1326     /// \brief The class of the item.
       
  1327     ///
       
  1328     /// \return The alive class id of the item, which is not nested into
       
  1329     /// other classes.
       
  1330     ///
       
  1331     /// The time complexity is O(log(n)).
       
  1332     int find(const Item& item) const {
       
  1333       return findClass(index[item]);
       
  1334     }
       
  1335     
       
  1336     /// \brief Joins the classes.
       
  1337     ///
       
  1338     /// The current function joins the given classes. The parameter is
       
  1339     /// an STL range which should be contains valid class ids. The
       
  1340     /// time complexity is O(log(n)*k) where n is the overall number
       
  1341     /// of the joined nodes and k is the number of classes.
       
  1342     /// \return The class of the joined classes.
       
  1343     /// \pre The range should contain at least two class ids.
       
  1344     template <typename Iterator>
       
  1345     int join(Iterator begin, Iterator end) {
       
  1346       std::vector<int> cs;
       
  1347       for (Iterator it = begin; it != end; ++it) {
       
  1348         cs.push_back(*it);
       
  1349       }
       
  1350 
       
  1351       int class_id = newClass();
       
  1352       { // creation union-find
       
  1353 
       
  1354         if (first_class != -1) {
       
  1355           classes[first_class].prev = class_id;
       
  1356         }
       
  1357         classes[class_id].next = first_class;
       
  1358         classes[class_id].prev = -1;
       
  1359         first_class = class_id;
       
  1360 
       
  1361         classes[class_id].depth = classes[cs[0]].depth;
       
  1362         classes[class_id].parent = classes[cs[0]].parent;
       
  1363         nodes[~(classes[class_id].parent)].parent = ~class_id;
       
  1364 
       
  1365         int l = cs[0];
       
  1366 
       
  1367         classes[class_id].left = l;
       
  1368         classes[class_id].right = l;
       
  1369 
       
  1370         if (classes[l].next != -1) {
       
  1371           classes[classes[l].next].prev = classes[l].prev;
       
  1372         }
       
  1373         classes[classes[l].prev].next = classes[l].next;
       
  1374         
       
  1375         classes[l].prev = -1;
       
  1376         classes[l].next = -1;
       
  1377 
       
  1378         classes[l].depth = leftNode(l);
       
  1379         classes[l].parent = class_id;
       
  1380         
       
  1381       }
       
  1382 
       
  1383       { // merging of heap
       
  1384         int l = class_id;
       
  1385         for (int ci = 1; ci < int(cs.size()); ++ci) {
       
  1386           int r = cs[ci];
       
  1387           int rln = leftNode(r);
       
  1388           if (classes[l].depth > classes[r].depth) {
       
  1389             int id = ~(classes[l].parent);
       
  1390             for (int i = classes[r].depth + 1; i < classes[l].depth; ++i) {
       
  1391               id = nodes[id].right;
       
  1392             }
       
  1393             while (id >= 0 && nodes[id].size == cmax) {
       
  1394               int new_id = newNode();
       
  1395               int right_id = nodes[id].right;
       
  1396 
       
  1397               popRight(id);
       
  1398               if (nodes[id].item == nodes[right_id].item) {
       
  1399                 setPrio(id);
       
  1400               }
       
  1401               push(new_id, right_id);
       
  1402               pushRight(new_id, ~(classes[r].parent));
       
  1403               setPrio(new_id);
       
  1404 
       
  1405               id = nodes[id].parent;
       
  1406               classes[r].parent = ~new_id;
       
  1407             }
       
  1408             if (id < 0) {
       
  1409               int new_parent = newNode();
       
  1410               nodes[new_parent].next = -1;
       
  1411               nodes[new_parent].prev = -1;
       
  1412               nodes[new_parent].parent = ~l;
       
  1413 
       
  1414               push(new_parent, ~(classes[l].parent));
       
  1415               pushRight(new_parent, ~(classes[r].parent));
       
  1416               setPrio(new_parent);
       
  1417 
       
  1418               classes[l].parent = ~new_parent;
       
  1419               classes[l].depth += 1;
       
  1420             } else {
       
  1421               pushRight(id, ~(classes[r].parent));
       
  1422               while (id >= 0 && less(~(classes[r].parent), id)) {
       
  1423                 nodes[id].prio = nodes[~(classes[r].parent)].prio;
       
  1424                 nodes[id].item = nodes[~(classes[r].parent)].item;
       
  1425                 id = nodes[id].parent;
       
  1426               }
       
  1427             }
       
  1428           } else if (classes[r].depth > classes[l].depth) {
       
  1429             int id = ~(classes[r].parent);
       
  1430             for (int i = classes[l].depth + 1; i < classes[r].depth; ++i) {
       
  1431               id = nodes[id].left;
       
  1432             }
       
  1433             while (id >= 0 && nodes[id].size == cmax) {
       
  1434               int new_id = newNode();
       
  1435               int left_id = nodes[id].left;
       
  1436 
       
  1437               popLeft(id);
       
  1438               if (nodes[id].prio == nodes[left_id].prio) {
       
  1439                 setPrio(id);
       
  1440               }
       
  1441               push(new_id, left_id);
       
  1442               pushLeft(new_id, ~(classes[l].parent));
       
  1443               setPrio(new_id);
       
  1444 
       
  1445               id = nodes[id].parent;
       
  1446               classes[l].parent = ~new_id;
       
  1447 
       
  1448             }
       
  1449             if (id < 0) {
       
  1450               int new_parent = newNode();
       
  1451               nodes[new_parent].next = -1;
       
  1452               nodes[new_parent].prev = -1;
       
  1453               nodes[new_parent].parent = ~l;
       
  1454 
       
  1455               push(new_parent, ~(classes[r].parent));
       
  1456               pushLeft(new_parent, ~(classes[l].parent));
       
  1457               setPrio(new_parent);
       
  1458               
       
  1459               classes[r].parent = ~new_parent;
       
  1460               classes[r].depth += 1;
       
  1461             } else {
       
  1462               pushLeft(id, ~(classes[l].parent));
       
  1463               while (id >= 0 && less(~(classes[l].parent), id)) {
       
  1464                 nodes[id].prio = nodes[~(classes[l].parent)].prio;
       
  1465                 nodes[id].item = nodes[~(classes[l].parent)].item;
       
  1466                 id = nodes[id].parent;
       
  1467               }
       
  1468             }
       
  1469             nodes[~(classes[r].parent)].parent = ~l;
       
  1470             classes[l].parent = classes[r].parent;
       
  1471             classes[l].depth = classes[r].depth;
       
  1472           } else {
       
  1473             if (classes[l].depth != 0 && 
       
  1474                 nodes[~(classes[l].parent)].size + 
       
  1475                 nodes[~(classes[r].parent)].size <= cmax) {
       
  1476               splice(~(classes[l].parent), ~(classes[r].parent));
       
  1477               deleteNode(~(classes[r].parent));
       
  1478               if (less(~(classes[r].parent), ~(classes[l].parent))) {
       
  1479                 nodes[~(classes[l].parent)].prio = 
       
  1480                   nodes[~(classes[r].parent)].prio;
       
  1481                 nodes[~(classes[l].parent)].item = 
       
  1482                   nodes[~(classes[r].parent)].item;
       
  1483               }
       
  1484             } else {
       
  1485               int new_parent = newNode();
       
  1486               nodes[new_parent].next = nodes[new_parent].prev = -1;
       
  1487               push(new_parent, ~(classes[l].parent));
       
  1488               pushRight(new_parent, ~(classes[r].parent));
       
  1489               setPrio(new_parent);
       
  1490             
       
  1491               classes[l].parent = ~new_parent;
       
  1492               classes[l].depth += 1;
       
  1493               nodes[new_parent].parent = ~l;
       
  1494             }
       
  1495           }
       
  1496           if (classes[r].next != -1) {
       
  1497             classes[classes[r].next].prev = classes[r].prev;
       
  1498           }
       
  1499           classes[classes[r].prev].next = classes[r].next;
       
  1500 
       
  1501           classes[r].prev = classes[l].right;
       
  1502           classes[classes[l].right].next = r;
       
  1503           classes[l].right = r;
       
  1504           classes[r].parent = l;
       
  1505 
       
  1506           classes[r].next = -1;
       
  1507           classes[r].depth = rln;
       
  1508         }
       
  1509       }
       
  1510       return class_id;
       
  1511     }
       
  1512 
       
  1513     /// \brief Split the class to subclasses.
       
  1514     ///
       
  1515     /// The current function splits the given class. The join, which
       
  1516     /// made the current class, stored a reference to the
       
  1517     /// subclasses. The \c splitClass() member restores the classes
       
  1518     /// and creates the heaps. The parameter is an STL output iterator
       
  1519     /// which will be filled with the subclass ids. The time
       
  1520     /// complexity is O(log(n)*k) where n is the overall number of
       
  1521     /// nodes in the splitted classes and k is the number of the
       
  1522     /// classes.
       
  1523     template <typename Iterator>
       
  1524     void split(int cls, Iterator out) {
       
  1525       std::vector<int> cs;
       
  1526       { // splitting union-find
       
  1527         int id = cls;
       
  1528         int l = classes[id].left;
       
  1529 
       
  1530         classes[l].parent = classes[id].parent;
       
  1531         classes[l].depth = classes[id].depth;
       
  1532 
       
  1533         nodes[~(classes[l].parent)].parent = ~l;
       
  1534 
       
  1535         *out++ = l;
       
  1536 
       
  1537         while (l != -1) {
       
  1538           cs.push_back(l);
       
  1539           l = classes[l].next;
       
  1540         }
       
  1541 
       
  1542         classes[classes[id].right].next = first_class;
       
  1543         classes[first_class].prev = classes[id].right;
       
  1544         first_class = classes[id].left;
       
  1545         
       
  1546         if (classes[id].next != -1) {
       
  1547           classes[classes[id].next].prev = classes[id].prev;
       
  1548         }
       
  1549         classes[classes[id].prev].next = classes[id].next;
       
  1550         
       
  1551         deleteClass(id);
       
  1552       }
       
  1553 
       
  1554       {
       
  1555         for (int i = 1; i < int(cs.size()); ++i) {
       
  1556           int l = classes[cs[i]].depth;
       
  1557           while (nodes[nodes[l].parent].left == l) {
       
  1558             l = nodes[l].parent;
       
  1559           }
       
  1560           int r = l;    
       
  1561           while (nodes[l].parent >= 0) {
       
  1562             l = nodes[l].parent;
       
  1563             int new_node = newNode();
       
  1564 
       
  1565             nodes[new_node].prev = -1;
       
  1566             nodes[new_node].next = -1;
       
  1567 
       
  1568             split(r, new_node);
       
  1569             pushAfter(l, new_node);
       
  1570             setPrio(l);
       
  1571             setPrio(new_node);
       
  1572             r = new_node;
       
  1573           }
       
  1574           classes[cs[i]].parent = ~r;
       
  1575           classes[cs[i]].depth = classes[~(nodes[l].parent)].depth;
       
  1576           nodes[r].parent = ~cs[i];
       
  1577 
       
  1578           nodes[l].next = -1;
       
  1579           nodes[r].prev = -1;
       
  1580 
       
  1581           repairRight(~(nodes[l].parent));
       
  1582           repairLeft(cs[i]);
       
  1583           
       
  1584           *out++ = cs[i];
       
  1585         }
       
  1586       }
       
  1587     }
       
  1588 
       
  1589     /// \brief Gives back the priority of the current item.
       
  1590     ///
       
  1591     /// \return Gives back the priority of the current item.
       
  1592     const Value& operator[](const Item& item) const {
       
  1593       return nodes[index[item]].prio;
       
  1594     }
       
  1595 
       
  1596     /// \brief Sets the priority of the current item.
       
  1597     ///
       
  1598     /// Sets the priority of the current item.
       
  1599     void set(const Item& item, const Value& prio) {
       
  1600       if (comp(prio, nodes[index[item]].prio)) {
       
  1601         decrease(item, prio);
       
  1602       } else if (!comp(prio, nodes[index[item]].prio)) {
       
  1603         increase(item, prio);
       
  1604       }
       
  1605     }
       
  1606       
       
  1607     /// \brief Increase the priority of the current item.
       
  1608     ///
       
  1609     /// Increase the priority of the current item.
       
  1610     void increase(const Item& item, const Value& prio) {
       
  1611       int id = index[item];
       
  1612       int kd = nodes[id].parent;
       
  1613       nodes[id].prio = prio;
       
  1614       while (kd >= 0 && nodes[kd].item == item) {
       
  1615         setPrio(kd);
       
  1616         kd = nodes[kd].parent;
       
  1617       }
       
  1618     }
       
  1619 
       
  1620     /// \brief Increase the priority of the current item.
       
  1621     ///
       
  1622     /// Increase the priority of the current item.
       
  1623     void decrease(const Item& item, const Value& prio) {
       
  1624       int id = index[item];
       
  1625       int kd = nodes[id].parent;
       
  1626       nodes[id].prio = prio;
       
  1627       while (kd >= 0 && less(id, kd)) {
       
  1628         nodes[kd].prio = prio;
       
  1629         nodes[kd].item = item;
       
  1630         kd = nodes[kd].parent;
       
  1631       }
       
  1632     }
       
  1633     
       
  1634     /// \brief Gives back the minimum priority of the class.
       
  1635     ///
       
  1636     /// \return Gives back the minimum priority of the class.
       
  1637     const Value& classPrio(int cls) const {
       
  1638       return nodes[~(classes[cls].parent)].prio;
       
  1639     }
       
  1640 
       
  1641     /// \brief Gives back the minimum priority item of the class.
       
  1642     ///
       
  1643     /// \return Gives back the minimum priority item of the class.
       
  1644     const Item& classTop(int cls) const {
       
  1645       return nodes[~(classes[cls].parent)].item;
       
  1646     }
       
  1647 
       
  1648     /// \brief Gives back a representant item of the class.
       
  1649     /// 
       
  1650     /// The representant is indpendent from the priorities of the
       
  1651     /// items. 
       
  1652     /// \return Gives back a representant item of the class.
       
  1653     const Item& classRep(int id) const {
       
  1654       int parent = classes[id].parent;
       
  1655       return nodes[parent >= 0 ? classes[id].depth : leftNode(id)].item;
       
  1656     }
       
  1657 
       
  1658     /// \brief Lemon style iterator for the items of a class.
       
  1659     ///
       
  1660     /// ClassIt is a lemon style iterator for the components. It iterates
       
  1661     /// on the items of a class. By example if you want to iterate on
       
  1662     /// each items of each classes then you may write the next code.
       
  1663     ///\code
       
  1664     /// for (ClassIt cit(huf); cit != INVALID; ++cit) {
       
  1665     ///   std::cout << "Class: ";
       
  1666     ///   for (ItemIt iit(huf, cit); iit != INVALID; ++iit) {
       
  1667     ///     std::cout << toString(iit) << ' ' << std::endl;
       
  1668     ///   }
       
  1669     ///   std::cout << std::endl;
       
  1670     /// }
       
  1671     ///\endcode
       
  1672     class ItemIt {
       
  1673     private:
       
  1674 
       
  1675       const HeapUnionFind* _huf;
       
  1676       int _id, _lid;
       
  1677       
       
  1678     public:
       
  1679 
       
  1680       /// \brief Default constructor 
       
  1681       ///
       
  1682       /// Default constructor 
       
  1683       ItemIt() {}
       
  1684 
       
  1685       ItemIt(const HeapUnionFind& huf, int cls) : _huf(&huf) {
       
  1686         int id = cls;
       
  1687         int parent = _huf->classes[id].parent;
       
  1688         if (parent >= 0) {
       
  1689           _id = _huf->classes[id].depth;
       
  1690           if (_huf->classes[id].next != -1) {
       
  1691             _lid = _huf->classes[_huf->classes[id].next].depth;
       
  1692           } else {
       
  1693             _lid = -1;
       
  1694           }
       
  1695         } else {
       
  1696           _id = _huf->leftNode(id);
       
  1697           _lid = -1;
       
  1698         } 
       
  1699       }
       
  1700       
       
  1701       /// \brief Increment operator
       
  1702       ///
       
  1703       /// It steps to the next item in the class.
       
  1704       ItemIt& operator++() {
       
  1705         _id = _huf->nextNode(_id);
       
  1706         return *this;
       
  1707       }
       
  1708 
       
  1709       /// \brief Conversion operator
       
  1710       ///
       
  1711       /// It converts the iterator to the current item.
       
  1712       operator const Item&() const {
       
  1713         return _huf->nodes[_id].item;
       
  1714       }
       
  1715       
       
  1716       /// \brief Equality operator
       
  1717       ///
       
  1718       /// Equality operator
       
  1719       bool operator==(const ItemIt& i) { 
       
  1720         return i._id == _id;
       
  1721       }
       
  1722 
       
  1723       /// \brief Inequality operator
       
  1724       ///
       
  1725       /// Inequality operator
       
  1726       bool operator!=(const ItemIt& i) { 
       
  1727         return i._id != _id;
       
  1728       }
       
  1729 
       
  1730       /// \brief Equality operator
       
  1731       ///
       
  1732       /// Equality operator
       
  1733       bool operator==(Invalid) { 
       
  1734         return _id == _lid;
       
  1735       }
       
  1736 
       
  1737       /// \brief Inequality operator
       
  1738       ///
       
  1739       /// Inequality operator
       
  1740       bool operator!=(Invalid) { 
       
  1741         return _id != _lid;
       
  1742       }
       
  1743       
       
  1744     };
       
  1745 
       
  1746     /// \brief Class iterator
       
  1747     ///
       
  1748     /// The iterator stores 
       
  1749     class ClassIt {
       
  1750     private:
       
  1751 
       
  1752       const HeapUnionFind* _huf;
       
  1753       int _id;
       
  1754 
       
  1755     public:
       
  1756 
       
  1757       ClassIt(const HeapUnionFind& huf) 
       
  1758         : _huf(&huf), _id(huf.first_class) {}
       
  1759 
       
  1760       ClassIt(const HeapUnionFind& huf, int cls) 
       
  1761         : _huf(&huf), _id(huf.classes[cls].left) {}
       
  1762 
       
  1763       ClassIt(Invalid) : _huf(0), _id(-1) {}
       
  1764       
       
  1765       const ClassIt& operator++() {
       
  1766         _id = _huf->classes[_id].next;
       
  1767 	return *this;
       
  1768       }
       
  1769 
       
  1770       /// \brief Equality operator
       
  1771       ///
       
  1772       /// Equality operator
       
  1773       bool operator==(const ClassIt& i) { 
       
  1774         return i._id == _id;
       
  1775       }
       
  1776 
       
  1777       /// \brief Inequality operator
       
  1778       ///
       
  1779       /// Inequality operator
       
  1780       bool operator!=(const ClassIt& i) { 
       
  1781         return i._id != _id;
       
  1782       }      
       
  1783       
       
  1784       operator int() const {
       
  1785 	return _id;
       
  1786       }
       
  1787             
       
  1788     };
       
  1789 
       
  1790   };
       
  1791 
       
  1792   //! @}
       
  1793 
       
  1794 } //namespace lemon
       
  1795 
       
  1796 #endif //LEMON_UNION_FIND_H