/* -*- mode: C++; indent-tabs-mode: nil; -*-
 *
 * This file is a part of LEMON, a generic C++ optimization library.
 *
 * Copyright (C) 2003-2010
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
 *
 * Permission to use, modify and distribute this software is granted
 * provided that this copyright notice appears in all copies. For
 * precise terms see the accompanying LICENSE file.
 *
 * This software is provided "AS IS" with no warranty of any kind,
 * express or implied, and with no claim as to its suitability for any
 * purpose.
 *
 */

#ifndef LEMON_GROSSO_LOCATELLI_PULLAN_MC_H
#define LEMON_GROSSO_LOCATELLI_PULLAN_MC_H

/// \ingroup approx_algs
///
/// \file
/// \brief The iterated local search algorithm of Grosso, Locatelli, and Pullan
/// for the maximum clique problem

#include <vector>
#include <limits>
#include <lemon/core.h>
#include <lemon/random.h>

namespace lemon {

  /// \addtogroup approx_algs
  /// @{

  /// \brief Implementation of the iterated local search algorithm of Grosso,
  /// Locatelli, and Pullan for the maximum clique problem
  ///
  /// \ref GrossoLocatelliPullanMc implements the iterated local search
  /// algorithm of Grosso, Locatelli, and Pullan for solving the \e maximum
  /// \e clique \e problem \cite grosso08maxclique.
  /// It is to find the largest complete subgraph (\e clique) in an
  /// undirected graph, i.e., the largest set of nodes where each
  /// pair of nodes is connected.
  ///
  /// This class provides a simple but highly efficient and robust heuristic
  /// method that quickly finds a quite large clique, but not necessarily the
  /// largest one.
  /// The algorithm performs a certain number of iterations to find several
  /// cliques and selects the largest one among them. Various limits can be
  /// specified to control the running time and the effectiveness of the
  /// search process.
  ///
  /// \tparam GR The undirected graph type the algorithm runs on.
  ///
  /// \note %GrossoLocatelliPullanMc provides three different node selection
  /// rules, from which the most powerful one is used by default.
  /// For more information, see \ref SelectionRule.
  template <typename GR>
  class GrossoLocatelliPullanMc
  {
  public:

    /// \brief Constants for specifying the node selection rule.
    ///
    /// Enum type containing constants for specifying the node selection rule
    /// for the \ref run() function.
    ///
    /// During the algorithm, nodes are selected for addition to the current
    /// clique according to the applied rule.
    /// In general, the PENALTY_BASED rule turned out to be the most powerful
    /// and the most robust, thus it is the default option.
    /// However, another selection rule can be specified using the \ref run()
    /// function with the proper parameter.
    enum SelectionRule {

      /// A node is selected randomly without any evaluation at each step.
      RANDOM,

      /// A node of maximum degree is selected randomly at each step.
      DEGREE_BASED,

      /// A node of minimum penalty is selected randomly at each step.
      /// The node penalties are updated adaptively after each stage of the
      /// search process.
      PENALTY_BASED
    };

    /// \brief Constants for the causes of search termination.
    ///
    /// Enum type containing constants for the different causes of search
    /// termination. The \ref run() function returns one of these values.
    enum TerminationCause {

      /// The iteration count limit is reached.
      ITERATION_LIMIT,

      /// The step count limit is reached.
      STEP_LIMIT,

      /// The clique size limit is reached.
      SIZE_LIMIT
    };

  private:

    TEMPLATE_GRAPH_TYPEDEFS(GR);

    typedef std::vector<int> IntVector;
    typedef std::vector<char> BoolVector;
    typedef std::vector<BoolVector> BoolMatrix;
    // Note: vector<char> is used instead of vector<bool> for efficiency reasons

    // The underlying graph
    const GR &_graph;
    IntNodeMap _id;

    // Internal matrix representation of the graph
    BoolMatrix _gr;
    int _n;
    
    // Search options
    bool _delta_based_restart;
    int _restart_delta_limit;
 
    // Search limits
    int _iteration_limit;
    int _step_limit;
    int _size_limit;

