1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/preflow.h Thu Dec 10 17:05:35 2009 +0100
1.3 @@ -0,0 +1,965 @@
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_PREFLOW_H
1.23 +#define LEMON_PREFLOW_H
1.24 +
1.25 +#include <lemon/tolerance.h>
1.26 +#include <lemon/elevator.h>
1.27 +
1.28 +/// \file
1.29 +/// \ingroup max_flow
1.30 +/// \brief Implementation of the preflow algorithm.
1.31 +
1.32 +namespace lemon {
1.33 +
1.34 + /// \brief Default traits class of Preflow class.
1.35 + ///
1.36 + /// Default traits class of Preflow class.
1.37 + /// \tparam GR Digraph type.
1.38 + /// \tparam CAP Capacity map type.
1.39 + template <typename GR, typename CAP>
1.40 + struct PreflowDefaultTraits {
1.41 +
1.42 + /// \brief The type of the digraph the algorithm runs on.
1.43 + typedef GR Digraph;
1.44 +
1.45 + /// \brief The type of the map that stores the arc capacities.
1.46 + ///
1.47 + /// The type of the map that stores the arc capacities.
1.48 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
1.49 + typedef CAP CapacityMap;
1.50 +
1.51 + /// \brief The type of the flow values.
1.52 + typedef typename CapacityMap::Value Value;
1.53 +
1.54 + /// \brief The type of the map that stores the flow values.
1.55 + ///
1.56 + /// The type of the map that stores the flow values.
1.57 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1.58 + typedef typename Digraph::template ArcMap<Value> FlowMap;
1.59 +
1.60 + /// \brief Instantiates a FlowMap.
1.61 + ///
1.62 + /// This function instantiates a \ref FlowMap.
1.63 + /// \param digraph The digraph for which we would like to define
1.64 + /// the flow map.
1.65 + static FlowMap* createFlowMap(const Digraph& digraph) {
1.66 + return new FlowMap(digraph);
1.67 + }
1.68 +
1.69 + /// \brief The elevator type used by Preflow algorithm.
1.70 + ///
1.71 + /// The elevator type used by Preflow algorithm.
1.72 + ///
1.73 + /// \sa Elevator
1.74 + /// \sa LinkedElevator
1.75 + typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;
1.76 +
1.77 + /// \brief Instantiates an Elevator.
1.78 + ///
1.79 + /// This function instantiates an \ref Elevator.
1.80 + /// \param digraph The digraph for which we would like to define
1.81 + /// the elevator.
1.82 + /// \param max_level The maximum level of the elevator.
1.83 + static Elevator* createElevator(const Digraph& digraph, int max_level) {
1.84 + return new Elevator(digraph, max_level);
1.85 + }
1.86 +
1.87 + /// \brief The tolerance used by the algorithm
1.88 + ///
1.89 + /// The tolerance used by the algorithm to handle inexact computation.
1.90 + typedef lemon::Tolerance<Value> Tolerance;
1.91 +
1.92 + };
1.93 +
1.94 +
1.95 + /// \ingroup max_flow
1.96 + ///
1.97 + /// \brief %Preflow algorithm class.
1.98 + ///
1.99 + /// This class provides an implementation of Goldberg-Tarjan's \e preflow
1.100 + /// \e push-relabel algorithm producing a \ref max_flow
1.101 + /// "flow of maximum value" in a digraph.
1.102 + /// The preflow algorithms are the fastest known maximum
1.103 + /// flow algorithms. The current implementation use a mixture of the
1.104 + /// \e "highest label" and the \e "bound decrease" heuristics.
1.105 + /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
1.106 + ///
1.107 + /// The algorithm consists of two phases. After the first phase
1.108 + /// the maximum flow value and the minimum cut is obtained. The
1.109 + /// second phase constructs a feasible maximum flow on each arc.
1.110 + ///
1.111 + /// \tparam GR The type of the digraph the algorithm runs on.
1.112 + /// \tparam CAP The type of the capacity map. The default map
1.113 + /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
1.114 +#ifdef DOXYGEN
1.115 + template <typename GR, typename CAP, typename TR>
1.116 +#else
1.117 + template <typename GR,
1.118 + typename CAP = typename GR::template ArcMap<int>,
1.119 + typename TR = PreflowDefaultTraits<GR, CAP> >
1.120 +#endif
1.121 + class Preflow {
1.122 + public:
1.123 +
1.124 + ///The \ref PreflowDefaultTraits "traits class" of the algorithm.
1.125 + typedef TR Traits;
1.126 + ///The type of the digraph the algorithm runs on.
1.127 + typedef typename Traits::Digraph Digraph;
1.128 + ///The type of the capacity map.
1.129 + typedef typename Traits::CapacityMap CapacityMap;
1.130 + ///The type of the flow values.
1.131 + typedef typename Traits::Value Value;
1.132 +
1.133 + ///The type of the flow map.
1.134 + typedef typename Traits::FlowMap FlowMap;
1.135 + ///The type of the elevator.
1.136 + typedef typename Traits::Elevator Elevator;
1.137 + ///The type of the tolerance.
1.138 + typedef typename Traits::Tolerance Tolerance;
1.139 +
1.140 + private:
1.141 +
1.142 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.143 +
1.144 + const Digraph& _graph;
1.145 + const CapacityMap* _capacity;
1.146 +
1.147 + int _node_num;
1.148 +
1.149 + Node _source, _target;
1.150 +
1.151 + FlowMap* _flow;
1.152 + bool _local_flow;
1.153 +
1.154 + Elevator* _level;
1.155 + bool _local_level;
1.156 +
1.157 + typedef typename Digraph::template NodeMap<Value> ExcessMap;
1.158 + ExcessMap* _excess;
1.159 +
1.160 + Tolerance _tolerance;
1.161 +
1.162 + bool _phase;
1.163 +
1.164 +
1.165 + void createStructures() {
1.166 + _node_num = countNodes(_graph);
1.167 +
1.168 + if (!_flow) {
1.169 + _flow = Traits::createFlowMap(_graph);
1.170 + _local_flow = true;
1.171 + }
1.172 + if (!_level) {
1.173 + _level = Traits::createElevator(_graph, _node_num);
1.174 + _local_level = true;
1.175 + }
1.176 + if (!_excess) {
1.177 + _excess = new ExcessMap(_graph);
1.178 + }
1.179 + }
1.180 +
1.181 + void destroyStructures() {
1.182 + if (_local_flow) {
1.183 + delete _flow;
1.184 + }
1.185 + if (_local_level) {
1.186 + delete _level;
1.187 + }
1.188 + if (_excess) {
1.189 + delete _excess;
1.190 + }
1.191 + }
1.192 +
1.193 + public:
1.194 +
1.195 + typedef Preflow Create;
1.196 +
1.197 + ///\name Named Template Parameters
1.198 +
1.199 + ///@{
1.200 +
1.201 + template <typename T>
1.202 + struct SetFlowMapTraits : public Traits {
1.203 + typedef T FlowMap;
1.204 + static FlowMap *createFlowMap(const Digraph&) {
1.205 + LEMON_ASSERT(false, "FlowMap is not initialized");
1.206 + return 0; // ignore warnings
1.207 + }
1.208 + };
1.209 +
1.210 + /// \brief \ref named-templ-param "Named parameter" for setting
1.211 + /// FlowMap type
1.212 + ///
1.213 + /// \ref named-templ-param "Named parameter" for setting FlowMap
1.214 + /// type.
1.215 + template <typename T>
1.216 + struct SetFlowMap
1.217 + : public Preflow<Digraph, CapacityMap, SetFlowMapTraits<T> > {
1.218 + typedef Preflow<Digraph, CapacityMap,
1.219 + SetFlowMapTraits<T> > Create;
1.220 + };
1.221 +
1.222 + template <typename T>
1.223 + struct SetElevatorTraits : public Traits {
1.224 + typedef T Elevator;
1.225 + static Elevator *createElevator(const Digraph&, int) {
1.226 + LEMON_ASSERT(false, "Elevator is not initialized");
1.227 + return 0; // ignore warnings
1.228 + }
1.229 + };
1.230 +
1.231 + /// \brief \ref named-templ-param "Named parameter" for setting
1.232 + /// Elevator type
1.233 + ///
1.234 + /// \ref named-templ-param "Named parameter" for setting Elevator
1.235 + /// type. If this named parameter is used, then an external
1.236 + /// elevator object must be passed to the algorithm using the
1.237 + /// \ref elevator(Elevator&) "elevator()" function before calling
1.238 + /// \ref run() or \ref init().
1.239 + /// \sa SetStandardElevator
1.240 + template <typename T>
1.241 + struct SetElevator
1.242 + : public Preflow<Digraph, CapacityMap, SetElevatorTraits<T> > {
1.243 + typedef Preflow<Digraph, CapacityMap,
1.244 + SetElevatorTraits<T> > Create;
1.245 + };
1.246 +
1.247 + template <typename T>
1.248 + struct SetStandardElevatorTraits : public Traits {
1.249 + typedef T Elevator;
1.250 + static Elevator *createElevator(const Digraph& digraph, int max_level) {
1.251 + return new Elevator(digraph, max_level);
1.252 + }
1.253 + };
1.254 +
1.255 + /// \brief \ref named-templ-param "Named parameter" for setting
1.256 + /// Elevator type with automatic allocation
1.257 + ///
1.258 + /// \ref named-templ-param "Named parameter" for setting Elevator
1.259 + /// type with automatic allocation.
1.260 + /// The Elevator should have standard constructor interface to be
1.261 + /// able to automatically created by the algorithm (i.e. the
1.262 + /// digraph and the maximum level should be passed to it).
1.263 + /// However an external elevator object could also be passed to the
1.264 + /// algorithm with the \ref elevator(Elevator&) "elevator()" function
1.265 + /// before calling \ref run() or \ref init().
1.266 + /// \sa SetElevator
1.267 + template <typename T>
1.268 + struct SetStandardElevator
1.269 + : public Preflow<Digraph, CapacityMap,
1.270 + SetStandardElevatorTraits<T> > {
1.271 + typedef Preflow<Digraph, CapacityMap,
1.272 + SetStandardElevatorTraits<T> > Create;
1.273 + };
1.274 +
1.275 + /// @}
1.276 +
1.277 + protected:
1.278 +
1.279 + Preflow() {}
1.280 +
1.281 + public:
1.282 +
1.283 +
1.284 + /// \brief The constructor of the class.
1.285 + ///
1.286 + /// The constructor of the class.
1.287 + /// \param digraph The digraph the algorithm runs on.
1.288 + /// \param capacity The capacity of the arcs.
1.289 + /// \param source The source node.
1.290 + /// \param target The target node.
1.291 + Preflow(const Digraph& digraph, const CapacityMap& capacity,
1.292 + Node source, Node target)
1.293 + : _graph(digraph), _capacity(&capacity),
1.294 + _node_num(0), _source(source), _target(target),
1.295 + _flow(0), _local_flow(false),
1.296 + _level(0), _local_level(false),
1.297 + _excess(0), _tolerance(), _phase() {}
1.298 +
1.299 + /// \brief Destructor.
1.300 + ///
1.301 + /// Destructor.
1.302 + ~Preflow() {
1.303 + destroyStructures();
1.304 + }
1.305 +
1.306 + /// \brief Sets the capacity map.
1.307 + ///
1.308 + /// Sets the capacity map.
1.309 + /// \return <tt>(*this)</tt>
1.310 + Preflow& capacityMap(const CapacityMap& map) {
1.311 + _capacity = ↦
1.312 + return *this;
1.313 + }
1.314 +
1.315 + /// \brief Sets the flow map.
1.316 + ///
1.317 + /// Sets the flow map.
1.318 + /// If you don't use this function before calling \ref run() or
1.319 + /// \ref init(), an instance will be allocated automatically.
1.320 + /// The destructor deallocates this automatically allocated map,
1.321 + /// of course.
1.322 + /// \return <tt>(*this)</tt>
1.323 + Preflow& flowMap(FlowMap& map) {
1.324 + if (_local_flow) {
1.325 + delete _flow;
1.326 + _local_flow = false;
1.327 + }
1.328 + _flow = ↦
1.329 + return *this;
1.330 + }
1.331 +
1.332 + /// \brief Sets the source node.
1.333 + ///
1.334 + /// Sets the source node.
1.335 + /// \return <tt>(*this)</tt>
1.336 + Preflow& source(const Node& node) {
1.337 + _source = node;
1.338 + return *this;
1.339 + }
1.340 +
1.341 + /// \brief Sets the target node.
1.342 + ///
1.343 + /// Sets the target node.
1.344 + /// \return <tt>(*this)</tt>
1.345 + Preflow& target(const Node& node) {
1.346 + _target = node;
1.347 + return *this;
1.348 + }
1.349 +
1.350 + /// \brief Sets the elevator used by algorithm.
1.351 + ///
1.352 + /// Sets the elevator used by algorithm.
1.353 + /// If you don't use this function before calling \ref run() or
1.354 + /// \ref init(), an instance will be allocated automatically.
1.355 + /// The destructor deallocates this automatically allocated elevator,
1.356 + /// of course.
1.357 + /// \return <tt>(*this)</tt>
1.358 + Preflow& elevator(Elevator& elevator) {
1.359 + if (_local_level) {
1.360 + delete _level;
1.361 + _local_level = false;
1.362 + }
1.363 + _level = &elevator;
1.364 + return *this;
1.365 + }
1.366 +
1.367 + /// \brief Returns a const reference to the elevator.
1.368 + ///
1.369 + /// Returns a const reference to the elevator.
1.370 + ///
1.371 + /// \pre Either \ref run() or \ref init() must be called before
1.372 + /// using this function.
1.373 + const Elevator& elevator() const {
1.374 + return *_level;
1.375 + }
1.376 +
1.377 + /// \brief Sets the tolerance used by algorithm.
1.378 + ///
1.379 + /// Sets the tolerance used by algorithm.
1.380 + Preflow& tolerance(const Tolerance& tolerance) const {
1.381 + _tolerance = tolerance;
1.382 + return *this;
1.383 + }
1.384 +
1.385 + /// \brief Returns a const reference to the tolerance.
1.386 + ///
1.387 + /// Returns a const reference to the tolerance.
1.388 + const Tolerance& tolerance() const {
1.389 + return tolerance;
1.390 + }
1.391 +
1.392 + /// \name Execution Control
1.393 + /// The simplest way to execute the preflow algorithm is to use
1.394 + /// \ref run() or \ref runMinCut().\n
1.395 + /// If you need more control on the initial solution or the execution,
1.396 + /// first you have to call one of the \ref init() functions, then
1.397 + /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
1.398 +
1.399 + ///@{
1.400 +
1.401 + /// \brief Initializes the internal data structures.
1.402 + ///
1.403 + /// Initializes the internal data structures and sets the initial
1.404 + /// flow to zero on each arc.
1.405 + void init() {
1.406 + createStructures();
1.407 +
1.408 + _phase = true;
1.409 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.410 + (*_excess)[n] = 0;
1.411 + }
1.412 +
1.413 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.414 + _flow->set(e, 0);
1.415 + }
1.416 +
1.417 + typename Digraph::template NodeMap<bool> reached(_graph, false);
1.418 +
1.419 + _level->initStart();
1.420 + _level->initAddItem(_target);
1.421 +
1.422 + std::vector<Node> queue;
1.423 + reached[_source] = true;
1.424 +
1.425 + queue.push_back(_target);
1.426 + reached[_target] = true;
1.427 + while (!queue.empty()) {
1.428 + _level->initNewLevel();
1.429 + std::vector<Node> nqueue;
1.430 + for (int i = 0; i < int(queue.size()); ++i) {
1.431 + Node n = queue[i];
1.432 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.433 + Node u = _graph.source(e);
1.434 + if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
1.435 + reached[u] = true;
1.436 + _level->initAddItem(u);
1.437 + nqueue.push_back(u);
1.438 + }
1.439 + }
1.440 + }
1.441 + queue.swap(nqueue);
1.442 + }
1.443 + _level->initFinish();
1.444 +
1.445 + for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
1.446 + if (_tolerance.positive((*_capacity)[e])) {
1.447 + Node u = _graph.target(e);
1.448 + if ((*_level)[u] == _level->maxLevel()) continue;
1.449 + _flow->set(e, (*_capacity)[e]);
1.450 + (*_excess)[u] += (*_capacity)[e];
1.451 + if (u != _target && !_level->active(u)) {
1.452 + _level->activate(u);
1.453 + }
1.454 + }
1.455 + }
1.456 + }
1.457 +
1.458 + /// \brief Initializes the internal data structures using the
1.459 + /// given flow map.
1.460 + ///
1.461 + /// Initializes the internal data structures and sets the initial
1.462 + /// flow to the given \c flowMap. The \c flowMap should contain a
1.463 + /// flow or at least a preflow, i.e. at each node excluding the
1.464 + /// source node the incoming flow should greater or equal to the
1.465 + /// outgoing flow.
1.466 + /// \return \c false if the given \c flowMap is not a preflow.
1.467 + template <typename FlowMap>
1.468 + bool init(const FlowMap& flowMap) {
1.469 + createStructures();
1.470 +
1.471 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.472 + _flow->set(e, flowMap[e]);
1.473 + }
1.474 +
1.475 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.476 + Value excess = 0;
1.477 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.478 + excess += (*_flow)[e];
1.479 + }
1.480 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.481 + excess -= (*_flow)[e];
1.482 + }
1.483 + if (excess < 0 && n != _source) return false;
1.484 + (*_excess)[n] = excess;
1.485 + }
1.486 +
1.487 + typename Digraph::template NodeMap<bool> reached(_graph, false);
1.488 +
1.489 + _level->initStart();
1.490 + _level->initAddItem(_target);
1.491 +
1.492 + std::vector<Node> queue;
1.493 + reached[_source] = true;
1.494 +
1.495 + queue.push_back(_target);
1.496 + reached[_target] = true;
1.497 + while (!queue.empty()) {
1.498 + _level->initNewLevel();
1.499 + std::vector<Node> nqueue;
1.500 + for (int i = 0; i < int(queue.size()); ++i) {
1.501 + Node n = queue[i];
1.502 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.503 + Node u = _graph.source(e);
1.504 + if (!reached[u] &&
1.505 + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
1.506 + reached[u] = true;
1.507 + _level->initAddItem(u);
1.508 + nqueue.push_back(u);
1.509 + }
1.510 + }
1.511 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.512 + Node v = _graph.target(e);
1.513 + if (!reached[v] && _tolerance.positive((*_flow)[e])) {
1.514 + reached[v] = true;
1.515 + _level->initAddItem(v);
1.516 + nqueue.push_back(v);
1.517 + }
1.518 + }
1.519 + }
1.520 + queue.swap(nqueue);
1.521 + }
1.522 + _level->initFinish();
1.523 +
1.524 + for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
1.525 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.526 + if (_tolerance.positive(rem)) {
1.527 + Node u = _graph.target(e);
1.528 + if ((*_level)[u] == _level->maxLevel()) continue;
1.529 + _flow->set(e, (*_capacity)[e]);
1.530 + (*_excess)[u] += rem;
1.531 + if (u != _target && !_level->active(u)) {
1.532 + _level->activate(u);
1.533 + }
1.534 + }
1.535 + }
1.536 + for (InArcIt e(_graph, _source); e != INVALID; ++e) {
1.537 + Value rem = (*_flow)[e];
1.538 + if (_tolerance.positive(rem)) {
1.539 + Node v = _graph.source(e);
1.540 + if ((*_level)[v] == _level->maxLevel()) continue;
1.541 + _flow->set(e, 0);
1.542 + (*_excess)[v] += rem;
1.543 + if (v != _target && !_level->active(v)) {
1.544 + _level->activate(v);
1.545 + }
1.546 + }
1.547 + }
1.548 + return true;
1.549 + }
1.550 +
1.551 + /// \brief Starts the first phase of the preflow algorithm.
1.552 + ///
1.553 + /// The preflow algorithm consists of two phases, this method runs
1.554 + /// the first phase. After the first phase the maximum flow value
1.555 + /// and a minimum value cut can already be computed, although a
1.556 + /// maximum flow is not yet obtained. So after calling this method
1.557 + /// \ref flowValue() returns the value of a maximum flow and \ref
1.558 + /// minCut() returns a minimum cut.
1.559 + /// \pre One of the \ref init() functions must be called before
1.560 + /// using this function.
1.561 + void startFirstPhase() {
1.562 + _phase = true;
1.563 +
1.564 + Node n = _level->highestActive();
1.565 + int level = _level->highestActiveLevel();
1.566 + while (n != INVALID) {
1.567 + int num = _node_num;
1.568 +
1.569 + while (num > 0 && n != INVALID) {
1.570 + Value excess = (*_excess)[n];
1.571 + int new_level = _level->maxLevel();
1.572 +
1.573 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.574 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.575 + if (!_tolerance.positive(rem)) continue;
1.576 + Node v = _graph.target(e);
1.577 + if ((*_level)[v] < level) {
1.578 + if (!_level->active(v) && v != _target) {
1.579 + _level->activate(v);
1.580 + }
1.581 + if (!_tolerance.less(rem, excess)) {
1.582 + _flow->set(e, (*_flow)[e] + excess);
1.583 + (*_excess)[v] += excess;
1.584 + excess = 0;
1.585 + goto no_more_push_1;
1.586 + } else {
1.587 + excess -= rem;
1.588 + (*_excess)[v] += rem;
1.589 + _flow->set(e, (*_capacity)[e]);
1.590 + }
1.591 + } else if (new_level > (*_level)[v]) {
1.592 + new_level = (*_level)[v];
1.593 + }
1.594 + }
1.595 +
1.596 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.597 + Value rem = (*_flow)[e];
1.598 + if (!_tolerance.positive(rem)) continue;
1.599 + Node v = _graph.source(e);
1.600 + if ((*_level)[v] < level) {
1.601 + if (!_level->active(v) && v != _target) {
1.602 + _level->activate(v);
1.603 + }
1.604 + if (!_tolerance.less(rem, excess)) {
1.605 + _flow->set(e, (*_flow)[e] - excess);
1.606 + (*_excess)[v] += excess;
1.607 + excess = 0;
1.608 + goto no_more_push_1;
1.609 + } else {
1.610 + excess -= rem;
1.611 + (*_excess)[v] += rem;
1.612 + _flow->set(e, 0);
1.613 + }
1.614 + } else if (new_level > (*_level)[v]) {
1.615 + new_level = (*_level)[v];
1.616 + }
1.617 + }
1.618 +
1.619 + no_more_push_1:
1.620 +
1.621 + (*_excess)[n] = excess;
1.622 +
1.623 + if (excess != 0) {
1.624 + if (new_level + 1 < _level->maxLevel()) {
1.625 + _level->liftHighestActive(new_level + 1);
1.626 + } else {
1.627 + _level->liftHighestActiveToTop();
1.628 + }
1.629 + if (_level->emptyLevel(level)) {
1.630 + _level->liftToTop(level);
1.631 + }
1.632 + } else {
1.633 + _level->deactivate(n);
1.634 + }
1.635 +
1.636 + n = _level->highestActive();
1.637 + level = _level->highestActiveLevel();
1.638 + --num;
1.639 + }
1.640 +
1.641 + num = _node_num * 20;
1.642 + while (num > 0 && n != INVALID) {
1.643 + Value excess = (*_excess)[n];
1.644 + int new_level = _level->maxLevel();
1.645 +
1.646 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.647 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.648 + if (!_tolerance.positive(rem)) continue;
1.649 + Node v = _graph.target(e);
1.650 + if ((*_level)[v] < level) {
1.651 + if (!_level->active(v) && v != _target) {
1.652 + _level->activate(v);
1.653 + }
1.654 + if (!_tolerance.less(rem, excess)) {
1.655 + _flow->set(e, (*_flow)[e] + excess);
1.656 + (*_excess)[v] += excess;
1.657 + excess = 0;
1.658 + goto no_more_push_2;
1.659 + } else {
1.660 + excess -= rem;
1.661 + (*_excess)[v] += rem;
1.662 + _flow->set(e, (*_capacity)[e]);
1.663 + }
1.664 + } else if (new_level > (*_level)[v]) {
1.665 + new_level = (*_level)[v];
1.666 + }
1.667 + }
1.668 +
1.669 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.670 + Value rem = (*_flow)[e];
1.671 + if (!_tolerance.positive(rem)) continue;
1.672 + Node v = _graph.source(e);
1.673 + if ((*_level)[v] < level) {
1.674 + if (!_level->active(v) && v != _target) {
1.675 + _level->activate(v);
1.676 + }
1.677 + if (!_tolerance.less(rem, excess)) {
1.678 + _flow->set(e, (*_flow)[e] - excess);
1.679 + (*_excess)[v] += excess;
1.680 + excess = 0;
1.681 + goto no_more_push_2;
1.682 + } else {
1.683 + excess -= rem;
1.684 + (*_excess)[v] += rem;
1.685 + _flow->set(e, 0);
1.686 + }
1.687 + } else if (new_level > (*_level)[v]) {
1.688 + new_level = (*_level)[v];
1.689 + }
1.690 + }
1.691 +
1.692 + no_more_push_2:
1.693 +
1.694 + (*_excess)[n] = excess;
1.695 +
1.696 + if (excess != 0) {
1.697 + if (new_level + 1 < _level->maxLevel()) {
1.698 + _level->liftActiveOn(level, new_level + 1);
1.699 + } else {
1.700 + _level->liftActiveToTop(level);
1.701 + }
1.702 + if (_level->emptyLevel(level)) {
1.703 + _level->liftToTop(level);
1.704 + }
1.705 + } else {
1.706 + _level->deactivate(n);
1.707 + }
1.708 +
1.709 + while (level >= 0 && _level->activeFree(level)) {
1.710 + --level;
1.711 + }
1.712 + if (level == -1) {
1.713 + n = _level->highestActive();
1.714 + level = _level->highestActiveLevel();
1.715 + } else {
1.716 + n = _level->activeOn(level);
1.717 + }
1.718 + --num;
1.719 + }
1.720 + }
1.721 + }
1.722 +
1.723 + /// \brief Starts the second phase of the preflow algorithm.
1.724 + ///
1.725 + /// The preflow algorithm consists of two phases, this method runs
1.726 + /// the second phase. After calling one of the \ref init() functions
1.727 + /// and \ref startFirstPhase() and then \ref startSecondPhase(),
1.728 + /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
1.729 + /// value of a maximum flow, \ref minCut() returns a minimum cut
1.730 + /// \pre One of the \ref init() functions and \ref startFirstPhase()
1.731 + /// must be called before using this function.
1.732 + void startSecondPhase() {
1.733 + _phase = false;
1.734 +
1.735 + typename Digraph::template NodeMap<bool> reached(_graph);
1.736 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.737 + reached[n] = (*_level)[n] < _level->maxLevel();
1.738 + }
1.739 +
1.740 + _level->initStart();
1.741 + _level->initAddItem(_source);
1.742 +
1.743 + std::vector<Node> queue;
1.744 + queue.push_back(_source);
1.745 + reached[_source] = true;
1.746 +
1.747 + while (!queue.empty()) {
1.748 + _level->initNewLevel();
1.749 + std::vector<Node> nqueue;
1.750 + for (int i = 0; i < int(queue.size()); ++i) {
1.751 + Node n = queue[i];
1.752 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.753 + Node v = _graph.target(e);
1.754 + if (!reached[v] && _tolerance.positive((*_flow)[e])) {
1.755 + reached[v] = true;
1.756 + _level->initAddItem(v);
1.757 + nqueue.push_back(v);
1.758 + }
1.759 + }
1.760 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.761 + Node u = _graph.source(e);
1.762 + if (!reached[u] &&
1.763 + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
1.764 + reached[u] = true;
1.765 + _level->initAddItem(u);
1.766 + nqueue.push_back(u);
1.767 + }
1.768 + }
1.769 + }
1.770 + queue.swap(nqueue);
1.771 + }
1.772 + _level->initFinish();
1.773 +
1.774 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.775 + if (!reached[n]) {
1.776 + _level->dirtyTopButOne(n);
1.777 + } else if ((*_excess)[n] > 0 && _target != n) {
1.778 + _level->activate(n);
1.779 + }
1.780 + }
1.781 +
1.782 + Node n;
1.783 + while ((n = _level->highestActive()) != INVALID) {
1.784 + Value excess = (*_excess)[n];
1.785 + int level = _level->highestActiveLevel();
1.786 + int new_level = _level->maxLevel();
1.787 +
1.788 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.789 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.790 + if (!_tolerance.positive(rem)) continue;
1.791 + Node v = _graph.target(e);
1.792 + if ((*_level)[v] < level) {
1.793 + if (!_level->active(v) && v != _source) {
1.794 + _level->activate(v);
1.795 + }
1.796 + if (!_tolerance.less(rem, excess)) {
1.797 + _flow->set(e, (*_flow)[e] + excess);
1.798 + (*_excess)[v] += excess;
1.799 + excess = 0;
1.800 + goto no_more_push;
1.801 + } else {
1.802 + excess -= rem;
1.803 + (*_excess)[v] += rem;
1.804 + _flow->set(e, (*_capacity)[e]);
1.805 + }
1.806 + } else if (new_level > (*_level)[v]) {
1.807 + new_level = (*_level)[v];
1.808 + }
1.809 + }
1.810 +
1.811 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.812 + Value rem = (*_flow)[e];
1.813 + if (!_tolerance.positive(rem)) continue;
1.814 + Node v = _graph.source(e);
1.815 + if ((*_level)[v] < level) {
1.816 + if (!_level->active(v) && v != _source) {
1.817 + _level->activate(v);
1.818 + }
1.819 + if (!_tolerance.less(rem, excess)) {
1.820 + _flow->set(e, (*_flow)[e] - excess);
1.821 + (*_excess)[v] += excess;
1.822 + excess = 0;
1.823 + goto no_more_push;
1.824 + } else {
1.825 + excess -= rem;
1.826 + (*_excess)[v] += rem;
1.827 + _flow->set(e, 0);
1.828 + }
1.829 + } else if (new_level > (*_level)[v]) {
1.830 + new_level = (*_level)[v];
1.831 + }
1.832 + }
1.833 +
1.834 + no_more_push:
1.835 +
1.836 + (*_excess)[n] = excess;
1.837 +
1.838 + if (excess != 0) {
1.839 + if (new_level + 1 < _level->maxLevel()) {
1.840 + _level->liftHighestActive(new_level + 1);
1.841 + } else {
1.842 + // Calculation error
1.843 + _level->liftHighestActiveToTop();
1.844 + }
1.845 + if (_level->emptyLevel(level)) {
1.846 + // Calculation error
1.847 + _level->liftToTop(level);
1.848 + }
1.849 + } else {
1.850 + _level->deactivate(n);
1.851 + }
1.852 +
1.853 + }
1.854 + }
1.855 +
1.856 + /// \brief Runs the preflow algorithm.
1.857 + ///
1.858 + /// Runs the preflow algorithm.
1.859 + /// \note pf.run() is just a shortcut of the following code.
1.860 + /// \code
1.861 + /// pf.init();
1.862 + /// pf.startFirstPhase();
1.863 + /// pf.startSecondPhase();
1.864 + /// \endcode
1.865 + void run() {
1.866 + init();
1.867 + startFirstPhase();
1.868 + startSecondPhase();
1.869 + }
1.870 +
1.871 + /// \brief Runs the preflow algorithm to compute the minimum cut.
1.872 + ///
1.873 + /// Runs the preflow algorithm to compute the minimum cut.
1.874 + /// \note pf.runMinCut() is just a shortcut of the following code.
1.875 + /// \code
1.876 + /// pf.init();
1.877 + /// pf.startFirstPhase();
1.878 + /// \endcode
1.879 + void runMinCut() {
1.880 + init();
1.881 + startFirstPhase();
1.882 + }
1.883 +
1.884 + /// @}
1.885 +
1.886 + /// \name Query Functions
1.887 + /// The results of the preflow algorithm can be obtained using these
1.888 + /// functions.\n
1.889 + /// Either one of the \ref run() "run*()" functions or one of the
1.890 + /// \ref startFirstPhase() "start*()" functions should be called
1.891 + /// before using them.
1.892 +
1.893 + ///@{
1.894 +
1.895 + /// \brief Returns the value of the maximum flow.
1.896 + ///
1.897 + /// Returns the value of the maximum flow by returning the excess
1.898 + /// of the target node. This value equals to the value of
1.899 + /// the maximum flow already after the first phase of the algorithm.
1.900 + ///
1.901 + /// \pre Either \ref run() or \ref init() must be called before
1.902 + /// using this function.
1.903 + Value flowValue() const {
1.904 + return (*_excess)[_target];
1.905 + }
1.906 +
1.907 + /// \brief Returns the flow value on the given arc.
1.908 + ///
1.909 + /// Returns the flow value on the given arc. This method can
1.910 + /// be called after the second phase of the algorithm.
1.911 + ///
1.912 + /// \pre Either \ref run() or \ref init() must be called before
1.913 + /// using this function.
1.914 + Value flow(const Arc& arc) const {
1.915 + return (*_flow)[arc];
1.916 + }
1.917 +
1.918 + /// \brief Returns a const reference to the flow map.
1.919 + ///
1.920 + /// Returns a const reference to the arc map storing the found flow.
1.921 + /// This method can be called after the second phase of the algorithm.
1.922 + ///
1.923 + /// \pre Either \ref run() or \ref init() must be called before
1.924 + /// using this function.
1.925 + const FlowMap& flowMap() const {
1.926 + return *_flow;
1.927 + }
1.928 +
1.929 + /// \brief Returns \c true when the node is on the source side of the
1.930 + /// minimum cut.
1.931 + ///
1.932 + /// Returns true when the node is on the source side of the found
1.933 + /// minimum cut. This method can be called both after running \ref
1.934 + /// startFirstPhase() and \ref startSecondPhase().
1.935 + ///
1.936 + /// \pre Either \ref run() or \ref init() must be called before
1.937 + /// using this function.
1.938 + bool minCut(const Node& node) const {
1.939 + return ((*_level)[node] == _level->maxLevel()) == _phase;
1.940 + }
1.941 +
1.942 + /// \brief Gives back a minimum value cut.
1.943 + ///
1.944 + /// Sets \c cutMap to the characteristic vector of a minimum value
1.945 + /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
1.946 + /// node map with \c bool (or convertible) value type.
1.947 + ///
1.948 + /// This method can be called both after running \ref startFirstPhase()
1.949 + /// and \ref startSecondPhase(). The result after the second phase
1.950 + /// could be slightly different if inexact computation is used.
1.951 + ///
1.952 + /// \note This function calls \ref minCut() for each node, so it runs in
1.953 + /// O(n) time.
1.954 + ///
1.955 + /// \pre Either \ref run() or \ref init() must be called before
1.956 + /// using this function.
1.957 + template <typename CutMap>
1.958 + void minCutMap(CutMap& cutMap) const {
1.959 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.960 + cutMap.set(n, minCut(n));
1.961 + }
1.962 + }
1.963 +
1.964 + /// @}
1.965 + };
1.966 +}
1.967 +
1.968 +#endif