lemon/grosso_locatelli_pullan_mc.h
author Alpar Juttner <alpar@cs.elte.hu>
Sat, 25 May 2013 06:59:31 +0200
changeset 1064 fc3854d936f7
parent 904 c279b19abc62
child 1053 1c978b5bcc65
permissions -rw-r--r--
Enable/disable options for LP/MIP backends (#465)
kpeter@904
     1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
kpeter@904
     2
 *
kpeter@904
     3
 * This file is a part of LEMON, a generic C++ optimization library.
kpeter@904
     4
 *
kpeter@904
     5
 * Copyright (C) 2003-2010
kpeter@904
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
kpeter@904
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
kpeter@904
     8
 *
kpeter@904
     9
 * Permission to use, modify and distribute this software is granted
kpeter@904
    10
 * provided that this copyright notice appears in all copies. For
kpeter@904
    11
 * precise terms see the accompanying LICENSE file.
kpeter@904
    12
 *
kpeter@904
    13
 * This software is provided "AS IS" with no warranty of any kind,
kpeter@904
    14
 * express or implied, and with no claim as to its suitability for any
kpeter@904
    15
 * purpose.
kpeter@904
    16
 *
kpeter@904
    17
 */
kpeter@904
    18
kpeter@904
    19
#ifndef LEMON_GROSSO_LOCATELLI_PULLAN_MC_H
kpeter@904
    20
#define LEMON_GROSSO_LOCATELLI_PULLAN_MC_H
kpeter@904
    21
kpeter@904
    22
/// \ingroup approx_algs
kpeter@904
    23
///
kpeter@904
    24
/// \file
kpeter@904
    25
/// \brief The iterated local search algorithm of Grosso, Locatelli, and Pullan
kpeter@904
    26
/// for the maximum clique problem
kpeter@904
    27
kpeter@904
    28
#include <vector>
kpeter@904
    29
#include <limits>
kpeter@904
    30
#include <lemon/core.h>
kpeter@904
    31
#include <lemon/random.h>
kpeter@904
    32
kpeter@904
    33
namespace lemon {
kpeter@904
    34
kpeter@904
    35
  /// \addtogroup approx_algs
kpeter@904
    36
  /// @{
kpeter@904
    37
kpeter@904
    38
  /// \brief Implementation of the iterated local search algorithm of Grosso,
kpeter@904
    39
  /// Locatelli, and Pullan for the maximum clique problem
kpeter@904
    40
  ///
kpeter@904
    41
  /// \ref GrossoLocatelliPullanMc implements the iterated local search
kpeter@904
    42
  /// algorithm of Grosso, Locatelli, and Pullan for solving the \e maximum
kpeter@904
    43
  /// \e clique \e problem \ref grosso08maxclique.
kpeter@904
    44
  /// It is to find the largest complete subgraph (\e clique) in an
kpeter@904
    45
  /// undirected graph, i.e., the largest set of nodes where each
kpeter@904
    46
  /// pair of nodes is connected.
kpeter@904
    47
  ///
kpeter@904
    48
  /// This class provides a simple but highly efficient and robust heuristic
kpeter@918
    49
  /// method that quickly finds a quite large clique, but not necessarily the
kpeter@904
    50
  /// largest one.
kpeter@918
    51
  /// The algorithm performs a certain number of iterations to find several
kpeter@918
    52
  /// cliques and selects the largest one among them. Various limits can be
kpeter@918
    53
  /// specified to control the running time and the effectiveness of the
kpeter@918
    54
  /// search process.
kpeter@904
    55
  ///
kpeter@904
    56
  /// \tparam GR The undirected graph type the algorithm runs on.
kpeter@904
    57
  ///
kpeter@904
    58
  /// \note %GrossoLocatelliPullanMc provides three different node selection
kpeter@904
    59
  /// rules, from which the most powerful one is used by default.
kpeter@904
    60
  /// For more information, see \ref SelectionRule.
kpeter@904
    61
  template <typename GR>
kpeter@904
    62
  class GrossoLocatelliPullanMc
kpeter@904
    63
  {
kpeter@904
    64
  public:
kpeter@904
    65
kpeter@904
    66
    /// \brief Constants for specifying the node selection rule.
kpeter@904
    67
    ///
kpeter@904
    68
    /// Enum type containing constants for specifying the node selection rule
kpeter@904
    69
    /// for the \ref run() function.
kpeter@904
    70
    ///
kpeter@904
    71
    /// During the algorithm, nodes are selected for addition to the current
kpeter@904
    72
    /// clique according to the applied rule.
kpeter@904
    73
    /// In general, the PENALTY_BASED rule turned out to be the most powerful
kpeter@904
    74
    /// and the most robust, thus it is the default option.
kpeter@904
    75
    /// However, another selection rule can be specified using the \ref run()
kpeter@904
    76
    /// function with the proper parameter.
kpeter@904
    77
    enum SelectionRule {
kpeter@904
    78
kpeter@904
    79
      /// A node is selected randomly without any evaluation at each step.
kpeter@904
    80
      RANDOM,
kpeter@904
    81
kpeter@904
    82
      /// A node of maximum degree is selected randomly at each step.
kpeter@904
    83
      DEGREE_BASED,
kpeter@904
    84
kpeter@904
    85
      /// A node of minimum penalty is selected randomly at each step.
kpeter@904
    86
      /// The node penalties are updated adaptively after each stage of the
kpeter@904
    87
      /// search process.
kpeter@904
    88
      PENALTY_BASED
kpeter@904
    89
    };
kpeter@904
    90
kpeter@918
    91
    /// \brief Constants for the causes of search termination.
kpeter@918
    92
    ///
kpeter@918
    93
    /// Enum type containing constants for the different causes of search
kpeter@918
    94
    /// termination. The \ref run() function returns one of these values.
kpeter@918
    95
    enum TerminationCause {
kpeter@918
    96
kpeter@918
    97
      /// The iteration count limit is reached.
kpeter@918
    98
      ITERATION_LIMIT,
kpeter@918
    99
kpeter@918
   100
      /// The step count limit is reached.
kpeter@918
   101
      STEP_LIMIT,
kpeter@918
   102
kpeter@918
   103
      /// The clique size limit is reached.
kpeter@918
   104
      SIZE_LIMIT
kpeter@918
   105
    };
kpeter@918
   106
kpeter@904
   107
  private:
kpeter@904
   108
kpeter@904
   109
    TEMPLATE_GRAPH_TYPEDEFS(GR);
kpeter@904
   110
kpeter@904
   111
    typedef std::vector<int> IntVector;
kpeter@904
   112
    typedef std::vector<char> BoolVector;
kpeter@904
   113
    typedef std::vector<BoolVector> BoolMatrix;
kpeter@904
   114
    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
kpeter@904
   115
kpeter@918
   116
    // The underlying graph
kpeter@904
   117
    const GR &_graph;
kpeter@904
   118
    IntNodeMap _id;
kpeter@904
   119
kpeter@904
   120
    // Internal matrix representation of the graph
kpeter@904
   121
    BoolMatrix _gr;
kpeter@904
   122
    int _n;
kpeter@918
   123
    
kpeter@918
   124
    // Search options
kpeter@918
   125
    bool _delta_based_restart;
kpeter@918
   126
    int _restart_delta_limit;
kpeter@918
   127
 
kpeter@918
   128
    // Search limits
kpeter@918
   129
    int _iteration_limit;
kpeter@918
   130
    int _step_limit;
kpeter@918
   131
    int _size_limit;
kpeter@904
   132
kpeter@904
   133
    // The current clique
kpeter@904
   134
    BoolVector _clique;
kpeter@904
   135
    int _size;
kpeter@904
   136
kpeter@904
   137
    // The best clique found so far
kpeter@904
   138
    BoolVector _best_clique;
kpeter@904
   139
    int _best_size;
kpeter@904
   140
kpeter@904
   141
    // The "distances" of the nodes from the current clique.
kpeter@904
   142
    // _delta[u] is the number of nodes in the clique that are
kpeter@904
   143
    // not connected with u.
kpeter@904
   144
    IntVector _delta;
kpeter@904
   145
kpeter@904
   146
    // The current tabu set
kpeter@904
   147
    BoolVector _tabu;
kpeter@904
   148
kpeter@904
   149
    // Random number generator
kpeter@904
   150
    Random _rnd;
kpeter@904
   151
kpeter@904
   152
  private:
kpeter@904
   153
kpeter@904
   154
    // Implementation of the RANDOM node selection rule.
kpeter@904
   155
    class RandomSelectionRule
kpeter@904
   156
    {
kpeter@904
   157
    private:
kpeter@904
   158
kpeter@904
   159
      // References to the algorithm instance
kpeter@904
   160
      const BoolVector &_clique;
kpeter@904
   161
      const IntVector  &_delta;
kpeter@904
   162
      const BoolVector &_tabu;
kpeter@904
   163
      Random &_rnd;
kpeter@904
   164
kpeter@904
   165
      // Pivot rule data
kpeter@904
   166
      int _n;
kpeter@904
   167
kpeter@904
   168
    public:
kpeter@904
   169
kpeter@904
   170
      // Constructor
kpeter@904
   171
      RandomSelectionRule(GrossoLocatelliPullanMc &mc) :
kpeter@904
   172
        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
kpeter@904
   173
        _rnd(mc._rnd), _n(mc._n)
kpeter@904
   174
      {}
kpeter@904
   175
kpeter@904
   176
      // Return a node index for a feasible add move or -1 if no one exists
kpeter@904
   177
      int nextFeasibleAddNode() const {
kpeter@904
   178
        int start_node = _rnd[_n];
kpeter@904
   179
        for (int i = start_node; i != _n; i++) {
kpeter@904
   180
          if (_delta[i] == 0 && !_tabu[i]) return i;
kpeter@904
   181
        }
kpeter@904
   182
        for (int i = 0; i != start_node; i++) {
kpeter@904
   183
          if (_delta[i] == 0 && !_tabu[i]) return i;
kpeter@904
   184
        }
kpeter@904
   185
        return -1;
kpeter@904
   186
      }
kpeter@904
   187
kpeter@904
   188
      // Return a node index for a feasible swap move or -1 if no one exists
kpeter@904
   189
      int nextFeasibleSwapNode() const {
kpeter@904
   190
        int start_node = _rnd[_n];
kpeter@904
   191
        for (int i = start_node; i != _n; i++) {
kpeter@904
   192
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i]) return i;
kpeter@904
   193
        }
kpeter@904
   194
        for (int i = 0; i != start_node; i++) {
kpeter@904
   195
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i]) return i;
kpeter@904
   196
        }
kpeter@904
   197
        return -1;
kpeter@904
   198
      }
