COIN-OR::LEMON - Graph Library

Ticket #179: 179-8-add-karp-f7ddbf131179.patch

File 179-8-add-karp-f7ddbf131179.patch, 19.8 KB (added by Peter Kovacs, 10 years ago)
  • lemon/Makefile.am

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1250016940 -7200
    # Node ID f7ddbf131179186bdcd76016762a1fc82c6b5a45
    # Parent  f2f32caf7d145c66cfbddb3d5e9999806cbd98fb
    Add Karp algorithm class (#179)
    based on the MinMeanCycle implementation in SVN -r3436.
    The interface is reworked to be the same as Howard's interface.
    
    diff --git a/lemon/Makefile.am b/lemon/Makefile.am
    a b  
    8585        lemon/grid_graph.h \
    8686        lemon/howard.h \
    8787        lemon/hypercube_graph.h \
     88        lemon/karp.h \
    8889        lemon/kruskal.h \
    8990        lemon/hao_orlin.h \
    9091        lemon/lgf_reader.h \
  • new file lemon/karp.h

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

    diff --git a/test/min_mean_cycle_test.cc b/test/min_mean_cycle_test.cc
    a b  
    2121
    2222#include <lemon/smart_graph.h>
    2323#include <lemon/lgf_reader.h>
    24 #include <lemon/howard.h>
    2524#include <lemon/path.h>
    2625#include <lemon/concepts/digraph.h>
    2726#include <lemon/concept_check.h>
    2827
     28#include <lemon/karp.h>
     29#include <lemon/howard.h>
     30
    2931#include "test_tools.h"
    3032
    3133using namespace lemon;
     
    141143  // Check the interface
    142144  {
    143145    typedef concepts::Digraph GR;
    144     typedef Howard<GR, concepts::ReadMap<GR::Arc, int> > IntMmcAlg;
    145     typedef Howard<GR, concepts::ReadMap<GR::Arc, float> > FloatMmcAlg;
     146
     147    // Karp
     148    checkConcept< MmcClassConcept<GR, int>,
     149                  Karp<GR, concepts::ReadMap<GR::Arc, int> > >();
     150    checkConcept< MmcClassConcept<GR, float>,
     151                  Karp<GR, concepts::ReadMap<GR::Arc, float> > >();
    146152   
    147     checkConcept<MmcClassConcept<GR, int>, IntMmcAlg>();
    148     checkConcept<MmcClassConcept<GR, float>, FloatMmcAlg>();
    149  
    150     if (IsSameType<IntMmcAlg::LargeValue, long_int>::result == 0)
    151       check(false, "Wrong LargeValue type");
    152     if (IsSameType<FloatMmcAlg::LargeValue, double>::result == 0)
    153       check(false, "Wrong LargeValue type");
     153    // Howard
     154    checkConcept< MmcClassConcept<GR, int>,
     155                  Howard<GR, concepts::ReadMap<GR::Arc, int> > >();
     156    checkConcept< MmcClassConcept<GR, float>,
     157                  Howard<GR, concepts::ReadMap<GR::Arc, float> > >();
     158
     159    if (IsSameType<Howard<GR, concepts::ReadMap<GR::Arc, int> >::LargeValue,
     160          long_int>::result == 0) check(false, "Wrong LargeValue type");
     161    if (IsSameType<Howard<GR, concepts::ReadMap<GR::Arc, float> >::LargeValue,
     162          double>::result == 0) check(false, "Wrong LargeValue type");
    154163  }
    155164
    156165  // Run various tests
     
    174183      arcMap("c4", c4).
    175184      run();
    176185
     186    // Karp
     187    checkMmcAlg<Karp<GR, IntArcMap> >(gr, l1, c1,  6, 3);
     188    checkMmcAlg<Karp<GR, IntArcMap> >(gr, l2, c2,  5, 2);
     189    checkMmcAlg<Karp<GR, IntArcMap> >(gr, l3, c3,  0, 1);
     190    checkMmcAlg<Karp<GR, IntArcMap> >(gr, l4, c4, -1, 1);
     191
     192    // Howard
    177193    checkMmcAlg<Howard<GR, IntArcMap> >(gr, l1, c1,  6, 3);
    178194    checkMmcAlg<Howard<GR, IntArcMap> >(gr, l2, c2,  5, 2);
    179195    checkMmcAlg<Howard<GR, IntArcMap> >(gr, l3, c3,  0, 1);