lemon/hao_orlin.h
author Alpar Juttner <alpar@cs.elte.hu>
Mon, 23 Feb 2009 11:48:47 +0000
changeset 525 635a8375227d
parent 412 7030149efed2
child 559 c5fd2d996909
permissions -rw-r--r--
dimacs.h reads MAT files to both dir and undir graphs (#231)
     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_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 = 0;
   229       std::vector<Node> queue(_node_num);
   230       int qfirst = 0, qlast = 0, qsep = 0;
   231 
   232       {
   233         typename Digraph::template NodeMap<bool> reached(_graph, false);
   234 
   235         reached.set(_source, true);
   236         bool first_set = true;
   237 
   238         for (NodeIt t(_graph); t != INVALID; ++t) {
   239           if (reached[t]) continue;
   240           _sets.push_front(std::list<int>());
   241 
   242           queue[qlast++] = t;
   243           reached.set(t, true);
   244 
   245           while (qfirst != qlast) {
   246             if (qsep == qfirst) {
   247               ++bucket_num;
   248               _sets.front().push_front(bucket_num);
   249               _dormant[bucket_num] = !first_set;
   250               _first[bucket_num] = _last[bucket_num] = INVALID;
   251               qsep = qlast;
   252             }
   253 
   254             Node n = queue[qfirst++];
   255             addItem(n, bucket_num);
   256 
   257             for (InArcIt a(_graph, n); a != INVALID; ++a) {
   258               Node u = _graph.source(a);
   259               if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
   260                 reached.set(u, true);
   261                 queue[qlast++] = u;
   262               }
   263             }
   264           }
   265           first_set = false;
   266         }
   267 
   268         ++bucket_num;
   269         _bucket->set(_source, 0);
   270         _dormant[0] = true;
   271       }
   272       _source_set->set(_source, true);
   273 
   274       Node target = _last[_sets.back().back()];
   275       {
   276         for (OutArcIt a(_graph, _source); a != INVALID; ++a) {
   277           if (_tolerance.positive((*_capacity)[a])) {
   278             Node u = _graph.target(a);
   279             _flow->set(a, (*_capacity)[a]);
   280             _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
   281             if (!(*_active)[u] && u != _source) {
   282               activate(u);
   283             }
   284           }
   285         }
   286 
   287         if ((*_active)[target]) {
   288           deactivate(target);
   289         }
   290 
   291         _highest = _sets.back().begin();
   292         while (_highest != _sets.back().end() &&
   293                !(*_active)[_first[*_highest]]) {
   294           ++_highest;
   295         }
   296       }
   297 
   298       while (true) {
   299         while (_highest != _sets.back().end()) {
   300           Node n = _first[*_highest];
   301           Value excess = (*_excess)[n];
   302           int next_bucket = _node_num;
   303 
   304           int under_bucket;
   305           if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
   306             under_bucket = -1;
   307           } else {
   308             under_bucket = *(++std::list<int>::iterator(_highest));
   309           }
   310 
   311           for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   312             Node v = _graph.target(a);
   313             if (_dormant[(*_bucket)[v]]) continue;
   314             Value rem = (*_capacity)[a] - (*_flow)[a];
   315             if (!_tolerance.positive(rem)) continue;
   316             if ((*_bucket)[v] == under_bucket) {
   317               if (!(*_active)[v] && v != target) {
   318                 activate(v);
   319               }
   320               if (!_tolerance.less(rem, excess)) {
   321                 _flow->set(a, (*_flow)[a] + excess);
   322                 _excess->set(v, (*_excess)[v] + excess);
   323                 excess = 0;
   324                 goto no_more_push;
   325               } else {
   326                 excess -= rem;
   327                 _excess->set(v, (*_excess)[v] + rem);
   328                 _flow->set(a, (*_capacity)[a]);
   329               }
   330             } else if (next_bucket > (*_bucket)[v]) {
   331               next_bucket = (*_bucket)[v];
   332             }
   333           }
   334 
   335           for (InArcIt a(_graph, n); a != INVALID; ++a) {
   336             Node v = _graph.source(a);
   337             if (_dormant[(*_bucket)[v]]) continue;
   338             Value rem = (*_flow)[a];
   339             if (!_tolerance.positive(rem)) continue;
   340             if ((*_bucket)[v] == under_bucket) {
   341               if (!(*_active)[v] && v != target) {
   342                 activate(v);
   343               }
   344               if (!_tolerance.less(rem, excess)) {
   345                 _flow->set(a, (*_flow)[a] - excess);
   346                 _excess->set(v, (*_excess)[v] + excess);
   347                 excess = 0;
   348                 goto no_more_push;
   349               } else {
   350                 excess -= rem;
   351                 _excess->set(v, (*_excess)[v] + rem);
   352                 _flow->set(a, 0);
   353               }
   354             } else if (next_bucket > (*_bucket)[v]) {
   355               next_bucket = (*_bucket)[v];
   356             }
   357           }
   358 
   359         no_more_push:
   360 
   361           _excess->set(n, excess);
   362 
   363           if (excess != 0) {
   364             if ((*_next)[n] == INVALID) {
   365               typename std::list<std::list<int> >::iterator new_set =
   366                 _sets.insert(--_sets.end(), std::list<int>());
   367               new_set->splice(new_set->end(), _sets.back(),
   368                               _sets.back().begin(), ++_highest);
   369               for (std::list<int>::iterator it = new_set->begin();
   370                    it != new_set->end(); ++it) {
   371                 _dormant[*it] = true;
   372               }
   373               while (_highest != _sets.back().end() &&
   374                      !(*_active)[_first[*_highest]]) {
   375                 ++_highest;
   376               }
   377             } else if (next_bucket == _node_num) {
   378               _first[(*_bucket)[n]] = (*_next)[n];
   379               _prev->set((*_next)[n], INVALID);
   380 
   381               std::list<std::list<int> >::iterator new_set =
   382                 _sets.insert(--_sets.end(), std::list<int>());
   383 
   384               new_set->push_front(bucket_num);
   385               _bucket->set(n, bucket_num);
   386               _first[bucket_num] = _last[bucket_num] = n;
   387               _next->set(n, INVALID);
   388               _prev->set(n, INVALID);
   389               _dormant[bucket_num] = true;
   390               ++bucket_num;
   391 
   392               while (_highest != _sets.back().end() &&
   393                      !(*_active)[_first[*_highest]]) {
   394                 ++_highest;
   395               }
   396             } else {
   397               _first[*_highest] = (*_next)[n];
   398               _prev->set((*_next)[n], INVALID);
   399 
   400               while (next_bucket != *_highest) {
   401                 --_highest;
   402               }
   403 
   404               if (_highest == _sets.back().begin()) {
   405                 _sets.back().push_front(bucket_num);
   406                 _dormant[bucket_num] = false;
   407                 _first[bucket_num] = _last[bucket_num] = INVALID;
   408                 ++bucket_num;
   409               }
   410               --_highest;
   411 
   412               _bucket->set(n, *_highest);
   413               _next->set(n, _first[*_highest]);
   414               if (_first[*_highest] != INVALID) {
   415                 _prev->set(_first[*_highest], n);
   416               } else {
   417                 _last[*_highest] = n;
   418               }
   419               _first[*_highest] = n;
   420             }
   421           } else {
   422 
   423             deactivate(n);
   424             if (!(*_active)[_first[*_highest]]) {
   425               ++_highest;
   426               if (_highest != _sets.back().end() &&
   427                   !(*_active)[_first[*_highest]]) {
   428                 _highest = _sets.back().end();
   429               }
   430             }
   431           }
   432         }
   433 
   434         if ((*_excess)[target] < _min_cut) {
   435           _min_cut = (*_excess)[target];
   436           for (NodeIt i(_graph); i != INVALID; ++i) {
   437             _min_cut_map->set(i, true);
   438           }
   439           for (std::list<int>::iterator it = _sets.back().begin();
   440                it != _sets.back().end(); ++it) {
   441             Node n = _first[*it];
   442             while (n != INVALID) {
   443               _min_cut_map->set(n, false);
   444               n = (*_next)[n];
   445             }
   446           }
   447         }
   448 
   449         {
   450           Node new_target;
   451           if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
   452             if ((*_next)[target] == INVALID) {
   453               _last[(*_bucket)[target]] = (*_prev)[target];
   454               new_target = (*_prev)[target];
   455             } else {
   456               _prev->set((*_next)[target], (*_prev)[target]);
   457               new_target = (*_next)[target];
   458             }
   459             if ((*_prev)[target] == INVALID) {
   460               _first[(*_bucket)[target]] = (*_next)[target];
   461             } else {
   462               _next->set((*_prev)[target], (*_next)[target]);
   463             }
   464           } else {
   465             _sets.back().pop_back();
   466             if (_sets.back().empty()) {
   467               _sets.pop_back();
   468               if (_sets.empty())
   469                 break;
   470               for (std::list<int>::iterator it = _sets.back().begin();
   471                    it != _sets.back().end(); ++it) {
   472                 _dormant[*it] = false;
   473               }
   474             }
   475             new_target = _last[_sets.back().back()];
   476           }
   477 
   478           _bucket->set(target, 0);
   479 
   480           _source_set->set(target, true);
   481           for (OutArcIt a(_graph, target); a != INVALID; ++a) {
   482             Value rem = (*_capacity)[a] - (*_flow)[a];
   483             if (!_tolerance.positive(rem)) continue;
   484             Node v = _graph.target(a);
   485             if (!(*_active)[v] && !(*_source_set)[v]) {
   486               activate(v);
   487             }
   488             _excess->set(v, (*_excess)[v] + rem);
   489             _flow->set(a, (*_capacity)[a]);
   490           }
   491 
   492           for (InArcIt a(_graph, target); a != INVALID; ++a) {
   493             Value rem = (*_flow)[a];
   494             if (!_tolerance.positive(rem)) continue;
   495             Node v = _graph.source(a);
   496             if (!(*_active)[v] && !(*_source_set)[v]) {
   497               activate(v);
   498             }
   499             _excess->set(v, (*_excess)[v] + rem);
   500             _flow->set(a, 0);
   501           }
   502 
   503           target = new_target;
   504           if ((*_active)[target]) {
   505             deactivate(target);
   506           }
   507 
   508           _highest = _sets.back().begin();
   509           while (_highest != _sets.back().end() &&
   510                  !(*_active)[_first[*_highest]]) {
   511             ++_highest;
   512           }
   513         }
   514       }
   515     }
   516 
   517     void findMinCutIn() {
   518 
   519       for (NodeIt n(_graph); n != INVALID; ++n) {
   520         _excess->set(n, 0);
   521       }
   522 
   523       for (ArcIt a(_graph); a != INVALID; ++a) {
   524         _flow->set(a, 0);
   525       }
   526 
   527       int bucket_num = 0;
   528       std::vector<Node> queue(_node_num);
   529       int qfirst = 0, qlast = 0, qsep = 0;
   530 
   531       {
   532         typename Digraph::template NodeMap<bool> reached(_graph, false);
   533 
   534         reached.set(_source, true);
   535 
   536         bool first_set = true;
   537 
   538         for (NodeIt t(_graph); t != INVALID; ++t) {
   539           if (reached[t]) continue;
   540           _sets.push_front(std::list<int>());
   541 
   542           queue[qlast++] = t;
   543           reached.set(t, true);
   544 
   545           while (qfirst != qlast) {
   546             if (qsep == qfirst) {
   547               ++bucket_num;
   548               _sets.front().push_front(bucket_num);
   549               _dormant[bucket_num] = !first_set;
   550               _first[bucket_num] = _last[bucket_num] = INVALID;
   551               qsep = qlast;
   552             }
   553 
   554             Node n = queue[qfirst++];
   555             addItem(n, bucket_num);
   556 
   557             for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   558               Node u = _graph.target(a);
   559               if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
   560                 reached.set(u, true);
   561                 queue[qlast++] = u;
   562               }
   563             }
   564           }
   565           first_set = false;
   566         }
   567 
   568         ++bucket_num;
   569         _bucket->set(_source, 0);
   570         _dormant[0] = true;
   571       }
   572       _source_set->set(_source, true);
   573 
   574       Node target = _last[_sets.back().back()];
   575       {
   576         for (InArcIt a(_graph, _source); a != INVALID; ++a) {
   577           if (_tolerance.positive((*_capacity)[a])) {
   578             Node u = _graph.source(a);
   579             _flow->set(a, (*_capacity)[a]);
   580             _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
   581             if (!(*_active)[u] && u != _source) {
   582               activate(u);
   583             }
   584           }
   585         }
   586         if ((*_active)[target]) {
   587           deactivate(target);
   588         }
   589 
   590         _highest = _sets.back().begin();
   591         while (_highest != _sets.back().end() &&
   592                !(*_active)[_first[*_highest]]) {
   593           ++_highest;
   594         }
   595       }
   596 
   597 
   598       while (true) {
   599         while (_highest != _sets.back().end()) {
   600           Node n = _first[*_highest];
   601           Value excess = (*_excess)[n];
   602           int next_bucket = _node_num;
   603 
   604           int under_bucket;
   605           if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
   606             under_bucket = -1;
   607           } else {
   608             under_bucket = *(++std::list<int>::iterator(_highest));
   609           }
   610 
   611           for (InArcIt a(_graph, n); a != INVALID; ++a) {
   612             Node v = _graph.source(a);
   613             if (_dormant[(*_bucket)[v]]) continue;
   614             Value rem = (*_capacity)[a] - (*_flow)[a];
   615             if (!_tolerance.positive(rem)) continue;
   616             if ((*_bucket)[v] == under_bucket) {
   617               if (!(*_active)[v] && v != target) {
   618                 activate(v);
   619               }
   620               if (!_tolerance.less(rem, excess)) {
   621                 _flow->set(a, (*_flow)[a] + excess);
   622                 _excess->set(v, (*_excess)[v] + excess);
   623                 excess = 0;
   624                 goto no_more_push;
   625               } else {
   626                 excess -= rem;
   627                 _excess->set(v, (*_excess)[v] + rem);
   628                 _flow->set(a, (*_capacity)[a]);
   629               }
   630             } else if (next_bucket > (*_bucket)[v]) {
   631               next_bucket = (*_bucket)[v];
   632             }
   633           }
   634 
   635           for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   636             Node v = _graph.target(a);
   637             if (_dormant[(*_bucket)[v]]) continue;
   638             Value rem = (*_flow)[a];
   639             if (!_tolerance.positive(rem)) continue;
   640             if ((*_bucket)[v] == under_bucket) {
   641               if (!(*_active)[v] && v != target) {
   642                 activate(v);
   643               }
   644               if (!_tolerance.less(rem, excess)) {
   645                 _flow->set(a, (*_flow)[a] - excess);
   646                 _excess->set(v, (*_excess)[v] + excess);
   647                 excess = 0;
   648                 goto no_more_push;
   649               } else {
   650                 excess -= rem;
   651                 _excess->set(v, (*_excess)[v] + rem);
   652                 _flow->set(a, 0);
   653               }
   654             } else if (next_bucket > (*_bucket)[v]) {
   655               next_bucket = (*_bucket)[v];
   656             }
   657           }
   658 
   659         no_more_push:
   660 
   661           _excess->set(n, excess);
   662 
   663           if (excess != 0) {
   664             if ((*_next)[n] == INVALID) {
   665               typename std::list<std::list<int> >::iterator new_set =
   666                 _sets.insert(--_sets.end(), std::list<int>());
   667               new_set->splice(new_set->end(), _sets.back(),
   668                               _sets.back().begin(), ++_highest);
   669               for (std::list<int>::iterator it = new_set->begin();
   670                    it != new_set->end(); ++it) {
   671                 _dormant[*it] = true;
   672               }
   673               while (_highest != _sets.back().end() &&
   674                      !(*_active)[_first[*_highest]]) {
   675                 ++_highest;
   676               }
   677             } else if (next_bucket == _node_num) {
   678               _first[(*_bucket)[n]] = (*_next)[n];
   679               _prev->set((*_next)[n], INVALID);
   680 
   681               std::list<std::list<int> >::iterator new_set =
   682                 _sets.insert(--_sets.end(), std::list<int>());
   683 
   684               new_set->push_front(bucket_num);
   685               _bucket->set(n, bucket_num);
   686               _first[bucket_num] = _last[bucket_num] = n;
   687               _next->set(n, INVALID);
   688               _prev->set(n, INVALID);
   689               _dormant[bucket_num] = true;
   690               ++bucket_num;
   691 
   692               while (_highest != _sets.back().end() &&
   693                      !(*_active)[_first[*_highest]]) {
   694                 ++_highest;
   695               }
   696             } else {
   697               _first[*_highest] = (*_next)[n];
   698               _prev->set((*_next)[n], INVALID);
   699 
   700               while (next_bucket != *_highest) {
   701                 --_highest;
   702               }
   703               if (_highest == _sets.back().begin()) {
   704                 _sets.back().push_front(bucket_num);
   705                 _dormant[bucket_num] = false;
   706                 _first[bucket_num] = _last[bucket_num] = INVALID;
   707                 ++bucket_num;
   708               }
   709               --_highest;
   710 
   711               _bucket->set(n, *_highest);
   712               _next->set(n, _first[*_highest]);
   713               if (_first[*_highest] != INVALID) {
   714                 _prev->set(_first[*_highest], n);
   715               } else {
   716                 _last[*_highest] = n;
   717               }
   718               _first[*_highest] = n;
   719             }
   720           } else {
   721 
   722             deactivate(n);
   723             if (!(*_active)[_first[*_highest]]) {
   724               ++_highest;
   725               if (_highest != _sets.back().end() &&
   726                   !(*_active)[_first[*_highest]]) {
   727                 _highest = _sets.back().end();
   728               }
   729             }
   730           }
   731         }
   732 
   733         if ((*_excess)[target] < _min_cut) {
   734           _min_cut = (*_excess)[target];
   735           for (NodeIt i(_graph); i != INVALID; ++i) {
   736             _min_cut_map->set(i, false);
   737           }
   738           for (std::list<int>::iterator it = _sets.back().begin();
   739                it != _sets.back().end(); ++it) {
   740             Node n = _first[*it];
   741             while (n != INVALID) {
   742               _min_cut_map->set(n, true);
   743               n = (*_next)[n];
   744             }
   745           }
   746         }
   747 
   748         {
   749           Node new_target;
   750           if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
   751             if ((*_next)[target] == INVALID) {
   752               _last[(*_bucket)[target]] = (*_prev)[target];
   753               new_target = (*_prev)[target];
   754             } else {
   755               _prev->set((*_next)[target], (*_prev)[target]);
   756               new_target = (*_next)[target];
   757             }
   758             if ((*_prev)[target] == INVALID) {
   759               _first[(*_bucket)[target]] = (*_next)[target];
   760             } else {
   761               _next->set((*_prev)[target], (*_next)[target]);
   762             }
   763           } else {
   764             _sets.back().pop_back();
   765             if (_sets.back().empty()) {
   766               _sets.pop_back();
   767               if (_sets.empty())
   768                 break;
   769               for (std::list<int>::iterator it = _sets.back().begin();
   770                    it != _sets.back().end(); ++it) {
   771                 _dormant[*it] = false;
   772               }
   773             }
   774             new_target = _last[_sets.back().back()];
   775           }
   776 
   777           _bucket->set(target, 0);
   778 
   779           _source_set->set(target, true);
   780           for (InArcIt a(_graph, target); a != INVALID; ++a) {
   781             Value rem = (*_capacity)[a] - (*_flow)[a];
   782             if (!_tolerance.positive(rem)) continue;
   783             Node v = _graph.source(a);
   784             if (!(*_active)[v] && !(*_source_set)[v]) {
   785               activate(v);
   786             }
   787             _excess->set(v, (*_excess)[v] + rem);
   788             _flow->set(a, (*_capacity)[a]);
   789           }
   790 
   791           for (OutArcIt a(_graph, target); a != INVALID; ++a) {
   792             Value rem = (*_flow)[a];
   793             if (!_tolerance.positive(rem)) continue;
   794             Node v = _graph.target(a);
   795             if (!(*_active)[v] && !(*_source_set)[v]) {
   796               activate(v);
   797             }
   798             _excess->set(v, (*_excess)[v] + rem);
   799             _flow->set(a, 0);
   800           }
   801 
   802           target = new_target;
   803           if ((*_active)[target]) {
   804             deactivate(target);
   805           }
   806 
   807           _highest = _sets.back().begin();
   808           while (_highest != _sets.back().end() &&
   809                  !(*_active)[_first[*_highest]]) {
   810             ++_highest;
   811           }
   812         }
   813       }
   814     }
   815 
   816   public:
   817 
   818     /// \name Execution control
   819     /// The simplest way to execute the algorithm is to use
   820     /// one of the member functions called \c run(...).
   821     /// \n
   822     /// If you need more control on the execution,
   823     /// first you must call \ref init(), then the \ref calculateIn() or
   824     /// \ref calculateOut() functions.
   825 
   826     /// @{
   827 
   828     /// \brief Initializes the internal data structures.
   829     ///
   830     /// Initializes the internal data structures. It creates
   831     /// the maps, residual graph adaptors and some bucket structures
   832     /// for the algorithm.
   833     void init() {
   834       init(NodeIt(_graph));
   835     }
   836 
   837     /// \brief Initializes the internal data structures.
   838     ///
   839     /// Initializes the internal data structures. It creates
   840     /// the maps, residual graph adaptor and some bucket structures
   841     /// for the algorithm. Node \c source  is used as the push-relabel
   842     /// algorithm's source.
   843     void init(const Node& source) {
   844       _source = source;
   845 
   846       _node_num = countNodes(_graph);
   847 
   848       _first.resize(_node_num);
   849       _last.resize(_node_num);
   850 
   851       _dormant.resize(_node_num);
   852 
   853       if (!_flow) {
   854         _flow = new FlowMap(_graph);
   855       }
   856       if (!_next) {
   857         _next = new typename Digraph::template NodeMap<Node>(_graph);
   858       }
   859       if (!_prev) {
   860         _prev = new typename Digraph::template NodeMap<Node>(_graph);
   861       }
   862       if (!_active) {
   863         _active = new typename Digraph::template NodeMap<bool>(_graph);
   864       }
   865       if (!_bucket) {
   866         _bucket = new typename Digraph::template NodeMap<int>(_graph);
   867       }
   868       if (!_excess) {
   869         _excess = new ExcessMap(_graph);
   870       }
   871       if (!_source_set) {
   872         _source_set = new SourceSetMap(_graph);
   873       }
   874       if (!_min_cut_map) {
   875         _min_cut_map = new MinCutMap(_graph);
   876       }
   877 
   878       _min_cut = std::numeric_limits<Value>::max();
   879     }
   880 
   881 
   882     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   883     /// source-side.
   884     ///
   885     /// Calculates a minimum cut with \f$ source \f$ on the
   886     /// source-side (i.e. a set \f$ X\subsetneq V \f$ with
   887     /// \f$ source \in X \f$ and minimal out-degree).
   888     void calculateOut() {
   889       findMinCutOut();
   890     }
   891 
   892     /// \brief Calculates a minimum cut with \f$ source \f$ on the
   893     /// target-side.
   894     ///
   895     /// Calculates a minimum cut with \f$ source \f$ on the
   896     /// target-side (i.e. a set \f$ X\subsetneq V \f$ with
   897     /// \f$ source \in X \f$ and minimal out-degree).
   898     void calculateIn() {
   899       findMinCutIn();
   900     }
   901 
   902 
   903     /// \brief Runs the algorithm.
   904     ///
   905     /// Runs the algorithm. It finds nodes \c source and \c target
   906     /// arbitrarily and then calls \ref init(), \ref calculateOut()
   907     /// and \ref calculateIn().
   908     void run() {
   909       init();
   910       calculateOut();
   911       calculateIn();
   912     }
   913 
   914     /// \brief Runs the algorithm.
   915     ///
   916     /// Runs the algorithm. It uses the given \c source node, finds a
   917     /// proper \c target and then calls the \ref init(), \ref
   918     /// calculateOut() and \ref calculateIn().
   919     void run(const Node& s) {
   920       init(s);
   921       calculateOut();
   922       calculateIn();
   923     }
   924 
   925     /// @}
   926 
   927     /// \name Query Functions
   928     /// The result of the %HaoOrlin algorithm
   929     /// can be obtained using these functions.
   930     /// \n
   931     /// Before using these functions, either \ref run(), \ref
   932     /// calculateOut() or \ref calculateIn() must be called.
   933 
   934     /// @{
   935 
   936     /// \brief Returns the value of the minimum value cut.
   937     ///
   938     /// Returns the value of the minimum value cut.
   939     Value minCutValue() const {
   940       return _min_cut;
   941     }
   942 
   943 
   944     /// \brief Returns a minimum cut.
   945     ///
   946     /// Sets \c nodeMap to the characteristic vector of a minimum
   947     /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
   948     /// with minimal out-degree (i.e. \c nodeMap will be true exactly
   949     /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
   950     /// bool-valued node-map.
   951     template <typename NodeMap>
   952     Value minCutMap(NodeMap& nodeMap) const {
   953       for (NodeIt it(_graph); it != INVALID; ++it) {
   954         nodeMap.set(it, (*_min_cut_map)[it]);
   955       }
   956       return _min_cut;
   957     }
   958 
   959     /// @}
   960 
   961   }; //class HaoOrlin
   962 
   963 
   964 } //namespace lemon
   965 
   966 #endif //LEMON_HAO_ORLIN_H