lemon/howard.h
author Peter Kovacs <kpeter@inf.elte.hu>
Tue, 11 Aug 2009 21:53:39 +0200
changeset 766 97744b6dabf8
parent 763 93cd93e82f9b
child 767 11c946fa8d13
permissions -rw-r--r--
Add HartmannOrlin algorithm class (#179)
This algorithm is an improved version of Karp's original method,
it applies an efficient early termination scheme.
The interface is the same as Karp's and Howard's interface.
     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_HOWARD_H
    20 #define LEMON_HOWARD_H
    21 
    22 /// \ingroup shortest_path
    23 ///
    24 /// \file
    25 /// \brief Howard's algorithm for finding a minimum mean cycle.
    26 
    27 #include <vector>
    28 #include <limits>
    29 #include <lemon/core.h>
    30 #include <lemon/path.h>
    31 #include <lemon/tolerance.h>
    32 #include <lemon/connectivity.h>
    33 
    34 namespace lemon {
    35 
    36   /// \brief Default traits class of Howard class.
    37   ///
    38   /// Default traits class of Howard class.
    39   /// \tparam GR The type of the digraph.
    40   /// \tparam LEN The type of the length map.
    41   /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
    42 #ifdef DOXYGEN
    43   template <typename GR, typename LEN>
    44 #else
    45   template <typename GR, typename LEN,
    46     bool integer = std::numeric_limits<typename LEN::Value>::is_integer>
    47 #endif
    48   struct HowardDefaultTraits
    49   {
    50     /// The type of the digraph
    51     typedef GR Digraph;
    52     /// The type of the length map
    53     typedef LEN LengthMap;
    54     /// The type of the arc lengths
    55     typedef typename LengthMap::Value Value;
    56 
    57     /// \brief The large value type used for internal computations
    58     ///
    59     /// The large value type used for internal computations.
    60     /// It is \c long \c long if the \c Value type is integer,
    61     /// otherwise it is \c double.
    62     /// \c Value must be convertible to \c LargeValue.
    63     typedef double LargeValue;
    64 
    65     /// The tolerance type used for internal computations
    66     typedef lemon::Tolerance<LargeValue> Tolerance;
    67 
    68     /// \brief The path type of the found cycles
    69     ///
    70     /// The path type of the found cycles.
    71     /// It must conform to the \ref lemon::concepts::Path "Path" concept
    72     /// and it must have an \c addBack() function.
    73     typedef lemon::Path<Digraph> Path;
    74   };
    75 
    76   // Default traits class for integer value types
    77   template <typename GR, typename LEN>
    78   struct HowardDefaultTraits<GR, LEN, true>
    79   {
    80     typedef GR Digraph;
    81     typedef LEN LengthMap;
    82     typedef typename LengthMap::Value Value;
    83 #ifdef LEMON_HAVE_LONG_LONG
    84     typedef long long LargeValue;
    85 #else
    86     typedef long LargeValue;
    87 #endif
    88     typedef lemon::Tolerance<LargeValue> Tolerance;
    89     typedef lemon::Path<Digraph> Path;
    90   };
    91 
    92 
    93   /// \addtogroup shortest_path
    94   /// @{
    95 
    96   /// \brief Implementation of Howard's algorithm for finding a minimum
    97   /// mean cycle.
    98   ///
    99   /// This class implements Howard's policy iteration algorithm for finding
   100   /// a directed cycle of minimum mean length (cost) in a digraph.
   101   ///
   102   /// \tparam GR The type of the digraph the algorithm runs on.
   103   /// \tparam LEN The type of the length map. The default
   104   /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
   105 #ifdef DOXYGEN
   106   template <typename GR, typename LEN, typename TR>
   107 #else
   108   template < typename GR,
   109              typename LEN = typename GR::template ArcMap<int>,
   110              typename TR = HowardDefaultTraits<GR, LEN> >
   111 #endif
   112   class Howard
   113   {
   114   public:
   115   
   116     /// The type of the digraph
   117     typedef typename TR::Digraph Digraph;
   118     /// The type of the length map
   119     typedef typename TR::LengthMap LengthMap;
   120     /// The type of the arc lengths
   121     typedef typename TR::Value Value;
   122 
   123     /// \brief The large value type
   124     ///
   125     /// The large value type used for internal computations.
   126     /// Using the \ref HowardDefaultTraits "default traits class",
   127     /// it is \c long \c long if the \c Value type is integer,
   128     /// otherwise it is \c double.
   129     typedef typename TR::LargeValue LargeValue;
   130 
   131     /// The tolerance type
   132     typedef typename TR::Tolerance Tolerance;
   133 
   134     /// \brief The path type of the found cycles
   135     ///
   136     /// The path type of the found cycles.
   137     /// Using the \ref HowardDefaultTraits "default traits class",
   138     /// it is \ref lemon::Path "Path<Digraph>".
   139     typedef typename TR::Path Path;
   140 
   141     /// The \ref HowardDefaultTraits "traits class" of the algorithm
   142     typedef TR Traits;
   143 
   144   private:
   145 
   146     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
   147   
   148     // The digraph the algorithm runs on
   149     const Digraph &_gr;
   150     // The length of the arcs
   151     const LengthMap &_length;
   152 
   153     // Data for the found cycles
   154     bool _curr_found, _best_found;
   155     LargeValue _curr_length, _best_length;
   156     int _curr_size, _best_size;
   157     Node _curr_node, _best_node;
   158 
   159     Path *_cycle_path;
   160     bool _local_path;
   161 
   162     // Internal data used by the algorithm
   163     typename Digraph::template NodeMap<Arc> _policy;
   164     typename Digraph::template NodeMap<bool> _reached;
   165     typename Digraph::template NodeMap<int> _level;
   166     typename Digraph::template NodeMap<LargeValue> _dist;
   167 
   168     // Data for storing the strongly connected components
   169     int _comp_num;
   170     typename Digraph::template NodeMap<int> _comp;
   171     std::vector<std::vector<Node> > _comp_nodes;
   172     std::vector<Node>* _nodes;
   173     typename Digraph::template NodeMap<std::vector<Arc> > _in_arcs;
   174     
   175     // Queue used for BFS search
   176     std::vector<Node> _queue;
   177     int _qfront, _qback;
   178 
   179     Tolerance _tolerance;
   180   
   181   public:
   182   
   183     /// \name Named Template Parameters
   184     /// @{
   185 
   186     template <typename T>
   187     struct SetLargeValueTraits : public Traits {
   188       typedef T LargeValue;
   189       typedef lemon::Tolerance<T> Tolerance;
   190     };
   191 
   192     /// \brief \ref named-templ-param "Named parameter" for setting
   193     /// \c LargeValue type.
   194     ///
   195     /// \ref named-templ-param "Named parameter" for setting \c LargeValue
   196     /// type. It is used for internal computations in the algorithm.
   197     template <typename T>
   198     struct SetLargeValue
   199       : public Howard<GR, LEN, SetLargeValueTraits<T> > {
   200       typedef Howard<GR, LEN, SetLargeValueTraits<T> > Create;
   201     };
   202 
   203     template <typename T>
   204     struct SetPathTraits : public Traits {
   205       typedef T Path;
   206     };
   207 
   208     /// \brief \ref named-templ-param "Named parameter" for setting
   209     /// \c %Path type.
   210     ///
   211     /// \ref named-templ-param "Named parameter" for setting the \c %Path
   212     /// type of the found cycles.
   213     /// It must conform to the \ref lemon::concepts::Path "Path" concept
   214     /// and it must have an \c addBack() function.
   215     template <typename T>
   216     struct SetPath
   217       : public Howard<GR, LEN, SetPathTraits<T> > {
   218       typedef Howard<GR, LEN, SetPathTraits<T> > Create;
   219     };
   220     
   221     /// @}
   222 
   223   public:
   224 
   225     /// \brief Constructor.
   226     ///
   227     /// The constructor of the class.
   228     ///
   229     /// \param digraph The digraph the algorithm runs on.
   230     /// \param length The lengths (costs) of the arcs.
   231     Howard( const Digraph &digraph,
   232             const LengthMap &length ) :
   233       _gr(digraph), _length(length), _cycle_path(NULL), _local_path(false),
   234       _policy(digraph), _reached(digraph), _level(digraph), _dist(digraph),
   235       _comp(digraph), _in_arcs(digraph)
   236     {}
   237 
   238     /// Destructor.
   239     ~Howard() {
   240       if (_local_path) delete _cycle_path;
   241     }
   242 
   243     /// \brief Set the path structure for storing the found cycle.
   244     ///
   245     /// This function sets an external path structure for storing the
   246     /// found cycle.
   247     ///
   248     /// If you don't call this function before calling \ref run() or
   249     /// \ref findMinMean(), it will allocate a local \ref Path "path"
   250     /// structure. The destuctor deallocates this automatically
   251     /// allocated object, of course.
   252     ///
   253     /// \note The algorithm calls only the \ref lemon::Path::addBack()
   254     /// "addBack()" function of the given path structure.
   255     ///
   256     /// \return <tt>(*this)</tt>
   257     Howard& cycle(Path &path) {
   258       if (_local_path) {
   259         delete _cycle_path;
   260         _local_path = false;
   261       }
   262       _cycle_path = &path;
   263       return *this;
   264     }
   265 
   266     /// \name Execution control
   267     /// The simplest way to execute the algorithm is to call the \ref run()
   268     /// function.\n
   269     /// If you only need the minimum mean length, you may call
   270     /// \ref findMinMean().
   271 
   272     /// @{
   273 
   274     /// \brief Run the algorithm.
   275     ///
   276     /// This function runs the algorithm.
   277     /// It can be called more than once (e.g. if the underlying digraph
   278     /// and/or the arc lengths have been modified).
   279     ///
   280     /// \return \c true if a directed cycle exists in the digraph.
   281     ///
   282     /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
   283     /// \code
   284     ///   return mmc.findMinMean() && mmc.findCycle();
   285     /// \endcode
   286     bool run() {
   287       return findMinMean() && findCycle();
   288     }
   289 
   290     /// \brief Find the minimum cycle mean.
   291     ///
   292     /// This function finds the minimum mean length of the directed
   293     /// cycles in the digraph.
   294     ///
   295     /// \return \c true if a directed cycle exists in the digraph.
   296     bool findMinMean() {
   297       // Initialize and find strongly connected components
   298       init();
   299       findComponents();
   300       
   301       // Find the minimum cycle mean in the components
   302       for (int comp = 0; comp < _comp_num; ++comp) {
   303         // Find the minimum mean cycle in the current component
   304         if (!buildPolicyGraph(comp)) continue;
   305         while (true) {
   306           findPolicyCycle();
   307           if (!computeNodeDistances()) break;
   308         }
   309         // Update the best cycle (global minimum mean cycle)
   310         if ( !_best_found || (_curr_found &&
   311              _curr_length * _best_size < _best_length * _curr_size) ) {
   312           _best_found = true;
   313           _best_length = _curr_length;
   314           _best_size = _curr_size;
   315           _best_node = _curr_node;
   316         }
   317       }
   318       return _best_found;
   319     }
   320 
   321     /// \brief Find a minimum mean directed cycle.
   322     ///
   323     /// This function finds a directed cycle of minimum mean length
   324     /// in the digraph using the data computed by findMinMean().
   325     ///
   326     /// \return \c true if a directed cycle exists in the digraph.
   327     ///
   328     /// \pre \ref findMinMean() must be called before using this function.
   329     bool findCycle() {
   330       if (!_best_found) return false;
   331       _cycle_path->addBack(_policy[_best_node]);
   332       for ( Node v = _best_node;
   333             (v = _gr.target(_policy[v])) != _best_node; ) {
   334         _cycle_path->addBack(_policy[v]);
   335       }
   336       return true;
   337     }
   338 
   339     /// @}
   340 
   341     /// \name Query Functions
   342     /// The results of the algorithm can be obtained using these
   343     /// functions.\n
   344     /// The algorithm should be executed before using them.
   345 
   346     /// @{
   347 
   348     /// \brief Return the total length of the found cycle.
   349     ///
   350     /// This function returns the total length of the found cycle.
   351     ///
   352     /// \pre \ref run() or \ref findMinMean() must be called before
   353     /// using this function.
   354     LargeValue cycleLength() const {
   355       return _best_length;
   356     }
   357 
   358     /// \brief Return the number of arcs on the found cycle.
   359     ///
   360     /// This function returns the number of arcs on the found cycle.
   361     ///
   362     /// \pre \ref run() or \ref findMinMean() must be called before
   363     /// using this function.
   364     int cycleArcNum() const {
   365       return _best_size;
   366     }
   367 
   368     /// \brief Return the mean length of the found cycle.
   369     ///
   370     /// This function returns the mean length of the found cycle.
   371     ///
   372     /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
   373     /// following code.
   374     /// \code
   375     ///   return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum();
   376     /// \endcode
   377     ///
   378     /// \pre \ref run() or \ref findMinMean() must be called before
   379     /// using this function.
   380     double cycleMean() const {
   381       return static_cast<double>(_best_length) / _best_size;
   382     }
   383 
   384     /// \brief Return the found cycle.
   385     ///
   386     /// This function returns a const reference to the path structure
   387     /// storing the found cycle.
   388     ///
   389     /// \pre \ref run() or \ref findCycle() must be called before using
   390     /// this function.
   391     const Path& cycle() const {
   392       return *_cycle_path;
   393     }
   394 
   395     ///@}
   396 
   397   private:
   398 
   399     // Initialize
   400     void init() {
   401       if (!_cycle_path) {
   402         _local_path = true;
   403         _cycle_path = new Path;
   404       }
   405       _queue.resize(countNodes(_gr));
   406       _best_found = false;
   407       _best_length = 0;
   408       _best_size = 1;
   409       _cycle_path->clear();
   410     }
   411     
   412     // Find strongly connected components and initialize _comp_nodes
   413     // and _in_arcs
   414     void findComponents() {
   415       _comp_num = stronglyConnectedComponents(_gr, _comp);
   416       _comp_nodes.resize(_comp_num);
   417       if (_comp_num == 1) {
   418         _comp_nodes[0].clear();
   419         for (NodeIt n(_gr); n != INVALID; ++n) {
   420           _comp_nodes[0].push_back(n);
   421           _in_arcs[n].clear();
   422           for (InArcIt a(_gr, n); a != INVALID; ++a) {
   423             _in_arcs[n].push_back(a);
   424           }
   425         }
   426       } else {
   427         for (int i = 0; i < _comp_num; ++i)
   428           _comp_nodes[i].clear();
   429         for (NodeIt n(_gr); n != INVALID; ++n) {
   430           int k = _comp[n];
   431           _comp_nodes[k].push_back(n);
   432           _in_arcs[n].clear();
   433           for (InArcIt a(_gr, n); a != INVALID; ++a) {
   434             if (_comp[_gr.source(a)] == k) _in_arcs[n].push_back(a);
   435           }
   436         }
   437       }
   438     }
   439 
   440     // Build the policy graph in the given strongly connected component
   441     // (the out-degree of every node is 1)
   442     bool buildPolicyGraph(int comp) {
   443       _nodes = &(_comp_nodes[comp]);
   444       if (_nodes->size() < 1 ||
   445           (_nodes->size() == 1 && _in_arcs[(*_nodes)[0]].size() == 0)) {
   446         return false;
   447       }
   448       for (int i = 0; i < int(_nodes->size()); ++i) {
   449         _dist[(*_nodes)[i]] = std::numeric_limits<LargeValue>::max();
   450       }
   451       Node u, v;
   452       Arc e;
   453       for (int i = 0; i < int(_nodes->size()); ++i) {
   454         v = (*_nodes)[i];
   455         for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
   456           e = _in_arcs[v][j];
   457           u = _gr.source(e);
   458           if (_length[e] < _dist[u]) {
   459             _dist[u] = _length[e];
   460             _policy[u] = e;
   461           }
   462         }
   463       }
   464       return true;
   465     }
   466 
   467     // Find the minimum mean cycle in the policy graph
   468     void findPolicyCycle() {
   469       for (int i = 0; i < int(_nodes->size()); ++i) {
   470         _level[(*_nodes)[i]] = -1;
   471       }
   472       LargeValue clength;
   473       int csize;
   474       Node u, v;
   475       _curr_found = false;
   476       for (int i = 0; i < int(_nodes->size()); ++i) {
   477         u = (*_nodes)[i];
   478         if (_level[u] >= 0) continue;
   479         for (; _level[u] < 0; u = _gr.target(_policy[u])) {
   480           _level[u] = i;
   481         }
   482         if (_level[u] == i) {
   483           // A cycle is found
   484           clength = _length[_policy[u]];
   485           csize = 1;
   486           for (v = u; (v = _gr.target(_policy[v])) != u; ) {
   487             clength += _length[_policy[v]];
   488             ++csize;
   489           }
   490           if ( !_curr_found ||
   491                (clength * _curr_size < _curr_length * csize) ) {
   492             _curr_found = true;
   493             _curr_length = clength;
   494             _curr_size = csize;
   495             _curr_node = u;
   496           }
   497         }
   498       }
   499     }
   500 
   501     // Contract the policy graph and compute node distances
   502     bool computeNodeDistances() {
   503       // Find the component of the main cycle and compute node distances
   504       // using reverse BFS
   505       for (int i = 0; i < int(_nodes->size()); ++i) {
   506         _reached[(*_nodes)[i]] = false;
   507       }
   508       _qfront = _qback = 0;
   509       _queue[0] = _curr_node;
   510       _reached[_curr_node] = true;
   511       _dist[_curr_node] = 0;
   512       Node u, v;
   513       Arc e;
   514       while (_qfront <= _qback) {
   515         v = _queue[_qfront++];
   516         for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
   517           e = _in_arcs[v][j];
   518           u = _gr.source(e);
   519           if (_policy[u] == e && !_reached[u]) {
   520             _reached[u] = true;
   521             _dist[u] = _dist[v] + _length[e] * _curr_size - _curr_length;
   522             _queue[++_qback] = u;
   523           }
   524         }
   525       }
   526 
   527       // Connect all other nodes to this component and compute node
   528       // distances using reverse BFS
   529       _qfront = 0;
   530       while (_qback < int(_nodes->size())-1) {
   531         v = _queue[_qfront++];
   532         for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
   533           e = _in_arcs[v][j];
   534           u = _gr.source(e);
   535           if (!_reached[u]) {
   536             _reached[u] = true;
   537             _policy[u] = e;
   538             _dist[u] = _dist[v] + _length[e] * _curr_size - _curr_length;
   539             _queue[++_qback] = u;
   540           }
   541         }
   542       }
   543 
   544       // Improve node distances
   545       bool improved = false;
   546       for (int i = 0; i < int(_nodes->size()); ++i) {
   547         v = (*_nodes)[i];
   548         for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
   549           e = _in_arcs[v][j];
   550           u = _gr.source(e);
   551           LargeValue delta = _dist[v] + _length[e] * _curr_size - _curr_length;
   552           if (_tolerance.less(delta, _dist[u])) {
   553             _dist[u] = delta;
   554             _policy[u] = e;
   555             improved = true;
   556           }
   557         }
   558       }
   559       return improved;
   560     }
   561 
   562   }; //class Howard
   563 
   564   ///@}
   565 
   566 } //namespace lemon
   567 
   568 #endif //LEMON_HOWARD_H