kpeter@904
   199
kpeter@904
   200
      // Return a node index for an add move or -1 if no one exists
kpeter@904
   201
      int nextAddNode() const {
kpeter@904
   202
        int start_node = _rnd[_n];
kpeter@904
   203
        for (int i = start_node; i != _n; i++) {
kpeter@904
   204
          if (_delta[i] == 0) return i;
kpeter@904
   205
        }
kpeter@904
   206
        for (int i = 0; i != start_node; i++) {
kpeter@904
   207
          if (_delta[i] == 0) return i;
kpeter@904
   208
        }
kpeter@904
   209
        return -1;
kpeter@904
   210
      }
kpeter@904
   211
kpeter@904
   212
      // Update internal data structures between stages (if necessary)
kpeter@904
   213
      void update() {}
kpeter@904
   214
kpeter@904
   215
    }; //class RandomSelectionRule
kpeter@904
   216
kpeter@904
   217
kpeter@904
   218
    // Implementation of the DEGREE_BASED node selection rule.
kpeter@904
   219
    class DegreeBasedSelectionRule
kpeter@904
   220
    {
kpeter@904
   221
    private:
kpeter@904
   222
kpeter@904
   223
      // References to the algorithm instance
kpeter@904
   224
      const BoolVector &_clique;
kpeter@904
   225
      const IntVector  &_delta;
kpeter@904
   226
      const BoolVector &_tabu;
kpeter@904
   227
      Random &_rnd;
kpeter@904
   228
kpeter@904
   229
      // Pivot rule data
kpeter@904
   230
      int _n;
kpeter@904
   231
      IntVector _deg;
kpeter@904
   232
kpeter@904
   233
    public:
kpeter@904
   234
kpeter@904
   235
      // Constructor
kpeter@904
   236
      DegreeBasedSelectionRule(GrossoLocatelliPullanMc &mc) :
kpeter@904
   237
        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
kpeter@904
   238
        _rnd(mc._rnd), _n(mc._n), _deg(_n)
kpeter@904
   239
      {
kpeter@904
   240
        for (int i = 0; i != _n; i++) {
kpeter@904
   241
          int d = 0;
kpeter@904
   242
          BoolVector &row = mc._gr[i];
kpeter@904
   243
          for (int j = 0; j != _n; j++) {
kpeter@904
   244
            if (row[j]) d++;
kpeter@904
   245
          }
kpeter@904
   246
          _deg[i] = d;
kpeter@904
   247
        }
kpeter@904
   248
      }
kpeter@904
   249
kpeter@904
   250
      // Return a node index for a feasible add move or -1 if no one exists
kpeter@904
   251
      int nextFeasibleAddNode() const {
kpeter@904
   252
        int start_node = _rnd[_n];
kpeter@904
   253
        int node = -1, max_deg = -1;
kpeter@904
   254
        for (int i = start_node; i != _n; i++) {
kpeter@904
   255
          if (_delta[i] == 0 && !_tabu[i] && _deg[i] > max_deg) {
kpeter@904
   256
            node = i;
kpeter@904
   257
            max_deg = _deg[i];
kpeter@904
   258
          }
kpeter@904
   259
        }
kpeter@904
   260
        for (int i = 0; i != start_node; i++) {
kpeter@904
   261
          if (_delta[i] == 0 && !_tabu[i] && _deg[i] > max_deg) {
kpeter@904
   262
            node = i;
kpeter@904
   263
            max_deg = _deg[i];
kpeter@904
   264
          }
kpeter@904
   265
        }
kpeter@904
   266
        return node;
kpeter@904
   267
      }
kpeter@904
   268
kpeter@904
   269
      // Return a node index for a feasible swap move or -1 if no one exists
kpeter@904
   270
      int nextFeasibleSwapNode() const {
kpeter@904
   271
        int start_node = _rnd[_n];
kpeter@904
   272
        int node = -1, max_deg = -1;
kpeter@904
   273
        for (int i = start_node; i != _n; i++) {
kpeter@904
   274
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
kpeter@904
   275
              _deg[i] > max_deg) {
kpeter@904
   276
            node = i;
kpeter@904
   277
            max_deg = _deg[i];
kpeter@904
   278
          }
kpeter@904
   279
        }
kpeter@904
   280
        for (int i = 0; i != start_node; i++) {
kpeter@904
   281
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
kpeter@904
   282
              _deg[i] > max_deg) {
kpeter@904
   283
            node = i;
kpeter@904
   284
            max_deg = _deg[i];
kpeter@904
   285
          }
kpeter@904
   286
        }
