kpeter@904: /* -*- mode: C++; indent-tabs-mode: nil; -*-
kpeter@904:  *
kpeter@904:  * This file is a part of LEMON, a generic C++ optimization library.
kpeter@904:  *
kpeter@904:  * Copyright (C) 2003-2010
kpeter@904:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
kpeter@904:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
kpeter@904:  *
kpeter@904:  * Permission to use, modify and distribute this software is granted
kpeter@904:  * provided that this copyright notice appears in all copies. For
kpeter@904:  * precise terms see the accompanying LICENSE file.
kpeter@904:  *
kpeter@904:  * This software is provided "AS IS" with no warranty of any kind,
kpeter@904:  * express or implied, and with no claim as to its suitability for any
kpeter@904:  * purpose.
kpeter@904:  *
kpeter@904:  */
kpeter@904: 
kpeter@904: #ifndef LEMON_GROSSO_LOCATELLI_PULLAN_MC_H
kpeter@904: #define LEMON_GROSSO_LOCATELLI_PULLAN_MC_H
kpeter@904: 
kpeter@904: /// \ingroup approx_algs
kpeter@904: ///
kpeter@904: /// \file
kpeter@904: /// \brief The iterated local search algorithm of Grosso, Locatelli, and Pullan
kpeter@904: /// for the maximum clique problem
kpeter@904: 
kpeter@904: #include <vector>
kpeter@904: #include <limits>
kpeter@904: #include <lemon/core.h>
kpeter@904: #include <lemon/random.h>
kpeter@904: 
kpeter@904: namespace lemon {
kpeter@904: 
kpeter@904:   /// \addtogroup approx_algs
kpeter@904:   /// @{
kpeter@904: 
kpeter@904:   /// \brief Implementation of the iterated local search algorithm of Grosso,
kpeter@904:   /// Locatelli, and Pullan for the maximum clique problem
kpeter@904:   ///
kpeter@904:   /// \ref GrossoLocatelliPullanMc implements the iterated local search
kpeter@904:   /// algorithm of Grosso, Locatelli, and Pullan for solving the \e maximum
kpeter@904:   /// \e clique \e problem \ref grosso08maxclique.
kpeter@904:   /// It is to find the largest complete subgraph (\e clique) in an
kpeter@904:   /// undirected graph, i.e., the largest set of nodes where each
kpeter@904:   /// pair of nodes is connected.
kpeter@904:   ///
kpeter@904:   /// This class provides a simple but highly efficient and robust heuristic
kpeter@904:   /// method that quickly finds a large clique, but not necessarily the
kpeter@904:   /// largest one.
kpeter@904:   ///
kpeter@904:   /// \tparam GR The undirected graph type the algorithm runs on.
kpeter@904:   ///
kpeter@904:   /// \note %GrossoLocatelliPullanMc provides three different node selection
kpeter@904:   /// rules, from which the most powerful one is used by default.
kpeter@904:   /// For more information, see \ref SelectionRule.
kpeter@904:   template <typename GR>
kpeter@904:   class GrossoLocatelliPullanMc
kpeter@904:   {
kpeter@904:   public:
kpeter@904: 
kpeter@904:     /// \brief Constants for specifying the node selection rule.
kpeter@904:     ///
kpeter@904:     /// Enum type containing constants for specifying the node selection rule
kpeter@904:     /// for the \ref run() function.
kpeter@904:     ///
kpeter@904:     /// During the algorithm, nodes are selected for addition to the current
kpeter@904:     /// clique according to the applied rule.
kpeter@904:     /// In general, the PENALTY_BASED rule turned out to be the most powerful
kpeter@904:     /// and the most robust, thus it is the default option.
kpeter@904:     /// However, another selection rule can be specified using the \ref run()
kpeter@904:     /// function with the proper parameter.
kpeter@904:     enum SelectionRule {
kpeter@904: 
kpeter@904:       /// A node is selected randomly without any evaluation at each step.
kpeter@904:       RANDOM,
kpeter@904: 
kpeter@904:       /// A node of maximum degree is selected randomly at each step.
kpeter@904:       DEGREE_BASED,
kpeter@904: 
kpeter@904:       /// A node of minimum penalty is selected randomly at each step.
kpeter@904:       /// The node penalties are updated adaptively after each stage of the
kpeter@904:       /// search process.
kpeter@904:       PENALTY_BASED
kpeter@904:     };
kpeter@904: 
kpeter@904:   private:
kpeter@904: 
kpeter@904:     TEMPLATE_GRAPH_TYPEDEFS(GR);
kpeter@904: 
kpeter@904:     typedef std::vector<int> IntVector;
kpeter@904:     typedef std::vector<char> BoolVector;
kpeter@904:     typedef std::vector<BoolVector> BoolMatrix;
kpeter@904:     // Note: vector<char> is used instead of vector<bool> for efficiency reasons
kpeter@904: 
kpeter@904:     const GR &_graph;
kpeter@904:     IntNodeMap _id;
kpeter@904: 
kpeter@904:     // Internal matrix representation of the graph
kpeter@904:     BoolMatrix _gr;
kpeter@904:     int _n;
kpeter@904: 
kpeter@904:     // The current clique
kpeter@904:     BoolVector _clique;
kpeter@904:     int _size;
kpeter@904: 
kpeter@904:     // The best clique found so far
kpeter@904:     BoolVector _best_clique;
kpeter@904:     int _best_size;
kpeter@904: 
kpeter@904:     // The "distances" of the nodes from the current clique.
kpeter@904:     // _delta[u] is the number of nodes in the clique that are
kpeter@904:     // not connected with u.
kpeter@904:     IntVector _delta;
kpeter@904: 
kpeter@904:     // The current tabu set
kpeter@904:     BoolVector _tabu;
kpeter@904: 
kpeter@904:     // Random number generator
kpeter@904:     Random _rnd;
kpeter@904: 
kpeter@904:   private:
kpeter@904: 
kpeter@904:     // Implementation of the RANDOM node selection rule.
kpeter@904:     class RandomSelectionRule
kpeter@904:     {
kpeter@904:     private:
kpeter@904: 
kpeter@904:       // References to the algorithm instance
kpeter@904:       const BoolVector &_clique;
kpeter@904:       const IntVector  &_delta;
kpeter@904:       const BoolVector &_tabu;
kpeter@904:       Random &_rnd;
kpeter@904: 
kpeter@904:       // Pivot rule data
kpeter@904:       int _n;
kpeter@904: 
kpeter@904:     public:
kpeter@904: 
kpeter@904:       // Constructor
kpeter@904:       RandomSelectionRule(GrossoLocatelliPullanMc &mc) :
kpeter@904:         _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
kpeter@904:         _rnd(mc._rnd), _n(mc._n)
kpeter@904:       {}
kpeter@904: 
kpeter@904:       // Return a node index for a feasible add move or -1 if no one exists
kpeter@904:       int nextFeasibleAddNode() const {
kpeter@904:         int start_node = _rnd[_n];
kpeter@904:         for (int i = start_node; i != _n; i++) {
kpeter@904:           if (_delta[i] == 0 && !_tabu[i]) return i;
kpeter@904:         }
kpeter@904:         for (int i = 0; i != start_node; i++) {
kpeter@904:           if (_delta[i] == 0 && !_tabu[i]) return i;
kpeter@904:         }
kpeter@904:         return -1;
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Return a node index for a feasible swap move or -1 if no one exists
kpeter@904:       int nextFeasibleSwapNode() const {
kpeter@904:         int start_node = _rnd[_n];
kpeter@904:         for (int i = start_node; i != _n; i++) {
kpeter@904:           if (!_clique[i] && _delta[i] == 1 && !_tabu[i]) return i;
kpeter@904:         }
kpeter@904:         for (int i = 0; i != start_node; i++) {
kpeter@904:           if (!_clique[i] && _delta[i] == 1 && !_tabu[i]) return i;
kpeter@904:         }
kpeter@904:         return -1;
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Return a node index for an add move or -1 if no one exists
kpeter@904:       int nextAddNode() const {
kpeter@904:         int start_node = _rnd[_n];
kpeter@904:         for (int i = start_node; i != _n; i++) {
kpeter@904:           if (_delta[i] == 0) return i;
kpeter@904:         }
kpeter@904:         for (int i = 0; i != start_node; i++) {
kpeter@904:           if (_delta[i] == 0) return i;
kpeter@904:         }
kpeter@904:         return -1;
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Update internal data structures between stages (if necessary)
kpeter@904:       void update() {}
kpeter@904: 
kpeter@904:     }; //class RandomSelectionRule
kpeter@904: 
kpeter@904: 
kpeter@904:     // Implementation of the DEGREE_BASED node selection rule.
kpeter@904:     class DegreeBasedSelectionRule
kpeter@904:     {
kpeter@904:     private:
kpeter@904: 
kpeter@904:       // References to the algorithm instance
kpeter@904:       const BoolVector &_clique;
kpeter@904:       const IntVector  &_delta;
kpeter@904:       const BoolVector &_tabu;
kpeter@904:       Random &_rnd;
kpeter@904: 
kpeter@904:       // Pivot rule data
kpeter@904:       int _n;
kpeter@904:       IntVector _deg;
kpeter@904: 
kpeter@904:     public:
kpeter@904: 
kpeter@904:       // Constructor
kpeter@904:       DegreeBasedSelectionRule(GrossoLocatelliPullanMc &mc) :
kpeter@904:         _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
kpeter@904:         _rnd(mc._rnd), _n(mc._n), _deg(_n)
kpeter@904:       {
kpeter@904:         for (int i = 0; i != _n; i++) {
kpeter@904:           int d = 0;
kpeter@904:           BoolVector &row = mc._gr[i];
kpeter@904:           for (int j = 0; j != _n; j++) {
kpeter@904:             if (row[j]) d++;
kpeter@904:           }
kpeter@904:           _deg[i] = d;
kpeter@904:         }
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Return a node index for a feasible add move or -1 if no one exists
kpeter@904:       int nextFeasibleAddNode() const {
kpeter@904:         int start_node = _rnd[_n];
kpeter@904:         int node = -1, max_deg = -1;
kpeter@904:         for (int i = start_node; i != _n; i++) {
kpeter@904:           if (_delta[i] == 0 && !_tabu[i] && _deg[i] > max_deg) {
kpeter@904:             node = i;
kpeter@904:             max_deg = _deg[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         for (int i = 0; i != start_node; i++) {
kpeter@904:           if (_delta[i] == 0 && !_tabu[i] && _deg[i] > max_deg) {
kpeter@904:             node = i;
kpeter@904:             max_deg = _deg[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         return node;
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Return a node index for a feasible swap move or -1 if no one exists
kpeter@904:       int nextFeasibleSwapNode() const {
kpeter@904:         int start_node = _rnd[_n];
kpeter@904:         int node = -1, max_deg = -1;
kpeter@904:         for (int i = start_node; i != _n; i++) {
kpeter@904:           if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
kpeter@904:               _deg[i] > max_deg) {
kpeter@904:             node = i;
kpeter@904:             max_deg = _deg[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         for (int i = 0; i != start_node; i++) {
kpeter@904:           if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
kpeter@904:               _deg[i] > max_deg) {
kpeter@904:             node = i;
kpeter@904:             max_deg = _deg[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         return node;
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Return a node index for an add move or -1 if no one exists
kpeter@904:       int nextAddNode() const {
kpeter@904:         int start_node = _rnd[_n];
kpeter@904:         int node = -1, max_deg = -1;
kpeter@904:         for (int i = start_node; i != _n; i++) {
kpeter@904:           if (_delta[i] == 0 && _deg[i] > max_deg) {
kpeter@904:             node = i;
kpeter@904:             max_deg = _deg[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         for (int i = 0; i != start_node; i++) {
kpeter@904:           if (_delta[i] == 0 && _deg[i] > max_deg) {
kpeter@904:             node = i;
kpeter@904:             max_deg = _deg[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         return node;
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Update internal data structures between stages (if necessary)
kpeter@904:       void update() {}
kpeter@904: 
kpeter@904:     }; //class DegreeBasedSelectionRule
kpeter@904: 
kpeter@904: 
kpeter@904:     // Implementation of the PENALTY_BASED node selection rule.
kpeter@904:     class PenaltyBasedSelectionRule
kpeter@904:     {
kpeter@904:     private:
kpeter@904: 
kpeter@904:       // References to the algorithm instance
kpeter@904:       const BoolVector &_clique;
kpeter@904:       const IntVector  &_delta;
kpeter@904:       const BoolVector &_tabu;
kpeter@904:       Random &_rnd;
kpeter@904: 
kpeter@904:       // Pivot rule data
kpeter@904:       int _n;
kpeter@904:       IntVector _penalty;
kpeter@904: 
kpeter@904:     public:
kpeter@904: 
kpeter@904:       // Constructor
kpeter@904:       PenaltyBasedSelectionRule(GrossoLocatelliPullanMc &mc) :
kpeter@904:         _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
kpeter@904:         _rnd(mc._rnd), _n(mc._n), _penalty(_n, 0)
kpeter@904:       {}
kpeter@904: 
kpeter@904:       // Return a node index for a feasible add move or -1 if no one exists
kpeter@904:       int nextFeasibleAddNode() const {
kpeter@904:         int start_node = _rnd[_n];
kpeter@904:         int node = -1, min_p = std::numeric_limits<int>::max();
kpeter@904:         for (int i = start_node; i != _n; i++) {
kpeter@904:           if (_delta[i] == 0 && !_tabu[i] && _penalty[i] < min_p) {
kpeter@904:             node = i;
kpeter@904:             min_p = _penalty[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         for (int i = 0; i != start_node; i++) {
kpeter@904:           if (_delta[i] == 0 && !_tabu[i] && _penalty[i] < min_p) {
kpeter@904:             node = i;
kpeter@904:             min_p = _penalty[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         return node;
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Return a node index for a feasible swap move or -1 if no one exists
kpeter@904:       int nextFeasibleSwapNode() const {
kpeter@904:         int start_node = _rnd[_n];
kpeter@904:         int node = -1, min_p = std::numeric_limits<int>::max();
kpeter@904:         for (int i = start_node; i != _n; i++) {
kpeter@904:           if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
kpeter@904:               _penalty[i] < min_p) {
kpeter@904:             node = i;
kpeter@904:             min_p = _penalty[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         for (int i = 0; i != start_node; i++) {
kpeter@904:           if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
kpeter@904:               _penalty[i] < min_p) {
kpeter@904:             node = i;
kpeter@904:             min_p = _penalty[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         return node;
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Return a node index for an add move or -1 if no one exists
kpeter@904:       int nextAddNode() const {
kpeter@904:         int start_node = _rnd[_n];
kpeter@904:         int node = -1, min_p = std::numeric_limits<int>::max();
kpeter@904:         for (int i = start_node; i != _n; i++) {
kpeter@904:           if (_delta[i] == 0 && _penalty[i] < min_p) {
kpeter@904:             node = i;
kpeter@904:             min_p = _penalty[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         for (int i = 0; i != start_node; i++) {
kpeter@904:           if (_delta[i] == 0 && _penalty[i] < min_p) {
kpeter@904:             node = i;
kpeter@904:             min_p = _penalty[i];
kpeter@904:           }
kpeter@904:         }
kpeter@904:         return node;
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Update internal data structures between stages (if necessary)
kpeter@904:       void update() {}
kpeter@904: 
kpeter@904:     }; //class PenaltyBasedSelectionRule
kpeter@904: 
kpeter@904:   public:
kpeter@904: 
kpeter@904:     /// \brief Constructor.
kpeter@904:     ///
kpeter@904:     /// Constructor.
kpeter@904:     /// The global \ref rnd "random number generator instance" is used
kpeter@904:     /// during the algorithm.
kpeter@904:     ///
kpeter@904:     /// \param graph The undirected graph the algorithm runs on.
kpeter@904:     GrossoLocatelliPullanMc(const GR& graph) :
kpeter@904:       _graph(graph), _id(_graph), _rnd(rnd)
kpeter@904:     {}
kpeter@904: 
kpeter@904:     /// \brief Constructor with random seed.
kpeter@904:     ///
kpeter@904:     /// Constructor with random seed.
kpeter@904:     ///
kpeter@904:     /// \param graph The undirected graph the algorithm runs on.
kpeter@904:     /// \param seed Seed value for the internal random number generator
kpeter@904:     /// that is used during the algorithm.
kpeter@904:     GrossoLocatelliPullanMc(const GR& graph, int seed) :
kpeter@904:       _graph(graph), _id(_graph), _rnd(seed)
kpeter@904:     {}
kpeter@904: 
kpeter@904:     /// \brief Constructor with random number generator.
kpeter@904:     ///
kpeter@904:     /// Constructor with random number generator.
kpeter@904:     ///
kpeter@904:     /// \param graph The undirected graph the algorithm runs on.
kpeter@904:     /// \param random A random number generator that is used during the
kpeter@904:     /// algorithm.
kpeter@904:     GrossoLocatelliPullanMc(const GR& graph, const Random& random) :
kpeter@904:       _graph(graph), _id(_graph), _rnd(random)
kpeter@904:     {}
kpeter@904: 
kpeter@904:     /// \name Execution Control
kpeter@904:     /// @{
kpeter@904: 
kpeter@904:     /// \brief Runs the algorithm.
kpeter@904:     ///
kpeter@904:     /// This function runs the algorithm.
kpeter@904:     ///
kpeter@904:     /// \param step_num The maximum number of node selections (steps)
kpeter@904:     /// during the search process.
kpeter@904:     /// This parameter controls the running time and the success of the
kpeter@904:     /// algorithm. For larger values, the algorithm runs slower but it more
kpeter@904:     /// likely finds larger cliques. For smaller values, the algorithm is
kpeter@904:     /// faster but probably gives worse results.
kpeter@904:     /// \param rule The node selection rule. For more information, see
kpeter@904:     /// \ref SelectionRule.
kpeter@904:     ///
kpeter@904:     /// \return The size of the found clique.
kpeter@904:     int run(int step_num = 100000,
kpeter@904:             SelectionRule rule = PENALTY_BASED)
kpeter@904:     {
kpeter@904:       init();
kpeter@904:       switch (rule) {
kpeter@904:         case RANDOM:
kpeter@904:           return start<RandomSelectionRule>(step_num);
kpeter@904:         case DEGREE_BASED:
kpeter@904:           return start<DegreeBasedSelectionRule>(step_num);
kpeter@904:         case PENALTY_BASED:
kpeter@904:           return start<PenaltyBasedSelectionRule>(step_num);
kpeter@904:       }
kpeter@904:       return 0; // avoid warning
kpeter@904:     }
kpeter@904: 
kpeter@904:     /// @}
kpeter@904: 
kpeter@904:     /// \name Query Functions
kpeter@904:     /// @{
kpeter@904: 
kpeter@904:     /// \brief The size of the found clique
kpeter@904:     ///
kpeter@904:     /// This function returns the size of the found clique.
kpeter@904:     ///
kpeter@904:     /// \pre run() must be called before using this function.
kpeter@904:     int cliqueSize() const {
kpeter@904:       return _best_size;
kpeter@904:     }
kpeter@904: 
kpeter@904:     /// \brief Gives back the found clique in a \c bool node map
kpeter@904:     ///
kpeter@904:     /// This function gives back the characteristic vector of the found
kpeter@904:     /// clique in the given node map.
kpeter@904:     /// It must be a \ref concepts::WriteMap "writable" node map with
kpeter@904:     /// \c bool (or convertible) value type.
kpeter@904:     ///
kpeter@904:     /// \pre run() must be called before using this function.
kpeter@904:     template <typename CliqueMap>
kpeter@904:     void cliqueMap(CliqueMap &map) const {
kpeter@904:       for (NodeIt n(_graph); n != INVALID; ++n) {
kpeter@904:         map[n] = static_cast<bool>(_best_clique[_id[n]]);
kpeter@904:       }
kpeter@904:     }
kpeter@904: 
kpeter@904:     /// \brief Iterator to list the nodes of the found clique
kpeter@904:     ///
kpeter@904:     /// This iterator class lists the nodes of the found clique.
kpeter@904:     /// Before using it, you must allocate a GrossoLocatelliPullanMc instance
kpeter@904:     /// and call its \ref GrossoLocatelliPullanMc::run() "run()" method.
kpeter@904:     ///
kpeter@904:     /// The following example prints out the IDs of the nodes in the found
kpeter@904:     /// clique.
kpeter@904:     /// \code
kpeter@904:     ///   GrossoLocatelliPullanMc<Graph> mc(g);
kpeter@904:     ///   mc.run();
kpeter@904:     ///   for (GrossoLocatelliPullanMc<Graph>::CliqueNodeIt n(mc);
kpeter@904:     ///        n != INVALID; ++n)
kpeter@904:     ///   {
kpeter@904:     ///     std::cout << g.id(n) << std::endl;
kpeter@904:     ///   }
kpeter@904:     /// \endcode
kpeter@904:     class CliqueNodeIt
kpeter@904:     {
kpeter@904:     private:
kpeter@904:       NodeIt _it;
kpeter@904:       BoolNodeMap _map;
kpeter@904: 
kpeter@904:     public:
kpeter@904: 
kpeter@904:       /// Constructor
kpeter@904: 
kpeter@904:       /// Constructor.
kpeter@904:       /// \param mc The algorithm instance.
kpeter@904:       CliqueNodeIt(const GrossoLocatelliPullanMc &mc)
kpeter@904:        : _map(mc._graph)
kpeter@904:       {
kpeter@904:         mc.cliqueMap(_map);
kpeter@904:         for (_it = NodeIt(mc._graph); _it != INVALID && !_map[_it]; ++_it) ;
kpeter@904:       }
kpeter@904: 
kpeter@904:       /// Conversion to \c Node
kpeter@904:       operator Node() const { return _it; }
kpeter@904: 
kpeter@904:       bool operator==(Invalid) const { return _it == INVALID; }
kpeter@904:       bool operator!=(Invalid) const { return _it != INVALID; }
kpeter@904: 
kpeter@904:       /// Next node
kpeter@904:       CliqueNodeIt &operator++() {
kpeter@904:         for (++_it; _it != INVALID && !_map[_it]; ++_it) ;
kpeter@904:         return *this;
kpeter@904:       }
kpeter@904: 
kpeter@904:       /// Postfix incrementation
kpeter@904: 
kpeter@904:       /// Postfix incrementation.
kpeter@904:       ///
kpeter@904:       /// \warning This incrementation returns a \c Node, not a
kpeter@904:       /// \c CliqueNodeIt as one may expect.
kpeter@904:       typename GR::Node operator++(int) {
kpeter@904:         Node n=*this;
kpeter@904:         ++(*this);
kpeter@904:         return n;
kpeter@904:       }
kpeter@904: 
kpeter@904:     };
kpeter@904: 
kpeter@904:     /// @}
kpeter@904: 
kpeter@904:   private:
kpeter@904: 
kpeter@904:     // Adds a node to the current clique
kpeter@904:     void addCliqueNode(int u) {
kpeter@904:       if (_clique[u]) return;
kpeter@904:       _clique[u] = true;
kpeter@904:       _size++;
kpeter@904:       BoolVector &row = _gr[u];
kpeter@904:       for (int i = 0; i != _n; i++) {
kpeter@904:         if (!row[i]) _delta[i]++;
kpeter@904:       }
kpeter@904:     }
kpeter@904: 
kpeter@904:     // Removes a node from the current clique
kpeter@904:     void delCliqueNode(int u) {
kpeter@904:       if (!_clique[u]) return;
kpeter@904:       _clique[u] = false;
kpeter@904:       _size--;
kpeter@904:       BoolVector &row = _gr[u];
kpeter@904:       for (int i = 0; i != _n; i++) {
kpeter@904:         if (!row[i]) _delta[i]--;
kpeter@904:       }
kpeter@904:     }
kpeter@904: 
kpeter@904:     // Initialize data structures
kpeter@904:     void init() {
kpeter@904:       _n = countNodes(_graph);
kpeter@904:       int ui = 0;
kpeter@904:       for (NodeIt u(_graph); u != INVALID; ++u) {
kpeter@904:         _id[u] = ui++;
kpeter@904:       }
kpeter@904:       _gr.clear();
kpeter@904:       _gr.resize(_n, BoolVector(_n, false));
kpeter@904:       ui = 0;
kpeter@904:       for (NodeIt u(_graph); u != INVALID; ++u) {
kpeter@904:         for (IncEdgeIt e(_graph, u); e != INVALID; ++e) {
kpeter@904:           int vi = _id[_graph.runningNode(e)];
kpeter@904:           _gr[ui][vi] = true;
kpeter@904:           _gr[vi][ui] = true;
kpeter@904:         }
kpeter@904:         ++ui;
kpeter@904:       }
kpeter@904: 
kpeter@904:       _clique.clear();
kpeter@904:       _clique.resize(_n, false);
kpeter@904:       _size = 0;
kpeter@904:       _best_clique.clear();
kpeter@904:       _best_clique.resize(_n, false);
kpeter@904:       _best_size = 0;
kpeter@904:       _delta.clear();
kpeter@904:       _delta.resize(_n, 0);
kpeter@904:       _tabu.clear();
kpeter@904:       _tabu.resize(_n, false);
kpeter@904:     }
kpeter@904: 
kpeter@904:     // Executes the algorithm
kpeter@904:     template <typename SelectionRuleImpl>
kpeter@904:     int start(int max_select) {
kpeter@904:       // Options for the restart rule
kpeter@904:       const bool delta_based_restart = true;
kpeter@904:       const int restart_delta_limit = 4;
kpeter@904: 
kpeter@904:       if (_n == 0) return 0;
kpeter@904:       if (_n == 1) {
kpeter@904:         _best_clique[0] = true;
kpeter@904:         _best_size = 1;
kpeter@904:         return _best_size;
kpeter@904:       }
kpeter@904: 
kpeter@904:       // Iterated local search
kpeter@904:       SelectionRuleImpl sel_method(*this);
kpeter@904:       int select = 0;
kpeter@904:       IntVector restart_nodes;
kpeter@904: 
kpeter@904:       while (select < max_select) {
kpeter@904: 
kpeter@904:         // Perturbation/restart
kpeter@904:         if (delta_based_restart) {
kpeter@904:           restart_nodes.clear();
kpeter@904:           for (int i = 0; i != _n; i++) {
kpeter@904:             if (_delta[i] >= restart_delta_limit)
kpeter@904:               restart_nodes.push_back(i);
kpeter@904:           }
kpeter@904:         }
kpeter@904:         int rs_node = -1;
kpeter@904:         if (restart_nodes.size() > 0) {
kpeter@904:           rs_node = restart_nodes[_rnd[restart_nodes.size()]];
kpeter@904:         } else {
kpeter@904:           rs_node = _rnd[_n];
kpeter@904:         }
kpeter@904:         BoolVector &row = _gr[rs_node];
kpeter@904:         for (int i = 0; i != _n; i++) {
kpeter@904:           if (_clique[i] && !row[i]) delCliqueNode(i);
kpeter@904:         }
kpeter@904:         addCliqueNode(rs_node);
kpeter@904: 
kpeter@904:         // Local search
kpeter@904:         _tabu.clear();
kpeter@904:         _tabu.resize(_n, false);
kpeter@904:         bool tabu_empty = true;
kpeter@904:         int max_swap = _size;
kpeter@904:         while (select < max_select) {
kpeter@904:           select++;
kpeter@904:           int u;
kpeter@904:           if ((u = sel_method.nextFeasibleAddNode()) != -1) {
kpeter@904:             // Feasible add move
kpeter@904:             addCliqueNode(u);
kpeter@904:             if (tabu_empty) max_swap = _size;
kpeter@904:           }
kpeter@904:           else if ((u = sel_method.nextFeasibleSwapNode()) != -1) {
kpeter@904:             // Feasible swap move
kpeter@904:             int v = -1;
kpeter@904:             BoolVector &row = _gr[u];
kpeter@904:             for (int i = 0; i != _n; i++) {
kpeter@904:               if (_clique[i] && !row[i]) {
kpeter@904:                 v = i;
kpeter@904:                 break;
kpeter@904:               }
kpeter@904:             }
kpeter@904:             addCliqueNode(u);
kpeter@904:             delCliqueNode(v);
kpeter@904:             _tabu[v] = true;
kpeter@904:             tabu_empty = false;
kpeter@904:             if (--max_swap <= 0) break;
kpeter@904:           }
kpeter@904:           else if ((u = sel_method.nextAddNode()) != -1) {
kpeter@904:             // Non-feasible add move
kpeter@904:             addCliqueNode(u);
kpeter@904:           }
kpeter@904:           else break;
kpeter@904:         }
kpeter@904:         if (_size > _best_size) {
kpeter@904:           _best_clique = _clique;
kpeter@904:           _best_size = _size;
kpeter@904:           if (_best_size == _n) return _best_size;
kpeter@904:         }
kpeter@904:         sel_method.update();
kpeter@904:       }
kpeter@904: 
kpeter@904:       return _best_size;
kpeter@904:     }
kpeter@904: 
kpeter@904:   }; //class GrossoLocatelliPullanMc
kpeter@904: 
kpeter@904:   ///@}
kpeter@904: 
kpeter@904: } //namespace lemon
kpeter@904: 
kpeter@904: #endif //LEMON_GROSSO_LOCATELLI_PULLAN_MC_H