RhodeCode

namespace lemon {

  /// \addtogroup shortest_path

  /// @{

  ///

  /// \ref lemon::Suurballe "Suurballe" implements an algorithm for

  /// finding arc-disjoint paths having minimum total length (cost)

  ///

  /// In fact, this implementation is the specialization of the

  /// \ref CapacityScaling "successive shortest path" algorithm.

  ///

  /// \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 SplitDigraphAdaptor.

             typename LengthMap = typename Digraph::template ArcMap<int> >

  class Suurballe

  {

    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);

    typedef typename LengthMap::Value Length;

    typedef ConstMap<Arc, int> ConstArcMap;

    typedef typename Digraph::template NodeMap<Length> PotentialMap;

    /// The type of the path structures.

    typedef SimplePath<Digraph> Path;

  private:

    /// 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 digraph with respect to the reduced arc

    /// lengths and modifying the node potentials according to the

    {

      typedef typename Digraph::template NodeMap<int> HeapCrossRef;

      typedef BinHeap<Length, HeapCrossRef> Heap;

    private:

      const Digraph &_graph;

      // The main maps

      const FlowMap &_flow;

      const LengthMap &_length;

      PotentialMap &_potential;

                        PotentialMap &potential,

                        PredMap &pred,

                        Node s, Node t ) :

        _graph(digraph), _flow(flow), _length(length), _potential(potential),

        _dist(digraph), _pred(pred), _s(s), _t(t) {}

      /// 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();

        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);

          for (OutArcIt 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]);

              case Heap::POST_HEAP:

                break;

              }

            }

          }

          for (InArcIt 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]);

              }

            }

          }

        }

        if (heap.empty()) return false;

        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:

    const Digraph &_graph;

    // The length map

    const LengthMap &_length;

    // Arc map of the current flow

    FlowMap *_flow;

  public:

    /// \brief Constructor.

    ///

    /// Constructor.

    ///

    /// \param length The length (cost) values of the arcs.

    /// \param s The source node.

    /// \param t The target node.

    Suurballe( const Digraph &digraph,

               const LengthMap &length,

               Node s, Node t ) :

    ~Suurballe() {

      if (_local_flow) delete _flow;

      if (_local_potential) delete _potential;

      delete _dijkstra;

    }

    ///

    ///

    /// The found flow contains only 0 and 1 values. It is the union of

    /// the found arc-disjoint paths.

    ///

    /// \return \c (*this)

    Suurballe& flowMap(FlowMap &map) {

        _local_flow = false;

      }

      _flow = ↦

      return *this;

    }

    ///

    ///

    /// The potentials provide the dual solution of the underlying 

    /// minimum cost flow problem.

    ///

    /// \return \c (*this)

    Suurballe& potentialMap(PotentialMap &map) {

    /// \n

    /// If you only need the flow that is the union of the found

    /// arc-disjoint paths, you may call init() and findFlow().

    /// @{

    ///

    ///

    /// \param k The number of paths to be found.

    ///

    /// arc-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();

      init();

      findFlow(k);

      findPaths();

      return _path_num;

    }

    ///

    void init() {

      if (!_flow) {

        _flow = new FlowMap(_graph);

        _local_flow = true;

      }

      if (!_potential) {

        _potential = new PotentialMap(_graph);

      _dijkstra = new ResidualDijkstra( _graph, *_flow, _length, 

                                        *_potential, _pred,

                                        _source, _target );

    }

    /// an optimal flow.

    ///

    /// arc-disjoint paths.

    ///

    /// arc-disjoint paths found.

    ///

    /// \pre \ref init() must be called before using this function.

    int findFlow(int k = 2) {

      _path_num = 0;

      while (_path_num < k) {

        if (!_dijkstra->run()) break;

        ++_path_num;

        Node u = _target;

        Arc e;

        while ((e = _pred[u]) != INVALID) {

          if (u == _graph.target(e)) {

            (*_flow)[e] = 1;

            u = _graph.source(e);

          }

        }

      }

      return _path_num;

    }

    ///

    ///

    /// \pre \ref init() and \ref findFlow() must be called before using

    /// this function.

    void findPaths() {

      FlowMap res_flow(_graph);

      paths.clear();

      paths.resize(_path_num);

      for (int i = 0; i < _path_num; ++i) {

        Node n = _source;

        while (n != _target) {

      }

    }

    /// @}

    /// \name Query Functions

    /// functions.

    /// \n The algorithm should be executed before using them.

    /// @{

    /// found flow.

    ///

    ///

    const FlowMap& flowMap() const {

      return *_flow;

    }

    /// found potentials (the dual solution).

    ///

    ///

    const PotentialMap& potentialMap() const {

      return *_potential;

    }

    ///

    /// It is \c 1 if the arc is involved in one of the found paths,

    /// otherwise it is \c 0.

    ///

    int flow(const Arc& arc) const {

      return (*_flow)[arc];

    }

    ///

    ///

    Length potential(const Node& node) const {

      return (*_potential)[node];

    }

    ///

    ///

    Length totalLength() const {

      Length c = 0;

      for (ArcIt e(_graph); e != INVALID; ++e)

        c += (*_flow)[e] * _length[e];

      return c;

    }

    ///

    ///

    int pathNum() const {

      return _path_num;

    }

    ///

    ///

    /// \param i The function returns the \c i-th path.

    /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>.

    ///

    Path path(int i) const {

      return paths[i];

    }

    /// @}

#include <lemon/suurballe.h>

#include "test_tools.h"

using namespace lemon;

template <typename Digraph, typename FlowMap>

bool checkFlow( const Digraph& gr, const FlowMap& flow, 

                typename Digraph::Node s, typename Digraph::Node t,

                int value )

{

  TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);

    if (n != s && n != t && sum != 0) return false;

  }

  return true;

}