kpeter@904
   287
        return node;
kpeter@904
   288
      }
kpeter@904
   289
kpeter@904
   290
      // Return a node index for an add move or -1 if no one exists
kpeter@904
   291
      int nextAddNode() const {
kpeter@904
   292
        int start_node = _rnd[_n];
kpeter@904
   293
        int node = -1, max_deg = -1;
kpeter@904
   294
        for (int i = start_node; i != _n; i++) {
kpeter@904
   295
          if (_delta[i] == 0 && _deg[i] > max_deg) {
kpeter@904
   296
            node = i;
kpeter@904
   297
            max_deg = _deg[i];
kpeter@904
   298
          }
kpeter@904
   299
        }
kpeter@904
   300
        for (int i = 0; i != start_node; i++) {
kpeter@904
   301
          if (_delta[i] == 0 && _deg[i] > max_deg) {
kpeter@904
   302
            node = i;
kpeter@904
   303
            max_deg = _deg[i];
kpeter@904
   304
          }
kpeter@904
   305
        }
kpeter@904
   306
        return node;
kpeter@904
   307
      }
kpeter@904
   308
kpeter@904
   309
      // Update internal data structures between stages (if necessary)
kpeter@904
   310
      void update() {}
kpeter@904
   311
kpeter@904
   312
    }; //class DegreeBasedSelectionRule
kpeter@904
   313
kpeter@904
   314
kpeter@904
   315
    // Implementation of the PENALTY_BASED node selection rule.
kpeter@904
   316
    class PenaltyBasedSelectionRule
