# HG changeset patch # User Alpar Juttner # Date 2008-11-21 15:11:29 # Node ID 660db48f324f9becef9aee3407289d4ce146a79e # Parent fc6954b4fce49fa09a571eb41b2d16ab7c463569 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 diff --git a/lemon/Makefile.am b/lemon/Makefile.am --- a/lemon/Makefile.am +++ b/lemon/Makefile.am @@ -42,6 +42,7 @@ lemon/max_matching.h \ lemon/nauty_reader.h \ lemon/path.h \ + lemon/preflow.h \ lemon/random.h \ lemon/smart_graph.h \ lemon/suurballe.h \ diff --git a/lemon/preflow.h b/lemon/preflow.h new file mode 100644 --- /dev/null +++ b/lemon/preflow.h @@ -0,0 +1,927 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_PREFLOW_H +#define LEMON_PREFLOW_H + +#include +#include +#include + +/// \file +/// \ingroup max_flow +/// \brief Implementation of the preflow algorithm. + +namespace lemon { + + /// \brief Default traits class of Preflow class. + /// + /// Default traits class of Preflow class. + /// \param _Graph Digraph type. + /// \param _CapacityMap Type of capacity map. + template + struct PreflowDefaultTraits { + + /// \brief The digraph type the algorithm runs on. + typedef _Graph Digraph; + + /// \brief The type of the map that stores the arc capacities. + /// + /// The type of the map that stores the arc capacities. + /// It must meet the \ref concepts::ReadMap "ReadMap" concept. + typedef _CapacityMap CapacityMap; + + /// \brief The type of the length of the arcs. + typedef typename CapacityMap::Value Value; + + /// \brief The map type that stores the flow values. + /// + /// The map type that stores the flow values. + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept. + typedef typename Digraph::template ArcMap FlowMap; + + /// \brief Instantiates a FlowMap. + /// + /// This function instantiates a \ref FlowMap. + /// \param digraph The digraph, to which we would like to define + /// the flow map. + static FlowMap* createFlowMap(const Digraph& digraph) { + return new FlowMap(digraph); + } + + /// \brief The eleavator type used by Preflow algorithm. + /// + /// The elevator type used by Preflow algorithm. + /// + /// \sa Elevator + /// \sa LinkedElevator + typedef LinkedElevator Elevator; + + /// \brief Instantiates an Elevator. + /// + /// This function instantiates a \ref Elevator. + /// \param digraph The digraph, to which we would like to define + /// the elevator. + /// \param max_level The maximum level of the elevator. + static Elevator* createElevator(const Digraph& digraph, int max_level) { + return new Elevator(digraph, max_level); + } + + /// \brief The tolerance used by the algorithm + /// + /// The tolerance used by the algorithm to handle inexact computation. + typedef lemon::Tolerance Tolerance; + + }; + + + /// \ingroup max_flow + /// + /// \brief %Preflow algorithms class. + /// + /// This class provides an implementation of the Goldberg's \e + /// preflow \e algorithm producing a flow of maximum value in a + /// digraph. The preflow algorithms are the fastest known max + /// flow algorithms. The current implementation use a mixture of the + /// \e "highest label" and the \e "bound decrease" heuristics. + /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$. + /// + /// The algorithm consists from two phases. After the first phase + /// the maximal flow value and the minimum cut can be obtained. The + /// second phase constructs the feasible maximum flow on each arc. + /// + /// \param _Graph The digraph type the algorithm runs on. + /// \param _CapacityMap The flow map type. + /// \param _Traits Traits class to set various data types used by + /// the algorithm. The default traits class is \ref + /// PreflowDefaultTraits. See \ref PreflowDefaultTraits for the + /// documentation of a %Preflow traits class. + /// + ///\author Jacint Szabo and Balazs Dezso +#ifdef DOXYGEN + template +#else + template , + typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> > +#endif + class Preflow { + public: + + typedef _Traits Traits; + typedef typename Traits::Digraph Digraph; + typedef typename Traits::CapacityMap CapacityMap; + typedef typename Traits::Value Value; + + typedef typename Traits::FlowMap FlowMap; + typedef typename Traits::Elevator Elevator; + typedef typename Traits::Tolerance Tolerance; + + /// \brief \ref Exception for uninitialized parameters. + /// + /// This error represents problems in the initialization + /// of the parameters of the algorithms. + class UninitializedParameter : public lemon::Exception { + public: + virtual const char* what() const throw() { + return "lemon::Preflow::UninitializedParameter"; + } + }; + + private: + + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); + + const Digraph& _graph; + const CapacityMap* _capacity; + + int _node_num; + + Node _source, _target; + + FlowMap* _flow; + bool _local_flow; + + Elevator* _level; + bool _local_level; + + typedef typename Digraph::template NodeMap ExcessMap; + ExcessMap* _excess; + + Tolerance _tolerance; + + bool _phase; + + + void createStructures() { + _node_num = countNodes(_graph); + + if (!_flow) { + _flow = Traits::createFlowMap(_graph); + _local_flow = true; + } + if (!_level) { + _level = Traits::createElevator(_graph, _node_num); + _local_level = true; + } + if (!_excess) { + _excess = new ExcessMap(_graph); + } + } + + void destroyStructures() { + if (_local_flow) { + delete _flow; + } + if (_local_level) { + delete _level; + } + if (_excess) { + delete _excess; + } + } + + public: + + typedef Preflow Create; + + ///\name Named template parameters + + ///@{ + + template + struct DefFlowMapTraits : public Traits { + typedef _FlowMap FlowMap; + static FlowMap *createFlowMap(const Digraph&) { + throw UninitializedParameter(); + } + }; + + /// \brief \ref named-templ-param "Named parameter" for setting + /// FlowMap type + /// + /// \ref named-templ-param "Named parameter" for setting FlowMap + /// type + template + struct DefFlowMap + : public Preflow > { + typedef Preflow > Create; + }; + + template + struct DefElevatorTraits : public Traits { + typedef _Elevator Elevator; + static Elevator *createElevator(const Digraph&, int) { + throw UninitializedParameter(); + } + }; + + /// \brief \ref named-templ-param "Named parameter" for setting + /// Elevator type + /// + /// \ref named-templ-param "Named parameter" for setting Elevator + /// type + template + struct DefElevator + : public Preflow > { + typedef Preflow > Create; + }; + + template + struct DefStandardElevatorTraits : public Traits { + typedef _Elevator Elevator; + static Elevator *createElevator(const Digraph& digraph, int max_level) { + return new Elevator(digraph, max_level); + } + }; + + /// \brief \ref named-templ-param "Named parameter" for setting + /// Elevator type + /// + /// \ref named-templ-param "Named parameter" for setting Elevator + /// type. The Elevator should be standard constructor interface, ie. + /// the digraph and the maximum level should be passed to it. + template + struct DefStandardElevator + : public Preflow > { + typedef Preflow > Create; + }; + + /// @} + + protected: + + Preflow() {} + + public: + + + /// \brief The constructor of the class. + /// + /// The constructor of the class. + /// \param digraph The digraph the algorithm runs on. + /// \param capacity The capacity of the arcs. + /// \param source The source node. + /// \param target The target node. + Preflow(const Digraph& digraph, const CapacityMap& capacity, + Node source, Node target) + : _graph(digraph), _capacity(&capacity), + _node_num(0), _source(source), _target(target), + _flow(0), _local_flow(false), + _level(0), _local_level(false), + _excess(0), _tolerance(), _phase() {} + + /// \brief Destrcutor. + /// + /// Destructor. + ~Preflow() { + destroyStructures(); + } + + /// \brief Sets the capacity map. + /// + /// Sets the capacity map. + /// \return \c (*this) + Preflow& capacityMap(const CapacityMap& map) { + _capacity = ↦ + return *this; + } + + /// \brief Sets the flow map. + /// + /// Sets the flow map. + /// \return \c (*this) + Preflow& flowMap(FlowMap& map) { + if (_local_flow) { + delete _flow; + _local_flow = false; + } + _flow = ↦ + return *this; + } + + /// \brief Returns the flow map. + /// + /// \return The flow map. + const FlowMap& flowMap() { + return *_flow; + } + + /// \brief Sets the elevator. + /// + /// Sets the elevator. + /// \return \c (*this) + Preflow& elevator(Elevator& elevator) { + if (_local_level) { + delete _level; + _local_level = false; + } + _level = &elevator; + return *this; + } + + /// \brief Returns the elevator. + /// + /// \return The elevator. + const Elevator& elevator() { + return *_level; + } + + /// \brief Sets the source node. + /// + /// Sets the source node. + /// \return \c (*this) + Preflow& source(const Node& node) { + _source = node; + return *this; + } + + /// \brief Sets the target node. + /// + /// Sets the target node. + /// \return \c (*this) + Preflow& target(const Node& node) { + _target = node; + return *this; + } + + /// \brief Sets the tolerance used by algorithm. + /// + /// Sets the tolerance used by algorithm. + Preflow& tolerance(const Tolerance& tolerance) const { + _tolerance = tolerance; + return *this; + } + + /// \brief Returns the tolerance used by algorithm. + /// + /// Returns the tolerance used by algorithm. + const Tolerance& tolerance() const { + return tolerance; + } + + /// \name Execution control The simplest way to execute the + /// algorithm is to use one of the member functions called \c + /// run(). + /// \n + /// If you need more control on initial solution or + /// execution then you have to call one \ref init() function and then + /// the startFirstPhase() and if you need the startSecondPhase(). + + ///@{ + + /// \brief Initializes the internal data structures. + /// + /// Initializes the internal data structures. + /// + void init() { + createStructures(); + + _phase = true; + for (NodeIt n(_graph); n != INVALID; ++n) { + _excess->set(n, 0); + } + + for (ArcIt e(_graph); e != INVALID; ++e) { + _flow->set(e, 0); + } + + typename Digraph::template NodeMap reached(_graph, false); + + _level->initStart(); + _level->initAddItem(_target); + + std::vector queue; + reached.set(_source, true); + + queue.push_back(_target); + reached.set(_target, true); + while (!queue.empty()) { + _level->initNewLevel(); + std::vector nqueue; + for (int i = 0; i < int(queue.size()); ++i) { + Node n = queue[i]; + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node u = _graph.source(e); + if (!reached[u] && _tolerance.positive((*_capacity)[e])) { + reached.set(u, true); + _level->initAddItem(u); + nqueue.push_back(u); + } + } + } + queue.swap(nqueue); + } + _level->initFinish(); + + for (OutArcIt e(_graph, _source); e != INVALID; ++e) { + if (_tolerance.positive((*_capacity)[e])) { + Node u = _graph.target(e); + if ((*_level)[u] == _level->maxLevel()) continue; + _flow->set(e, (*_capacity)[e]); + _excess->set(u, (*_excess)[u] + (*_capacity)[e]); + if (u != _target && !_level->active(u)) { + _level->activate(u); + } + } + } + } + + /// \brief Initializes the internal data structures. + /// + /// Initializes the internal data structures and sets the initial + /// flow to the given \c flowMap. The \c flowMap should contain a + /// flow or at least a preflow, ie. in each node excluding the + /// target the incoming flow should greater or equal to the + /// outgoing flow. + /// \return %False when the given \c flowMap is not a preflow. + template + bool flowInit(const FlowMap& flowMap) { + createStructures(); + + for (ArcIt e(_graph); e != INVALID; ++e) { + _flow->set(e, flowMap[e]); + } + + for (NodeIt n(_graph); n != INVALID; ++n) { + Value excess = 0; + for (InArcIt e(_graph, n); e != INVALID; ++e) { + excess += (*_flow)[e]; + } + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + excess -= (*_flow)[e]; + } + if (excess < 0 && n != _source) return false; + _excess->set(n, excess); + } + + typename Digraph::template NodeMap reached(_graph, false); + + _level->initStart(); + _level->initAddItem(_target); + + std::vector queue; + reached.set(_source, true); + + queue.push_back(_target); + reached.set(_target, true); + while (!queue.empty()) { + _level->initNewLevel(); + std::vector nqueue; + for (int i = 0; i < int(queue.size()); ++i) { + Node n = queue[i]; + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node u = _graph.source(e); + if (!reached[u] && + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) { + reached.set(u, true); + _level->initAddItem(u); + nqueue.push_back(u); + } + } + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.target(e); + if (!reached[v] && _tolerance.positive((*_flow)[e])) { + reached.set(v, true); + _level->initAddItem(v); + nqueue.push_back(v); + } + } + } + queue.swap(nqueue); + } + _level->initFinish(); + + for (OutArcIt e(_graph, _source); e != INVALID; ++e) { + Value rem = (*_capacity)[e] - (*_flow)[e]; + if (_tolerance.positive(rem)) { + Node u = _graph.target(e); + if ((*_level)[u] == _level->maxLevel()) continue; + _flow->set(e, (*_capacity)[e]); + _excess->set(u, (*_excess)[u] + rem); + if (u != _target && !_level->active(u)) { + _level->activate(u); + } + } + } + for (InArcIt e(_graph, _source); e != INVALID; ++e) { + Value rem = (*_flow)[e]; + if (_tolerance.positive(rem)) { + Node v = _graph.source(e); + if ((*_level)[v] == _level->maxLevel()) continue; + _flow->set(e, 0); + _excess->set(v, (*_excess)[v] + rem); + if (v != _target && !_level->active(v)) { + _level->activate(v); + } + } + } + return true; + } + + /// \brief Starts the first phase of the preflow algorithm. + /// + /// The preflow algorithm consists of two phases, this method runs + /// the first phase. After the first phase the maximum flow value + /// and a minimum value cut can already be computed, although a + /// maximum flow is not yet obtained. So after calling this method + /// \ref flowValue() returns the value of a maximum flow and \ref + /// minCut() returns a minimum cut. + /// \pre One of the \ref init() functions should be called. + void startFirstPhase() { + _phase = true; + + Node n = _level->highestActive(); + int level = _level->highestActiveLevel(); + while (n != INVALID) { + int num = _node_num; + + while (num > 0 && n != INVALID) { + Value excess = (*_excess)[n]; + int new_level = _level->maxLevel(); + + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + Value rem = (*_capacity)[e] - (*_flow)[e]; + if (!_tolerance.positive(rem)) continue; + Node v = _graph.target(e); + if ((*_level)[v] < level) { + if (!_level->active(v) && v != _target) { + _level->activate(v); + } + if (!_tolerance.less(rem, excess)) { + _flow->set(e, (*_flow)[e] + excess); + _excess->set(v, (*_excess)[v] + excess); + excess = 0; + goto no_more_push_1; + } else { + excess -= rem; + _excess->set(v, (*_excess)[v] + rem); + _flow->set(e, (*_capacity)[e]); + } + } else if (new_level > (*_level)[v]) { + new_level = (*_level)[v]; + } + } + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Value rem = (*_flow)[e]; + if (!_tolerance.positive(rem)) continue; + Node v = _graph.source(e); + if ((*_level)[v] < level) { + if (!_level->active(v) && v != _target) { + _level->activate(v); + } + if (!_tolerance.less(rem, excess)) { + _flow->set(e, (*_flow)[e] - excess); + _excess->set(v, (*_excess)[v] + excess); + excess = 0; + goto no_more_push_1; + } else { + excess -= rem; + _excess->set(v, (*_excess)[v] + rem); + _flow->set(e, 0); + } + } else if (new_level > (*_level)[v]) { + new_level = (*_level)[v]; + } + } + + no_more_push_1: + + _excess->set(n, excess); + + if (excess != 0) { + if (new_level + 1 < _level->maxLevel()) { + _level->liftHighestActive(new_level + 1); + } else { + _level->liftHighestActiveToTop(); + } + if (_level->emptyLevel(level)) { + _level->liftToTop(level); + } + } else { + _level->deactivate(n); + } + + n = _level->highestActive(); + level = _level->highestActiveLevel(); + --num; + } + + num = _node_num * 20; + while (num > 0 && n != INVALID) { + Value excess = (*_excess)[n]; + int new_level = _level->maxLevel(); + + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + Value rem = (*_capacity)[e] - (*_flow)[e]; + if (!_tolerance.positive(rem)) continue; + Node v = _graph.target(e); + if ((*_level)[v] < level) { + if (!_level->active(v) && v != _target) { + _level->activate(v); + } + if (!_tolerance.less(rem, excess)) { + _flow->set(e, (*_flow)[e] + excess); + _excess->set(v, (*_excess)[v] + excess); + excess = 0; + goto no_more_push_2; + } else { + excess -= rem; + _excess->set(v, (*_excess)[v] + rem); + _flow->set(e, (*_capacity)[e]); + } + } else if (new_level > (*_level)[v]) { + new_level = (*_level)[v]; + } + } + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Value rem = (*_flow)[e]; + if (!_tolerance.positive(rem)) continue; + Node v = _graph.source(e); + if ((*_level)[v] < level) { + if (!_level->active(v) && v != _target) { + _level->activate(v); + } + if (!_tolerance.less(rem, excess)) { + _flow->set(e, (*_flow)[e] - excess); + _excess->set(v, (*_excess)[v] + excess); + excess = 0; + goto no_more_push_2; + } else { + excess -= rem; + _excess->set(v, (*_excess)[v] + rem); + _flow->set(e, 0); + } + } else if (new_level > (*_level)[v]) { + new_level = (*_level)[v]; + } + } + + no_more_push_2: + + _excess->set(n, excess); + + if (excess != 0) { + if (new_level + 1 < _level->maxLevel()) { + _level->liftActiveOn(level, new_level + 1); + } else { + _level->liftActiveToTop(level); + } + if (_level->emptyLevel(level)) { + _level->liftToTop(level); + } + } else { + _level->deactivate(n); + } + + while (level >= 0 && _level->activeFree(level)) { + --level; + } + if (level == -1) { + n = _level->highestActive(); + level = _level->highestActiveLevel(); + } else { + n = _level->activeOn(level); + } + --num; + } + } + } + + /// \brief Starts the second phase of the preflow algorithm. + /// + /// The preflow algorithm consists of two phases, this method runs + /// the second phase. After calling \ref init() and \ref + /// startFirstPhase() and then \ref startSecondPhase(), \ref + /// flowMap() return a maximum flow, \ref flowValue() returns the + /// value of a maximum flow, \ref minCut() returns a minimum cut + /// \pre The \ref init() and startFirstPhase() functions should be + /// called before. + void startSecondPhase() { + _phase = false; + + typename Digraph::template NodeMap reached(_graph); + for (NodeIt n(_graph); n != INVALID; ++n) { + reached.set(n, (*_level)[n] < _level->maxLevel()); + } + + _level->initStart(); + _level->initAddItem(_source); + + std::vector queue; + queue.push_back(_source); + reached.set(_source, true); + + while (!queue.empty()) { + _level->initNewLevel(); + std::vector nqueue; + for (int i = 0; i < int(queue.size()); ++i) { + Node n = queue[i]; + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + Node v = _graph.target(e); + if (!reached[v] && _tolerance.positive((*_flow)[e])) { + reached.set(v, true); + _level->initAddItem(v); + nqueue.push_back(v); + } + } + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Node u = _graph.source(e); + if (!reached[u] && + _tolerance.positive((*_capacity)[e] - (*_flow)[e])) { + reached.set(u, true); + _level->initAddItem(u); + nqueue.push_back(u); + } + } + } + queue.swap(nqueue); + } + _level->initFinish(); + + for (NodeIt n(_graph); n != INVALID; ++n) { + if (!reached[n]) { + _level->dirtyTopButOne(n); + } else if ((*_excess)[n] > 0 && _target != n) { + _level->activate(n); + } + } + + Node n; + while ((n = _level->highestActive()) != INVALID) { + Value excess = (*_excess)[n]; + int level = _level->highestActiveLevel(); + int new_level = _level->maxLevel(); + + for (OutArcIt e(_graph, n); e != INVALID; ++e) { + Value rem = (*_capacity)[e] - (*_flow)[e]; + if (!_tolerance.positive(rem)) continue; + Node v = _graph.target(e); + if ((*_level)[v] < level) { + if (!_level->active(v) && v != _source) { + _level->activate(v); + } + if (!_tolerance.less(rem, excess)) { + _flow->set(e, (*_flow)[e] + excess); + _excess->set(v, (*_excess)[v] + excess); + excess = 0; + goto no_more_push; + } else { + excess -= rem; + _excess->set(v, (*_excess)[v] + rem); + _flow->set(e, (*_capacity)[e]); + } + } else if (new_level > (*_level)[v]) { + new_level = (*_level)[v]; + } + } + + for (InArcIt e(_graph, n); e != INVALID; ++e) { + Value rem = (*_flow)[e]; + if (!_tolerance.positive(rem)) continue; + Node v = _graph.source(e); + if ((*_level)[v] < level) { + if (!_level->active(v) && v != _source) { + _level->activate(v); + } + if (!_tolerance.less(rem, excess)) { + _flow->set(e, (*_flow)[e] - excess); + _excess->set(v, (*_excess)[v] + excess); + excess = 0; + goto no_more_push; + } else { + excess -= rem; + _excess->set(v, (*_excess)[v] + rem); + _flow->set(e, 0); + } + } else if (new_level > (*_level)[v]) { + new_level = (*_level)[v]; + } + } + + no_more_push: + + _excess->set(n, excess); + + if (excess != 0) { + if (new_level + 1 < _level->maxLevel()) { + _level->liftHighestActive(new_level + 1); + } else { + // Calculation error + _level->liftHighestActiveToTop(); + } + if (_level->emptyLevel(level)) { + // Calculation error + _level->liftToTop(level); + } + } else { + _level->deactivate(n); + } + + } + } + + /// \brief Runs the preflow algorithm. + /// + /// Runs the preflow algorithm. + /// \note pf.run() is just a shortcut of the following code. + /// \code + /// pf.init(); + /// pf.startFirstPhase(); + /// pf.startSecondPhase(); + /// \endcode + void run() { + init(); + startFirstPhase(); + startSecondPhase(); + } + + /// \brief Runs the preflow algorithm to compute the minimum cut. + /// + /// Runs the preflow algorithm to compute the minimum cut. + /// \note pf.runMinCut() is just a shortcut of the following code. + /// \code + /// pf.init(); + /// pf.startFirstPhase(); + /// \endcode + void runMinCut() { + init(); + startFirstPhase(); + } + + /// @} + + /// \name Query Functions + /// The result of the %Preflow algorithm can be obtained using these + /// functions.\n + /// Before the use of these functions, + /// either run() or start() must be called. + + ///@{ + + /// \brief Returns the value of the maximum flow. + /// + /// Returns the value of the maximum flow by returning the excess + /// of the target node \c t. This value equals to the value of + /// the maximum flow already after the first phase. + Value flowValue() const { + return (*_excess)[_target]; + } + + /// \brief Returns true when the node is on the source side of minimum cut. + /// + /// Returns true when the node is on the source side of minimum + /// cut. This method can be called both after running \ref + /// startFirstPhase() and \ref startSecondPhase(). + bool minCut(const Node& node) const { + return ((*_level)[node] == _level->maxLevel()) == _phase; + } + + /// \brief Returns a minimum value cut. + /// + /// Sets the \c cutMap to the characteristic vector of a minimum value + /// cut. This method can be called both after running \ref + /// startFirstPhase() and \ref startSecondPhase(). The result after second + /// phase could be changed slightly if inexact computation is used. + /// \pre The \c cutMap should be a bool-valued node-map. + template + void minCutMap(CutMap& cutMap) const { + for (NodeIt n(_graph); n != INVALID; ++n) { + cutMap.set(n, minCut(n)); + } + } + + /// \brief Returns the flow on the arc. + /// + /// Sets the \c flowMap to the flow on the arcs. This method can + /// be called after the second phase of algorithm. + Value flow(const Arc& arc) const { + return (*_flow)[arc]; + } + + /// @} + }; +} + +#endif diff --git a/test/Makefile.am b/test/Makefile.am --- a/test/Makefile.am +++ b/test/Makefile.am @@ -1,6 +1,7 @@ EXTRA_DIST += \ test/CMakeLists.txt \ - test/min_cost_flow_test.lgf + test/min_cost_flow_test.lgf \ + test/preflow_graph.lgf noinst_HEADERS += \ test/graph_test.h \ @@ -23,6 +24,7 @@ test/max_matching_test \ test/random_test \ test/path_test \ + test/preflow_test \ test/suurballe_test \ test/test_tools_fail \ test/test_tools_pass \ @@ -47,6 +49,7 @@ test_maps_test_SOURCES = test/maps_test.cc test_max_matching_test_SOURCES = test/max_matching_test.cc test_path_test_SOURCES = test/path_test.cc +test_preflow_test_SOURCES = test/preflow_test.cc test_suurballe_test_SOURCES = test/suurballe_test.cc test_random_test_SOURCES = test/random_test.cc test_test_tools_fail_SOURCES = test/test_tools_fail.cc diff --git a/test/preflow_graph.lgf b/test/preflow_graph.lgf new file mode 100644 --- /dev/null +++ b/test/preflow_graph.lgf @@ -0,0 +1,35 @@ +@nodes +label +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +@edges + label capacity +0 1 0 20 +0 2 1 0 +1 1 2 3 +1 2 3 8 +1 3 4 8 +2 5 5 5 +3 2 6 5 +3 5 7 5 +3 6 8 5 +4 3 9 3 +5 7 10 3 +5 6 11 10 +5 8 12 10 +6 8 13 8 +8 9 14 20 +8 1 15 5 +9 5 16 5 +@attributes +source 1 +target 8 +@end diff --git a/test/preflow_test.cc b/test/preflow_test.cc new file mode 100644 --- /dev/null +++ b/test/preflow_test.cc @@ -0,0 +1,205 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#include +#include + +#include "test_tools.h" +#include +#include +#include +#include +#include + +using namespace lemon; + +void checkPreflow() +{ + typedef int VType; + typedef concepts::Digraph Digraph; + + typedef Digraph::Node Node; + typedef Digraph::Arc Arc; + typedef concepts::ReadMap CapMap; + typedef concepts::ReadWriteMap FlowMap; + typedef concepts::WriteMap CutMap; + + Digraph g; + Node n; + Arc e; + CapMap cap; + FlowMap flow; + CutMap cut; + + Preflow::DefFlowMap::Create preflow_test(g,cap,n,n); + + preflow_test.capacityMap(cap); + flow = preflow_test.flowMap(); + preflow_test.flowMap(flow); + preflow_test.source(n); + preflow_test.target(n); + + preflow_test.init(); + preflow_test.flowInit(cap); + preflow_test.startFirstPhase(); + preflow_test.startSecondPhase(); + preflow_test.run(); + preflow_test.runMinCut(); + + preflow_test.flowValue(); + preflow_test.minCut(n); + preflow_test.minCutMap(cut); + preflow_test.flow(e); + +} + +int cutValue (const SmartDigraph& g, + const SmartDigraph::NodeMap& cut, + const SmartDigraph::ArcMap& cap) { + + int c=0; + for(SmartDigraph::ArcIt e(g); e!=INVALID; ++e) { + if (cut[g.source(e)] && !cut[g.target(e)]) c+=cap[e]; + } + return c; +} + +bool checkFlow(const SmartDigraph& g, + const SmartDigraph::ArcMap& flow, + const SmartDigraph::ArcMap& cap, + SmartDigraph::Node s, SmartDigraph::Node t) { + + for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) { + if (flow[e] < 0 || flow[e] > cap[e]) return false; + } + + for (SmartDigraph::NodeIt n(g); n != INVALID; ++n) { + if (n == s || n == t) continue; + int sum = 0; + for (SmartDigraph::OutArcIt e(g, n); e != INVALID; ++e) { + sum += flow[e]; + } + for (SmartDigraph::InArcIt e(g, n); e != INVALID; ++e) { + sum -= flow[e]; + } + if (sum != 0) return false; + } + return true; +} + +int main() { + + typedef SmartDigraph Digraph; + + typedef Digraph::Node Node; + typedef Digraph::NodeIt NodeIt; + typedef Digraph::ArcIt ArcIt; + typedef Digraph::ArcMap CapMap; + typedef Digraph::ArcMap FlowMap; + typedef Digraph::NodeMap CutMap; + + typedef Preflow PType; + + std::string f_name; + if( getenv("srcdir") ) + f_name = std::string(getenv("srcdir")); + else f_name = "."; + f_name += "/test/preflow_graph.lgf"; + + std::ifstream file(f_name.c_str()); + + check(file, "Input file '" << f_name << "' not found."); + + Digraph g; + Node s, t; + CapMap cap(g); + DigraphReader(g,file). + arcMap("capacity", cap). + node("source",s). + node("target",t). + run(); + + PType preflow_test(g, cap, s, t); + preflow_test.run(); + + check(checkFlow(g, preflow_test.flowMap(), cap, s, t), + "The flow is not feasible."); + + CutMap min_cut(g); + preflow_test.minCutMap(min_cut); + int min_cut_value=cutValue(g,min_cut,cap); + + check(preflow_test.flowValue() == min_cut_value, + "The max flow value is not equal to the three min cut values."); + + FlowMap flow(g); + for(ArcIt e(g); e!=INVALID; ++e) flow[e] = preflow_test.flowMap()[e]; + + int flow_value=preflow_test.flowValue(); + + for(ArcIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e]; + preflow_test.flowInit(flow); + preflow_test.startFirstPhase(); + + CutMap min_cut1(g); + preflow_test.minCutMap(min_cut1); + min_cut_value=cutValue(g,min_cut1,cap); + + check(preflow_test.flowValue() == min_cut_value && + min_cut_value == 2*flow_value, + "The max flow value or the min cut value is wrong."); + + preflow_test.startSecondPhase(); + + check(checkFlow(g, preflow_test.flowMap(), cap, s, t), + "The flow is not feasible."); + + CutMap min_cut2(g); + preflow_test.minCutMap(min_cut2); + min_cut_value=cutValue(g,min_cut2,cap); + + check(preflow_test.flowValue() == min_cut_value && + min_cut_value == 2*flow_value, + "The max flow value or the three min cut values were not doubled"); + + + preflow_test.flowMap(flow); + + NodeIt tmp1(g,s); + ++tmp1; + if ( tmp1 != INVALID ) s=tmp1; + + NodeIt tmp2(g,t); + ++tmp2; + if ( tmp2 != INVALID ) t=tmp2; + + preflow_test.source(s); + preflow_test.target(t); + + preflow_test.run(); + + CutMap min_cut3(g); + preflow_test.minCutMap(min_cut3); + min_cut_value=cutValue(g,min_cut3,cap); + + + check(preflow_test.flowValue() == min_cut_value, + "The max flow value or the three min cut values are incorrect."); + + return 0; +}