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