lemon/hao_orlin.h
author Balazs Dezso <deba@inf.elte.hu>
Mon, 01 Dec 2008 23:15:15 +0100
changeset 410 eac19fb31a09
child 411 01c443515ad2
permissions -rw-r--r--
Simple test for HaoOrlin algorithm class (#58)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2008
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_HAO_ORLIN_H
    20 #define LEMON_HAO_ORLIN_H
    21 
    22 #include <vector>
    23 #include <list>
    24 #include <limits>
    25 
    26 #include <lemon/maps.h>
    27 #include <lemon/core.h>
    28 #include <lemon/tolerance.h>
    29 
    30 /// \file
    31 /// \ingroup min_cut
    32 /// \brief Implementation of the Hao-Orlin algorithm.
    33 ///
    34 /// Implementation of the Hao-Orlin algorithm class for testing network
    35 /// reliability.
    36 
    37 namespace lemon {
    38 
    39   /// \ingroup min_cut
    40   ///
    41   /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
    42   ///
    43   /// Hao-Orlin calculates a minimum cut in a directed graph
    44   /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and
    45   /// consists of two phases: in the first phase it determines a
    46   /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
    47   /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal
    48   /// out-degree) and in the second phase it determines a minimum cut
    49   /// with \f$ source \f$ on the sink-side (i.e. a set
    50   /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal
    51   /// out-degree). Obviously, the smaller of these two cuts will be a
    52   /// minimum cut of \f$ D \f$. The algorithm is a modified
    53   /// push-relabel preflow algorithm and our implementation calculates
    54   /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
    55   /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
    56   /// purpose of such algorithm is testing network reliability. For an
    57   /// undirected graph you can run just the first phase of the
    58   /// algorithm or you can use the algorithm of Nagamochi and Ibaraki
    59   /// which solves the undirected problem in
    60   /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the
    61   /// NagamochiIbaraki algorithm class.
    62   ///
    63   /// \param _Digraph is the graph type of the algorithm.
    64   /// \param _CapacityMap is an edge map of capacities which should
    65   /// be any numreric type. The default type is _Digraph::ArcMap<int>.
    66   /// \param _Tolerance is the handler of the inexact computation. The
    67   /// default type for this is Tolerance<CapacityMap::Value>.
    68 #ifdef DOXYGEN
    69   template <typename _Digraph, typename _CapacityMap, typename _Tolerance>
    70 #else
    71   template <typename _Digraph,
    72             typename _CapacityMap = typename _Digraph::template ArcMap<int>,
    73             typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
    74 #endif
    75   class HaoOrlin {
    76   private:
    77 
    78     typedef _Digraph Digraph;
    79     typedef _CapacityMap CapacityMap;
    80     typedef _Tolerance Tolerance;
    81 
    82     typedef typename CapacityMap::Value Value;
    83 
    84     TEMPLATE_GRAPH_TYPEDEFS(Digraph);
    85 
    86     const Digraph& _graph;
    87     const CapacityMap* _capacity;
    88 
    89     typedef typename Digraph::template ArcMap<Value> FlowMap;
    90     FlowMap* _flow;
    91 
    92     Node _source;
    93 
    94     int _node_num;
    95 
    96     // Bucketing structure
    97     std::vector<Node> _first, _last;
    98     typename Digraph::template NodeMap<Node>* _next;
    99     typename Digraph::template NodeMap<Node>* _prev;
   100     typename Digraph::template NodeMap<bool>* _active;
   101     typename Digraph::template NodeMap<int>* _bucket;
   102 
   103     std::vector<bool> _dormant;
   104 
   105     std::list<std::list<int> > _sets;
   106     std::list<int>::iterator _highest;
   107 
   108     typedef typename Digraph::template NodeMap<Value> ExcessMap;
   109     ExcessMap* _excess;
   110 
   111     typedef typename Digraph::template NodeMap<bool> SourceSetMap;
   112     SourceSetMap* _source_set;
   113 
   114     Value _min_cut;
   115 
   116     typedef typename Digraph::template NodeMap<bool> MinCutMap;
   117     MinCutMap* _min_cut_map;
   118 
   119     Tolerance _tolerance;
   120 
   121   public:
   122 
   123     /// \brief Constructor
   124     ///
   125     /// Constructor of the algorithm class.
   126     HaoOrlin(const Digraph& graph, const CapacityMap& capacity,
   127              const Tolerance& tolerance = Tolerance()) :
   128       _graph(graph), _capacity(&capacity), _flow(0), _source(),
   129       _node_num(), _first(), _last(), _next(0), _prev(0),
   130       _active(0), _bucket(0), _dormant(), _sets(), _highest(),
   131       _excess(0), _source_set(0), _min_cut(), _min_cut_map(0),
   132       _tolerance(tolerance) {}
   133 
   134     ~HaoOrlin() {
   135       if (_min_cut_map) {
   136         delete _min_cut_map;
   137       }
   138       if (_source_set) {
   139         delete _source_set;
   140       }
   141       if (_excess) {
   142         delete _excess;
   143       }
   144       if (_next) {
   145         delete _next;
   146       }
   147       if (_prev) {
   148         delete _prev;
   149       }
   150       if (_active) {
   151         delete _active;
   152       }
   153       if (_bucket) {
   154         delete _bucket;
   155       }
   156       if (_flow) {
   157         delete _flow;
   158       }
   159     }
   160 
   161   private:
   162 
   163     void activate(const Node& i) {
   164       _active->set(i, true);
   165 
   166       int bucket = (*_bucket)[i];
   167 
   168       if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;
   169       //unlace
   170       _next->set((*_prev)[i], (*_next)[i]);
   171       if ((*_next)[i] != INVALID) {
   172         _prev->set((*_next)[i], (*_prev)[i]);
   173       } else {
   174         _last[bucket] = (*_prev)[i];
   175       }
   176       //lace
   177       _next->set(i, _first[bucket]);
   178       _prev->set(_first[bucket], i);
   179       _prev->set(i, INVALID);
   180       _first[bucket] = i;
   181     }
   182 
   183     void deactivate(const Node& i) {
   184       _active->set(i, false);
   185       int bucket = (*_bucket)[i];
   186 
   187       if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
   188 
   189       //unlace
   190       _prev->set((*_next)[i], (*_prev)[i]);
   191       if ((*_prev)[i] != INVALID) {
   192         _next->set((*_prev)[i], (*_next)[i]);
   193       } else {
   194         _first[bucket] = (*_next)[i];
   195       }
   196       //lace
   197       _prev->set(i, _last[bucket]);
   198       _next->set(_last[bucket], i);
   199       _next->set(i, INVALID);
   200       _last[bucket] = i;
   201     }
   202 
   203     void addItem(const Node& i, int bucket) {
   204       (*_bucket)[i] = bucket;
   205       if (_last[bucket] != INVALID) {
   206         _prev->set(i, _last[bucket]);
   207         _next->set(_last[bucket], i);
   208         _next->set(i, INVALID);
   209         _last[bucket] = i;
   210       } else {
   211         _prev->set(i, INVALID);
   212         _first[bucket] = i;
   213         _next->set(i, INVALID);
   214         _last[bucket] = i;
   215       }
   216     }
   217 
   218     void findMinCutOut() {
   219 
   220       for (NodeIt n(_graph); n != INVALID; ++n) {
   221         _excess->set(n, 0);
   222       }
   223 
   224       for (ArcIt a(_graph); a != INVALID; ++a) {
   225         _flow->set(a, 0);
   226       }
   227 
   228       int bucket_num = 1;
   229 
   230       {
   231         typename Digraph::template NodeMap<bool> reached(_graph, false);
   232 
   233         reached.set(_source, true);
   234 
   235         bool first_set = true;
   236 
   237         for (NodeIt t(_graph); t != INVALID; ++t) {
   238           if (reached[t]) continue;
   239           _sets.push_front(std::list<int>());
   240           _sets.front().push_front(bucket_num);
   241           _dormant[bucket_num] = !first_set;
   242 
   243           _bucket->set(t, bucket_num);
   244           _first[bucket_num] = _last[bucket_num] = t;
   245           _next->set(t, INVALID);
   246           _prev->set(t, INVALID);
   247 
   248           ++bucket_num;
   249 
   250           std::vector<Node> queue;
   251           queue.push_back(t);
   252           reached.set(t, true);
   253 
   254           while (!queue.empty()) {
   255             _sets.front().push_front(bucket_num);
   256             _dormant[bucket_num] = !first_set;
   257             _first[bucket_num] = _last[bucket_num] = INVALID;
   258 
   259             std::vector<Node> nqueue;
   260             for (int i = 0; i < int(queue.size()); ++i) {
   261               Node n = queue[i];
   262               for (InArcIt a(_graph, n); a != INVALID; ++a) {
   263                 Node u = _graph.source(a);
   264                 if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
   265                   reached.set(u, true);
   266                   addItem(u, bucket_num);
   267                   nqueue.push_back(u);
   268                 }
   269               }
   270             }
   271             queue.swap(nqueue);
   272             ++bucket_num;
   273           }
   274           _sets.front().pop_front();
   275           --bucket_num;
   276           first_set = false;
   277         }
   278 
   279         _bucket->set(_source, 0);
   280         _dormant[0] = true;
   281       }
   282       _source_set->set(_source, true);
   283 
   284       Node target = _last[_sets.back().back()];
   285       {
   286         for (OutArcIt a(_graph, _source); a != INVALID; ++a) {
   287           if (_tolerance.positive((*_capacity)[a])) {
   288             Node u = _graph.target(a);
   289             _flow->set(a, (*_capacity)[a]);
   290             _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
   291             if (!(*_active)[u] && u != _source) {
   292               activate(u);
   293             }
   294           }
   295         }
   296 
   297         if ((*_active)[target]) {
   298           deactivate(target);
   299         }
   300 
   301         _highest = _sets.back().begin();
   302         while (_highest != _sets.back().end() &&
   303                !(*_active)[_first[*_highest]]) {
   304           ++_highest;
   305         }
   306       }
   307 
   308       while (true) {
   309         while (_highest != _sets.back().end()) {
   310           Node n = _first[*_highest];
   311           Value excess = (*_excess)[n];
   312           int next_bucket = _node_num;
   313 
   314           int under_bucket;
   315           if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
   316             under_bucket = -1;
   317           } else {
   318             under_bucket = *(++std::list<int>::iterator(_highest));
   319           }
   320 
   321           for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   322             Node v = _graph.target(a);
   323             if (_dormant[(*_bucket)[v]]) continue;
   324             Value rem = (*_capacity)[a] - (*_flow)[a];
   325             if (!_tolerance.positive(rem)) continue;
   326             if ((*_bucket)[v] == under_bucket) {
   327               if (!(*_active)[v] && v != target) {
   328                 activate(v);
   329               }
   330               if (!_tolerance.less(rem, excess)) {
   331                 _flow->set(a, (*_flow)[a] + excess);
   332                 _excess->set(v, (*_excess)[v] + excess);
   333                 excess = 0;
   334                 goto no_more_push;
   335               } else {
   336                 excess -= rem;
   337                 _excess->set(v, (*_excess)[v] + rem);
   338                 _flow->set(a, (*_capacity)[a]);
   339               }
   340             } else if (next_bucket > (*_bucket)[v]) {
   341               next_bucket = (*_bucket)[v];
   342             }
   343           }
   344 
   345           for (InArcIt a(_graph, n); a != INVALID; ++a) {
   346             Node v = _graph.source(a);
   347             if (_dormant[(*_bucket)[v]]) continue;
   348             Value rem = (*_flow)[a];
   349             if (!_tolerance.positive(rem)) continue;
   350             if ((*_bucket)[v] == under_bucket) {
   351               if (!(*_active)[v] && v != target) {
   352                 activate(v);
   353               }
   354               if (!_tolerance.less(rem, excess)) {
   355                 _flow->set(a, (*_flow)[a] - excess);
   356                 _excess->set(v, (*_excess)[v] + excess);
   357                 excess = 0;
   358                 goto no_more_push;
   359               } else {
   360                 excess -= rem;
   361                 _excess->set(v, (*_excess)[v] + rem);
   362                 _flow->set(a, 0);
   363               }
   364             } else if (next_bucket > (*_bucket)[v]) {
   365               next_bucket = (*_bucket)[v];
   366             }
   367           }
   368 
   369         no_more_push:
   370 
   371           _excess->set(n, excess);
   372 
   373           if (excess != 0) {
   374             if ((*_next)[n] == INVALID) {
   375               typename std::list<std::list<int> >::iterator new_set =
   376                 _sets.insert(--_sets.end(), std::list<int>());
   377               new_set->splice(new_set->end(), _sets.back(),
   378                               _sets.back().begin(), ++_highest);
   379               for (std::list<int>::iterator it = new_set->begin();
   380                    it != new_set->end(); ++it) {
   381                 _dormant[*it] = true;
   382               }
   383               while (_highest != _sets.back().end() &&
   384                      !(*_active)[_first[*_highest]]) {
   385                 ++_highest;
   386               }
   387             } else if (next_bucket == _node_num) {
   388               _first[(*_bucket)[n]] = (*_next)[n];
   389               _prev->set((*_next)[n], INVALID);
   390 
   391               std::list<std::list<int> >::iterator new_set =
   392                 _sets.insert(--_sets.end(), std::list<int>());
   393 
   394               new_set->push_front(bucket_num);
   395               _bucket->set(n, bucket_num);
   396               _first[bucket_num] = _last[bucket_num] = n;
   397               _next->set(n, INVALID);
   398               _prev->set(n, INVALID);
   399               _dormant[bucket_num] = true;
   400               ++bucket_num;
   401 
   402               while (_highest != _sets.back().end() &&
   403                      !(*_active)[_first[*_highest]]) {
   404                 ++_highest;
   405               }
   406             } else {
   407               _first[*_highest] = (*_next)[n];
   408               _prev->set((*_next)[n], INVALID);
   409 
   410               while (next_bucket != *_highest) {
   411                 --_highest;
   412               }
   413 
   414               if (_highest == _sets.back().begin()) {
   415                 _sets.back().push_front(bucket_num);
   416                 _dormant[bucket_num] = false;
   417                 _first[bucket_num] = _last[bucket_num] = INVALID;
   418                 ++bucket_num;
   419               }
   420               --_highest;
   421 
   422               _bucket->set(n, *_highest);
   423               _next->set(n, _first[*_highest]);
   424               if (_first[*_highest] != INVALID) {
   425                 _prev->set(_first[*_highest], n);
   426               } else {
   427                 _last[*_highest] = n;
   428               }
   429               _first[*_highest] = n;
   430             }
   431           } else {
   432 
   433             deactivate(n);
   434             if (!(*_active)[_first[*_highest]]) {
   435               ++_highest;
   436               if (_highest != _sets.back().end() &&
   437                   !(*_active)[_first[*_highest]]) {
   438                 _highest = _sets.back().end();
   439               }
   440             }
   441           }
   442         }
   443 
   444         if ((*_excess)[target] < _min_cut) {
   445           _min_cut = (*_excess)[target];
   446           for (NodeIt i(_graph); i != INVALID; ++i) {
   447             _min_cut_map->set(i, true);
   448           }
   449           for (std::list<int>::iterator it = _sets.back().begin();
   450                it != _sets.back().end(); ++it) {
   451             Node n = _first[*it];
   452             while (n != INVALID) {
   453               _min_cut_map->set(n, false);
   454               n = (*_next)[n];
   455             }
   456           }
   457         }
   458 
   459         {
   460           Node new_target;
   461           if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
   462             if ((*_next)[target] == INVALID) {
   463               _last[(*_bucket)[target]] = (*_prev)[target];
   464               new_target = (*_prev)[target];
   465             } else {
   466               _prev->set((*_next)[target], (*_prev)[target]);
   467               new_target = (*_next)[target];
   468             }
   469             if ((*_prev)[target] == INVALID) {
   470               _first[(*_bucket)[target]] = (*_next)[target];
   471             } else {
   472               _next->set((*_prev)[target], (*_next)[target]);
   473             }
   474           } else {
   475             _sets.back().pop_back();
   476             if (_sets.back().empty()) {
   477               _sets.pop_back();
   478               if (_sets.empty())
   479                 break;
   480               for (std::list<int>::iterator it = _sets.back().begin();
   481                    it != _sets.back().end(); ++it) {
   482                 _dormant[*it] = false;
   483               }
   484             }
   485             new_target = _last[_sets.back().back()];
   486           }
   487 
   488           _bucket->set(target, 0);
   489 
   490           _source_set->set(target, true);
   491           for (OutArcIt a(_graph, target); a != INVALID; ++a) {
   492             Value rem = (*_capacity)[a] - (*_flow)[a];
   493             if (!_tolerance.positive(rem)) continue;
   494             Node v = _graph.target(a);
   495             if (!(*_active)[v] && !(*_source_set)[v]) {
   496               activate(v);
   497             }
   498             _excess->set(v, (*_excess)[v] + rem);
   499             _flow->set(a, (*_capacity)[a]);
   500           }
   501 
   502           for (InArcIt a(_graph, target); a != INVALID; ++a) {
   503             Value rem = (*_flow)[a];
   504             if (!_tolerance.positive(rem)) continue;
   505             Node v = _graph.source(a);
   506             if (!(*_active)[v] && !(*_source_set)[v]) {
   507               activate(v);
   508             }
   509             _excess->set(v, (*_excess)[v] + rem);
   510             _flow->set(a, 0);
   511           }
   512 
   513           target = new_target;
   514           if ((*_active)[target]) {
   515             deactivate(target);
   516           }
   517 
   518           _highest = _sets.back().begin();
   519           while (_highest != _sets.back().end() &&
   520                  !(*_active)[_first[*_highest]]) {
   521             ++_highest;
   522           }
   523         }
   524       }
   525     }
   526 
   527     void findMinCutIn() {
   528 
   529       for (NodeIt n(_graph); n != INVALID; ++n) {
   530         _excess->set(n, 0);
   531       }
   532 
   533       for (ArcIt a(_graph); a != INVALID; ++a) {
   534         _flow->set(a, 0);
   535       }
   536 
   537       int bucket_num = 1;
   538 
   539       {
   540         typename Digraph::template NodeMap<bool> reached(_graph, false);
   541 
   542         reached.set(_source, true);
   543 
   544         bool first_set = true;
   545 
   546         for (NodeIt t(_graph); t != INVALID; ++t) {
   547           if (reached[t]) continue;
   548           _sets.push_front(std::list<int>());
   549           _sets.front().push_front(bucket_num);
   550           _dormant[bucket_num] = !first_set;
   551 
   552           _bucket->set(t, bucket_num);
   553           _first[bucket_num] = _last[bucket_num] = t;
   554           _next->set(t, INVALID);
   555           _prev->set(t, INVALID);
   556 
   557           ++bucket_num;
   558 
   559           std::vector<Node> queue;
   560           queue.push_back(t);
   561           reached.set(t, true);
   562 
   563           while (!queue.empty()) {
   564             _sets.front().push_front(bucket_num);
   565             _dormant[bucket_num] = !first_set;
   566             _first[bucket_num] = _last[bucket_num] = INVALID;
   567 
   568             std::vector<Node> nqueue;
   569             for (int i = 0; i < int(queue.size()); ++i) {
   570               Node n = queue[i];
   571               for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   572                 Node u = _graph.target(a);
   573                 if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
   574                   reached.set(u, true);
   575                   addItem(u, bucket_num);
   576                   nqueue.push_back(u);
   577                 }
   578               }
   579             }
   580             queue.swap(nqueue);
   581             ++bucket_num;
   582           }
   583           _sets.front().pop_front();
   584           --bucket_num;
   585           first_set = false;
   586         }
   587 
   588         _bucket->set(_source, 0);
   589         _dormant[0] = true;
   590       }
   591       _source_set->set(_source, true);
   592 
   593       Node target = _last[_sets.back().back()];
   594       {
   595         for (InArcIt a(_graph, _source); a != INVALID; ++a) {
   596           if (_tolerance.positive((*_capacity)[a])) {
   597             Node u = _graph.source(a);
   598             _flow->set(a, (*_capacity)[a]);
   599             _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
   600             if (!(*_active)[u] && u != _source) {
   601               activate(u);
   602             }
   603           }
   604         }
   605         if ((*_active)[target]) {
   606           deactivate(target);
   607         }
   608 
   609         _highest = _sets.back().begin();
   610         while (_highest != _sets.back().end() &&
   611                !(*_active)[_first[*_highest]]) {
   612           ++_highest;
   613         }
   614       }
   615 
   616 
   617       while (true) {
   618         while (_highest != _sets.back().end()) {
   619           Node n = _first[*_highest];
   620           Value excess = (*_excess)[n];
   621           int next_bucket = _node_num;
   622 
   623           int under_bucket;
   624           if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
   625             under_bucket = -1;
   626           } else {
   627             under_bucket = *(++std::list<int>::iterator(_highest));
   628           }
   629 
   630           for (InArcIt a(_graph, n); a != INVALID; ++a) {
   631             Node v = _graph.source(a);
   632             if (_dormant[(*_bucket)[v]]) continue;
   633             Value rem = (*_capacity)[a] - (*_flow)[a];
   634             if (!_tolerance.positive(rem)) continue;
   635             if ((*_bucket)[v] == under_bucket) {
   636               if (!(*_active)[v] && v != target) {
   637                 activate(v);
   638               }
   639               if (!_tolerance.less(rem, excess)) {
   640                 _flow->set(a, (*_flow)[a] + excess);
   641                 _excess->set(v, (*_excess)[v] + excess);
   642                 excess = 0;
   643                 goto no_more_push;
   644               } else {
   645                 excess -= rem;
   646                 _excess->set(v, (*_excess)[v] + rem);
   647                 _flow->set(a, (*_capacity)[a]);
   648               }
   649             } else if (next_bucket > (*_bucket)[v]) {
   650               next_bucket = (*_bucket)[v];
   651             }
   652           }
   653 
   654           for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   655             Node v = _graph.target(a);
   656             if (_dormant[(*_bucket)[v]]) continue;
   657             Value rem = (*_flow)[a];
   658             if (!_tolerance.positive(rem)) continue;
   659             if ((*_bucket)[v] == under_bucket) {
   660               if (!(*_active)[v] && v != target) {
   661                 activate(v);
   662               }
   663               if (!_tolerance.less(rem, excess)) {
   664                 _flow->set(a, (*_flow)[a] - excess);
   665                 _excess->set(v, (*_excess)[v] + excess);
   666                 excess = 0;
   667                 goto no_more_push;
   668               } else {
   669                 excess -= rem;
   670                 _excess->set(v, (*_excess)[v] + rem);
   671                 _flow->set(a, 0);
   672               }
   673             } else if (next_bucket > (*_bucket)[v]) {
   674               next_bucket = (*_bucket)[v];
   675             }
   676           }
   677 
   678         no_more_push:
   679 
   680           _excess->set(n, excess);
   681 
   682           if (excess != 0) {
   683             if ((*_next)[n] == INVALID) {
   684               typename std::list<std::list<int> >::iterator new_set =
   685                 _sets.insert(--_sets.end(), std::list<int>());
   686               new_set->splice(new_set->end(), _sets.back(),
   687                               _sets.back().begin(), ++_highest);
   688               for (std::list<int>::iterator it = new_set->begin();
   689                    it != new_set->end(); ++it) {
   690                 _dormant[*it] = true;
   691               }
   692               while (_highest != _sets.back().end() &&
   693                      !(*_active)[_first[*_highest]]) {
   694                 ++_highest;
   695               }
   696             } else if (next_bucket == _node_num) {
   697               _first[(*_bucket)[n]] = (*_next)[n];
   698               _prev->set((*_next)[n], INVALID);
   699 
   700               std::list<std::list<int> >::iterator new_set =
   701                 _sets.insert(--_sets.end(), std::list<int>());
   702 
   703               new_set->push_front(bucket_num);
   704               _bucket->set(n, bucket_num);
   705               _first[bucket_num] = _last[bucket_num] = n;
   706               _next->set(n, INVALID);
   707               _prev->set(n, INVALID);
   708               _dormant[bucket_num] = true;
   709               ++bucket_num;
   710 
   711               while (_highest != _sets.back().end() &&
   712                      !(*_active)[_first[*_highest]]) {
   713                 ++_highest;
   714               }
   715             } else {
   716               _first[*_highest] = (*_next)[n];
   717               _prev->set((*_next)[n], INVALID);
   718 
   719               while (next_bucket != *_highest) {
   720                 --_highest;
   721               }
   722               if (_highest == _sets.back().begin()) {
   723                 _sets.back().push_front(bucket_num);
   724                 _dormant[bucket_num] = false;
   725                 _first[bucket_num] = _last[bucket_num] = INVALID;
   726                 ++bucket_num;
   727               }
   728               --_highest;
   729 
   730               _bucket->set(n, *_highest);
   731               _next->set(n, _first[*_highest]);
   732               if (_first[*_highest] != INVALID) {
   733                 _prev->set(_first[*_highest], n);
   734               } else {
   735                 _last[*_highest] = n;
   736               }
   737               _first[*_highest] = n;
   738             }
   739           } else {
   740 
   741             deactivate(n);
   742             if (!(*_active)[_first[*_highest]]) {
   743               ++_highest;
   744               if (_highest != _sets.back().end() &&
   745                   !(*_active)[_first[*_highest]]) {
   746                 _highest = _sets.back().end();
   747               }
   748             }
   749           }
   750         }
   751 
   752         if ((*_excess)[target] < _min_cut) {
   753           _min_cut = (*_excess)[target];
   754           for (NodeIt i(_graph); i != INVALID; ++i) {
   755             _min_cut_map->set(i, false);
   756           }
   757           for (std::list<int>::iterator it = _sets.back().begin();
   758                it != _sets.back().end(); ++it) {
   759             Node n = _first[*it];
   760             while (n != INVALID) {
   761               _min_cut_map->set(n, true);
   762               n = (*_next)[n];
   763             }
   764           }
   765         }
   766 
   767         {
   768           Node new_target;
   769           if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
   770             if ((*_next)[target] == INVALID) {
   771               _last[(*_bucket)[target]] = (*_prev)[target];
   772               new_target = (*_prev)[target];
   773             } else {
   774               _prev->set((*_next)[target], (*_prev)[target]);
   775               new_target = (*_next)[target];
   776             }
   777             if ((*_prev)[target] == INVALID) {
   778               _first[(*_bucket)[target]] = (*_next)[target];
   779             } else {
   780               _next->set((*_prev)[target], (*_next)[target]);
   781             }
   782           } else {
   783             _sets.back().pop_back();
   784             if (_sets.back().empty()) {
   785               _sets.pop_back();
   786               if (_sets.empty())
   787                 break;
   788               for (std::list<int>::iterator it = _sets.back().begin();
   789                    it != _sets.back().end(); ++it) {
   790                 _dormant[*it] = false;
   791               }
   792             }
   793             new_target = _last[_sets.back().back()];
   794           }
   795 
   796           _bucket->set(target, 0);
   797 
   798           _source_set->set(target, true);
   799           for (InArcIt a(_graph, target); a != INVALID; ++a) {
   800             Value rem = (*_capacity)[a] - (*_flow)[a];
   801             if (!_tolerance.positive(rem)) continue;
   802             Node v = _graph.source(a);
   803             if (!(*_active)[v] && !(*_source_set)[v]) {
   804               activate(v);
   805             }
   806             _excess->set(v, (*_excess)[v] + rem);
   807             _flow->set(a, (*_capacity)[a]);
   808           }
   809 
   810           for (OutArcIt a(_graph, target); a != INVALID; ++a) {
   811             Value rem = (*_flow)[a];
   812             if (!_tolerance.positive(rem)) continue;
   813             Node v = _graph.target(a);
   814             if (!(*_active)[v] && !(*_source_set)[v]) {
   815               activate(v);
   816             }
   817             _excess->set(v, (*_excess)[v] + rem);
   818             _flow->set(a, 0);
   819           }
   820 
   821           target = new_target;
   822           if ((*_active)[target]) {
   823             deactivate(target);
   824           }
   825 
   826           _highest = _sets.back().begin();
   827           while (_highest != _sets.back().end() &&
   828                  !(*_active)[_first[*_highest]]) {
   829             ++_highest;
   830           }
   831         }
   832       }
   833     }
   834 
   835   public:
   836 
   837     /// \name Execution control
   838     /// The simplest way to execute the algorithm is to use
   839     /// one of the member functions called \c run(...).
   840     /// \n
   841     /// If you need more control on the execution,
   842     /// first you must call \ref init(), then the \ref calculateIn() or
   843     /// \ref calculateIn() functions.
   844 
   845     /// @{
   846 
   847     /// \brief Initializes the internal data structures.
   848     ///
   849     /// Initializes the internal data structures. It creates
   850     /// the maps, residual graph adaptors and some bucket structures
   851     /// for the algorithm.
   852     void init() {
   853       init(NodeIt(_graph));
   854     }
   855 
   856     /// \brief Initializes the internal data structures.
   857     ///
   858     /// Initializes the internal data structures. It creates
   859     /// the maps, residual graph adaptor and some bucket structures
   860     /// for the algorithm. Node \c source  is used as the push-relabel
   861     /// algorithm's source.
   862     void init(const Node& source) {
   863       _source = source;
   864 
   865       _node_num = countNodes(_graph);
   866 
   867       _first.resize(_node_num + 1);
   868       _last.resize(_node_num + 1);
   869 
   870       _dormant.resize(_node_num + 1);
   871 
   872       if (!_flow) {
   873         _flow = new FlowMap(_graph);
   874       }
   875       if (!_next) {
   876         _next = new typename Digraph::template NodeMap<Node>(_graph);
   877       }
   878       if (!_prev) {
   879         _prev = new typename Digraph::template NodeMap<Node>(_graph);
   880       }
   881       if (!_active) {
   882         _active = new typename Digraph::template NodeMap<bool>(_graph);
   883       }
   884       if (!_bucket) {
   885         _bucket = new typename Digraph::template NodeMap<int>(_graph);
   886       }
   887       if (!_excess) {
   888         _excess = new ExcessMap(_graph);
   889       }
   890       if (!_source_set) {
   891         _source_set = new SourceSetMap(_graph);
   892       }
   893       if (!_min_cut_map) {
   894         _min_cut_map = new MinCutMap(_graph);
   895       }
   896 
   897       _min_cut = std::numeric_limits<Value>::max();
   898     }
   899 
   900 
   901     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   902     /// source-side.
   903     ///
   904     /// Calculates a minimum cut with \f$ source \f$ on the
   905     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
   906     /// \in X \f$ and minimal out-degree).
   907     void calculateOut() {
   908       findMinCutOut();
   909     }
   910 
   911     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   912     /// target-side.
   913     ///
   914     /// Calculates a minimum cut with \f$ source \f$ on the
   915     /// target-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
   916     /// \in X \f$ and minimal out-degree).
   917     void calculateIn() {
   918       findMinCutIn();
   919     }
   920 
   921 
   922     /// \brief Runs the algorithm.
   923     ///
   924     /// Runs the algorithm. It finds nodes \c source and \c target
   925     /// arbitrarily and then calls \ref init(), \ref calculateOut()
   926     /// and \ref calculateIn().
   927     void run() {
   928       init();
   929       calculateOut();
   930       calculateIn();
   931     }
   932 
   933     /// \brief Runs the algorithm.
   934     ///
   935     /// Runs the algorithm. It uses the given \c source node, finds a
   936     /// proper \c target and then calls the \ref init(), \ref
   937     /// calculateOut() and \ref calculateIn().
   938     void run(const Node& s) {
   939       init(s);
   940       calculateOut();
   941       calculateIn();
   942     }
   943 
   944     /// @}
   945 
   946     /// \name Query Functions
   947     /// The result of the %HaoOrlin algorithm
   948     /// can be obtained using these functions.
   949     /// \n
   950     /// Before using these functions, either \ref run(), \ref
   951     /// calculateOut() or \ref calculateIn() must be called.
   952 
   953     /// @{
   954 
   955     /// \brief Returns the value of the minimum value cut.
   956     ///
   957     /// Returns the value of the minimum value cut.
   958     Value minCutValue() const {
   959       return _min_cut;
   960     }
   961 
   962 
   963     /// \brief Returns a minimum cut.
   964     ///
   965     /// Sets \c nodeMap to the characteristic vector of a minimum
   966     /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
   967     /// with minimal out-degree (i.e. \c nodeMap will be true exactly
   968     /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
   969     /// bool-valued node-map.
   970     template <typename NodeMap>
   971     Value minCutMap(NodeMap& nodeMap) const {
   972       for (NodeIt it(_graph); it != INVALID; ++it) {
   973         nodeMap.set(it, (*_min_cut_map)[it]);
   974       }
   975       return _min_cut;
   976     }
   977 
   978     /// @}
   979 
   980   }; //class HaoOrlin
   981 
   982 
   983 } //namespace lemon
   984 
   985 #endif //LEMON_HAO_ORLIN_H