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