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