lemon/karp.h
changeset 766 97744b6dabf8
child 767 11c946fa8d13
equal deleted inserted replaced
-1:000000000000 0:5970d96a49ae
       
     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_KARP_H
       
    20 #define LEMON_KARP_H
       
    21 
       
    22 /// \ingroup shortest_path
       
    23 ///
       
    24 /// \file
       
    25 /// \brief Karp'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 Karp algorithm.
       
    37   ///
       
    38   /// Default traits class of Karp algorithm.
       
    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 KarpDefaultTraits
       
    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 KarpDefaultTraits<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 Karp's algorithm for finding a minimum
       
    97   /// mean cycle.
       
    98   ///
       
    99   /// This class implements Karp's algorithm for finding a directed
       
   100   /// 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 = KarpDefaultTraits<GR, LEN> >
       
   111 #endif
       
   112   class Karp
       
   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 KarpDefaultTraits "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 KarpDefaultTraits "default traits class",
       
   138     /// it is \ref lemon::Path "Path<Digraph>".
       
   139     typedef typename TR::Path Path;
       
   140 
       
   141     /// The \ref KarpDefaultTraits "traits class" of the algorithm
       
   142     typedef TR Traits;
       
   143 
       
   144   private:
       
   145 
       
   146     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
       
   147 
       
   148     // Data sturcture for path data
       
   149     struct PathData
       
   150     {
       
   151       bool found;
       
   152       LargeValue dist;
       
   153       Arc pred;
       
   154       PathData(bool f = false, LargeValue d = 0, Arc p = INVALID) :
       
   155         found(f), dist(d), pred(p) {}
       
   156     };
       
   157 
       
   158     typedef typename Digraph::template NodeMap<std::vector<PathData> >
       
   159       PathDataNodeMap;
       
   160 
       
   161   private:
       
   162 
       
   163     // The digraph the algorithm runs on
       
   164     const Digraph &_gr;
       
   165     // The length of the arcs
       
   166     const LengthMap &_length;
       
   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> > _out_arcs;
       
   174 
       
   175     // Data for the found cycle
       
   176     LargeValue _cycle_length;
       
   177     int _cycle_size;
       
   178     Node _cycle_node;
       
   179 
       
   180     Path *_cycle_path;
       
   181     bool _local_path;
       
   182 
       
   183     // Node map for storing path data
       
   184     PathDataNodeMap _data;
       
   185     // The processed nodes in the last round
       
   186     std::vector<Node> _process;
       
   187 
       
   188     Tolerance _tolerance;
       
   189 
       
   190   public:
       
   191 
       
   192     /// \name Named Template Parameters
       
   193     /// @{
       
   194 
       
   195     template <typename T>
       
   196     struct SetLargeValueTraits : public Traits {
       
   197       typedef T LargeValue;
       
   198       typedef lemon::Tolerance<T> Tolerance;
       
   199     };
       
   200 
       
   201     /// \brief \ref named-templ-param "Named parameter" for setting
       
   202     /// \c LargeValue type.
       
   203     ///
       
   204     /// \ref named-templ-param "Named parameter" for setting \c LargeValue
       
   205     /// type. It is used for internal computations in the algorithm.
       
   206     template <typename T>
       
   207     struct SetLargeValue
       
   208       : public Karp<GR, LEN, SetLargeValueTraits<T> > {
       
   209       typedef Karp<GR, LEN, SetLargeValueTraits<T> > Create;
       
   210     };
       
   211 
       
   212     template <typename T>
       
   213     struct SetPathTraits : public Traits {
       
   214       typedef T Path;
       
   215     };
       
   216 
       
   217     /// \brief \ref named-templ-param "Named parameter" for setting
       
   218     /// \c %Path type.
       
   219     ///
       
   220     /// \ref named-templ-param "Named parameter" for setting the \c %Path
       
   221     /// type of the found cycles.
       
   222     /// It must conform to the \ref lemon::concepts::Path "Path" concept
       
   223     /// and it must have an \c addFront() function.
       
   224     template <typename T>
       
   225     struct SetPath
       
   226       : public Karp<GR, LEN, SetPathTraits<T> > {
       
   227       typedef Karp<GR, LEN, SetPathTraits<T> > Create;
       
   228     };
       
   229 
       
   230     /// @}
       
   231 
       
   232   public:
       
   233 
       
   234     /// \brief Constructor.
       
   235     ///
       
   236     /// The constructor of the class.
       
   237     ///
       
   238     /// \param digraph The digraph the algorithm runs on.
       
   239     /// \param length The lengths (costs) of the arcs.
       
   240     Karp( const Digraph &digraph,
       
   241           const LengthMap &length ) :
       
   242       _gr(digraph), _length(length), _comp(digraph), _out_arcs(digraph),
       
   243       _cycle_length(0), _cycle_size(1), _cycle_node(INVALID),
       
   244       _cycle_path(NULL), _local_path(false), _data(digraph)
       
   245     {}
       
   246 
       
   247     /// Destructor.
       
   248     ~Karp() {
       
   249       if (_local_path) delete _cycle_path;
       
   250     }
       
   251 
       
   252     /// \brief Set the path structure for storing the found cycle.
       
   253     ///
       
   254     /// This function sets an external path structure for storing the
       
   255     /// found cycle.
       
   256     ///
       
   257     /// If you don't call this function before calling \ref run() or
       
   258     /// \ref findMinMean(), it will allocate a local \ref Path "path"
       
   259     /// structure. The destuctor deallocates this automatically
       
   260     /// allocated object, of course.
       
   261     ///
       
   262     /// \note The algorithm calls only the \ref lemon::Path::addFront()
       
   263     /// "addFront()" function of the given path structure.
       
   264     ///
       
   265     /// \return <tt>(*this)</tt>
       
   266     Karp& cycle(Path &path) {
       
   267       if (_local_path) {
       
   268         delete _cycle_path;
       
   269         _local_path = false;
       
   270       }
       
   271       _cycle_path = &path;
       
   272       return *this;
       
   273     }
       
   274 
       
   275     /// \name Execution control
       
   276     /// The simplest way to execute the algorithm is to call the \ref run()
       
   277     /// function.\n
       
   278     /// If you only need the minimum mean length, you may call
       
   279     /// \ref findMinMean().
       
   280 
       
   281     /// @{
       
   282 
       
   283     /// \brief Run the algorithm.
       
   284     ///
       
   285     /// This function runs the algorithm.
       
   286     /// It can be called more than once (e.g. if the underlying digraph
       
   287     /// and/or the arc lengths have been modified).
       
   288     ///
       
   289     /// \return \c true if a directed cycle exists in the digraph.
       
   290     ///
       
   291     /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
       
   292     /// \code
       
   293     ///   return mmc.findMinMean() && mmc.findCycle();
       
   294     /// \endcode
       
   295     bool run() {
       
   296       return findMinMean() && findCycle();
       
   297     }
       
   298 
       
   299     /// \brief Find the minimum cycle mean.
       
   300     ///
       
   301     /// This function finds the minimum mean length of the directed
       
   302     /// cycles in the digraph.
       
   303     ///
       
   304     /// \return \c true if a directed cycle exists in the digraph.
       
   305     bool findMinMean() {
       
   306       // Initialization and find strongly connected components
       
   307       init();
       
   308       findComponents();
       
   309       
       
   310       // Find the minimum cycle mean in the components
       
   311       for (int comp = 0; comp < _comp_num; ++comp) {
       
   312         if (!initComponent(comp)) continue;
       
   313         processRounds();
       
   314         updateMinMean();
       
   315       }
       
   316       return (_cycle_node != INVALID);
       
   317     }
       
   318 
       
   319     /// \brief Find a minimum mean directed cycle.
       
   320     ///
       
   321     /// This function finds a directed cycle of minimum mean length
       
   322     /// in the digraph using the data computed by findMinMean().
       
   323     ///
       
   324     /// \return \c true if a directed cycle exists in the digraph.
       
   325     ///
       
   326     /// \pre \ref findMinMean() must be called before using this function.
       
   327     bool findCycle() {
       
   328       if (_cycle_node == INVALID) return false;
       
   329       IntNodeMap reached(_gr, -1);
       
   330       int r = _data[_cycle_node].size();
       
   331       Node u = _cycle_node;
       
   332       while (reached[u] < 0) {
       
   333         reached[u] = --r;
       
   334         u = _gr.source(_data[u][r].pred);
       
   335       }
       
   336       r = reached[u];
       
   337       Arc e = _data[u][r].pred;
       
   338       _cycle_path->addFront(e);
       
   339       _cycle_length = _length[e];
       
   340       _cycle_size = 1;
       
   341       Node v;
       
   342       while ((v = _gr.source(e)) != u) {
       
   343         e = _data[v][--r].pred;
       
   344         _cycle_path->addFront(e);
       
   345         _cycle_length += _length[e];
       
   346         ++_cycle_size;
       
   347       }
       
   348       return true;
       
   349     }
       
   350 
       
   351     /// @}
       
   352 
       
   353     /// \name Query Functions
       
   354     /// The results of the algorithm can be obtained using these
       
   355     /// functions.\n
       
   356     /// The algorithm should be executed before using them.
       
   357 
       
   358     /// @{
       
   359 
       
   360     /// \brief Return the total length of the found cycle.
       
   361     ///
       
   362     /// This function returns the total length of the found cycle.
       
   363     ///
       
   364     /// \pre \ref run() or \ref findMinMean() must be called before
       
   365     /// using this function.
       
   366     LargeValue cycleLength() const {
       
   367       return _cycle_length;
       
   368     }
       
   369 
       
   370     /// \brief Return the number of arcs on the found cycle.
       
   371     ///
       
   372     /// This function returns the number of arcs on the found cycle.
       
   373     ///
       
   374     /// \pre \ref run() or \ref findMinMean() must be called before
       
   375     /// using this function.
       
   376     int cycleArcNum() const {
       
   377       return _cycle_size;
       
   378     }
       
   379 
       
   380     /// \brief Return the mean length of the found cycle.
       
   381     ///
       
   382     /// This function returns the mean length of the found cycle.
       
   383     ///
       
   384     /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
       
   385     /// following code.
       
   386     /// \code
       
   387     ///   return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum();
       
   388     /// \endcode
       
   389     ///
       
   390     /// \pre \ref run() or \ref findMinMean() must be called before
       
   391     /// using this function.
       
   392     double cycleMean() const {
       
   393       return static_cast<double>(_cycle_length) / _cycle_size;
       
   394     }
       
   395 
       
   396     /// \brief Return the found cycle.
       
   397     ///
       
   398     /// This function returns a const reference to the path structure
       
   399     /// storing the found cycle.
       
   400     ///
       
   401     /// \pre \ref run() or \ref findCycle() must be called before using
       
   402     /// this function.
       
   403     const Path& cycle() const {
       
   404       return *_cycle_path;
       
   405     }
       
   406 
       
   407     ///@}
       
   408 
       
   409   private:
       
   410 
       
   411     // Initialization
       
   412     void init() {
       
   413       if (!_cycle_path) {
       
   414         _local_path = true;
       
   415         _cycle_path = new Path;
       
   416       }
       
   417       _cycle_path->clear();
       
   418       _cycle_length = 0;
       
   419       _cycle_size = 1;
       
   420       _cycle_node = INVALID;
       
   421       for (NodeIt u(_gr); u != INVALID; ++u)
       
   422         _data[u].clear();
       
   423     }
       
   424 
       
   425     // Find strongly connected components and initialize _comp_nodes
       
   426     // and _out_arcs
       
   427     void findComponents() {
       
   428       _comp_num = stronglyConnectedComponents(_gr, _comp);
       
   429       _comp_nodes.resize(_comp_num);
       
   430       if (_comp_num == 1) {
       
   431         _comp_nodes[0].clear();
       
   432         for (NodeIt n(_gr); n != INVALID; ++n) {
       
   433           _comp_nodes[0].push_back(n);
       
   434           _out_arcs[n].clear();
       
   435           for (OutArcIt a(_gr, n); a != INVALID; ++a) {
       
   436             _out_arcs[n].push_back(a);
       
   437           }
       
   438         }
       
   439       } else {
       
   440         for (int i = 0; i < _comp_num; ++i)
       
   441           _comp_nodes[i].clear();
       
   442         for (NodeIt n(_gr); n != INVALID; ++n) {
       
   443           int k = _comp[n];
       
   444           _comp_nodes[k].push_back(n);
       
   445           _out_arcs[n].clear();
       
   446           for (OutArcIt a(_gr, n); a != INVALID; ++a) {
       
   447             if (_comp[_gr.target(a)] == k) _out_arcs[n].push_back(a);
       
   448           }
       
   449         }
       
   450       }
       
   451     }
       
   452 
       
   453     // Initialize path data for the current component
       
   454     bool initComponent(int comp) {
       
   455       _nodes = &(_comp_nodes[comp]);
       
   456       int n = _nodes->size();
       
   457       if (n < 1 || (n == 1 && _out_arcs[(*_nodes)[0]].size() == 0)) {
       
   458         return false;
       
   459       }      
       
   460       for (int i = 0; i < n; ++i) {
       
   461         _data[(*_nodes)[i]].resize(n + 1);
       
   462       }
       
   463       return true;
       
   464     }
       
   465 
       
   466     // Process all rounds of computing path data for the current component.
       
   467     // _data[v][k] is the length of a shortest directed walk from the root
       
   468     // node to node v containing exactly k arcs.
       
   469     void processRounds() {
       
   470       Node start = (*_nodes)[0];
       
   471       _data[start][0] = PathData(true, 0);
       
   472       _process.clear();
       
   473       _process.push_back(start);
       
   474 
       
   475       int k, n = _nodes->size();
       
   476       for (k = 1; k <= n && int(_process.size()) < n; ++k) {
       
   477         processNextBuildRound(k);
       
   478       }
       
   479       for ( ; k <= n; ++k) {
       
   480         processNextFullRound(k);
       
   481       }
       
   482     }
       
   483 
       
   484     // Process one round and rebuild _process
       
   485     void processNextBuildRound(int k) {
       
   486       std::vector<Node> next;
       
   487       Node u, v;
       
   488       Arc e;
       
   489       LargeValue d;
       
   490       for (int i = 0; i < int(_process.size()); ++i) {
       
   491         u = _process[i];
       
   492         for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
       
   493           e = _out_arcs[u][j];
       
   494           v = _gr.target(e);
       
   495           d = _data[u][k-1].dist + _length[e];
       
   496           if (!_data[v][k].found) {
       
   497             next.push_back(v);
       
   498             _data[v][k] = PathData(true, _data[u][k-1].dist + _length[e], e);
       
   499           }
       
   500           else if (_tolerance.less(d, _data[v][k].dist)) {
       
   501             _data[v][k] = PathData(true, d, e);
       
   502           }
       
   503         }
       
   504       }
       
   505       _process.swap(next);
       
   506     }
       
   507 
       
   508     // Process one round using _nodes instead of _process
       
   509     void processNextFullRound(int k) {
       
   510       Node u, v;
       
   511       Arc e;
       
   512       LargeValue d;
       
   513       for (int i = 0; i < int(_nodes->size()); ++i) {
       
   514         u = (*_nodes)[i];
       
   515         for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
       
   516           e = _out_arcs[u][j];
       
   517           v = _gr.target(e);
       
   518           d = _data[u][k-1].dist + _length[e];
       
   519           if (!_data[v][k].found || _tolerance.less(d, _data[v][k].dist)) {
       
   520             _data[v][k] = PathData(true, d, e);
       
   521           }
       
   522         }
       
   523       }
       
   524     }
       
   525 
       
   526     // Update the minimum cycle mean
       
   527     void updateMinMean() {
       
   528       int n = _nodes->size();
       
   529       for (int i = 0; i < n; ++i) {
       
   530         Node u = (*_nodes)[i];
       
   531         if (!_data[u][n].found) continue;
       
   532         LargeValue length, max_length = 0;
       
   533         int size, max_size = 1;
       
   534         bool found_curr = false;
       
   535         for (int k = 0; k < n; ++k) {
       
   536           if (!_data[u][k].found) continue;
       
   537           length = _data[u][n].dist - _data[u][k].dist;
       
   538           size = n - k;
       
   539           if (!found_curr || length * max_size > max_length * size) {
       
   540             found_curr = true;
       
   541             max_length = length;
       
   542             max_size = size;
       
   543           }
       
   544         }
       
   545         if ( found_curr && (_cycle_node == INVALID ||
       
   546              max_length * _cycle_size < _cycle_length * max_size) ) {
       
   547           _cycle_length = max_length;
       
   548           _cycle_size = max_size;
       
   549           _cycle_node = u;
       
   550         }
       
   551       }
       
   552     }
       
   553 
       
   554   }; //class Karp
       
   555 
       
   556   ///@}
       
   557 
       
   558 } //namespace lemon
       
   559 
       
   560 #endif //LEMON_KARP_H