    // The current clique
    BoolVector _clique;
    int _size;

    // The best clique found so far
    BoolVector _best_clique;
    int _best_size;

    // The "distances" of the nodes from the current clique.
    // _delta[u] is the number of nodes in the clique that are
    // not connected with u.
    IntVector _delta;

    // The current tabu set
    BoolVector _tabu;

    // Random number generator
    Random _rnd;

  private:

    // Implementation of the RANDOM node selection rule.
    class RandomSelectionRule
    {
    private:

      // References to the algorithm instance
      const BoolVector &_clique;
      const IntVector  &_delta;
      const BoolVector &_tabu;
      Random &_rnd;

      // Pivot rule data
      int _n;

    public:

      // Constructor
      RandomSelectionRule(GrossoLocatelliPullanMc &mc) :
        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
        _rnd(mc._rnd), _n(mc._n)
      {}

      // Return a node index for a feasible add move or -1 if no one exists
      int nextFeasibleAddNode() const {
        int start_node = _rnd[_n];
        for (int i = start_node; i != _n; i++) {
          if (_delta[i] == 0 && !_tabu[i]) return i;
        }
        for (int i = 0; i != start_node; i++) {
          if (_delta[i] == 0 && !_tabu[i]) return i;
        }
        return -1;
      }

      // Return a node index for a feasible swap move or -1 if no one exists
      int nextFeasibleSwapNode() const {
        int start_node = _rnd[_n];
        for (int i = start_node; i != _n; i++) {
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i]) return i;
        }
        for (int i = 0; i != start_node; i++) {
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i]) return i;
        }
        return -1;
      }

      // Return a node index for an add move or -1 if no one exists
      int nextAddNode() const {
        int start_node = _rnd[_n];
        for (int i = start_node; i != _n; i++) {
          if (_delta[i] == 0) return i;
        }
        for (int i = 0; i != start_node; i++) {
          if (_delta[i] == 0) return i;
        }
        return -1;
      }

      // Update internal data structures between stages (if necessary)
      void update() {}

    }; //class RandomSelectionRule


    // Implementation of the DEGREE_BASED node selection rule.
    class DegreeBasedSelectionRule
    {
    private:

      // References to the algorithm instance
      const BoolVector &_clique;
      const IntVector  &_delta;
      const BoolVector &_tabu;
      Random &_rnd;

      // Pivot rule data
      int _n;
      IntVector _deg;

    public:

      // Constructor
      DegreeBasedSelectionRule(GrossoLocatelliPullanMc &mc) :
        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
        _rnd(mc._rnd), _n(mc._n), _deg(_n)
      {
        for (int i = 0; i != _n; i++) {
          int d = 0;
          BoolVector &row = mc._gr[i];
          for (int j = 0; j != _n; j++) {
            if (row[j]) d++;
          }
          _deg[i] = d;
        }
      }

      // Return a node index for a feasible add move or -1 if no one exists
      int nextFeasibleAddNode() const {
        int start_node = _rnd[_n];
        int node = -1, max_deg = -1;
        for (int i = start_node; i != _n; i++) {
          if (_delta[i] == 0 && !_tabu[i] && _deg[i] > max_deg) {
            node = i;
            max_deg = _deg[i];
          }
        }
        for (int i = 0; i != start_node; i++) {
          if (_delta[i] == 0 && !_tabu[i] && _deg[i] > max_deg) {
            node = i;
            max_deg = _deg[i];
          }
        }
        return node;
      }

      // Return a node index for a feasible swap move or -1 if no one exists
      int nextFeasibleSwapNode() const {
        int start_node = _rnd[_n];
        int node = -1, max_deg = -1;
        for (int i = start_node; i != _n; i++) {
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
              _deg[i] > max_deg) {
            node = i;
            max_deg = _deg[i];
          }
        }
        for (int i = 0; i != start_node; i++) {
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
              _deg[i] > max_deg) {
            node = i;
            max_deg = _deg[i];
          }
        }
        return node;
      }

      // Return a node index for an add move or -1 if no one exists
      int nextAddNode() const {
        int start_node = _rnd[_n];
        int node = -1, max_deg = -1;
        for (int i = start_node; i != _n; i++) {
          if (_delta[i] == 0 && _deg[i] > max_deg) {
            node = i;
            max_deg = _deg[i];
          }
        }
        for (int i = 0; i != start_node; i++) {
          if (_delta[i] == 0 && _deg[i] > max_deg) {
            node = i;
            max_deg = _deg[i];
          }
        }
        return node;
      }

      // Update internal data structures between stages (if necessary)
      void update() {}

    }; //class DegreeBasedSelectionRule


    // Implementation of the PENALTY_BASED node selection rule.
    class PenaltyBasedSelectionRule
    {
    private:

      // References to the algorithm instance
      const BoolVector &_clique;
      const IntVector  &_delta;
      const BoolVector &_tabu;
      Random &_rnd;

      // Pivot rule data
      int _n;
      IntVector _penalty;

    public:

      // Constructor
      PenaltyBasedSelectionRule(GrossoLocatelliPullanMc &mc) :
        _clique(mc._clique), _delta(mc._delta), _tabu(mc._tabu),
        _rnd(mc._rnd), _n(mc._n), _penalty(_n, 0)
      {}

      // Return a node index for a feasible add move or -1 if no one exists
      int nextFeasibleAddNode() const {
        int start_node = _rnd[_n];
        int node = -1, min_p = std::numeric_limits<int>::max();
        for (int i = start_node; i != _n; i++) {
          if (_delta[i] == 0 && !_tabu[i] && _penalty[i] < min_p) {
            node = i;
            min_p = _penalty[i];
          }
        }
        for (int i = 0; i != start_node; i++) {
          if (_delta[i] == 0 && !_tabu[i] && _penalty[i] < min_p) {
            node = i;
            min_p = _penalty[i];
          }
        }
        return node;
      }

      // Return a node index for a feasible swap move or -1 if no one exists
      int nextFeasibleSwapNode() const {
        int start_node = _rnd[_n];
        int node = -1, min_p = std::numeric_limits<int>::max();
        for (int i = start_node; i != _n; i++) {
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
              _penalty[i] < min_p) {
            node = i;
            min_p = _penalty[i];
          }
        }
        for (int i = 0; i != start_node; i++) {
          if (!_clique[i] && _delta[i] == 1 && !_tabu[i] &&
              _penalty[i] < min_p) {
            node = i;
            min_p = _penalty[i];
          }
        }
        return node;
      }

      // Return a node index for an add move or -1 if no one exists
      int nextAddNode() const {
        int start_node = _rnd[_n];
        int node = -1, min_p = std::numeric_limits<int>::max();
        for (int i = start_node; i != _n; i++) {
          if (_delta[i] == 0 && _penalty[i] < min_p) {
            node = i;
            min_p = _penalty[i];
          }
        }
        for (int i = 0; i != start_node; i++) {
          if (_delta[i] == 0 && _penalty[i] < min_p) {
            node = i;
            min_p = _penalty[i];
          }
        }
        return node;
      }

      // Update internal data structures between stages (if necessary)
      void update() {}

    }; //class PenaltyBasedSelectionRule

  public:

    /// \brief Constructor.
    ///
    /// Constructor.
    /// The global \ref rnd "random number generator instance" is used
    /// during the algorithm.
    ///
    /// \param graph The undirected graph the algorithm runs on.
    GrossoLocatelliPullanMc(const GR& graph) :
      _graph(graph), _id(_graph), _rnd(rnd)
    {
      initOptions();
    }

    /// \brief Constructor with random seed.
    ///
    /// Constructor with random seed.
    ///
    /// \param graph The undirected graph the algorithm runs on.
    /// \param seed Seed value for the internal random number generator
    /// that is used during the algorithm.
    GrossoLocatelliPullanMc(const GR& graph, int seed) :
      _graph(graph), _id(_graph), _rnd(seed)
    {
      initOptions();
    }

    /// \brief Constructor with random number generator.
    ///
    /// Constructor with random number generator.
    ///
    /// \param graph The undirected graph the algorithm runs on.
    /// \param random A random number generator that is used during the
    /// algorithm.
    GrossoLocatelliPullanMc(const GR& graph, const Random& random) :
      _graph(graph), _id(_graph), _rnd(random)
    {
      initOptions();
    }

    /// \name Execution Control
    /// The \ref run() function can be used to execute the algorithm.\n
    /// The functions \ref iterationLimit(int), \ref stepLimit(int), and 
    /// \ref sizeLimit(int) can be used to specify various limits for the
    /// search process.
    
    /// @{
    
    /// \brief Sets the maximum number of iterations.
    ///
    /// This function sets the maximum number of iterations.
    /// Each iteration of the algorithm finds a maximal clique (but not
    /// necessarily the largest one) by performing several search steps
    /// (node selections).
    ///
    /// This limit controls the running time and the success of the
    /// algorithm. For larger values, the algorithm runs slower, but it more
    /// likely finds larger cliques. For smaller values, the algorithm is
    /// faster but probably gives worse results.
    /// 
    /// The default value is \c 1000.
    /// \c -1 means that number of iterations is not limited.
    ///
    /// \warning You should specify a reasonable limit for the number of
    /// iterations and/or the number of search steps.
    ///
    /// \return <tt>(*this)</tt>
    ///
    /// \sa stepLimit(int)
    /// \sa sizeLimit(int)
    GrossoLocatelliPullanMc& iterationLimit(int limit) {
      _iteration_limit = limit;
      return *this;
    }
    
    /// \brief Sets the maximum number of search steps.
    ///
    /// This function sets the maximum number of elementary search steps.
    /// Each iteration of the algorithm finds a maximal clique (but not
    /// necessarily the largest one) by performing several search steps
    /// (node selections).
    ///
    /// This limit controls the running time and the success of the
    /// algorithm. For larger values, the algorithm runs slower, but it more
    /// likely finds larger cliques. For smaller values, the algorithm is
    /// faster but probably gives worse results.
    /// 
    /// The default value is \c -1, which means that number of steps
    /// is not limited explicitly. However, the number of iterations is
    /// limited and each iteration performs a finite number of search steps.
    ///
    /// \warning You should specify a reasonable limit for the number of
    /// iterations and/or the number of search steps.
    ///
    /// \return <tt>(*this)</tt>
    ///
    /// \sa iterationLimit(int)
    /// \sa sizeLimit(int)
    GrossoLocatelliPullanMc& stepLimit(int limit) {
      _step_limit = limit;
      return *this;
    }
    
    /// \brief Sets the desired clique size.
    ///
    /// This function sets the desired clique size that serves as a search
    /// limit. If a clique of this size (or a larger one) is found, then the
    /// algorithm terminates.
    /// 
    /// This function is especially useful if you know an exact upper bound
    /// for the size of the cliques in the graph or if any clique above 
    /// a certain size limit is sufficient for your application.
    /// 
    /// The default value is \c -1, which means that the size limit is set to
    /// the number of nodes in the graph.
    ///
    /// \return <tt>(*this)</tt>
    ///
    /// \sa iterationLimit(int)
    /// \sa stepLimit(int)
    GrossoLocatelliPullanMc& sizeLimit(int limit) {
      _size_limit = limit;
      return *this;
    }
    
    /// \brief The maximum number of iterations.
    ///
    /// This function gives back the maximum number of iterations.
    /// \c -1 means that no limit is specified.
    ///
    /// \sa iterationLimit(int)
    int iterationLimit() const {
      return _iteration_limit;
    }
    
    /// \brief The maximum number of search steps.
    ///
    /// This function gives back the maximum number of search steps.
    /// \c -1 means that no limit is specified.
    ///
    /// \sa stepLimit(int)
    int stepLimit() const {
      return _step_limit;
    }
    
    /// \brief The desired clique size.
    ///
    /// This function gives back the desired clique size that serves as a
    /// search limit. \c -1 means that this limit is set to the number of
    /// nodes in the graph.
    ///
    /// \sa sizeLimit(int)
    int sizeLimit() const {
      return _size_limit;
    }

    /// \brief Runs the algorithm.
    ///
    /// This function runs the algorithm. If one of the specified limits
    /// is reached, the search process terminates.
    ///
    /// \param rule The node selection rule. For more information, see
    /// \ref SelectionRule.
    ///
    /// \return The termination cause of the search. For more information,
    /// see \ref TerminationCause.
    TerminationCause run(SelectionRule rule = PENALTY_BASED)
    {
      init();
      switch (rule) {
        case RANDOM:
          return start<RandomSelectionRule>();
        case DEGREE_BASED:
          return start<DegreeBasedSelectionRule>();
        default:
          return start<PenaltyBasedSelectionRule>();
      }
    }

    /// @}

    /// \name Query Functions
    /// The results of the algorithm can be obtained using these functions.\n
    /// The run() function must be called before using them. 

    /// @{

    /// \brief The size of the found clique
    ///
    /// This function returns the size of the found clique.
    ///
    /// \pre run() must be called before using this function.
    int cliqueSize() const {
      return _best_size;
    }

    /// \brief Gives back the found clique in a \c bool node map
    ///
    /// This function gives back the characteristic vector of the found
    /// clique in the given node map.
    /// It must be a \ref concepts::WriteMap "writable" node map with
    /// \c bool (or convertible) value type.
    ///
    /// \pre run() must be called before using this function.
    template <typename CliqueMap>
    void cliqueMap(CliqueMap &map) const {
      for (NodeIt n(_graph); n != INVALID; ++n) {
        map[n] = static_cast<bool>(_best_clique[_id[n]]);
      }
    }

    /// \brief Iterator to list the nodes of the found clique
    ///
    /// This iterator class lists the nodes of the found clique.
    /// Before using it, you must allocate a GrossoLocatelliPullanMc instance
    /// and call its \ref GrossoLocatelliPullanMc::run() "run()" method.
    ///
    /// The following example prints out the IDs of the nodes in the found
    /// clique.
    /// \code
    ///   GrossoLocatelliPullanMc<Graph> mc(g);
    ///   mc.run();
    ///   for (GrossoLocatelliPullanMc<Graph>::CliqueNodeIt n(mc);
    ///        n != INVALID; ++n)
    ///   {
    ///     std::cout << g.id(n) << std::endl;
    ///   }
    /// \endcode
    class CliqueNodeIt
    {
    private:
      NodeIt _it;
      BoolNodeMap _map;

    public:

      /// Constructor

      /// Constructor.
      /// \param mc The algorithm instance.
      CliqueNodeIt(const GrossoLocatelliPullanMc &mc)
       : _map(mc._graph)
      {
        mc.cliqueMap(_map);
        for (_it = NodeIt(mc._graph); _it != INVALID && !_map[_it]; ++_it) ;
      }

      /// Conversion to \c Node
      operator Node() const { return _it; }

      bool operator==(Invalid) const { return _it == INVALID; }
      bool operator!=(Invalid) const { return _it != INVALID; }

      /// Next node
      CliqueNodeIt &operator++() {
        for (++_it; _it != INVALID && !_map[_it]; ++_it) ;
        return *this;
      }

      /// Postfix incrementation

      /// Postfix incrementation.
      ///
      /// \warning This incrementation returns a \c Node, not a
      /// \c CliqueNodeIt as one may expect.
      typename GR::Node operator++(int) {
        Node n=*this;
        ++(*this);
        return n;
      }

    };

    /// @}

  private:
  
    // Initialize search options and limits
    void initOptions() {
      // Search options
      _delta_based_restart = true;
      _restart_delta_limit = 4;
     
      // Search limits
      _iteration_limit = 1000;
      _step_limit = -1;             // this is disabled by default
      _size_limit = -1;             // this is disabled by default
    }

    // Adds a node to the current clique
    void addCliqueNode(int u) {
      if (_clique[u]) return;
      _clique[u] = true;
      _size++;
      BoolVector &row = _gr[u];
      for (int i = 0; i != _n; i++) {
        if (!row[i]) _delta[i]++;
      }
    }

    // Removes a node from the current clique
    void delCliqueNode(int u) {
      if (!_clique[u]) return;
      _clique[u] = false;
      _size--;
      BoolVector &row = _gr[u];
      for (int i = 0; i != _n; i++) {
        if (!row[i]) _delta[i]--;
      }
    }

    // Initialize data structures
    void init() {
      _n = countNodes(_graph);
      int ui = 0;
      for (NodeIt u(_graph); u != INVALID; ++u) {
        _id[u] = ui++;
      }
      _gr.clear();
      _gr.resize(_n, BoolVector(_n, false));
      ui = 0;
      for (NodeIt u(_graph); u != INVALID; ++u) {
        for (IncEdgeIt e(_graph, u); e != INVALID; ++e) {
          int vi = _id[_graph.runningNode(e)];
          _gr[ui][vi] = true;
          _gr[vi][ui] = true;
        }
        ++ui;
      }

      _clique.clear();
      _clique.resize(_n, false);
      _size = 0;
      _best_clique.clear();
      _best_clique.resize(_n, false);
      _best_size = 0;
      _delta.clear();
      _delta.resize(_n, 0);
      _tabu.clear();
      _tabu.resize(_n, false);
    }

    // Executes the algorithm
    template <typename SelectionRuleImpl>
    TerminationCause start() {
      if (_n == 0) return SIZE_LIMIT;
      if (_n == 1) {
        _best_clique[0] = true;
        _best_size = 1;
        return SIZE_LIMIT;
      }

      // Iterated local search algorithm
      const int max_size = _size_limit >= 0 ? _size_limit : _n;
      const int max_restart = _iteration_limit >= 0 ?
        _iteration_limit : std::numeric_limits<int>::max();
      const int max_select = _step_limit >= 0 ?
        _step_limit : std::numeric_limits<int>::max();

      SelectionRuleImpl sel_method(*this);
      int select = 0, restart = 0;
      IntVector restart_nodes;
      while (select < max_select && restart < max_restart) {

        // Perturbation/restart
        restart++;
        if (_delta_based_restart) {
          restart_nodes.clear();
          for (int i = 0; i != _n; i++) {
            if (_delta[i] >= _restart_delta_limit)
              restart_nodes.push_back(i);
          }
        }
        int rs_node = -1;
        if (restart_nodes.size() > 0) {
          rs_node = restart_nodes[_rnd[restart_nodes.size()]];
        } else {
          rs_node = _rnd[_n];
        }
        BoolVector &row = _gr[rs_node];
        for (int i = 0; i != _n; i++) {
          if (_clique[i] && !row[i]) delCliqueNode(i);
        }
        addCliqueNode(rs_node);

        // Local search
        _tabu.clear();
        _tabu.resize(_n, false);
        bool tabu_empty = true;
        int max_swap = _size;
        while (select < max_select) {
          select++;
          int u;
          if ((u = sel_method.nextFeasibleAddNode()) != -1) {
            // Feasible add move
            addCliqueNode(u);
            if (tabu_empty) max_swap = _size;
          }
          else if ((u = sel_method.nextFeasibleSwapNode()) != -1) {
            // Feasible swap move
            int v = -1;
            BoolVector &row = _gr[u];
            for (int i = 0; i != _n; i++) {
              if (_clique[i] && !row[i]) {
                v = i;
                break;
              }
            }
            addCliqueNode(u);
            delCliqueNode(v);
            _tabu[v] = true;
            tabu_empty = false;
            if (--max_swap <= 0) break;
          }
          else if ((u = sel_method.nextAddNode()) != -1) {
            // Non-feasible add move
            addCliqueNode(u);
          }
          else break;
        }
        if (_size > _best_size) {
          _best_clique = _clique;
          _best_size = _size;
          if (_best_size >= max_size) return SIZE_LIMIT;
        }
        sel_method.update();
      }

      return (restart >= max_restart ? ITERATION_LIMIT : STEP_LIMIT);
    }

  }; //class GrossoLocatelliPullanMc

  ///@}

} //namespace lemon

#endif //LEMON_GROSSO_LOCATELLI_PULLAN_MC_H