1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/hao_orlin.h Thu Nov 05 15:50:01 2009 +0100
1.3 @@ -0,0 +1,988 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2009
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_HAO_ORLIN_H
1.23 +#define LEMON_HAO_ORLIN_H
1.24 +
1.25 +#include <vector>
1.26 +#include <list>
1.27 +#include <limits>
1.28 +
1.29 +#include <lemon/maps.h>
1.30 +#include <lemon/core.h>
1.31 +#include <lemon/tolerance.h>
1.32 +
1.33 +/// \file
1.34 +/// \ingroup min_cut
1.35 +/// \brief Implementation of the Hao-Orlin algorithm.
1.36 +///
1.37 +/// Implementation of the Hao-Orlin algorithm for finding a minimum cut
1.38 +/// in a digraph.
1.39 +
1.40 +namespace lemon {
1.41 +
1.42 + /// \ingroup min_cut
1.43 + ///
1.44 + /// \brief Hao-Orlin algorithm for finding a minimum cut in a digraph.
1.45 + ///
1.46 + /// This class implements the Hao-Orlin algorithm for finding a minimum
1.47 + /// value cut in a directed graph \f$D=(V,A)\f$.
1.48 + /// It takes a fixed node \f$ source \in V \f$ and
1.49 + /// consists of two phases: in the first phase it determines a
1.50 + /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
1.51 + /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal outgoing
1.52 + /// capacity) and in the second phase it determines a minimum cut
1.53 + /// with \f$ source \f$ on the sink-side (i.e. a set
1.54 + /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal outgoing
1.55 + /// capacity). Obviously, the smaller of these two cuts will be a
1.56 + /// minimum cut of \f$ D \f$. The algorithm is a modified
1.57 + /// preflow push-relabel algorithm. Our implementation calculates
1.58 + /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
1.59 + /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
1.60 + /// purpose of such algorithm is e.g. testing network reliability.
1.61 + ///
1.62 + /// For an undirected graph you can run just the first phase of the
1.63 + /// algorithm or you can use the algorithm of Nagamochi and Ibaraki,
1.64 + /// which solves the undirected problem in \f$ O(nm + n^2 \log n) \f$
1.65 + /// time. It is implemented in the NagamochiIbaraki algorithm class.
1.66 + ///
1.67 + /// \tparam GR The type of the digraph the algorithm runs on.
1.68 + /// \tparam CAP The type of the arc map containing the capacities,
1.69 + /// which can be any numreric type. The default map type is
1.70 + /// \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
1.71 + /// \tparam TOL Tolerance class for handling inexact computations. The
1.72 + /// default tolerance type is \ref Tolerance "Tolerance<CAP::Value>".
1.73 +#ifdef DOXYGEN
1.74 + template <typename GR, typename CAP, typename TOL>
1.75 +#else
1.76 + template <typename GR,
1.77 + typename CAP = typename GR::template ArcMap<int>,
1.78 + typename TOL = Tolerance<typename CAP::Value> >
1.79 +#endif
1.80 + class HaoOrlin {
1.81 + public:
1.82 +
1.83 + /// The digraph type of the algorithm
1.84 + typedef GR Digraph;
1.85 + /// The capacity map type of the algorithm
1.86 + typedef CAP CapacityMap;
1.87 + /// The tolerance type of the algorithm
1.88 + typedef TOL Tolerance;
1.89 +
1.90 + private:
1.91 +
1.92 + typedef typename CapacityMap::Value Value;
1.93 +
1.94 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.95 +
1.96 + const Digraph& _graph;
1.97 + const CapacityMap* _capacity;
1.98 +
1.99 + typedef typename Digraph::template ArcMap<Value> FlowMap;
1.100 + FlowMap* _flow;
1.101 +
1.102 + Node _source;
1.103 +
1.104 + int _node_num;
1.105 +
1.106 + // Bucketing structure
1.107 + std::vector<Node> _first, _last;
1.108 + typename Digraph::template NodeMap<Node>* _next;
1.109 + typename Digraph::template NodeMap<Node>* _prev;
1.110 + typename Digraph::template NodeMap<bool>* _active;
1.111 + typename Digraph::template NodeMap<int>* _bucket;
1.112 +
1.113 + std::vector<bool> _dormant;
1.114 +
1.115 + std::list<std::list<int> > _sets;
1.116 + std::list<int>::iterator _highest;
1.117 +
1.118 + typedef typename Digraph::template NodeMap<Value> ExcessMap;
1.119 + ExcessMap* _excess;
1.120 +
1.121 + typedef typename Digraph::template NodeMap<bool> SourceSetMap;
1.122 + SourceSetMap* _source_set;
1.123 +
1.124 + Value _min_cut;
1.125 +
1.126 + typedef typename Digraph::template NodeMap<bool> MinCutMap;
1.127 + MinCutMap* _min_cut_map;
1.128 +
1.129 + Tolerance _tolerance;
1.130 +
1.131 + public:
1.132 +
1.133 + /// \brief Constructor
1.134 + ///
1.135 + /// Constructor of the algorithm class.
1.136 + HaoOrlin(const Digraph& graph, const CapacityMap& capacity,
1.137 + const Tolerance& tolerance = Tolerance()) :
1.138 + _graph(graph), _capacity(&capacity), _flow(0), _source(),
1.139 + _node_num(), _first(), _last(), _next(0), _prev(0),
1.140 + _active(0), _bucket(0), _dormant(), _sets(), _highest(),
1.141 + _excess(0), _source_set(0), _min_cut(), _min_cut_map(0),
1.142 + _tolerance(tolerance) {}
1.143 +
1.144 + ~HaoOrlin() {
1.145 + if (_min_cut_map) {
1.146 + delete _min_cut_map;
1.147 + }
1.148 + if (_source_set) {
1.149 + delete _source_set;
1.150 + }
1.151 + if (_excess) {
1.152 + delete _excess;
1.153 + }
1.154 + if (_next) {
1.155 + delete _next;
1.156 + }
1.157 + if (_prev) {
1.158 + delete _prev;
1.159 + }
1.160 + if (_active) {
1.161 + delete _active;
1.162 + }
1.163 + if (_bucket) {
1.164 + delete _bucket;
1.165 + }
1.166 + if (_flow) {
1.167 + delete _flow;
1.168 + }
1.169 + }
1.170 +
1.171 + private:
1.172 +
1.173 + void activate(const Node& i) {
1.174 + (*_active)[i] = true;
1.175 +
1.176 + int bucket = (*_bucket)[i];
1.177 +
1.178 + if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;
1.179 + //unlace
1.180 + (*_next)[(*_prev)[i]] = (*_next)[i];
1.181 + if ((*_next)[i] != INVALID) {
1.182 + (*_prev)[(*_next)[i]] = (*_prev)[i];
1.183 + } else {
1.184 + _last[bucket] = (*_prev)[i];
1.185 + }
1.186 + //lace
1.187 + (*_next)[i] = _first[bucket];
1.188 + (*_prev)[_first[bucket]] = i;
1.189 + (*_prev)[i] = INVALID;
1.190 + _first[bucket] = i;
1.191 + }
1.192 +
1.193 + void deactivate(const Node& i) {
1.194 + (*_active)[i] = false;
1.195 + int bucket = (*_bucket)[i];
1.196 +
1.197 + if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
1.198 +
1.199 + //unlace
1.200 + (*_prev)[(*_next)[i]] = (*_prev)[i];
1.201 + if ((*_prev)[i] != INVALID) {
1.202 + (*_next)[(*_prev)[i]] = (*_next)[i];
1.203 + } else {
1.204 + _first[bucket] = (*_next)[i];
1.205 + }
1.206 + //lace
1.207 + (*_prev)[i] = _last[bucket];
1.208 + (*_next)[_last[bucket]] = i;
1.209 + (*_next)[i] = INVALID;
1.210 + _last[bucket] = i;
1.211 + }
1.212 +
1.213 + void addItem(const Node& i, int bucket) {
1.214 + (*_bucket)[i] = bucket;
1.215 + if (_last[bucket] != INVALID) {
1.216 + (*_prev)[i] = _last[bucket];
1.217 + (*_next)[_last[bucket]] = i;
1.218 + (*_next)[i] = INVALID;
1.219 + _last[bucket] = i;
1.220 + } else {
1.221 + (*_prev)[i] = INVALID;
1.222 + _first[bucket] = i;
1.223 + (*_next)[i] = INVALID;
1.224 + _last[bucket] = i;
1.225 + }
1.226 + }
1.227 +
1.228 + void findMinCutOut() {
1.229 +
1.230 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.231 + (*_excess)[n] = 0;
1.232 + (*_source_set)[n] = false;
1.233 + }
1.234 +
1.235 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.236 + (*_flow)[a] = 0;
1.237 + }
1.238 +
1.239 + int bucket_num = 0;
1.240 + std::vector<Node> queue(_node_num);
1.241 + int qfirst = 0, qlast = 0, qsep = 0;
1.242 +
1.243 + {
1.244 + typename Digraph::template NodeMap<bool> reached(_graph, false);
1.245 +
1.246 + reached[_source] = true;
1.247 + bool first_set = true;
1.248 +
1.249 + for (NodeIt t(_graph); t != INVALID; ++t) {
1.250 + if (reached[t]) continue;
1.251 + _sets.push_front(std::list<int>());
1.252 +
1.253 + queue[qlast++] = t;
1.254 + reached[t] = true;
1.255 +
1.256 + while (qfirst != qlast) {
1.257 + if (qsep == qfirst) {
1.258 + ++bucket_num;
1.259 + _sets.front().push_front(bucket_num);
1.260 + _dormant[bucket_num] = !first_set;
1.261 + _first[bucket_num] = _last[bucket_num] = INVALID;
1.262 + qsep = qlast;
1.263 + }
1.264 +
1.265 + Node n = queue[qfirst++];
1.266 + addItem(n, bucket_num);
1.267 +
1.268 + for (InArcIt a(_graph, n); a != INVALID; ++a) {
1.269 + Node u = _graph.source(a);
1.270 + if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
1.271 + reached[u] = true;
1.272 + queue[qlast++] = u;
1.273 + }
1.274 + }
1.275 + }
1.276 + first_set = false;
1.277 + }
1.278 +
1.279 + ++bucket_num;
1.280 + (*_bucket)[_source] = 0;
1.281 + _dormant[0] = true;
1.282 + }
1.283 + (*_source_set)[_source] = true;
1.284 +
1.285 + Node target = _last[_sets.back().back()];
1.286 + {
1.287 + for (OutArcIt a(_graph, _source); a != INVALID; ++a) {
1.288 + if (_tolerance.positive((*_capacity)[a])) {
1.289 + Node u = _graph.target(a);
1.290 + (*_flow)[a] = (*_capacity)[a];
1.291 + (*_excess)[u] += (*_capacity)[a];
1.292 + if (!(*_active)[u] && u != _source) {
1.293 + activate(u);
1.294 + }
1.295 + }
1.296 + }
1.297 +
1.298 + if ((*_active)[target]) {
1.299 + deactivate(target);
1.300 + }
1.301 +
1.302 + _highest = _sets.back().begin();
1.303 + while (_highest != _sets.back().end() &&
1.304 + !(*_active)[_first[*_highest]]) {
1.305 + ++_highest;
1.306 + }
1.307 + }
1.308 +
1.309 + while (true) {
1.310 + while (_highest != _sets.back().end()) {
1.311 + Node n = _first[*_highest];
1.312 + Value excess = (*_excess)[n];
1.313 + int next_bucket = _node_num;
1.314 +
1.315 + int under_bucket;
1.316 + if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
1.317 + under_bucket = -1;
1.318 + } else {
1.319 + under_bucket = *(++std::list<int>::iterator(_highest));
1.320 + }
1.321 +
1.322 + for (OutArcIt a(_graph, n); a != INVALID; ++a) {
1.323 + Node v = _graph.target(a);
1.324 + if (_dormant[(*_bucket)[v]]) continue;
1.325 + Value rem = (*_capacity)[a] - (*_flow)[a];
1.326 + if (!_tolerance.positive(rem)) continue;
1.327 + if ((*_bucket)[v] == under_bucket) {
1.328 + if (!(*_active)[v] && v != target) {
1.329 + activate(v);
1.330 + }
1.331 + if (!_tolerance.less(rem, excess)) {
1.332 + (*_flow)[a] += excess;
1.333 + (*_excess)[v] += excess;
1.334 + excess = 0;
1.335 + goto no_more_push;
1.336 + } else {
1.337 + excess -= rem;
1.338 + (*_excess)[v] += rem;
1.339 + (*_flow)[a] = (*_capacity)[a];
1.340 + }
1.341 + } else if (next_bucket > (*_bucket)[v]) {
1.342 + next_bucket = (*_bucket)[v];
1.343 + }
1.344 + }
1.345 +
1.346 + for (InArcIt a(_graph, n); a != INVALID; ++a) {
1.347 + Node v = _graph.source(a);
1.348 + if (_dormant[(*_bucket)[v]]) continue;
1.349 + Value rem = (*_flow)[a];
1.350 + if (!_tolerance.positive(rem)) continue;
1.351 + if ((*_bucket)[v] == under_bucket) {
1.352 + if (!(*_active)[v] && v != target) {
1.353 + activate(v);
1.354 + }
1.355 + if (!_tolerance.less(rem, excess)) {
1.356 + (*_flow)[a] -= excess;
1.357 + (*_excess)[v] += excess;
1.358 + excess = 0;
1.359 + goto no_more_push;
1.360 + } else {
1.361 + excess -= rem;
1.362 + (*_excess)[v] += rem;
1.363 + (*_flow)[a] = 0;
1.364 + }
1.365 + } else if (next_bucket > (*_bucket)[v]) {
1.366 + next_bucket = (*_bucket)[v];
1.367 + }
1.368 + }
1.369 +
1.370 + no_more_push:
1.371 +
1.372 + (*_excess)[n] = excess;
1.373 +
1.374 + if (excess != 0) {
1.375 + if ((*_next)[n] == INVALID) {
1.376 + typename std::list<std::list<int> >::iterator new_set =
1.377 + _sets.insert(--_sets.end(), std::list<int>());
1.378 + new_set->splice(new_set->end(), _sets.back(),
1.379 + _sets.back().begin(), ++_highest);
1.380 + for (std::list<int>::iterator it = new_set->begin();
1.381 + it != new_set->end(); ++it) {
1.382 + _dormant[*it] = true;
1.383 + }
1.384 + while (_highest != _sets.back().end() &&
1.385 + !(*_active)[_first[*_highest]]) {
1.386 + ++_highest;
1.387 + }
1.388 + } else if (next_bucket == _node_num) {
1.389 + _first[(*_bucket)[n]] = (*_next)[n];
1.390 + (*_prev)[(*_next)[n]] = INVALID;
1.391 +
1.392 + std::list<std::list<int> >::iterator new_set =
1.393 + _sets.insert(--_sets.end(), std::list<int>());
1.394 +
1.395 + new_set->push_front(bucket_num);
1.396 + (*_bucket)[n] = bucket_num;
1.397 + _first[bucket_num] = _last[bucket_num] = n;
1.398 + (*_next)[n] = INVALID;
1.399 + (*_prev)[n] = INVALID;
1.400 + _dormant[bucket_num] = true;
1.401 + ++bucket_num;
1.402 +
1.403 + while (_highest != _sets.back().end() &&
1.404 + !(*_active)[_first[*_highest]]) {
1.405 + ++_highest;
1.406 + }
1.407 + } else {
1.408 + _first[*_highest] = (*_next)[n];
1.409 + (*_prev)[(*_next)[n]] = INVALID;
1.410 +
1.411 + while (next_bucket != *_highest) {
1.412 + --_highest;
1.413 + }
1.414 +
1.415 + if (_highest == _sets.back().begin()) {
1.416 + _sets.back().push_front(bucket_num);
1.417 + _dormant[bucket_num] = false;
1.418 + _first[bucket_num] = _last[bucket_num] = INVALID;
1.419 + ++bucket_num;
1.420 + }
1.421 + --_highest;
1.422 +
1.423 + (*_bucket)[n] = *_highest;
1.424 + (*_next)[n] = _first[*_highest];
1.425 + if (_first[*_highest] != INVALID) {
1.426 + (*_prev)[_first[*_highest]] = n;
1.427 + } else {
1.428 + _last[*_highest] = n;
1.429 + }
1.430 + _first[*_highest] = n;
1.431 + }
1.432 + } else {
1.433 +
1.434 + deactivate(n);
1.435 + if (!(*_active)[_first[*_highest]]) {
1.436 + ++_highest;
1.437 + if (_highest != _sets.back().end() &&
1.438 + !(*_active)[_first[*_highest]]) {
1.439 + _highest = _sets.back().end();
1.440 + }
1.441 + }
1.442 + }
1.443 + }
1.444 +
1.445 + if ((*_excess)[target] < _min_cut) {
1.446 + _min_cut = (*_excess)[target];
1.447 + for (NodeIt i(_graph); i != INVALID; ++i) {
1.448 + (*_min_cut_map)[i] = true;
1.449 + }
1.450 + for (std::list<int>::iterator it = _sets.back().begin();
1.451 + it != _sets.back().end(); ++it) {
1.452 + Node n = _first[*it];
1.453 + while (n != INVALID) {
1.454 + (*_min_cut_map)[n] = false;
1.455 + n = (*_next)[n];
1.456 + }
1.457 + }
1.458 + }
1.459 +
1.460 + {
1.461 + Node new_target;
1.462 + if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
1.463 + if ((*_next)[target] == INVALID) {
1.464 + _last[(*_bucket)[target]] = (*_prev)[target];
1.465 + new_target = (*_prev)[target];
1.466 + } else {
1.467 + (*_prev)[(*_next)[target]] = (*_prev)[target];
1.468 + new_target = (*_next)[target];
1.469 + }
1.470 + if ((*_prev)[target] == INVALID) {
1.471 + _first[(*_bucket)[target]] = (*_next)[target];
1.472 + } else {
1.473 + (*_next)[(*_prev)[target]] = (*_next)[target];
1.474 + }
1.475 + } else {
1.476 + _sets.back().pop_back();
1.477 + if (_sets.back().empty()) {
1.478 + _sets.pop_back();
1.479 + if (_sets.empty())
1.480 + break;
1.481 + for (std::list<int>::iterator it = _sets.back().begin();
1.482 + it != _sets.back().end(); ++it) {
1.483 + _dormant[*it] = false;
1.484 + }
1.485 + }
1.486 + new_target = _last[_sets.back().back()];
1.487 + }
1.488 +
1.489 + (*_bucket)[target] = 0;
1.490 +
1.491 + (*_source_set)[target] = true;
1.492 + for (OutArcIt a(_graph, target); a != INVALID; ++a) {
1.493 + Value rem = (*_capacity)[a] - (*_flow)[a];
1.494 + if (!_tolerance.positive(rem)) continue;
1.495 + Node v = _graph.target(a);
1.496 + if (!(*_active)[v] && !(*_source_set)[v]) {
1.497 + activate(v);
1.498 + }
1.499 + (*_excess)[v] += rem;
1.500 + (*_flow)[a] = (*_capacity)[a];
1.501 + }
1.502 +
1.503 + for (InArcIt a(_graph, target); a != INVALID; ++a) {
1.504 + Value rem = (*_flow)[a];
1.505 + if (!_tolerance.positive(rem)) continue;
1.506 + Node v = _graph.source(a);
1.507 + if (!(*_active)[v] && !(*_source_set)[v]) {
1.508 + activate(v);
1.509 + }
1.510 + (*_excess)[v] += rem;
1.511 + (*_flow)[a] = 0;
1.512 + }
1.513 +
1.514 + target = new_target;
1.515 + if ((*_active)[target]) {
1.516 + deactivate(target);
1.517 + }
1.518 +
1.519 + _highest = _sets.back().begin();
1.520 + while (_highest != _sets.back().end() &&
1.521 + !(*_active)[_first[*_highest]]) {
1.522 + ++_highest;
1.523 + }
1.524 + }
1.525 + }
1.526 + }
1.527 +
1.528 + void findMinCutIn() {
1.529 +
1.530 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.531 + (*_excess)[n] = 0;
1.532 + (*_source_set)[n] = false;
1.533 + }
1.534 +
1.535 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.536 + (*_flow)[a] = 0;
1.537 + }
1.538 +
1.539 + int bucket_num = 0;
1.540 + std::vector<Node> queue(_node_num);
1.541 + int qfirst = 0, qlast = 0, qsep = 0;
1.542 +
1.543 + {
1.544 + typename Digraph::template NodeMap<bool> reached(_graph, false);
1.545 +
1.546 + reached[_source] = true;
1.547 +
1.548 + bool first_set = true;
1.549 +
1.550 + for (NodeIt t(_graph); t != INVALID; ++t) {
1.551 + if (reached[t]) continue;
1.552 + _sets.push_front(std::list<int>());
1.553 +
1.554 + queue[qlast++] = t;
1.555 + reached[t] = true;
1.556 +
1.557 + while (qfirst != qlast) {
1.558 + if (qsep == qfirst) {
1.559 + ++bucket_num;
1.560 + _sets.front().push_front(bucket_num);
1.561 + _dormant[bucket_num] = !first_set;
1.562 + _first[bucket_num] = _last[bucket_num] = INVALID;
1.563 + qsep = qlast;
1.564 + }
1.565 +
1.566 + Node n = queue[qfirst++];
1.567 + addItem(n, bucket_num);
1.568 +
1.569 + for (OutArcIt a(_graph, n); a != INVALID; ++a) {
1.570 + Node u = _graph.target(a);
1.571 + if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
1.572 + reached[u] = true;
1.573 + queue[qlast++] = u;
1.574 + }
1.575 + }
1.576 + }
1.577 + first_set = false;
1.578 + }
1.579 +
1.580 + ++bucket_num;
1.581 + (*_bucket)[_source] = 0;
1.582 + _dormant[0] = true;
1.583 + }
1.584 + (*_source_set)[_source] = true;
1.585 +
1.586 + Node target = _last[_sets.back().back()];
1.587 + {
1.588 + for (InArcIt a(_graph, _source); a != INVALID; ++a) {
1.589 + if (_tolerance.positive((*_capacity)[a])) {
1.590 + Node u = _graph.source(a);
1.591 + (*_flow)[a] = (*_capacity)[a];
1.592 + (*_excess)[u] += (*_capacity)[a];
1.593 + if (!(*_active)[u] && u != _source) {
1.594 + activate(u);
1.595 + }
1.596 + }
1.597 + }
1.598 + if ((*_active)[target]) {
1.599 + deactivate(target);
1.600 + }
1.601 +
1.602 + _highest = _sets.back().begin();
1.603 + while (_highest != _sets.back().end() &&
1.604 + !(*_active)[_first[*_highest]]) {
1.605 + ++_highest;
1.606 + }
1.607 + }
1.608 +
1.609 +
1.610 + while (true) {
1.611 + while (_highest != _sets.back().end()) {
1.612 + Node n = _first[*_highest];
1.613 + Value excess = (*_excess)[n];
1.614 + int next_bucket = _node_num;
1.615 +
1.616 + int under_bucket;
1.617 + if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
1.618 + under_bucket = -1;
1.619 + } else {
1.620 + under_bucket = *(++std::list<int>::iterator(_highest));
1.621 + }
1.622 +
1.623 + for (InArcIt a(_graph, n); a != INVALID; ++a) {
1.624 + Node v = _graph.source(a);
1.625 + if (_dormant[(*_bucket)[v]]) continue;
1.626 + Value rem = (*_capacity)[a] - (*_flow)[a];
1.627 + if (!_tolerance.positive(rem)) continue;
1.628 + if ((*_bucket)[v] == under_bucket) {
1.629 + if (!(*_active)[v] && v != target) {
1.630 + activate(v);
1.631 + }
1.632 + if (!_tolerance.less(rem, excess)) {
1.633 + (*_flow)[a] += excess;
1.634 + (*_excess)[v] += excess;
1.635 + excess = 0;
1.636 + goto no_more_push;
1.637 + } else {
1.638 + excess -= rem;
1.639 + (*_excess)[v] += rem;
1.640 + (*_flow)[a] = (*_capacity)[a];
1.641 + }
1.642 + } else if (next_bucket > (*_bucket)[v]) {
1.643 + next_bucket = (*_bucket)[v];
1.644 + }
1.645 + }
1.646 +
1.647 + for (OutArcIt a(_graph, n); a != INVALID; ++a) {
1.648 + Node v = _graph.target(a);
1.649 + if (_dormant[(*_bucket)[v]]) continue;
1.650 + Value rem = (*_flow)[a];
1.651 + if (!_tolerance.positive(rem)) continue;
1.652 + if ((*_bucket)[v] == under_bucket) {
1.653 + if (!(*_active)[v] && v != target) {
1.654 + activate(v);
1.655 + }
1.656 + if (!_tolerance.less(rem, excess)) {
1.657 + (*_flow)[a] -= excess;
1.658 + (*_excess)[v] += excess;
1.659 + excess = 0;
1.660 + goto no_more_push;
1.661 + } else {
1.662 + excess -= rem;
1.663 + (*_excess)[v] += rem;
1.664 + (*_flow)[a] = 0;
1.665 + }
1.666 + } else if (next_bucket > (*_bucket)[v]) {
1.667 + next_bucket = (*_bucket)[v];
1.668 + }
1.669 + }
1.670 +
1.671 + no_more_push:
1.672 +
1.673 + (*_excess)[n] = excess;
1.674 +
1.675 + if (excess != 0) {
1.676 + if ((*_next)[n] == INVALID) {
1.677 + typename std::list<std::list<int> >::iterator new_set =
1.678 + _sets.insert(--_sets.end(), std::list<int>());
1.679 + new_set->splice(new_set->end(), _sets.back(),
1.680 + _sets.back().begin(), ++_highest);
1.681 + for (std::list<int>::iterator it = new_set->begin();
1.682 + it != new_set->end(); ++it) {
1.683 + _dormant[*it] = true;
1.684 + }
1.685 + while (_highest != _sets.back().end() &&
1.686 + !(*_active)[_first[*_highest]]) {
1.687 + ++_highest;
1.688 + }
1.689 + } else if (next_bucket == _node_num) {
1.690 + _first[(*_bucket)[n]] = (*_next)[n];
1.691 + (*_prev)[(*_next)[n]] = INVALID;
1.692 +
1.693 + std::list<std::list<int> >::iterator new_set =
1.694 + _sets.insert(--_sets.end(), std::list<int>());
1.695 +
1.696 + new_set->push_front(bucket_num);
1.697 + (*_bucket)[n] = bucket_num;
1.698 + _first[bucket_num] = _last[bucket_num] = n;
1.699 + (*_next)[n] = INVALID;
1.700 + (*_prev)[n] = INVALID;
1.701 + _dormant[bucket_num] = true;
1.702 + ++bucket_num;
1.703 +
1.704 + while (_highest != _sets.back().end() &&
1.705 + !(*_active)[_first[*_highest]]) {
1.706 + ++_highest;
1.707 + }
1.708 + } else {
1.709 + _first[*_highest] = (*_next)[n];
1.710 + (*_prev)[(*_next)[n]] = INVALID;
1.711 +
1.712 + while (next_bucket != *_highest) {
1.713 + --_highest;
1.714 + }
1.715 + if (_highest == _sets.back().begin()) {
1.716 + _sets.back().push_front(bucket_num);
1.717 + _dormant[bucket_num] = false;
1.718 + _first[bucket_num] = _last[bucket_num] = INVALID;
1.719 + ++bucket_num;
1.720 + }
1.721 + --_highest;
1.722 +
1.723 + (*_bucket)[n] = *_highest;
1.724 + (*_next)[n] = _first[*_highest];
1.725 + if (_first[*_highest] != INVALID) {
1.726 + (*_prev)[_first[*_highest]] = n;
1.727 + } else {
1.728 + _last[*_highest] = n;
1.729 + }
1.730 + _first[*_highest] = n;
1.731 + }
1.732 + } else {
1.733 +
1.734 + deactivate(n);
1.735 + if (!(*_active)[_first[*_highest]]) {
1.736 + ++_highest;
1.737 + if (_highest != _sets.back().end() &&
1.738 + !(*_active)[_first[*_highest]]) {
1.739 + _highest = _sets.back().end();
1.740 + }
1.741 + }
1.742 + }
1.743 + }
1.744 +
1.745 + if ((*_excess)[target] < _min_cut) {
1.746 + _min_cut = (*_excess)[target];
1.747 + for (NodeIt i(_graph); i != INVALID; ++i) {
1.748 + (*_min_cut_map)[i] = false;
1.749 + }
1.750 + for (std::list<int>::iterator it = _sets.back().begin();
1.751 + it != _sets.back().end(); ++it) {
1.752 + Node n = _first[*it];
1.753 + while (n != INVALID) {
1.754 + (*_min_cut_map)[n] = true;
1.755 + n = (*_next)[n];
1.756 + }
1.757 + }
1.758 + }
1.759 +
1.760 + {
1.761 + Node new_target;
1.762 + if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
1.763 + if ((*_next)[target] == INVALID) {
1.764 + _last[(*_bucket)[target]] = (*_prev)[target];
1.765 + new_target = (*_prev)[target];
1.766 + } else {
1.767 + (*_prev)[(*_next)[target]] = (*_prev)[target];
1.768 + new_target = (*_next)[target];
1.769 + }
1.770 + if ((*_prev)[target] == INVALID) {
1.771 + _first[(*_bucket)[target]] = (*_next)[target];
1.772 + } else {
1.773 + (*_next)[(*_prev)[target]] = (*_next)[target];
1.774 + }
1.775 + } else {
1.776 + _sets.back().pop_back();
1.777 + if (_sets.back().empty()) {
1.778 + _sets.pop_back();
1.779 + if (_sets.empty())
1.780 + break;
1.781 + for (std::list<int>::iterator it = _sets.back().begin();
1.782 + it != _sets.back().end(); ++it) {
1.783 + _dormant[*it] = false;
1.784 + }
1.785 + }
1.786 + new_target = _last[_sets.back().back()];
1.787 + }
1.788 +
1.789 + (*_bucket)[target] = 0;
1.790 +
1.791 + (*_source_set)[target] = true;
1.792 + for (InArcIt a(_graph, target); a != INVALID; ++a) {
1.793 + Value rem = (*_capacity)[a] - (*_flow)[a];
1.794 + if (!_tolerance.positive(rem)) continue;
1.795 + Node v = _graph.source(a);
1.796 + if (!(*_active)[v] && !(*_source_set)[v]) {
1.797 + activate(v);
1.798 + }
1.799 + (*_excess)[v] += rem;
1.800 + (*_flow)[a] = (*_capacity)[a];
1.801 + }
1.802 +
1.803 + for (OutArcIt a(_graph, target); a != INVALID; ++a) {
1.804 + Value rem = (*_flow)[a];
1.805 + if (!_tolerance.positive(rem)) continue;
1.806 + Node v = _graph.target(a);
1.807 + if (!(*_active)[v] && !(*_source_set)[v]) {
1.808 + activate(v);
1.809 + }
1.810 + (*_excess)[v] += rem;
1.811 + (*_flow)[a] = 0;
1.812 + }
1.813 +
1.814 + target = new_target;
1.815 + if ((*_active)[target]) {
1.816 + deactivate(target);
1.817 + }
1.818 +
1.819 + _highest = _sets.back().begin();
1.820 + while (_highest != _sets.back().end() &&
1.821 + !(*_active)[_first[*_highest]]) {
1.822 + ++_highest;
1.823 + }
1.824 + }
1.825 + }
1.826 + }
1.827 +
1.828 + public:
1.829 +
1.830 + /// \name Execution Control
1.831 + /// The simplest way to execute the algorithm is to use
1.832 + /// one of the member functions called \ref run().
1.833 + /// \n
1.834 + /// If you need better control on the execution,
1.835 + /// you have to call one of the \ref init() functions first, then
1.836 + /// \ref calculateOut() and/or \ref calculateIn().
1.837 +
1.838 + /// @{
1.839 +
1.840 + /// \brief Initialize the internal data structures.
1.841 + ///
1.842 + /// This function initializes the internal data structures. It creates
1.843 + /// the maps and some bucket structures for the algorithm.
1.844 + /// The first node is used as the source node for the push-relabel
1.845 + /// algorithm.
1.846 + void init() {
1.847 + init(NodeIt(_graph));
1.848 + }
1.849 +
1.850 + /// \brief Initialize the internal data structures.
1.851 + ///
1.852 + /// This function initializes the internal data structures. It creates
1.853 + /// the maps and some bucket structures for the algorithm.
1.854 + /// The given node is used as the source node for the push-relabel
1.855 + /// algorithm.
1.856 + void init(const Node& source) {
1.857 + _source = source;
1.858 +
1.859 + _node_num = countNodes(_graph);
1.860 +
1.861 + _first.resize(_node_num);
1.862 + _last.resize(_node_num);
1.863 +
1.864 + _dormant.resize(_node_num);
1.865 +
1.866 + if (!_flow) {
1.867 + _flow = new FlowMap(_graph);
1.868 + }
1.869 + if (!_next) {
1.870 + _next = new typename Digraph::template NodeMap<Node>(_graph);
1.871 + }
1.872 + if (!_prev) {
1.873 + _prev = new typename Digraph::template NodeMap<Node>(_graph);
1.874 + }
1.875 + if (!_active) {
1.876 + _active = new typename Digraph::template NodeMap<bool>(_graph);
1.877 + }
1.878 + if (!_bucket) {
1.879 + _bucket = new typename Digraph::template NodeMap<int>(_graph);
1.880 + }
1.881 + if (!_excess) {
1.882 + _excess = new ExcessMap(_graph);
1.883 + }
1.884 + if (!_source_set) {
1.885 + _source_set = new SourceSetMap(_graph);
1.886 + }
1.887 + if (!_min_cut_map) {
1.888 + _min_cut_map = new MinCutMap(_graph);
1.889 + }
1.890 +
1.891 + _min_cut = std::numeric_limits<Value>::max();
1.892 + }
1.893 +
1.894 +
1.895 + /// \brief Calculate a minimum cut with \f$ source \f$ on the
1.896 + /// source-side.
1.897 + ///
1.898 + /// This function calculates a minimum cut with \f$ source \f$ on the
1.899 + /// source-side (i.e. a set \f$ X\subsetneq V \f$ with
1.900 + /// \f$ source \in X \f$ and minimal outgoing capacity).
1.901 + ///
1.902 + /// \pre \ref init() must be called before using this function.
1.903 + void calculateOut() {
1.904 + findMinCutOut();
1.905 + }
1.906 +
1.907 + /// \brief Calculate a minimum cut with \f$ source \f$ on the
1.908 + /// sink-side.
1.909 + ///
1.910 + /// This function calculates a minimum cut with \f$ source \f$ on the
1.911 + /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with
1.912 + /// \f$ source \notin X \f$ and minimal outgoing capacity).
1.913 + ///
1.914 + /// \pre \ref init() must be called before using this function.
1.915 + void calculateIn() {
1.916 + findMinCutIn();
1.917 + }
1.918 +
1.919 +
1.920 + /// \brief Run the algorithm.
1.921 + ///
1.922 + /// This function runs the algorithm. It finds nodes \c source and
1.923 + /// \c target arbitrarily and then calls \ref init(), \ref calculateOut()
1.924 + /// and \ref calculateIn().
1.925 + void run() {
1.926 + init();
1.927 + calculateOut();
1.928 + calculateIn();
1.929 + }
1.930 +
1.931 + /// \brief Run the algorithm.
1.932 + ///
1.933 + /// This function runs the algorithm. It uses the given \c source node,
1.934 + /// finds a proper \c target node and then calls the \ref init(),
1.935 + /// \ref calculateOut() and \ref calculateIn().
1.936 + void run(const Node& s) {
1.937 + init(s);
1.938 + calculateOut();
1.939 + calculateIn();
1.940 + }
1.941 +
1.942 + /// @}
1.943 +
1.944 + /// \name Query Functions
1.945 + /// The result of the %HaoOrlin algorithm
1.946 + /// can be obtained using these functions.\n
1.947 + /// \ref run(), \ref calculateOut() or \ref calculateIn()
1.948 + /// should be called before using them.
1.949 +
1.950 + /// @{
1.951 +
1.952 + /// \brief Return the value of the minimum cut.
1.953 + ///
1.954 + /// This function returns the value of the minimum cut.
1.955 + ///
1.956 + /// \pre \ref run(), \ref calculateOut() or \ref calculateIn()
1.957 + /// must be called before using this function.
1.958 + Value minCutValue() const {
1.959 + return _min_cut;
1.960 + }
1.961 +
1.962 +
1.963 + /// \brief Return a minimum cut.
1.964 + ///
1.965 + /// This function sets \c cutMap to the characteristic vector of a
1.966 + /// minimum value cut: it will give a non-empty set \f$ X\subsetneq V \f$
1.967 + /// with minimal outgoing capacity (i.e. \c cutMap will be \c true exactly
1.968 + /// for the nodes of \f$ X \f$).
1.969 + ///
1.970 + /// \param cutMap A \ref concepts::WriteMap "writable" node map with
1.971 + /// \c bool (or convertible) value type.
1.972 + ///
1.973 + /// \return The value of the minimum cut.
1.974 + ///
1.975 + /// \pre \ref run(), \ref calculateOut() or \ref calculateIn()
1.976 + /// must be called before using this function.
1.977 + template <typename CutMap>
1.978 + Value minCutMap(CutMap& cutMap) const {
1.979 + for (NodeIt it(_graph); it != INVALID; ++it) {
1.980 + cutMap.set(it, (*_min_cut_map)[it]);
1.981 + }
1.982 + return _min_cut;
1.983 + }
1.984 +
1.985 + /// @}
1.986 +
1.987 + }; //class HaoOrlin
1.988 +
1.989 +} //namespace lemon
1.990 +
1.991 +#endif //LEMON_HAO_ORLIN_H