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