kpeter@904
   317
    {
kpeter@904
   318
    private:
kpeter@904
   319
kpeter@904
   320
      // References to the algorithm instance
kpeter@904
   321
      const BoolVector &_clique;
kpeter@904
   322
      const IntVector  &_delta;
kpeter@904
   323
      const BoolVector &_tabu;
kpeter@904
   324
      Random &_rnd;
kpeter@904
   325
kpeter@904
   326
      // Pivot rule data
kpeter@904
   327
      int _n;
kpeter@904
   328
      IntVector _penalty;
kpeter@904
   329
kpeter@904
   330
    public:
kpeter@904
   331
kpeter@904
   332
      // Constructor
kpeter@904
   333
      PenaltyBasedSelectionRule(GrossoLocatelliPullanMc &mc) :
kpeter@904
   334
        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
kpeter@904
   335
        _rnd(mc._rnd), _n(mc._n), _penalty(_n, 0)
kpeter@904
   336
      {}
kpeter@904
   337
kpeter@904
   338
      // Return a node index for a feasible add move or -1 if no one exists
kpeter@904
   339
      int nextFeasibleAddNode() const {
kpeter@904
   340
        int start_node = _rnd[_n];
kpeter@904
   341
        int node = -1, min_p = std::numeric_limits<int>::max();
kpeter@904
   342
        for (int i = start_node; i != _n; i++) {
kpeter@904
   343
          if (_delta[i] == 0 && !_tabu[i] && _penalty[i] < min_p) {
kpeter@904
   344
            node = i;
kpeter@904
   345
            min_p = _penalty[i];
kpeter@904
   346
          }
kpeter@904
   347
        }
kpeter@904
   348
        for (int i = 0; i != start_node; i++) {
kpeter@904
   349
          if (_delta[i] == 0 && !_tabu[i] && _penalty[i] < min_p) {
kpeter@904
   350
            node = i;
kpeter@904
   351
            min_p = _penalty[i];
kpeter@904
   352
          }
kpeter@904
   353
        }
kpeter@904
   354
        return node;
kpeter@904
   355
      }
kpeter@904
   356
kpeter@904
   357
      // Return a node index for a feasible swap move or -1 if no one exists
kpeter@904
   358
      int nextFeasibleSwapNode() const {
kpeter@904
   359
        int start_node = _rnd[_n];
kpeter@904
   360
        int node = -1, min_p = std::numeric_limits<int>::max();
kpeter@904
   361
        for (int i = start_node; i != _n; i++) {
kpeter@904
   362
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
kpeter@904
   363
              _penalty[i] < min_p) {
kpeter@904
   364
            node = i;
kpeter@904
   365
            min_p = _penalty[i];
kpeter@904
   366
          }
kpeter@904
   367
        }
kpeter@904
   368
        for (int i = 0; i != start_node; i++) {
kpeter@904
   369
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
kpeter@904
   370
              _penalty[i] < min_p) {
kpeter@904
   371
            node = i;
kpeter@904
   372
            min_p = _penalty[i];
kpeter@904
   373
          }
kpeter@904
   374
        }
kpeter@904
   375
        return node;
kpeter@904
   376
      }
kpeter@904
   377
kpeter@904
   378
      // Return a node index for an add move or -1 if no one exists
kpeter@904
   379
      int nextAddNode() const {
kpeter@904
   380
        int start_node = _rnd[_n];
kpeter@904
   381
        int node = -1, min_p = std::numeric_limits<int>::max();
kpeter@904
   382
        for (int i = start_node; i != _n; i++) {
kpeter@904
   383
          if (_delta[i] == 0 && _penalty[i] < min_p) {
kpeter@904
   384
            node = i;
kpeter@904
   385
            min_p = _penalty[i];
kpeter@904
   386
          }
kpeter@904
   387
        }
kpeter@904
   388
        for (int i = 0; i != start_node; i++) {
kpeter@904
   389
          if (_delta[i] == 0 && _penalty[i] < min_p) {
kpeter@904
   390
            node = i;
kpeter@904
   391
            min_p = _penalty[i];
kpeter@904
   392
          }
kpeter@904
   393
        }
kpeter@904
   394
        return node;
kpeter@904
   395
      }
kpeter@904
   396
kpeter@904
   397
      // Update internal data structures between stages (if necessary)
kpeter@904
   398
      void update() {}
kpeter@904
   399
kpeter@904
   400
    }; //class PenaltyBasedSelectionRule
kpeter@904
   401
kpeter@904
   402
  public:
kpeter@904
   403
kpeter@904
   404
    /// \brief Constructor.
kpeter@904
   405
    ///
kpeter@904
   406
    /// Constructor.
kpeter@904
   407
    /// The global \ref rnd "random number generator instance" is used
kpeter@904
   408
    /// during the algorithm.
kpeter@904
   409
    ///
kpeter@904
   410
    /// \param graph The undirected graph the algorithm runs on.
kpeter@904
   411
    GrossoLocatelliPullanMc(const GR& graph) :
kpeter@904
   412
      _graph(graph), _id(_graph), _rnd(rnd)
kpeter@918
   413
    {
kpeter@918
   414
      initOptions();
kpeter@918
   415
    }
kpeter@904
   416
kpeter@904
   417
    /// \brief Constructor with random seed.
kpeter@904
   418
    ///
kpeter@904
   419
    /// Constructor with random seed.
kpeter@904
   420
    ///
kpeter@904
   421
    /// \param graph The undirected graph the algorithm runs on.
kpeter@904
   422
    /// \param seed Seed value for the internal random number generator
kpeter@904
   423
    /// that is used during the algorithm.
kpeter@904
   424
    GrossoLocatelliPullanMc(const GR& graph, int seed) :
kpeter@904
   425
      _graph(graph), _id(_graph), _rnd(seed)
