lemon/unionfind.h
author Balazs Dezso <deba@inf.elte.hu>
Sun, 18 Jan 2009 17:49:08 +0100
changeset 471 fb5af0411793
parent 438 81d40f1c850c
child 550 c5fd2d996909
permissions -rw-r--r--
Fix lp indexing bug (#205)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     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/core.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(jd, nodes[jd].next) ||
  1181                   nodes[jd].item == nodes[pd].item) {
  1182                 nodes[nodes[jd].next].prio = nodes[jd].prio;
  1183                 nodes[nodes[jd].next].item = nodes[jd].item;
  1184               }
  1185               popLeft(pd);
  1186               deleteNode(jd);
  1187               jd = pd;
  1188             } else {
  1189               int ld = nodes[nodes[jd].next].left;
  1190               popLeft(nodes[jd].next);
  1191               pushRight(jd, ld);
  1192               if (less(ld, nodes[jd].left) ||
  1193                   nodes[ld].item == nodes[pd].item) {
  1194                 nodes[jd].item = nodes[ld].item;
  1195                 nodes[jd].prio = nodes[ld].prio;
  1196               }
  1197               if (nodes[nodes[jd].next].item == nodes[ld].item) {
  1198                 setPrio(nodes[jd].next);
  1199               }
  1200               jd = nodes[jd].left;
  1201             }
  1202           }
  1203         } else {
  1204           jd = nodes[jd].left;
  1205         }
  1206       }
  1207     }
  1208 
  1209     void repairRight(int id) {
  1210       int jd = ~(classes[id].parent);
  1211       while (nodes[jd].right != -1) {
  1212         int kd = nodes[jd].right;
  1213         if (nodes[jd].size == 1) {
  1214           if (nodes[jd].parent < 0) {
  1215             classes[id].parent = ~kd;
  1216             classes[id].depth -= 1;
  1217             nodes[kd].parent = ~id;
  1218             deleteNode(jd);
  1219             jd = kd;
  1220           } else {
  1221             int pd = nodes[jd].parent;
  1222             if (nodes[nodes[jd].prev].size < cmax) {
  1223               pushRight(nodes[jd].prev, nodes[jd].right);
  1224               if (less(jd, nodes[jd].prev) ||
  1225                   nodes[jd].item == nodes[pd].item) {
  1226                 nodes[nodes[jd].prev].prio = nodes[jd].prio;
  1227                 nodes[nodes[jd].prev].item = nodes[jd].item;
  1228               }
  1229               popRight(pd);
  1230               deleteNode(jd);
  1231               jd = pd;
  1232             } else {
  1233               int ld = nodes[nodes[jd].prev].right;
  1234               popRight(nodes[jd].prev);
  1235               pushLeft(jd, ld);
  1236               if (less(ld, nodes[jd].right) ||
  1237                   nodes[ld].item == nodes[pd].item) {
  1238                 nodes[jd].item = nodes[ld].item;
  1239                 nodes[jd].prio = nodes[ld].prio;
  1240               }
  1241               if (nodes[nodes[jd].prev].item == nodes[ld].item) {
  1242                 setPrio(nodes[jd].prev);
  1243               }
  1244               jd = nodes[jd].right;
  1245             }
  1246           }
  1247         } else {
  1248           jd = nodes[jd].right;
  1249         }
  1250       }
  1251     }
  1252 
  1253 
  1254     bool less(int id, int jd) const {
  1255       return comp(nodes[id].prio, nodes[jd].prio);
  1256     }
  1257 
  1258   public:
  1259 
  1260     /// \brief Returns true when the given class is alive.
  1261     ///
  1262     /// Returns true when the given class is alive, ie. the class is
  1263     /// not nested into other class.
  1264     bool alive(int cls) const {
  1265       return classes[cls].parent < 0;
  1266     }
  1267 
  1268     /// \brief Returns true when the given class is trivial.
  1269     ///
  1270     /// Returns true when the given class is trivial, ie. the class
  1271     /// contains just one item directly.
  1272     bool trivial(int cls) const {
  1273       return classes[cls].left == -1;
  1274     }
  1275 
  1276     /// \brief Constructs the union-find.
  1277     ///
  1278     /// Constructs the union-find.
  1279     /// \brief _index The index map of the union-find. The data
  1280     /// structure uses internally for store references.
  1281     HeapUnionFind(ItemIntMap& _index)
  1282       : index(_index), first_class(-1),
  1283         first_free_class(-1), first_free_node(-1) {}
  1284 
  1285     /// \brief Insert a new node into a new component.
  1286     ///
  1287     /// Insert a new node into a new component.
  1288     /// \param item The item of the new node.
  1289     /// \param prio The priority of the new node.
  1290     /// \return The class id of the one-item-heap.
  1291     int insert(const Item& item, const Value& prio) {
  1292       int id = newNode();
  1293       nodes[id].item = item;
  1294       nodes[id].prio = prio;
  1295       nodes[id].size = 0;
  1296 
  1297       nodes[id].prev = -1;
  1298       nodes[id].next = -1;
  1299 
  1300       nodes[id].left = -1;
  1301       nodes[id].right = -1;
  1302 
  1303       nodes[id].item = item;
  1304       index[item] = id;
  1305 
  1306       int class_id = newClass();
  1307       classes[class_id].parent = ~id;
  1308       classes[class_id].depth = 0;
  1309 
  1310       classes[class_id].left = -1;
  1311       classes[class_id].right = -1;
  1312 
  1313       if (first_class != -1) {
  1314         classes[first_class].prev = class_id;
  1315       }
  1316       classes[class_id].next = first_class;
  1317       classes[class_id].prev = -1;
  1318       first_class = class_id;
  1319 
  1320       nodes[id].parent = ~class_id;
  1321 
  1322       return class_id;
  1323     }
  1324 
  1325     /// \brief The class of the item.
  1326     ///
  1327     /// \return The alive class id of the item, which is not nested into
  1328     /// other classes.
  1329     ///
  1330     /// The time complexity is O(log(n)).
  1331     int find(const Item& item) const {
  1332       return findClass(index[item]);
  1333     }
  1334 
  1335     /// \brief Joins the classes.
  1336     ///
  1337     /// The current function joins the given classes. The parameter is
  1338     /// an STL range which should be contains valid class ids. The
  1339     /// time complexity is O(log(n)*k) where n is the overall number
  1340     /// of the joined nodes and k is the number of classes.
  1341     /// \return The class of the joined classes.
  1342     /// \pre The range should contain at least two class ids.
  1343     template <typename Iterator>
  1344     int join(Iterator begin, Iterator end) {
  1345       std::vector<int> cs;
  1346       for (Iterator it = begin; it != end; ++it) {
  1347         cs.push_back(*it);
  1348       }
  1349 
  1350       int class_id = newClass();
  1351       { // creation union-find
  1352 
  1353         if (first_class != -1) {
  1354           classes[first_class].prev = class_id;
  1355         }
  1356         classes[class_id].next = first_class;
  1357         classes[class_id].prev = -1;
  1358         first_class = class_id;
  1359 
  1360         classes[class_id].depth = classes[cs[0]].depth;
  1361         classes[class_id].parent = classes[cs[0]].parent;
  1362         nodes[~(classes[class_id].parent)].parent = ~class_id;
  1363 
  1364         int l = cs[0];
  1365 
  1366         classes[class_id].left = l;
  1367         classes[class_id].right = l;
  1368 
  1369         if (classes[l].next != -1) {
  1370           classes[classes[l].next].prev = classes[l].prev;
  1371         }
  1372         classes[classes[l].prev].next = classes[l].next;
  1373 
  1374         classes[l].prev = -1;
  1375         classes[l].next = -1;
  1376 
  1377         classes[l].depth = leftNode(l);
  1378         classes[l].parent = class_id;
  1379 
  1380       }
  1381 
  1382       { // merging of heap
  1383         int l = class_id;
  1384         for (int ci = 1; ci < int(cs.size()); ++ci) {
  1385           int r = cs[ci];
  1386           int rln = leftNode(r);
  1387           if (classes[l].depth > classes[r].depth) {
  1388             int id = ~(classes[l].parent);
  1389             for (int i = classes[r].depth + 1; i < classes[l].depth; ++i) {
  1390               id = nodes[id].right;
  1391             }
  1392             while (id >= 0 && nodes[id].size == cmax) {
  1393               int new_id = newNode();
  1394               int right_id = nodes[id].right;
  1395 
  1396               popRight(id);
  1397               if (nodes[id].item == nodes[right_id].item) {
  1398                 setPrio(id);
  1399               }
  1400               push(new_id, right_id);
  1401               pushRight(new_id, ~(classes[r].parent));
  1402 
  1403               if (less(~classes[r].parent, right_id)) {
  1404                 nodes[new_id].item = nodes[~classes[r].parent].item;
  1405                 nodes[new_id].prio = nodes[~classes[r].parent].prio;
  1406               } else {
  1407                 nodes[new_id].item = nodes[right_id].item;
  1408                 nodes[new_id].prio = nodes[right_id].prio;
  1409               }
  1410 
  1411               id = nodes[id].parent;
  1412               classes[r].parent = ~new_id;
  1413             }
  1414             if (id < 0) {
  1415               int new_parent = newNode();
  1416               nodes[new_parent].next = -1;
  1417               nodes[new_parent].prev = -1;
  1418               nodes[new_parent].parent = ~l;
  1419 
  1420               push(new_parent, ~(classes[l].parent));
  1421               pushRight(new_parent, ~(classes[r].parent));
  1422               setPrio(new_parent);
  1423 
  1424               classes[l].parent = ~new_parent;
  1425               classes[l].depth += 1;
  1426             } else {
  1427               pushRight(id, ~(classes[r].parent));
  1428               while (id >= 0 && less(~(classes[r].parent), id)) {
  1429                 nodes[id].prio = nodes[~(classes[r].parent)].prio;
  1430                 nodes[id].item = nodes[~(classes[r].parent)].item;
  1431                 id = nodes[id].parent;
  1432               }
  1433             }
  1434           } else if (classes[r].depth > classes[l].depth) {
  1435             int id = ~(classes[r].parent);
  1436             for (int i = classes[l].depth + 1; i < classes[r].depth; ++i) {
  1437               id = nodes[id].left;
  1438             }
  1439             while (id >= 0 && nodes[id].size == cmax) {
  1440               int new_id = newNode();
  1441               int left_id = nodes[id].left;
  1442 
  1443               popLeft(id);
  1444               if (nodes[id].prio == nodes[left_id].prio) {
  1445                 setPrio(id);
  1446               }
  1447               push(new_id, left_id);
  1448               pushLeft(new_id, ~(classes[l].parent));
  1449 
  1450               if (less(~classes[l].parent, left_id)) {
  1451                 nodes[new_id].item = nodes[~classes[l].parent].item;
  1452                 nodes[new_id].prio = nodes[~classes[l].parent].prio;
  1453               } else {
  1454                 nodes[new_id].item = nodes[left_id].item;
  1455                 nodes[new_id].prio = nodes[left_id].prio;
  1456               }
  1457 
  1458               id = nodes[id].parent;
  1459               classes[l].parent = ~new_id;
  1460 
  1461             }
  1462             if (id < 0) {
  1463               int new_parent = newNode();
  1464               nodes[new_parent].next = -1;
  1465               nodes[new_parent].prev = -1;
  1466               nodes[new_parent].parent = ~l;
  1467 
  1468               push(new_parent, ~(classes[r].parent));
  1469               pushLeft(new_parent, ~(classes[l].parent));
  1470               setPrio(new_parent);
  1471 
  1472               classes[r].parent = ~new_parent;
  1473               classes[r].depth += 1;
  1474             } else {
  1475               pushLeft(id, ~(classes[l].parent));
  1476               while (id >= 0 && less(~(classes[l].parent), id)) {
  1477                 nodes[id].prio = nodes[~(classes[l].parent)].prio;
  1478                 nodes[id].item = nodes[~(classes[l].parent)].item;
  1479                 id = nodes[id].parent;
  1480               }
  1481             }
  1482             nodes[~(classes[r].parent)].parent = ~l;
  1483             classes[l].parent = classes[r].parent;
  1484             classes[l].depth = classes[r].depth;
  1485           } else {
  1486             if (classes[l].depth != 0 &&
  1487                 nodes[~(classes[l].parent)].size +
  1488                 nodes[~(classes[r].parent)].size <= cmax) {
  1489               splice(~(classes[l].parent), ~(classes[r].parent));
  1490               deleteNode(~(classes[r].parent));
  1491               if (less(~(classes[r].parent), ~(classes[l].parent))) {
  1492                 nodes[~(classes[l].parent)].prio =
  1493                   nodes[~(classes[r].parent)].prio;
  1494                 nodes[~(classes[l].parent)].item =
  1495                   nodes[~(classes[r].parent)].item;
  1496               }
  1497             } else {
  1498               int new_parent = newNode();
  1499               nodes[new_parent].next = nodes[new_parent].prev = -1;
  1500               push(new_parent, ~(classes[l].parent));
  1501               pushRight(new_parent, ~(classes[r].parent));
  1502               setPrio(new_parent);
  1503 
  1504               classes[l].parent = ~new_parent;
  1505               classes[l].depth += 1;
  1506               nodes[new_parent].parent = ~l;
  1507             }
  1508           }
  1509           if (classes[r].next != -1) {
  1510             classes[classes[r].next].prev = classes[r].prev;
  1511           }
  1512           classes[classes[r].prev].next = classes[r].next;
  1513 
  1514           classes[r].prev = classes[l].right;
  1515           classes[classes[l].right].next = r;
  1516           classes[l].right = r;
  1517           classes[r].parent = l;
  1518 
  1519           classes[r].next = -1;
  1520           classes[r].depth = rln;
  1521         }
  1522       }
  1523       return class_id;
  1524     }
  1525 
  1526     /// \brief Split the class to subclasses.
  1527     ///
  1528     /// The current function splits the given class. The join, which
  1529     /// made the current class, stored a reference to the
  1530     /// subclasses. The \c splitClass() member restores the classes
  1531     /// and creates the heaps. The parameter is an STL output iterator
  1532     /// which will be filled with the subclass ids. The time
  1533     /// complexity is O(log(n)*k) where n is the overall number of
  1534     /// nodes in the splitted classes and k is the number of the
  1535     /// classes.
  1536     template <typename Iterator>
  1537     void split(int cls, Iterator out) {
  1538       std::vector<int> cs;
  1539       { // splitting union-find
  1540         int id = cls;
  1541         int l = classes[id].left;
  1542 
  1543         classes[l].parent = classes[id].parent;
  1544         classes[l].depth = classes[id].depth;
  1545 
  1546         nodes[~(classes[l].parent)].parent = ~l;
  1547 
  1548         *out++ = l;
  1549 
  1550         while (l != -1) {
  1551           cs.push_back(l);
  1552           l = classes[l].next;
  1553         }
  1554 
  1555         classes[classes[id].right].next = first_class;
  1556         classes[first_class].prev = classes[id].right;
  1557         first_class = classes[id].left;
  1558 
  1559         if (classes[id].next != -1) {
  1560           classes[classes[id].next].prev = classes[id].prev;
  1561         }
  1562         classes[classes[id].prev].next = classes[id].next;
  1563 
  1564         deleteClass(id);
  1565       }
  1566 
  1567       {
  1568         for (int i = 1; i < int(cs.size()); ++i) {
  1569           int l = classes[cs[i]].depth;
  1570           while (nodes[nodes[l].parent].left == l) {
  1571             l = nodes[l].parent;
  1572           }
  1573           int r = l;
  1574           while (nodes[l].parent >= 0) {
  1575             l = nodes[l].parent;
  1576             int new_node = newNode();
  1577 
  1578             nodes[new_node].prev = -1;
  1579             nodes[new_node].next = -1;
  1580 
  1581             split(r, new_node);
  1582             pushAfter(l, new_node);
  1583             setPrio(l);
  1584             setPrio(new_node);
  1585             r = new_node;
  1586           }
  1587           classes[cs[i]].parent = ~r;
  1588           classes[cs[i]].depth = classes[~(nodes[l].parent)].depth;
  1589           nodes[r].parent = ~cs[i];
  1590 
  1591           nodes[l].next = -1;
  1592           nodes[r].prev = -1;
  1593 
  1594           repairRight(~(nodes[l].parent));
  1595           repairLeft(cs[i]);
  1596 
  1597           *out++ = cs[i];
  1598         }
  1599       }
  1600     }
  1601 
  1602     /// \brief Gives back the priority of the current item.
  1603     ///
  1604     /// \return Gives back the priority of the current item.
  1605     const Value& operator[](const Item& item) const {
  1606       return nodes[index[item]].prio;
  1607     }
  1608 
  1609     /// \brief Sets the priority of the current item.
  1610     ///
  1611     /// Sets the priority of the current item.
  1612     void set(const Item& item, const Value& prio) {
  1613       if (comp(prio, nodes[index[item]].prio)) {
  1614         decrease(item, prio);
  1615       } else if (!comp(prio, nodes[index[item]].prio)) {
  1616         increase(item, prio);
  1617       }
  1618     }
  1619 
  1620     /// \brief Increase the priority of the current item.
  1621     ///
  1622     /// Increase the priority of the current item.
  1623     void increase(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 && nodes[kd].item == item) {
  1628         setPrio(kd);
  1629         kd = nodes[kd].parent;
  1630       }
  1631     }
  1632 
  1633     /// \brief Increase the priority of the current item.
  1634     ///
  1635     /// Increase the priority of the current item.
  1636     void decrease(const Item& item, const Value& prio) {
  1637       int id = index[item];
  1638       int kd = nodes[id].parent;
  1639       nodes[id].prio = prio;
  1640       while (kd >= 0 && less(id, kd)) {
  1641         nodes[kd].prio = prio;
  1642         nodes[kd].item = item;
  1643         kd = nodes[kd].parent;
  1644       }
  1645     }
  1646 
  1647     /// \brief Gives back the minimum priority of the class.
  1648     ///
  1649     /// \return Gives back the minimum priority of the class.
  1650     const Value& classPrio(int cls) const {
  1651       return nodes[~(classes[cls].parent)].prio;
  1652     }
  1653 
  1654     /// \brief Gives back the minimum priority item of the class.
  1655     ///
  1656     /// \return Gives back the minimum priority item of the class.
  1657     const Item& classTop(int cls) const {
  1658       return nodes[~(classes[cls].parent)].item;
  1659     }
  1660 
  1661     /// \brief Gives back a representant item of the class.
  1662     ///
  1663     /// The representant is indpendent from the priorities of the
  1664     /// items.
  1665     /// \return Gives back a representant item of the class.
  1666     const Item& classRep(int id) const {
  1667       int parent = classes[id].parent;
  1668       return nodes[parent >= 0 ? classes[id].depth : leftNode(id)].item;
  1669     }
  1670 
  1671     /// \brief LEMON style iterator for the items of a class.
  1672     ///
  1673     /// ClassIt is a lemon style iterator for the components. It iterates
  1674     /// on the items of a class. By example if you want to iterate on
  1675     /// each items of each classes then you may write the next code.
  1676     ///\code
  1677     /// for (ClassIt cit(huf); cit != INVALID; ++cit) {
  1678     ///   std::cout << "Class: ";
  1679     ///   for (ItemIt iit(huf, cit); iit != INVALID; ++iit) {
  1680     ///     std::cout << toString(iit) << ' ' << std::endl;
  1681     ///   }
  1682     ///   std::cout << std::endl;
  1683     /// }
  1684     ///\endcode
  1685     class ItemIt {
  1686     private:
  1687 
  1688       const HeapUnionFind* _huf;
  1689       int _id, _lid;
  1690 
  1691     public:
  1692 
  1693       /// \brief Default constructor
  1694       ///
  1695       /// Default constructor
  1696       ItemIt() {}
  1697 
  1698       ItemIt(const HeapUnionFind& huf, int cls) : _huf(&huf) {
  1699         int id = cls;
  1700         int parent = _huf->classes[id].parent;
  1701         if (parent >= 0) {
  1702           _id = _huf->classes[id].depth;
  1703           if (_huf->classes[id].next != -1) {
  1704             _lid = _huf->classes[_huf->classes[id].next].depth;
  1705           } else {
  1706             _lid = -1;
  1707           }
  1708         } else {
  1709           _id = _huf->leftNode(id);
  1710           _lid = -1;
  1711         }
  1712       }
  1713 
  1714       /// \brief Increment operator
  1715       ///
  1716       /// It steps to the next item in the class.
  1717       ItemIt& operator++() {
  1718         _id = _huf->nextNode(_id);
  1719         return *this;
  1720       }
  1721 
  1722       /// \brief Conversion operator
  1723       ///
  1724       /// It converts the iterator to the current item.
  1725       operator const Item&() const {
  1726         return _huf->nodes[_id].item;
  1727       }
  1728 
  1729       /// \brief Equality operator
  1730       ///
  1731       /// Equality operator
  1732       bool operator==(const ItemIt& i) {
  1733         return i._id == _id;
  1734       }
  1735 
  1736       /// \brief Inequality operator
  1737       ///
  1738       /// Inequality operator
  1739       bool operator!=(const ItemIt& i) {
  1740         return i._id != _id;
  1741       }
  1742 
  1743       /// \brief Equality operator
  1744       ///
  1745       /// Equality operator
  1746       bool operator==(Invalid) {
  1747         return _id == _lid;
  1748       }
  1749 
  1750       /// \brief Inequality operator
  1751       ///
  1752       /// Inequality operator
  1753       bool operator!=(Invalid) {
  1754         return _id != _lid;
  1755       }
  1756 
  1757     };
  1758 
  1759     /// \brief Class iterator
  1760     ///
  1761     /// The iterator stores
  1762     class ClassIt {
  1763     private:
  1764 
  1765       const HeapUnionFind* _huf;
  1766       int _id;
  1767 
  1768     public:
  1769 
  1770       ClassIt(const HeapUnionFind& huf)
  1771         : _huf(&huf), _id(huf.first_class) {}
  1772 
  1773       ClassIt(const HeapUnionFind& huf, int cls)
  1774         : _huf(&huf), _id(huf.classes[cls].left) {}
  1775 
  1776       ClassIt(Invalid) : _huf(0), _id(-1) {}
  1777 
  1778       const ClassIt& operator++() {
  1779         _id = _huf->classes[_id].next;
  1780         return *this;
  1781       }
  1782 
  1783       /// \brief Equality operator
  1784       ///
  1785       /// Equality operator
  1786       bool operator==(const ClassIt& i) {
  1787         return i._id == _id;
  1788       }
  1789 
  1790       /// \brief Inequality operator
  1791       ///
  1792       /// Inequality operator
  1793       bool operator!=(const ClassIt& i) {
  1794         return i._id != _id;
  1795       }
  1796 
  1797       operator int() const {
  1798         return _id;
  1799       }
  1800 
  1801     };
  1802 
  1803   };
  1804 
  1805   //! @}
  1806 
  1807 } //namespace lemon
  1808 
  1809 #endif //LEMON_UNION_FIND_H