template < typename Digraph, typename CostMap, 

           typename FlowMap, typename PotentialMap >

bool checkOptimality( const Digraph& gr, const CostMap& cost,

                      const FlowMap& flow, const PotentialMap& pi )

{

  TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);

  bool opt = true;

  for (ArcIt 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;

}

bool checkPath( const Digraph& gr, const Path& path,

                typename Digraph::Node s, typename Digraph::Node t)

{

  TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);

  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));

  }

int main()

{

  DIGRAPH_TYPEDEFS(ListDigraph);

  ListDigraph digraph;

  ListDigraph::ArcMap<int> length(digraph);

  Node source, target;

  std::string fname;

  if(getenv("srcdir"))

    arcMap("cost", length).

    node("source", source).

    node("target", target).

    run();

  input.close();

  {

    Suurballe<ListDigraph> suurballe(digraph, length, source, target);

    check(suurballe.run(2) == 2, "Wrong number of paths");

    check(checkFlow(digraph, suurballe.flowMap(), source, target, 2),

          "The flow is not feasible");

    check(suurballe.totalLength() == 510, "The flow is not optimal");

          "Wrong potentials");

    for (int i = 0; i < suurballe.pathNum(); ++i)

      check(checkPath(digraph, suurballe.path(i), source, target),

            "Wrong path");

  }

  {

    Suurballe<ListDigraph> suurballe(digraph, length, source, target);

    check(suurballe.run(3) == 3, "Wrong number of paths");

    check(checkFlow(digraph, suurballe.flowMap(), source, target, 3),

          "The flow is not feasible");

    check(suurballe.totalLength() == 1040, "The flow is not optimal");

          "Wrong potentials");

    for (int i = 0; i < suurballe.pathNum(); ++i)

      check(checkPath(digraph, suurballe.path(i), source, target),

            "Wrong path");

  }

  {

    Suurballe<ListDigraph> suurballe(digraph, length, source, target);

    check(suurballe.run(5) == 3, "Wrong number of paths");

    check(checkFlow(digraph, suurballe.flowMap(), source, target, 3),

          "The flow is not feasible");

    check(suurballe.totalLength() == 1040, "The flow is not optimal");