diff -r 20d42311e344 -r 37fb2c384c78 lemon/suurballe.h
--- a/lemon/suurballe.h	Thu Feb 28 16:33:40 2008 +0000
+++ b/lemon/suurballe.h	Fri Feb 29 15:55:13 2008 +0000
@@ -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).
+  /// \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.
   ///
-  ///\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).
+  /// In fact, this implementation is the specialization of the
+  /// \ref CapacityScaling "successive shortest path" algorithm.
   ///
-  ///\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{
-
+  /// \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 ConstMap<Edge, int> ConstEdgeMap;
+    typedef typename Graph::template NodeMap<Edge> PredMap;
+
+  public:
+
+    /// 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;
+
+    public:
+
+      /// 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;
+      }
+
+    }; //class ResidualDijkstra
+
+  private:
+
+    // The directed graph the algorithm runs on
+    const Graph &_graph;
+    // The length map
+    const LengthMap &_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;
+    // Edge map of the current flow
+    FlowMap *_flow;
+    bool _local_flow;
+    // Node map of the current potentials
+    PotentialMap *_potential;
+    bool _local_potential;
 
-    typedef ConstMap<Edge,int> ConstMap;
+    // The source node
+    Node _source;
+    // The target node
+    Node _target;
 
-    const Graph& G;
+    // Container to store the found paths
+    std::vector< SimplePath<Graph> > paths;
+    int _path_num;
 
-    Node s;
-    Node t;
+    // The pred edge map
+    PredMap _pred;
+    // Implementation of the Dijkstra algorithm for finding augmenting
+    // shortest paths in the residual network
+    ResidualDijkstra *_dijkstra;
 
-    //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;
+  public:
 
-    //Container to store found paths
-    std::vector<SimplePath<Graph> > paths;
+    /// \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) {}
 
-  public :
+    /// 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 = &map;
+      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 = &map;
+      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?
+    /// \param k The number of paths to be found.
     ///
-    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); 
-
-      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];
-	}
-	
-      }
-      return i;
+    /// \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;
     }
 
-    
-    /// \brief Returns the total length of the paths.
+    /// \brief Initializes the algorithm.
     ///
-    /// This function gives back the total length of the found paths.
-    Length totalLength(){
-      return min_cost_flow.totalLength();
+    /// Initializes the algorithm.
+    void init() {
+      // Initializing maps
+      if (!_flow) {
+        _flow = new FlowMap(_graph);
+        _local_flow = true;
+      }
+      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 Returns the found flow.
+    /// \brief Executes the successive shortest path algorithm to find
+    /// an optimal flow.
     ///
-    /// This function returns a const reference to the EdgeMap \c flow.
-    const EdgeIntMap &getFlow() const { return min_cost_flow.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;
 
-    /// \brief Returns the optimal dual solution
+        // 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 Computes the paths from the flow.
     ///
-    /// This function returns a const reference to the NodeMap \c
-    /// potential (the dual solution).
-    const EdgeIntMap &getPotential() const { return min_cost_flow.potential;}
+    /// 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);
 
-    /// \brief Checks whether the complementary slackness holds.
-    ///
-    /// 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();
+      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;
+        }
+      }
     }
 
-    typedef SimplePath<Graph> Path; 
+    /// @}
 
-    /// \brief Read the found paths.
+    /// \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 a const reference to the edge map storing the
+    /// found flow.
     ///
-    /// 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 a const reference to the edge map storing the flow that
+    /// is the union of the found edge-disjoint paths.
     ///
-    /// \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).
-    ///
-    /// \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.
+    const FlowMap& flowMap() const {
+      return *_flow;
     }
 
-    /// \brief Gives back the number of the paths.
+    /// \brief Returns a const reference to the node map storing the
+    /// found potentials (the dual solution).
     ///
-    /// Gives back the number of the constructed paths.
+    /// 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;
+    }
+
+    /// \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 Returns the potential of the given node.
+    ///
+    /// 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).
+    ///
+    /// Returns the total length (cost) of the found paths (flow).
+    /// The complexity of the function is \f$ O(e) \f$.
+    ///
+    /// \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 Returns the number of the found 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
 
   ///@}