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