Port Edmonds-Karp algorithm from svn -r3524 (#177)
authorAntal Nemes <thoneyvazul@gmail.com>
Tue, 30 Nov 2010 20:21:52 +0100
changeset 105692a884824429
parent 1055 05b34170866b
child 1057 6a8a688eacf6
Port Edmonds-Karp algorithm from svn -r3524 (#177)
lemon/edmonds_karp.h
test/CMakeLists.txt
test/edmonds_karp_test.cc
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lemon/edmonds_karp.h	Tue Nov 30 20:21:52 2010 +0100
     1.3 @@ -0,0 +1,515 @@
     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-2010
     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 length of the arcs.
    1.52 +    typedef typename CapacityMap::Value Value;
    1.53 +
    1.54 +    /// \brief The map type that stores the flow values.
    1.55 +    ///
    1.56 +    /// The map type that stores the flow values. 
    1.57 +    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
    1.58 +    typedef typename Digraph::template ArcMap<Value> FlowMap;
    1.59 +
    1.60 +    /// \brief Instantiates a FlowMap.
    1.61 +    ///
    1.62 +    /// This function instantiates a \ref FlowMap. 
    1.63 +    /// \param digraph The digraph, to which we would like to define the flow map.
    1.64 +    static FlowMap* createFlowMap(const Digraph& digraph) {
    1.65 +      return new FlowMap(digraph);
    1.66 +    }
    1.67 +
    1.68 +    /// \brief The tolerance used by the algorithm
    1.69 +    ///
    1.70 +    /// The tolerance used by the algorithm to handle inexact computation.
    1.71 +    typedef lemon::Tolerance<Value> Tolerance;
    1.72 +
    1.73 +  };
    1.74 +
    1.75 +  /// \ingroup max_flow
    1.76 +  ///
    1.77 +  /// \brief Edmonds-Karp algorithms class.
    1.78 +  ///
    1.79 +  /// This class provides an implementation of the \e Edmonds-Karp \e
    1.80 +  /// algorithm producing a flow of maximum value in directed
    1.81 +  /// digraphs. The Edmonds-Karp algorithm is slower than the Preflow
    1.82 +  /// algorithm but it has an advantage of the step-by-step execution
    1.83 +  /// control with feasible flow solutions. The \e source node, the \e
    1.84 +  /// target node, the \e capacity of the arcs and the \e starting \e
    1.85 +  /// flow value of the arcs should be passed to the algorithm
    1.86 +  /// through the constructor.
    1.87 +  ///
    1.88 +  /// The time complexity of the algorithm is \f$ O(nm^2) \f$ in
    1.89 +  /// worst case.  Always try the preflow algorithm instead of this if
    1.90 +  /// you just want to compute the optimal flow.
    1.91 +  ///
    1.92 +  /// \param GR The digraph type the algorithm runs on.
    1.93 +  /// \param CAP The capacity map type.
    1.94 +  /// \param TR Traits class to set various data types used by
    1.95 +  /// the algorithm.  The default traits class is \ref
    1.96 +  /// EdmondsKarpDefaultTraits.  See \ref EdmondsKarpDefaultTraits for the
    1.97 +  /// documentation of a Edmonds-Karp traits class. 
    1.98 +
    1.99 +#ifdef DOXYGEN
   1.100 +  template <typename GR, typename CAP, typename TR>
   1.101 +#else 
   1.102 +  template <typename GR,
   1.103 +	    typename CAP = typename GR::template ArcMap<int>,
   1.104 +            typename TR = EdmondsKarpDefaultTraits<GR, CAP> >
   1.105 +#endif
   1.106 +  class EdmondsKarp {
   1.107 +  public:
   1.108 +
   1.109 +    typedef TR Traits;
   1.110 +    typedef typename Traits::Digraph Digraph;
   1.111 +    typedef typename Traits::CapacityMap CapacityMap;
   1.112 +    typedef typename Traits::Value Value; 
   1.113 +
   1.114 +    typedef typename Traits::FlowMap FlowMap;
   1.115 +    typedef typename Traits::Tolerance Tolerance;
   1.116 +
   1.117 +  private:
   1.118 +
   1.119 +    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
   1.120 +    typedef typename Digraph::template NodeMap<Arc> PredMap;
   1.121 +    
   1.122 +    const Digraph& _graph;
   1.123 +    const CapacityMap* _capacity;
   1.124 +
   1.125 +    Node _source, _target;
   1.126 +
   1.127 +    FlowMap* _flow;
   1.128 +    bool _local_flow;
   1.129 +
   1.130 +    PredMap* _pred;
   1.131 +    std::vector<Node> _queue;
   1.132 +    
   1.133 +    Tolerance _tolerance;
   1.134 +    Value _flow_value;
   1.135 +
   1.136 +    void createStructures() {
   1.137 +      if (!_flow) {
   1.138 +	_flow = Traits::createFlowMap(_graph);
   1.139 +	_local_flow = true;
   1.140 +      }
   1.141 +      if (!_pred) {
   1.142 +	_pred = new PredMap(_graph);
   1.143 +      }
   1.144 +      _queue.resize(countNodes(_graph));
   1.145 +    }
   1.146 +
   1.147 +    void destroyStructures() {
   1.148 +      if (_local_flow) {
   1.149 +	delete _flow;
   1.150 +      }
   1.151 +      if (_pred) {
   1.152 +	delete _pred;
   1.153 +      }
   1.154 +    }
   1.155 +    
   1.156 +  public:
   1.157 +
   1.158 +    ///\name Named template parameters
   1.159 +
   1.160 +    ///@{
   1.161 +
   1.162 +    template <typename T>
   1.163 +    struct DefFlowMapTraits : public Traits {
   1.164 +      typedef T FlowMap;
   1.165 +      static FlowMap *createFlowMap(const Digraph&) {
   1.166 +	LEMON_ASSERT(false,"Uninitialized parameter.");
   1.167 +        return 0;
   1.168 +      }
   1.169 +    };
   1.170 +
   1.171 +    /// \brief \ref named-templ-param "Named parameter" for setting
   1.172 +    /// FlowMap type
   1.173 +    ///
   1.174 +    /// \ref named-templ-param "Named parameter" for setting FlowMap
   1.175 +    /// type
   1.176 +    template <typename T>
   1.177 +    struct DefFlowMap 
   1.178 +      : public EdmondsKarp<Digraph, CapacityMap, DefFlowMapTraits<T> > {
   1.179 +      typedef EdmondsKarp<Digraph, CapacityMap, DefFlowMapTraits<T> > 
   1.180 +      Create;
   1.181 +    };
   1.182 +
   1.183 +
   1.184 +    /// @}
   1.185 +
   1.186 +  protected:
   1.187 +    
   1.188 +    EdmondsKarp() {}
   1.189 +
   1.190 +  public:
   1.191 +
   1.192 +    /// \brief The constructor of the class.
   1.193 +    ///
   1.194 +    /// The constructor of the class. 
   1.195 +    /// \param digraph The digraph the algorithm runs on. 
   1.196 +    /// \param capacity The capacity of the arcs. 
   1.197 +    /// \param source The source node.
   1.198 +    /// \param target The target node.
   1.199 +    EdmondsKarp(const Digraph& digraph, const CapacityMap& capacity,
   1.200 +		Node source, Node target)
   1.201 +      : _graph(digraph), _capacity(&capacity), _source(source), _target(target),
   1.202 +	_flow(0), _local_flow(false), _pred(0), _tolerance(), _flow_value()
   1.203 +    {
   1.204 +      LEMON_ASSERT(_source != _target,"Flow source and target are the same nodes.");
   1.205 +    }
   1.206 +
   1.207 +    /// \brief Destructor.
   1.208 +    ///
   1.209 +    /// Destructor.
   1.210 +    ~EdmondsKarp() {
   1.211 +      destroyStructures();
   1.212 +    }
   1.213 +
   1.214 +    /// \brief Sets the capacity map.
   1.215 +    ///
   1.216 +    /// Sets the capacity map.
   1.217 +    /// \return \c (*this)
   1.218 +    EdmondsKarp& capacityMap(const CapacityMap& map) {
   1.219 +      _capacity = &map;
   1.220 +      return *this;
   1.221 +    }
   1.222 +
   1.223 +    /// \brief Sets the flow map.
   1.224 +    ///
   1.225 +    /// Sets the flow map.
   1.226 +    /// \return \c (*this)
   1.227 +    EdmondsKarp& flowMap(FlowMap& map) {
   1.228 +      if (_local_flow) {
   1.229 +	delete _flow;
   1.230 +	_local_flow = false;
   1.231 +      }
   1.232 +      _flow = &map;
   1.233 +      return *this;
   1.234 +    }
   1.235 +
   1.236 +    /// \brief Returns the flow map.
   1.237 +    ///
   1.238 +    /// \return The flow map.
   1.239 +    const FlowMap& flowMap() const {
   1.240 +      return *_flow;
   1.241 +    }
   1.242 +
   1.243 +    /// \brief Sets the source node.
   1.244 +    ///
   1.245 +    /// Sets the source node.
   1.246 +    /// \return \c (*this)
   1.247 +    EdmondsKarp& source(const Node& node) {
   1.248 +      _source = node;
   1.249 +      return *this;
   1.250 +    }
   1.251 +
   1.252 +    /// \brief Sets the target node.
   1.253 +    ///
   1.254 +    /// Sets the target node.
   1.255 +    /// \return \c (*this)
   1.256 +    EdmondsKarp& target(const Node& node) {
   1.257 +      _target = node;
   1.258 +      return *this;
   1.259 +    }
   1.260 +
   1.261 +    /// \brief Sets the tolerance used by algorithm.
   1.262 +    ///
   1.263 +    /// Sets the tolerance used by algorithm.
   1.264 +    EdmondsKarp& tolerance(const Tolerance& tolerance) {
   1.265 +      _tolerance = tolerance;
   1.266 +      return *this;
   1.267 +    } 
   1.268 +
   1.269 +    /// \brief Returns the tolerance used by algorithm.
   1.270 +    ///
   1.271 +    /// Returns the tolerance used by algorithm.
   1.272 +    const Tolerance& tolerance() const {
   1.273 +      return _tolerance;
   1.274 +    } 
   1.275 +
   1.276 +    /// \name Execution control
   1.277 +    /// The simplest way to execute the
   1.278 +    /// algorithm is to use the \c run() member functions.
   1.279 +    /// \n
   1.280 +    /// If you need more control on initial solution or
   1.281 +    /// execution then you have to call one \ref init() function and then
   1.282 +    /// the start() or multiple times the \c augment() member function.  
   1.283 +    
   1.284 +    ///@{
   1.285 +
   1.286 +    /// \brief Initializes the algorithm
   1.287 +    /// 
   1.288 +    /// Sets the flow to empty flow.
   1.289 +    void init() {
   1.290 +      createStructures();
   1.291 +      for (ArcIt it(_graph); it != INVALID; ++it) {
   1.292 +        _flow->set(it, 0);
   1.293 +      }
   1.294 +      _flow_value = 0;
   1.295 +    }
   1.296 +    
   1.297 +    /// \brief Initializes the algorithm
   1.298 +    /// 
   1.299 +    /// Initializes the flow to the \c flowMap. The \c flowMap should
   1.300 +    /// contain a feasible flow, ie. in each node excluding the source
   1.301 +    /// and the target the incoming flow should be equal to the
   1.302 +    /// outgoing flow.
   1.303 +    template <typename FlowMap>
   1.304 +    void flowInit(const FlowMap& flowMap) {
   1.305 +      createStructures();
   1.306 +      for (ArcIt e(_graph); e != INVALID; ++e) {
   1.307 +	_flow->set(e, flowMap[e]);
   1.308 +      }
   1.309 +      _flow_value = 0;
   1.310 +      for (OutArcIt jt(_graph, _source); jt != INVALID; ++jt) {
   1.311 +        _flow_value += (*_flow)[jt];
   1.312 +      }
   1.313 +      for (InArcIt jt(_graph, _source); jt != INVALID; ++jt) {
   1.314 +        _flow_value -= (*_flow)[jt];
   1.315 +      }
   1.316 +    }
   1.317 +
   1.318 +    /// \brief Initializes the algorithm
   1.319 +    /// 
   1.320 +    /// Initializes the flow to the \c flowMap. The \c flowMap should
   1.321 +    /// contain a feasible flow, ie. in each node excluding the source
   1.322 +    /// and the target the incoming flow should be equal to the
   1.323 +    /// outgoing flow.  
   1.324 +    /// \return %False when the given flowMap does not contain
   1.325 +    /// feasible flow.
   1.326 +    template <typename FlowMap>
   1.327 +    bool checkedFlowInit(const FlowMap& flowMap) {
   1.328 +      createStructures();
   1.329 +      for (ArcIt e(_graph); e != INVALID; ++e) {
   1.330 +	_flow->set(e, flowMap[e]);
   1.331 +      }
   1.332 +      for (NodeIt it(_graph); it != INVALID; ++it) {
   1.333 +        if (it == _source || it == _target) continue;
   1.334 +        Value outFlow = 0;
   1.335 +        for (OutArcIt jt(_graph, it); jt != INVALID; ++jt) {
   1.336 +          outFlow += (*_flow)[jt];
   1.337 +        }
   1.338 +        Value inFlow = 0;
   1.339 +        for (InArcIt jt(_graph, it); jt != INVALID; ++jt) {
   1.340 +          inFlow += (*_flow)[jt];
   1.341 +        }
   1.342 +        if (_tolerance.different(outFlow, inFlow)) {
   1.343 +          return false;
   1.344 +        }
   1.345 +      }
   1.346 +      for (ArcIt it(_graph); it != INVALID; ++it) {
   1.347 +        if (_tolerance.less((*_flow)[it], 0)) return false;
   1.348 +        if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
   1.349 +      }
   1.350 +      _flow_value = 0;
   1.351 +      for (OutArcIt jt(_graph, _source); jt != INVALID; ++jt) {
   1.352 +        _flow_value += (*_flow)[jt];
   1.353 +      }
   1.354 +      for (InArcIt jt(_graph, _source); jt != INVALID; ++jt) {
   1.355 +        _flow_value -= (*_flow)[jt];
   1.356 +      }
   1.357 +      return true;
   1.358 +    }
   1.359 +
   1.360 +    /// \brief Augment the solution on an arc shortest path.
   1.361 +    /// 
   1.362 +    /// Augment the solution on an arc shortest path. It searches an
   1.363 +    /// arc shortest path between the source and the target
   1.364 +    /// in the residual digraph by the bfs algoritm.
   1.365 +    /// Then it increases the flow on this path with the minimal residual
   1.366 +    /// capacity on the path. If there is no such path it gives back
   1.367 +    /// false.
   1.368 +    /// \return %False when the augmenting didn't success so the
   1.369 +    /// current flow is a feasible and optimal solution.
   1.370 +    bool augment() {
   1.371 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.372 +	_pred->set(n, INVALID);
   1.373 +      }
   1.374 +      
   1.375 +      int first = 0, last = 1;
   1.376 +      
   1.377 +      _queue[0] = _source;
   1.378 +      _pred->set(_source, OutArcIt(_graph, _source));
   1.379 +
   1.380 +      while (first != last && (*_pred)[_target] == INVALID) {
   1.381 +	Node n = _queue[first++];
   1.382 +	
   1.383 +	for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   1.384 +	  Value rem = (*_capacity)[e] - (*_flow)[e];
   1.385 +	  Node t = _graph.target(e);
   1.386 +	  if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
   1.387 +	    _pred->set(t, e);
   1.388 +	    _queue[last++] = t;
   1.389 +	  }
   1.390 +	}
   1.391 +	for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.392 +	  Value rem = (*_flow)[e];
   1.393 +	  Node t = _graph.source(e);
   1.394 +	  if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {
   1.395 +	    _pred->set(t, e);
   1.396 +	    _queue[last++] = t;
   1.397 +	  }
   1.398 +	}
   1.399 +      }
   1.400 +
   1.401 +      if ((*_pred)[_target] != INVALID) {
   1.402 +	Node n = _target;
   1.403 +	Arc e = (*_pred)[n];
   1.404 +
   1.405 +	Value prem = (*_capacity)[e] - (*_flow)[e];
   1.406 +	n = _graph.source(e);
   1.407 +	while (n != _source) {
   1.408 +	  e = (*_pred)[n];
   1.409 +	  if (_graph.target(e) == n) {
   1.410 +	    Value rem = (*_capacity)[e] - (*_flow)[e];
   1.411 +	    if (rem < prem) prem = rem;
   1.412 +	    n = _graph.source(e);
   1.413 +	  } else {
   1.414 +	    Value rem = (*_flow)[e];
   1.415 +	    if (rem < prem) prem = rem;
   1.416 +	    n = _graph.target(e);   
   1.417 +	  } 
   1.418 +	}
   1.419 +
   1.420 +	n = _target;
   1.421 +	e = (*_pred)[n];
   1.422 +
   1.423 +	_flow->set(e, (*_flow)[e] + prem);
   1.424 +	n = _graph.source(e);
   1.425 +	while (n != _source) {
   1.426 +	  e = (*_pred)[n];
   1.427 +	  if (_graph.target(e) == n) {
   1.428 +	    _flow->set(e, (*_flow)[e] + prem);
   1.429 +	    n = _graph.source(e);
   1.430 +	  } else {
   1.431 +	    _flow->set(e, (*_flow)[e] - prem);
   1.432 +	    n = _graph.target(e);   
   1.433 +	  } 
   1.434 +	}
   1.435 +
   1.436 +	_flow_value += prem;	
   1.437 +	return true;
   1.438 +      } else {
   1.439 +	return false;
   1.440 +      }
   1.441 +    }
   1.442 +
   1.443 +    /// \brief Executes the algorithm
   1.444 +    ///
   1.445 +    /// It runs augmenting phases until the optimal solution is reached. 
   1.446 +    void start() {
   1.447 +      while (augment()) {}
   1.448 +    }
   1.449 +
   1.450 +    /// \brief Runs the algorithm.
   1.451 +    /// 
   1.452 +    /// It is just a shorthand for:
   1.453 +    ///
   1.454 +    ///\code 
   1.455 +    /// ek.init();
   1.456 +    /// ek.start();
   1.457 +    ///\endcode
   1.458 +    void run() {
   1.459 +      init();
   1.460 +      start();
   1.461 +    }
   1.462 +
   1.463 +    /// @}
   1.464 +
   1.465 +    /// \name Query Functions
   1.466 +    /// The result of the Edmonds-Karp algorithm can be obtained using these
   1.467 +    /// functions.\n
   1.468 +    /// Before the use of these functions,
   1.469 +    /// either run() or start() must be called.
   1.470 +    
   1.471 +    ///@{
   1.472 +
   1.473 +    /// \brief Returns the value of the maximum flow.
   1.474 +    ///
   1.475 +    /// Returns the value of the maximum flow by returning the excess
   1.476 +    /// of the target node \c t.
   1.477 +
   1.478 +    Value flowValue() const {
   1.479 +      return _flow_value;
   1.480 +    }
   1.481 +
   1.482 +
   1.483 +    /// \brief Returns the flow on the arc.
   1.484 +    ///
   1.485 +    /// Sets the \c flowMap to the flow on the arcs.
   1.486 +    Value flow(const Arc& arc) const {
   1.487 +      return (*_flow)[arc];
   1.488 +    }
   1.489 +
   1.490 +    /// \brief Returns true when the node is on the source side of minimum cut.
   1.491 +    ///
   1.492 +
   1.493 +    /// Returns true when the node is on the source side of minimum
   1.494 +    /// cut.
   1.495 +
   1.496 +    bool minCut(const Node& node) const {
   1.497 +      return ((*_pred)[node] != INVALID) or node == _source;
   1.498 +    }
   1.499 +
   1.500 +    /// \brief Returns a minimum value cut.
   1.501 +    ///
   1.502 +    /// Sets \c cutMap to the characteristic vector of a minimum value cut.
   1.503 +
   1.504 +    template <typename CutMap>
   1.505 +    void minCutMap(CutMap& cutMap) const {
   1.506 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.507 +	cutMap.set(n, (*_pred)[n] != INVALID);
   1.508 +      }
   1.509 +      cutMap.set(_source, true);
   1.510 +    }    
   1.511 +
   1.512 +    /// @}
   1.513 +
   1.514 +  };
   1.515 +
   1.516 +}
   1.517 +
   1.518 +#endif
     2.1 --- a/test/CMakeLists.txt	Mon Mar 18 18:48:58 2013 +0100
     2.2 +++ b/test/CMakeLists.txt	Tue Nov 30 20:21:52 2010 +0100
     2.3 @@ -25,6 +25,7 @@
     2.4    dijkstra_test
     2.5    dim_test
     2.6    edge_set_test
     2.7 +  edmonds_karp_test
     2.8    error_test
     2.9    euler_test
    2.10    fractional_matching_test
     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/test/edmonds_karp_test.cc	Tue Nov 30 20:21:52 2010 +0100
     3.3 @@ -0,0 +1,236 @@
     3.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     3.5 + *
     3.6 + * This file is a part of LEMON, a generic C++ optimization library.
     3.7 + *
     3.8 + * Copyright (C) 2003-2010
     3.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    3.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    3.11 + *
    3.12 + * Permission to use, modify and distribute this software is granted
    3.13 + * provided that this copyright notice appears in all copies. For
    3.14 + * precise terms see the accompanying LICENSE file.
    3.15 + *
    3.16 + * This software is provided "AS IS" with no warranty of any kind,
    3.17 + * express or implied, and with no claim as to its suitability for any
    3.18 + * purpose.
    3.19 + *
    3.20 + */
    3.21 +
    3.22 +#include<iostream>
    3.23 +
    3.24 +#include "test_tools.h"
    3.25 +#include<lemon/smart_graph.h>
    3.26 +#include<lemon/edmonds_karp.h>
    3.27 +#include <lemon/concepts/digraph.h>
    3.28 +#include <lemon/concepts/maps.h>
    3.29 +#include <lemon/lgf_reader.h>
    3.30 +
    3.31 +using namespace lemon;
    3.32 +
    3.33 +char test_lgf[] =
    3.34 +  "@nodes\n"
    3.35 +  "label\n"
    3.36 +  "0\n"
    3.37 +  "1\n"
    3.38 +  "2\n"
    3.39 +  "3\n"
    3.40 +  "4\n"
    3.41 +  "5\n"
    3.42 +  "6\n"
    3.43 +  "7\n"
    3.44 +  "8\n"
    3.45 +  "9\n"
    3.46 +  "@arcs\n"
    3.47 +  "    label capacity\n"
    3.48 +  "0 1 0     20\n"
    3.49 +  "0 2 1     0\n"
    3.50 +  "1 1 2     3\n"
    3.51 +  "1 2 3     8\n"
    3.52 +  "1 3 4     8\n"
    3.53 +  "2 5 5     5\n"
    3.54 +  "3 2 6     5\n"
    3.55 +  "3 5 7     5\n"
    3.56 +  "3 6 8     5\n"
    3.57 +  "4 3 9     3\n"
    3.58 +  "5 7 10    3\n"
    3.59 +  "5 6 11    10\n"
    3.60 +  "5 8 12    10\n"
    3.61 +  "6 8 13    8\n"
    3.62 +  "8 9 14    20\n"
    3.63 +  "8 1 15    5\n"
    3.64 +  "9 5 16    5\n"
    3.65 +  "@attributes\n"
    3.66 +  "source 1\n"
    3.67 +  "target 8\n";
    3.68 +
    3.69 +void checkEdmondKarpCompile() {
    3.70 +  typedef int VType;
    3.71 +  typedef concepts::Digraph Digraph;
    3.72 +
    3.73 +  typedef Digraph::Node Node;
    3.74 +  typedef Digraph::Arc Arc;
    3.75 +  typedef concepts::ReadMap<Arc,VType> CapMap;
    3.76 +  typedef concepts::ReadWriteMap<Arc,VType> FlowMap;
    3.77 +  typedef concepts::WriteMap<Node,bool> CutMap;
    3.78 +
    3.79 +  Digraph g;
    3.80 +  Node n;
    3.81 +  Arc e;
    3.82 +  CapMap cap;
    3.83 +  FlowMap flow;
    3.84 +  CutMap cut;
    3.85 +  VType v;
    3.86 +  bool b;
    3.87 +  ignore_unused_variable_warning(v,b);
    3.88 +  typedef EdmondsKarp<Digraph, CapMap>
    3.89 +              ::DefFlowMap<FlowMap>
    3.90 +              ::Create EKType;
    3.91 +  EKType ek_test(g, cap, n, n);
    3.92 +  const EKType& const_ek_test = ek_test;
    3.93 +
    3.94 +  EKType::Tolerance tol = const_ek_test.tolerance();
    3.95 +  ek_test.tolerance(tol);
    3.96 +
    3.97 +  ek_test
    3.98 +    .capacityMap(cap)
    3.99 +    .flowMap(flow)
   3.100 +    .source(n)
   3.101 +    .target(n);
   3.102 +
   3.103 +  ek_test.init();
   3.104 +  ek_test.start();
   3.105 +
   3.106 +  v = const_ek_test.flowValue();
   3.107 +  v = const_ek_test.flow(e);
   3.108 +
   3.109 +  const FlowMap& fm = const_ek_test.flowMap();
   3.110 +  b = const_ek_test.minCut(n);
   3.111 +  const_ek_test.minCutMap(cut);
   3.112 +
   3.113 +  ignore_unused_variable_warning(fm);
   3.114 +}
   3.115 +
   3.116 +int cutValue (const SmartDigraph& g,
   3.117 +              const SmartDigraph::NodeMap<bool>& cut,
   3.118 +              const SmartDigraph::ArcMap<int>& cap) {
   3.119 +
   3.120 +  int c=0;
   3.121 +  for(SmartDigraph::ArcIt e(g); e!=INVALID; ++e) {
   3.122 +    if (cut[g.source(e)] && !cut[g.target(e)]) c+=cap[e];
   3.123 +  }
   3.124 +  return c;
   3.125 +}
   3.126 +
   3.127 +bool checkFlow(const SmartDigraph& g,
   3.128 +               const SmartDigraph::ArcMap<int>& flow,
   3.129 +               const SmartDigraph::ArcMap<int>& cap,
   3.130 +               SmartDigraph::Node s, SmartDigraph::Node t) {
   3.131 +
   3.132 +  for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) {
   3.133 +    if (flow[e] < 0 || flow[e] > cap[e]) return false;
   3.134 +  }
   3.135 +
   3.136 +  for (SmartDigraph::NodeIt n(g); n != INVALID; ++n) {
   3.137 +    if (n == s || n == t) continue;
   3.138 +    int sum = 0;
   3.139 +    for (SmartDigraph::OutArcIt e(g, n); e != INVALID; ++e) {
   3.140 +      sum += flow[e];
   3.141 +    }
   3.142 +    for (SmartDigraph::InArcIt e(g, n); e != INVALID; ++e) {
   3.143 +      sum -= flow[e];
   3.144 +    }
   3.145 +    if (sum != 0) return false;
   3.146 +  }
   3.147 +  return true;
   3.148 +}
   3.149 +
   3.150 +int main() {
   3.151 +
   3.152 +  typedef SmartDigraph Digraph;
   3.153 +
   3.154 +  typedef Digraph::Node Node;
   3.155 +  typedef Digraph::NodeIt NodeIt;
   3.156 +  typedef Digraph::ArcIt ArcIt;
   3.157 +  typedef Digraph::ArcMap<int> CapMap;
   3.158 +  typedef Digraph::ArcMap<int> FlowMap;
   3.159 +  typedef Digraph::NodeMap<bool> CutMap;
   3.160 +
   3.161 +  typedef EdmondsKarp<Digraph, CapMap> EKType;
   3.162 +
   3.163 +  Digraph g;
   3.164 +  Node s, t;
   3.165 +  CapMap cap(g);
   3.166 +  std::istringstream input(test_lgf);
   3.167 +  DigraphReader<Digraph>(g,input).
   3.168 +    arcMap("capacity", cap).
   3.169 +    node("source",s).
   3.170 +    node("target",t).
   3.171 +    run();
   3.172 +
   3.173 +  EKType ek_test(g, cap, s, t);
   3.174 +  ek_test.run();
   3.175 +
   3.176 +  check(checkFlow(g, ek_test.flowMap(), cap, s, t),
   3.177 +        "The flow is not feasible.");
   3.178 +
   3.179 +  CutMap min_cut(g);
   3.180 +  ek_test.minCutMap(min_cut);
   3.181 +  int min_cut_value=cutValue(g,min_cut,cap);
   3.182 +
   3.183 +  check(ek_test.flowValue() == min_cut_value,
   3.184 +        "The max flow value is not equal to the three min cut values.");
   3.185 +
   3.186 +  FlowMap flow(g);
   3.187 +  for(ArcIt e(g); e!=INVALID; ++e) flow[e] = ek_test.flowMap()[e];
   3.188 +
   3.189 +  int flow_value=ek_test.flowValue();
   3.190 +
   3.191 +  for(ArcIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e];
   3.192 +  ek_test.flowInit(flow);
   3.193 +  ek_test.start();
   3.194 +
   3.195 +  CutMap min_cut1(g);
   3.196 +  ek_test.minCutMap(min_cut1);
   3.197 +  min_cut_value=cutValue(g,min_cut1,cap);
   3.198 +
   3.199 +  check(ek_test.flowValue() == min_cut_value &&
   3.200 +        min_cut_value == 2*flow_value,
   3.201 +        "The max flow value or the min cut value is wrong.");
   3.202 +
   3.203 +  check(checkFlow(g, ek_test.flowMap(), cap, s, t),
   3.204 +        "The flow is not feasible.");
   3.205 +
   3.206 +  CutMap min_cut2(g);
   3.207 +  ek_test.minCutMap(min_cut2);
   3.208 +  min_cut_value=cutValue(g,min_cut2,cap);
   3.209 +
   3.210 +  check(ek_test.flowValue() == min_cut_value &&
   3.211 +        min_cut_value == 2*flow_value,
   3.212 +        "The max flow value or the three min cut values were not doubled.");
   3.213 +
   3.214 +
   3.215 +  ek_test.flowMap(flow);
   3.216 +
   3.217 +  NodeIt tmp1(g,s);
   3.218 +  ++tmp1;
   3.219 +  if ( tmp1 != INVALID ) s=tmp1;
   3.220 +
   3.221 +  NodeIt tmp2(g,t);
   3.222 +  ++tmp2;
   3.223 +  if ( tmp2 != INVALID ) t=tmp2;
   3.224 +
   3.225 +  ek_test.source(s);
   3.226 +  ek_test.target(t);
   3.227 +
   3.228 +  ek_test.run();
   3.229 +
   3.230 +  CutMap min_cut3(g);
   3.231 +  ek_test.minCutMap(min_cut3);
   3.232 +  min_cut_value=cutValue(g,min_cut3,cap);
   3.233 +
   3.234 +
   3.235 +  check(ek_test.flowValue() == min_cut_value,
   3.236 +        "The max flow value or the three min cut values are incorrect.");
   3.237 +
   3.238 +  return 0;
   3.239 +}