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