Port preflow push max flow alg. from svn -r3516 (#176)
authorAlpar Juttner <alpar@cs.elte.hu>
Fri, 21 Nov 2008 14:11:29 +0000
changeset 389660db48f324f
parent 384 fc6954b4fce4
child 390 53c5277ba294
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
lemon/Makefile.am
lemon/preflow.h
test/Makefile.am
test/preflow_graph.lgf
test/preflow_test.cc
     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 = &map;
   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 = &map;
   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 +}