lemon/suurballe.h
author Alpar Juttner <alpar@cs.elte.hu>
Fri, 15 Apr 2011 09:37:47 +0200
changeset 945 2024bf2ecdd6
parent 851 c67e235c832f
child 853 ec0b1b423b8b
permissions -rw-r--r--
Build of mip_test and lp_test precede the running of the tests (#419)
     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.
    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 
   257     /// Destructor.
   258     ~Suurballe() {
   259       if (_local_flow) delete _flow;
   260       if (_local_potential) delete _potential;
   261       delete _dijkstra;
   262     }
   263 
   264     /// \brief Set the flow map.
   265     ///
   266     /// This function sets the flow map.
   267     /// If it is not used before calling \ref run() or \ref init(),
   268     /// an instance will be allocated automatically. The destructor
   269     /// deallocates this automatically allocated map, of course.
   270     ///
   271     /// The found flow contains only 0 and 1 values, since it is the
   272     /// union of the found arc-disjoint paths.
   273     ///
   274     /// \return <tt>(*this)</tt>
   275     Suurballe& flowMap(FlowMap &map) {
   276       if (_local_flow) {
   277         delete _flow;
   278         _local_flow = false;
   279       }
   280       _flow = &map;
   281       return *this;
   282     }
   283 
   284     /// \brief Set the potential map.
   285     ///
   286     /// This function sets the potential map.
   287     /// If it is not used before calling \ref run() or \ref init(),
   288     /// an instance will be allocated automatically. The destructor
   289     /// deallocates this automatically allocated map, of course.
   290     ///
   291     /// The node potentials provide the dual solution of the underlying
   292     /// \ref min_cost_flow "minimum cost flow problem".
   293     ///
   294     /// \return <tt>(*this)</tt>
   295     Suurballe& potentialMap(PotentialMap &map) {
   296       if (_local_potential) {
   297         delete _potential;
   298         _local_potential = false;
   299       }
   300       _potential = &map;
   301       return *this;
   302     }
   303 
   304     /// \name Execution Control
   305     /// The simplest way to execute the algorithm is to call the run()
   306     /// function.
   307     /// \n
   308     /// If you only need the flow that is the union of the found
   309     /// arc-disjoint paths, you may call init() and findFlow().
   310 
   311     /// @{
   312 
   313     /// \brief Run the algorithm.
   314     ///
   315     /// This function runs the algorithm.
   316     ///
   317     /// \param s The source node.
   318     /// \param t The target node.
   319     /// \param k The number of paths to be found.
   320     ///
   321     /// \return \c k if there are at least \c k arc-disjoint paths from
   322     /// \c s to \c t in the digraph. Otherwise it returns the number of
   323     /// arc-disjoint paths found.
   324     ///
   325     /// \note Apart from the return value, <tt>s.run(s, t, k)</tt> is
   326     /// just a shortcut of the following code.
   327     /// \code
   328     ///   s.init(s);
   329     ///   s.findFlow(t, k);
   330     ///   s.findPaths();
   331     /// \endcode
   332     int run(const Node& s, const Node& t, int k = 2) {
   333       init(s);
   334       findFlow(t, k);
   335       findPaths();
   336       return _path_num;
   337     }
   338 
   339     /// \brief Initialize the algorithm.
   340     ///
   341     /// This function initializes the algorithm.
   342     ///
   343     /// \param s The source node.
   344     void init(const Node& s) {
   345       _source = s;
   346 
   347       // Initialize maps
   348       if (!_flow) {
   349         _flow = new FlowMap(_graph);
   350         _local_flow = true;
   351       }
   352       if (!_potential) {
   353         _potential = new PotentialMap(_graph);
   354         _local_potential = true;
   355       }
   356       for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
   357       for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
   358     }
   359 
   360     /// \brief Execute the algorithm to find an optimal flow.
   361     ///
   362     /// This function executes the successive shortest path algorithm to
   363     /// find a minimum cost flow, which is the union of \c k (or less)
   364     /// arc-disjoint paths.
   365     ///
   366     /// \param t The target node.
   367     /// \param k The number of paths to be found.
   368     ///
   369     /// \return \c k if there are at least \c k arc-disjoint paths from
   370     /// the source node to the given node \c t in the digraph.
   371     /// Otherwise it returns the number of arc-disjoint paths found.
   372     ///
   373     /// \pre \ref init() must be called before using this function.
   374     int findFlow(const Node& t, int k = 2) {
   375       _target = t;
   376       _dijkstra =
   377         new ResidualDijkstra( _graph, *_flow, _length, *_potential, _pred,
   378                               _source, _target );
   379 
   380       // Find shortest paths
   381       _path_num = 0;
   382       while (_path_num < k) {
   383         // Run Dijkstra
   384         if (!_dijkstra->run()) break;
   385         ++_path_num;
   386 
   387         // Set the flow along the found shortest path
   388         Node u = _target;
   389         Arc e;
   390         while ((e = _pred[u]) != INVALID) {
   391           if (u == _graph.target(e)) {
   392             (*_flow)[e] = 1;
   393             u = _graph.source(e);
   394           } else {
   395             (*_flow)[e] = 0;
   396             u = _graph.target(e);
   397           }
   398         }
   399       }
   400       return _path_num;
   401     }
   402 
   403     /// \brief Compute the paths from the flow.
   404     ///
   405     /// This function computes the paths from the found minimum cost flow,
   406     /// which is the union of some arc-disjoint paths.
   407     ///
   408     /// \pre \ref init() and \ref findFlow() must be called before using
   409     /// this function.
   410     void findPaths() {
   411       FlowMap res_flow(_graph);
   412       for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a];
   413 
   414       paths.clear();
   415       paths.resize(_path_num);
   416       for (int i = 0; i < _path_num; ++i) {
   417         Node n = _source;
   418         while (n != _target) {
   419           OutArcIt e(_graph, n);
   420           for ( ; res_flow[e] == 0; ++e) ;
   421           n = _graph.target(e);
   422           paths[i].addBack(e);
   423           res_flow[e] = 0;
   424         }
   425       }
   426     }
   427 
   428     /// @}
   429 
   430     /// \name Query Functions
   431     /// The results of the algorithm can be obtained using these
   432     /// functions.
   433     /// \n The algorithm should be executed before using them.
   434 
   435     /// @{
   436 
   437     /// \brief Return the total length of the found paths.
   438     ///
   439     /// This function returns the total length of the found paths, i.e.
   440     /// the total cost of the found flow.
   441     /// The complexity of the function is O(e).
   442     ///
   443     /// \pre \ref run() or \ref findFlow() must be called before using
   444     /// this function.
   445     Length totalLength() const {
   446       Length c = 0;
   447       for (ArcIt e(_graph); e != INVALID; ++e)
   448         c += (*_flow)[e] * _length[e];
   449       return c;
   450     }
   451 
   452     /// \brief Return the flow value on the given arc.
   453     ///
   454     /// This function returns the flow value on the given arc.
   455     /// It is \c 1 if the arc is involved in one of the found arc-disjoint
   456     /// paths, otherwise it is \c 0.
   457     ///
   458     /// \pre \ref run() or \ref findFlow() must be called before using
   459     /// this function.
   460     int flow(const Arc& arc) const {
   461       return (*_flow)[arc];
   462     }
   463 
   464     /// \brief Return a const reference to an arc map storing the
   465     /// found flow.
   466     ///
   467     /// This function returns a const reference to an arc map storing
   468     /// the flow that is the union of the found arc-disjoint paths.
   469     ///
   470     /// \pre \ref run() or \ref findFlow() must be called before using
   471     /// this function.
   472     const FlowMap& flowMap() const {
   473       return *_flow;
   474     }
   475 
   476     /// \brief Return the potential of the given node.
   477     ///
   478     /// This function returns the potential of the given node.
   479     /// The node potentials provide the dual solution of the
   480     /// underlying \ref min_cost_flow "minimum cost flow problem".
   481     ///
   482     /// \pre \ref run() or \ref findFlow() must be called before using
   483     /// this function.
   484     Length potential(const Node& node) const {
   485       return (*_potential)[node];
   486     }
   487 
   488     /// \brief Return a const reference to a node map storing the
   489     /// found potentials (the dual solution).
   490     ///
   491     /// This function returns a const reference to a node map storing
   492     /// the found potentials that provide the dual solution of the
   493     /// underlying \ref min_cost_flow "minimum cost flow problem".
   494     ///
   495     /// \pre \ref run() or \ref findFlow() must be called before using
   496     /// this function.
   497     const PotentialMap& potentialMap() const {
   498       return *_potential;
   499     }
   500 
   501     /// \brief Return the number of the found paths.
   502     ///
   503     /// This function returns the number of the found paths.
   504     ///
   505     /// \pre \ref run() or \ref findFlow() must be called before using
   506     /// this function.
   507     int pathNum() const {
   508       return _path_num;
   509     }
   510 
   511     /// \brief Return a const reference to the specified path.
   512     ///
   513     /// This function returns a const reference to the specified path.
   514     ///
   515     /// \param i The function returns the <tt>i</tt>-th path.
   516     /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>.
   517     ///
   518     /// \pre \ref run() or \ref findPaths() must be called before using
   519     /// this function.
   520     const Path& path(int i) const {
   521       return paths[i];
   522     }
   523 
   524     /// @}
   525 
   526   }; //class Suurballe
   527 
   528   ///@}
   529 
   530 } //namespace lemon
   531 
   532 #endif //LEMON_SUURBALLE_H