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