[Lemon-commits] kpeter: r3469 - in lemon/trunk: lemon test

Lemon SVN svn at lemon.cs.elte.hu
Fri Feb 29 16:55:14 CET 2008


Author: kpeter
Date: Fri Feb 29 16:55:13 2008
New Revision: 3469

Modified:
   lemon/trunk/lemon/suurballe.h
   lemon/trunk/test/suurballe_test.cc

Log:
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.



Modified: lemon/trunk/lemon/suurballe.h
==============================================================================
--- lemon/trunk/lemon/suurballe.h	(original)
+++ lemon/trunk/lemon/suurballe.h	Fri Feb 29 16:55:13 2008
@@ -21,182 +21,474 @@
 
 ///\ingroup shortest_path
 ///\file
-///\brief An algorithm for finding k paths of minimal total length.
+///\brief An algorithm for finding edge-disjoint paths between two
+/// nodes having minimum total length.
 
-
-#include <lemon/maps.h>
 #include <vector>
+#include <lemon/bin_heap.h>
 #include <lemon/path.h>
-#include <lemon/ssp_min_cost_flow.h>
 
 namespace lemon {
 
-/// \addtogroup shortest_path
-/// @{
+  /// \addtogroup shortest_path
+  /// @{
 
-  ///\brief Implementation of an algorithm for finding k edge-disjoint
-  /// paths between 2 nodes of minimal total length
+  /// \brief Implementation of an algorithm for finding edge-disjoint
+  /// paths between two nodes having minimum total length.
   ///
-  /// The class \ref lemon::Suurballe implements
-  /// an algorithm for finding k edge-disjoint paths
-  /// from a given source node to a given target node in an
-  /// edge-weighted directed graph having minimal total weight (length).
-  ///
-  ///\warning Length values should be nonnegative!
-  /// 
-  ///\param Graph The directed graph type the algorithm runs on.
-  ///\param LengthMap The type of the length map (values should be nonnegative).
-  ///
-  ///\note It it questionable whether it is correct to call this method after
-  ///%Suurballe for it is just a special case of Edmonds' and Karp's algorithm
-  ///for finding minimum cost flows. In fact, this implementation just
-  ///wraps the SspMinCostFlow algorithms. The paper of both %Suurballe and
-  ///Edmonds-Karp published in 1972, therefore it is possibly right to
-  ///state that they are
-  ///independent results. Most frequently this special case is referred as
-  ///%Suurballe method in the literature, especially in communication
-  ///network context.
-  ///\author Attila Bernath
-  template <typename Graph, typename LengthMap>
-  class Suurballe{
-
+  /// \ref lemon::Suurballe "Suurballe" implements an algorithm for
+  /// finding edge-disjoint paths having minimum total length (cost)
+  /// from a given source node to a given target node in a directed
+  /// graph.
+  ///
+  /// In fact, this implementation is the specialization of the
+  /// \ref CapacityScaling "successive shortest path" algorithm.
+  ///
+  /// \tparam Graph The directed graph type the algorithm runs on.
+  /// \tparam LengthMap The type of the length (cost) map.
+  ///
+  /// \warning Length values should be \e non-negative \e integers.
+  ///
+  /// \note For finding node-disjoint paths this algorithm can be used
+  /// with \ref SplitGraphAdaptor.
+  ///
+  /// \author Attila Bernath and Peter Kovacs
+  
+  template < typename Graph, 
+             typename LengthMap = typename Graph::template EdgeMap<int> >
+  class Suurballe
+  {
+    GRAPH_TYPEDEFS(typename Graph);
 
     typedef typename LengthMap::Value Length;
-    
-    typedef typename Graph::Node Node;
-    typedef typename Graph::NodeIt NodeIt;
-    typedef typename Graph::Edge Edge;
-    typedef typename Graph::OutEdgeIt OutEdgeIt;
-    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
+    typedef ConstMap<Edge, int> ConstEdgeMap;
+    typedef typename Graph::template NodeMap<Edge> PredMap;
 
-    typedef ConstMap<Edge,int> ConstMap;
+  public:
 
-    const Graph& G;
+    /// The type of the flow map.
+    typedef typename Graph::template EdgeMap<int> FlowMap;
+    /// The type of the potential map.
+    typedef typename Graph::template NodeMap<Length> PotentialMap;
+    /// The type of the path structures.
+    typedef SimplePath<Graph> Path;
+
+  private:
+  
+    /// \brief Special implementation of the \ref Dijkstra algorithm
+    /// for finding shortest paths in the residual network.
+    ///
+    /// \ref ResidualDijkstra is a special implementation of the
+    /// \ref Dijkstra algorithm for finding shortest paths in the
+    /// residual network of the graph with respect to the reduced edge
+    /// lengths and modifying the node potentials according to the
+    /// distance of the nodes.
+    class ResidualDijkstra
+    {
+      typedef typename Graph::template NodeMap<int> HeapCrossRef;
+      typedef BinHeap<Length, HeapCrossRef> Heap;
+
+    private:
+
+      // The directed graph the algorithm runs on
+      const Graph &_graph;
+
+      // The main maps
+      const FlowMap &_flow;
+      const LengthMap &_length;
+      PotentialMap &_potential;
+
+      // The distance map
+      PotentialMap _dist;
+      // The pred edge map
+      PredMap &_pred;
+      // The processed (i.e. permanently labeled) nodes
+      std::vector<Node> _proc_nodes;
+      
+      Node _s;
+      Node _t;
 
-    Node s;
-    Node t;
+    public:
 
-    //Auxiliary variables
-    //This is the capacity map for the mincostflow problem
-    ConstMap const1map;
-    //This MinCostFlow instance will actually solve the problem
-    SspMinCostFlow<Graph, LengthMap, ConstMap> min_cost_flow;
+      /// Constructor.
+      ResidualDijkstra( const Graph &graph,
+                        const FlowMap &flow,
+                        const LengthMap &length,
+                        PotentialMap &potential,
+                        PredMap &pred,
+                        Node s, Node t ) :
+        _graph(graph), _flow(flow), _length(length), _potential(potential),
+        _dist(graph), _pred(pred), _s(s), _t(t) {}
+
+      /// \brief Runs the algorithm. Returns \c true if a path is found
+      /// from the source node to the target node.
+      bool run() {
+        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
+        Heap heap(heap_cross_ref);
+        heap.push(_s, 0);
+        _pred[_s] = INVALID;
+        _proc_nodes.clear();
+
+        // Processing nodes
+        while (!heap.empty() && heap.top() != _t) {
+          Node u = heap.top(), v;
+          Length d = heap.prio() + _potential[u], nd;
+          _dist[u] = heap.prio();
+          heap.pop();
+          _proc_nodes.push_back(u);
+
+          // Traversing outgoing edges
+          for (OutEdgeIt e(_graph, u); e != INVALID; ++e) {
+            if (_flow[e] == 0) {
+              v = _graph.target(e);
+              switch(heap.state(v)) {
+              case Heap::PRE_HEAP:
+                heap.push(v, d + _length[e] - _potential[v]);
+                _pred[v] = e;
+                break;
+              case Heap::IN_HEAP:
+                nd = d + _length[e] - _potential[v];
+                if (nd < heap[v]) {
+                  heap.decrease(v, nd);
+                  _pred[v] = e;
+                }
+                break;
+              case Heap::POST_HEAP:
+                break;
+              }
+            }
+          }
+
+          // Traversing incoming edges
+          for (InEdgeIt e(_graph, u); e != INVALID; ++e) {
+            if (_flow[e] == 1) {
+              v = _graph.source(e);
+              switch(heap.state(v)) {
+              case Heap::PRE_HEAP:
+                heap.push(v, d - _length[e] - _potential[v]);
+                _pred[v] = e;
+                break;
+              case Heap::IN_HEAP:
+                nd = d - _length[e] - _potential[v];
+                if (nd < heap[v]) {
+                  heap.decrease(v, nd);
+                  _pred[v] = e;
+                }
+                break;
+              case Heap::POST_HEAP:
+                break;
+              }
+            }
+          }
+        }
+        if (heap.empty()) return false;
+
+        // Updating potentials of processed nodes
+        Length t_dist = heap.prio();
+        for (int i = 0; i < int(_proc_nodes.size()); ++i)
+          _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
+        return true;
+      }
 
-    //Container to store found paths
-    std::vector<SimplePath<Graph> > paths;
+    }; //class ResidualDijkstra
 
-  public :
+  private:
 
+    // The directed graph the algorithm runs on
+    const Graph &_graph;
+    // The length map
+    const LengthMap &_length;
+    
+    // Edge map of the current flow
+    FlowMap *_flow;
+    bool _local_flow;
+    // Node map of the current potentials
+    PotentialMap *_potential;
+    bool _local_potential;
+
+    // The source node
+    Node _source;
+    // The target node
+    Node _target;
+
+    // Container to store the found paths
+    std::vector< SimplePath<Graph> > paths;
+    int _path_num;
+
+    // The pred edge map
+    PredMap _pred;
+    // Implementation of the Dijkstra algorithm for finding augmenting
+    // shortest paths in the residual network
+    ResidualDijkstra *_dijkstra;
+
+  public:
+
+    /// \brief Constructor.
+    ///
+    /// Constructor.
+    ///
+    /// \param graph The directed graph the algorithm runs on.
+    /// \param length The length (cost) values of the edges.
+    /// \param s The source node.
+    /// \param t The target node.
+    Suurballe( const Graph &graph,
+               const LengthMap &length,
+               Node s, Node t ) :
+      _graph(graph), _length(length), _flow(0), _local_flow(false),
+      _potential(0), _local_potential(false), _source(s), _target(t),
+      _pred(graph) {}
+
+    /// Destructor.
+    ~Suurballe() {
+      if (_local_flow) delete _flow;
+      if (_local_potential) delete _potential;
+      delete _dijkstra;
+    }
+
+    /// \brief Sets the flow map.
+    ///
+    /// Sets the flow map.
+    ///
+    /// The found flow contains only 0 and 1 values. It is the union of
+    /// the found edge-disjoint paths.
+    ///
+    /// \return \c (*this)
+    Suurballe& flowMap(FlowMap &map) {
+      if (_local_flow) {
+        delete _flow;
+        _local_flow = false;
+      }
+      _flow = ↦
+      return *this;
+    }
 
-    /// \brief The constructor of the class.
+    /// \brief Sets the potential map.
     ///
-    /// \param _G The directed graph the algorithm runs on. 
-    /// \param _length The length (weight or cost) of the edges. 
-    /// \param _s Source node.
-    /// \param _t Target node.
-    Suurballe(Graph& _G, LengthMap& _length, Node _s, Node _t) : 
-      G(_G), s(_s), t(_t), const1map(1), 
-      min_cost_flow(_G, _length, const1map, _s, _t) { }
+    /// Sets the potential map.
+    ///
+    /// The potentials provide the dual solution of the underlying 
+    /// minimum cost flow problem.
+    ///
+    /// \return \c (*this)
+    Suurballe& potentialMap(PotentialMap &map) {
+      if (_local_potential) {
+        delete _potential;
+        _local_potential = false;
+      }
+      _potential = ↦
+      return *this;
+    }
+
+    /// \name Execution control
+    /// The simplest way to execute the algorithm is to call the run()
+    /// function.
+    /// \n
+    /// If you only need the flow that is the union of the found
+    /// edge-disjoint paths, you may call init() and findFlow().
+
+    /// @{
 
     /// \brief Runs the algorithm.
     ///
     /// Runs the algorithm.
-    /// Returns k if there are at least k edge-disjoint paths from s to t.
-    /// Otherwise it returns the number of edge-disjoint paths found 
-    /// from s to t.
-    ///
-    /// \param k How many paths are we looking for?
-    ///
-    int run(int k) {
-      int i = min_cost_flow.run(k);
-
-      //Let's find the paths
-      //We put the paths into stl vectors (as an inner representation). 
-      //In the meantime we lose the information stored in 'reversed'.
-      //We suppose the lengths to be positive now.
-
-      //We don't want to change the flow of min_cost_flow, so we make a copy
-      //The name here suggests that the flow has only 0/1 values.
-      EdgeIntMap reversed(G); 
+    ///
+    /// \param k The number of paths to be found.
+    ///
+    /// \return \c k if there are at least \c k edge-disjoint paths
+    /// from \c s to \c t. Otherwise it returns the number of
+    /// edge-disjoint paths found.
+    ///
+    /// \note Apart from the return value, <tt>s.run(k)</tt> is just a
+    /// shortcut of the following code.
+    /// \code
+    ///   s.init();
+    ///   s.findFlow(k);
+    ///   s.findPaths();
+    /// \endcode
+    int run(int k = 2) {
+      init();
+      findFlow(k);
+      findPaths();
+      return _path_num;
+    }
 
-      for(typename Graph::EdgeIt e(G); e!=INVALID; ++e) 
-	reversed[e] = min_cost_flow.getFlow()[e];
-      
-      paths.clear();
-      paths.resize(k);
-      for (int j=0; j<i; ++j){
-	Node n=s;
-
-	while (n!=t){
-
-	  OutEdgeIt e(G, n);
-	  
-	  while (!reversed[e]){
-	    ++e;
-	  }
-	  n = G.target(e);
-	  paths[j].addBack(e);
-	  reversed[e] = 1-reversed[e];
-	}
-	
+    /// \brief Initializes the algorithm.
+    ///
+    /// Initializes the algorithm.
+    void init() {
+      // Initializing maps
+      if (!_flow) {
+        _flow = new FlowMap(_graph);
+        _local_flow = true;
       }
-      return i;
+      if (!_potential) {
+        _potential = new PotentialMap(_graph);
+        _local_potential = true;
+      }
+      for (EdgeIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
+
+      _dijkstra = new ResidualDijkstra( _graph, *_flow, _length, 
+                                        *_potential, _pred,
+                                        _source, _target );
     }
 
+    /// \brief Executes the successive shortest path algorithm to find
+    /// an optimal flow.
+    ///
+    /// Executes the successive shortest path algorithm to find a
+    /// minimum cost flow, which is the union of \c k or less
+    /// edge-disjoint paths.
+    ///
+    /// \return \c k if there are at least \c k edge-disjoint paths
+    /// from \c s to \c t. Otherwise it returns the number of
+    /// edge-disjoint paths found.
+    ///
+    /// \pre \ref init() must be called before using this function.
+    int findFlow(int k = 2) {
+      // Finding shortest paths
+      _path_num = 0;
+      while (_path_num < k) {
+        // Running Dijkstra
+        if (!_dijkstra->run()) break;
+        ++_path_num;
+
+        // Setting the flow along the found shortest path
+        Node u = _target;
+        Edge e;
+        while ((e = _pred[u]) != INVALID) {
+          if (u == _graph.target(e)) {
+            (*_flow)[e] = 1;
+            u = _graph.source(e);
+          } else {
+            (*_flow)[e] = 0;
+            u = _graph.target(e);
+          }
+        }
+      }
+      return _path_num;
+    }
     
-    /// \brief Returns the total length of the paths.
+    /// \brief Computes the paths from the flow.
     ///
-    /// This function gives back the total length of the found paths.
-    Length totalLength(){
-      return min_cost_flow.totalLength();
+    /// Computes the paths from the flow.
+    ///
+    /// \pre \ref init() and \ref findFlow() must be called before using
+    /// this function.
+    void findPaths() {
+      // Creating the residual flow map (the union of the paths not
+      // found so far)
+      FlowMap res_flow(*_flow);
+
+      paths.clear();
+      paths.resize(_path_num);
+      for (int i = 0; i < _path_num; ++i) {
+        Node n = _source;
+        while (n != _target) {
+          OutEdgeIt e(_graph, n);
+          for ( ; res_flow[e] == 0; ++e) ;
+          n = _graph.target(e);
+          paths[i].addBack(e);
+          res_flow[e] = 0;
+        }
+      }
     }
 
-    /// \brief Returns the found flow.
-    ///
-    /// This function returns a const reference to the EdgeMap \c flow.
-    const EdgeIntMap &getFlow() const { return min_cost_flow.flow;}
+    /// @}
+
+    /// \name Query Functions
+    /// The result of the algorithm can be obtained using these
+    /// functions.
+    /// \n The algorithm should be executed before using them.
+
+    /// @{
 
-    /// \brief Returns the optimal dual solution
+    /// \brief Returns a const reference to the edge map storing the
+    /// found flow.
     ///
-    /// This function returns a const reference to the NodeMap \c
-    /// potential (the dual solution).
-    const EdgeIntMap &getPotential() const { return min_cost_flow.potential;}
+    /// Returns a const reference to the edge map storing the flow that
+    /// is the union of the found edge-disjoint paths.
+    ///
+    /// \pre \ref run() or findFlow() must be called before using this
+    /// function.
+    const FlowMap& flowMap() const {
+      return *_flow;
+    }
 
-    /// \brief Checks whether the complementary slackness holds.
+    /// \brief Returns a const reference to the node map storing the
+    /// found potentials (the dual solution).
     ///
-    /// This function checks, whether the given solution is optimal.
-    /// Currently this function only checks optimality, doesn't bother
-    /// with feasibility.  It is meant for testing purposes.
-    bool checkComplementarySlackness(){
-      return min_cost_flow.checkComplementarySlackness();
+    /// Returns a const reference to the node map storing the found
+    /// potentials that provide the dual solution of the underlying 
+    /// minimum cost flow problem.
+    ///
+    /// \pre \ref run() or findFlow() must be called before using this
+    /// function.
+    const PotentialMap& potentialMap() const {
+      return *_potential;
     }
 
-    typedef SimplePath<Graph> Path; 
+    /// \brief Returns the flow on the given edge.
+    ///
+    /// Returns the flow on the given edge.
+    /// It is \c 1 if the edge is involved in one of the found paths,
+    /// otherwise it is \c 0.
+    ///
+    /// \pre \ref run() or findFlow() must be called before using this
+    /// function.
+    int flow(const Edge& edge) const {
+      return (*_flow)[edge];
+    }
 
-    /// \brief Read the found paths.
+    /// \brief Returns the potential of the given node.
     ///
-    /// This function gives back the \c j-th path in argument p.
-    /// Assumes that \c run() has been run and nothing has changed
-    /// since then.
+    /// Returns the potential of the given node.
+    ///
+    /// \pre \ref run() or findFlow() must be called before using this
+    /// function.
+    Length potential(const Node& node) const {
+      return (*_potential)[node];
+    }
+
+    /// \brief Returns the total length (cost) of the found paths (flow).
     ///
-    /// \warning It is assumed that \c p is constructed to be a path
-    /// of graph \c G.  If \c j is not less than the result of
-    /// previous \c run, then the result here will be an empty path
-    /// (\c j can be 0 as well).
+    /// Returns the total length (cost) of the found paths (flow).
+    /// The complexity of the function is \f$ O(e) \f$.
     ///
-    /// \param j Which path you want to get from the found paths (in a
-    /// real application you would get the found paths iteratively).
-    Path path(int j) const {
-      return paths[j];
+    /// \pre \ref run() or findFlow() must be called before using this
+    /// function.
+    Length totalLength() const {
+      Length c = 0;
+      for (EdgeIt e(_graph); e != INVALID; ++e)
+        c += (*_flow)[e] * _length[e];
+      return c;
     }
 
-    /// \brief Gives back the number of the paths.
+    /// \brief Returns the number of the found paths.
     ///
-    /// Gives back the number of the constructed paths.
+    /// Returns the number of the found paths.
+    ///
+    /// \pre \ref run() or findFlow() must be called before using this
+    /// function.
     int pathNum() const {
-      return paths.size();
+      return _path_num;
     }
 
+    /// \brief Returns a const reference to the specified path.
+    ///
+    /// Returns a const reference to the specified path.
+    ///
+    /// \param i The function returns the \c i-th path.
+    /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>.
+    ///
+    /// \pre \ref run() or findPaths() must be called before using this
+    /// function.
+    Path path(int i) const {
+      return paths[i];
+    }
+
+    /// @}
+
   }; //class Suurballe
 
   ///@}

Modified: lemon/trunk/test/suurballe_test.cc
==============================================================================
--- lemon/trunk/test/suurballe_test.cc	(original)
+++ lemon/trunk/test/suurballe_test.cc	Fri Feb 29 16:55:13 2008
@@ -17,94 +17,144 @@
  */
 
 #include <iostream>
+#include <fstream>
+
 #include <lemon/list_graph.h>
+#include <lemon/graph_reader.h>
+#include <lemon/path.h>
 #include <lemon/suurballe.h>
-//#include <path.h>
+
 #include "test_tools.h"
 
 using namespace lemon;
 
+// Checks the feasibility of the flow
+template <typename Graph, typename FlowMap>
+bool checkFlow( const Graph& gr, const FlowMap& flow, 
+                typename Graph::Node s, typename Graph::Node t,
+                int value )
+{
+  GRAPH_TYPEDEFS(typename Graph);
+  for (EdgeIt e(gr); e != INVALID; ++e)
+    if (!(flow[e] == 0 || flow[e] == 1)) return false;
+
+  for (NodeIt n(gr); n != INVALID; ++n) {
+    int sum = 0;
+    for (OutEdgeIt e(gr, n); e != INVALID; ++e)
+      sum += flow[e];
+    for (InEdgeIt e(gr, n); e != INVALID; ++e)
+      sum -= flow[e];
+    if (n == s && sum != value) return false;
+    if (n == t && sum != -value) return false;
+    if (n != s && n != t && sum != 0) return false;
+  }
 
-bool passed = true;
-
+  return true;
+}
 
-int main()
+// Checks the optimalitiy of the flow
+template < typename Graph, typename CostMap, 
+           typename FlowMap, typename PotentialMap >
+bool checkOptimality( const Graph& gr, const CostMap& cost,
+                      const FlowMap& flow, const PotentialMap& pi )
 {
-  typedef ListGraph Graph;
-  typedef Graph::Node Node;
-  typedef Graph::Edge Edge;
-
-  Graph graph;
-
-  //Ahuja könyv példája
-
-  Node s=graph.addNode();
-  Node v1=graph.addNode();  
-  Node v2=graph.addNode();
-  Node v3=graph.addNode();
-  Node v4=graph.addNode();
-  Node v5=graph.addNode();
-  Node t=graph.addNode();
-
-  Edge s_v1=graph.addEdge(s, v1);
-  Edge v1_v2=graph.addEdge(v1, v2);
-  Edge s_v3=graph.addEdge(s, v3);
-  Edge v2_v4=graph.addEdge(v2, v4);
-  Edge v2_v5=graph.addEdge(v2, v5);
-  Edge v3_v5=graph.addEdge(v3, v5);
-  Edge v4_t=graph.addEdge(v4, t);
-  Edge v5_t=graph.addEdge(v5, t);
-  
+  // Checking the Complementary Slackness optimality condition
+  GRAPH_TYPEDEFS(typename Graph);
+  bool opt = true;
+  for (EdgeIt e(gr); e != INVALID; ++e) {
+    typename CostMap::Value red_cost =
+      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
+    opt = (flow[e] == 0 && red_cost >= 0) ||
+          (flow[e] == 1 && red_cost <= 0);
+    if (!opt) break;
+  }
+  return opt;
+}
 
-  Graph::EdgeMap<int> length(graph);
+// Checks a path
+template < typename Graph, typename Path >
+bool checkPath( const Graph& gr, const Path& path,
+                typename Graph::Node s, typename Graph::Node t)
+{
+  // Checking the Complementary Slackness optimality condition
+  GRAPH_TYPEDEFS(typename Graph);
+  Node n = s;
+  for (int i = 0; i < path.length(); ++i) {
+    if (gr.source(path.nth(i)) != n) return false;
+    n = gr.target(path.nth(i));
+  }
+  return n == t;
+}
 
-  length.set(s_v1, 6);
-  length.set(v1_v2, 4);
-  length.set(s_v3, 10);
-  length.set(v2_v4, 5);
-  length.set(v2_v5, 1);
-  length.set(v3_v5, 5);
-  length.set(v4_t, 8);
-  length.set(v5_t, 8);
 
-  std::cout << "Minlengthpaths algorithm test..." << std::endl;
+int main()
+{
+  GRAPH_TYPEDEFS(ListGraph);
 
+  // Reading the test graph
+  ListGraph graph;
+  ListGraph::EdgeMap<int> length(graph);
+  Node source, target;
+
+  std::string fname;
+  if(getenv("srcdir"))
+    fname = std::string(getenv("srcdir"));
+  else fname = ".";
+  fname += "/test/min_cost_flow_test.lgf";
+
+  std::ifstream input(fname.c_str());
+  check(input, "Input file '" << fname << "' not found");
+  GraphReader<ListGraph>(input, graph).
+    readEdgeMap("cost", length).
+    readNode("source", source).
+    readNode("target", target).
+    run();
+  input.close();
   
-  int k=3;
-  Suurballe< Graph, Graph::EdgeMap<int> >
-    surb_test(graph, length, s, t);
-
-  check(  surb_test.run(k) == 2 && surb_test.totalLength() == 46,
-	  "Two paths, total length should be 46");
-
-  check(  surb_test.checkComplementarySlackness(),
-	  "Complementary slackness conditions are not met.");
-
-  //  typedef DirPath<Graph> DPath;
-  //  DPath P(graph);
-
-  /*
-  surb_test.getPath(P,0);
-  check(P.length() == 4, "First path should contain 4 edges.");  
-  std::cout<<P.length()<<std::endl;
-  surb_test.getPath(P,1);
-  check(P.length() == 3, "Second path: 3 edges.");
-  std::cout<<P.length()<<std::endl;
-  */  
-
-  k=1;
-  check(  surb_test.run(k) == 1 && surb_test.totalLength() == 19,
-	  "One path, total length should be 19");
-
-  check(  surb_test.checkComplementarySlackness(),
-	  "Complementary slackness conditions are not met.");
- 
-  //  surb_test.getPath(P,0);
-  //  check(P.length() == 4, "First path should contain 4 edges.");  
-
-  std::cout << (passed ? "All tests passed." : "Some of the tests failed!!!")
-	    << std::endl;
-
-  return passed ? 0 : 1;
+  // Finding 2 paths
+  {
+    Suurballe<ListGraph> suurballe(graph, length, source, target);
+    check(suurballe.run(2) == 2, "Wrong number of paths");
+    check(checkFlow(graph, suurballe.flowMap(), source, target, 2),
+          "The flow is not feasible");
+    check(suurballe.totalLength() == 510, "The flow is not optimal");
+    check(checkOptimality(graph, length, suurballe.flowMap(), 
+                          suurballe.potentialMap()),
+          "Wrong potentials");
+    for (int i = 0; i < suurballe.pathNum(); ++i)
+      check(checkPath(graph, suurballe.path(i), source, target),
+            "Wrong path");
+  }
+
+  // Finding 3 paths
+  {
+    Suurballe<ListGraph> suurballe(graph, length, source, target);
+    check(suurballe.run(3) == 3, "Wrong number of paths");
+    check(checkFlow(graph, suurballe.flowMap(), source, target, 3),
+          "The flow is not feasible");
+    check(suurballe.totalLength() == 1040, "The flow is not optimal");
+    check(checkOptimality(graph, length, suurballe.flowMap(), 
+                          suurballe.potentialMap()),
+          "Wrong potentials");
+    for (int i = 0; i < suurballe.pathNum(); ++i)
+      check(checkPath(graph, suurballe.path(i), source, target),
+            "Wrong path");
+  }
+
+  // Finding 5 paths (only 3 can be found)
+  {
+    Suurballe<ListGraph> suurballe(graph, length, source, target);
+    check(suurballe.run(5) == 3, "Wrong number of paths");
+    check(checkFlow(graph, suurballe.flowMap(), source, target, 3),
+          "The flow is not feasible");
+    check(suurballe.totalLength() == 1040, "The flow is not optimal");
+    check(checkOptimality(graph, length, suurballe.flowMap(), 
+                          suurballe.potentialMap()),
+          "Wrong potentials");
+    for (int i = 0; i < suurballe.pathNum(); ++i)
+      check(checkPath(graph, suurballe.path(i), source, target),
+            "Wrong path");
+  }
 
+  return 0;
 }



More information about the Lemon-commits mailing list