1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/edmonds_karp.h Wed Oct 17 19:14:07 2018 +0200
1.3 @@ -0,0 +1,556 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2013
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_EDMONDS_KARP_H
1.23 +#define LEMON_EDMONDS_KARP_H
1.24 +
1.25 +/// \file
1.26 +/// \ingroup max_flow
1.27 +/// \brief Implementation of the Edmonds-Karp algorithm.
1.28 +
1.29 +#include <lemon/tolerance.h>
1.30 +#include <vector>
1.31 +
1.32 +namespace lemon {
1.33 +
1.34 + /// \brief Default traits class of EdmondsKarp class.
1.35 + ///
1.36 + /// Default traits class of EdmondsKarp class.
1.37 + /// \param GR Digraph type.
1.38 + /// \param CAP Type of capacity map.
1.39 + template <typename GR, typename CAP>
1.40 + struct EdmondsKarpDefaultTraits {
1.41 +
1.42 + /// \brief The digraph type the algorithm runs on.
1.43 + typedef GR Digraph;
1.44 +
1.45 + /// \brief The type of the map that stores the arc capacities.
1.46 + ///
1.47 + /// The type of the map that stores the arc capacities.
1.48 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
1.49 + typedef CAP CapacityMap;
1.50 +
1.51 + /// \brief The type of the flow values.
1.52 + typedef typename CapacityMap::Value Value;
1.53 +
1.54 + /// \brief The type of the map that stores the flow values.
1.55 + ///
1.56 + /// The type of the map that stores the flow values.
1.57 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1.58 +#ifdef DOXYGEN
1.59 + typedef GR::ArcMap<Value> FlowMap;
1.60 +#else
1.61 + typedef typename Digraph::template ArcMap<Value> FlowMap;
1.62 +#endif
1.63 +
1.64 + /// \brief Instantiates a FlowMap.
1.65 + ///
1.66 + /// This function instantiates a \ref FlowMap.
1.67 + /// \param digraph The digraph for which we would like to define
1.68 + /// the flow map.
1.69 + static FlowMap* createFlowMap(const Digraph& digraph) {
1.70 + return new FlowMap(digraph);
1.71 + }
1.72 +
1.73 + /// \brief The tolerance used by the algorithm
1.74 + ///
1.75 + /// The tolerance used by the algorithm to handle inexact computation.
1.76 + typedef lemon::Tolerance<Value> Tolerance;
1.77 +
1.78 + };
1.79 +
1.80 + /// \ingroup max_flow
1.81 + ///
1.82 + /// \brief Edmonds-Karp algorithms class.
1.83 + ///
1.84 + /// This class provides an implementation of the \e Edmonds-Karp \e
1.85 + /// algorithm producing a \ref max_flow "flow of maximum value" in a
1.86 + /// digraph \cite clrs01algorithms, \cite amo93networkflows,
1.87 + /// \cite edmondskarp72theoretical.
1.88 + /// The Edmonds-Karp algorithm is slower than the Preflow
1.89 + /// algorithm, but it has an advantage of the step-by-step execution
1.90 + /// control with feasible flow solutions. The \e source node, the \e
1.91 + /// target node, the \e capacity of the arcs and the \e starting \e
1.92 + /// flow value of the arcs should be passed to the algorithm
1.93 + /// through the constructor.
1.94 + ///
1.95 + /// The time complexity of the algorithm is \f$ O(nm^2) \f$ in
1.96 + /// worst case. Always try the Preflow algorithm instead of this if
1.97 + /// you just want to compute the optimal flow.
1.98 + ///
1.99 + /// \tparam GR The type of the digraph the algorithm runs on.
1.100 + /// \tparam CAP The type of the capacity map. The default map
1.101 + /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
1.102 + /// \tparam TR The traits class that defines various types used by the
1.103 + /// algorithm. By default, it is \ref EdmondsKarpDefaultTraits
1.104 + /// "EdmondsKarpDefaultTraits<GR, CAP>".
1.105 + /// In most cases, this parameter should not be set directly,
1.106 + /// consider to use the named template parameters instead.
1.107 +
1.108 +#ifdef DOXYGEN
1.109 + template <typename GR, typename CAP, typename TR>
1.110 +#else
1.111 + template <typename GR,
1.112 + typename CAP = typename GR::template ArcMap<int>,
1.113 + typename TR = EdmondsKarpDefaultTraits<GR, CAP> >
1.114 +#endif
1.115 + class EdmondsKarp {
1.116 + public:
1.117 +
1.118 + /// \brief The \ref lemon::EdmondsKarpDefaultTraits "traits class"
1.119 + /// of the algorithm.
1.120 + typedef TR Traits;
1.121 + /// The type of the digraph the algorithm runs on.
1.122 + typedef typename Traits::Digraph Digraph;
1.123 + /// The type of the capacity map.
1.124 + typedef typename Traits::CapacityMap CapacityMap;
1.125 + /// The type of the flow values.
1.126 + typedef typename Traits::Value Value;
1.127 +
1.128 + /// The type of the flow map.
1.129 + typedef typename Traits::FlowMap FlowMap;
1.130 + /// The type of the tolerance.
1.131 + typedef typename Traits::Tolerance Tolerance;
1.132 +
1.133 + private:
1.134 +
1.135 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.136 + typedef typename Digraph::template NodeMap<Arc> PredMap;
1.137 +
1.138 + const Digraph& _graph;
1.139 + const CapacityMap* _capacity;
1.140 +
1.141 + Node _source, _target;
1.142 +
1.143 + FlowMap* _flow;
1.144 + bool _local_flow;
1.145 +
1.146 + PredMap* _pred;
1.147 + std::vector<Node> _queue;
1.148 +
1.149 + Tolerance _tolerance;
1.150 + Value _flow_value;
1.151 +
1.152 + void createStructures() {
1.153 + if (!_flow) {
1.154 + _flow = Traits::createFlowMap(_graph);
1.155 + _local_flow = true;
1.156 + }
1.157 + if (!_pred) {
1.158 + _pred = new PredMap(_graph);
1.159 + }
1.160 + _queue.resize(countNodes(_graph));
1.161 + }
1.162 +
1.163 + void destroyStructures() {
1.164 + if (_local_flow) {
1.165 + delete _flow;
1.166 + }
1.167 + if (_pred) {
1.168 + delete _pred;
1.169 + }
1.170 + }
1.171 +
1.172 + public:
1.173 +
1.174 + typedef EdmondsKarp Create;
1.175 +
1.176 + ///\name Named template parameters
1.177 +
1.178 + ///@{
1.179 +
1.180 + template <typename T>
1.181 + struct SetFlowMapTraits : public Traits {
1.182 + typedef T FlowMap;
1.183 + static FlowMap *createFlowMap(const Digraph&) {
1.184 + LEMON_ASSERT(false, "FlowMap is not initialized");
1.185 + return 0;
1.186 + }
1.187 + };
1.188 +
1.189 + /// \brief \ref named-templ-param "Named parameter" for setting
1.190 + /// FlowMap type
1.191 + ///
1.192 + /// \ref named-templ-param "Named parameter" for setting FlowMap
1.193 + /// type
1.194 + template <typename T>
1.195 + struct SetFlowMap
1.196 + : public EdmondsKarp<Digraph, CapacityMap, SetFlowMapTraits<T> > {
1.197 + typedef EdmondsKarp<Digraph, CapacityMap, SetFlowMapTraits<T> > Create;
1.198 + };
1.199 +
1.200 + /// @}
1.201 +
1.202 + protected:
1.203 +
1.204 + EdmondsKarp() {}
1.205 +
1.206 + public:
1.207 +
1.208 + /// \brief The constructor of the class.
1.209 + ///
1.210 + /// The constructor of the class.
1.211 + /// \param digraph The digraph the algorithm runs on.
1.212 + /// \param capacity The capacity of the arcs.
1.213 + /// \param source The source node.
1.214 + /// \param target The target node.
1.215 + EdmondsKarp(const Digraph& digraph, const CapacityMap& capacity,
1.216 + Node source, Node target)
1.217 + : _graph(digraph), _capacity(&capacity), _source(source), _target(target),
1.218 + _flow(0), _local_flow(false), _pred(0), _tolerance(), _flow_value()
1.219 + {
1.220 + LEMON_ASSERT(_source != _target,
1.221 + "Flow source and target are the same nodes.");
1.222 + }
1.223 +
1.224 + /// \brief Destructor.
1.225 + ///
1.226 + /// Destructor.
1.227 + ~EdmondsKarp() {
1.228 + destroyStructures();
1.229 + }
1.230 +
1.231 + /// \brief Sets the capacity map.
1.232 + ///
1.233 + /// Sets the capacity map.
1.234 + /// \return <tt>(*this)</tt>
1.235 + EdmondsKarp& capacityMap(const CapacityMap& map) {
1.236 + _capacity = ↦
1.237 + return *this;
1.238 + }
1.239 +
1.240 + /// \brief Sets the flow map.
1.241 + ///
1.242 + /// Sets the flow map.
1.243 + /// If you don't use this function before calling \ref run() or
1.244 + /// \ref init(), an instance will be allocated automatically.
1.245 + /// The destructor deallocates this automatically allocated map,
1.246 + /// of course.
1.247 + /// \return <tt>(*this)</tt>
1.248 + EdmondsKarp& flowMap(FlowMap& map) {
1.249 + if (_local_flow) {
1.250 + delete _flow;
1.251 + _local_flow = false;
1.252 + }
1.253 + _flow = ↦
1.254 + return *this;
1.255 + }
1.256 +
1.257 + /// \brief Sets the source node.
1.258 + ///
1.259 + /// Sets the source node.
1.260 + /// \return <tt>(*this)</tt>
1.261 + EdmondsKarp& source(const Node& node) {
1.262 + _source = node;
1.263 + return *this;
1.264 + }
1.265 +
1.266 + /// \brief Sets the target node.
1.267 + ///
1.268 + /// Sets the target node.
1.269 + /// \return <tt>(*this)</tt>
1.270 + EdmondsKarp& target(const Node& node) {
1.271 + _target = node;
1.272 + return *this;
1.273 + }
1.274 +
1.275 + /// \brief Sets the tolerance used by algorithm.
1.276 + ///
1.277 + /// Sets the tolerance used by algorithm.
1.278 + /// \return <tt>(*this)</tt>
1.279 + EdmondsKarp& tolerance(const Tolerance& tolerance) {
1.280 + _tolerance = tolerance;
1.281 + return *this;
1.282 + }
1.283 +
1.284 + /// \brief Returns a const reference to the tolerance.
1.285 + ///
1.286 + /// Returns a const reference to the tolerance object used by
1.287 + /// the algorithm.
1.288 + const Tolerance& tolerance() const {
1.289 + return _tolerance;
1.290 + }
1.291 +
1.292 + /// \name Execution control
1.293 + /// The simplest way to execute the algorithm is to use \ref run().\n
1.294 + /// If you need better control on the initial solution or the execution,
1.295 + /// you have to call one of the \ref init() functions first, then
1.296 + /// \ref start() or multiple times the \ref augment() function.
1.297 +
1.298 + ///@{
1.299 +
1.300 + /// \brief Initializes the algorithm.
1.301 + ///
1.302 + /// Initializes the internal data structures and sets the initial
1.303 + /// flow to zero on each arc.
1.304 + void init() {
1.305 + createStructures();
1.306 + for (ArcIt it(_graph); it != INVALID; ++it) {
1.307 + _flow->set(it, 0);
1.308 + }
1.309 + _flow_value = 0;
1.310 + }
1.311 +
1.312 + /// \brief Initializes the algorithm using the given flow map.
1.313 + ///
1.314 + /// Initializes the internal data structures and sets the initial
1.315 + /// flow to the given \c flowMap. The \c flowMap should
1.316 + /// contain a feasible flow, i.e. at each node excluding the source
1.317 + /// and the target, the incoming flow should be equal to the
1.318 + /// outgoing flow.
1.319 + template <typename FlowMap>
1.320 + void init(const FlowMap& flowMap) {
1.321 + createStructures();
1.322 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.323 + _flow->set(e, flowMap[e]);
1.324 + }
1.325 + _flow_value = 0;
1.326 + for (OutArcIt jt(_graph, _source); jt != INVALID; ++jt) {
1.327 + _flow_value += (*_flow)[jt];
1.328 + }
1.329 + for (InArcIt jt(_graph, _source); jt != INVALID; ++jt) {
1.330 + _flow_value -= (*_flow)[jt];
1.331 + }
1.332 + }
1.333 +
1.334 + /// \brief Initializes the algorithm using the given flow map.
1.335 + ///
1.336 + /// Initializes the internal data structures and sets the initial
1.337 + /// flow to the given \c flowMap. The \c flowMap should
1.338 + /// contain a feasible flow, i.e. at each node excluding the source
1.339 + /// and the target, the incoming flow should be equal to the
1.340 + /// outgoing flow.
1.341 + /// \return \c false when the given \c flowMap does not contain a
1.342 + /// feasible flow.
1.343 + template <typename FlowMap>
1.344 + bool checkedInit(const FlowMap& flowMap) {
1.345 + createStructures();
1.346 + for (ArcIt e(_graph); e != INVALID; ++e) {
1.347 + _flow->set(e, flowMap[e]);
1.348 + }
1.349 + for (NodeIt it(_graph); it != INVALID; ++it) {
1.350 + if (it == _source || it == _target) continue;
1.351 + Value outFlow = 0;
1.352 + for (OutArcIt jt(_graph, it); jt != INVALID; ++jt) {
1.353 + outFlow += (*_flow)[jt];
1.354 + }
1.355 + Value inFlow = 0;
1.356 + for (InArcIt jt(_graph, it); jt != INVALID; ++jt) {
1.357 + inFlow += (*_flow)[jt];
1.358 + }
1.359 + if (_tolerance.different(outFlow, inFlow)) {
1.360 + return false;
1.361 + }
1.362 + }
1.363 + for (ArcIt it(_graph); it != INVALID; ++it) {
1.364 + if (_tolerance.less((*_flow)[it], 0)) return false;
1.365 + if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
1.366 + }
1.367 + _flow_value = 0;
1.368 + for (OutArcIt jt(_graph, _source); jt != INVALID; ++jt) {
1.369 + _flow_value += (*_flow)[jt];
1.370 + }
1.371 + for (InArcIt jt(_graph, _source); jt != INVALID; ++jt) {
1.372 + _flow_value -= (*_flow)[jt];
1.373 + }
1.374 + return true;
1.375 + }
1.376 +
1.377 + /// \brief Augments the solution along a shortest path.
1.378 + ///
1.379 + /// Augments the solution along a shortest path. This function searches a
1.380 + /// shortest path between the source and the target
1.381 + /// in the residual digraph by the Bfs algoritm.
1.382 + /// Then it increases the flow on this path with the minimal residual
1.383 + /// capacity on the path. If there is no such path, it gives back
1.384 + /// false.
1.385 + /// \return \c false when the augmenting did not success, i.e. the
1.386 + /// current flow is a feasible and optimal solution.
1.387 + bool augment() {
1.388 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.389 + _pred->set(n, INVALID);
1.390 + }
1.391 +
1.392 + int first = 0, last = 1;
1.393 +
1.394 + _queue[0] = _source;
1.395 + _pred->set(_source, OutArcIt(_graph, _source));
1.396 +
1.397 + while (first != last && (*_pred)[_target] == INVALID) {
1.398 + Node n = _queue[first++];
1.399 +
1.400 + for (OutArcIt e(_graph, n); e != INVALID; ++e) {
1.401 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.402 + Node t = _graph.target(e);
1.403 + if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
1.404 + _pred->set(t, e);
1.405 + _queue[last++] = t;
1.406 + }
1.407 + }
1.408 + for (InArcIt e(_graph, n); e != INVALID; ++e) {
1.409 + Value rem = (*_flow)[e];
1.410 + Node t = _graph.source(e);
1.411 + if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
1.412 + _pred->set(t, e);
1.413 + _queue[last++] = t;
1.414 + }
1.415 + }
1.416 + }
1.417 +
1.418 + if ((*_pred)[_target] != INVALID) {
1.419 + Node n = _target;
1.420 + Arc e = (*_pred)[n];
1.421 +
1.422 + Value prem = (*_capacity)[e] - (*_flow)[e];
1.423 + n = _graph.source(e);
1.424 + while (n != _source) {
1.425 + e = (*_pred)[n];
1.426 + if (_graph.target(e) == n) {
1.427 + Value rem = (*_capacity)[e] - (*_flow)[e];
1.428 + if (rem < prem) prem = rem;
1.429 + n = _graph.source(e);
1.430 + } else {
1.431 + Value rem = (*_flow)[e];
1.432 + if (rem < prem) prem = rem;
1.433 + n = _graph.target(e);
1.434 + }
1.435 + }
1.436 +
1.437 + n = _target;
1.438 + e = (*_pred)[n];
1.439 +
1.440 + _flow->set(e, (*_flow)[e] + prem);
1.441 + n = _graph.source(e);
1.442 + while (n != _source) {
1.443 + e = (*_pred)[n];
1.444 + if (_graph.target(e) == n) {
1.445 + _flow->set(e, (*_flow)[e] + prem);
1.446 + n = _graph.source(e);
1.447 + } else {
1.448 + _flow->set(e, (*_flow)[e] - prem);
1.449 + n = _graph.target(e);
1.450 + }
1.451 + }
1.452 +
1.453 + _flow_value += prem;
1.454 + return true;
1.455 + } else {
1.456 + return false;
1.457 + }
1.458 + }
1.459 +
1.460 + /// \brief Executes the algorithm
1.461 + ///
1.462 + /// Executes the algorithm by performing augmenting phases until the
1.463 + /// optimal solution is reached.
1.464 + /// \pre One of the \ref init() functions must be called before
1.465 + /// using this function.
1.466 + void start() {
1.467 + while (augment()) {}
1.468 + }
1.469 +
1.470 + /// \brief Runs the algorithm.
1.471 + ///
1.472 + /// Runs the Edmonds-Karp algorithm.
1.473 + /// \note ek.run() is just a shortcut of the following code.
1.474 + ///\code
1.475 + /// ek.init();
1.476 + /// ek.start();
1.477 + ///\endcode
1.478 + void run() {
1.479 + init();
1.480 + start();
1.481 + }
1.482 +
1.483 + /// @}
1.484 +
1.485 + /// \name Query Functions
1.486 + /// The result of the Edmonds-Karp algorithm can be obtained using these
1.487 + /// functions.\n
1.488 + /// Either \ref run() or \ref start() should be called before using them.
1.489 +
1.490 + ///@{
1.491 +
1.492 + /// \brief Returns the value of the maximum flow.
1.493 + ///
1.494 + /// Returns the value of the maximum flow found by the algorithm.
1.495 + ///
1.496 + /// \pre Either \ref run() or \ref init() must be called before
1.497 + /// using this function.
1.498 + Value flowValue() const {
1.499 + return _flow_value;
1.500 + }
1.501 +
1.502 + /// \brief Returns the flow value on the given arc.
1.503 + ///
1.504 + /// Returns the flow value on the given arc.
1.505 + ///
1.506 + /// \pre Either \ref run() or \ref init() must be called before
1.507 + /// using this function.
1.508 + Value flow(const Arc& arc) const {
1.509 + return (*_flow)[arc];
1.510 + }
1.511 +
1.512 + /// \brief Returns a const reference to the flow map.
1.513 + ///
1.514 + /// Returns a const reference to the arc map storing the found flow.
1.515 + ///
1.516 + /// \pre Either \ref run() or \ref init() must be called before
1.517 + /// using this function.
1.518 + const FlowMap& flowMap() const {
1.519 + return *_flow;
1.520 + }
1.521 +
1.522 + /// \brief Returns \c true when the node is on the source side of the
1.523 + /// minimum cut.
1.524 + ///
1.525 + /// Returns true when the node is on the source side of the found
1.526 + /// minimum cut.
1.527 + ///
1.528 + /// \pre Either \ref run() or \ref init() must be called before
1.529 + /// using this function.
1.530 + bool minCut(const Node& node) const {
1.531 + return ((*_pred)[node] != INVALID) || node == _source;
1.532 + }
1.533 +
1.534 + /// \brief Gives back a minimum value cut.
1.535 + ///
1.536 + /// Sets \c cutMap to the characteristic vector of a minimum value
1.537 + /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
1.538 + /// node map with \c bool (or convertible) value type.
1.539 + ///
1.540 + /// \note This function calls \ref minCut() for each node, so it runs in
1.541 + /// O(n) time.
1.542 + ///
1.543 + /// \pre Either \ref run() or \ref init() must be called before
1.544 + /// using this function.
1.545 + template <typename CutMap>
1.546 + void minCutMap(CutMap& cutMap) const {
1.547 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.548 + cutMap.set(n, (*_pred)[n] != INVALID);
1.549 + }
1.550 + cutMap.set(_source, true);
1.551 + }
1.552 +
1.553 + /// @}
1.554 +
1.555 + };
1.556 +
1.557 +}
1.558 +
1.559 +#endif