lemon/suurballe.h
author Peter Kovacs <kpeter@inf.elte.hu>
Sun, 30 Nov 2008 00:51:20 +0100
changeset 394 e7707b3069f1
parent 345 2f64c4a692a8
child 425 cace3206223b
permissions -rw-r--r--
Better test files for Preflow (#176)

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