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