kpeter@918
   426
    {
kpeter@918
   427
      initOptions();
kpeter@918
   428
    }
kpeter@904
   429
kpeter@904
   430
    /// \brief Constructor with random number generator.
kpeter@904
   431
    ///
kpeter@904
   432
    /// Constructor with random number generator.
kpeter@904
   433
    ///
kpeter@904
   434
    /// \param graph The undirected graph the algorithm runs on.
kpeter@904
   435
    /// \param random A random number generator that is used during the
kpeter@904
   436
    /// algorithm.
kpeter@904
   437
    GrossoLocatelliPullanMc(const GR& graph, const Random& random) :
kpeter@904
   438
      _graph(graph), _id(_graph), _rnd(random)
kpeter@918
   439
    {
kpeter@918
   440
      initOptions();
kpeter@918
   441
    }
kpeter@904
   442
kpeter@904
   443
    /// \name Execution Control
kpeter@918
   444
    /// The \ref run() function can be used to execute the algorithm.\n
kpeter@918
   445
    /// The functions \ref iterationLimit(int), \ref stepLimit(int), and 
kpeter@918
   446
    /// \ref sizeLimit(int) can be used to specify various limits for the
kpeter@918
   447
    /// search process.
kpeter@918
   448
    
kpeter@904
   449
    /// @{
kpeter@918
   450
    
kpeter@918
   451
    /// \brief Sets the maximum number of iterations.
kpeter@918
   452
    ///
kpeter@918
   453
    /// This function sets the maximum number of iterations.
kpeter@918
   454
    /// Each iteration of the algorithm finds a maximal clique (but not
kpeter@918
   455
    /// necessarily the largest one) by performing several search steps
kpeter@918
   456
    /// (node selections).
kpeter@918
   457
    ///
kpeter@918
   458
    /// This limit controls the running time and the success of the
kpeter@918
   459
    /// algorithm. For larger values, the algorithm runs slower, but it more
kpeter@918
   460
    /// likely finds larger cliques. For smaller values, the algorithm is
kpeter@918
   461
    /// faster but probably gives worse results.
kpeter@918
   462
    /// 
kpeter@918
   463
    /// The default value is \c 1000.
kpeter@918
   464
    /// \c -1 means that number of iterations is not limited.
kpeter@918
   465
    ///
kpeter@918
   466
    /// \warning You should specify a reasonable limit for the number of
kpeter@918
   467
    /// iterations and/or the number of search steps.
kpeter@918
   468
    ///
kpeter@918
   469
    /// \return <tt>(*this)</tt>
kpeter@918
   470
    ///
kpeter@918
   471
    /// \sa stepLimit(int)
kpeter@918
   472
    /// \sa sizeLimit(int)
kpeter@918
   473
    GrossoLocatelliPullanMc& iterationLimit(int limit) {
kpeter@918
   474
      _iteration_limit = limit;
kpeter@918
   475
      return *this;
kpeter@918
   476
    }
kpeter@918
   477
    
kpeter@918
   478
    /// \brief Sets the maximum number of search steps.
kpeter@918
   479
    ///
kpeter@918
   480
    /// This function sets the maximum number of elementary search steps.
kpeter@918
   481
    /// Each iteration of the algorithm finds a maximal clique (but not
kpeter@918
   482
    /// necessarily the largest one) by performing several search steps
kpeter@918
   483
    /// (node selections).
kpeter@918
   484
    ///
kpeter@918
   485
    /// This limit controls the running time and the success of the
kpeter@918
   486
    /// algorithm. For larger values, the algorithm runs slower, but it more
kpeter@918
   487
    /// likely finds larger cliques. For smaller values, the algorithm is
kpeter@918
   488
    /// faster but probably gives worse results.
kpeter@918
   489
    /// 
kpeter@918
   490
    /// The default value is \c -1, which means that number of steps
kpeter@918
   491
    /// is not limited explicitly. However, the number of iterations is
kpeter@918
   492
    /// limited and each iteration performs a finite number of search steps.
kpeter@918
   493
    ///
kpeter@918
   494
    /// \warning You should specify a reasonable limit for the number of
kpeter@918
   495
    /// iterations and/or the number of search steps.
kpeter@918
   496
    ///
kpeter@918
   497
    /// \return <tt>(*this)</tt>
kpeter@918
   498
    ///
kpeter@918
   499
    /// \sa iterationLimit(int)
kpeter@918
   500
    /// \sa sizeLimit(int)
kpeter@918
   501
    GrossoLocatelliPullanMc& stepLimit(int limit) {
kpeter@918
   502
      _step_limit = limit;
kpeter@918
   503
      return *this;
kpeter@918
   504
    }
kpeter@918
   505
    
kpeter@918
   506
    /// \brief Sets the desired clique size.
kpeter@918
   507
    ///
kpeter@918
   508
    /// This function sets the desired clique size that serves as a search
kpeter@918
   509
    /// limit. If a clique of this size (or a larger one) is found, then the
kpeter@918
   510
    /// algorithm terminates.
kpeter@918
   511
    /// 
kpeter@918
   512
    /// This function is especially useful if you know an exact upper bound
kpeter@918
   513
    /// for the size of the cliques in the graph or if any clique above 
kpeter@918
   514
    /// a certain size limit is sufficient for your application.
kpeter@918
   515
    /// 
kpeter@918
   516
    /// The default value is \c -1, which means that the size limit is set to
kpeter@918
   517
    /// the number of nodes in the graph.
kpeter@918
   518
    ///
kpeter@918
   519
    /// \return <tt>(*this)</tt>
kpeter@918
   520
    ///
kpeter@918
   521
    /// \sa iterationLimit(int)
kpeter@918
   522
    /// \sa stepLimit(int)
