1.1 --- a/lemon/Makefile.am Tue Oct 28 18:33:51 2008 +0100
1.2 +++ b/lemon/Makefile.am Tue Oct 28 18:39:53 2008 +0000
1.3 @@ -40,6 +40,7 @@
1.4 lemon/path.h \
1.5 lemon/random.h \
1.6 lemon/smart_graph.h \
1.7 + lemon/suurballe.h \
1.8 lemon/time_measure.h \
1.9 lemon/tolerance.h \
1.10 lemon/unionfind.h
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/lemon/suurballe.h Tue Oct 28 18:39:53 2008 +0000
2.3 @@ -0,0 +1,499 @@
2.4 +/* -*- C++ -*-
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_SUURBALLE_H
2.23 +#define LEMON_SUURBALLE_H
2.24 +
2.25 +///\ingroup shortest_path
2.26 +///\file
2.27 +///\brief An algorithm for finding arc-disjoint paths between two
2.28 +/// nodes having minimum total length.
2.29 +
2.30 +#include <vector>
2.31 +#include <lemon/bin_heap.h>
2.32 +#include <lemon/path.h>
2.33 +
2.34 +namespace lemon {
2.35 +
2.36 + /// \addtogroup shortest_path
2.37 + /// @{
2.38 +
2.39 + /// \brief Implementation of an algorithm for finding arc-disjoint
2.40 + /// paths between two nodes having minimum total length.
2.41 + ///
2.42 + /// \ref lemon::Suurballe "Suurballe" implements an algorithm for
2.43 + /// finding arc-disjoint paths having minimum total length (cost)
2.44 + /// from a given source node to a given target node in a directed
2.45 + /// digraph.
2.46 + ///
2.47 + /// In fact, this implementation is the specialization of the
2.48 + /// \ref CapacityScaling "successive shortest path" algorithm.
2.49 + ///
2.50 + /// \tparam Digraph The directed digraph type the algorithm runs on.
2.51 + /// \tparam LengthMap The type of the length (cost) map.
2.52 + ///
2.53 + /// \warning Length values should be \e non-negative \e integers.
2.54 + ///
2.55 + /// \note For finding node-disjoint paths this algorithm can be used
2.56 + /// with \ref SplitDigraphAdaptor.
2.57 + ///
2.58 + /// \author Attila Bernath and Peter Kovacs
2.59 +
2.60 + template < typename Digraph,
2.61 + typename LengthMap = typename Digraph::template ArcMap<int> >
2.62 + class Suurballe
2.63 + {
2.64 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
2.65 +
2.66 + typedef typename LengthMap::Value Length;
2.67 + typedef ConstMap<Arc, int> ConstArcMap;
2.68 + typedef typename Digraph::template NodeMap<Arc> PredMap;
2.69 +
2.70 + public:
2.71 +
2.72 + /// The type of the flow map.
2.73 + typedef typename Digraph::template ArcMap<int> FlowMap;
2.74 + /// The type of the potential map.
2.75 + typedef typename Digraph::template NodeMap<Length> PotentialMap;
2.76 + /// The type of the path structures.
2.77 + typedef SimplePath<Digraph> Path;
2.78 +
2.79 + private:
2.80 +
2.81 + /// \brief Special implementation of the \ref Dijkstra algorithm
2.82 + /// for finding shortest paths in the residual network.
2.83 + ///
2.84 + /// \ref ResidualDijkstra is a special implementation of the
2.85 + /// \ref Dijkstra algorithm for finding shortest paths in the
2.86 + /// residual network of the digraph with respect to the reduced arc
2.87 + /// lengths and modifying the node potentials according to the
2.88 + /// distance of the nodes.
2.89 + class ResidualDijkstra
2.90 + {
2.91 + typedef typename Digraph::template NodeMap<int> HeapCrossRef;
2.92 + typedef BinHeap<Length, HeapCrossRef> Heap;
2.93 +
2.94 + private:
2.95 +
2.96 + // The directed digraph the algorithm runs on
2.97 + const Digraph &_graph;
2.98 +
2.99 + // The main maps
2.100 + const FlowMap &_flow;
2.101 + const LengthMap &_length;
2.102 + PotentialMap &_potential;
2.103 +
2.104 + // The distance map
2.105 + PotentialMap _dist;
2.106 + // The pred arc map
2.107 + PredMap &_pred;
2.108 + // The processed (i.e. permanently labeled) nodes
2.109 + std::vector<Node> _proc_nodes;
2.110 +
2.111 + Node _s;
2.112 + Node _t;
2.113 +
2.114 + public:
2.115 +
2.116 + /// Constructor.
2.117 + ResidualDijkstra( const Digraph &digraph,
2.118 + const FlowMap &flow,
2.119 + const LengthMap &length,
2.120 + PotentialMap &potential,
2.121 + PredMap &pred,
2.122 + Node s, Node t ) :
2.123 + _graph(digraph), _flow(flow), _length(length), _potential(potential),
2.124 + _dist(digraph), _pred(pred), _s(s), _t(t) {}
2.125 +
2.126 + /// \brief Runs the algorithm. Returns \c true if a path is found
2.127 + /// from the source node to the target node.
2.128 + bool run() {
2.129 + HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
2.130 + Heap heap(heap_cross_ref);
2.131 + heap.push(_s, 0);
2.132 + _pred[_s] = INVALID;
2.133 + _proc_nodes.clear();
2.134 +
2.135 + // Processing nodes
2.136 + while (!heap.empty() && heap.top() != _t) {
2.137 + Node u = heap.top(), v;
2.138 + Length d = heap.prio() + _potential[u], nd;
2.139 + _dist[u] = heap.prio();
2.140 + heap.pop();
2.141 + _proc_nodes.push_back(u);
2.142 +
2.143 + // Traversing outgoing arcs
2.144 + for (OutArcIt e(_graph, u); e != INVALID; ++e) {
2.145 + if (_flow[e] == 0) {
2.146 + v = _graph.target(e);
2.147 + switch(heap.state(v)) {
2.148 + case Heap::PRE_HEAP:
2.149 + heap.push(v, d + _length[e] - _potential[v]);
2.150 + _pred[v] = e;
2.151 + break;
2.152 + case Heap::IN_HEAP:
2.153 + nd = d + _length[e] - _potential[v];
2.154 + if (nd < heap[v]) {
2.155 + heap.decrease(v, nd);
2.156 + _pred[v] = e;
2.157 + }
2.158 + break;
2.159 + case Heap::POST_HEAP:
2.160 + break;
2.161 + }
2.162 + }
2.163 + }
2.164 +
2.165 + // Traversing incoming arcs
2.166 + for (InArcIt e(_graph, u); e != INVALID; ++e) {
2.167 + if (_flow[e] == 1) {
2.168 + v = _graph.source(e);
2.169 + switch(heap.state(v)) {
2.170 + case Heap::PRE_HEAP:
2.171 + heap.push(v, d - _length[e] - _potential[v]);
2.172 + _pred[v] = e;
2.173 + break;
2.174 + case Heap::IN_HEAP:
2.175 + nd = d - _length[e] - _potential[v];
2.176 + if (nd < heap[v]) {
2.177 + heap.decrease(v, nd);
2.178 + _pred[v] = e;
2.179 + }
2.180 + break;
2.181 + case Heap::POST_HEAP:
2.182 + break;
2.183 + }
2.184 + }
2.185 + }
2.186 + }
2.187 + if (heap.empty()) return false;
2.188 +
2.189 + // Updating potentials of processed nodes
2.190 + Length t_dist = heap.prio();
2.191 + for (int i = 0; i < int(_proc_nodes.size()); ++i)
2.192 + _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
2.193 + return true;
2.194 + }
2.195 +
2.196 + }; //class ResidualDijkstra
2.197 +
2.198 + private:
2.199 +
2.200 + // The directed digraph the algorithm runs on
2.201 + const Digraph &_graph;
2.202 + // The length map
2.203 + const LengthMap &_length;
2.204 +
2.205 + // Arc map of the current flow
2.206 + FlowMap *_flow;
2.207 + bool _local_flow;
2.208 + // Node map of the current potentials
2.209 + PotentialMap *_potential;
2.210 + bool _local_potential;
2.211 +
2.212 + // The source node
2.213 + Node _source;
2.214 + // The target node
2.215 + Node _target;
2.216 +
2.217 + // Container to store the found paths
2.218 + std::vector< SimplePath<Digraph> > paths;
2.219 + int _path_num;
2.220 +
2.221 + // The pred arc map
2.222 + PredMap _pred;
2.223 + // Implementation of the Dijkstra algorithm for finding augmenting
2.224 + // shortest paths in the residual network
2.225 + ResidualDijkstra *_dijkstra;
2.226 +
2.227 + public:
2.228 +
2.229 + /// \brief Constructor.
2.230 + ///
2.231 + /// Constructor.
2.232 + ///
2.233 + /// \param digraph The directed digraph the algorithm runs on.
2.234 + /// \param length The length (cost) values of the arcs.
2.235 + /// \param s The source node.
2.236 + /// \param t The target node.
2.237 + Suurballe( const Digraph &digraph,
2.238 + const LengthMap &length,
2.239 + Node s, Node t ) :
2.240 + _graph(digraph), _length(length), _flow(0), _local_flow(false),
2.241 + _potential(0), _local_potential(false), _source(s), _target(t),
2.242 + _pred(digraph) {}
2.243 +
2.244 + /// Destructor.
2.245 + ~Suurballe() {
2.246 + if (_local_flow) delete _flow;
2.247 + if (_local_potential) delete _potential;
2.248 + delete _dijkstra;
2.249 + }
2.250 +
2.251 + /// \brief Sets the flow map.
2.252 + ///
2.253 + /// Sets the flow map.
2.254 + ///
2.255 + /// The found flow contains only 0 and 1 values. It is the union of
2.256 + /// the found arc-disjoint paths.
2.257 + ///
2.258 + /// \return \c (*this)
2.259 + Suurballe& flowMap(FlowMap &map) {
2.260 + if (_local_flow) {
2.261 + delete _flow;
2.262 + _local_flow = false;
2.263 + }
2.264 + _flow = ↦
2.265 + return *this;
2.266 + }
2.267 +
2.268 + /// \brief Sets the potential map.
2.269 + ///
2.270 + /// Sets the potential map.
2.271 + ///
2.272 + /// The potentials provide the dual solution of the underlying
2.273 + /// minimum cost flow problem.
2.274 + ///
2.275 + /// \return \c (*this)
2.276 + Suurballe& potentialMap(PotentialMap &map) {
2.277 + if (_local_potential) {
2.278 + delete _potential;
2.279 + _local_potential = false;
2.280 + }
2.281 + _potential = ↦
2.282 + return *this;
2.283 + }
2.284 +
2.285 + /// \name Execution control
2.286 + /// The simplest way to execute the algorithm is to call the run()
2.287 + /// function.
2.288 + /// \n
2.289 + /// If you only need the flow that is the union of the found
2.290 + /// arc-disjoint paths, you may call init() and findFlow().
2.291 +
2.292 + /// @{
2.293 +
2.294 + /// \brief Runs the algorithm.
2.295 + ///
2.296 + /// Runs the algorithm.
2.297 + ///
2.298 + /// \param k The number of paths to be found.
2.299 + ///
2.300 + /// \return \c k if there are at least \c k arc-disjoint paths
2.301 + /// from \c s to \c t. Otherwise it returns the number of
2.302 + /// arc-disjoint paths found.
2.303 + ///
2.304 + /// \note Apart from the return value, <tt>s.run(k)</tt> is just a
2.305 + /// shortcut of the following code.
2.306 + /// \code
2.307 + /// s.init();
2.308 + /// s.findFlow(k);
2.309 + /// s.findPaths();
2.310 + /// \endcode
2.311 + int run(int k = 2) {
2.312 + init();
2.313 + findFlow(k);
2.314 + findPaths();
2.315 + return _path_num;
2.316 + }
2.317 +
2.318 + /// \brief Initializes the algorithm.
2.319 + ///
2.320 + /// Initializes the algorithm.
2.321 + void init() {
2.322 + // Initializing maps
2.323 + if (!_flow) {
2.324 + _flow = new FlowMap(_graph);
2.325 + _local_flow = true;
2.326 + }
2.327 + if (!_potential) {
2.328 + _potential = new PotentialMap(_graph);
2.329 + _local_potential = true;
2.330 + }
2.331 + for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
2.332 + for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
2.333 +
2.334 + _dijkstra = new ResidualDijkstra( _graph, *_flow, _length,
2.335 + *_potential, _pred,
2.336 + _source, _target );
2.337 + }
2.338 +
2.339 + /// \brief Executes the successive shortest path algorithm to find
2.340 + /// an optimal flow.
2.341 + ///
2.342 + /// Executes the successive shortest path algorithm to find a
2.343 + /// minimum cost flow, which is the union of \c k or less
2.344 + /// arc-disjoint paths.
2.345 + ///
2.346 + /// \return \c k if there are at least \c k arc-disjoint paths
2.347 + /// from \c s to \c t. Otherwise it returns the number of
2.348 + /// arc-disjoint paths found.
2.349 + ///
2.350 + /// \pre \ref init() must be called before using this function.
2.351 + int findFlow(int k = 2) {
2.352 + // Finding shortest paths
2.353 + _path_num = 0;
2.354 + while (_path_num < k) {
2.355 + // Running Dijkstra
2.356 + if (!_dijkstra->run()) break;
2.357 + ++_path_num;
2.358 +
2.359 + // Setting the flow along the found shortest path
2.360 + Node u = _target;
2.361 + Arc e;
2.362 + while ((e = _pred[u]) != INVALID) {
2.363 + if (u == _graph.target(e)) {
2.364 + (*_flow)[e] = 1;
2.365 + u = _graph.source(e);
2.366 + } else {
2.367 + (*_flow)[e] = 0;
2.368 + u = _graph.target(e);
2.369 + }
2.370 + }
2.371 + }
2.372 + return _path_num;
2.373 + }
2.374 +
2.375 + /// \brief Computes the paths from the flow.
2.376 + ///
2.377 + /// Computes the paths from the flow.
2.378 + ///
2.379 + /// \pre \ref init() and \ref findFlow() must be called before using
2.380 + /// this function.
2.381 + void findPaths() {
2.382 + // Creating the residual flow map (the union of the paths not
2.383 + // found so far)
2.384 + FlowMap res_flow(_graph);
2.385 + for(ArcIt a(_graph);a!=INVALID;++a) res_flow[a]=(*_flow)[a];
2.386 +
2.387 + paths.clear();
2.388 + paths.resize(_path_num);
2.389 + for (int i = 0; i < _path_num; ++i) {
2.390 + Node n = _source;
2.391 + while (n != _target) {
2.392 + OutArcIt e(_graph, n);
2.393 + for ( ; res_flow[e] == 0; ++e) ;
2.394 + n = _graph.target(e);
2.395 + paths[i].addBack(e);
2.396 + res_flow[e] = 0;
2.397 + }
2.398 + }
2.399 + }
2.400 +
2.401 + /// @}
2.402 +
2.403 + /// \name Query Functions
2.404 + /// The result of the algorithm can be obtained using these
2.405 + /// functions.
2.406 + /// \n The algorithm should be executed before using them.
2.407 +
2.408 + /// @{
2.409 +
2.410 + /// \brief Returns a const reference to the arc map storing the
2.411 + /// found flow.
2.412 + ///
2.413 + /// Returns a const reference to the arc map storing the flow that
2.414 + /// is the union of the found arc-disjoint paths.
2.415 + ///
2.416 + /// \pre \ref run() or findFlow() must be called before using this
2.417 + /// function.
2.418 + const FlowMap& flowMap() const {
2.419 + return *_flow;
2.420 + }
2.421 +
2.422 + /// \brief Returns a const reference to the node map storing the
2.423 + /// found potentials (the dual solution).
2.424 + ///
2.425 + /// Returns a const reference to the node map storing the found
2.426 + /// potentials that provide the dual solution of the underlying
2.427 + /// minimum cost flow problem.
2.428 + ///
2.429 + /// \pre \ref run() or findFlow() must be called before using this
2.430 + /// function.
2.431 + const PotentialMap& potentialMap() const {
2.432 + return *_potential;
2.433 + }
2.434 +
2.435 + /// \brief Returns the flow on the given arc.
2.436 + ///
2.437 + /// Returns the flow on the given arc.
2.438 + /// It is \c 1 if the arc is involved in one of the found paths,
2.439 + /// otherwise it is \c 0.
2.440 + ///
2.441 + /// \pre \ref run() or findFlow() must be called before using this
2.442 + /// function.
2.443 + int flow(const Arc& arc) const {
2.444 + return (*_flow)[arc];
2.445 + }
2.446 +
2.447 + /// \brief Returns the potential of the given node.
2.448 + ///
2.449 + /// Returns the potential of the given node.
2.450 + ///
2.451 + /// \pre \ref run() or findFlow() must be called before using this
2.452 + /// function.
2.453 + Length potential(const Node& node) const {
2.454 + return (*_potential)[node];
2.455 + }
2.456 +
2.457 + /// \brief Returns the total length (cost) of the found paths (flow).
2.458 + ///
2.459 + /// Returns the total length (cost) of the found paths (flow).
2.460 + /// The complexity of the function is \f$ O(e) \f$.
2.461 + ///
2.462 + /// \pre \ref run() or findFlow() must be called before using this
2.463 + /// function.
2.464 + Length totalLength() const {
2.465 + Length c = 0;
2.466 + for (ArcIt e(_graph); e != INVALID; ++e)
2.467 + c += (*_flow)[e] * _length[e];
2.468 + return c;
2.469 + }
2.470 +
2.471 + /// \brief Returns the number of the found paths.
2.472 + ///
2.473 + /// Returns the number of the found paths.
2.474 + ///
2.475 + /// \pre \ref run() or findFlow() must be called before using this
2.476 + /// function.
2.477 + int pathNum() const {
2.478 + return _path_num;
2.479 + }
2.480 +
2.481 + /// \brief Returns a const reference to the specified path.
2.482 + ///
2.483 + /// Returns a const reference to the specified path.
2.484 + ///
2.485 + /// \param i The function returns the \c i-th path.
2.486 + /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>.
2.487 + ///
2.488 + /// \pre \ref run() or findPaths() must be called before using this
2.489 + /// function.
2.490 + Path path(int i) const {
2.491 + return paths[i];
2.492 + }
2.493 +
2.494 + /// @}
2.495 +
2.496 + }; //class Suurballe
2.497 +
2.498 + ///@}
2.499 +
2.500 +} //namespace lemon
2.501 +
2.502 +#endif //LEMON_SUURBALLE_H
3.1 --- a/test/Makefile.am Tue Oct 28 18:33:51 2008 +0100
3.2 +++ b/test/Makefile.am Tue Oct 28 18:39:53 2008 +0000
3.3 @@ -1,5 +1,6 @@
3.4 EXTRA_DIST += \
3.5 - test/CMakeLists.txt
3.6 + test/CMakeLists.txt \
3.7 + test/min_cost_flow_test.lgf
3.8
3.9 noinst_HEADERS += \
3.10 test/graph_test.h \
3.11 @@ -22,6 +23,7 @@
3.12 test/max_matching_test \
3.13 test/random_test \
3.14 test/path_test \
3.15 + test/suurballe_test \
3.16 test/test_tools_fail \
3.17 test/test_tools_pass \
3.18 test/time_measure_test \
3.19 @@ -45,6 +47,7 @@
3.20 test_maps_test_SOURCES = test/maps_test.cc
3.21 test_max_matching_test_SOURCES = test/max_matching_test.cc
3.22 test_path_test_SOURCES = test/path_test.cc
3.23 +test_suurballe_test_SOURCES = test/suurballe_test.cc
3.24 test_random_test_SOURCES = test/random_test.cc
3.25 test_test_tools_fail_SOURCES = test/test_tools_fail.cc
3.26 test_test_tools_pass_SOURCES = test/test_tools_pass.cc
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/test/min_cost_flow_test.lgf Tue Oct 28 18:39:53 2008 +0000
4.3 @@ -0,0 +1,44 @@
4.4 +@nodes
4.5 +label supply1 supply2 supply3
4.6 +1 0 20 27
4.7 +2 0 -4 0
4.8 +3 0 0 0
4.9 +4 0 0 0
4.10 +5 0 9 0
4.11 +6 0 -6 0
4.12 +7 0 0 0
4.13 +8 0 0 0
4.14 +9 0 3 0
4.15 +10 0 -2 0
4.16 +11 0 0 0
4.17 +12 0 -20 -27
4.18 +
4.19 +@arcs
4.20 + cost capacity lower1 lower2
4.21 +1 2 70 11 0 8
4.22 +1 3 150 3 0 1
4.23 +1 4 80 15 0 2
4.24 +2 8 80 12 0 0
4.25 +3 5 140 5 0 3
4.26 +4 6 60 10 0 1
4.27 +4 7 80 2 0 0
4.28 +4 8 110 3 0 0
4.29 +5 7 60 14 0 0
4.30 +5 11 120 12 0 0
4.31 +6 3 0 3 0 0
4.32 +6 9 140 4 0 0
4.33 +6 10 90 8 0 0
4.34 +7 1 30 5 0 0
4.35 +8 12 60 16 0 4
4.36 +9 12 50 6 0 0
4.37 +10 12 70 13 0 5
4.38 +10 2 100 7 0 0
4.39 +10 7 60 10 0 0
4.40 +11 10 20 14 0 6
4.41 +12 11 30 10 0 0
4.42 +
4.43 +@attributes
4.44 +source 1
4.45 +target 12
4.46 +
4.47 +@end
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/test/suurballe_test.cc Tue Oct 28 18:39:53 2008 +0000
5.3 @@ -0,0 +1,160 @@
5.4 +/* -*- C++ -*-
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 <iostream>
5.23 +#include <fstream>
5.24 +
5.25 +#include <lemon/list_graph.h>
5.26 +#include <lemon/lgf_reader.h>
5.27 +#include <lemon/path.h>
5.28 +#include <lemon/suurballe.h>
5.29 +
5.30 +#include "test_tools.h"
5.31 +
5.32 +using namespace lemon;
5.33 +
5.34 +// Checks the feasibility of the flow
5.35 +template <typename Digraph, typename FlowMap>
5.36 +bool checkFlow( const Digraph& gr, const FlowMap& flow,
5.37 + typename Digraph::Node s, typename Digraph::Node t,
5.38 + int value )
5.39 +{
5.40 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
5.41 + for (ArcIt e(gr); e != INVALID; ++e)
5.42 + if (!(flow[e] == 0 || flow[e] == 1)) return false;
5.43 +
5.44 + for (NodeIt n(gr); n != INVALID; ++n) {
5.45 + int sum = 0;
5.46 + for (OutArcIt e(gr, n); e != INVALID; ++e)
5.47 + sum += flow[e];
5.48 + for (InArcIt e(gr, n); e != INVALID; ++e)
5.49 + sum -= flow[e];
5.50 + if (n == s && sum != value) return false;
5.51 + if (n == t && sum != -value) return false;
5.52 + if (n != s && n != t && sum != 0) return false;
5.53 + }
5.54 +
5.55 + return true;
5.56 +}
5.57 +
5.58 +// Checks the optimalitiy of the flow
5.59 +template < typename Digraph, typename CostMap,
5.60 + typename FlowMap, typename PotentialMap >
5.61 +bool checkOptimality( const Digraph& gr, const CostMap& cost,
5.62 + const FlowMap& flow, const PotentialMap& pi )
5.63 +{
5.64 + // Checking the Complementary Slackness optimality condition
5.65 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
5.66 + bool opt = true;
5.67 + for (ArcIt e(gr); e != INVALID; ++e) {
5.68 + typename CostMap::Value red_cost =
5.69 + cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
5.70 + opt = (flow[e] == 0 && red_cost >= 0) ||
5.71 + (flow[e] == 1 && red_cost <= 0);
5.72 + if (!opt) break;
5.73 + }
5.74 + return opt;
5.75 +}
5.76 +
5.77 +// Checks a path
5.78 +template < typename Digraph, typename Path >
5.79 +bool checkPath( const Digraph& gr, const Path& path,
5.80 + typename Digraph::Node s, typename Digraph::Node t)
5.81 +{
5.82 + // Checking the Complementary Slackness optimality condition
5.83 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
5.84 + Node n = s;
5.85 + for (int i = 0; i < path.length(); ++i) {
5.86 + if (gr.source(path.nth(i)) != n) return false;
5.87 + n = gr.target(path.nth(i));
5.88 + }
5.89 + return n == t;
5.90 +}
5.91 +
5.92 +
5.93 +int main()
5.94 +{
5.95 + DIGRAPH_TYPEDEFS(ListDigraph);
5.96 +
5.97 + // Reading the test digraph
5.98 + ListDigraph digraph;
5.99 + ListDigraph::ArcMap<int> length(digraph);
5.100 + Node source, target;
5.101 +
5.102 + std::string fname;
5.103 + if(getenv("srcdir"))
5.104 + fname = std::string(getenv("srcdir"));
5.105 + else fname = ".";
5.106 + fname += "/test/min_cost_flow_test.lgf";
5.107 +
5.108 + std::ifstream input(fname.c_str());
5.109 + check(input, "Input file '" << fname << "' not found");
5.110 + DigraphReader<ListDigraph>(digraph, input).
5.111 + arcMap("cost", length).
5.112 + node("source", source).
5.113 + node("target", target).
5.114 + run();
5.115 + input.close();
5.116 +
5.117 + // Finding 2 paths
5.118 + {
5.119 + Suurballe<ListDigraph> suurballe(digraph, length, source, target);
5.120 + check(suurballe.run(2) == 2, "Wrong number of paths");
5.121 + check(checkFlow(digraph, suurballe.flowMap(), source, target, 2),
5.122 + "The flow is not feasible");
5.123 + check(suurballe.totalLength() == 510, "The flow is not optimal");
5.124 + check(checkOptimality(digraph, length, suurballe.flowMap(),
5.125 + suurballe.potentialMap()),
5.126 + "Wrong potentials");
5.127 + for (int i = 0; i < suurballe.pathNum(); ++i)
5.128 + check(checkPath(digraph, suurballe.path(i), source, target),
5.129 + "Wrong path");
5.130 + }
5.131 +
5.132 + // Finding 3 paths
5.133 + {
5.134 + Suurballe<ListDigraph> suurballe(digraph, length, source, target);
5.135 + check(suurballe.run(3) == 3, "Wrong number of paths");
5.136 + check(checkFlow(digraph, suurballe.flowMap(), source, target, 3),
5.137 + "The flow is not feasible");
5.138 + check(suurballe.totalLength() == 1040, "The flow is not optimal");
5.139 + check(checkOptimality(digraph, length, suurballe.flowMap(),
5.140 + suurballe.potentialMap()),
5.141 + "Wrong potentials");
5.142 + for (int i = 0; i < suurballe.pathNum(); ++i)
5.143 + check(checkPath(digraph, suurballe.path(i), source, target),
5.144 + "Wrong path");
5.145 + }
5.146 +
5.147 + // Finding 5 paths (only 3 can be found)
5.148 + {
5.149 + Suurballe<ListDigraph> suurballe(digraph, length, source, target);
5.150 + check(suurballe.run(5) == 3, "Wrong number of paths");
5.151 + check(checkFlow(digraph, suurballe.flowMap(), source, target, 3),
5.152 + "The flow is not feasible");
5.153 + check(suurballe.totalLength() == 1040, "The flow is not optimal");
5.154 + check(checkOptimality(digraph, length, suurballe.flowMap(),
5.155 + suurballe.potentialMap()),
5.156 + "Wrong potentials");
5.157 + for (int i = 0; i < suurballe.pathNum(); ++i)
5.158 + check(checkPath(digraph, suurballe.path(i), source, target),
5.159 + "Wrong path");
5.160 + }
5.161 +
5.162 + return 0;
5.163 +}