lemon/grosso_locatelli_pullan_mc.h
changeset 1177 3c00344f49c9
parent 1053 1c978b5bcc65
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lemon/grosso_locatelli_pullan_mc.h	Wed Oct 17 19:14:07 2018 +0200
     1.3 @@ -0,0 +1,840 @@
     1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     1.5 + *
     1.6 + * This file is a part of LEMON, a generic C++ optimization library.
     1.7 + *
     1.8 + * Copyright (C) 2003-2013
     1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    1.11 + *
    1.12 + * Permission to use, modify and distribute this software is granted
    1.13 + * provided that this copyright notice appears in all copies. For
    1.14 + * precise terms see the accompanying LICENSE file.
    1.15 + *
    1.16 + * This software is provided "AS IS" with no warranty of any kind,
    1.17 + * express or implied, and with no claim as to its suitability for any
    1.18 + * purpose.
    1.19 + *
    1.20 + */
    1.21 +
    1.22 +#ifndef LEMON_GROSSO_LOCATELLI_PULLAN_MC_H
    1.23 +#define LEMON_GROSSO_LOCATELLI_PULLAN_MC_H
    1.24 +
    1.25 +/// \ingroup approx_algs
    1.26 +///
    1.27 +/// \file
    1.28 +/// \brief The iterated local search algorithm of Grosso, Locatelli, and Pullan
    1.29 +/// for the maximum clique problem
    1.30 +
    1.31 +#include <vector>
    1.32 +#include <limits>
    1.33 +#include <lemon/core.h>
    1.34 +#include <lemon/random.h>
    1.35 +
    1.36 +namespace lemon {
    1.37 +
    1.38 +  /// \addtogroup approx_algs
    1.39 +  /// @{
    1.40 +
    1.41 +  /// \brief Implementation of the iterated local search algorithm of Grosso,
    1.42 +  /// Locatelli, and Pullan for the maximum clique problem
    1.43 +  ///
    1.44 +  /// \ref GrossoLocatelliPullanMc implements the iterated local search
    1.45 +  /// algorithm of Grosso, Locatelli, and Pullan for solving the \e maximum
    1.46 +  /// \e clique \e problem \cite grosso08maxclique.
    1.47 +  /// It is to find the largest complete subgraph (\e clique) in an
    1.48 +  /// undirected graph, i.e., the largest set of nodes where each
    1.49 +  /// pair of nodes is connected.
    1.50 +  ///
    1.51 +  /// This class provides a simple but highly efficient and robust heuristic
    1.52 +  /// method that quickly finds a quite large clique, but not necessarily the
    1.53 +  /// largest one.
    1.54 +  /// The algorithm performs a certain number of iterations to find several
    1.55 +  /// cliques and selects the largest one among them. Various limits can be
    1.56 +  /// specified to control the running time and the effectiveness of the
    1.57 +  /// search process.
    1.58 +  ///
    1.59 +  /// \tparam GR The undirected graph type the algorithm runs on.
    1.60 +  ///
    1.61 +  /// \note %GrossoLocatelliPullanMc provides three different node selection
    1.62 +  /// rules, from which the most powerful one is used by default.
    1.63 +  /// For more information, see \ref SelectionRule.
    1.64 +  template <typename GR>
    1.65 +  class GrossoLocatelliPullanMc
    1.66 +  {
    1.67 +  public:
    1.68 +
    1.69 +    /// \brief Constants for specifying the node selection rule.
    1.70 +    ///
    1.71 +    /// Enum type containing constants for specifying the node selection rule
    1.72 +    /// for the \ref run() function.
    1.73 +    ///
    1.74 +    /// During the algorithm, nodes are selected for addition to the current
    1.75 +    /// clique according to the applied rule.
    1.76 +    /// In general, the PENALTY_BASED rule turned out to be the most powerful
    1.77 +    /// and the most robust, thus it is the default option.
    1.78 +    /// However, another selection rule can be specified using the \ref run()
    1.79 +    /// function with the proper parameter.
    1.80 +    enum SelectionRule {
    1.81 +
    1.82 +      /// A node is selected randomly without any evaluation at each step.
    1.83 +      RANDOM,
    1.84 +
    1.85 +      /// A node of maximum degree is selected randomly at each step.
    1.86 +      DEGREE_BASED,
    1.87 +
    1.88 +      /// A node of minimum penalty is selected randomly at each step.
    1.89 +      /// The node penalties are updated adaptively after each stage of the
    1.90 +      /// search process.
    1.91 +      PENALTY_BASED
    1.92 +    };
    1.93 +
    1.94 +    /// \brief Constants for the causes of search termination.
    1.95 +    ///
    1.96 +    /// Enum type containing constants for the different causes of search
    1.97 +    /// termination. The \ref run() function returns one of these values.
    1.98 +    enum TerminationCause {
    1.99 +
   1.100 +      /// The iteration count limit is reached.
   1.101 +      ITERATION_LIMIT,
   1.102 +
   1.103 +      /// The step count limit is reached.
   1.104 +      STEP_LIMIT,
   1.105 +
   1.106 +      /// The clique size limit is reached.
   1.107 +      SIZE_LIMIT
   1.108 +    };
   1.109 +
   1.110 +  private:
   1.111 +
   1.112 +    TEMPLATE_GRAPH_TYPEDEFS(GR);
   1.113 +
   1.114 +    typedef std::vector<int> IntVector;
   1.115 +    typedef std::vector<char> BoolVector;
   1.116 +    typedef std::vector<BoolVector> BoolMatrix;
   1.117 +    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
   1.118 +
   1.119 +    // The underlying graph
   1.120 +    const GR &_graph;
   1.121 +    IntNodeMap _id;
   1.122 +
   1.123 +    // Internal matrix representation of the graph
   1.124 +    BoolMatrix _gr;
   1.125 +    int _n;
   1.126 +
   1.127 +    // Search options
   1.128 +    bool _delta_based_restart;
   1.129 +    int _restart_delta_limit;
   1.130 +
   1.131 +    // Search limits
   1.132 +    int _iteration_limit;
   1.133 +    int _step_limit;
   1.134 +    int _size_limit;
   1.135 +
   1.136 +    // The current clique
   1.137 +    BoolVector _clique;
   1.138 +    int _size;
   1.139 +
   1.140 +    // The best clique found so far
   1.141 +    BoolVector _best_clique;
   1.142 +    int _best_size;
   1.143 +
   1.144 +    // The "distances" of the nodes from the current clique.
   1.145 +    // _delta[u] is the number of nodes in the clique that are
   1.146 +    // not connected with u.
   1.147 +    IntVector _delta;
   1.148 +
   1.149 +    // The current tabu set
   1.150 +    BoolVector _tabu;
   1.151 +
   1.152 +    // Random number generator
   1.153 +    Random _rnd;
   1.154 +
   1.155 +  private:
   1.156 +
   1.157 +    // Implementation of the RANDOM node selection rule.
   1.158 +    class RandomSelectionRule
   1.159 +    {
   1.160 +    private:
   1.161 +
   1.162 +      // References to the algorithm instance
   1.163 +      const BoolVector &_clique;
   1.164 +      const IntVector  &_delta;
   1.165 +      const BoolVector &_tabu;
   1.166 +      Random &_rnd;
   1.167 +
   1.168 +      // Pivot rule data
   1.169 +      int _n;
   1.170 +
   1.171 +    public:
   1.172 +
   1.173 +      // Constructor
   1.174 +      RandomSelectionRule(GrossoLocatelliPullanMc &mc) :
   1.175 +        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
   1.176 +        _rnd(mc._rnd), _n(mc._n)
   1.177 +      {}
   1.178 +
   1.179 +      // Return a node index for a feasible add move or -1 if no one exists
   1.180 +      int nextFeasibleAddNode() const {
   1.181 +        int start_node = _rnd[_n];
   1.182 +        for (int i = start_node; i != _n; i++) {
   1.183 +          if (_delta[i] == 0 && !_tabu[i]) return i;
   1.184 +        }
   1.185 +        for (int i = 0; i != start_node; i++) {
   1.186 +          if (_delta[i] == 0 && !_tabu[i]) return i;
   1.187 +        }
   1.188 +        return -1;
   1.189 +      }
   1.190 +
   1.191 +      // Return a node index for a feasible swap move or -1 if no one exists
   1.192 +      int nextFeasibleSwapNode() const {
   1.193 +        int start_node = _rnd[_n];
   1.194 +        for (int i = start_node; i != _n; i++) {
   1.195 +          if (!_clique[i] && _delta[i] == 1 && !_tabu[i]) return i;
   1.196 +        }
   1.197 +        for (int i = 0; i != start_node; i++) {
   1.198 +          if (!_clique[i] && _delta[i] == 1 && !_tabu[i]) return i;
   1.199 +        }
   1.200 +        return -1;
   1.201 +      }
   1.202 +
   1.203 +      // Return a node index for an add move or -1 if no one exists
   1.204 +      int nextAddNode() const {
   1.205 +        int start_node = _rnd[_n];
   1.206 +        for (int i = start_node; i != _n; i++) {
   1.207 +          if (_delta[i] == 0) return i;
   1.208 +        }
   1.209 +        for (int i = 0; i != start_node; i++) {
   1.210 +          if (_delta[i] == 0) return i;
   1.211 +        }
   1.212 +        return -1;
   1.213 +      }
   1.214 +
   1.215 +      // Update internal data structures between stages (if necessary)
   1.216 +      void update() {}
   1.217 +
   1.218 +    }; //class RandomSelectionRule
   1.219 +
   1.220 +
   1.221 +    // Implementation of the DEGREE_BASED node selection rule.
   1.222 +    class DegreeBasedSelectionRule
   1.223 +    {
   1.224 +    private:
   1.225 +
   1.226 +      // References to the algorithm instance
   1.227 +      const BoolVector &_clique;
   1.228 +      const IntVector  &_delta;
   1.229 +      const BoolVector &_tabu;
   1.230 +      Random &_rnd;
   1.231 +
   1.232 +      // Pivot rule data
   1.233 +      int _n;
   1.234 +      IntVector _deg;
   1.235 +
   1.236 +    public:
   1.237 +
   1.238 +      // Constructor
   1.239 +      DegreeBasedSelectionRule(GrossoLocatelliPullanMc &mc) :
   1.240 +        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
   1.241 +        _rnd(mc._rnd), _n(mc._n), _deg(_n)
   1.242 +      {
   1.243 +        for (int i = 0; i != _n; i++) {
   1.244 +          int d = 0;
   1.245 +          BoolVector &row = mc._gr[i];
   1.246 +          for (int j = 0; j != _n; j++) {
   1.247 +            if (row[j]) d++;
   1.248 +          }
   1.249 +          _deg[i] = d;
   1.250 +        }
   1.251 +      }
   1.252 +
   1.253 +      // Return a node index for a feasible add move or -1 if no one exists
   1.254 +      int nextFeasibleAddNode() const {
   1.255 +        int start_node = _rnd[_n];
   1.256 +        int node = -1, max_deg = -1;
   1.257 +        for (int i = start_node; i != _n; i++) {
   1.258 +          if (_delta[i] == 0 && !_tabu[i] && _deg[i] > max_deg) {
   1.259 +            node = i;
   1.260 +            max_deg = _deg[i];
   1.261 +          }
   1.262 +        }
   1.263 +        for (int i = 0; i != start_node; i++) {
   1.264 +          if (_delta[i] == 0 && !_tabu[i] && _deg[i] > max_deg) {
   1.265 +            node = i;
   1.266 +            max_deg = _deg[i];
   1.267 +          }
   1.268 +        }
   1.269 +        return node;
   1.270 +      }
   1.271 +
   1.272 +      // Return a node index for a feasible swap move or -1 if no one exists
   1.273 +      int nextFeasibleSwapNode() const {
   1.274 +        int start_node = _rnd[_n];
   1.275 +        int node = -1, max_deg = -1;
   1.276 +        for (int i = start_node; i != _n; i++) {
   1.277 +          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
   1.278 +              _deg[i] > max_deg) {
   1.279 +            node = i;
   1.280 +            max_deg = _deg[i];
   1.281 +          }
   1.282 +        }
   1.283 +        for (int i = 0; i != start_node; i++) {
   1.284 +          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
   1.285 +              _deg[i] > max_deg) {
   1.286 +            node = i;
   1.287 +            max_deg = _deg[i];
   1.288 +          }
   1.289 +        }
   1.290 +        return node;
   1.291 +      }
   1.292 +
   1.293 +      // Return a node index for an add move or -1 if no one exists
   1.294 +      int nextAddNode() const {
   1.295 +        int start_node = _rnd[_n];
   1.296 +        int node = -1, max_deg = -1;
   1.297 +        for (int i = start_node; i != _n; i++) {
   1.298 +          if (_delta[i] == 0 && _deg[i] > max_deg) {
   1.299 +            node = i;
   1.300 +            max_deg = _deg[i];
   1.301 +          }
   1.302 +        }
   1.303 +        for (int i = 0; i != start_node; i++) {
   1.304 +          if (_delta[i] == 0 && _deg[i] > max_deg) {
   1.305 +            node = i;
   1.306 +            max_deg = _deg[i];
   1.307 +          }
   1.308 +        }
   1.309 +        return node;
   1.310 +      }
   1.311 +
   1.312 +      // Update internal data structures between stages (if necessary)
   1.313 +      void update() {}
   1.314 +
   1.315 +    }; //class DegreeBasedSelectionRule
   1.316 +
   1.317 +
   1.318 +    // Implementation of the PENALTY_BASED node selection rule.
   1.319 +    class PenaltyBasedSelectionRule
   1.320 +    {
   1.321 +    private:
   1.322 +
   1.323 +      // References to the algorithm instance
   1.324 +      const BoolVector &_clique;
   1.325 +      const IntVector  &_delta;
   1.326 +      const BoolVector &_tabu;
   1.327 +      Random &_rnd;
   1.328 +
   1.329 +      // Pivot rule data
   1.330 +      int _n;
   1.331 +      IntVector _penalty;
   1.332 +
   1.333 +    public:
   1.334 +
   1.335 +      // Constructor
   1.336 +      PenaltyBasedSelectionRule(GrossoLocatelliPullanMc &mc) :
   1.337 +        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
   1.338 +        _rnd(mc._rnd), _n(mc._n), _penalty(_n, 0)
   1.339 +      {}
   1.340 +
   1.341 +      // Return a node index for a feasible add move or -1 if no one exists
   1.342 +      int nextFeasibleAddNode() const {
   1.343 +        int start_node = _rnd[_n];
   1.344 +        int node = -1, min_p = std::numeric_limits<int>::max();
   1.345 +        for (int i = start_node; i != _n; i++) {
   1.346 +          if (_delta[i] == 0 && !_tabu[i] && _penalty[i] < min_p) {
   1.347 +            node = i;
   1.348 +            min_p = _penalty[i];
   1.349 +          }
   1.350 +        }
   1.351 +        for (int i = 0; i != start_node; i++) {
   1.352 +          if (_delta[i] == 0 && !_tabu[i] && _penalty[i] < min_p) {
   1.353 +            node = i;
   1.354 +            min_p = _penalty[i];
   1.355 +          }
   1.356 +        }
   1.357 +        return node;
   1.358 +      }
   1.359 +
   1.360 +      // Return a node index for a feasible swap move or -1 if no one exists
   1.361 +      int nextFeasibleSwapNode() const {
   1.362 +        int start_node = _rnd[_n];
   1.363 +        int node = -1, min_p = std::numeric_limits<int>::max();
   1.364 +        for (int i = start_node; i != _n; i++) {
   1.365 +          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
   1.366 +              _penalty[i] < min_p) {
   1.367 +            node = i;
   1.368 +            min_p = _penalty[i];
   1.369 +          }
   1.370 +        }
   1.371 +        for (int i = 0; i != start_node; i++) {
   1.372 +          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
   1.373 +              _penalty[i] < min_p) {
   1.374 +            node = i;
   1.375 +            min_p = _penalty[i];
   1.376 +          }
   1.377 +        }
   1.378 +        return node;
   1.379 +      }
   1.380 +
   1.381 +      // Return a node index for an add move or -1 if no one exists
   1.382 +      int nextAddNode() const {
   1.383 +        int start_node = _rnd[_n];
   1.384 +        int node = -1, min_p = std::numeric_limits<int>::max();
   1.385 +        for (int i = start_node; i != _n; i++) {
   1.386 +          if (_delta[i] == 0 && _penalty[i] < min_p) {
   1.387 +            node = i;
   1.388 +            min_p = _penalty[i];
   1.389 +          }
   1.390 +        }
   1.391 +        for (int i = 0; i != start_node; i++) {
   1.392 +          if (_delta[i] == 0 && _penalty[i] < min_p) {
   1.393 +            node = i;
   1.394 +            min_p = _penalty[i];
   1.395 +          }
   1.396 +        }
   1.397 +        return node;
   1.398 +      }
   1.399 +
   1.400 +      // Update internal data structures between stages (if necessary)
   1.401 +      void update() {}
   1.402 +
   1.403 +    }; //class PenaltyBasedSelectionRule
   1.404 +
   1.405 +  public:
   1.406 +
   1.407 +    /// \brief Constructor.
   1.408 +    ///
   1.409 +    /// Constructor.
   1.410 +    /// The global \ref rnd "random number generator instance" is used
   1.411 +    /// during the algorithm.
   1.412 +    ///
   1.413 +    /// \param graph The undirected graph the algorithm runs on.
   1.414 +    GrossoLocatelliPullanMc(const GR& graph) :
   1.415 +      _graph(graph), _id(_graph), _rnd(rnd)
   1.416 +    {
   1.417 +      initOptions();
   1.418 +    }
   1.419 +
   1.420 +    /// \brief Constructor with random seed.
   1.421 +    ///
   1.422 +    /// Constructor with random seed.
   1.423 +    ///
   1.424 +    /// \param graph The undirected graph the algorithm runs on.
   1.425 +    /// \param seed Seed value for the internal random number generator
   1.426 +    /// that is used during the algorithm.
   1.427 +    GrossoLocatelliPullanMc(const GR& graph, int seed) :
   1.428 +      _graph(graph), _id(_graph), _rnd(seed)
   1.429 +    {
   1.430 +      initOptions();
   1.431 +    }
   1.432 +
   1.433 +    /// \brief Constructor with random number generator.
   1.434 +    ///
   1.435 +    /// Constructor with random number generator.
   1.436 +    ///
   1.437 +    /// \param graph The undirected graph the algorithm runs on.
   1.438 +    /// \param random A random number generator that is used during the
   1.439 +    /// algorithm.
   1.440 +    GrossoLocatelliPullanMc(const GR& graph, const Random& random) :
   1.441 +      _graph(graph), _id(_graph), _rnd(random)
   1.442 +    {
   1.443 +      initOptions();
   1.444 +    }
   1.445 +
   1.446 +    /// \name Execution Control
   1.447 +    /// The \ref run() function can be used to execute the algorithm.\n
   1.448 +    /// The functions \ref iterationLimit(int), \ref stepLimit(int), and
   1.449 +    /// \ref sizeLimit(int) can be used to specify various limits for the
   1.450 +    /// search process.
   1.451 +
   1.452 +    /// @{
   1.453 +
   1.454 +    /// \brief Sets the maximum number of iterations.
   1.455 +    ///
   1.456 +    /// This function sets the maximum number of iterations.
   1.457 +    /// Each iteration of the algorithm finds a maximal clique (but not
   1.458 +    /// necessarily the largest one) by performing several search steps
   1.459 +    /// (node selections).
   1.460 +    ///
   1.461 +    /// This limit controls the running time and the success of the
   1.462 +    /// algorithm. For larger values, the algorithm runs slower, but it more
   1.463 +    /// likely finds larger cliques. For smaller values, the algorithm is
   1.464 +    /// faster but probably gives worse results.
   1.465 +    ///
   1.466 +    /// The default value is \c 1000.
   1.467 +    /// \c -1 means that number of iterations is not limited.
   1.468 +    ///
   1.469 +    /// \warning You should specify a reasonable limit for the number of
   1.470 +    /// iterations and/or the number of search steps.
   1.471 +    ///
   1.472 +    /// \return <tt>(*this)</tt>
   1.473 +    ///
   1.474 +    /// \sa stepLimit(int)
   1.475 +    /// \sa sizeLimit(int)
   1.476 +    GrossoLocatelliPullanMc& iterationLimit(int limit) {
   1.477 +      _iteration_limit = limit;
   1.478 +      return *this;
   1.479 +    }
   1.480 +
   1.481 +    /// \brief Sets the maximum number of search steps.
   1.482 +    ///
   1.483 +    /// This function sets the maximum number of elementary search steps.
   1.484 +    /// Each iteration of the algorithm finds a maximal clique (but not
   1.485 +    /// necessarily the largest one) by performing several search steps
   1.486 +    /// (node selections).
   1.487 +    ///
   1.488 +    /// This limit controls the running time and the success of the
   1.489 +    /// algorithm. For larger values, the algorithm runs slower, but it more
   1.490 +    /// likely finds larger cliques. For smaller values, the algorithm is
   1.491 +    /// faster but probably gives worse results.
   1.492 +    ///
   1.493 +    /// The default value is \c -1, which means that number of steps
   1.494 +    /// is not limited explicitly. However, the number of iterations is
   1.495 +    /// limited and each iteration performs a finite number of search steps.
   1.496 +    ///
   1.497 +    /// \warning You should specify a reasonable limit for the number of
   1.498 +    /// iterations and/or the number of search steps.
   1.499 +    ///
   1.500 +    /// \return <tt>(*this)</tt>
   1.501 +    ///
   1.502 +    /// \sa iterationLimit(int)
   1.503 +    /// \sa sizeLimit(int)
   1.504 +    GrossoLocatelliPullanMc& stepLimit(int limit) {
   1.505 +      _step_limit = limit;
   1.506 +      return *this;
   1.507 +    }
   1.508 +
   1.509 +    /// \brief Sets the desired clique size.
   1.510 +    ///
   1.511 +    /// This function sets the desired clique size that serves as a search
   1.512 +    /// limit. If a clique of this size (or a larger one) is found, then the
   1.513 +    /// algorithm terminates.
   1.514 +    ///
   1.515 +    /// This function is especially useful if you know an exact upper bound
   1.516 +    /// for the size of the cliques in the graph or if any clique above
   1.517 +    /// a certain size limit is sufficient for your application.
   1.518 +    ///
   1.519 +    /// The default value is \c -1, which means that the size limit is set to
   1.520 +    /// the number of nodes in the graph.
   1.521 +    ///
   1.522 +    /// \return <tt>(*this)</tt>
   1.523 +    ///
   1.524 +    /// \sa iterationLimit(int)
   1.525 +    /// \sa stepLimit(int)
   1.526 +    GrossoLocatelliPullanMc& sizeLimit(int limit) {
   1.527 +      _size_limit = limit;
   1.528 +      return *this;
   1.529 +    }
   1.530 +
   1.531 +    /// \brief The maximum number of iterations.
   1.532 +    ///
   1.533 +    /// This function gives back the maximum number of iterations.
   1.534 +    /// \c -1 means that no limit is specified.
   1.535 +    ///
   1.536 +    /// \sa iterationLimit(int)
   1.537 +    int iterationLimit() const {
   1.538 +      return _iteration_limit;
   1.539 +    }
   1.540 +
   1.541 +    /// \brief The maximum number of search steps.
   1.542 +    ///
   1.543 +    /// This function gives back the maximum number of search steps.
   1.544 +    /// \c -1 means that no limit is specified.
   1.545 +    ///
   1.546 +    /// \sa stepLimit(int)
   1.547 +    int stepLimit() const {
   1.548 +      return _step_limit;
   1.549 +    }
   1.550 +
   1.551 +    /// \brief The desired clique size.
   1.552 +    ///
   1.553 +    /// This function gives back the desired clique size that serves as a
   1.554 +    /// search limit. \c -1 means that this limit is set to the number of
   1.555 +    /// nodes in the graph.
   1.556 +    ///
   1.557 +    /// \sa sizeLimit(int)
   1.558 +    int sizeLimit() const {
   1.559 +      return _size_limit;
   1.560 +    }
   1.561 +
   1.562 +    /// \brief Runs the algorithm.
   1.563 +    ///
   1.564 +    /// This function runs the algorithm. If one of the specified limits
   1.565 +    /// is reached, the search process terminates.
   1.566 +    ///
   1.567 +    /// \param rule The node selection rule. For more information, see
   1.568 +    /// \ref SelectionRule.
   1.569 +    ///
   1.570 +    /// \return The termination cause of the search. For more information,
   1.571 +    /// see \ref TerminationCause.
   1.572 +    TerminationCause run(SelectionRule rule = PENALTY_BASED)
   1.573 +    {
   1.574 +      init();
   1.575 +      switch (rule) {
   1.576 +        case RANDOM:
   1.577 +          return start<RandomSelectionRule>();
   1.578 +        case DEGREE_BASED:
   1.579 +          return start<DegreeBasedSelectionRule>();
   1.580 +        default:
   1.581 +          return start<PenaltyBasedSelectionRule>();
   1.582 +      }
   1.583 +    }
   1.584 +
   1.585 +    /// @}
   1.586 +
   1.587 +    /// \name Query Functions
   1.588 +    /// The results of the algorithm can be obtained using these functions.\n
   1.589 +    /// The run() function must be called before using them.
   1.590 +
   1.591 +    /// @{
   1.592 +
   1.593 +    /// \brief The size of the found clique
   1.594 +    ///
   1.595 +    /// This function returns the size of the found clique.
   1.596 +    ///
   1.597 +    /// \pre run() must be called before using this function.
   1.598 +    int cliqueSize() const {
   1.599 +      return _best_size;
   1.600 +    }
   1.601 +
   1.602 +    /// \brief Gives back the found clique in a \c bool node map
   1.603 +    ///
   1.604 +    /// This function gives back the characteristic vector of the found
   1.605 +    /// clique in the given node map.
   1.606 +    /// It must be a \ref concepts::WriteMap "writable" node map with
   1.607 +    /// \c bool (or convertible) value type.
   1.608 +    ///
   1.609 +    /// \pre run() must be called before using this function.
   1.610 +    template <typename CliqueMap>
   1.611 +    void cliqueMap(CliqueMap &map) const {
   1.612 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.613 +        map[n] = static_cast<bool>(_best_clique[_id[n]]);
   1.614 +      }
   1.615 +    }
   1.616 +
   1.617 +    /// \brief Iterator to list the nodes of the found clique
   1.618 +    ///
   1.619 +    /// This iterator class lists the nodes of the found clique.
   1.620 +    /// Before using it, you must allocate a GrossoLocatelliPullanMc instance
   1.621 +    /// and call its \ref GrossoLocatelliPullanMc::run() "run()" method.
   1.622 +    ///
   1.623 +    /// The following example prints out the IDs of the nodes in the found
   1.624 +    /// clique.
   1.625 +    /// \code
   1.626 +    ///   GrossoLocatelliPullanMc<Graph> mc(g);
   1.627 +    ///   mc.run();
   1.628 +    ///   for (GrossoLocatelliPullanMc<Graph>::CliqueNodeIt n(mc);
   1.629 +    ///        n != INVALID; ++n)
   1.630 +    ///   {
   1.631 +    ///     std::cout << g.id(n) << std::endl;
   1.632 +    ///   }
   1.633 +    /// \endcode
   1.634 +    class CliqueNodeIt
   1.635 +    {
   1.636 +    private:
   1.637 +      NodeIt _it;
   1.638 +      BoolNodeMap _map;
   1.639 +
   1.640 +    public:
   1.641 +
   1.642 +      /// Constructor
   1.643 +
   1.644 +      /// Constructor.
   1.645 +      /// \param mc The algorithm instance.
   1.646 +      CliqueNodeIt(const GrossoLocatelliPullanMc &mc)
   1.647 +       : _map(mc._graph)
   1.648 +      {
   1.649 +        mc.cliqueMap(_map);
   1.650 +        for (_it = NodeIt(mc._graph); _it != INVALID && !_map[_it]; ++_it) ;
   1.651 +      }
   1.652 +
   1.653 +      /// Conversion to \c Node
   1.654 +      operator Node() const { return _it; }
   1.655 +
   1.656 +      bool operator==(Invalid) const { return _it == INVALID; }
   1.657 +      bool operator!=(Invalid) const { return _it != INVALID; }
   1.658 +
   1.659 +      /// Next node
   1.660 +      CliqueNodeIt &operator++() {
   1.661 +        for (++_it; _it != INVALID && !_map[_it]; ++_it) ;
   1.662 +        return *this;
   1.663 +      }
   1.664 +
   1.665 +      /// Postfix incrementation
   1.666 +
   1.667 +      /// Postfix incrementation.
   1.668 +      ///
   1.669 +      /// \warning This incrementation returns a \c Node, not a
   1.670 +      /// \c CliqueNodeIt as one may expect.
   1.671 +      typename GR::Node operator++(int) {
   1.672 +        Node n=*this;
   1.673 +        ++(*this);
   1.674 +        return n;
   1.675 +      }
   1.676 +
   1.677 +    };
   1.678 +
   1.679 +    /// @}
   1.680 +
   1.681 +  private:
   1.682 +
   1.683 +    // Initialize search options and limits
   1.684 +    void initOptions() {
   1.685 +      // Search options
   1.686 +      _delta_based_restart = true;
   1.687 +      _restart_delta_limit = 4;
   1.688 +
   1.689 +      // Search limits
   1.690 +      _iteration_limit = 1000;
   1.691 +      _step_limit = -1;             // this is disabled by default
   1.692 +      _size_limit = -1;             // this is disabled by default
   1.693 +    }
   1.694 +
   1.695 +    // Adds a node to the current clique
   1.696 +    void addCliqueNode(int u) {
   1.697 +      if (_clique[u]) return;
   1.698 +      _clique[u] = true;
   1.699 +      _size++;
   1.700 +      BoolVector &row = _gr[u];
   1.701 +      for (int i = 0; i != _n; i++) {
   1.702 +        if (!row[i]) _delta[i]++;
   1.703 +      }
   1.704 +    }
   1.705 +
   1.706 +    // Removes a node from the current clique
   1.707 +    void delCliqueNode(int u) {
   1.708 +      if (!_clique[u]) return;
   1.709 +      _clique[u] = false;
   1.710 +      _size--;
   1.711 +      BoolVector &row = _gr[u];
   1.712 +      for (int i = 0; i != _n; i++) {
   1.713 +        if (!row[i]) _delta[i]--;
   1.714 +      }
   1.715 +    }
   1.716 +
   1.717 +    // Initialize data structures
   1.718 +    void init() {
   1.719 +      _n = countNodes(_graph);
   1.720 +      int ui = 0;
   1.721 +      for (NodeIt u(_graph); u != INVALID; ++u) {
   1.722 +        _id[u] = ui++;
   1.723 +      }
   1.724 +      _gr.clear();
   1.725 +      _gr.resize(_n, BoolVector(_n, false));
   1.726 +      ui = 0;
   1.727 +      for (NodeIt u(_graph); u != INVALID; ++u) {
   1.728 +        for (IncEdgeIt e(_graph, u); e != INVALID; ++e) {
   1.729 +          int vi = _id[_graph.runningNode(e)];
   1.730 +          _gr[ui][vi] = true;
   1.731 +          _gr[vi][ui] = true;
   1.732 +        }
   1.733 +        ++ui;
   1.734 +      }
   1.735 +
   1.736 +      _clique.clear();
   1.737 +      _clique.resize(_n, false);
   1.738 +      _size = 0;
   1.739 +      _best_clique.clear();
   1.740 +      _best_clique.resize(_n, false);
   1.741 +      _best_size = 0;
   1.742 +      _delta.clear();
   1.743 +      _delta.resize(_n, 0);
   1.744 +      _tabu.clear();
   1.745 +      _tabu.resize(_n, false);
   1.746 +    }
   1.747 +
   1.748 +    // Executes the algorithm
   1.749 +    template <typename SelectionRuleImpl>
   1.750 +    TerminationCause start() {
   1.751 +      if (_n == 0) return SIZE_LIMIT;
   1.752 +      if (_n == 1) {
   1.753 +        _best_clique[0] = true;
   1.754 +        _best_size = 1;
   1.755 +        return SIZE_LIMIT;
   1.756 +      }
   1.757 +
   1.758 +      // Iterated local search algorithm
   1.759 +      const int max_size = _size_limit >= 0 ? _size_limit : _n;
   1.760 +      const int max_restart = _iteration_limit >= 0 ?
   1.761 +        _iteration_limit : std::numeric_limits<int>::max();
   1.762 +      const int max_select = _step_limit >= 0 ?
   1.763 +        _step_limit : std::numeric_limits<int>::max();
   1.764 +
   1.765 +      SelectionRuleImpl sel_method(*this);
   1.766 +      int select = 0, restart = 0;
   1.767 +      IntVector restart_nodes;
   1.768 +      while (select < max_select && restart < max_restart) {
   1.769 +
   1.770 +        // Perturbation/restart
   1.771 +        restart++;
   1.772 +        if (_delta_based_restart) {
   1.773 +          restart_nodes.clear();
   1.774 +          for (int i = 0; i != _n; i++) {
   1.775 +            if (_delta[i] >= _restart_delta_limit)
   1.776 +              restart_nodes.push_back(i);
   1.777 +          }
   1.778 +        }
   1.779 +        int rs_node = -1;
   1.780 +        if (restart_nodes.size() > 0) {
   1.781 +          rs_node = restart_nodes[_rnd[restart_nodes.size()]];
   1.782 +        } else {
   1.783 +          rs_node = _rnd[_n];
   1.784 +        }
   1.785 +        BoolVector &row = _gr[rs_node];
   1.786 +        for (int i = 0; i != _n; i++) {
   1.787 +          if (_clique[i] && !row[i]) delCliqueNode(i);
   1.788 +        }
   1.789 +        addCliqueNode(rs_node);
   1.790 +
   1.791 +        // Local search
   1.792 +        _tabu.clear();
   1.793 +        _tabu.resize(_n, false);
   1.794 +        bool tabu_empty = true;
   1.795 +        int max_swap = _size;
   1.796 +        while (select < max_select) {
   1.797 +          select++;
   1.798 +          int u;
   1.799 +          if ((u = sel_method.nextFeasibleAddNode()) != -1) {
   1.800 +            // Feasible add move
   1.801 +            addCliqueNode(u);
   1.802 +            if (tabu_empty) max_swap = _size;
   1.803 +          }
   1.804 +          else if ((u = sel_method.nextFeasibleSwapNode()) != -1) {
   1.805 +            // Feasible swap move
   1.806 +            int v = -1;
   1.807 +            BoolVector &row = _gr[u];
   1.808 +            for (int i = 0; i != _n; i++) {
   1.809 +              if (_clique[i] && !row[i]) {
   1.810 +                v = i;
   1.811 +                break;
   1.812 +              }
   1.813 +            }
   1.814 +            addCliqueNode(u);
   1.815 +            delCliqueNode(v);
   1.816 +            _tabu[v] = true;
   1.817 +            tabu_empty = false;
   1.818 +            if (--max_swap <= 0) break;
   1.819 +          }
   1.820 +          else if ((u = sel_method.nextAddNode()) != -1) {
   1.821 +            // Non-feasible add move
   1.822 +            addCliqueNode(u);
   1.823 +          }
   1.824 +          else break;
   1.825 +        }
   1.826 +        if (_size > _best_size) {
   1.827 +          _best_clique = _clique;
   1.828 +          _best_size = _size;
   1.829 +          if (_best_size >= max_size) return SIZE_LIMIT;
   1.830 +        }
   1.831 +        sel_method.update();
   1.832 +      }
   1.833 +
   1.834 +      return (restart >= max_restart ? ITERATION_LIMIT : STEP_LIMIT);
   1.835 +    }
   1.836 +
   1.837 +  }; //class GrossoLocatelliPullanMc
   1.838 +
   1.839 +  ///@}
   1.840 +
   1.841 +} //namespace lemon
   1.842 +
   1.843 +#endif //LEMON_GROSSO_LOCATELLI_PULLAN_MC_H