Port preflow push max flow alg. from svn -r3516 (#176)
Namely,
- port the files
- apply the migrate script
- apply the unify script
- break the long lines in lemon/preflow.h
- convert the .dim test file to .lgf
- fix compilation problems
1.1 --- a/lemon/Makefile.am Fri Nov 21 10:49:39 2008 +0000
1.2 +++ b/lemon/Makefile.am Fri Nov 21 14:11:29 2008 +0000
1.3 @@ -42,6 +42,7 @@
1.4 lemon/max_matching.h \
1.5 lemon/nauty_reader.h \
1.6 lemon/path.h \
1.7 + lemon/preflow.h \
1.8 lemon/random.h \
1.9 lemon/smart_graph.h \
1.10 lemon/suurballe.h \
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/lemon/preflow.h Fri Nov 21 14:11:29 2008 +0000
2.3 @@ -0,0 +1,927 @@
2.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
2.5 + *
2.6 + * This file is a part of LEMON, a generic C++ optimization library.
2.7 + *
2.8 + * Copyright (C) 2003-2008
2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
2.11 + *
2.12 + * Permission to use, modify and distribute this software is granted
2.13 + * provided that this copyright notice appears in all copies. For
2.14 + * precise terms see the accompanying LICENSE file.
2.15 + *
2.16 + * This software is provided "AS IS" with no warranty of any kind,
2.17 + * express or implied, and with no claim as to its suitability for any
2.18 + * purpose.
2.19 + *
2.20 + */
2.21 +
2.22 +#ifndef LEMON_PREFLOW_H
2.23 +#define LEMON_PREFLOW_H
2.24 +
2.25 +#include <lemon/error.h>
2.26 +#include <lemon/tolerance.h>
2.27 +#include <lemon/elevator.h>
2.28 +
2.29 +/// \file
2.30 +/// \ingroup max_flow
2.31 +/// \brief Implementation of the preflow algorithm.
2.32 +
2.33 +namespace lemon {
2.34 +
2.35 + /// \brief Default traits class of Preflow class.
2.36 + ///
2.37 + /// Default traits class of Preflow class.
2.38 + /// \param _Graph Digraph type.
2.39 + /// \param _CapacityMap Type of capacity map.
2.40 + template <typename _Graph, typename _CapacityMap>
2.41 + struct PreflowDefaultTraits {
2.42 +
2.43 + /// \brief The digraph type the algorithm runs on.
2.44 + typedef _Graph Digraph;
2.45 +
2.46 + /// \brief The type of the map that stores the arc capacities.
2.47 + ///
2.48 + /// The type of the map that stores the arc capacities.
2.49 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
2.50 + typedef _CapacityMap CapacityMap;
2.51 +
2.52 + /// \brief The type of the length of the arcs.
2.53 + typedef typename CapacityMap::Value Value;
2.54 +
2.55 + /// \brief The map type that stores the flow values.
2.56 + ///
2.57 + /// The map type that stores the flow values.
2.58 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
2.59 + typedef typename Digraph::template ArcMap<Value> FlowMap;
2.60 +
2.61 + /// \brief Instantiates a FlowMap.
2.62 + ///
2.63 + /// This function instantiates a \ref FlowMap.
2.64 + /// \param digraph The digraph, to which we would like to define
2.65 + /// the flow map.
2.66 + static FlowMap* createFlowMap(const Digraph& digraph) {
2.67 + return new FlowMap(digraph);
2.68 + }
2.69 +
2.70 + /// \brief The eleavator type used by Preflow algorithm.
2.71 + ///
2.72 + /// The elevator type used by Preflow algorithm.
2.73 + ///
2.74 + /// \sa Elevator
2.75 + /// \sa LinkedElevator
2.76 + typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;
2.77 +
2.78 + /// \brief Instantiates an Elevator.
2.79 + ///
2.80 + /// This function instantiates a \ref Elevator.
2.81 + /// \param digraph The digraph, to which we would like to define
2.82 + /// the elevator.
2.83 + /// \param max_level The maximum level of the elevator.
2.84 + static Elevator* createElevator(const Digraph& digraph, int max_level) {
2.85 + return new Elevator(digraph, max_level);
2.86 + }
2.87 +
2.88 + /// \brief The tolerance used by the algorithm
2.89 + ///
2.90 + /// The tolerance used by the algorithm to handle inexact computation.
2.91 + typedef lemon::Tolerance<Value> Tolerance;
2.92 +
2.93 + };
2.94 +
2.95 +
2.96 + /// \ingroup max_flow
2.97 + ///
2.98 + /// \brief %Preflow algorithms class.
2.99 + ///
2.100 + /// This class provides an implementation of the Goldberg's \e
2.101 + /// preflow \e algorithm producing a flow of maximum value in a
2.102 + /// digraph. The preflow algorithms are the fastest known max
2.103 + /// flow algorithms. The current implementation use a mixture of the
2.104 + /// \e "highest label" and the \e "bound decrease" heuristics.
2.105 + /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
2.106 + ///
2.107 + /// The algorithm consists from two phases. After the first phase
2.108 + /// the maximal flow value and the minimum cut can be obtained. The
2.109 + /// second phase constructs the feasible maximum flow on each arc.
2.110 + ///
2.111 + /// \param _Graph The digraph type the algorithm runs on.
2.112 + /// \param _CapacityMap The flow map type.
2.113 + /// \param _Traits Traits class to set various data types used by
2.114 + /// the algorithm. The default traits class is \ref
2.115 + /// PreflowDefaultTraits. See \ref PreflowDefaultTraits for the
2.116 + /// documentation of a %Preflow traits class.
2.117 + ///
2.118 + ///\author Jacint Szabo and Balazs Dezso
2.119 +#ifdef DOXYGEN
2.120 + template <typename _Graph, typename _CapacityMap, typename _Traits>
2.121 +#else
2.122 + template <typename _Graph,
2.123 + typename _CapacityMap = typename _Graph::template ArcMap<int>,
2.124 + typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
2.125 +#endif
2.126 + class Preflow {
2.127 + public:
2.128 +
2.129 + typedef _Traits Traits;
2.130 + typedef typename Traits::Digraph Digraph;
2.131 + typedef typename Traits::CapacityMap CapacityMap;
2.132 + typedef typename Traits::Value Value;
2.133 +
2.134 + typedef typename Traits::FlowMap FlowMap;
2.135 + typedef typename Traits::Elevator Elevator;
2.136 + typedef typename Traits::Tolerance Tolerance;
2.137 +
2.138 + /// \brief \ref Exception for uninitialized parameters.
2.139 + ///
2.140 + /// This error represents problems in the initialization
2.141 + /// of the parameters of the algorithms.
2.142 + class UninitializedParameter : public lemon::Exception {
2.143 + public:
2.144 + virtual const char* what() const throw() {
2.145 + return "lemon::Preflow::UninitializedParameter";
2.146 + }
2.147 + };
2.148 +
2.149 + private:
2.150 +
2.151 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
2.152 +
2.153 + const Digraph& _graph;
2.154 + const CapacityMap* _capacity;
2.155 +
2.156 + int _node_num;
2.157 +
2.158 + Node _source, _target;
2.159 +
2.160 + FlowMap* _flow;
2.161 + bool _local_flow;
2.162 +
2.163 + Elevator* _level;
2.164 + bool _local_level;
2.165 +
2.166 + typedef typename Digraph::template NodeMap<Value> ExcessMap;
2.167 + ExcessMap* _excess;
2.168 +
2.169 + Tolerance _tolerance;
2.170 +
2.171 + bool _phase;
2.172 +
2.173 +
2.174 + void createStructures() {
2.175 + _node_num = countNodes(_graph);
2.176 +
2.177 + if (!_flow) {
2.178 + _flow = Traits::createFlowMap(_graph);
2.179 + _local_flow = true;
2.180 + }
2.181 + if (!_level) {
2.182 + _level = Traits::createElevator(_graph, _node_num);
2.183 + _local_level = true;
2.184 + }
2.185 + if (!_excess) {
2.186 + _excess = new ExcessMap(_graph);
2.187 + }
2.188 + }
2.189 +
2.190 + void destroyStructures() {
2.191 + if (_local_flow) {
2.192 + delete _flow;
2.193 + }
2.194 + if (_local_level) {
2.195 + delete _level;
2.196 + }
2.197 + if (_excess) {
2.198 + delete _excess;
2.199 + }
2.200 + }
2.201 +
2.202 + public:
2.203 +
2.204 + typedef Preflow Create;
2.205 +
2.206 + ///\name Named template parameters
2.207 +
2.208 + ///@{
2.209 +
2.210 + template <typename _FlowMap>
2.211 + struct DefFlowMapTraits : public Traits {
2.212 + typedef _FlowMap FlowMap;
2.213 + static FlowMap *createFlowMap(const Digraph&) {
2.214 + throw UninitializedParameter();
2.215 + }
2.216 + };
2.217 +
2.218 + /// \brief \ref named-templ-param "Named parameter" for setting
2.219 + /// FlowMap type
2.220 + ///
2.221 + /// \ref named-templ-param "Named parameter" for setting FlowMap
2.222 + /// type
2.223 + template <typename _FlowMap>
2.224 + struct DefFlowMap
2.225 + : public Preflow<Digraph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
2.226 + typedef Preflow<Digraph, CapacityMap,
2.227 + DefFlowMapTraits<_FlowMap> > Create;
2.228 + };
2.229 +
2.230 + template <typename _Elevator>
2.231 + struct DefElevatorTraits : public Traits {
2.232 + typedef _Elevator Elevator;
2.233 + static Elevator *createElevator(const Digraph&, int) {
2.234 + throw UninitializedParameter();
2.235 + }
2.236 + };
2.237 +
2.238 + /// \brief \ref named-templ-param "Named parameter" for setting
2.239 + /// Elevator type
2.240 + ///
2.241 + /// \ref named-templ-param "Named parameter" for setting Elevator
2.242 + /// type
2.243 + template <typename _Elevator>
2.244 + struct DefElevator
2.245 + : public Preflow<Digraph, CapacityMap, DefElevatorTraits<_Elevator> > {
2.246 + typedef Preflow<Digraph, CapacityMap,
2.247 + DefElevatorTraits<_Elevator> > Create;
2.248 + };
2.249 +
2.250 + template <typename _Elevator>
2.251 + struct DefStandardElevatorTraits : public Traits {
2.252 + typedef _Elevator Elevator;
2.253 + static Elevator *createElevator(const Digraph& digraph, int max_level) {
2.254 + return new Elevator(digraph, max_level);
2.255 + }
2.256 + };
2.257 +
2.258 + /// \brief \ref named-templ-param "Named parameter" for setting
2.259 + /// Elevator type
2.260 + ///
2.261 + /// \ref named-templ-param "Named parameter" for setting Elevator
2.262 + /// type. The Elevator should be standard constructor interface, ie.
2.263 + /// the digraph and the maximum level should be passed to it.
2.264 + template <typename _Elevator>
2.265 + struct DefStandardElevator
2.266 + : public Preflow<Digraph, CapacityMap,
2.267 + DefStandardElevatorTraits<_Elevator> > {
2.268 + typedef Preflow<Digraph, CapacityMap,
2.269 + DefStandardElevatorTraits<_Elevator> > Create;
2.270 + };
2.271 +
2.272 + /// @}
2.273 +
2.274 + protected:
2.275 +
2.276 + Preflow() {}
2.277 +
2.278 + public:
2.279 +
2.280 +
2.281 + /// \brief The constructor of the class.
2.282 + ///
2.283 + /// The constructor of the class.
2.284 + /// \param digraph The digraph the algorithm runs on.
2.285 + /// \param capacity The capacity of the arcs.
2.286 + /// \param source The source node.
2.287 + /// \param target The target node.
2.288 + Preflow(const Digraph& digraph, const CapacityMap& capacity,
2.289 + Node source, Node target)
2.290 + : _graph(digraph), _capacity(&capacity),
2.291 + _node_num(0), _source(source), _target(target),
2.292 + _flow(0), _local_flow(false),
2.293 + _level(0), _local_level(false),
2.294 + _excess(0), _tolerance(), _phase() {}
2.295 +
2.296 + /// \brief Destrcutor.
2.297 + ///
2.298 + /// Destructor.
2.299 + ~Preflow() {
2.300 + destroyStructures();
2.301 + }
2.302 +
2.303 + /// \brief Sets the capacity map.
2.304 + ///
2.305 + /// Sets the capacity map.
2.306 + /// \return \c (*this)
2.307 + Preflow& capacityMap(const CapacityMap& map) {
2.308 + _capacity = ↦
2.309 + return *this;
2.310 + }
2.311 +
2.312 + /// \brief Sets the flow map.
2.313 + ///
2.314 + /// Sets the flow map.
2.315 + /// \return \c (*this)
2.316 + Preflow& flowMap(FlowMap& map) {
2.317 + if (_local_flow) {
2.318 + delete _flow;
2.319 + _local_flow = false;
2.320 + }
2.321 + _flow = ↦
2.322 + return *this;
2.323 + }
2.324 +
2.325 + /// \brief Returns the flow map.
2.326 + ///
2.327 + /// \return The flow map.
2.328 + const FlowMap& flowMap() {
2.329 + return *_flow;
2.330 + }
2.331 +
2.332 + /// \brief Sets the elevator.
2.333 + ///
2.334 + /// Sets the elevator.
2.335 + /// \return \c (*this)
2.336 + Preflow& elevator(Elevator& elevator) {
2.337 + if (_local_level) {
2.338 + delete _level;
2.339 + _local_level = false;
2.340 + }
2.341 + _level = &elevator;
2.342 + return *this;
2.343 + }
2.344 +
2.345 + /// \brief Returns the elevator.
2.346 + ///
2.347 + /// \return The elevator.
2.348 + const Elevator& elevator() {
2.349 + return *_level;
2.350 + }
2.351 +
2.352 + /// \brief Sets the source node.
2.353 + ///
2.354 + /// Sets the source node.
2.355 + /// \return \c (*this)
2.356 + Preflow& source(const Node& node) {
2.357 + _source = node;
2.358 + return *this;
2.359 + }
2.360 +
2.361 + /// \brief Sets the target node.
2.362 + ///
2.363 + /// Sets the target node.
2.364 + /// \return \c (*this)
2.365 + Preflow& target(const Node& node) {
2.366 + _target = node;
2.367 + return *this;
2.368 + }
2.369 +
2.370 + /// \brief Sets the tolerance used by algorithm.
2.371 + ///
2.372 + /// Sets the tolerance used by algorithm.
2.373 + Preflow& tolerance(const Tolerance& tolerance) const {
2.374 + _tolerance = tolerance;
2.375 + return *this;
2.376 + }
2.377 +
2.378 + /// \brief Returns the tolerance used by algorithm.
2.379 + ///
2.380 + /// Returns the tolerance used by algorithm.
2.381 + const Tolerance& tolerance() const {
2.382 + return tolerance;
2.383 + }
2.384 +
2.385 + /// \name Execution control The simplest way to execute the
2.386 + /// algorithm is to use one of the member functions called \c
2.387 + /// run().
2.388 + /// \n
2.389 + /// If you need more control on initial solution or
2.390 + /// execution then you have to call one \ref init() function and then
2.391 + /// the startFirstPhase() and if you need the startSecondPhase().
2.392 +
2.393 + ///@{
2.394 +
2.395 + /// \brief Initializes the internal data structures.
2.396 + ///
2.397 + /// Initializes the internal data structures.
2.398 + ///
2.399 + void init() {
2.400 + createStructures();
2.401 +
2.402 + _phase = true;
2.403 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.404 + _excess->set(n, 0);
2.405 + }
2.406 +
2.407 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.408 + _flow->set(e, 0);
2.409 + }
2.410 +
2.411 + typename Digraph::template NodeMap<bool> reached(_graph, false);
2.412 +
2.413 + _level->initStart();
2.414 + _level->initAddItem(_target);
2.415 +
2.416 + std::vector<Node> queue;
2.417 + reached.set(_source, true);
2.418 +
2.419 + queue.push_back(_target);
2.420 + reached.set(_target, true);
2.421 + while (!queue.empty()) {
2.422 + _level->initNewLevel();
2.423 + std::vector<Node> nqueue;
2.424 + for (int i = 0; i < int(queue.size()); ++i) {
2.425 + Node n = queue[i];
2.426 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.427 + Node u = _graph.source(e);
2.428 + if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
2.429 + reached.set(u, true);
2.430 + _level->initAddItem(u);
2.431 + nqueue.push_back(u);
2.432 + }
2.433 + }
2.434 + }
2.435 + queue.swap(nqueue);
2.436 + }
2.437 + _level->initFinish();
2.438 +
2.439 + for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
2.440 + if (_tolerance.positive((*_capacity)[e])) {
2.441 + Node u = _graph.target(e);
2.442 + if ((*_level)[u] == _level->maxLevel()) continue;
2.443 + _flow->set(e, (*_capacity)[e]);
2.444 + _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
2.445 + if (u != _target && !_level->active(u)) {
2.446 + _level->activate(u);
2.447 + }
2.448 + }
2.449 + }
2.450 + }
2.451 +
2.452 + /// \brief Initializes the internal data structures.
2.453 + ///
2.454 + /// Initializes the internal data structures and sets the initial
2.455 + /// flow to the given \c flowMap. The \c flowMap should contain a
2.456 + /// flow or at least a preflow, ie. in each node excluding the
2.457 + /// target the incoming flow should greater or equal to the
2.458 + /// outgoing flow.
2.459 + /// \return %False when the given \c flowMap is not a preflow.
2.460 + template <typename FlowMap>
2.461 + bool flowInit(const FlowMap& flowMap) {
2.462 + createStructures();
2.463 +
2.464 + for (ArcIt e(_graph); e != INVALID; ++e) {
2.465 + _flow->set(e, flowMap[e]);
2.466 + }
2.467 +
2.468 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.469 + Value excess = 0;
2.470 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.471 + excess += (*_flow)[e];
2.472 + }
2.473 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2.474 + excess -= (*_flow)[e];
2.475 + }
2.476 + if (excess < 0 && n != _source) return false;
2.477 + _excess->set(n, excess);
2.478 + }
2.479 +
2.480 + typename Digraph::template NodeMap<bool> reached(_graph, false);
2.481 +
2.482 + _level->initStart();
2.483 + _level->initAddItem(_target);
2.484 +
2.485 + std::vector<Node> queue;
2.486 + reached.set(_source, true);
2.487 +
2.488 + queue.push_back(_target);
2.489 + reached.set(_target, true);
2.490 + while (!queue.empty()) {
2.491 + _level->initNewLevel();
2.492 + std::vector<Node> nqueue;
2.493 + for (int i = 0; i < int(queue.size()); ++i) {
2.494 + Node n = queue[i];
2.495 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.496 + Node u = _graph.source(e);
2.497 + if (!reached[u] &&
2.498 + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
2.499 + reached.set(u, true);
2.500 + _level->initAddItem(u);
2.501 + nqueue.push_back(u);
2.502 + }
2.503 + }
2.504 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2.505 + Node v = _graph.target(e);
2.506 + if (!reached[v] && _tolerance.positive((*_flow)[e])) {
2.507 + reached.set(v, true);
2.508 + _level->initAddItem(v);
2.509 + nqueue.push_back(v);
2.510 + }
2.511 + }
2.512 + }
2.513 + queue.swap(nqueue);
2.514 + }
2.515 + _level->initFinish();
2.516 +
2.517 + for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
2.518 + Value rem = (*_capacity)[e] - (*_flow)[e];
2.519 + if (_tolerance.positive(rem)) {
2.520 + Node u = _graph.target(e);
2.521 + if ((*_level)[u] == _level->maxLevel()) continue;
2.522 + _flow->set(e, (*_capacity)[e]);
2.523 + _excess->set(u, (*_excess)[u] + rem);
2.524 + if (u != _target && !_level->active(u)) {
2.525 + _level->activate(u);
2.526 + }
2.527 + }
2.528 + }
2.529 + for (InArcIt e(_graph, _source); e != INVALID; ++e) {
2.530 + Value rem = (*_flow)[e];
2.531 + if (_tolerance.positive(rem)) {
2.532 + Node v = _graph.source(e);
2.533 + if ((*_level)[v] == _level->maxLevel()) continue;
2.534 + _flow->set(e, 0);
2.535 + _excess->set(v, (*_excess)[v] + rem);
2.536 + if (v != _target && !_level->active(v)) {
2.537 + _level->activate(v);
2.538 + }
2.539 + }
2.540 + }
2.541 + return true;
2.542 + }
2.543 +
2.544 + /// \brief Starts the first phase of the preflow algorithm.
2.545 + ///
2.546 + /// The preflow algorithm consists of two phases, this method runs
2.547 + /// the first phase. After the first phase the maximum flow value
2.548 + /// and a minimum value cut can already be computed, although a
2.549 + /// maximum flow is not yet obtained. So after calling this method
2.550 + /// \ref flowValue() returns the value of a maximum flow and \ref
2.551 + /// minCut() returns a minimum cut.
2.552 + /// \pre One of the \ref init() functions should be called.
2.553 + void startFirstPhase() {
2.554 + _phase = true;
2.555 +
2.556 + Node n = _level->highestActive();
2.557 + int level = _level->highestActiveLevel();
2.558 + while (n != INVALID) {
2.559 + int num = _node_num;
2.560 +
2.561 + while (num > 0 && n != INVALID) {
2.562 + Value excess = (*_excess)[n];
2.563 + int new_level = _level->maxLevel();
2.564 +
2.565 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2.566 + Value rem = (*_capacity)[e] - (*_flow)[e];
2.567 + if (!_tolerance.positive(rem)) continue;
2.568 + Node v = _graph.target(e);
2.569 + if ((*_level)[v] < level) {
2.570 + if (!_level->active(v) && v != _target) {
2.571 + _level->activate(v);
2.572 + }
2.573 + if (!_tolerance.less(rem, excess)) {
2.574 + _flow->set(e, (*_flow)[e] + excess);
2.575 + _excess->set(v, (*_excess)[v] + excess);
2.576 + excess = 0;
2.577 + goto no_more_push_1;
2.578 + } else {
2.579 + excess -= rem;
2.580 + _excess->set(v, (*_excess)[v] + rem);
2.581 + _flow->set(e, (*_capacity)[e]);
2.582 + }
2.583 + } else if (new_level > (*_level)[v]) {
2.584 + new_level = (*_level)[v];
2.585 + }
2.586 + }
2.587 +
2.588 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.589 + Value rem = (*_flow)[e];
2.590 + if (!_tolerance.positive(rem)) continue;
2.591 + Node v = _graph.source(e);
2.592 + if ((*_level)[v] < level) {
2.593 + if (!_level->active(v) && v != _target) {
2.594 + _level->activate(v);
2.595 + }
2.596 + if (!_tolerance.less(rem, excess)) {
2.597 + _flow->set(e, (*_flow)[e] - excess);
2.598 + _excess->set(v, (*_excess)[v] + excess);
2.599 + excess = 0;
2.600 + goto no_more_push_1;
2.601 + } else {
2.602 + excess -= rem;
2.603 + _excess->set(v, (*_excess)[v] + rem);
2.604 + _flow->set(e, 0);
2.605 + }
2.606 + } else if (new_level > (*_level)[v]) {
2.607 + new_level = (*_level)[v];
2.608 + }
2.609 + }
2.610 +
2.611 + no_more_push_1:
2.612 +
2.613 + _excess->set(n, excess);
2.614 +
2.615 + if (excess != 0) {
2.616 + if (new_level + 1 < _level->maxLevel()) {
2.617 + _level->liftHighestActive(new_level + 1);
2.618 + } else {
2.619 + _level->liftHighestActiveToTop();
2.620 + }
2.621 + if (_level->emptyLevel(level)) {
2.622 + _level->liftToTop(level);
2.623 + }
2.624 + } else {
2.625 + _level->deactivate(n);
2.626 + }
2.627 +
2.628 + n = _level->highestActive();
2.629 + level = _level->highestActiveLevel();
2.630 + --num;
2.631 + }
2.632 +
2.633 + num = _node_num * 20;
2.634 + while (num > 0 && n != INVALID) {
2.635 + Value excess = (*_excess)[n];
2.636 + int new_level = _level->maxLevel();
2.637 +
2.638 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2.639 + Value rem = (*_capacity)[e] - (*_flow)[e];
2.640 + if (!_tolerance.positive(rem)) continue;
2.641 + Node v = _graph.target(e);
2.642 + if ((*_level)[v] < level) {
2.643 + if (!_level->active(v) && v != _target) {
2.644 + _level->activate(v);
2.645 + }
2.646 + if (!_tolerance.less(rem, excess)) {
2.647 + _flow->set(e, (*_flow)[e] + excess);
2.648 + _excess->set(v, (*_excess)[v] + excess);
2.649 + excess = 0;
2.650 + goto no_more_push_2;
2.651 + } else {
2.652 + excess -= rem;
2.653 + _excess->set(v, (*_excess)[v] + rem);
2.654 + _flow->set(e, (*_capacity)[e]);
2.655 + }
2.656 + } else if (new_level > (*_level)[v]) {
2.657 + new_level = (*_level)[v];
2.658 + }
2.659 + }
2.660 +
2.661 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.662 + Value rem = (*_flow)[e];
2.663 + if (!_tolerance.positive(rem)) continue;
2.664 + Node v = _graph.source(e);
2.665 + if ((*_level)[v] < level) {
2.666 + if (!_level->active(v) && v != _target) {
2.667 + _level->activate(v);
2.668 + }
2.669 + if (!_tolerance.less(rem, excess)) {
2.670 + _flow->set(e, (*_flow)[e] - excess);
2.671 + _excess->set(v, (*_excess)[v] + excess);
2.672 + excess = 0;
2.673 + goto no_more_push_2;
2.674 + } else {
2.675 + excess -= rem;
2.676 + _excess->set(v, (*_excess)[v] + rem);
2.677 + _flow->set(e, 0);
2.678 + }
2.679 + } else if (new_level > (*_level)[v]) {
2.680 + new_level = (*_level)[v];
2.681 + }
2.682 + }
2.683 +
2.684 + no_more_push_2:
2.685 +
2.686 + _excess->set(n, excess);
2.687 +
2.688 + if (excess != 0) {
2.689 + if (new_level + 1 < _level->maxLevel()) {
2.690 + _level->liftActiveOn(level, new_level + 1);
2.691 + } else {
2.692 + _level->liftActiveToTop(level);
2.693 + }
2.694 + if (_level->emptyLevel(level)) {
2.695 + _level->liftToTop(level);
2.696 + }
2.697 + } else {
2.698 + _level->deactivate(n);
2.699 + }
2.700 +
2.701 + while (level >= 0 && _level->activeFree(level)) {
2.702 + --level;
2.703 + }
2.704 + if (level == -1) {
2.705 + n = _level->highestActive();
2.706 + level = _level->highestActiveLevel();
2.707 + } else {
2.708 + n = _level->activeOn(level);
2.709 + }
2.710 + --num;
2.711 + }
2.712 + }
2.713 + }
2.714 +
2.715 + /// \brief Starts the second phase of the preflow algorithm.
2.716 + ///
2.717 + /// The preflow algorithm consists of two phases, this method runs
2.718 + /// the second phase. After calling \ref init() and \ref
2.719 + /// startFirstPhase() and then \ref startSecondPhase(), \ref
2.720 + /// flowMap() return a maximum flow, \ref flowValue() returns the
2.721 + /// value of a maximum flow, \ref minCut() returns a minimum cut
2.722 + /// \pre The \ref init() and startFirstPhase() functions should be
2.723 + /// called before.
2.724 + void startSecondPhase() {
2.725 + _phase = false;
2.726 +
2.727 + typename Digraph::template NodeMap<bool> reached(_graph);
2.728 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.729 + reached.set(n, (*_level)[n] < _level->maxLevel());
2.730 + }
2.731 +
2.732 + _level->initStart();
2.733 + _level->initAddItem(_source);
2.734 +
2.735 + std::vector<Node> queue;
2.736 + queue.push_back(_source);
2.737 + reached.set(_source, true);
2.738 +
2.739 + while (!queue.empty()) {
2.740 + _level->initNewLevel();
2.741 + std::vector<Node> nqueue;
2.742 + for (int i = 0; i < int(queue.size()); ++i) {
2.743 + Node n = queue[i];
2.744 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2.745 + Node v = _graph.target(e);
2.746 + if (!reached[v] && _tolerance.positive((*_flow)[e])) {
2.747 + reached.set(v, true);
2.748 + _level->initAddItem(v);
2.749 + nqueue.push_back(v);
2.750 + }
2.751 + }
2.752 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.753 + Node u = _graph.source(e);
2.754 + if (!reached[u] &&
2.755 + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
2.756 + reached.set(u, true);
2.757 + _level->initAddItem(u);
2.758 + nqueue.push_back(u);
2.759 + }
2.760 + }
2.761 + }
2.762 + queue.swap(nqueue);
2.763 + }
2.764 + _level->initFinish();
2.765 +
2.766 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.767 + if (!reached[n]) {
2.768 + _level->dirtyTopButOne(n);
2.769 + } else if ((*_excess)[n] > 0 && _target != n) {
2.770 + _level->activate(n);
2.771 + }
2.772 + }
2.773 +
2.774 + Node n;
2.775 + while ((n = _level->highestActive()) != INVALID) {
2.776 + Value excess = (*_excess)[n];
2.777 + int level = _level->highestActiveLevel();
2.778 + int new_level = _level->maxLevel();
2.779 +
2.780 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
2.781 + Value rem = (*_capacity)[e] - (*_flow)[e];
2.782 + if (!_tolerance.positive(rem)) continue;
2.783 + Node v = _graph.target(e);
2.784 + if ((*_level)[v] < level) {
2.785 + if (!_level->active(v) && v != _source) {
2.786 + _level->activate(v);
2.787 + }
2.788 + if (!_tolerance.less(rem, excess)) {
2.789 + _flow->set(e, (*_flow)[e] + excess);
2.790 + _excess->set(v, (*_excess)[v] + excess);
2.791 + excess = 0;
2.792 + goto no_more_push;
2.793 + } else {
2.794 + excess -= rem;
2.795 + _excess->set(v, (*_excess)[v] + rem);
2.796 + _flow->set(e, (*_capacity)[e]);
2.797 + }
2.798 + } else if (new_level > (*_level)[v]) {
2.799 + new_level = (*_level)[v];
2.800 + }
2.801 + }
2.802 +
2.803 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
2.804 + Value rem = (*_flow)[e];
2.805 + if (!_tolerance.positive(rem)) continue;
2.806 + Node v = _graph.source(e);
2.807 + if ((*_level)[v] < level) {
2.808 + if (!_level->active(v) && v != _source) {
2.809 + _level->activate(v);
2.810 + }
2.811 + if (!_tolerance.less(rem, excess)) {
2.812 + _flow->set(e, (*_flow)[e] - excess);
2.813 + _excess->set(v, (*_excess)[v] + excess);
2.814 + excess = 0;
2.815 + goto no_more_push;
2.816 + } else {
2.817 + excess -= rem;
2.818 + _excess->set(v, (*_excess)[v] + rem);
2.819 + _flow->set(e, 0);
2.820 + }
2.821 + } else if (new_level > (*_level)[v]) {
2.822 + new_level = (*_level)[v];
2.823 + }
2.824 + }
2.825 +
2.826 + no_more_push:
2.827 +
2.828 + _excess->set(n, excess);
2.829 +
2.830 + if (excess != 0) {
2.831 + if (new_level + 1 < _level->maxLevel()) {
2.832 + _level->liftHighestActive(new_level + 1);
2.833 + } else {
2.834 + // Calculation error
2.835 + _level->liftHighestActiveToTop();
2.836 + }
2.837 + if (_level->emptyLevel(level)) {
2.838 + // Calculation error
2.839 + _level->liftToTop(level);
2.840 + }
2.841 + } else {
2.842 + _level->deactivate(n);
2.843 + }
2.844 +
2.845 + }
2.846 + }
2.847 +
2.848 + /// \brief Runs the preflow algorithm.
2.849 + ///
2.850 + /// Runs the preflow algorithm.
2.851 + /// \note pf.run() is just a shortcut of the following code.
2.852 + /// \code
2.853 + /// pf.init();
2.854 + /// pf.startFirstPhase();
2.855 + /// pf.startSecondPhase();
2.856 + /// \endcode
2.857 + void run() {
2.858 + init();
2.859 + startFirstPhase();
2.860 + startSecondPhase();
2.861 + }
2.862 +
2.863 + /// \brief Runs the preflow algorithm to compute the minimum cut.
2.864 + ///
2.865 + /// Runs the preflow algorithm to compute the minimum cut.
2.866 + /// \note pf.runMinCut() is just a shortcut of the following code.
2.867 + /// \code
2.868 + /// pf.init();
2.869 + /// pf.startFirstPhase();
2.870 + /// \endcode
2.871 + void runMinCut() {
2.872 + init();
2.873 + startFirstPhase();
2.874 + }
2.875 +
2.876 + /// @}
2.877 +
2.878 + /// \name Query Functions
2.879 + /// The result of the %Preflow algorithm can be obtained using these
2.880 + /// functions.\n
2.881 + /// Before the use of these functions,
2.882 + /// either run() or start() must be called.
2.883 +
2.884 + ///@{
2.885 +
2.886 + /// \brief Returns the value of the maximum flow.
2.887 + ///
2.888 + /// Returns the value of the maximum flow by returning the excess
2.889 + /// of the target node \c t. This value equals to the value of
2.890 + /// the maximum flow already after the first phase.
2.891 + Value flowValue() const {
2.892 + return (*_excess)[_target];
2.893 + }
2.894 +
2.895 + /// \brief Returns true when the node is on the source side of minimum cut.
2.896 + ///
2.897 + /// Returns true when the node is on the source side of minimum
2.898 + /// cut. This method can be called both after running \ref
2.899 + /// startFirstPhase() and \ref startSecondPhase().
2.900 + bool minCut(const Node& node) const {
2.901 + return ((*_level)[node] == _level->maxLevel()) == _phase;
2.902 + }
2.903 +
2.904 + /// \brief Returns a minimum value cut.
2.905 + ///
2.906 + /// Sets the \c cutMap to the characteristic vector of a minimum value
2.907 + /// cut. This method can be called both after running \ref
2.908 + /// startFirstPhase() and \ref startSecondPhase(). The result after second
2.909 + /// phase could be changed slightly if inexact computation is used.
2.910 + /// \pre The \c cutMap should be a bool-valued node-map.
2.911 + template <typename CutMap>
2.912 + void minCutMap(CutMap& cutMap) const {
2.913 + for (NodeIt n(_graph); n != INVALID; ++n) {
2.914 + cutMap.set(n, minCut(n));
2.915 + }
2.916 + }
2.917 +
2.918 + /// \brief Returns the flow on the arc.
2.919 + ///
2.920 + /// Sets the \c flowMap to the flow on the arcs. This method can
2.921 + /// be called after the second phase of algorithm.
2.922 + Value flow(const Arc& arc) const {
2.923 + return (*_flow)[arc];
2.924 + }
2.925 +
2.926 + /// @}
2.927 + };
2.928 +}
2.929 +
2.930 +#endif
3.1 --- a/test/Makefile.am Fri Nov 21 10:49:39 2008 +0000
3.2 +++ b/test/Makefile.am Fri Nov 21 14:11:29 2008 +0000
3.3 @@ -1,6 +1,7 @@
3.4 EXTRA_DIST += \
3.5 test/CMakeLists.txt \
3.6 - test/min_cost_flow_test.lgf
3.7 + test/min_cost_flow_test.lgf \
3.8 + test/preflow_graph.lgf
3.9
3.10 noinst_HEADERS += \
3.11 test/graph_test.h \
3.12 @@ -23,6 +24,7 @@
3.13 test/max_matching_test \
3.14 test/random_test \
3.15 test/path_test \
3.16 + test/preflow_test \
3.17 test/suurballe_test \
3.18 test/test_tools_fail \
3.19 test/test_tools_pass \
3.20 @@ -47,6 +49,7 @@
3.21 test_maps_test_SOURCES = test/maps_test.cc
3.22 test_max_matching_test_SOURCES = test/max_matching_test.cc
3.23 test_path_test_SOURCES = test/path_test.cc
3.24 +test_preflow_test_SOURCES = test/preflow_test.cc
3.25 test_suurballe_test_SOURCES = test/suurballe_test.cc
3.26 test_random_test_SOURCES = test/random_test.cc
3.27 test_test_tools_fail_SOURCES = test/test_tools_fail.cc
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/test/preflow_graph.lgf Fri Nov 21 14:11:29 2008 +0000
4.3 @@ -0,0 +1,35 @@
4.4 +@nodes
4.5 +label
4.6 +0
4.7 +1
4.8 +2
4.9 +3
4.10 +4
4.11 +5
4.12 +6
4.13 +7
4.14 +8
4.15 +9
4.16 +@edges
4.17 + label capacity
4.18 +0 1 0 20
4.19 +0 2 1 0
4.20 +1 1 2 3
4.21 +1 2 3 8
4.22 +1 3 4 8
4.23 +2 5 5 5
4.24 +3 2 6 5
4.25 +3 5 7 5
4.26 +3 6 8 5
4.27 +4 3 9 3
4.28 +5 7 10 3
4.29 +5 6 11 10
4.30 +5 8 12 10
4.31 +6 8 13 8
4.32 +8 9 14 20
4.33 +8 1 15 5
4.34 +9 5 16 5
4.35 +@attributes
4.36 +source 1
4.37 +target 8
4.38 +@end
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/test/preflow_test.cc Fri Nov 21 14:11:29 2008 +0000
5.3 @@ -0,0 +1,205 @@
5.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
5.5 + *
5.6 + * This file is a part of LEMON, a generic C++ optimization library.
5.7 + *
5.8 + * Copyright (C) 2003-2008
5.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
5.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
5.11 + *
5.12 + * Permission to use, modify and distribute this software is granted
5.13 + * provided that this copyright notice appears in all copies. For
5.14 + * precise terms see the accompanying LICENSE file.
5.15 + *
5.16 + * This software is provided "AS IS" with no warranty of any kind,
5.17 + * express or implied, and with no claim as to its suitability for any
5.18 + * purpose.
5.19 + *
5.20 + */
5.21 +
5.22 +#include <fstream>
5.23 +#include <string>
5.24 +
5.25 +#include "test_tools.h"
5.26 +#include <lemon/smart_graph.h>
5.27 +#include <lemon/preflow.h>
5.28 +#include <lemon/concepts/digraph.h>
5.29 +#include <lemon/concepts/maps.h>
5.30 +#include <lemon/lgf_reader.h>
5.31 +
5.32 +using namespace lemon;
5.33 +
5.34 +void checkPreflow()
5.35 +{
5.36 + typedef int VType;
5.37 + typedef concepts::Digraph Digraph;
5.38 +
5.39 + typedef Digraph::Node Node;
5.40 + typedef Digraph::Arc Arc;
5.41 + typedef concepts::ReadMap<Arc,VType> CapMap;
5.42 + typedef concepts::ReadWriteMap<Arc,VType> FlowMap;
5.43 + typedef concepts::WriteMap<Node,bool> CutMap;
5.44 +
5.45 + Digraph g;
5.46 + Node n;
5.47 + Arc e;
5.48 + CapMap cap;
5.49 + FlowMap flow;
5.50 + CutMap cut;
5.51 +
5.52 + Preflow<Digraph, CapMap>::DefFlowMap<FlowMap>::Create preflow_test(g,cap,n,n);
5.53 +
5.54 + preflow_test.capacityMap(cap);
5.55 + flow = preflow_test.flowMap();
5.56 + preflow_test.flowMap(flow);
5.57 + preflow_test.source(n);
5.58 + preflow_test.target(n);
5.59 +
5.60 + preflow_test.init();
5.61 + preflow_test.flowInit(cap);
5.62 + preflow_test.startFirstPhase();
5.63 + preflow_test.startSecondPhase();
5.64 + preflow_test.run();
5.65 + preflow_test.runMinCut();
5.66 +
5.67 + preflow_test.flowValue();
5.68 + preflow_test.minCut(n);
5.69 + preflow_test.minCutMap(cut);
5.70 + preflow_test.flow(e);
5.71 +
5.72 +}
5.73 +
5.74 +int cutValue (const SmartDigraph& g,
5.75 + const SmartDigraph::NodeMap<bool>& cut,
5.76 + const SmartDigraph::ArcMap<int>& cap) {
5.77 +
5.78 + int c=0;
5.79 + for(SmartDigraph::ArcIt e(g); e!=INVALID; ++e) {
5.80 + if (cut[g.source(e)] && !cut[g.target(e)]) c+=cap[e];
5.81 + }
5.82 + return c;
5.83 +}
5.84 +
5.85 +bool checkFlow(const SmartDigraph& g,
5.86 + const SmartDigraph::ArcMap<int>& flow,
5.87 + const SmartDigraph::ArcMap<int>& cap,
5.88 + SmartDigraph::Node s, SmartDigraph::Node t) {
5.89 +
5.90 + for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) {
5.91 + if (flow[e] < 0 || flow[e] > cap[e]) return false;
5.92 + }
5.93 +
5.94 + for (SmartDigraph::NodeIt n(g); n != INVALID; ++n) {
5.95 + if (n == s || n == t) continue;
5.96 + int sum = 0;
5.97 + for (SmartDigraph::OutArcIt e(g, n); e != INVALID; ++e) {
5.98 + sum += flow[e];
5.99 + }
5.100 + for (SmartDigraph::InArcIt e(g, n); e != INVALID; ++e) {
5.101 + sum -= flow[e];
5.102 + }
5.103 + if (sum != 0) return false;
5.104 + }
5.105 + return true;
5.106 +}
5.107 +
5.108 +int main() {
5.109 +
5.110 + typedef SmartDigraph Digraph;
5.111 +
5.112 + typedef Digraph::Node Node;
5.113 + typedef Digraph::NodeIt NodeIt;
5.114 + typedef Digraph::ArcIt ArcIt;
5.115 + typedef Digraph::ArcMap<int> CapMap;
5.116 + typedef Digraph::ArcMap<int> FlowMap;
5.117 + typedef Digraph::NodeMap<bool> CutMap;
5.118 +
5.119 + typedef Preflow<Digraph, CapMap> PType;
5.120 +
5.121 + std::string f_name;
5.122 + if( getenv("srcdir") )
5.123 + f_name = std::string(getenv("srcdir"));
5.124 + else f_name = ".";
5.125 + f_name += "/test/preflow_graph.lgf";
5.126 +
5.127 + std::ifstream file(f_name.c_str());
5.128 +
5.129 + check(file, "Input file '" << f_name << "' not found.");
5.130 +
5.131 + Digraph g;
5.132 + Node s, t;
5.133 + CapMap cap(g);
5.134 + DigraphReader<Digraph>(g,file).
5.135 + arcMap("capacity", cap).
5.136 + node("source",s).
5.137 + node("target",t).
5.138 + run();
5.139 +
5.140 + PType preflow_test(g, cap, s, t);
5.141 + preflow_test.run();
5.142 +
5.143 + check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
5.144 + "The flow is not feasible.");
5.145 +
5.146 + CutMap min_cut(g);
5.147 + preflow_test.minCutMap(min_cut);
5.148 + int min_cut_value=cutValue(g,min_cut,cap);
5.149 +
5.150 + check(preflow_test.flowValue() == min_cut_value,
5.151 + "The max flow value is not equal to the three min cut values.");
5.152 +
5.153 + FlowMap flow(g);
5.154 + for(ArcIt e(g); e!=INVALID; ++e) flow[e] = preflow_test.flowMap()[e];
5.155 +
5.156 + int flow_value=preflow_test.flowValue();
5.157 +
5.158 + for(ArcIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e];
5.159 + preflow_test.flowInit(flow);
5.160 + preflow_test.startFirstPhase();
5.161 +
5.162 + CutMap min_cut1(g);
5.163 + preflow_test.minCutMap(min_cut1);
5.164 + min_cut_value=cutValue(g,min_cut1,cap);
5.165 +
5.166 + check(preflow_test.flowValue() == min_cut_value &&
5.167 + min_cut_value == 2*flow_value,
5.168 + "The max flow value or the min cut value is wrong.");
5.169 +
5.170 + preflow_test.startSecondPhase();
5.171 +
5.172 + check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
5.173 + "The flow is not feasible.");
5.174 +
5.175 + CutMap min_cut2(g);
5.176 + preflow_test.minCutMap(min_cut2);
5.177 + min_cut_value=cutValue(g,min_cut2,cap);
5.178 +
5.179 + check(preflow_test.flowValue() == min_cut_value &&
5.180 + min_cut_value == 2*flow_value,
5.181 + "The max flow value or the three min cut values were not doubled");
5.182 +
5.183 +
5.184 + preflow_test.flowMap(flow);
5.185 +
5.186 + NodeIt tmp1(g,s);
5.187 + ++tmp1;
5.188 + if ( tmp1 != INVALID ) s=tmp1;
5.189 +
5.190 + NodeIt tmp2(g,t);
5.191 + ++tmp2;
5.192 + if ( tmp2 != INVALID ) t=tmp2;
5.193 +
5.194 + preflow_test.source(s);
5.195 + preflow_test.target(t);
5.196 +
5.197 + preflow_test.run();
5.198 +
5.199 + CutMap min_cut3(g);
5.200 + preflow_test.minCutMap(min_cut3);
5.201 + min_cut_value=cutValue(g,min_cut3,cap);
5.202 +
5.203 +
5.204 + check(preflow_test.flowValue() == min_cut_value,
5.205 + "The max flow value or the three min cut values are incorrect.");
5.206 +
5.207 + return 0;
5.208 +}