lemon/min_mean_cycle.h
 author Peter Kovacs Thu, 06 Aug 2009 20:31:04 +0200 changeset 809 03887b5e0f6f parent 808 5795860737f5 child 810 93cd93e82f9b permissions -rw-r--r--
Rename cyclePath() to cycle() in MinMeanCycle (#179)
     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_MIN_MEAN_CYCLE_H

    20 #define LEMON_MIN_MEAN_CYCLE_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 <lemon/core.h>

    29 #include <lemon/path.h>

    30 #include <lemon/tolerance.h>

    31 #include <lemon/connectivity.h>

    32

    33 namespace lemon {

    34

    35   /// \brief Default traits class of MinMeanCycle class.

    36   ///

    37   /// Default traits class of MinMeanCycle class.

    38   /// \tparam GR The type of the digraph.

    39   /// \tparam LEN The type of the length map.

    40   /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.

    41 #ifdef DOXYGEN

    42   template <typename GR, typename LEN>

    43 #else

    44   template <typename GR, typename LEN,

    45     bool integer = std::numeric_limits<typename LEN::Value>::is_integer>

    46 #endif

    47   struct MinMeanCycleDefaultTraits

    48   {

    49     /// The type of the digraph

    50     typedef GR Digraph;

    51     /// The type of the length map

    52     typedef LEN LengthMap;

    53     /// The type of the arc lengths

    54     typedef typename LengthMap::Value Value;

    55

    56     /// \brief The large value type used for internal computations

    57     ///

    58     /// The large value type used for internal computations.

    59     /// It is \c long \c long if the \c Value type is integer,

    60     /// otherwise it is \c double.

    61     /// \c Value must be convertible to \c LargeValue.

    62     typedef double LargeValue;

    63

    64     /// The tolerance type used for internal computations

    65     typedef lemon::Tolerance<LargeValue> Tolerance;

    66

    67     /// \brief The path type of the found cycles

    68     ///

    69     /// The path type of the found cycles.

    70     /// It must conform to the \ref lemon::concepts::Path "Path" concept

    71     /// and it must have an \c addBack() function.

    72     typedef lemon::Path<Digraph> Path;

    73   };

    74

    75   // Default traits class for integer value types

    76   template <typename GR, typename LEN>

    77   struct MinMeanCycleDefaultTraits<GR, LEN, true>

    78   {

    79     typedef GR Digraph;

    80     typedef LEN LengthMap;

    81     typedef typename LengthMap::Value Value;

    82 #ifdef LEMON_HAVE_LONG_LONG

    83     typedef long long LargeValue;

    84 #else

    85     typedef long LargeValue;

    86 #endif

    87     typedef lemon::Tolerance<LargeValue> Tolerance;

    88     typedef lemon::Path<Digraph> Path;

    89   };

    90

    91

    92   /// \addtogroup shortest_path

    93   /// @{

    94

    95   /// \brief Implementation of Howard's algorithm for finding a minimum

    96   /// mean cycle.

    97   ///

    98   /// \ref MinMeanCycle implements Howard's algorithm for finding a

    99   /// directed cycle of minimum mean length (cost) in a digraph.

   100   ///

   101   /// \tparam GR The type of the digraph the algorithm runs on.

   102   /// \tparam LEN The type of the length map. The default

   103   /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".

   104 #ifdef DOXYGEN

   105   template <typename GR, typename LEN, typename TR>

   106 #else

   107   template < typename GR,

   108              typename LEN = typename GR::template ArcMap<int>,

   109              typename TR = MinMeanCycleDefaultTraits<GR, LEN> >

   110 #endif

   111   class MinMeanCycle

   112   {

   113   public:

   114

   115     /// The type of the digraph

   116     typedef typename TR::Digraph Digraph;

   117     /// The type of the length map

   118     typedef typename TR::LengthMap LengthMap;

   119     /// The type of the arc lengths

   120     typedef typename TR::Value Value;

   121

   122     /// \brief The large value type

   123     ///

   124     /// The large value type used for internal computations.

   125     /// Using the \ref MinMeanCycleDefaultTraits "default traits class",

   126     /// it is \c long \c long if the \c Value type is integer,

   127     /// otherwise it is \c double.

   128     typedef typename TR::LargeValue LargeValue;

   129

   130     /// The tolerance type

   131     typedef typename TR::Tolerance Tolerance;

   132

   133     /// \brief The path type of the found cycles

   134     ///

   135     /// The path type of the found cycles.

   136     /// Using the \ref MinMeanCycleDefaultTraits "default traits class",

   137     /// it is \ref lemon::Path "Path<Digraph>".

   138     typedef typename TR::Path Path;

   139

   140     /// The \ref MinMeanCycleDefaultTraits "traits class" of the algorithm

   141     typedef TR Traits;

   142

   143   private:

   144

   145     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);

   146

   147     // The digraph the algorithm runs on

   148     const Digraph &_gr;

   149     // The length of the arcs

   150     const LengthMap &_length;

   151

   152     // Data for the found cycles

   153     bool _curr_found, _best_found;

   154     LargeValue _curr_length, _best_length;

   155     int _curr_size, _best_size;

   156     Node _curr_node, _best_node;

   157

   158     Path *_cycle_path;

   159     bool _local_path;

   160

   161     // Internal data used by the algorithm

   162     typename Digraph::template NodeMap<Arc> _policy;

   163     typename Digraph::template NodeMap<bool> _reached;

   164     typename Digraph::template NodeMap<int> _level;

   165     typename Digraph::template NodeMap<LargeValue> _dist;

   166

   167     // Data for storing the strongly connected components

   168     int _comp_num;

   169     typename Digraph::template NodeMap<int> _comp;

   170     std::vector<std::vector<Node> > _comp_nodes;

   171     std::vector<Node>* _nodes;

   172     typename Digraph::template NodeMap<std::vector<Arc> > _in_arcs;

   173

   174     // Queue used for BFS search

   175     std::vector<Node> _queue;

   176     int _qfront, _qback;

   177

   178     Tolerance _tolerance;

   179

   180   public:

   181

   182     /// \name Named Template Parameters

   183     /// @{

   184

   185     template <typename T>

   186     struct SetLargeValueTraits : public Traits {

   187       typedef T LargeValue;

   188       typedef lemon::Tolerance<T> Tolerance;

   189     };

   190

   191     /// \brief \ref named-templ-param "Named parameter" for setting

   192     /// \c LargeValue type.

   193     ///

   194     /// \ref named-templ-param "Named parameter" for setting \c LargeValue

   195     /// type. It is used for internal computations in the algorithm.

   196     template <typename T>

   197     struct SetLargeValue

   198       : public MinMeanCycle<GR, LEN, SetLargeValueTraits<T> > {

   199       typedef MinMeanCycle<GR, LEN, SetLargeValueTraits<T> > Create;

   200     };

   201

   202     template <typename T>

   203     struct SetPathTraits : public Traits {

   204       typedef T Path;

   205     };

   206

   207     /// \brief \ref named-templ-param "Named parameter" for setting

   208     /// \c %Path type.

   209     ///

   210     /// \ref named-templ-param "Named parameter" for setting the \c %Path

   211     /// type of the found cycles.

   212     /// It must conform to the \ref lemon::concepts::Path "Path" concept

   213     /// and it must have an \c addBack() function.

   214     template <typename T>

   215     struct SetPath

   216       : public MinMeanCycle<GR, LEN, SetPathTraits<T> > {

   217       typedef MinMeanCycle<GR, LEN, SetPathTraits<T> > Create;

   218     };

   219

   220     /// @}

   221

   222   public:

   223

   224     /// \brief Constructor.

   225     ///

   226     /// The constructor of the class.

   227     ///

   228     /// \param digraph The digraph the algorithm runs on.

   229     /// \param length The lengths (costs) of the arcs.

   230     MinMeanCycle( const Digraph &digraph,

   231                   const LengthMap &length ) :

   232       _gr(digraph), _length(length), _cycle_path(NULL), _local_path(false),

   233       _policy(digraph), _reached(digraph), _level(digraph), _dist(digraph),

   234       _comp(digraph), _in_arcs(digraph)

   235     {}

   236

   237     /// Destructor.

   238     ~MinMeanCycle() {

   239       if (_local_path) delete _cycle_path;

   240     }

   241

   242     /// \brief Set the path structure for storing the found cycle.

   243     ///

   244     /// This function sets an external path structure for storing the

   245     /// found cycle.

   246     ///

   247     /// If you don't call this function before calling \ref run() or

   248     /// \ref findMinMean(), it will allocate a local \ref Path "path"

   249     /// structure. The destuctor deallocates this automatically

   250     /// allocated object, of course.

   251     ///

   252     /// \note The algorithm calls only the \ref lemon::Path::addBack()

   253     /// "addBack()" function of the given path structure.

   254     ///

   255     /// \return <tt>(*this)</tt>

   256     MinMeanCycle& cycle(Path &path) {

   257       if (_local_path) {

   258         delete _cycle_path;

   259         _local_path = false;

   260       }

   261       _cycle_path = &path;

   262       return *this;

   263     }

   264

   265     /// \name Execution control

   266     /// The simplest way to execute the algorithm is to call the \ref run()

   267     /// function.\n

   268     /// If you only need the minimum mean length, you may call

   269     /// \ref findMinMean().

   270

   271     /// @{

   272

   273     /// \brief Run the algorithm.

   274     ///

   275     /// This function runs the algorithm.

   276     /// It can be called more than once (e.g. if the underlying digraph

   277     /// and/or the arc lengths have been modified).

   278     ///

   279     /// \return \c true if a directed cycle exists in the digraph.

   280     ///

   281     /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.

   282     /// \code

   283     ///   return mmc.findMinMean() && mmc.findCycle();

   284     /// \endcode

   285     bool run() {

   286       return findMinMean() && findCycle();

   287     }

   288

   289     /// \brief Find the minimum cycle mean.

   290     ///

   291     /// This function finds the minimum mean length of the directed

   292     /// cycles in the digraph.

   293     ///

   294     /// \return \c true if a directed cycle exists in the digraph.

   295     bool findMinMean() {

   296       // Initialize and find strongly connected components

   297       init();

   298       findComponents();

   299

   300       // Find the minimum cycle mean in the components

   301       for (int comp = 0; comp < _comp_num; ++comp) {

   302         // Find the minimum mean cycle in the current component

   303         if (!buildPolicyGraph(comp)) continue;

   304         while (true) {

   305           findPolicyCycle();

   306           if (!computeNodeDistances()) break;

   307         }

   308         // Update the best cycle (global minimum mean cycle)

   309         if ( !_best_found || (_curr_found &&

   310              _curr_length * _best_size < _best_length * _curr_size) ) {

   311           _best_found = true;

   312           _best_length = _curr_length;

   313           _best_size = _curr_size;

   314           _best_node = _curr_node;

   315         }

   316       }

   317       return _best_found;

   318     }

   319

   320     /// \brief Find a minimum mean directed cycle.

   321     ///

   322     /// This function finds a directed cycle of minimum mean length

   323     /// in the digraph using the data computed by findMinMean().

   324     ///

   325     /// \return \c true if a directed cycle exists in the digraph.

   326     ///

   327     /// \pre \ref findMinMean() must be called before using this function.

   328     bool findCycle() {

   329       if (!_best_found) return false;

   330       _cycle_path->addBack(_policy[_best_node]);

   331       for ( Node v = _best_node;

   332             (v = _gr.target(_policy[v])) != _best_node; ) {

   333         _cycle_path->addBack(_policy[v]);

   334       }

   335       return true;

   336     }

   337

   338     /// @}

   339

   340     /// \name Query Functions

   341     /// The results of the algorithm can be obtained using these

   342     /// functions.\n

   343     /// The algorithm should be executed before using them.

   344

   345     /// @{

   346

   347     /// \brief Return the total length of the found cycle.

   348     ///

   349     /// This function returns the total length of the found cycle.

   350     ///

   351     /// \pre \ref run() or \ref findMinMean() must be called before

   352     /// using this function.

   353     LargeValue cycleLength() const {

   354       return _best_length;

   355     }

   356

   357     /// \brief Return the number of arcs on the found cycle.

   358     ///

   359     /// This function returns the number of arcs on the found cycle.

   360     ///

   361     /// \pre \ref run() or \ref findMinMean() must be called before

   362     /// using this function.

   363     int cycleArcNum() const {

   364       return _best_size;

   365     }

   366

   367     /// \brief Return the mean length of the found cycle.

   368     ///

   369     /// This function returns the mean length of the found cycle.

   370     ///

   371     /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the

   372     /// following code.

   373     /// \code

   374     ///   return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum();

   375     /// \endcode

   376     ///

   377     /// \pre \ref run() or \ref findMinMean() must be called before

   378     /// using this function.

   379     double cycleMean() const {

   380       return static_cast<double>(_best_length) / _best_size;

   381     }

   382

   383     /// \brief Return the found cycle.

   384     ///

   385     /// This function returns a const reference to the path structure

   386     /// storing the found cycle.

   387     ///

   388     /// \pre \ref run() or \ref findCycle() must be called before using

   389     /// this function.

   390     const Path& cycle() const {

   391       return *_cycle_path;

   392     }

   393

   394     ///@}

   395

   396   private:

   397

   398     // Initialize

   399     void init() {

   400       if (!_cycle_path) {

   401         _local_path = true;

   402         _cycle_path = new Path;

   403       }

   404       _queue.resize(countNodes(_gr));

   405       _best_found = false;

   406       _best_length = 0;

   407       _best_size = 1;

   408       _cycle_path->clear();

   409     }

   410

   411     // Find strongly connected components and initialize _comp_nodes

   412     // and _in_arcs

   413     void findComponents() {

   414       _comp_num = stronglyConnectedComponents(_gr, _comp);

   415       _comp_nodes.resize(_comp_num);

   416       if (_comp_num == 1) {

   417         _comp_nodes[0].clear();

   418         for (NodeIt n(_gr); n != INVALID; ++n) {

   419           _comp_nodes[0].push_back(n);

   420           _in_arcs[n].clear();

   421           for (InArcIt a(_gr, n); a != INVALID; ++a) {

   422             _in_arcs[n].push_back(a);

   423           }

   424         }

   425       } else {

   426         for (int i = 0; i < _comp_num; ++i)

   427           _comp_nodes[i].clear();

   428         for (NodeIt n(_gr); n != INVALID; ++n) {

   429           int k = _comp[n];

   430           _comp_nodes[k].push_back(n);

   431           _in_arcs[n].clear();

   432           for (InArcIt a(_gr, n); a != INVALID; ++a) {

   433             if (_comp[_gr.source(a)] == k) _in_arcs[n].push_back(a);

   434           }

   435         }

   436       }

   437     }

   438

   439     // Build the policy graph in the given strongly connected component

   440     // (the out-degree of every node is 1)

   441     bool buildPolicyGraph(int comp) {

   442       _nodes = &(_comp_nodes[comp]);

   443       if (_nodes->size() < 1 ||

   444           (_nodes->size() == 1 && _in_arcs[(*_nodes)[0]].size() == 0)) {

   445         return false;

   446       }

   447       for (int i = 0; i < int(_nodes->size()); ++i) {

   448         _dist[(*_nodes)[i]] = std::numeric_limits<LargeValue>::max();

   449       }

   450       Node u, v;

   451       Arc e;

   452       for (int i = 0; i < int(_nodes->size()); ++i) {

   453         v = (*_nodes)[i];

   454         for (int j = 0; j < int(_in_arcs[v].size()); ++j) {

   455           e = _in_arcs[v][j];

   456           u = _gr.source(e);

   457           if (_length[e] < _dist[u]) {

   458             _dist[u] = _length[e];

   459             _policy[u] = e;

   460           }

   461         }

   462       }

   463       return true;

   464     }

   465

   466     // Find the minimum mean cycle in the policy graph

   467     void findPolicyCycle() {

   468       for (int i = 0; i < int(_nodes->size()); ++i) {

   469         _level[(*_nodes)[i]] = -1;

   470       }

   471       LargeValue clength;

   472       int csize;

   473       Node u, v;

   474       _curr_found = false;

   475       for (int i = 0; i < int(_nodes->size()); ++i) {

   476         u = (*_nodes)[i];

   477         if (_level[u] >= 0) continue;

   478         for (; _level[u] < 0; u = _gr.target(_policy[u])) {

   479           _level[u] = i;

   480         }

   481         if (_level[u] == i) {

   482           // A cycle is found

   483           clength = _length[_policy[u]];

   484           csize = 1;

   485           for (v = u; (v = _gr.target(_policy[v])) != u; ) {

   486             clength += _length[_policy[v]];

   487             ++csize;

   488           }

   489           if ( !_curr_found ||

   490                (clength * _curr_size < _curr_length * csize) ) {

   491             _curr_found = true;

   492             _curr_length = clength;

   493             _curr_size = csize;

   494             _curr_node = u;

   495           }

   496         }

   497       }

   498     }

   499

   500     // Contract the policy graph and compute node distances

   501     bool computeNodeDistances() {

   502       // Find the component of the main cycle and compute node distances

   503       // using reverse BFS

   504       for (int i = 0; i < int(_nodes->size()); ++i) {

   505         _reached[(*_nodes)[i]] = false;

   506       }

   507       _qfront = _qback = 0;

   508       _queue[0] = _curr_node;

   509       _reached[_curr_node] = true;

   510       _dist[_curr_node] = 0;

   511       Node u, v;

   512       Arc e;

   513       while (_qfront <= _qback) {

   514         v = _queue[_qfront++];

   515         for (int j = 0; j < int(_in_arcs[v].size()); ++j) {

   516           e = _in_arcs[v][j];

   517           u = _gr.source(e);

   518           if (_policy[u] == e && !_reached[u]) {

   519             _reached[u] = true;

   520             _dist[u] = _dist[v] + _length[e] * _curr_size - _curr_length;

   521             _queue[++_qback] = u;

   522           }

   523         }

   524       }

   525

   526       // Connect all other nodes to this component and compute node

   527       // distances using reverse BFS

   528       _qfront = 0;

   529       while (_qback < int(_nodes->size())-1) {

   530         v = _queue[_qfront++];

   531         for (int j = 0; j < int(_in_arcs[v].size()); ++j) {

   532           e = _in_arcs[v][j];

   533           u = _gr.source(e);

   534           if (!_reached[u]) {

   535             _reached[u] = true;

   536             _policy[u] = e;

   537             _dist[u] = _dist[v] + _length[e] * _curr_size - _curr_length;

   538             _queue[++_qback] = u;

   539           }

   540         }

   541       }

   542

   543       // Improve node distances

   544       bool improved = false;

   545       for (int i = 0; i < int(_nodes->size()); ++i) {

   546         v = (*_nodes)[i];

   547         for (int j = 0; j < int(_in_arcs[v].size()); ++j) {

   548           e = _in_arcs[v][j];

   549           u = _gr.source(e);

   550           LargeValue delta = _dist[v] + _length[e] * _curr_size - _curr_length;

   551           if (_tolerance.less(delta, _dist[u])) {

   552             _dist[u] = delta;

   553             _policy[u] = e;

   554             improved = true;

   555           }

   556         }

   557       }

   558       return improved;

   559     }

   560

   561   }; //class MinMeanCycle

   562

   563   ///@}

   564

   565 } //namespace lemon

   566

   567 #endif //LEMON_MIN_MEAN_CYCLE_H