lemon/hartmann_orlin_mmc.h
author Alpar Juttner <alpar@cs.elte.hu>
Thu, 21 Jan 2021 08:36:53 +0100
changeset 1204 736a341e604b
parent 1080 c5cd8960df74
permissions -rw-r--r--
Update CMakeLists.txt to match current CMAKE versions (#641)
alpar@877
     1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
kpeter@766
     2
 *
alpar@877
     3
 * This file is a part of LEMON, a generic C++ optimization library.
kpeter@766
     4
 *
alpar@1092
     5
 * Copyright (C) 2003-2013
kpeter@766
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
kpeter@766
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
kpeter@766
     8
 *
kpeter@766
     9
 * Permission to use, modify and distribute this software is granted
kpeter@766
    10
 * provided that this copyright notice appears in all copies. For
kpeter@766
    11
 * precise terms see the accompanying LICENSE file.
kpeter@766
    12
 *
kpeter@766
    13
 * This software is provided "AS IS" with no warranty of any kind,
kpeter@766
    14
 * express or implied, and with no claim as to its suitability for any
kpeter@766
    15
 * purpose.
kpeter@766
    16
 *
kpeter@766
    17
 */
kpeter@766
    18
kpeter@864
    19
#ifndef LEMON_HARTMANN_ORLIN_MMC_H
kpeter@864
    20
#define LEMON_HARTMANN_ORLIN_MMC_H
kpeter@766
    21
kpeter@768
    22
/// \ingroup min_mean_cycle
kpeter@766
    23
///
kpeter@766
    24
/// \file
kpeter@766
    25
/// \brief Hartmann-Orlin's algorithm for finding a minimum mean cycle.
kpeter@766
    26
kpeter@766
    27
#include <vector>
kpeter@766
    28
#include <limits>
kpeter@766
    29
#include <lemon/core.h>
kpeter@766
    30
#include <lemon/path.h>
kpeter@766
    31
#include <lemon/tolerance.h>
kpeter@766
    32
#include <lemon/connectivity.h>
kpeter@766
    33
kpeter@766
    34
namespace lemon {
kpeter@766
    35
kpeter@864
    36
  /// \brief Default traits class of HartmannOrlinMmc class.
kpeter@766
    37
  ///
kpeter@864
    38
  /// Default traits class of HartmannOrlinMmc class.
kpeter@766
    39
  /// \tparam GR The type of the digraph.
kpeter@864
    40
  /// \tparam CM The type of the cost map.
kpeter@879
    41
  /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
kpeter@766
    42
#ifdef DOXYGEN
kpeter@864
    43
  template <typename GR, typename CM>
kpeter@766
    44
#else
kpeter@864
    45
  template <typename GR, typename CM,
kpeter@864
    46
    bool integer = std::numeric_limits<typename CM::Value>::is_integer>
kpeter@766
    47
#endif
kpeter@864
    48
  struct HartmannOrlinMmcDefaultTraits
kpeter@766
    49
  {
kpeter@766
    50
    /// The type of the digraph
kpeter@766
    51
    typedef GR Digraph;
kpeter@864
    52
    /// The type of the cost map
kpeter@864
    53
    typedef CM CostMap;
kpeter@864
    54
    /// The type of the arc costs
kpeter@864
    55
    typedef typename CostMap::Value Cost;
kpeter@766
    56
kpeter@864
    57
    /// \brief The large cost type used for internal computations
kpeter@766
    58
    ///
kpeter@864
    59
    /// The large cost type used for internal computations.
kpeter@864
    60
    /// It is \c long \c long if the \c Cost type is integer,
kpeter@766
    61
    /// otherwise it is \c double.
kpeter@864
    62
    /// \c Cost must be convertible to \c LargeCost.
kpeter@864
    63
    typedef double LargeCost;
kpeter@766
    64
kpeter@766
    65
    /// The tolerance type used for internal computations
kpeter@864
    66
    typedef lemon::Tolerance<LargeCost> Tolerance;
kpeter@766
    67
kpeter@766
    68
    /// \brief The path type of the found cycles
kpeter@766
    69
    ///
kpeter@766
    70
    /// The path type of the found cycles.
kpeter@766
    71
    /// It must conform to the \ref lemon::concepts::Path "Path" concept
kpeter@772
    72
    /// and it must have an \c addFront() function.
kpeter@766
    73
    typedef lemon::Path<Digraph> Path;
kpeter@766
    74
  };
kpeter@766
    75
kpeter@864
    76
  // Default traits class for integer cost types
kpeter@864
    77
  template <typename GR, typename CM>
kpeter@864
    78
  struct HartmannOrlinMmcDefaultTraits<GR, CM, true>
kpeter@766
    79
  {
kpeter@766
    80
    typedef GR Digraph;
kpeter@864
    81
    typedef CM CostMap;
kpeter@864
    82
    typedef typename CostMap::Value Cost;
kpeter@766
    83
#ifdef LEMON_HAVE_LONG_LONG
kpeter@864
    84
    typedef long long LargeCost;
kpeter@766
    85
#else
kpeter@864
    86
    typedef long LargeCost;
kpeter@766
    87
#endif
kpeter@864
    88
    typedef lemon::Tolerance<LargeCost> Tolerance;
kpeter@766
    89
    typedef lemon::Path<Digraph> Path;
kpeter@766
    90
  };
kpeter@766
    91
kpeter@766
    92
kpeter@768
    93
  /// \addtogroup min_mean_cycle
kpeter@766
    94
  /// @{
kpeter@766
    95
kpeter@766
    96
  /// \brief Implementation of the Hartmann-Orlin algorithm for finding
kpeter@766
    97
  /// a minimum mean cycle.
kpeter@766
    98
  ///
kpeter@766
    99
  /// This class implements the Hartmann-Orlin algorithm for finding
kpeter@864
   100
  /// a directed cycle of minimum mean cost in a digraph
alpar@1053
   101
  /// \cite hartmann93finding, \cite dasdan98minmeancycle.
kpeter@1049
   102
  /// This method is based on \ref KarpMmc "Karp"'s original algorithm, but
kpeter@1049
   103
  /// applies an early termination scheme. It makes the algorithm
kpeter@1049
   104
  /// significantly faster for some problem instances, but slower for others.
kpeter@1080
   105
  /// The algorithm runs in time O(nm) and uses space O(n<sup>2</sup>+m).
kpeter@766
   106
  ///
kpeter@766
   107
  /// \tparam GR The type of the digraph the algorithm runs on.
kpeter@864
   108
  /// \tparam CM The type of the cost map. The default
kpeter@766
   109
  /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
kpeter@825
   110
  /// \tparam TR The traits class that defines various types used by the
kpeter@864
   111
  /// algorithm. By default, it is \ref HartmannOrlinMmcDefaultTraits
kpeter@864
   112
  /// "HartmannOrlinMmcDefaultTraits<GR, CM>".
kpeter@825
   113
  /// In most cases, this parameter should not be set directly,
kpeter@825
   114
  /// consider to use the named template parameters instead.
kpeter@766
   115
#ifdef DOXYGEN
kpeter@864
   116
  template <typename GR, typename CM, typename TR>
kpeter@766
   117
#else
kpeter@766
   118
  template < typename GR,
kpeter@864
   119
             typename CM = typename GR::template ArcMap<int>,
kpeter@864
   120
             typename TR = HartmannOrlinMmcDefaultTraits<GR, CM> >
kpeter@766
   121
#endif
kpeter@864
   122
  class HartmannOrlinMmc
kpeter@766
   123
  {
kpeter@766
   124
  public:
kpeter@766
   125
kpeter@766
   126
    /// The type of the digraph
kpeter@766
   127
    typedef typename TR::Digraph Digraph;
kpeter@864
   128
    /// The type of the cost map
kpeter@864
   129
    typedef typename TR::CostMap CostMap;
kpeter@864
   130
    /// The type of the arc costs
kpeter@864
   131
    typedef typename TR::Cost Cost;
kpeter@766
   132
kpeter@864
   133
    /// \brief The large cost type
kpeter@766
   134
    ///
kpeter@864
   135
    /// The large cost type used for internal computations.
kpeter@864
   136
    /// By default, it is \c long \c long if the \c Cost type is integer,
kpeter@766
   137
    /// otherwise it is \c double.
kpeter@864
   138
    typedef typename TR::LargeCost LargeCost;
kpeter@766
   139
kpeter@766
   140
    /// The tolerance type
kpeter@766
   141
    typedef typename TR::Tolerance Tolerance;
kpeter@766
   142
kpeter@766
   143
    /// \brief The path type of the found cycles
kpeter@766
   144
    ///
kpeter@766
   145
    /// The path type of the found cycles.
alpar@1074
   146
    /// Using the \ref lemon::HartmannOrlinMmcDefaultTraits
alpar@1074
   147
    /// "default traits class",
kpeter@766
   148
    /// it is \ref lemon::Path "Path<Digraph>".
kpeter@766
   149
    typedef typename TR::Path Path;
kpeter@766
   150
alpar@1074
   151
    /// \brief The
alpar@1074
   152
    /// \ref lemon::HartmannOrlinMmcDefaultTraits "traits class"
alpar@1074
   153
    /// of the algorithm
kpeter@766
   154
    typedef TR Traits;
kpeter@766
   155
kpeter@766
   156
  private:
kpeter@766
   157
kpeter@766
   158
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
kpeter@766
   159
kpeter@766
   160
    // Data sturcture for path data
kpeter@766
   161
    struct PathData
kpeter@766
   162
    {
kpeter@864
   163
      LargeCost dist;
kpeter@766
   164
      Arc pred;
kpeter@864
   165
      PathData(LargeCost d, Arc p = INVALID) :
kpeter@767
   166
        dist(d), pred(p) {}
kpeter@766
   167
    };
kpeter@766
   168
kpeter@766
   169
    typedef typename Digraph::template NodeMap<std::vector<PathData> >
kpeter@766
   170
      PathDataNodeMap;
kpeter@766
   171
kpeter@766
   172
  private:
kpeter@766
   173
kpeter@766
   174
    // The digraph the algorithm runs on
kpeter@766
   175
    const Digraph &_gr;
kpeter@864
   176
    // The cost of the arcs
kpeter@864
   177
    const CostMap &_cost;
kpeter@766
   178
kpeter@766
   179
    // Data for storing the strongly connected components
kpeter@766
   180
    int _comp_num;
kpeter@766
   181
    typename Digraph::template NodeMap<int> _comp;
kpeter@766
   182
    std::vector<std::vector<Node> > _comp_nodes;
kpeter@766
   183
    std::vector<Node>* _nodes;
kpeter@766
   184
    typename Digraph::template NodeMap<std::vector<Arc> > _out_arcs;
kpeter@766
   185
kpeter@766
   186
    // Data for the found cycles
kpeter@766
   187
    bool _curr_found, _best_found;
kpeter@864
   188
    LargeCost _curr_cost, _best_cost;
kpeter@766
   189
    int _curr_size, _best_size;
kpeter@766
   190
    Node _curr_node, _best_node;
kpeter@766
   191
    int _curr_level, _best_level;
kpeter@766
   192
kpeter@766
   193
    Path *_cycle_path;
kpeter@766
   194
    bool _local_path;
kpeter@766
   195
kpeter@766
   196
    // Node map for storing path data
kpeter@766
   197
    PathDataNodeMap _data;
kpeter@766
   198
    // The processed nodes in the last round
kpeter@766
   199
    std::vector<Node> _process;
kpeter@766
   200
kpeter@766
   201
    Tolerance _tolerance;
kpeter@766
   202
kpeter@767
   203
    // Infinite constant
kpeter@864
   204
    const LargeCost INF;
kpeter@767
   205
kpeter@766
   206
  public:
kpeter@766
   207
kpeter@766
   208
    /// \name Named Template Parameters
kpeter@766
   209
    /// @{
kpeter@766
   210
kpeter@766
   211
    template <typename T>
kpeter@864
   212
    struct SetLargeCostTraits : public Traits {
kpeter@864
   213
      typedef T LargeCost;
kpeter@766
   214
      typedef lemon::Tolerance<T> Tolerance;
kpeter@766
   215
    };
kpeter@766
   216
kpeter@766
   217
    /// \brief \ref named-templ-param "Named parameter" for setting
kpeter@864
   218
    /// \c LargeCost type.
kpeter@766
   219
    ///
kpeter@864
   220
    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
kpeter@766
   221
    /// type. It is used for internal computations in the algorithm.
kpeter@766
   222
    template <typename T>
kpeter@864
   223
    struct SetLargeCost
kpeter@864
   224
      : public HartmannOrlinMmc<GR, CM, SetLargeCostTraits<T> > {
kpeter@864
   225
      typedef HartmannOrlinMmc<GR, CM, SetLargeCostTraits<T> > Create;
kpeter@766
   226
    };
kpeter@766
   227
kpeter@766
   228
    template <typename T>
kpeter@766
   229
    struct SetPathTraits : public Traits {
kpeter@766
   230
      typedef T Path;
kpeter@766
   231
    };
kpeter@766
   232
kpeter@766
   233
    /// \brief \ref named-templ-param "Named parameter" for setting
kpeter@766
   234
    /// \c %Path type.
kpeter@766
   235
    ///
kpeter@766
   236
    /// \ref named-templ-param "Named parameter" for setting the \c %Path
kpeter@766
   237
    /// type of the found cycles.
kpeter@766
   238
    /// It must conform to the \ref lemon::concepts::Path "Path" concept
kpeter@766
   239
    /// and it must have an \c addFront() function.
kpeter@766
   240
    template <typename T>
kpeter@766
   241
    struct SetPath
kpeter@864
   242
      : public HartmannOrlinMmc<GR, CM, SetPathTraits<T> > {
kpeter@864
   243
      typedef HartmannOrlinMmc<GR, CM, SetPathTraits<T> > Create;
kpeter@766
   244
    };
kpeter@766
   245
kpeter@766
   246
    /// @}
kpeter@766
   247
kpeter@863
   248
  protected:
kpeter@863
   249
kpeter@864
   250
    HartmannOrlinMmc() {}
kpeter@863
   251
kpeter@766
   252
  public:
kpeter@766
   253
kpeter@766
   254
    /// \brief Constructor.
kpeter@766
   255
    ///
kpeter@766
   256
    /// The constructor of the class.
kpeter@766
   257
    ///
kpeter@766
   258
    /// \param digraph The digraph the algorithm runs on.
kpeter@864
   259
    /// \param cost The costs of the arcs.
kpeter@864
   260
    HartmannOrlinMmc( const Digraph &digraph,
kpeter@864
   261
                      const CostMap &cost ) :
kpeter@864
   262
      _gr(digraph), _cost(cost), _comp(digraph), _out_arcs(digraph),
kpeter@864
   263
      _best_found(false), _best_cost(0), _best_size(1),
kpeter@767
   264
      _cycle_path(NULL), _local_path(false), _data(digraph),
kpeter@864
   265
      INF(std::numeric_limits<LargeCost>::has_infinity ?
kpeter@864
   266
          std::numeric_limits<LargeCost>::infinity() :
kpeter@864
   267
          std::numeric_limits<LargeCost>::max())
kpeter@766
   268
    {}
kpeter@766
   269
kpeter@766
   270
    /// Destructor.
kpeter@864
   271
    ~HartmannOrlinMmc() {
kpeter@766
   272
      if (_local_path) delete _cycle_path;
kpeter@766
   273
    }
kpeter@766
   274
kpeter@766
   275
    /// \brief Set the path structure for storing the found cycle.
kpeter@766
   276
    ///
kpeter@766
   277
    /// This function sets an external path structure for storing the
kpeter@766
   278
    /// found cycle.
kpeter@766
   279
    ///
kpeter@766
   280
    /// If you don't call this function before calling \ref run() or
kpeter@1049
   281
    /// \ref findCycleMean(), a local \ref Path "path" structure
kpeter@1049
   282
    /// will be allocated. The destuctor deallocates this automatically
kpeter@766
   283
    /// allocated object, of course.
kpeter@766
   284
    ///
kpeter@766
   285
    /// \note The algorithm calls only the \ref lemon::Path::addFront()
kpeter@766
   286
    /// "addFront()" function of the given path structure.
kpeter@766
   287
    ///
kpeter@766
   288
    /// \return <tt>(*this)</tt>
kpeter@864
   289
    HartmannOrlinMmc& cycle(Path &path) {
kpeter@766
   290
      if (_local_path) {
kpeter@766
   291
        delete _cycle_path;
kpeter@766
   292
        _local_path = false;
kpeter@766
   293
      }
kpeter@766
   294
      _cycle_path = &path;
kpeter@766
   295
      return *this;
kpeter@766
   296
    }
kpeter@766
   297
kpeter@769
   298
    /// \brief Set the tolerance used by the algorithm.
kpeter@769
   299
    ///
kpeter@769
   300
    /// This function sets the tolerance object used by the algorithm.
kpeter@769
   301
    ///
kpeter@769
   302
    /// \return <tt>(*this)</tt>
kpeter@864
   303
    HartmannOrlinMmc& tolerance(const Tolerance& tolerance) {
kpeter@769
   304
      _tolerance = tolerance;
kpeter@769
   305
      return *this;
kpeter@769
   306
    }
kpeter@769
   307
kpeter@769
   308
    /// \brief Return a const reference to the tolerance.
kpeter@769
   309
    ///
kpeter@769
   310
    /// This function returns a const reference to the tolerance object
kpeter@769
   311
    /// used by the algorithm.
kpeter@769
   312
    const Tolerance& tolerance() const {
kpeter@769
   313
      return _tolerance;
kpeter@769
   314
    }
kpeter@769
   315
kpeter@766
   316
    /// \name Execution control
kpeter@766
   317
    /// The simplest way to execute the algorithm is to call the \ref run()
kpeter@766
   318
    /// function.\n
kpeter@864
   319
    /// If you only need the minimum mean cost, you may call
kpeter@864
   320
    /// \ref findCycleMean().
kpeter@766
   321
kpeter@766
   322
    /// @{
kpeter@766
   323
kpeter@766
   324
    /// \brief Run the algorithm.
kpeter@766
   325
    ///
kpeter@766
   326
    /// This function runs the algorithm.
kpeter@766
   327
    /// It can be called more than once (e.g. if the underlying digraph
kpeter@864
   328
    /// and/or the arc costs have been modified).
kpeter@766
   329
    ///
kpeter@766
   330
    /// \return \c true if a directed cycle exists in the digraph.
kpeter@766
   331
    ///
kpeter@766
   332
    /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
kpeter@766
   333
    /// \code
kpeter@864
   334
    ///   return mmc.findCycleMean() && mmc.findCycle();
kpeter@766
   335
    /// \endcode
kpeter@766
   336
    bool run() {
kpeter@864
   337
      return findCycleMean() && findCycle();
kpeter@766
   338
    }
kpeter@766
   339
kpeter@766
   340
    /// \brief Find the minimum cycle mean.
kpeter@766
   341
    ///
kpeter@864
   342
    /// This function finds the minimum mean cost of the directed
kpeter@766
   343
    /// cycles in the digraph.
kpeter@766
   344
    ///
kpeter@766
   345
    /// \return \c true if a directed cycle exists in the digraph.
kpeter@864
   346
    bool findCycleMean() {
kpeter@766
   347
      // Initialization and find strongly connected components
kpeter@766
   348
      init();
kpeter@766
   349
      findComponents();
alpar@877
   350
kpeter@766
   351
      // Find the minimum cycle mean in the components
kpeter@766
   352
      for (int comp = 0; comp < _comp_num; ++comp) {
kpeter@766
   353
        if (!initComponent(comp)) continue;
kpeter@766
   354
        processRounds();
alpar@877
   355
kpeter@766
   356
        // Update the best cycle (global minimum mean cycle)
alpar@877
   357
        if ( _curr_found && (!_best_found ||
kpeter@864
   358
             _curr_cost * _best_size < _best_cost * _curr_size) ) {
kpeter@766
   359
          _best_found = true;
kpeter@864
   360
          _best_cost = _curr_cost;
kpeter@766
   361
          _best_size = _curr_size;
kpeter@766
   362
          _best_node = _curr_node;
kpeter@766
   363
          _best_level = _curr_level;
kpeter@766
   364
        }
kpeter@766
   365
      }
kpeter@766
   366
      return _best_found;
kpeter@766
   367
    }
kpeter@766
   368
kpeter@766
   369
    /// \brief Find a minimum mean directed cycle.
kpeter@766
   370
    ///
kpeter@864
   371
    /// This function finds a directed cycle of minimum mean cost
kpeter@864
   372
    /// in the digraph using the data computed by findCycleMean().
kpeter@766
   373
    ///
kpeter@766
   374
    /// \return \c true if a directed cycle exists in the digraph.
kpeter@766
   375
    ///
kpeter@864
   376
    /// \pre \ref findCycleMean() must be called before using this function.
kpeter@766
   377
    bool findCycle() {
kpeter@766
   378
      if (!_best_found) return false;
kpeter@766
   379
      IntNodeMap reached(_gr, -1);
kpeter@766
   380
      int r = _best_level + 1;
kpeter@766
   381
      Node u = _best_node;
kpeter@766
   382
      while (reached[u] < 0) {
kpeter@766
   383
        reached[u] = --r;
kpeter@766
   384
        u = _gr.source(_data[u][r].pred);
kpeter@766
   385
      }
kpeter@766
   386
      r = reached[u];
kpeter@766
   387
      Arc e = _data[u][r].pred;
kpeter@766
   388
      _cycle_path->addFront(e);
kpeter@864
   389
      _best_cost = _cost[e];
kpeter@766
   390
      _best_size = 1;
kpeter@766
   391
      Node v;
kpeter@766
   392
      while ((v = _gr.source(e)) != u) {
kpeter@766
   393
        e = _data[v][--r].pred;
kpeter@766
   394
        _cycle_path->addFront(e);
kpeter@864
   395
        _best_cost += _cost[e];
kpeter@766
   396
        ++_best_size;
kpeter@766
   397
      }
kpeter@766
   398
      return true;
kpeter@766
   399
    }
kpeter@766
   400
kpeter@766
   401
    /// @}
kpeter@766
   402
kpeter@766
   403
    /// \name Query Functions
kpeter@766
   404
    /// The results of the algorithm can be obtained using these
kpeter@766
   405
    /// functions.\n
kpeter@766
   406
    /// The algorithm should be executed before using them.
kpeter@766
   407
kpeter@766
   408
    /// @{
kpeter@766
   409
kpeter@864
   410
    /// \brief Return the total cost of the found cycle.
kpeter@766
   411
    ///
kpeter@864
   412
    /// This function returns the total cost of the found cycle.
kpeter@766
   413
    ///
kpeter@864
   414
    /// \pre \ref run() or \ref findCycleMean() must be called before
kpeter@766
   415
    /// using this function.
kpeter@864
   416
    Cost cycleCost() const {
kpeter@864
   417
      return static_cast<Cost>(_best_cost);
kpeter@766
   418
    }
kpeter@766
   419
kpeter@766
   420
    /// \brief Return the number of arcs on the found cycle.
kpeter@766
   421
    ///
kpeter@766
   422
    /// This function returns the number of arcs on the found cycle.
kpeter@766
   423
    ///
kpeter@864
   424
    /// \pre \ref run() or \ref findCycleMean() must be called before
kpeter@766
   425
    /// using this function.
kpeter@864
   426
    int cycleSize() const {
kpeter@766
   427
      return _best_size;
kpeter@766
   428
    }
kpeter@766
   429
kpeter@864
   430
    /// \brief Return the mean cost of the found cycle.
kpeter@766
   431
    ///
kpeter@864
   432
    /// This function returns the mean cost of the found cycle.
kpeter@766
   433
    ///
kpeter@766
   434
    /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
kpeter@766
   435
    /// following code.
kpeter@766
   436
    /// \code
kpeter@864
   437
    ///   return static_cast<double>(alg.cycleCost()) / alg.cycleSize();
kpeter@766
   438
    /// \endcode
kpeter@766
   439
    ///
kpeter@864
   440
    /// \pre \ref run() or \ref findCycleMean() must be called before
kpeter@766
   441
    /// using this function.
kpeter@766
   442
    double cycleMean() const {
kpeter@864
   443
      return static_cast<double>(_best_cost) / _best_size;
kpeter@766
   444
    }
kpeter@766
   445
kpeter@766
   446
    /// \brief Return the found cycle.
kpeter@766
   447
    ///
kpeter@766
   448
    /// This function returns a const reference to the path structure
kpeter@766
   449
    /// storing the found cycle.
kpeter@766
   450
    ///
kpeter@766
   451
    /// \pre \ref run() or \ref findCycle() must be called before using
kpeter@766
   452
    /// this function.
kpeter@766
   453
    const Path& cycle() const {
kpeter@766
   454
      return *_cycle_path;
kpeter@766
   455
    }
kpeter@766
   456
kpeter@766
   457
    ///@}
kpeter@766
   458
kpeter@766
   459
  private:
kpeter@766
   460
kpeter@766
   461
    // Initialization
kpeter@766
   462
    void init() {
kpeter@766
   463
      if (!_cycle_path) {
kpeter@766
   464
        _local_path = true;
kpeter@766
   465
        _cycle_path = new Path;
kpeter@766
   466
      }
kpeter@766
   467
      _cycle_path->clear();
kpeter@766
   468
      _best_found = false;
kpeter@864
   469
      _best_cost = 0;
kpeter@766
   470
      _best_size = 1;
kpeter@766
   471
      _cycle_path->clear();
kpeter@766
   472
      for (NodeIt u(_gr); u != INVALID; ++u)
kpeter@766
   473
        _data[u].clear();
kpeter@766
   474
    }
kpeter@766
   475
kpeter@766
   476
    // Find strongly connected components and initialize _comp_nodes
kpeter@766
   477
    // and _out_arcs
kpeter@766
   478
    void findComponents() {
kpeter@766
   479
      _comp_num = stronglyConnectedComponents(_gr, _comp);
kpeter@766
   480
      _comp_nodes.resize(_comp_num);
kpeter@766
   481
      if (_comp_num == 1) {
kpeter@766
   482
        _comp_nodes[0].clear();
kpeter@766
   483
        for (NodeIt n(_gr); n != INVALID; ++n) {
kpeter@766
   484
          _comp_nodes[0].push_back(n);
kpeter@766
   485
          _out_arcs[n].clear();
kpeter@766
   486
          for (OutArcIt a(_gr, n); a != INVALID; ++a) {
kpeter@766
   487
            _out_arcs[n].push_back(a);
kpeter@766
   488
          }
kpeter@766
   489
        }
kpeter@766
   490
      } else {
kpeter@766
   491
        for (int i = 0; i < _comp_num; ++i)
kpeter@766
   492
          _comp_nodes[i].clear();
kpeter@766
   493
        for (NodeIt n(_gr); n != INVALID; ++n) {
kpeter@766
   494
          int k = _comp[n];
kpeter@766
   495
          _comp_nodes[k].push_back(n);
kpeter@766
   496
          _out_arcs[n].clear();
kpeter@766
   497
          for (OutArcIt a(_gr, n); a != INVALID; ++a) {
kpeter@766
   498
            if (_comp[_gr.target(a)] == k) _out_arcs[n].push_back(a);
kpeter@766
   499
          }
kpeter@766
   500
        }
kpeter@766
   501
      }
kpeter@766
   502
    }
kpeter@766
   503
kpeter@766
   504
    // Initialize path data for the current component
kpeter@766
   505
    bool initComponent(int comp) {
kpeter@766
   506
      _nodes = &(_comp_nodes[comp]);
kpeter@766
   507
      int n = _nodes->size();
kpeter@766
   508
      if (n < 1 || (n == 1 && _out_arcs[(*_nodes)[0]].size() == 0)) {
kpeter@766
   509
        return false;
alpar@877
   510
      }
kpeter@766
   511
      for (int i = 0; i < n; ++i) {
kpeter@767
   512
        _data[(*_nodes)[i]].resize(n + 1, PathData(INF));
kpeter@766
   513
      }
kpeter@766
   514
      return true;
kpeter@766
   515
    }
kpeter@766
   516
kpeter@766
   517
    // Process all rounds of computing path data for the current component.
kpeter@864
   518
    // _data[v][k] is the cost of a shortest directed walk from the root
kpeter@766
   519
    // node to node v containing exactly k arcs.
kpeter@766
   520
    void processRounds() {
kpeter@766
   521
      Node start = (*_nodes)[0];
kpeter@767
   522
      _data[start][0] = PathData(0);
kpeter@766
   523
      _process.clear();
kpeter@766
   524
      _process.push_back(start);
kpeter@766
   525
kpeter@766
   526
      int k, n = _nodes->size();
kpeter@766
   527
      int next_check = 4;
kpeter@766
   528
      bool terminate = false;
kpeter@766
   529
      for (k = 1; k <= n && int(_process.size()) < n && !terminate; ++k) {
kpeter@766
   530
        processNextBuildRound(k);
kpeter@766
   531
        if (k == next_check || k == n) {
kpeter@766
   532
          terminate = checkTermination(k);
kpeter@766
   533
          next_check = next_check * 3 / 2;
kpeter@766
   534
        }
kpeter@766
   535
      }
kpeter@766
   536
      for ( ; k <= n && !terminate; ++k) {
kpeter@766
   537
        processNextFullRound(k);
kpeter@766
   538
        if (k == next_check || k == n) {
kpeter@766
   539
          terminate = checkTermination(k);
kpeter@766
   540
          next_check = next_check * 3 / 2;
kpeter@766
   541
        }
kpeter@766
   542
      }
kpeter@766
   543
    }
kpeter@766
   544
kpeter@766
   545
    // Process one round and rebuild _process
kpeter@766
   546
    void processNextBuildRound(int k) {
kpeter@766
   547
      std::vector<Node> next;
kpeter@766
   548
      Node u, v;
kpeter@766
   549
      Arc e;
kpeter@864
   550
      LargeCost d;
kpeter@766
   551
      for (int i = 0; i < int(_process.size()); ++i) {
kpeter@766
   552
        u = _process[i];
kpeter@766
   553
        for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
kpeter@766
   554
          e = _out_arcs[u][j];
kpeter@766
   555
          v = _gr.target(e);
kpeter@864
   556
          d = _data[u][k-1].dist + _cost[e];
kpeter@767
   557
          if (_tolerance.less(d, _data[v][k].dist)) {
kpeter@767
   558
            if (_data[v][k].dist == INF) next.push_back(v);
kpeter@767
   559
            _data[v][k] = PathData(d, e);
kpeter@766
   560
          }
kpeter@766
   561
        }
kpeter@766
   562
      }
kpeter@766
   563
      _process.swap(next);
kpeter@766
   564
    }
kpeter@766
   565
kpeter@766
   566
    // Process one round using _nodes instead of _process
kpeter@766
   567
    void processNextFullRound(int k) {
kpeter@766
   568
      Node u, v;
kpeter@766
   569
      Arc e;
kpeter@864
   570
      LargeCost d;
kpeter@766
   571
      for (int i = 0; i < int(_nodes->size()); ++i) {
kpeter@766
   572
        u = (*_nodes)[i];
kpeter@766
   573
        for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
kpeter@766
   574
          e = _out_arcs[u][j];
kpeter@766
   575
          v = _gr.target(e);
kpeter@864
   576
          d = _data[u][k-1].dist + _cost[e];
kpeter@767
   577
          if (_tolerance.less(d, _data[v][k].dist)) {
kpeter@767
   578
            _data[v][k] = PathData(d, e);
kpeter@766
   579
          }
kpeter@766
   580
        }
kpeter@766
   581
      }
kpeter@766
   582
    }
alpar@877
   583
kpeter@766
   584
    // Check early termination
kpeter@766
   585
    bool checkTermination(int k) {
kpeter@766
   586
      typedef std::pair<int, int> Pair;
kpeter@766
   587
      typename GR::template NodeMap<Pair> level(_gr, Pair(-1, 0));
kpeter@864
   588
      typename GR::template NodeMap<LargeCost> pi(_gr);
kpeter@766
   589
      int n = _nodes->size();
kpeter@864
   590
      LargeCost cost;
kpeter@766
   591
      int size;
kpeter@766
   592
      Node u;
alpar@877
   593
kpeter@766
   594
      // Search for cycles that are already found
kpeter@766
   595
      _curr_found = false;
kpeter@766
   596
      for (int i = 0; i < n; ++i) {
kpeter@766
   597
        u = (*_nodes)[i];
kpeter@767
   598
        if (_data[u][k].dist == INF) continue;
kpeter@766
   599
        for (int j = k; j >= 0; --j) {
kpeter@766
   600
          if (level[u].first == i && level[u].second > 0) {
kpeter@766
   601
            // A cycle is found
kpeter@864
   602
            cost = _data[u][level[u].second].dist - _data[u][j].dist;
kpeter@766
   603
            size = level[u].second - j;
kpeter@864
   604
            if (!_curr_found || cost * _curr_size < _curr_cost * size) {
kpeter@864
   605
              _curr_cost = cost;
kpeter@766
   606
              _curr_size = size;
kpeter@766
   607
              _curr_node = u;
kpeter@766
   608
              _curr_level = level[u].second;
kpeter@766
   609
              _curr_found = true;
kpeter@766
   610
            }
kpeter@766
   611
          }
kpeter@766
   612
          level[u] = Pair(i, j);
deba@795
   613
          if (j != 0) {
alpar@877
   614
            u = _gr.source(_data[u][j].pred);
alpar@877
   615
          }
kpeter@766
   616
        }
kpeter@766
   617
      }
kpeter@766
   618
kpeter@766
   619
      // If at least one cycle is found, check the optimality condition
kpeter@864
   620
      LargeCost d;
kpeter@766
   621
      if (_curr_found && k < n) {
kpeter@766
   622
        // Find node potentials
kpeter@766
   623
        for (int i = 0; i < n; ++i) {
kpeter@766
   624
          u = (*_nodes)[i];
kpeter@767
   625
          pi[u] = INF;
kpeter@766
   626
          for (int j = 0; j <= k; ++j) {
kpeter@767
   627
            if (_data[u][j].dist < INF) {
kpeter@864
   628
              d = _data[u][j].dist * _curr_size - j * _curr_cost;
kpeter@767
   629
              if (_tolerance.less(d, pi[u])) pi[u] = d;
kpeter@766
   630
            }
kpeter@766
   631
          }
kpeter@766
   632
        }
kpeter@766
   633
kpeter@766
   634
        // Check the optimality condition for all arcs
kpeter@766
   635
        bool done = true;
kpeter@766
   636
        for (ArcIt a(_gr); a != INVALID; ++a) {
kpeter@864
   637
          if (_tolerance.less(_cost[a] * _curr_size - _curr_cost,
kpeter@766
   638
                              pi[_gr.target(a)] - pi[_gr.source(a)]) ) {
kpeter@766
   639
            done = false;
kpeter@766
   640
            break;
kpeter@766
   641
          }
kpeter@766
   642
        }
kpeter@766
   643
        return done;
kpeter@766
   644
      }
kpeter@766
   645
      return (k == n);
kpeter@766
   646
    }
kpeter@766
   647
kpeter@864
   648
  }; //class HartmannOrlinMmc
kpeter@766
   649
kpeter@766
   650
  ///@}
kpeter@766
   651
kpeter@766
   652
} //namespace lemon
kpeter@766
   653
kpeter@864
   654
#endif //LEMON_HARTMANN_ORLIN_MMC_H