lemon/fib_heap.h
author Peter Kovacs <kpeter@inf.elte.hu>
Thu, 12 Nov 2009 23:26:13 +0100
changeset 806 fa6f37d7a25b
parent 710 f1fe0ddad6f7
permissions -rw-r--r--
Entirely rework CapacityScaling (#180)

- Use the new interface similarly to NetworkSimplex.
- Rework the implementation using an efficient internal structure
for handling the residual network. This improvement made the
code much faster (up to 2-5 times faster on large graphs).
- Handle GEQ supply type (LEQ is not supported).
- Handle negative costs for arcs of finite capacity.
(Note that this algorithm cannot handle arcs of negative cost
and infinite upper bound, thus it returns UNBOUNDED if such
an arc exists.)
- Extend the documentation.
     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_FIB_HEAP_H
    20 #define LEMON_FIB_HEAP_H
    21 
    22 ///\file
    23 ///\ingroup heaps
    24 ///\brief Fibonacci heap implementation.
    25 
    26 #include <vector>
    27 #include <utility>
    28 #include <functional>
    29 #include <lemon/math.h>
    30 
    31 namespace lemon {
    32 
    33   /// \ingroup heaps
    34   ///
    35   /// \brief Fibonacci heap data structure.
    36   ///
    37   /// This class implements the \e Fibonacci \e heap data structure.
    38   /// It fully conforms to the \ref concepts::Heap "heap concept".
    39   ///
    40   /// The methods \ref increase() and \ref erase() are not efficient in a
    41   /// Fibonacci heap. In case of many calls of these operations, it is
    42   /// better to use other heap structure, e.g. \ref BinHeap "binary heap".
    43   ///
    44   /// \tparam PR Type of the priorities of the items.
    45   /// \tparam IM A read-writable item map with \c int values, used
    46   /// internally to handle the cross references.
    47   /// \tparam CMP A functor class for comparing the priorities.
    48   /// The default is \c std::less<PR>.
    49 #ifdef DOXYGEN
    50   template <typename PR, typename IM, typename CMP>
    51 #else
    52   template <typename PR, typename IM, typename CMP = std::less<PR> >
    53 #endif
    54   class FibHeap {
    55   public:
    56 
    57     /// Type of the item-int map.
    58     typedef IM ItemIntMap;
    59     /// Type of the priorities.
    60     typedef PR Prio;
    61     /// Type of the items stored in the heap.
    62     typedef typename ItemIntMap::Key Item;
    63     /// Type of the item-priority pairs.
    64     typedef std::pair<Item,Prio> Pair;
    65     /// Functor type for comparing the priorities.
    66     typedef CMP Compare;
    67 
    68   private:
    69     class Store;
    70 
    71     std::vector<Store> _data;
    72     int _minimum;
    73     ItemIntMap &_iim;
    74     Compare _comp;
    75     int _num;
    76 
    77   public:
    78 
    79     /// \brief Type to represent the states of the items.
    80     ///
    81     /// Each item has a state associated to it. It can be "in heap",
    82     /// "pre-heap" or "post-heap". The latter two are indifferent from the
    83     /// heap's point of view, but may be useful to the user.
    84     ///
    85     /// The item-int map must be initialized in such way that it assigns
    86     /// \c PRE_HEAP (<tt>-1</tt>) to any element to be put in the heap.
    87     enum State {
    88       IN_HEAP = 0,    ///< = 0.
    89       PRE_HEAP = -1,  ///< = -1.
    90       POST_HEAP = -2  ///< = -2.
    91     };
    92 
    93     /// \brief Constructor.
    94     ///
    95     /// Constructor.
    96     /// \param map A map that assigns \c int values to the items.
    97     /// It is used internally to handle the cross references.
    98     /// The assigned value must be \c PRE_HEAP (<tt>-1</tt>) for each item.
    99     explicit FibHeap(ItemIntMap &map)
   100       : _minimum(0), _iim(map), _num() {}
   101 
   102     /// \brief Constructor.
   103     ///
   104     /// Constructor.
   105     /// \param map A map that assigns \c int values to the items.
   106     /// It is used internally to handle the cross references.
   107     /// The assigned value must be \c PRE_HEAP (<tt>-1</tt>) for each item.
   108     /// \param comp The function object used for comparing the priorities.
   109     FibHeap(ItemIntMap &map, const Compare &comp)
   110       : _minimum(0), _iim(map), _comp(comp), _num() {}
   111 
   112     /// \brief The number of items stored in the heap.
   113     ///
   114     /// This function returns the number of items stored in the heap.
   115     int size() const { return _num; }
   116 
   117     /// \brief Check if the heap is empty.
   118     ///
   119     /// This function returns \c true if the heap is empty.
   120     bool empty() const { return _num==0; }
   121 
   122     /// \brief Make the heap empty.
   123     ///
   124     /// This functon makes the heap empty.
   125     /// It does not change the cross reference map. If you want to reuse
   126     /// a heap that is not surely empty, you should first clear it and
   127     /// then you should set the cross reference map to \c PRE_HEAP
   128     /// for each item.
   129     void clear() {
   130       _data.clear(); _minimum = 0; _num = 0;
   131     }
   132 
   133     /// \brief Insert an item into the heap with the given priority.
   134     ///
   135     /// This function inserts the given item into the heap with the
   136     /// given priority.
   137     /// \param item The item to insert.
   138     /// \param prio The priority of the item.
   139     /// \pre \e item must not be stored in the heap.
   140     void push (const Item& item, const Prio& prio) {
   141       int i=_iim[item];
   142       if ( i < 0 ) {
   143         int s=_data.size();
   144         _iim.set( item, s );
   145         Store st;
   146         st.name=item;
   147         _data.push_back(st);
   148         i=s;
   149       } else {
   150         _data[i].parent=_data[i].child=-1;
   151         _data[i].degree=0;
   152         _data[i].in=true;
   153         _data[i].marked=false;
   154       }
   155 
   156       if ( _num ) {
   157         _data[_data[_minimum].right_neighbor].left_neighbor=i;
   158         _data[i].right_neighbor=_data[_minimum].right_neighbor;
   159         _data[_minimum].right_neighbor=i;
   160         _data[i].left_neighbor=_minimum;
   161         if ( _comp( prio, _data[_minimum].prio) ) _minimum=i;
   162       } else {
   163         _data[i].right_neighbor=_data[i].left_neighbor=i;
   164         _minimum=i;
   165       }
   166       _data[i].prio=prio;
   167       ++_num;
   168     }
   169 
   170     /// \brief Return the item having minimum priority.
   171     ///
   172     /// This function returns the item having minimum priority.
   173     /// \pre The heap must be non-empty.
   174     Item top() const { return _data[_minimum].name; }
   175 
   176     /// \brief The minimum priority.
   177     ///
   178     /// This function returns the minimum priority.
   179     /// \pre The heap must be non-empty.
   180     Prio prio() const { return _data[_minimum].prio; }
   181 
   182     /// \brief Remove the item having minimum priority.
   183     ///
   184     /// This function removes the item having minimum priority.
   185     /// \pre The heap must be non-empty.
   186     void pop() {
   187       /*The first case is that there are only one root.*/
   188       if ( _data[_minimum].left_neighbor==_minimum ) {
   189         _data[_minimum].in=false;
   190         if ( _data[_minimum].degree!=0 ) {
   191           makeRoot(_data[_minimum].child);
   192           _minimum=_data[_minimum].child;
   193           balance();
   194         }
   195       } else {
   196         int right=_data[_minimum].right_neighbor;
   197         unlace(_minimum);
   198         _data[_minimum].in=false;
   199         if ( _data[_minimum].degree > 0 ) {
   200           int left=_data[_minimum].left_neighbor;
   201           int child=_data[_minimum].child;
   202           int last_child=_data[child].left_neighbor;
   203 
   204           makeRoot(child);
   205 
   206           _data[left].right_neighbor=child;
   207           _data[child].left_neighbor=left;
   208           _data[right].left_neighbor=last_child;
   209           _data[last_child].right_neighbor=right;
   210         }
   211         _minimum=right;
   212         balance();
   213       } // the case where there are more roots
   214       --_num;
   215     }
   216 
   217     /// \brief Remove the given item from the heap.
   218     ///
   219     /// This function removes the given item from the heap if it is
   220     /// already stored.
   221     /// \param item The item to delete.
   222     /// \pre \e item must be in the heap.
   223     void erase (const Item& item) {
   224       int i=_iim[item];
   225 
   226       if ( i >= 0 && _data[i].in ) {
   227         if ( _data[i].parent!=-1 ) {
   228           int p=_data[i].parent;
   229           cut(i,p);
   230           cascade(p);
   231         }
   232         _minimum=i;     //As if its prio would be -infinity
   233         pop();
   234       }
   235     }
   236 
   237     /// \brief The priority of the given item.
   238     ///
   239     /// This function returns the priority of the given item.
   240     /// \param item The item.
   241     /// \pre \e item must be in the heap.
   242     Prio operator[](const Item& item) const {
   243       return _data[_iim[item]].prio;
   244     }
   245 
   246     /// \brief Set the priority of an item or insert it, if it is
   247     /// not stored in the heap.
   248     ///
   249     /// This method sets the priority of the given item if it is
   250     /// already stored in the heap. Otherwise it inserts the given
   251     /// item into the heap with the given priority.
   252     /// \param item The item.
   253     /// \param prio The priority.
   254     void set (const Item& item, const Prio& prio) {
   255       int i=_iim[item];
   256       if ( i >= 0 && _data[i].in ) {
   257         if ( _comp(prio, _data[i].prio) ) decrease(item, prio);
   258         if ( _comp(_data[i].prio, prio) ) increase(item, prio);
   259       } else push(item, prio);
   260     }
   261 
   262     /// \brief Decrease the priority of an item to the given value.
   263     ///
   264     /// This function decreases the priority of an item to the given value.
   265     /// \param item The item.
   266     /// \param prio The priority.
   267     /// \pre \e item must be stored in the heap with priority at least \e prio.
   268     void decrease (const Item& item, const Prio& prio) {
   269       int i=_iim[item];
   270       _data[i].prio=prio;
   271       int p=_data[i].parent;
   272 
   273       if ( p!=-1 && _comp(prio, _data[p].prio) ) {
   274         cut(i,p);
   275         cascade(p);
   276       }
   277       if ( _comp(prio, _data[_minimum].prio) ) _minimum=i;
   278     }
   279 
   280     /// \brief Increase the priority of an item to the given value.
   281     ///
   282     /// This function increases the priority of an item to the given value.
   283     /// \param item The item.
   284     /// \param prio The priority.
   285     /// \pre \e item must be stored in the heap with priority at most \e prio.
   286     void increase (const Item& item, const Prio& prio) {
   287       erase(item);
   288       push(item, prio);
   289     }
   290 
   291     /// \brief Return the state of an item.
   292     ///
   293     /// This method returns \c PRE_HEAP if the given item has never
   294     /// been in the heap, \c IN_HEAP if it is in the heap at the moment,
   295     /// and \c POST_HEAP otherwise.
   296     /// In the latter case it is possible that the item will get back
   297     /// to the heap again.
   298     /// \param item The item.
   299     State state(const Item &item) const {
   300       int i=_iim[item];
   301       if( i>=0 ) {
   302         if ( _data[i].in ) i=0;
   303         else i=-2;
   304       }
   305       return State(i);
   306     }
   307 
   308     /// \brief Set the state of an item in the heap.
   309     ///
   310     /// This function sets the state of the given item in the heap.
   311     /// It can be used to manually clear the heap when it is important
   312     /// to achive better time complexity.
   313     /// \param i The item.
   314     /// \param st The state. It should not be \c IN_HEAP.
   315     void state(const Item& i, State st) {
   316       switch (st) {
   317       case POST_HEAP:
   318       case PRE_HEAP:
   319         if (state(i) == IN_HEAP) {
   320           erase(i);
   321         }
   322         _iim[i] = st;
   323         break;
   324       case IN_HEAP:
   325         break;
   326       }
   327     }
   328 
   329   private:
   330 
   331     void balance() {
   332 
   333       int maxdeg=int( std::floor( 2.08*log(double(_data.size()))))+1;
   334 
   335       std::vector<int> A(maxdeg,-1);
   336 
   337       /*
   338        *Recall that now minimum does not point to the minimum prio element.
   339        *We set minimum to this during balance().
   340        */
   341       int anchor=_data[_minimum].left_neighbor;
   342       int next=_minimum;
   343       bool end=false;
   344 
   345       do {
   346         int active=next;
   347         if ( anchor==active ) end=true;
   348         int d=_data[active].degree;
   349         next=_data[active].right_neighbor;
   350 
   351         while (A[d]!=-1) {
   352           if( _comp(_data[active].prio, _data[A[d]].prio) ) {
   353             fuse(active,A[d]);
   354           } else {
   355             fuse(A[d],active);
   356             active=A[d];
   357           }
   358           A[d]=-1;
   359           ++d;
   360         }
   361         A[d]=active;
   362       } while ( !end );
   363 
   364 
   365       while ( _data[_minimum].parent >=0 )
   366         _minimum=_data[_minimum].parent;
   367       int s=_minimum;
   368       int m=_minimum;
   369       do {
   370         if ( _comp(_data[s].prio, _data[_minimum].prio) ) _minimum=s;
   371         s=_data[s].right_neighbor;
   372       } while ( s != m );
   373     }
   374 
   375     void makeRoot(int c) {
   376       int s=c;
   377       do {
   378         _data[s].parent=-1;
   379         s=_data[s].right_neighbor;
   380       } while ( s != c );
   381     }
   382 
   383     void cut(int a, int b) {
   384       /*
   385        *Replacing a from the children of b.
   386        */
   387       --_data[b].degree;
   388 
   389       if ( _data[b].degree !=0 ) {
   390         int child=_data[b].child;
   391         if ( child==a )
   392           _data[b].child=_data[child].right_neighbor;
   393         unlace(a);
   394       }
   395 
   396 
   397       /*Lacing a to the roots.*/
   398       int right=_data[_minimum].right_neighbor;
   399       _data[_minimum].right_neighbor=a;
   400       _data[a].left_neighbor=_minimum;
   401       _data[a].right_neighbor=right;
   402       _data[right].left_neighbor=a;
   403 
   404       _data[a].parent=-1;
   405       _data[a].marked=false;
   406     }
   407 
   408     void cascade(int a) {
   409       if ( _data[a].parent!=-1 ) {
   410         int p=_data[a].parent;
   411 
   412         if ( _data[a].marked==false ) _data[a].marked=true;
   413         else {
   414           cut(a,p);
   415           cascade(p);
   416         }
   417       }
   418     }
   419 
   420     void fuse(int a, int b) {
   421       unlace(b);
   422 
   423       /*Lacing b under a.*/
   424       _data[b].parent=a;
   425 
   426       if (_data[a].degree==0) {
   427         _data[b].left_neighbor=b;
   428         _data[b].right_neighbor=b;
   429         _data[a].child=b;
   430       } else {
   431         int child=_data[a].child;
   432         int last_child=_data[child].left_neighbor;
   433         _data[child].left_neighbor=b;
   434         _data[b].right_neighbor=child;
   435         _data[last_child].right_neighbor=b;
   436         _data[b].left_neighbor=last_child;
   437       }
   438 
   439       ++_data[a].degree;
   440 
   441       _data[b].marked=false;
   442     }
   443 
   444     /*
   445      *It is invoked only if a has siblings.
   446      */
   447     void unlace(int a) {
   448       int leftn=_data[a].left_neighbor;
   449       int rightn=_data[a].right_neighbor;
   450       _data[leftn].right_neighbor=rightn;
   451       _data[rightn].left_neighbor=leftn;
   452     }
   453 
   454 
   455     class Store {
   456       friend class FibHeap;
   457 
   458       Item name;
   459       int parent;
   460       int left_neighbor;
   461       int right_neighbor;
   462       int child;
   463       int degree;
   464       bool marked;
   465       bool in;
   466       Prio prio;
   467 
   468       Store() : parent(-1), child(-1), degree(), marked(false), in(true) {}
   469     };
   470   };
   471 
   472 } //namespace lemon
   473 
   474 #endif //LEMON_FIB_HEAP_H
   475