kpeter@918
   523
    GrossoLocatelliPullanMc& sizeLimit(int limit) {
kpeter@918
   524
      _size_limit = limit;
kpeter@918
   525
      return *this;
kpeter@918
   526
    }
kpeter@918
   527
    
kpeter@918
   528
    /// \brief The maximum number of iterations.
kpeter@918
   529
    ///
kpeter@918
   530
    /// This function gives back the maximum number of iterations.
kpeter@918
   531
    /// \c -1 means that no limit is specified.
kpeter@918
   532
    ///
kpeter@918
   533
    /// \sa iterationLimit(int)
kpeter@918
   534
    int iterationLimit() const {
kpeter@918
   535
      return _iteration_limit;
kpeter@918
   536
    }
kpeter@918
   537
    
kpeter@918
   538
    /// \brief The maximum number of search steps.
kpeter@918
   539
    ///
kpeter@918
   540
    /// This function gives back the maximum number of search steps.
kpeter@918
   541
    /// \c -1 means that no limit is specified.
kpeter@918
   542
    ///
kpeter@918
   543
    /// \sa stepLimit(int)
kpeter@918
   544
    int stepLimit() const {
kpeter@918
   545
      return _step_limit;
kpeter@918
   546
    }
kpeter@918
   547
    
kpeter@918
   548
    /// \brief The desired clique size.
kpeter@918
   549
    ///
kpeter@918
   550
    /// This function gives back the desired clique size that serves as a
kpeter@918
   551
    /// search limit. \c -1 means that this limit is set to the number of
kpeter@918
   552
    /// nodes in the graph.
kpeter@918
   553
    ///
kpeter@918
   554
    /// \sa sizeLimit(int)
kpeter@918
   555
    int sizeLimit() const {
kpeter@918
   556
      return _size_limit;
kpeter@918
   557
    }
kpeter@904
   558
kpeter@904
   559
    /// \brief Runs the algorithm.
kpeter@904
   560
    ///
kpeter@918
   561
    /// This function runs the algorithm. If one of the specified limits
kpeter@918
   562
    /// is reached, the search process terminates.
kpeter@904
   563
    ///
kpeter@904
   564
    /// \param rule The node selection rule. For more information, see
kpeter@904
   565
    /// \ref SelectionRule.
kpeter@904
   566
    ///
kpeter@918
   567
    /// \return The termination cause of the search. For more information,
kpeter@918
   568
    /// see \ref TerminationCause.
kpeter@918
   569
    TerminationCause run(SelectionRule rule = PENALTY_BASED)
kpeter@904
   570
    {
kpeter@904
   571
      init();
kpeter@904
   572
      switch (rule) {
kpeter@904
   573
        case RANDOM:
kpeter@918
   574
          return start<RandomSelectionRule>();
kpeter@904
   575
        case DEGREE_BASED:
kpeter@918
   576
          return start<DegreeBasedSelectionRule>();
kpeter@918
   577
        default:
kpeter@918
   578
          return start<PenaltyBasedSelectionRule>();
kpeter@904
   579
      }
kpeter@904
   580
    }
kpeter@904
   581
kpeter@904
   582
    /// @}
kpeter@904
   583
kpeter@904
   584
    /// \name Query Functions
kpeter@918
   585
    /// The results of the algorithm can be obtained using these functions.\n
kpeter@918
   586
    /// The run() function must be called before using them. 
kpeter@918
   587
kpeter@904
   588
    /// @{
kpeter@904
   589
kpeter@904
   590
    /// \brief The size of the found clique
kpeter@904
   591
    ///
kpeter@904
   592
    /// This function returns the size of the found clique.
kpeter@904
   593
    ///
kpeter@904
   594
    /// \pre run() must be called before using this function.
kpeter@904
   595
    int cliqueSize() const {
kpeter@904
   596
      return _best_size;
kpeter@904
   597
    }
kpeter@904
   598
kpeter@904
   599
    /// \brief Gives back the found clique in a \c bool node map
kpeter@904
   600
    ///
kpeter@904
   601
    /// This function gives back the characteristic vector of the found
kpeter@904
   602
    /// clique in the given node map.
kpeter@904
   603
    /// It must be a \ref concepts::WriteMap "writable" node map with
kpeter@904
   604
    /// \c bool (or convertible) value type.
kpeter@904
   605
    ///
kpeter@904
   606
    /// \pre run() must be called before using this function.
kpeter@904
   607
    template <typename CliqueMap>
kpeter@904
   608
    void cliqueMap(CliqueMap &map) const {
kpeter@904
   609
      for (NodeIt n(_graph); n != INVALID; ++n) {
kpeter@904
   610
        map[n] = static_cast<bool>(_best_clique[_id[n]]);
kpeter@904
   611
      }
kpeter@904
   612
    }
kpeter@904
   613
kpeter@904
   614
    /// \brief Iterator to list the nodes of the found clique
kpeter@904
   615
    ///
kpeter@904
   616
    /// This iterator class lists the nodes of the found clique.
kpeter@904
   617
    /// Before using it, you must allocate a GrossoLocatelliPullanMc instance
kpeter@904
   618
    /// and call its \ref GrossoLocatelliPullanMc::run() "run()" method.
kpeter@904
   619
    ///
kpeter@904
   620
    /// The following example prints out the IDs of the nodes in the found
kpeter@904
   621
    /// clique.
kpeter@904
   622
    /// \code
kpeter@904
   623
    ///   GrossoLocatelliPullanMc<Graph> mc(g);
kpeter@904
   624
    ///   mc.run();
kpeter@904
   625
    ///   for (GrossoLocatelliPullanMc<Graph>::CliqueNodeIt n(mc);
