lemon/hao_orlin.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 17 Apr 2009 18:04:36 +0200
changeset 609 e6927fe719e6
parent 412 7030149efed2
child 559 c5fd2d996909
permissions -rw-r--r--
Support >= and <= constraints in NetworkSimplex (#219, #234)

By default the same inequality constraints are supported as by
Circulation (the GEQ form), but the LEQ form can also be selected
using the problemType() function.

The documentation of the min. cost flow module is reworked and
extended with important notes and explanations about the different
variants of the problem and about the dual solution and optimality
conditions.
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
 *
alpar@440
     5
 * Copyright (C) 2003-2009
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>());
alpar@440
   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>());
alpar@440
   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