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