kpeter@904
   626
    ///        n != INVALID; ++n)
kpeter@904
   627
    ///   {
kpeter@904
   628
    ///     std::cout << g.id(n) << std::endl;
kpeter@904
   629
    ///   }
kpeter@904
   630
    /// \endcode
kpeter@904
   631
    class CliqueNodeIt
kpeter@904
   632
    {
kpeter@904
   633
    private:
kpeter@904
   634
      NodeIt _it;
kpeter@904
   635
      BoolNodeMap _map;
kpeter@904
   636
kpeter@904
   637
    public:
kpeter@904
   638
kpeter@904
   639
      /// Constructor
kpeter@904
   640
kpeter@904
   641
      /// Constructor.
kpeter@904
   642
      /// \param mc The algorithm instance.
kpeter@904
   643
      CliqueNodeIt(const GrossoLocatelliPullanMc &mc)
kpeter@904
   644
       : _map(mc._graph)
kpeter@904
   645
      {
kpeter@904
   646
        mc.cliqueMap(_map);
kpeter@904
   647
        for (_it = NodeIt(mc._graph); _it != INVALID && !_map[_it]; ++_it) ;
kpeter@904
   648
      }
kpeter@904
   649
kpeter@904
   650
      /// Conversion to \c Node
kpeter@904
   651
      operator Node() const { return _it; }
kpeter@904
   652
kpeter@904
   653
      bool operator==(Invalid) const { return _it == INVALID; }
kpeter@904
   654
      bool operator!=(Invalid) const { return _it != INVALID; }
kpeter@904
   655
kpeter@904
   656
      /// Next node
kpeter@904
   657
      CliqueNodeIt &operator++() {
kpeter@904
   658
        for (++_it; _it != INVALID && !_map[_it]; ++_it) ;
kpeter@904
   659
        return *this;
kpeter@904
   660
      }
kpeter@904
   661
kpeter@904
   662
      /// Postfix incrementation
kpeter@904
   663
kpeter@904
   664
      /// Postfix incrementation.
kpeter@904
   665
      ///
kpeter@904
   666
      /// \warning This incrementation returns a \c Node, not a
kpeter@904
   667
      /// \c CliqueNodeIt as one may expect.
kpeter@904
   668
      typename GR::Node operator++(int) {
kpeter@904
   669
        Node n=*this;
kpeter@904
   670
        ++(*this);
kpeter@904
   671
        return n;
kpeter@904
   672
      }
kpeter@904
   673
kpeter@904
   674
    };
kpeter@904
   675
kpeter@904
   676
    /// @}
kpeter@904
   677
kpeter@904
   678
  private:
kpeter@918
   679
  
kpeter@918
   680
    // Initialize search options and limits
kpeter@918
   681
    void initOptions() {
kpeter@918
   682
      // Search options
kpeter@918
   683
      _delta_based_restart = true;
kpeter@918
   684
      _restart_delta_limit = 4;
kpeter@918
   685
     
kpeter@918
   686
      // Search limits
kpeter@918
   687
      _iteration_limit = 1000;
kpeter@918
   688
      _step_limit = -1;             // this is disabled by default
kpeter@918
   689
      _size_limit = -1;             // this is disabled by default
kpeter@918
   690
    }
kpeter@904
   691
kpeter@904
   692
    // Adds a node to the current clique
kpeter@904
   693
    void addCliqueNode(int u) {
kpeter@904
   694
      if (_clique[u]) return;
kpeter@904
   695
      _clique[u] = true;
kpeter@904
   696
      _size++;
kpeter@904
   697
      BoolVector &row = _gr[u];
kpeter@904
   698
      for (int i = 0; i != _n; i++) {
kpeter@904
   699
        if (!row[i]) _delta[i]++;
kpeter@904
   700
      }
kpeter@904
   701
    }
kpeter@904
   702
kpeter@904
   703
    // Removes a node from the current clique
kpeter@904
   704
    void delCliqueNode(int u) {
kpeter@904
   705
      if (!_clique[u]) return;
kpeter@904
   706
      _clique[u] = false;
kpeter@904
   707
      _size--;
kpeter@904
   708
      BoolVector &row = _gr[u];
kpeter@904
   709
      for (int i = 0; i != _n; i++) {
kpeter@904
   710
        if (!row[i]) _delta[i]--;
kpeter@904
   711
      }
kpeter@904
   712
    }
kpeter@904
   713
kpeter@904
   714
    // Initialize data structures
kpeter@904
   715
    void init() {
kpeter@904
   716
      _n = countNodes(_graph);
kpeter@904
   717
      int ui = 0;
kpeter@904
   718
      for (NodeIt u(_graph); u != INVALID; ++u) {
kpeter@904
   719
        _id[u] = ui++;
kpeter@904
   720
      }
kpeter@904
   721
      _gr.clear();
kpeter@904
   722
      _gr.resize(_n, BoolVector(_n, false));
kpeter@904
   723
      ui = 0;
kpeter@904
   724
      for (NodeIt u(_graph); u != INVALID; ++u) {
kpeter@904
   725
        for (IncEdgeIt e(_graph, u); e != INVALID; ++e) {
kpeter@904
   726
          int vi = _id[_graph.runningNode(e)];
kpeter@904
   727
          _gr[ui][vi] = true;
kpeter@904
   728
          _gr[vi][ui] = true;
kpeter@904
   729
        }
kpeter@904
   730
        ++ui;
kpeter@904
   731
      }
kpeter@904
   732
kpeter@904
   733
      _clique.clear();
kpeter@904
   734
      _clique.resize(_n, false);
kpeter@904
   735
      _size = 0;
kpeter@904
   736
      _best_clique.clear();
kpeter@904
   737
      _best_clique.resize(_n, false);
kpeter@904
   738
      _best_size = 0;
kpeter@904
   739
      _delta.clear();
kpeter@904
   740
      _delta.resize(_n, 0);
kpeter@904
   741
      _tabu.clear();
kpeter@904
   742
      _tabu.resize(_n, false);
kpeter@904
   743
    }
