lemon/unionfind.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 03 Apr 2009 18:59:15 +0200
changeset 600 6ac5d9ae1d3d
parent 438 81d40f1c850c
child 550 c5fd2d996909
permissions -rw-r--r--
Support real types + numerical stability fix in NS (#254)

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