Port Hao-Orlin algorithm from SVN -r3509 (#58)
authorBalazs Dezso <deba@inf.elte.hu>
Mon, 01 Dec 2008 23:12:16 +0100
changeset 409b8ce15103485
parent 404 59d3aa4f921f
child 410 eac19fb31a09
Port Hao-Orlin algorithm from SVN -r3509 (#58)
lemon/Makefile.am
lemon/hao_orlin.h
     1.1 --- a/lemon/Makefile.am	Mon Dec 01 14:18:40 2008 +0000
     1.2 +++ b/lemon/Makefile.am	Mon Dec 01 23:12:16 2008 +0100
     1.3 @@ -36,6 +36,7 @@
     1.4          lemon/grid_graph.h \
     1.5  	lemon/hypercube_graph.h \
     1.6  	lemon/kruskal.h \
     1.7 +	lemon/hao_orlin.h \
     1.8  	lemon/lgf_reader.h \
     1.9  	lemon/lgf_writer.h \
    1.10  	lemon/list_graph.h \
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/lemon/hao_orlin.h	Mon Dec 01 23:12:16 2008 +0100
     2.3 @@ -0,0 +1,985 @@
     2.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2.5 + *
     2.6 + * This file is a part of LEMON, a generic C++ optimization library.
     2.7 + *
     2.8 + * Copyright (C) 2003-2008
     2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    2.11 + *
    2.12 + * Permission to use, modify and distribute this software is granted
    2.13 + * provided that this copyright notice appears in all copies. For
    2.14 + * precise terms see the accompanying LICENSE file.
    2.15 + *
    2.16 + * This software is provided "AS IS" with no warranty of any kind,
    2.17 + * express or implied, and with no claim as to its suitability for any
    2.18 + * purpose.
    2.19 + *
    2.20 + */
    2.21 +
    2.22 +#ifndef LEMON_HAO_ORLIN_H
    2.23 +#define LEMON_HAO_ORLIN_H
    2.24 +
    2.25 +#include <vector>
    2.26 +#include <list>
    2.27 +#include <limits>
    2.28 +
    2.29 +#include <lemon/maps.h>
    2.30 +#include <lemon/core.h>
    2.31 +#include <lemon/tolerance.h>
    2.32 +
    2.33 +/// \file
    2.34 +/// \ingroup min_cut
    2.35 +/// \brief Implementation of the Hao-Orlin algorithm.
    2.36 +///
    2.37 +/// Implementation of the Hao-Orlin algorithm class for testing network
    2.38 +/// reliability.
    2.39 +
    2.40 +namespace lemon {
    2.41 +
    2.42 +  /// \ingroup min_cut
    2.43 +  ///
    2.44 +  /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
    2.45 +  ///
    2.46 +  /// Hao-Orlin calculates a minimum cut in a directed graph
    2.47 +  /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and
    2.48 +  /// consists of two phases: in the first phase it determines a
    2.49 +  /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
    2.50 +  /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal
    2.51 +  /// out-degree) and in the second phase it determines a minimum cut
    2.52 +  /// with \f$ source \f$ on the sink-side (i.e. a set
    2.53 +  /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal
    2.54 +  /// out-degree). Obviously, the smaller of these two cuts will be a
    2.55 +  /// minimum cut of \f$ D \f$. The algorithm is a modified
    2.56 +  /// push-relabel preflow algorithm and our implementation calculates
    2.57 +  /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
    2.58 +  /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
    2.59 +  /// purpose of such algorithm is testing network reliability. For an
    2.60 +  /// undirected graph you can run just the first phase of the
    2.61 +  /// algorithm or you can use the algorithm of Nagamochi and Ibaraki
    2.62 +  /// which solves the undirected problem in
    2.63 +  /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the
    2.64 +  /// NagamochiIbaraki algorithm class.
    2.65 +  ///
    2.66 +  /// \param _Digraph is the graph type of the algorithm.
    2.67 +  /// \param _CapacityMap is an edge map of capacities which should
    2.68 +  /// be any numreric type. The default type is _Digraph::ArcMap<int>.
    2.69 +  /// \param _Tolerance is the handler of the inexact computation. The
    2.70 +  /// default type for this is Tolerance<CapacityMap::Value>.
    2.71 +#ifdef DOXYGEN
    2.72 +  template <typename _Digraph, typename _CapacityMap, typename _Tolerance>
    2.73 +#else
    2.74 +  template <typename _Digraph,
    2.75 +            typename _CapacityMap = typename _Digraph::template ArcMap<int>,
    2.76 +            typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
    2.77 +#endif
    2.78 +  class HaoOrlin {
    2.79 +  private:
    2.80 +
    2.81 +    typedef _Digraph Digraph;
    2.82 +    typedef _CapacityMap CapacityMap;
    2.83 +    typedef _Tolerance Tolerance;
    2.84 +
    2.85 +    typedef typename CapacityMap::Value Value;
    2.86 +
    2.87 +    TEMPLATE_GRAPH_TYPEDEFS(Digraph);
    2.88 +
    2.89 +    const Digraph& _graph;
    2.90 +    const CapacityMap* _capacity;
    2.91 +
    2.92 +    typedef typename Digraph::template ArcMap<Value> FlowMap;
    2.93 +    FlowMap* _flow;
    2.94 +
    2.95 +    Node _source;
    2.96 +
    2.97 +    int _node_num;
    2.98 +
    2.99 +    // Bucketing structure
   2.100 +    std::vector<Node> _first, _last;
   2.101 +    typename Digraph::template NodeMap<Node>* _next;
   2.102 +    typename Digraph::template NodeMap<Node>* _prev;
   2.103 +    typename Digraph::template NodeMap<bool>* _active;
   2.104 +    typename Digraph::template NodeMap<int>* _bucket;
   2.105 +
   2.106 +    std::vector<bool> _dormant;
   2.107 +
   2.108 +    std::list<std::list<int> > _sets;
   2.109 +    std::list<int>::iterator _highest;
   2.110 +
   2.111 +    typedef typename Digraph::template NodeMap<Value> ExcessMap;
   2.112 +    ExcessMap* _excess;
   2.113 +
   2.114 +    typedef typename Digraph::template NodeMap<bool> SourceSetMap;
   2.115 +    SourceSetMap* _source_set;
   2.116 +
   2.117 +    Value _min_cut;
   2.118 +
   2.119 +    typedef typename Digraph::template NodeMap<bool> MinCutMap;
   2.120 +    MinCutMap* _min_cut_map;
   2.121 +
   2.122 +    Tolerance _tolerance;
   2.123 +
   2.124 +  public:
   2.125 +
   2.126 +    /// \brief Constructor
   2.127 +    ///
   2.128 +    /// Constructor of the algorithm class.
   2.129 +    HaoOrlin(const Digraph& graph, const CapacityMap& capacity,
   2.130 +             const Tolerance& tolerance = Tolerance()) :
   2.131 +      _graph(graph), _capacity(&capacity), _flow(0), _source(),
   2.132 +      _node_num(), _first(), _last(), _next(0), _prev(0),
   2.133 +      _active(0), _bucket(0), _dormant(), _sets(), _highest(),
   2.134 +      _excess(0), _source_set(0), _min_cut(), _min_cut_map(0),
   2.135 +      _tolerance(tolerance) {}
   2.136 +
   2.137 +    ~HaoOrlin() {
   2.138 +      if (_min_cut_map) {
   2.139 +        delete _min_cut_map;
   2.140 +      }
   2.141 +      if (_source_set) {
   2.142 +        delete _source_set;
   2.143 +      }
   2.144 +      if (_excess) {
   2.145 +        delete _excess;
   2.146 +      }
   2.147 +      if (_next) {
   2.148 +        delete _next;
   2.149 +      }
   2.150 +      if (_prev) {
   2.151 +        delete _prev;
   2.152 +      }
   2.153 +      if (_active) {
   2.154 +        delete _active;
   2.155 +      }
   2.156 +      if (_bucket) {
   2.157 +        delete _bucket;
   2.158 +      }
   2.159 +      if (_flow) {
   2.160 +        delete _flow;
   2.161 +      }
   2.162 +    }
   2.163 +
   2.164 +  private:
   2.165 +
   2.166 +    void activate(const Node& i) {
   2.167 +      _active->set(i, true);
   2.168 +
   2.169 +      int bucket = (*_bucket)[i];
   2.170 +
   2.171 +      if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;
   2.172 +      //unlace
   2.173 +      _next->set((*_prev)[i], (*_next)[i]);
   2.174 +      if ((*_next)[i] != INVALID) {
   2.175 +        _prev->set((*_next)[i], (*_prev)[i]);
   2.176 +      } else {
   2.177 +        _last[bucket] = (*_prev)[i];
   2.178 +      }
   2.179 +      //lace
   2.180 +      _next->set(i, _first[bucket]);
   2.181 +      _prev->set(_first[bucket], i);
   2.182 +      _prev->set(i, INVALID);
   2.183 +      _first[bucket] = i;
   2.184 +    }
   2.185 +
   2.186 +    void deactivate(const Node& i) {
   2.187 +      _active->set(i, false);
   2.188 +      int bucket = (*_bucket)[i];
   2.189 +
   2.190 +      if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
   2.191 +
   2.192 +      //unlace
   2.193 +      _prev->set((*_next)[i], (*_prev)[i]);
   2.194 +      if ((*_prev)[i] != INVALID) {
   2.195 +        _next->set((*_prev)[i], (*_next)[i]);
   2.196 +      } else {
   2.197 +        _first[bucket] = (*_next)[i];
   2.198 +      }
   2.199 +      //lace
   2.200 +      _prev->set(i, _last[bucket]);
   2.201 +      _next->set(_last[bucket], i);
   2.202 +      _next->set(i, INVALID);
   2.203 +      _last[bucket] = i;
   2.204 +    }
   2.205 +
   2.206 +    void addItem(const Node& i, int bucket) {
   2.207 +      (*_bucket)[i] = bucket;
   2.208 +      if (_last[bucket] != INVALID) {
   2.209 +        _prev->set(i, _last[bucket]);
   2.210 +        _next->set(_last[bucket], i);
   2.211 +        _next->set(i, INVALID);
   2.212 +        _last[bucket] = i;
   2.213 +      } else {
   2.214 +        _prev->set(i, INVALID);
   2.215 +        _first[bucket] = i;
   2.216 +        _next->set(i, INVALID);
   2.217 +        _last[bucket] = i;
   2.218 +      }
   2.219 +    }
   2.220 +
   2.221 +    void findMinCutOut() {
   2.222 +
   2.223 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   2.224 +        _excess->set(n, 0);
   2.225 +      }
   2.226 +
   2.227 +      for (ArcIt a(_graph); a != INVALID; ++a) {
   2.228 +        _flow->set(a, 0);
   2.229 +      }
   2.230 +
   2.231 +      int bucket_num = 1;
   2.232 +
   2.233 +      {
   2.234 +        typename Digraph::template NodeMap<bool> reached(_graph, false);
   2.235 +
   2.236 +        reached.set(_source, true);
   2.237 +
   2.238 +        bool first_set = true;
   2.239 +
   2.240 +        for (NodeIt t(_graph); t != INVALID; ++t) {
   2.241 +          if (reached[t]) continue;
   2.242 +          _sets.push_front(std::list<int>());
   2.243 +          _sets.front().push_front(bucket_num);
   2.244 +          _dormant[bucket_num] = !first_set;
   2.245 +
   2.246 +          _bucket->set(t, bucket_num);
   2.247 +          _first[bucket_num] = _last[bucket_num] = t;
   2.248 +          _next->set(t, INVALID);
   2.249 +          _prev->set(t, INVALID);
   2.250 +
   2.251 +          ++bucket_num;
   2.252 +
   2.253 +          std::vector<Node> queue;
   2.254 +          queue.push_back(t);
   2.255 +          reached.set(t, true);
   2.256 +
   2.257 +          while (!queue.empty()) {
   2.258 +            _sets.front().push_front(bucket_num);
   2.259 +            _dormant[bucket_num] = !first_set;
   2.260 +            _first[bucket_num] = _last[bucket_num] = INVALID;
   2.261 +
   2.262 +            std::vector<Node> nqueue;
   2.263 +            for (int i = 0; i < int(queue.size()); ++i) {
   2.264 +              Node n = queue[i];
   2.265 +              for (InArcIt a(_graph, n); a != INVALID; ++a) {
   2.266 +                Node u = _graph.source(a);
   2.267 +                if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
   2.268 +                  reached.set(u, true);
   2.269 +                  addItem(u, bucket_num);
   2.270 +                  nqueue.push_back(u);
   2.271 +                }
   2.272 +              }
   2.273 +            }
   2.274 +            queue.swap(nqueue);
   2.275 +            ++bucket_num;
   2.276 +          }
   2.277 +          _sets.front().pop_front();
   2.278 +          --bucket_num;
   2.279 +          first_set = false;
   2.280 +        }
   2.281 +
   2.282 +        _bucket->set(_source, 0);
   2.283 +        _dormant[0] = true;
   2.284 +      }
   2.285 +      _source_set->set(_source, true);
   2.286 +
   2.287 +      Node target = _last[_sets.back().back()];
   2.288 +      {
   2.289 +        for (OutArcIt a(_graph, _source); a != INVALID; ++a) {
   2.290 +          if (_tolerance.positive((*_capacity)[a])) {
   2.291 +            Node u = _graph.target(a);
   2.292 +            _flow->set(a, (*_capacity)[a]);
   2.293 +            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
   2.294 +            if (!(*_active)[u] && u != _source) {
   2.295 +              activate(u);
   2.296 +            }
   2.297 +          }
   2.298 +        }
   2.299 +
   2.300 +        if ((*_active)[target]) {
   2.301 +          deactivate(target);
   2.302 +        }
   2.303 +
   2.304 +        _highest = _sets.back().begin();
   2.305 +        while (_highest != _sets.back().end() &&
   2.306 +               !(*_active)[_first[*_highest]]) {
   2.307 +          ++_highest;
   2.308 +        }
   2.309 +      }
   2.310 +
   2.311 +      while (true) {
   2.312 +        while (_highest != _sets.back().end()) {
   2.313 +          Node n = _first[*_highest];
   2.314 +          Value excess = (*_excess)[n];
   2.315 +          int next_bucket = _node_num;
   2.316 +
   2.317 +          int under_bucket;
   2.318 +          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
   2.319 +            under_bucket = -1;
   2.320 +          } else {
   2.321 +            under_bucket = *(++std::list<int>::iterator(_highest));
   2.322 +          }
   2.323 +
   2.324 +          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   2.325 +            Node v = _graph.target(a);
   2.326 +            if (_dormant[(*_bucket)[v]]) continue;
   2.327 +            Value rem = (*_capacity)[a] - (*_flow)[a];
   2.328 +            if (!_tolerance.positive(rem)) continue;
   2.329 +            if ((*_bucket)[v] == under_bucket) {
   2.330 +              if (!(*_active)[v] && v != target) {
   2.331 +                activate(v);
   2.332 +              }
   2.333 +              if (!_tolerance.less(rem, excess)) {
   2.334 +                _flow->set(a, (*_flow)[a] + excess);
   2.335 +                _excess->set(v, (*_excess)[v] + excess);
   2.336 +                excess = 0;
   2.337 +                goto no_more_push;
   2.338 +              } else {
   2.339 +                excess -= rem;
   2.340 +                _excess->set(v, (*_excess)[v] + rem);
   2.341 +                _flow->set(a, (*_capacity)[a]);
   2.342 +              }
   2.343 +            } else if (next_bucket > (*_bucket)[v]) {
   2.344 +              next_bucket = (*_bucket)[v];
   2.345 +            }
   2.346 +          }
   2.347 +
   2.348 +          for (InArcIt a(_graph, n); a != INVALID; ++a) {
   2.349 +            Node v = _graph.source(a);
   2.350 +            if (_dormant[(*_bucket)[v]]) continue;
   2.351 +            Value rem = (*_flow)[a];
   2.352 +            if (!_tolerance.positive(rem)) continue;
   2.353 +            if ((*_bucket)[v] == under_bucket) {
   2.354 +              if (!(*_active)[v] && v != target) {
   2.355 +                activate(v);
   2.356 +              }
   2.357 +              if (!_tolerance.less(rem, excess)) {
   2.358 +                _flow->set(a, (*_flow)[a] - excess);
   2.359 +                _excess->set(v, (*_excess)[v] + excess);
   2.360 +                excess = 0;
   2.361 +                goto no_more_push;
   2.362 +              } else {
   2.363 +                excess -= rem;
   2.364 +                _excess->set(v, (*_excess)[v] + rem);
   2.365 +                _flow->set(a, 0);
   2.366 +              }
   2.367 +            } else if (next_bucket > (*_bucket)[v]) {
   2.368 +              next_bucket = (*_bucket)[v];
   2.369 +            }
   2.370 +          }
   2.371 +
   2.372 +        no_more_push:
   2.373 +
   2.374 +          _excess->set(n, excess);
   2.375 +
   2.376 +          if (excess != 0) {
   2.377 +            if ((*_next)[n] == INVALID) {
   2.378 +              typename std::list<std::list<int> >::iterator new_set =
   2.379 +                _sets.insert(--_sets.end(), std::list<int>());
   2.380 +              new_set->splice(new_set->end(), _sets.back(),
   2.381 +                              _sets.back().begin(), ++_highest);
   2.382 +              for (std::list<int>::iterator it = new_set->begin();
   2.383 +                   it != new_set->end(); ++it) {
   2.384 +                _dormant[*it] = true;
   2.385 +              }
   2.386 +              while (_highest != _sets.back().end() &&
   2.387 +                     !(*_active)[_first[*_highest]]) {
   2.388 +                ++_highest;
   2.389 +              }
   2.390 +            } else if (next_bucket == _node_num) {
   2.391 +              _first[(*_bucket)[n]] = (*_next)[n];
   2.392 +              _prev->set((*_next)[n], INVALID);
   2.393 +
   2.394 +              std::list<std::list<int> >::iterator new_set =
   2.395 +                _sets.insert(--_sets.end(), std::list<int>());
   2.396 +
   2.397 +              new_set->push_front(bucket_num);
   2.398 +              _bucket->set(n, bucket_num);
   2.399 +              _first[bucket_num] = _last[bucket_num] = n;
   2.400 +              _next->set(n, INVALID);
   2.401 +              _prev->set(n, INVALID);
   2.402 +              _dormant[bucket_num] = true;
   2.403 +              ++bucket_num;
   2.404 +
   2.405 +              while (_highest != _sets.back().end() &&
   2.406 +                     !(*_active)[_first[*_highest]]) {
   2.407 +                ++_highest;
   2.408 +              }
   2.409 +            } else {
   2.410 +              _first[*_highest] = (*_next)[n];
   2.411 +              _prev->set((*_next)[n], INVALID);
   2.412 +
   2.413 +              while (next_bucket != *_highest) {
   2.414 +                --_highest;
   2.415 +              }
   2.416 +
   2.417 +              if (_highest == _sets.back().begin()) {
   2.418 +                _sets.back().push_front(bucket_num);
   2.419 +                _dormant[bucket_num] = false;
   2.420 +                _first[bucket_num] = _last[bucket_num] = INVALID;
   2.421 +                ++bucket_num;
   2.422 +              }
   2.423 +              --_highest;
   2.424 +
   2.425 +              _bucket->set(n, *_highest);
   2.426 +              _next->set(n, _first[*_highest]);
   2.427 +              if (_first[*_highest] != INVALID) {
   2.428 +                _prev->set(_first[*_highest], n);
   2.429 +              } else {
   2.430 +                _last[*_highest] = n;
   2.431 +              }
   2.432 +              _first[*_highest] = n;
   2.433 +            }
   2.434 +          } else {
   2.435 +
   2.436 +            deactivate(n);
   2.437 +            if (!(*_active)[_first[*_highest]]) {
   2.438 +              ++_highest;
   2.439 +              if (_highest != _sets.back().end() &&
   2.440 +                  !(*_active)[_first[*_highest]]) {
   2.441 +                _highest = _sets.back().end();
   2.442 +              }
   2.443 +            }
   2.444 +          }
   2.445 +        }
   2.446 +
   2.447 +        if ((*_excess)[target] < _min_cut) {
   2.448 +          _min_cut = (*_excess)[target];
   2.449 +          for (NodeIt i(_graph); i != INVALID; ++i) {
   2.450 +            _min_cut_map->set(i, true);
   2.451 +          }
   2.452 +          for (std::list<int>::iterator it = _sets.back().begin();
   2.453 +               it != _sets.back().end(); ++it) {
   2.454 +            Node n = _first[*it];
   2.455 +            while (n != INVALID) {
   2.456 +              _min_cut_map->set(n, false);
   2.457 +              n = (*_next)[n];
   2.458 +            }
   2.459 +          }
   2.460 +        }
   2.461 +
   2.462 +        {
   2.463 +          Node new_target;
   2.464 +          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
   2.465 +            if ((*_next)[target] == INVALID) {
   2.466 +              _last[(*_bucket)[target]] = (*_prev)[target];
   2.467 +              new_target = (*_prev)[target];
   2.468 +            } else {
   2.469 +              _prev->set((*_next)[target], (*_prev)[target]);
   2.470 +              new_target = (*_next)[target];
   2.471 +            }
   2.472 +            if ((*_prev)[target] == INVALID) {
   2.473 +              _first[(*_bucket)[target]] = (*_next)[target];
   2.474 +            } else {
   2.475 +              _next->set((*_prev)[target], (*_next)[target]);
   2.476 +            }
   2.477 +          } else {
   2.478 +            _sets.back().pop_back();
   2.479 +            if (_sets.back().empty()) {
   2.480 +              _sets.pop_back();
   2.481 +              if (_sets.empty())
   2.482 +                break;
   2.483 +              for (std::list<int>::iterator it = _sets.back().begin();
   2.484 +                   it != _sets.back().end(); ++it) {
   2.485 +                _dormant[*it] = false;
   2.486 +              }
   2.487 +            }
   2.488 +            new_target = _last[_sets.back().back()];
   2.489 +          }
   2.490 +
   2.491 +          _bucket->set(target, 0);
   2.492 +
   2.493 +          _source_set->set(target, true);
   2.494 +          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
   2.495 +            Value rem = (*_capacity)[a] - (*_flow)[a];
   2.496 +            if (!_tolerance.positive(rem)) continue;
   2.497 +            Node v = _graph.target(a);
   2.498 +            if (!(*_active)[v] && !(*_source_set)[v]) {
   2.499 +              activate(v);
   2.500 +            }
   2.501 +            _excess->set(v, (*_excess)[v] + rem);
   2.502 +            _flow->set(a, (*_capacity)[a]);
   2.503 +          }
   2.504 +
   2.505 +          for (InArcIt a(_graph, target); a != INVALID; ++a) {
   2.506 +            Value rem = (*_flow)[a];
   2.507 +            if (!_tolerance.positive(rem)) continue;
   2.508 +            Node v = _graph.source(a);
   2.509 +            if (!(*_active)[v] && !(*_source_set)[v]) {
   2.510 +              activate(v);
   2.511 +            }
   2.512 +            _excess->set(v, (*_excess)[v] + rem);
   2.513 +            _flow->set(a, 0);
   2.514 +          }
   2.515 +
   2.516 +          target = new_target;
   2.517 +          if ((*_active)[target]) {
   2.518 +            deactivate(target);
   2.519 +          }
   2.520 +
   2.521 +          _highest = _sets.back().begin();
   2.522 +          while (_highest != _sets.back().end() &&
   2.523 +                 !(*_active)[_first[*_highest]]) {
   2.524 +            ++_highest;
   2.525 +          }
   2.526 +        }
   2.527 +      }
   2.528 +    }
   2.529 +
   2.530 +    void findMinCutIn() {
   2.531 +
   2.532 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   2.533 +        _excess->set(n, 0);
   2.534 +      }
   2.535 +
   2.536 +      for (ArcIt a(_graph); a != INVALID; ++a) {
   2.537 +        _flow->set(a, 0);
   2.538 +      }
   2.539 +
   2.540 +      int bucket_num = 1;
   2.541 +
   2.542 +      {
   2.543 +        typename Digraph::template NodeMap<bool> reached(_graph, false);
   2.544 +
   2.545 +        reached.set(_source, true);
   2.546 +
   2.547 +        bool first_set = true;
   2.548 +
   2.549 +        for (NodeIt t(_graph); t != INVALID; ++t) {
   2.550 +          if (reached[t]) continue;
   2.551 +          _sets.push_front(std::list<int>());
   2.552 +          _sets.front().push_front(bucket_num);
   2.553 +          _dormant[bucket_num] = !first_set;
   2.554 +
   2.555 +          _bucket->set(t, bucket_num);
   2.556 +          _first[bucket_num] = _last[bucket_num] = t;
   2.557 +          _next->set(t, INVALID);
   2.558 +          _prev->set(t, INVALID);
   2.559 +
   2.560 +          ++bucket_num;
   2.561 +
   2.562 +          std::vector<Node> queue;
   2.563 +          queue.push_back(t);
   2.564 +          reached.set(t, true);
   2.565 +
   2.566 +          while (!queue.empty()) {
   2.567 +            _sets.front().push_front(bucket_num);
   2.568 +            _dormant[bucket_num] = !first_set;
   2.569 +            _first[bucket_num] = _last[bucket_num] = INVALID;
   2.570 +
   2.571 +            std::vector<Node> nqueue;
   2.572 +            for (int i = 0; i < int(queue.size()); ++i) {
   2.573 +              Node n = queue[i];
   2.574 +              for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   2.575 +                Node u = _graph.target(a);
   2.576 +                if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
   2.577 +                  reached.set(u, true);
   2.578 +                  addItem(u, bucket_num);
   2.579 +                  nqueue.push_back(u);
   2.580 +                }
   2.581 +              }
   2.582 +            }
   2.583 +            queue.swap(nqueue);
   2.584 +            ++bucket_num;
   2.585 +          }
   2.586 +          _sets.front().pop_front();
   2.587 +          --bucket_num;
   2.588 +          first_set = false;
   2.589 +        }
   2.590 +
   2.591 +        _bucket->set(_source, 0);
   2.592 +        _dormant[0] = true;
   2.593 +      }
   2.594 +      _source_set->set(_source, true);
   2.595 +
   2.596 +      Node target = _last[_sets.back().back()];
   2.597 +      {
   2.598 +        for (InArcIt a(_graph, _source); a != INVALID; ++a) {
   2.599 +          if (_tolerance.positive((*_capacity)[a])) {
   2.600 +            Node u = _graph.source(a);
   2.601 +            _flow->set(a, (*_capacity)[a]);
   2.602 +            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
   2.603 +            if (!(*_active)[u] && u != _source) {
   2.604 +              activate(u);
   2.605 +            }
   2.606 +          }
   2.607 +        }
   2.608 +        if ((*_active)[target]) {
   2.609 +          deactivate(target);
   2.610 +        }
   2.611 +
   2.612 +        _highest = _sets.back().begin();
   2.613 +        while (_highest != _sets.back().end() &&
   2.614 +               !(*_active)[_first[*_highest]]) {
   2.615 +          ++_highest;
   2.616 +        }
   2.617 +      }
   2.618 +
   2.619 +
   2.620 +      while (true) {
   2.621 +        while (_highest != _sets.back().end()) {
   2.622 +          Node n = _first[*_highest];
   2.623 +          Value excess = (*_excess)[n];
   2.624 +          int next_bucket = _node_num;
   2.625 +
   2.626 +          int under_bucket;
   2.627 +          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
   2.628 +            under_bucket = -1;
   2.629 +          } else {
   2.630 +            under_bucket = *(++std::list<int>::iterator(_highest));
   2.631 +          }
   2.632 +
   2.633 +          for (InArcIt a(_graph, n); a != INVALID; ++a) {
   2.634 +            Node v = _graph.source(a);
   2.635 +            if (_dormant[(*_bucket)[v]]) continue;
   2.636 +            Value rem = (*_capacity)[a] - (*_flow)[a];
   2.637 +            if (!_tolerance.positive(rem)) continue;
   2.638 +            if ((*_bucket)[v] == under_bucket) {
   2.639 +              if (!(*_active)[v] && v != target) {
   2.640 +                activate(v);
   2.641 +              }
   2.642 +              if (!_tolerance.less(rem, excess)) {
   2.643 +                _flow->set(a, (*_flow)[a] + excess);
   2.644 +                _excess->set(v, (*_excess)[v] + excess);
   2.645 +                excess = 0;
   2.646 +                goto no_more_push;
   2.647 +              } else {
   2.648 +                excess -= rem;
   2.649 +                _excess->set(v, (*_excess)[v] + rem);
   2.650 +                _flow->set(a, (*_capacity)[a]);
   2.651 +              }
   2.652 +            } else if (next_bucket > (*_bucket)[v]) {
   2.653 +              next_bucket = (*_bucket)[v];
   2.654 +            }
   2.655 +          }
   2.656 +
   2.657 +          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   2.658 +            Node v = _graph.target(a);
   2.659 +            if (_dormant[(*_bucket)[v]]) continue;
   2.660 +            Value rem = (*_flow)[a];
   2.661 +            if (!_tolerance.positive(rem)) continue;
   2.662 +            if ((*_bucket)[v] == under_bucket) {
   2.663 +              if (!(*_active)[v] && v != target) {
   2.664 +                activate(v);
   2.665 +              }
   2.666 +              if (!_tolerance.less(rem, excess)) {
   2.667 +                _flow->set(a, (*_flow)[a] - excess);
   2.668 +                _excess->set(v, (*_excess)[v] + excess);
   2.669 +                excess = 0;
   2.670 +                goto no_more_push;
   2.671 +              } else {
   2.672 +                excess -= rem;
   2.673 +                _excess->set(v, (*_excess)[v] + rem);
   2.674 +                _flow->set(a, 0);
   2.675 +              }
   2.676 +            } else if (next_bucket > (*_bucket)[v]) {
   2.677 +              next_bucket = (*_bucket)[v];
   2.678 +            }
   2.679 +          }
   2.680 +
   2.681 +        no_more_push:
   2.682 +
   2.683 +          _excess->set(n, excess);
   2.684 +
   2.685 +          if (excess != 0) {
   2.686 +            if ((*_next)[n] == INVALID) {
   2.687 +              typename std::list<std::list<int> >::iterator new_set =
   2.688 +                _sets.insert(--_sets.end(), std::list<int>());
   2.689 +              new_set->splice(new_set->end(), _sets.back(),
   2.690 +                              _sets.back().begin(), ++_highest);
   2.691 +              for (std::list<int>::iterator it = new_set->begin();
   2.692 +                   it != new_set->end(); ++it) {
   2.693 +                _dormant[*it] = true;
   2.694 +              }
   2.695 +              while (_highest != _sets.back().end() &&
   2.696 +                     !(*_active)[_first[*_highest]]) {
   2.697 +                ++_highest;
   2.698 +              }
   2.699 +            } else if (next_bucket == _node_num) {
   2.700 +              _first[(*_bucket)[n]] = (*_next)[n];
   2.701 +              _prev->set((*_next)[n], INVALID);
   2.702 +
   2.703 +              std::list<std::list<int> >::iterator new_set =
   2.704 +                _sets.insert(--_sets.end(), std::list<int>());
   2.705 +
   2.706 +              new_set->push_front(bucket_num);
   2.707 +              _bucket->set(n, bucket_num);
   2.708 +              _first[bucket_num] = _last[bucket_num] = n;
   2.709 +              _next->set(n, INVALID);
   2.710 +              _prev->set(n, INVALID);
   2.711 +              _dormant[bucket_num] = true;
   2.712 +              ++bucket_num;
   2.713 +
   2.714 +              while (_highest != _sets.back().end() &&
   2.715 +                     !(*_active)[_first[*_highest]]) {
   2.716 +                ++_highest;
   2.717 +              }
   2.718 +            } else {
   2.719 +              _first[*_highest] = (*_next)[n];
   2.720 +              _prev->set((*_next)[n], INVALID);
   2.721 +
   2.722 +              while (next_bucket != *_highest) {
   2.723 +                --_highest;
   2.724 +              }
   2.725 +              if (_highest == _sets.back().begin()) {
   2.726 +                _sets.back().push_front(bucket_num);
   2.727 +                _dormant[bucket_num] = false;
   2.728 +                _first[bucket_num] = _last[bucket_num] = INVALID;
   2.729 +                ++bucket_num;
   2.730 +              }
   2.731 +              --_highest;
   2.732 +
   2.733 +              _bucket->set(n, *_highest);
   2.734 +              _next->set(n, _first[*_highest]);
   2.735 +              if (_first[*_highest] != INVALID) {
   2.736 +                _prev->set(_first[*_highest], n);
   2.737 +              } else {
   2.738 +                _last[*_highest] = n;
   2.739 +              }
   2.740 +              _first[*_highest] = n;
   2.741 +            }
   2.742 +          } else {
   2.743 +
   2.744 +            deactivate(n);
   2.745 +            if (!(*_active)[_first[*_highest]]) {
   2.746 +              ++_highest;
   2.747 +              if (_highest != _sets.back().end() &&
   2.748 +                  !(*_active)[_first[*_highest]]) {
   2.749 +                _highest = _sets.back().end();
   2.750 +              }
   2.751 +            }
   2.752 +          }
   2.753 +        }
   2.754 +
   2.755 +        if ((*_excess)[target] < _min_cut) {
   2.756 +          _min_cut = (*_excess)[target];
   2.757 +          for (NodeIt i(_graph); i != INVALID; ++i) {
   2.758 +            _min_cut_map->set(i, false);
   2.759 +          }
   2.760 +          for (std::list<int>::iterator it = _sets.back().begin();
   2.761 +               it != _sets.back().end(); ++it) {
   2.762 +            Node n = _first[*it];
   2.763 +            while (n != INVALID) {
   2.764 +              _min_cut_map->set(n, true);
   2.765 +              n = (*_next)[n];
   2.766 +            }
   2.767 +          }
   2.768 +        }
   2.769 +
   2.770 +        {
   2.771 +          Node new_target;
   2.772 +          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
   2.773 +            if ((*_next)[target] == INVALID) {
   2.774 +              _last[(*_bucket)[target]] = (*_prev)[target];
   2.775 +              new_target = (*_prev)[target];
   2.776 +            } else {
   2.777 +              _prev->set((*_next)[target], (*_prev)[target]);
   2.778 +              new_target = (*_next)[target];
   2.779 +            }
   2.780 +            if ((*_prev)[target] == INVALID) {
   2.781 +              _first[(*_bucket)[target]] = (*_next)[target];
   2.782 +            } else {
   2.783 +              _next->set((*_prev)[target], (*_next)[target]);
   2.784 +            }
   2.785 +          } else {
   2.786 +            _sets.back().pop_back();
   2.787 +            if (_sets.back().empty()) {
   2.788 +              _sets.pop_back();
   2.789 +              if (_sets.empty())
   2.790 +                break;
   2.791 +              for (std::list<int>::iterator it = _sets.back().begin();
   2.792 +                   it != _sets.back().end(); ++it) {
   2.793 +                _dormant[*it] = false;
   2.794 +              }
   2.795 +            }
   2.796 +            new_target = _last[_sets.back().back()];
   2.797 +          }
   2.798 +
   2.799 +          _bucket->set(target, 0);
   2.800 +
   2.801 +          _source_set->set(target, true);
   2.802 +          for (InArcIt a(_graph, target); a != INVALID; ++a) {
   2.803 +            Value rem = (*_capacity)[a] - (*_flow)[a];
   2.804 +            if (!_tolerance.positive(rem)) continue;
   2.805 +            Node v = _graph.source(a);
   2.806 +            if (!(*_active)[v] && !(*_source_set)[v]) {
   2.807 +              activate(v);
   2.808 +            }
   2.809 +            _excess->set(v, (*_excess)[v] + rem);
   2.810 +            _flow->set(a, (*_capacity)[a]);
   2.811 +          }
   2.812 +
   2.813 +          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
   2.814 +            Value rem = (*_flow)[a];
   2.815 +            if (!_tolerance.positive(rem)) continue;
   2.816 +            Node v = _graph.target(a);
   2.817 +            if (!(*_active)[v] && !(*_source_set)[v]) {
   2.818 +              activate(v);
   2.819 +            }
   2.820 +            _excess->set(v, (*_excess)[v] + rem);
   2.821 +            _flow->set(a, 0);
   2.822 +          }
   2.823 +
   2.824 +          target = new_target;
   2.825 +          if ((*_active)[target]) {
   2.826 +            deactivate(target);
   2.827 +          }
   2.828 +
   2.829 +          _highest = _sets.back().begin();
   2.830 +          while (_highest != _sets.back().end() &&
   2.831 +                 !(*_active)[_first[*_highest]]) {
   2.832 +            ++_highest;
   2.833 +          }
   2.834 +        }
   2.835 +      }
   2.836 +    }
   2.837 +
   2.838 +  public:
   2.839 +
   2.840 +    /// \name Execution control
   2.841 +    /// The simplest way to execute the algorithm is to use
   2.842 +    /// one of the member functions called \c run(...).
   2.843 +    /// \n
   2.844 +    /// If you need more control on the execution,
   2.845 +    /// first you must call \ref init(), then the \ref calculateIn() or
   2.846 +    /// \ref calculateIn() functions.
   2.847 +
   2.848 +    /// @{
   2.849 +
   2.850 +    /// \brief Initializes the internal data structures.
   2.851 +    ///
   2.852 +    /// Initializes the internal data structures. It creates
   2.853 +    /// the maps, residual graph adaptors and some bucket structures
   2.854 +    /// for the algorithm.
   2.855 +    void init() {
   2.856 +      init(NodeIt(_graph));
   2.857 +    }
   2.858 +
   2.859 +    /// \brief Initializes the internal data structures.
   2.860 +    ///
   2.861 +    /// Initializes the internal data structures. It creates
   2.862 +    /// the maps, residual graph adaptor and some bucket structures
   2.863 +    /// for the algorithm. Node \c source  is used as the push-relabel
   2.864 +    /// algorithm's source.
   2.865 +    void init(const Node& source) {
   2.866 +      _source = source;
   2.867 +
   2.868 +      _node_num = countNodes(_graph);
   2.869 +
   2.870 +      _first.resize(_node_num + 1);
   2.871 +      _last.resize(_node_num + 1);
   2.872 +
   2.873 +      _dormant.resize(_node_num + 1);
   2.874 +
   2.875 +      if (!_flow) {
   2.876 +        _flow = new FlowMap(_graph);
   2.877 +      }
   2.878 +      if (!_next) {
   2.879 +        _next = new typename Digraph::template NodeMap<Node>(_graph);
   2.880 +      }
   2.881 +      if (!_prev) {
   2.882 +        _prev = new typename Digraph::template NodeMap<Node>(_graph);
   2.883 +      }
   2.884 +      if (!_active) {
   2.885 +        _active = new typename Digraph::template NodeMap<bool>(_graph);
   2.886 +      }
   2.887 +      if (!_bucket) {
   2.888 +        _bucket = new typename Digraph::template NodeMap<int>(_graph);
   2.889 +      }
   2.890 +      if (!_excess) {
   2.891 +        _excess = new ExcessMap(_graph);
   2.892 +      }
   2.893 +      if (!_source_set) {
   2.894 +        _source_set = new SourceSetMap(_graph);
   2.895 +      }
   2.896 +      if (!_min_cut_map) {
   2.897 +        _min_cut_map = new MinCutMap(_graph);
   2.898 +      }
   2.899 +
   2.900 +      _min_cut = std::numeric_limits<Value>::max();
   2.901 +    }
   2.902 +
   2.903 +
   2.904 +    /// \brief Calculates a minimum cut with \f$ source \f$ on the
   2.905 +    /// source-side.
   2.906 +    ///
   2.907 +    /// Calculates a minimum cut with \f$ source \f$ on the
   2.908 +    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
   2.909 +    /// \in X \f$ and minimal out-degree).
   2.910 +    void calculateOut() {
   2.911 +      findMinCutOut();
   2.912 +    }
   2.913 +
   2.914 +    /// \brief Calculates a minimum cut with \f$ source \f$ on the
   2.915 +    /// target-side.
   2.916 +    ///
   2.917 +    /// Calculates a minimum cut with \f$ source \f$ on the
   2.918 +    /// target-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
   2.919 +    /// \in X \f$ and minimal out-degree).
   2.920 +    void calculateIn() {
   2.921 +      findMinCutIn();
   2.922 +    }
   2.923 +
   2.924 +
   2.925 +    /// \brief Runs the algorithm.
   2.926 +    ///
   2.927 +    /// Runs the algorithm. It finds nodes \c source and \c target
   2.928 +    /// arbitrarily and then calls \ref init(), \ref calculateOut()
   2.929 +    /// and \ref calculateIn().
   2.930 +    void run() {
   2.931 +      init();
   2.932 +      calculateOut();
   2.933 +      calculateIn();
   2.934 +    }
   2.935 +
   2.936 +    /// \brief Runs the algorithm.
   2.937 +    ///
   2.938 +    /// Runs the algorithm. It uses the given \c source node, finds a
   2.939 +    /// proper \c target and then calls the \ref init(), \ref
   2.940 +    /// calculateOut() and \ref calculateIn().
   2.941 +    void run(const Node& s) {
   2.942 +      init(s);
   2.943 +      calculateOut();
   2.944 +      calculateIn();
   2.945 +    }
   2.946 +
   2.947 +    /// @}
   2.948 +
   2.949 +    /// \name Query Functions
   2.950 +    /// The result of the %HaoOrlin algorithm
   2.951 +    /// can be obtained using these functions.
   2.952 +    /// \n
   2.953 +    /// Before using these functions, either \ref run(), \ref
   2.954 +    /// calculateOut() or \ref calculateIn() must be called.
   2.955 +
   2.956 +    /// @{
   2.957 +
   2.958 +    /// \brief Returns the value of the minimum value cut.
   2.959 +    ///
   2.960 +    /// Returns the value of the minimum value cut.
   2.961 +    Value minCutValue() const {
   2.962 +      return _min_cut;
   2.963 +    }
   2.964 +
   2.965 +
   2.966 +    /// \brief Returns a minimum cut.
   2.967 +    ///
   2.968 +    /// Sets \c nodeMap to the characteristic vector of a minimum
   2.969 +    /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
   2.970 +    /// with minimal out-degree (i.e. \c nodeMap will be true exactly
   2.971 +    /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
   2.972 +    /// bool-valued node-map.
   2.973 +    template <typename NodeMap>
   2.974 +    Value minCutMap(NodeMap& nodeMap) const {
   2.975 +      for (NodeIt it(_graph); it != INVALID; ++it) {
   2.976 +        nodeMap.set(it, (*_min_cut_map)[it]);
   2.977 +      }
   2.978 +      return _min_cut;
   2.979 +    }
   2.980 +
   2.981 +    /// @}
   2.982 +
   2.983 +  }; //class HaoOrlin
   2.984 +
   2.985 +
   2.986 +} //namespace lemon
   2.987 +
   2.988 +#endif //LEMON_HAO_ORLIN_H