kpeter@904
   744
kpeter@904
   745
    // Executes the algorithm
kpeter@904
   746
    template <typename SelectionRuleImpl>
kpeter@918
   747
    TerminationCause start() {
kpeter@918
   748
      if (_n == 0) return SIZE_LIMIT;
kpeter@904
   749
      if (_n == 1) {
kpeter@904
   750
        _best_clique[0] = true;
kpeter@904
   751
        _best_size = 1;
kpeter@918
   752
        return SIZE_LIMIT;
kpeter@904
   753
      }
kpeter@904
   754
kpeter@918
   755
      // Iterated local search algorithm
kpeter@918
   756
      const int max_size = _size_limit >= 0 ? _size_limit : _n;
kpeter@918
   757
      const int max_restart = _iteration_limit >= 0 ?
kpeter@918
   758
        _iteration_limit : std::numeric_limits<int>::max();
kpeter@918
   759
      const int max_select = _step_limit >= 0 ?
kpeter@918
   760
        _step_limit : std::numeric_limits<int>::max();
kpeter@918
   761
kpeter@904
   762
      SelectionRuleImpl sel_method(*this);
kpeter@918
   763
      int select = 0, restart = 0;
kpeter@904
   764
      IntVector restart_nodes;
kpeter@918
   765
      while (select < max_select && restart < max_restart) {
kpeter@904
   766
kpeter@904
   767
        // Perturbation/restart
kpeter@918
   768
        restart++;
kpeter@918
   769
        if (_delta_based_restart) {
kpeter@904
   770
          restart_nodes.clear();
kpeter@904
   771
          for (int i = 0; i != _n; i++) {
kpeter@918
   772
            if (_delta[i] >= _restart_delta_limit)
kpeter@904
   773
              restart_nodes.push_back(i);
kpeter@904
   774
          }
kpeter@904
   775
        }
kpeter@904
   776
        int rs_node = -1;
kpeter@904
   777
        if (restart_nodes.size() > 0) {
kpeter@904
   778
          rs_node = restart_nodes[_rnd[restart_nodes.size()]];
kpeter@904
   779
        } else {
kpeter@904
   780
          rs_node = _rnd[_n];
kpeter@904
   781
        }
kpeter@904
   782
        BoolVector &row = _gr[rs_node];
kpeter@904
   783
        for (int i = 0; i != _n; i++) {
kpeter@904
   784
          if (_clique[i] && !row[i]) delCliqueNode(i);
kpeter@904
   785
        }
kpeter@904
   786
        addCliqueNode(rs_node);
kpeter@904
   787
kpeter@904
   788
        // Local search
kpeter@904
   789
        _tabu.clear();
kpeter@904
   790
        _tabu.resize(_n, false);
kpeter@904
   791
        bool tabu_empty = true;
kpeter@904
   792
        int max_swap = _size;
kpeter@904
   793
        while (select < max_select) {
kpeter@904
   794
          select++;
kpeter@904
   795
          int u;
kpeter@904
   796
          if ((u = sel_method.nextFeasibleAddNode()) != -1) {
kpeter@904
   797
            // Feasible add move
kpeter@904
   798
            addCliqueNode(u);
kpeter@904
   799
            if (tabu_empty) max_swap = _size;
kpeter@904
   800
          }
kpeter@904
   801
          else if ((u = sel_method.nextFeasibleSwapNode()) != -1) {
kpeter@904
   802
            // Feasible swap move
kpeter@904
   803
            int v = -1;
kpeter@904
   804
            BoolVector &row = _gr[u];
kpeter@904
   805
            for (int i = 0; i != _n; i++) {
kpeter@904
   806
              if (_clique[i] && !row[i]) {
kpeter@904
   807
                v = i;
kpeter@904
   808
                break;
kpeter@904
   809
              }
kpeter@904
   810
            }
kpeter@904
   811
            addCliqueNode(u);
kpeter@904
   812
            delCliqueNode(v);
kpeter@904
   813
            _tabu[v] = true;
kpeter@904
   814
            tabu_empty = false;
kpeter@904
   815
            if (--max_swap <= 0) break;
kpeter@904
   816
          }
kpeter@904
   817
          else if ((u = sel_method.nextAddNode()) != -1) {
kpeter@904
   818
            // Non-feasible add move
kpeter@904
   819
            addCliqueNode(u);
kpeter@904
   820
          }
kpeter@904
   821
          else break;
kpeter@904
   822
        }
kpeter@904
   823
        if (_size > _best_size) {
kpeter@904
   824
          _best_clique = _clique;
kpeter@904
   825
          _best_size = _size;
kpeter@918
   826
          if (_best_size >= max_size) return SIZE_LIMIT;
kpeter@904
   827
        }
kpeter@904
   828
        sel_method.update();
kpeter@904
   829
      }
kpeter@904
   830
kpeter@918
   831
      return (restart >= max_restart ? ITERATION_LIMIT : STEP_LIMIT);
kpeter@904
   832
    }
kpeter@904
   833
kpeter@904
   834
  }; //class GrossoLocatelliPullanMc
kpeter@904
   835
kpeter@904
   836
  ///@}
kpeter@904
   837
kpeter@904
   838
} //namespace lemon
kpeter@904
   839
kpeter@904
   840
#endif //LEMON_GROSSO_LOCATELLI_PULLAN_MC_H