lemon/unionfind.h
author Alpar Juttner <alpar@cs.elte.hu>
Fri, 12 Dec 2008 21:52:53 +0000
branch1.0
changeset 350 d5cfcff85dfd
parent 220 a5d8c039f218
child 351 561be42c4b99
permissions -rw-r--r--
Merge bugfix #195
     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-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/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(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