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