Port Suurballe algorithm from svn -r3512
authorAlpar Juttner <alpar@cs.elte.hu>
Tue, 28 Oct 2008 18:39:53 +0000
changeset 3572f64c4a692a8
parent 355 956a29f30887
child 358 7f26c4b32651
Port Suurballe algorithm from svn -r3512
lemon/Makefile.am
lemon/suurballe.h
test/Makefile.am
test/min_cost_flow_test.lgf
test/suurballe_test.cc
     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 = &map;
   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 = &map;
   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 +}