COIN-OR::LEMON - Graph Library

source: lemon/lemon/grosso_locatelli_pullan_mc.h @ 999:c279b19abc62

Last change on this file since 999:c279b19abc62 was 999:c279b19abc62, checked in by Peter Kovacs <kpeter@…>, 10 years ago

Add a heuristic algorithm for the max clique problem (#380)

File size: 19.3 KB
Line 
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
33namespace 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
Note: See TracBrowser for help on using the repository browser.