lemon/suurballe.h
author Alpar Juttner <alpar@cs.elte.hu>
Fri, 25 Sep 2009 09:13:03 +0200
changeset 717 684964884a2e
parent 584 33c6b6e755cd
child 851 c67e235c832f
permissions -rw-r--r--
Merge
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_SUURBALLE_H
    20 #define LEMON_SUURBALLE_H
    21 
    22 ///\ingroup shortest_path
    23 ///\file
    24 ///\brief An algorithm for finding arc-disjoint paths between two
    25 /// nodes having minimum total length.
    26 
    27 #include <vector>
    28 #include <limits>
    29 #include <lemon/bin_heap.h>
    30 #include <lemon/path.h>
    31 #include <lemon/list_graph.h>
    32 #include <lemon/maps.h>
    33 
    34 namespace lemon {
    35 
    36   /// \addtogroup shortest_path
    37   /// @{
    38 
    39   /// \brief Algorithm for finding arc-disjoint paths between two nodes
    40   /// having minimum total length.
    41   ///
    42   /// \ref lemon::Suurballe "Suurballe" implements an algorithm for
    43   /// finding arc-disjoint paths having minimum total length (cost)
    44   /// from a given source node to a given target node in a digraph.
    45   ///
    46   /// Note that this problem is a special case of the \ref min_cost_flow
    47   /// "minimum cost flow problem". This implementation is actually an
    48   /// efficient specialized version of the \ref CapacityScaling
    49   /// "Successive Shortest Path" algorithm directly for this problem.
    50   /// Therefore this class provides query functions for flow values and
    51   /// node potentials (the dual solution) just like the minimum cost flow
    52   /// algorithms.
    53   ///
    54   /// \tparam GR The digraph type the algorithm runs on.
    55   /// \tparam LEN The type of the length map.
    56   /// The default value is <tt>GR::ArcMap<int></tt>.
    57   ///
    58   /// \warning Length values should be \e non-negative \e integers.
    59   ///
    60   /// \note For finding node-disjoint paths this algorithm can be used
    61   /// along with the \ref SplitNodes adaptor.
    62 #ifdef DOXYGEN
    63   template <typename GR, typename LEN>
    64 #else
    65   template < typename GR,
    66              typename LEN = typename GR::template ArcMap<int> >
    67 #endif
    68   class Suurballe
    69   {
    70     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    71 
    72     typedef ConstMap<Arc, int> ConstArcMap;
    73     typedef typename GR::template NodeMap<Arc> PredMap;
    74 
    75   public:
    76 
    77     /// The type of the digraph the algorithm runs on.
    78     typedef GR Digraph;
    79     /// The type of the length map.
    80     typedef LEN LengthMap;
    81     /// The type of the lengths.
    82     typedef typename LengthMap::Value Length;
    83 #ifdef DOXYGEN
    84     /// The type of the flow map.
    85     typedef GR::ArcMap<int> FlowMap;
    86     /// The type of the potential map.
    87     typedef GR::NodeMap<Length> PotentialMap;
    88 #else
    89     /// The type of the flow map.
    90     typedef typename Digraph::template ArcMap<int> FlowMap;
    91     /// The type of the potential map.
    92     typedef typename Digraph::template NodeMap<Length> PotentialMap;
    93 #endif
    94 
    95     /// The type of the path structures.
    96     typedef SimplePath<GR> Path;
    97 
    98   private:
    99 
   100     // ResidualDijkstra is a special implementation of the
   101     // Dijkstra algorithm for finding shortest paths in the
   102     // residual network with respect to the reduced arc lengths
   103     // and modifying the node potentials according to the
   104     // distance of the nodes.
   105     class ResidualDijkstra
   106     {
   107       typedef typename Digraph::template NodeMap<int> HeapCrossRef;
   108       typedef BinHeap<Length, HeapCrossRef> Heap;
   109 
   110     private:
   111 
   112       // The digraph the algorithm runs on
   113       const Digraph &_graph;
   114 
   115       // The main maps
   116       const FlowMap &_flow;
   117       const LengthMap &_length;
   118       PotentialMap &_potential;
   119 
   120       // The distance map
   121       PotentialMap _dist;
   122       // The pred arc map
   123       PredMap &_pred;
   124       // The processed (i.e. permanently labeled) nodes
   125       std::vector<Node> _proc_nodes;
   126 
   127       Node _s;
   128       Node _t;
   129 
   130     public:
   131 
   132       /// Constructor.
   133       ResidualDijkstra( const Digraph &graph,
   134                         const FlowMap &flow,
   135                         const LengthMap &length,
   136                         PotentialMap &potential,
   137                         PredMap &pred,
   138                         Node s, Node t ) :
   139         _graph(graph), _flow(flow), _length(length), _potential(potential),
   140         _dist(graph), _pred(pred), _s(s), _t(t) {}
   141 
   142       /// \brief Run the algorithm. It returns \c true if a path is found
   143       /// from the source node to the target node.
   144       bool run() {
   145         HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
   146         Heap heap(heap_cross_ref);
   147         heap.push(_s, 0);
   148         _pred[_s] = INVALID;
   149         _proc_nodes.clear();
   150 
   151         // Process nodes
   152         while (!heap.empty() && heap.top() != _t) {
   153           Node u = heap.top(), v;
   154           Length d = heap.prio() + _potential[u], nd;
   155           _dist[u] = heap.prio();
   156           heap.pop();
   157           _proc_nodes.push_back(u);
   158 
   159           // Traverse outgoing arcs
   160           for (OutArcIt e(_graph, u); e != INVALID; ++e) {
   161             if (_flow[e] == 0) {
   162               v = _graph.target(e);
   163               switch(heap.state(v)) {
   164               case Heap::PRE_HEAP:
   165                 heap.push(v, d + _length[e] - _potential[v]);
   166                 _pred[v] = e;
   167                 break;
   168               case Heap::IN_HEAP:
   169                 nd = d + _length[e] - _potential[v];
   170                 if (nd < heap[v]) {
   171                   heap.decrease(v, nd);
   172                   _pred[v] = e;
   173                 }
   174                 break;
   175               case Heap::POST_HEAP:
   176                 break;
   177               }
   178             }
   179           }
   180 
   181           // Traverse incoming arcs
   182           for (InArcIt e(_graph, u); e != INVALID; ++e) {
   183             if (_flow[e] == 1) {
   184               v = _graph.source(e);
   185               switch(heap.state(v)) {
   186               case Heap::PRE_HEAP:
   187                 heap.push(v, d - _length[e] - _potential[v]);
   188                 _pred[v] = e;
   189                 break;
   190               case Heap::IN_HEAP:
   191                 nd = d - _length[e] - _potential[v];
   192                 if (nd < heap[v]) {
   193                   heap.decrease(v, nd);
   194                   _pred[v] = e;
   195                 }
   196                 break;
   197               case Heap::POST_HEAP:
   198                 break;
   199               }
   200             }
   201           }
   202         }
   203         if (heap.empty()) return false;
   204 
   205         // Update potentials of processed nodes
   206         Length t_dist = heap.prio();
   207         for (int i = 0; i < int(_proc_nodes.size()); ++i)
   208           _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
   209         return true;
   210       }
   211 
   212     }; //class ResidualDijkstra
   213 
   214   private:
   215 
   216     // The digraph the algorithm runs on
   217     const Digraph &_graph;
   218     // The length map
   219     const LengthMap &_length;
   220 
   221     // Arc map of the current flow
   222     FlowMap *_flow;
   223     bool _local_flow;
   224     // Node map of the current potentials
   225     PotentialMap *_potential;
   226     bool _local_potential;
   227 
   228     // The source node
   229     Node _source;
   230     // The target node
   231     Node _target;
   232 
   233     // Container to store the found paths
   234     std::vector< SimplePath<Digraph> > paths;
   235     int _path_num;
   236 
   237     // The pred arc map
   238     PredMap _pred;
   239     // Implementation of the Dijkstra algorithm for finding augmenting
   240     // shortest paths in the residual network
   241     ResidualDijkstra *_dijkstra;
   242 
   243   public:
   244 
   245     /// \brief Constructor.
   246     ///
   247     /// Constructor.
   248     ///
   249     /// \param graph The digraph the algorithm runs on.
   250     /// \param length The length (cost) values of the arcs.
   251     Suurballe( const Digraph &graph,
   252                const LengthMap &length ) :
   253       _graph(graph), _length(length), _flow(0), _local_flow(false),
   254       _potential(0), _local_potential(false), _pred(graph)
   255     {
   256       LEMON_ASSERT(std::numeric_limits<Length>::is_integer,
   257         "The length type of Suurballe must be integer");
   258     }
   259 
   260     /// Destructor.
   261     ~Suurballe() {
   262       if (_local_flow) delete _flow;
   263       if (_local_potential) delete _potential;
   264       delete _dijkstra;
   265     }
   266 
   267     /// \brief Set the flow map.
   268     ///
   269     /// This function sets the flow map.
   270     /// If it is not used before calling \ref run() or \ref init(),
   271     /// an instance will be allocated automatically. The destructor
   272     /// deallocates this automatically allocated map, of course.
   273     ///
   274     /// The found flow contains only 0 and 1 values, since it is the
   275     /// union of the found arc-disjoint paths.
   276     ///
   277     /// \return <tt>(*this)</tt>
   278     Suurballe& flowMap(FlowMap &map) {
   279       if (_local_flow) {
   280         delete _flow;
   281         _local_flow = false;
   282       }
   283       _flow = &map;
   284       return *this;
   285     }
   286 
   287     /// \brief Set the potential map.
   288     ///
   289     /// This function sets the potential map.
   290     /// If it is not used before calling \ref run() or \ref init(),
   291     /// an instance will be allocated automatically. The destructor
   292     /// deallocates this automatically allocated map, of course.
   293     ///
   294     /// The node potentials provide the dual solution of the underlying
   295     /// \ref min_cost_flow "minimum cost flow problem".
   296     ///
   297     /// \return <tt>(*this)</tt>
   298     Suurballe& potentialMap(PotentialMap &map) {
   299       if (_local_potential) {
   300         delete _potential;
   301         _local_potential = false;
   302       }
   303       _potential = &map;
   304       return *this;
   305     }
   306 
   307     /// \name Execution Control
   308     /// The simplest way to execute the algorithm is to call the run()
   309     /// function.
   310     /// \n
   311     /// If you only need the flow that is the union of the found
   312     /// arc-disjoint paths, you may call init() and findFlow().
   313 
   314     /// @{
   315 
   316     /// \brief Run the algorithm.
   317     ///
   318     /// This function runs the algorithm.
   319     ///
   320     /// \param s The source node.
   321     /// \param t The target node.
   322     /// \param k The number of paths to be found.
   323     ///
   324     /// \return \c k if there are at least \c k arc-disjoint paths from
   325     /// \c s to \c t in the digraph. Otherwise it returns the number of
   326     /// arc-disjoint paths found.
   327     ///
   328     /// \note Apart from the return value, <tt>s.run(s, t, k)</tt> is
   329     /// just a shortcut of the following code.
   330     /// \code
   331     ///   s.init(s);
   332     ///   s.findFlow(t, k);
   333     ///   s.findPaths();
   334     /// \endcode
   335     int run(const Node& s, const Node& t, int k = 2) {
   336       init(s);
   337       findFlow(t, k);
   338       findPaths();
   339       return _path_num;
   340     }
   341 
   342     /// \brief Initialize the algorithm.
   343     ///
   344     /// This function initializes the algorithm.
   345     ///
   346     /// \param s The source node.
   347     void init(const Node& s) {
   348       _source = s;
   349 
   350       // Initialize maps
   351       if (!_flow) {
   352         _flow = new FlowMap(_graph);
   353         _local_flow = true;
   354       }
   355       if (!_potential) {
   356         _potential = new PotentialMap(_graph);
   357         _local_potential = true;
   358       }
   359       for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
   360       for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
   361     }
   362 
   363     /// \brief Execute the algorithm to find an optimal flow.
   364     ///
   365     /// This function executes the successive shortest path algorithm to
   366     /// find a minimum cost flow, which is the union of \c k (or less)
   367     /// arc-disjoint paths.
   368     ///
   369     /// \param t The target node.
   370     /// \param k The number of paths to be found.
   371     ///
   372     /// \return \c k if there are at least \c k arc-disjoint paths from
   373     /// the source node to the given node \c t in the digraph.
   374     /// Otherwise it returns the number of arc-disjoint paths found.
   375     ///
   376     /// \pre \ref init() must be called before using this function.
   377     int findFlow(const Node& t, int k = 2) {
   378       _target = t;
   379       _dijkstra =
   380         new ResidualDijkstra( _graph, *_flow, _length, *_potential, _pred,
   381                               _source, _target );
   382 
   383       // Find shortest paths
   384       _path_num = 0;
   385       while (_path_num < k) {
   386         // Run Dijkstra
   387         if (!_dijkstra->run()) break;
   388         ++_path_num;
   389 
   390         // Set the flow along the found shortest path
   391         Node u = _target;
   392         Arc e;
   393         while ((e = _pred[u]) != INVALID) {
   394           if (u == _graph.target(e)) {
   395             (*_flow)[e] = 1;
   396             u = _graph.source(e);
   397           } else {
   398             (*_flow)[e] = 0;
   399             u = _graph.target(e);
   400           }
   401         }
   402       }
   403       return _path_num;
   404     }
   405 
   406     /// \brief Compute the paths from the flow.
   407     ///
   408     /// This function computes the paths from the found minimum cost flow,
   409     /// which is the union of some arc-disjoint paths.
   410     ///
   411     /// \pre \ref init() and \ref findFlow() must be called before using
   412     /// this function.
   413     void findPaths() {
   414       FlowMap res_flow(_graph);
   415       for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a];
   416 
   417       paths.clear();
   418       paths.resize(_path_num);
   419       for (int i = 0; i < _path_num; ++i) {
   420         Node n = _source;
   421         while (n != _target) {
   422           OutArcIt e(_graph, n);
   423           for ( ; res_flow[e] == 0; ++e) ;
   424           n = _graph.target(e);
   425           paths[i].addBack(e);
   426           res_flow[e] = 0;
   427         }
   428       }
   429     }
   430 
   431     /// @}
   432 
   433     /// \name Query Functions
   434     /// The results of the algorithm can be obtained using these
   435     /// functions.
   436     /// \n The algorithm should be executed before using them.
   437 
   438     /// @{
   439 
   440     /// \brief Return the total length of the found paths.
   441     ///
   442     /// This function returns the total length of the found paths, i.e.
   443     /// the total cost of the found flow.
   444     /// The complexity of the function is O(e).
   445     ///
   446     /// \pre \ref run() or \ref findFlow() must be called before using
   447     /// this function.
   448     Length totalLength() const {
   449       Length c = 0;
   450       for (ArcIt e(_graph); e != INVALID; ++e)
   451         c += (*_flow)[e] * _length[e];
   452       return c;
   453     }
   454 
   455     /// \brief Return the flow value on the given arc.
   456     ///
   457     /// This function returns the flow value on the given arc.
   458     /// It is \c 1 if the arc is involved in one of the found arc-disjoint
   459     /// paths, otherwise it is \c 0.
   460     ///
   461     /// \pre \ref run() or \ref findFlow() must be called before using
   462     /// this function.
   463     int flow(const Arc& arc) const {
   464       return (*_flow)[arc];
   465     }
   466 
   467     /// \brief Return a const reference to an arc map storing the
   468     /// found flow.
   469     ///
   470     /// This function returns a const reference to an arc map storing
   471     /// the flow that is the union of the found arc-disjoint paths.
   472     ///
   473     /// \pre \ref run() or \ref findFlow() must be called before using
   474     /// this function.
   475     const FlowMap& flowMap() const {
   476       return *_flow;
   477     }
   478 
   479     /// \brief Return the potential of the given node.
   480     ///
   481     /// This function returns the potential of the given node.
   482     /// The node potentials provide the dual solution of the
   483     /// underlying \ref min_cost_flow "minimum cost flow problem".
   484     ///
   485     /// \pre \ref run() or \ref findFlow() must be called before using
   486     /// this function.
   487     Length potential(const Node& node) const {
   488       return (*_potential)[node];
   489     }
   490 
   491     /// \brief Return a const reference to a node map storing the
   492     /// found potentials (the dual solution).
   493     ///
   494     /// This function returns a const reference to a node map storing
   495     /// the found potentials that provide the dual solution of the
   496     /// underlying \ref min_cost_flow "minimum cost flow problem".
   497     ///
   498     /// \pre \ref run() or \ref findFlow() must be called before using
   499     /// this function.
   500     const PotentialMap& potentialMap() const {
   501       return *_potential;
   502     }
   503 
   504     /// \brief Return the number of the found paths.
   505     ///
   506     /// This function returns the number of the found paths.
   507     ///
   508     /// \pre \ref run() or \ref findFlow() must be called before using
   509     /// this function.
   510     int pathNum() const {
   511       return _path_num;
   512     }
   513 
   514     /// \brief Return a const reference to the specified path.
   515     ///
   516     /// This function returns a const reference to the specified path.
   517     ///
   518     /// \param i The function returns the <tt>i</tt>-th path.
   519     /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>.
   520     ///
   521     /// \pre \ref run() or \ref findPaths() must be called before using
   522     /// this function.
   523     Path path(int i) const {
   524       return paths[i];
   525     }
   526 
   527     /// @}
   528 
   529   }; //class Suurballe
   530 
   531   ///@}
   532 
   533 } //namespace lemon
   534 
   535 #endif //LEMON_SUURBALLE_H