Reimplemented Suurballe class.
- The new version is the specialized version of CapacityScaling.
- It is about 10-20 times faster than the former Suurballe algorithm
and about 20-50 percent faster than CapacityScaling.
- Doc improvements.
- The test file is also replaced.
1.1 --- a/lemon/suurballe.h Thu Feb 28 16:33:40 2008 +0000
1.2 +++ b/lemon/suurballe.h Fri Feb 29 15:55:13 2008 +0000
1.3 @@ -21,182 +21,474 @@
1.4
1.5 ///\ingroup shortest_path
1.6 ///\file
1.7 -///\brief An algorithm for finding k paths of minimal total length.
1.8 +///\brief An algorithm for finding edge-disjoint paths between two
1.9 +/// nodes having minimum total length.
1.10
1.11 -
1.12 -#include <lemon/maps.h>
1.13 #include <vector>
1.14 +#include <lemon/bin_heap.h>
1.15 #include <lemon/path.h>
1.16 -#include <lemon/ssp_min_cost_flow.h>
1.17
1.18 namespace lemon {
1.19
1.20 -/// \addtogroup shortest_path
1.21 -/// @{
1.22 + /// \addtogroup shortest_path
1.23 + /// @{
1.24
1.25 - ///\brief Implementation of an algorithm for finding k edge-disjoint
1.26 - /// paths between 2 nodes of minimal total length
1.27 + /// \brief Implementation of an algorithm for finding edge-disjoint
1.28 + /// paths between two nodes having minimum total length.
1.29 ///
1.30 - /// The class \ref lemon::Suurballe implements
1.31 - /// an algorithm for finding k edge-disjoint paths
1.32 - /// from a given source node to a given target node in an
1.33 - /// edge-weighted directed graph having minimal total weight (length).
1.34 + /// \ref lemon::Suurballe "Suurballe" implements an algorithm for
1.35 + /// finding edge-disjoint paths having minimum total length (cost)
1.36 + /// from a given source node to a given target node in a directed
1.37 + /// graph.
1.38 ///
1.39 - ///\warning Length values should be nonnegative!
1.40 - ///
1.41 - ///\param Graph The directed graph type the algorithm runs on.
1.42 - ///\param LengthMap The type of the length map (values should be nonnegative).
1.43 + /// In fact, this implementation is the specialization of the
1.44 + /// \ref CapacityScaling "successive shortest path" algorithm.
1.45 ///
1.46 - ///\note It it questionable whether it is correct to call this method after
1.47 - ///%Suurballe for it is just a special case of Edmonds' and Karp's algorithm
1.48 - ///for finding minimum cost flows. In fact, this implementation just
1.49 - ///wraps the SspMinCostFlow algorithms. The paper of both %Suurballe and
1.50 - ///Edmonds-Karp published in 1972, therefore it is possibly right to
1.51 - ///state that they are
1.52 - ///independent results. Most frequently this special case is referred as
1.53 - ///%Suurballe method in the literature, especially in communication
1.54 - ///network context.
1.55 - ///\author Attila Bernath
1.56 - template <typename Graph, typename LengthMap>
1.57 - class Suurballe{
1.58 -
1.59 + /// \tparam Graph The directed graph type the algorithm runs on.
1.60 + /// \tparam LengthMap The type of the length (cost) map.
1.61 + ///
1.62 + /// \warning Length values should be \e non-negative \e integers.
1.63 + ///
1.64 + /// \note For finding node-disjoint paths this algorithm can be used
1.65 + /// with \ref SplitGraphAdaptor.
1.66 + ///
1.67 + /// \author Attila Bernath and Peter Kovacs
1.68 +
1.69 + template < typename Graph,
1.70 + typename LengthMap = typename Graph::template EdgeMap<int> >
1.71 + class Suurballe
1.72 + {
1.73 + GRAPH_TYPEDEFS(typename Graph);
1.74
1.75 typedef typename LengthMap::Value Length;
1.76 + typedef ConstMap<Edge, int> ConstEdgeMap;
1.77 + typedef typename Graph::template NodeMap<Edge> PredMap;
1.78 +
1.79 + public:
1.80 +
1.81 + /// The type of the flow map.
1.82 + typedef typename Graph::template EdgeMap<int> FlowMap;
1.83 + /// The type of the potential map.
1.84 + typedef typename Graph::template NodeMap<Length> PotentialMap;
1.85 + /// The type of the path structures.
1.86 + typedef SimplePath<Graph> Path;
1.87 +
1.88 + private:
1.89 +
1.90 + /// \brief Special implementation of the \ref Dijkstra algorithm
1.91 + /// for finding shortest paths in the residual network.
1.92 + ///
1.93 + /// \ref ResidualDijkstra is a special implementation of the
1.94 + /// \ref Dijkstra algorithm for finding shortest paths in the
1.95 + /// residual network of the graph with respect to the reduced edge
1.96 + /// lengths and modifying the node potentials according to the
1.97 + /// distance of the nodes.
1.98 + class ResidualDijkstra
1.99 + {
1.100 + typedef typename Graph::template NodeMap<int> HeapCrossRef;
1.101 + typedef BinHeap<Length, HeapCrossRef> Heap;
1.102 +
1.103 + private:
1.104 +
1.105 + // The directed graph the algorithm runs on
1.106 + const Graph &_graph;
1.107 +
1.108 + // The main maps
1.109 + const FlowMap &_flow;
1.110 + const LengthMap &_length;
1.111 + PotentialMap &_potential;
1.112 +
1.113 + // The distance map
1.114 + PotentialMap _dist;
1.115 + // The pred edge map
1.116 + PredMap &_pred;
1.117 + // The processed (i.e. permanently labeled) nodes
1.118 + std::vector<Node> _proc_nodes;
1.119 +
1.120 + Node _s;
1.121 + Node _t;
1.122 +
1.123 + public:
1.124 +
1.125 + /// Constructor.
1.126 + ResidualDijkstra( const Graph &graph,
1.127 + const FlowMap &flow,
1.128 + const LengthMap &length,
1.129 + PotentialMap &potential,
1.130 + PredMap &pred,
1.131 + Node s, Node t ) :
1.132 + _graph(graph), _flow(flow), _length(length), _potential(potential),
1.133 + _dist(graph), _pred(pred), _s(s), _t(t) {}
1.134 +
1.135 + /// \brief Runs the algorithm. Returns \c true if a path is found
1.136 + /// from the source node to the target node.
1.137 + bool run() {
1.138 + HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
1.139 + Heap heap(heap_cross_ref);
1.140 + heap.push(_s, 0);
1.141 + _pred[_s] = INVALID;
1.142 + _proc_nodes.clear();
1.143 +
1.144 + // Processing nodes
1.145 + while (!heap.empty() && heap.top() != _t) {
1.146 + Node u = heap.top(), v;
1.147 + Length d = heap.prio() + _potential[u], nd;
1.148 + _dist[u] = heap.prio();
1.149 + heap.pop();
1.150 + _proc_nodes.push_back(u);
1.151 +
1.152 + // Traversing outgoing edges
1.153 + for (OutEdgeIt e(_graph, u); e != INVALID; ++e) {
1.154 + if (_flow[e] == 0) {
1.155 + v = _graph.target(e);
1.156 + switch(heap.state(v)) {
1.157 + case Heap::PRE_HEAP:
1.158 + heap.push(v, d + _length[e] - _potential[v]);
1.159 + _pred[v] = e;
1.160 + break;
1.161 + case Heap::IN_HEAP:
1.162 + nd = d + _length[e] - _potential[v];
1.163 + if (nd < heap[v]) {
1.164 + heap.decrease(v, nd);
1.165 + _pred[v] = e;
1.166 + }
1.167 + break;
1.168 + case Heap::POST_HEAP:
1.169 + break;
1.170 + }
1.171 + }
1.172 + }
1.173 +
1.174 + // Traversing incoming edges
1.175 + for (InEdgeIt e(_graph, u); e != INVALID; ++e) {
1.176 + if (_flow[e] == 1) {
1.177 + v = _graph.source(e);
1.178 + switch(heap.state(v)) {
1.179 + case Heap::PRE_HEAP:
1.180 + heap.push(v, d - _length[e] - _potential[v]);
1.181 + _pred[v] = e;
1.182 + break;
1.183 + case Heap::IN_HEAP:
1.184 + nd = d - _length[e] - _potential[v];
1.185 + if (nd < heap[v]) {
1.186 + heap.decrease(v, nd);
1.187 + _pred[v] = e;
1.188 + }
1.189 + break;
1.190 + case Heap::POST_HEAP:
1.191 + break;
1.192 + }
1.193 + }
1.194 + }
1.195 + }
1.196 + if (heap.empty()) return false;
1.197 +
1.198 + // Updating potentials of processed nodes
1.199 + Length t_dist = heap.prio();
1.200 + for (int i = 0; i < int(_proc_nodes.size()); ++i)
1.201 + _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
1.202 + return true;
1.203 + }
1.204 +
1.205 + }; //class ResidualDijkstra
1.206 +
1.207 + private:
1.208 +
1.209 + // The directed graph the algorithm runs on
1.210 + const Graph &_graph;
1.211 + // The length map
1.212 + const LengthMap &_length;
1.213
1.214 - typedef typename Graph::Node Node;
1.215 - typedef typename Graph::NodeIt NodeIt;
1.216 - typedef typename Graph::Edge Edge;
1.217 - typedef typename Graph::OutEdgeIt OutEdgeIt;
1.218 - typedef typename Graph::template EdgeMap<int> EdgeIntMap;
1.219 + // Edge map of the current flow
1.220 + FlowMap *_flow;
1.221 + bool _local_flow;
1.222 + // Node map of the current potentials
1.223 + PotentialMap *_potential;
1.224 + bool _local_potential;
1.225
1.226 - typedef ConstMap<Edge,int> ConstMap;
1.227 + // The source node
1.228 + Node _source;
1.229 + // The target node
1.230 + Node _target;
1.231
1.232 - const Graph& G;
1.233 + // Container to store the found paths
1.234 + std::vector< SimplePath<Graph> > paths;
1.235 + int _path_num;
1.236
1.237 - Node s;
1.238 - Node t;
1.239 + // The pred edge map
1.240 + PredMap _pred;
1.241 + // Implementation of the Dijkstra algorithm for finding augmenting
1.242 + // shortest paths in the residual network
1.243 + ResidualDijkstra *_dijkstra;
1.244
1.245 - //Auxiliary variables
1.246 - //This is the capacity map for the mincostflow problem
1.247 - ConstMap const1map;
1.248 - //This MinCostFlow instance will actually solve the problem
1.249 - SspMinCostFlow<Graph, LengthMap, ConstMap> min_cost_flow;
1.250 + public:
1.251
1.252 - //Container to store found paths
1.253 - std::vector<SimplePath<Graph> > paths;
1.254 + /// \brief Constructor.
1.255 + ///
1.256 + /// Constructor.
1.257 + ///
1.258 + /// \param graph The directed graph the algorithm runs on.
1.259 + /// \param length The length (cost) values of the edges.
1.260 + /// \param s The source node.
1.261 + /// \param t The target node.
1.262 + Suurballe( const Graph &graph,
1.263 + const LengthMap &length,
1.264 + Node s, Node t ) :
1.265 + _graph(graph), _length(length), _flow(0), _local_flow(false),
1.266 + _potential(0), _local_potential(false), _source(s), _target(t),
1.267 + _pred(graph) {}
1.268
1.269 - public :
1.270 + /// Destructor.
1.271 + ~Suurballe() {
1.272 + if (_local_flow) delete _flow;
1.273 + if (_local_potential) delete _potential;
1.274 + delete _dijkstra;
1.275 + }
1.276
1.277 + /// \brief Sets the flow map.
1.278 + ///
1.279 + /// Sets the flow map.
1.280 + ///
1.281 + /// The found flow contains only 0 and 1 values. It is the union of
1.282 + /// the found edge-disjoint paths.
1.283 + ///
1.284 + /// \return \c (*this)
1.285 + Suurballe& flowMap(FlowMap &map) {
1.286 + if (_local_flow) {
1.287 + delete _flow;
1.288 + _local_flow = false;
1.289 + }
1.290 + _flow = ↦
1.291 + return *this;
1.292 + }
1.293
1.294 - /// \brief The constructor of the class.
1.295 + /// \brief Sets the potential map.
1.296 ///
1.297 - /// \param _G The directed graph the algorithm runs on.
1.298 - /// \param _length The length (weight or cost) of the edges.
1.299 - /// \param _s Source node.
1.300 - /// \param _t Target node.
1.301 - Suurballe(Graph& _G, LengthMap& _length, Node _s, Node _t) :
1.302 - G(_G), s(_s), t(_t), const1map(1),
1.303 - min_cost_flow(_G, _length, const1map, _s, _t) { }
1.304 + /// Sets the potential map.
1.305 + ///
1.306 + /// The potentials provide the dual solution of the underlying
1.307 + /// minimum cost flow problem.
1.308 + ///
1.309 + /// \return \c (*this)
1.310 + Suurballe& potentialMap(PotentialMap &map) {
1.311 + if (_local_potential) {
1.312 + delete _potential;
1.313 + _local_potential = false;
1.314 + }
1.315 + _potential = ↦
1.316 + return *this;
1.317 + }
1.318 +
1.319 + /// \name Execution control
1.320 + /// The simplest way to execute the algorithm is to call the run()
1.321 + /// function.
1.322 + /// \n
1.323 + /// If you only need the flow that is the union of the found
1.324 + /// edge-disjoint paths, you may call init() and findFlow().
1.325 +
1.326 + /// @{
1.327
1.328 /// \brief Runs the algorithm.
1.329 ///
1.330 /// Runs the algorithm.
1.331 - /// Returns k if there are at least k edge-disjoint paths from s to t.
1.332 - /// Otherwise it returns the number of edge-disjoint paths found
1.333 - /// from s to t.
1.334 ///
1.335 - /// \param k How many paths are we looking for?
1.336 + /// \param k The number of paths to be found.
1.337 ///
1.338 - int run(int k) {
1.339 - int i = min_cost_flow.run(k);
1.340 -
1.341 - //Let's find the paths
1.342 - //We put the paths into stl vectors (as an inner representation).
1.343 - //In the meantime we lose the information stored in 'reversed'.
1.344 - //We suppose the lengths to be positive now.
1.345 -
1.346 - //We don't want to change the flow of min_cost_flow, so we make a copy
1.347 - //The name here suggests that the flow has only 0/1 values.
1.348 - EdgeIntMap reversed(G);
1.349 -
1.350 - for(typename Graph::EdgeIt e(G); e!=INVALID; ++e)
1.351 - reversed[e] = min_cost_flow.getFlow()[e];
1.352 -
1.353 - paths.clear();
1.354 - paths.resize(k);
1.355 - for (int j=0; j<i; ++j){
1.356 - Node n=s;
1.357 -
1.358 - while (n!=t){
1.359 -
1.360 - OutEdgeIt e(G, n);
1.361 -
1.362 - while (!reversed[e]){
1.363 - ++e;
1.364 - }
1.365 - n = G.target(e);
1.366 - paths[j].addBack(e);
1.367 - reversed[e] = 1-reversed[e];
1.368 - }
1.369 -
1.370 - }
1.371 - return i;
1.372 + /// \return \c k if there are at least \c k edge-disjoint paths
1.373 + /// from \c s to \c t. Otherwise it returns the number of
1.374 + /// edge-disjoint paths found.
1.375 + ///
1.376 + /// \note Apart from the return value, <tt>s.run(k)</tt> is just a
1.377 + /// shortcut of the following code.
1.378 + /// \code
1.379 + /// s.init();
1.380 + /// s.findFlow(k);
1.381 + /// s.findPaths();
1.382 + /// \endcode
1.383 + int run(int k = 2) {
1.384 + init();
1.385 + findFlow(k);
1.386 + findPaths();
1.387 + return _path_num;
1.388 }
1.389
1.390 -
1.391 - /// \brief Returns the total length of the paths.
1.392 + /// \brief Initializes the algorithm.
1.393 ///
1.394 - /// This function gives back the total length of the found paths.
1.395 - Length totalLength(){
1.396 - return min_cost_flow.totalLength();
1.397 + /// Initializes the algorithm.
1.398 + void init() {
1.399 + // Initializing maps
1.400 + if (!_flow) {
1.401 + _flow = new FlowMap(_graph);
1.402 + _local_flow = true;
1.403 + }
1.404 + if (!_potential) {
1.405 + _potential = new PotentialMap(_graph);
1.406 + _local_potential = true;
1.407 + }
1.408 + for (EdgeIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
1.409 + for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
1.410 +
1.411 + _dijkstra = new ResidualDijkstra( _graph, *_flow, _length,
1.412 + *_potential, _pred,
1.413 + _source, _target );
1.414 }
1.415
1.416 - /// \brief Returns the found flow.
1.417 + /// \brief Executes the successive shortest path algorithm to find
1.418 + /// an optimal flow.
1.419 ///
1.420 - /// This function returns a const reference to the EdgeMap \c flow.
1.421 - const EdgeIntMap &getFlow() const { return min_cost_flow.flow;}
1.422 + /// Executes the successive shortest path algorithm to find a
1.423 + /// minimum cost flow, which is the union of \c k or less
1.424 + /// edge-disjoint paths.
1.425 + ///
1.426 + /// \return \c k if there are at least \c k edge-disjoint paths
1.427 + /// from \c s to \c t. Otherwise it returns the number of
1.428 + /// edge-disjoint paths found.
1.429 + ///
1.430 + /// \pre \ref init() must be called before using this function.
1.431 + int findFlow(int k = 2) {
1.432 + // Finding shortest paths
1.433 + _path_num = 0;
1.434 + while (_path_num < k) {
1.435 + // Running Dijkstra
1.436 + if (!_dijkstra->run()) break;
1.437 + ++_path_num;
1.438
1.439 - /// \brief Returns the optimal dual solution
1.440 + // Setting the flow along the found shortest path
1.441 + Node u = _target;
1.442 + Edge e;
1.443 + while ((e = _pred[u]) != INVALID) {
1.444 + if (u == _graph.target(e)) {
1.445 + (*_flow)[e] = 1;
1.446 + u = _graph.source(e);
1.447 + } else {
1.448 + (*_flow)[e] = 0;
1.449 + u = _graph.target(e);
1.450 + }
1.451 + }
1.452 + }
1.453 + return _path_num;
1.454 + }
1.455 +
1.456 + /// \brief Computes the paths from the flow.
1.457 ///
1.458 - /// This function returns a const reference to the NodeMap \c
1.459 - /// potential (the dual solution).
1.460 - const EdgeIntMap &getPotential() const { return min_cost_flow.potential;}
1.461 + /// Computes the paths from the flow.
1.462 + ///
1.463 + /// \pre \ref init() and \ref findFlow() must be called before using
1.464 + /// this function.
1.465 + void findPaths() {
1.466 + // Creating the residual flow map (the union of the paths not
1.467 + // found so far)
1.468 + FlowMap res_flow(*_flow);
1.469
1.470 - /// \brief Checks whether the complementary slackness holds.
1.471 - ///
1.472 - /// This function checks, whether the given solution is optimal.
1.473 - /// Currently this function only checks optimality, doesn't bother
1.474 - /// with feasibility. It is meant for testing purposes.
1.475 - bool checkComplementarySlackness(){
1.476 - return min_cost_flow.checkComplementarySlackness();
1.477 + paths.clear();
1.478 + paths.resize(_path_num);
1.479 + for (int i = 0; i < _path_num; ++i) {
1.480 + Node n = _source;
1.481 + while (n != _target) {
1.482 + OutEdgeIt e(_graph, n);
1.483 + for ( ; res_flow[e] == 0; ++e) ;
1.484 + n = _graph.target(e);
1.485 + paths[i].addBack(e);
1.486 + res_flow[e] = 0;
1.487 + }
1.488 + }
1.489 }
1.490
1.491 - typedef SimplePath<Graph> Path;
1.492 + /// @}
1.493
1.494 - /// \brief Read the found paths.
1.495 + /// \name Query Functions
1.496 + /// The result of the algorithm can be obtained using these
1.497 + /// functions.
1.498 + /// \n The algorithm should be executed before using them.
1.499 +
1.500 + /// @{
1.501 +
1.502 + /// \brief Returns a const reference to the edge map storing the
1.503 + /// found flow.
1.504 ///
1.505 - /// This function gives back the \c j-th path in argument p.
1.506 - /// Assumes that \c run() has been run and nothing has changed
1.507 - /// since then.
1.508 + /// Returns a const reference to the edge map storing the flow that
1.509 + /// is the union of the found edge-disjoint paths.
1.510 ///
1.511 - /// \warning It is assumed that \c p is constructed to be a path
1.512 - /// of graph \c G. If \c j is not less than the result of
1.513 - /// previous \c run, then the result here will be an empty path
1.514 - /// (\c j can be 0 as well).
1.515 - ///
1.516 - /// \param j Which path you want to get from the found paths (in a
1.517 - /// real application you would get the found paths iteratively).
1.518 - Path path(int j) const {
1.519 - return paths[j];
1.520 + /// \pre \ref run() or findFlow() must be called before using this
1.521 + /// function.
1.522 + const FlowMap& flowMap() const {
1.523 + return *_flow;
1.524 }
1.525
1.526 - /// \brief Gives back the number of the paths.
1.527 + /// \brief Returns a const reference to the node map storing the
1.528 + /// found potentials (the dual solution).
1.529 ///
1.530 - /// Gives back the number of the constructed paths.
1.531 + /// Returns a const reference to the node map storing the found
1.532 + /// potentials that provide the dual solution of the underlying
1.533 + /// minimum cost flow problem.
1.534 + ///
1.535 + /// \pre \ref run() or findFlow() must be called before using this
1.536 + /// function.
1.537 + const PotentialMap& potentialMap() const {
1.538 + return *_potential;
1.539 + }
1.540 +
1.541 + /// \brief Returns the flow on the given edge.
1.542 + ///
1.543 + /// Returns the flow on the given edge.
1.544 + /// It is \c 1 if the edge is involved in one of the found paths,
1.545 + /// otherwise it is \c 0.
1.546 + ///
1.547 + /// \pre \ref run() or findFlow() must be called before using this
1.548 + /// function.
1.549 + int flow(const Edge& edge) const {
1.550 + return (*_flow)[edge];
1.551 + }
1.552 +
1.553 + /// \brief Returns the potential of the given node.
1.554 + ///
1.555 + /// Returns the potential of the given node.
1.556 + ///
1.557 + /// \pre \ref run() or findFlow() must be called before using this
1.558 + /// function.
1.559 + Length potential(const Node& node) const {
1.560 + return (*_potential)[node];
1.561 + }
1.562 +
1.563 + /// \brief Returns the total length (cost) of the found paths (flow).
1.564 + ///
1.565 + /// Returns the total length (cost) of the found paths (flow).
1.566 + /// The complexity of the function is \f$ O(e) \f$.
1.567 + ///
1.568 + /// \pre \ref run() or findFlow() must be called before using this
1.569 + /// function.
1.570 + Length totalLength() const {
1.571 + Length c = 0;
1.572 + for (EdgeIt e(_graph); e != INVALID; ++e)
1.573 + c += (*_flow)[e] * _length[e];
1.574 + return c;
1.575 + }
1.576 +
1.577 + /// \brief Returns the number of the found paths.
1.578 + ///
1.579 + /// Returns the number of the found paths.
1.580 + ///
1.581 + /// \pre \ref run() or findFlow() must be called before using this
1.582 + /// function.
1.583 int pathNum() const {
1.584 - return paths.size();
1.585 + return _path_num;
1.586 }
1.587
1.588 + /// \brief Returns a const reference to the specified path.
1.589 + ///
1.590 + /// Returns a const reference to the specified path.
1.591 + ///
1.592 + /// \param i The function returns the \c i-th path.
1.593 + /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>.
1.594 + ///
1.595 + /// \pre \ref run() or findPaths() must be called before using this
1.596 + /// function.
1.597 + Path path(int i) const {
1.598 + return paths[i];
1.599 + }
1.600 +
1.601 + /// @}
1.602 +
1.603 }; //class Suurballe
1.604
1.605 ///@}
2.1 --- a/test/suurballe_test.cc Thu Feb 28 16:33:40 2008 +0000
2.2 +++ b/test/suurballe_test.cc Fri Feb 29 15:55:13 2008 +0000
2.3 @@ -17,94 +17,144 @@
2.4 */
2.5
2.6 #include <iostream>
2.7 +#include <fstream>
2.8 +
2.9 #include <lemon/list_graph.h>
2.10 +#include <lemon/graph_reader.h>
2.11 +#include <lemon/path.h>
2.12 #include <lemon/suurballe.h>
2.13 -//#include <path.h>
2.14 +
2.15 #include "test_tools.h"
2.16
2.17 using namespace lemon;
2.18
2.19 +// Checks the feasibility of the flow
2.20 +template <typename Graph, typename FlowMap>
2.21 +bool checkFlow( const Graph& gr, const FlowMap& flow,
2.22 + typename Graph::Node s, typename Graph::Node t,
2.23 + int value )
2.24 +{
2.25 + GRAPH_TYPEDEFS(typename Graph);
2.26 + for (EdgeIt e(gr); e != INVALID; ++e)
2.27 + if (!(flow[e] == 0 || flow[e] == 1)) return false;
2.28
2.29 -bool passed = true;
2.30 + for (NodeIt n(gr); n != INVALID; ++n) {
2.31 + int sum = 0;
2.32 + for (OutEdgeIt e(gr, n); e != INVALID; ++e)
2.33 + sum += flow[e];
2.34 + for (InEdgeIt e(gr, n); e != INVALID; ++e)
2.35 + sum -= flow[e];
2.36 + if (n == s && sum != value) return false;
2.37 + if (n == t && sum != -value) return false;
2.38 + if (n != s && n != t && sum != 0) return false;
2.39 + }
2.40 +
2.41 + return true;
2.42 +}
2.43 +
2.44 +// Checks the optimalitiy of the flow
2.45 +template < typename Graph, typename CostMap,
2.46 + typename FlowMap, typename PotentialMap >
2.47 +bool checkOptimality( const Graph& gr, const CostMap& cost,
2.48 + const FlowMap& flow, const PotentialMap& pi )
2.49 +{
2.50 + // Checking the Complementary Slackness optimality condition
2.51 + GRAPH_TYPEDEFS(typename Graph);
2.52 + bool opt = true;
2.53 + for (EdgeIt e(gr); e != INVALID; ++e) {
2.54 + typename CostMap::Value red_cost =
2.55 + cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
2.56 + opt = (flow[e] == 0 && red_cost >= 0) ||
2.57 + (flow[e] == 1 && red_cost <= 0);
2.58 + if (!opt) break;
2.59 + }
2.60 + return opt;
2.61 +}
2.62 +
2.63 +// Checks a path
2.64 +template < typename Graph, typename Path >
2.65 +bool checkPath( const Graph& gr, const Path& path,
2.66 + typename Graph::Node s, typename Graph::Node t)
2.67 +{
2.68 + // Checking the Complementary Slackness optimality condition
2.69 + GRAPH_TYPEDEFS(typename Graph);
2.70 + Node n = s;
2.71 + for (int i = 0; i < path.length(); ++i) {
2.72 + if (gr.source(path.nth(i)) != n) return false;
2.73 + n = gr.target(path.nth(i));
2.74 + }
2.75 + return n == t;
2.76 +}
2.77
2.78
2.79 int main()
2.80 {
2.81 - typedef ListGraph Graph;
2.82 - typedef Graph::Node Node;
2.83 - typedef Graph::Edge Edge;
2.84 + GRAPH_TYPEDEFS(ListGraph);
2.85
2.86 - Graph graph;
2.87 + // Reading the test graph
2.88 + ListGraph graph;
2.89 + ListGraph::EdgeMap<int> length(graph);
2.90 + Node source, target;
2.91
2.92 - //Ahuja könyv példája
2.93 + std::string fname;
2.94 + if(getenv("srcdir"))
2.95 + fname = std::string(getenv("srcdir"));
2.96 + else fname = ".";
2.97 + fname += "/test/min_cost_flow_test.lgf";
2.98
2.99 - Node s=graph.addNode();
2.100 - Node v1=graph.addNode();
2.101 - Node v2=graph.addNode();
2.102 - Node v3=graph.addNode();
2.103 - Node v4=graph.addNode();
2.104 - Node v5=graph.addNode();
2.105 - Node t=graph.addNode();
2.106 + std::ifstream input(fname.c_str());
2.107 + check(input, "Input file '" << fname << "' not found");
2.108 + GraphReader<ListGraph>(input, graph).
2.109 + readEdgeMap("cost", length).
2.110 + readNode("source", source).
2.111 + readNode("target", target).
2.112 + run();
2.113 + input.close();
2.114 +
2.115 + // Finding 2 paths
2.116 + {
2.117 + Suurballe<ListGraph> suurballe(graph, length, source, target);
2.118 + check(suurballe.run(2) == 2, "Wrong number of paths");
2.119 + check(checkFlow(graph, suurballe.flowMap(), source, target, 2),
2.120 + "The flow is not feasible");
2.121 + check(suurballe.totalLength() == 510, "The flow is not optimal");
2.122 + check(checkOptimality(graph, length, suurballe.flowMap(),
2.123 + suurballe.potentialMap()),
2.124 + "Wrong potentials");
2.125 + for (int i = 0; i < suurballe.pathNum(); ++i)
2.126 + check(checkPath(graph, suurballe.path(i), source, target),
2.127 + "Wrong path");
2.128 + }
2.129
2.130 - Edge s_v1=graph.addEdge(s, v1);
2.131 - Edge v1_v2=graph.addEdge(v1, v2);
2.132 - Edge s_v3=graph.addEdge(s, v3);
2.133 - Edge v2_v4=graph.addEdge(v2, v4);
2.134 - Edge v2_v5=graph.addEdge(v2, v5);
2.135 - Edge v3_v5=graph.addEdge(v3, v5);
2.136 - Edge v4_t=graph.addEdge(v4, t);
2.137 - Edge v5_t=graph.addEdge(v5, t);
2.138 -
2.139 + // Finding 3 paths
2.140 + {
2.141 + Suurballe<ListGraph> suurballe(graph, length, source, target);
2.142 + check(suurballe.run(3) == 3, "Wrong number of paths");
2.143 + check(checkFlow(graph, suurballe.flowMap(), source, target, 3),
2.144 + "The flow is not feasible");
2.145 + check(suurballe.totalLength() == 1040, "The flow is not optimal");
2.146 + check(checkOptimality(graph, length, suurballe.flowMap(),
2.147 + suurballe.potentialMap()),
2.148 + "Wrong potentials");
2.149 + for (int i = 0; i < suurballe.pathNum(); ++i)
2.150 + check(checkPath(graph, suurballe.path(i), source, target),
2.151 + "Wrong path");
2.152 + }
2.153
2.154 - Graph::EdgeMap<int> length(graph);
2.155 + // Finding 5 paths (only 3 can be found)
2.156 + {
2.157 + Suurballe<ListGraph> suurballe(graph, length, source, target);
2.158 + check(suurballe.run(5) == 3, "Wrong number of paths");
2.159 + check(checkFlow(graph, suurballe.flowMap(), source, target, 3),
2.160 + "The flow is not feasible");
2.161 + check(suurballe.totalLength() == 1040, "The flow is not optimal");
2.162 + check(checkOptimality(graph, length, suurballe.flowMap(),
2.163 + suurballe.potentialMap()),
2.164 + "Wrong potentials");
2.165 + for (int i = 0; i < suurballe.pathNum(); ++i)
2.166 + check(checkPath(graph, suurballe.path(i), source, target),
2.167 + "Wrong path");
2.168 + }
2.169
2.170 - length.set(s_v1, 6);
2.171 - length.set(v1_v2, 4);
2.172 - length.set(s_v3, 10);
2.173 - length.set(v2_v4, 5);
2.174 - length.set(v2_v5, 1);
2.175 - length.set(v3_v5, 5);
2.176 - length.set(v4_t, 8);
2.177 - length.set(v5_t, 8);
2.178 -
2.179 - std::cout << "Minlengthpaths algorithm test..." << std::endl;
2.180 -
2.181 -
2.182 - int k=3;
2.183 - Suurballe< Graph, Graph::EdgeMap<int> >
2.184 - surb_test(graph, length, s, t);
2.185 -
2.186 - check( surb_test.run(k) == 2 && surb_test.totalLength() == 46,
2.187 - "Two paths, total length should be 46");
2.188 -
2.189 - check( surb_test.checkComplementarySlackness(),
2.190 - "Complementary slackness conditions are not met.");
2.191 -
2.192 - // typedef DirPath<Graph> DPath;
2.193 - // DPath P(graph);
2.194 -
2.195 - /*
2.196 - surb_test.getPath(P,0);
2.197 - check(P.length() == 4, "First path should contain 4 edges.");
2.198 - std::cout<<P.length()<<std::endl;
2.199 - surb_test.getPath(P,1);
2.200 - check(P.length() == 3, "Second path: 3 edges.");
2.201 - std::cout<<P.length()<<std::endl;
2.202 - */
2.203 -
2.204 - k=1;
2.205 - check( surb_test.run(k) == 1 && surb_test.totalLength() == 19,
2.206 - "One path, total length should be 19");
2.207 -
2.208 - check( surb_test.checkComplementarySlackness(),
2.209 - "Complementary slackness conditions are not met.");
2.210 -
2.211 - // surb_test.getPath(P,0);
2.212 - // check(P.length() == 4, "First path should contain 4 edges.");
2.213 -
2.214 - std::cout << (passed ? "All tests passed." : "Some of the tests failed!!!")
2.215 - << std::endl;
2.216 -
2.217 - return passed ? 0 : 1;
2.218 -
2.219 + return 0;
2.220 }