alpar@9: /* glpnet08.c */
alpar@9: 
alpar@9: /***********************************************************************
alpar@9: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@9: *
alpar@9: *  Two subroutines sub() and wclique() below are intended to find a
alpar@9: *  maximum weight clique in a given undirected graph. These subroutines
alpar@9: *  are slightly modified version of the program WCLIQUE developed by
alpar@9: *  Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based
alpar@9: *  on ideas from the article "P. R. J. Ostergard, A new algorithm for
alpar@9: *  the maximum-weight clique problem, submitted for publication", which
alpar@9: *  in turn is a generalization of the algorithm for unweighted graphs
alpar@9: *  presented in "P. R. J. Ostergard, A fast algorithm for the maximum
alpar@9: *  clique problem, submitted for publication".
alpar@9: *
alpar@9: *  USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE.
alpar@9: *
alpar@9: *  Changes were made by Andrew Makhorin <mao@gnu.org>.
alpar@9: *
alpar@9: *  GLPK is free software: you can redistribute it and/or modify it
alpar@9: *  under the terms of the GNU General Public License as published by
alpar@9: *  the Free Software Foundation, either version 3 of the License, or
alpar@9: *  (at your option) any later version.
alpar@9: *
alpar@9: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9: *  License for more details.
alpar@9: *
alpar@9: *  You should have received a copy of the GNU General Public License
alpar@9: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9: ***********************************************************************/
alpar@9: 
alpar@9: #include "glpenv.h"
alpar@9: #include "glpnet.h"
alpar@9: 
alpar@9: /***********************************************************************
alpar@9: *  NAME
alpar@9: *
alpar@9: *  wclique - find maximum weight clique with Ostergard's algorithm
alpar@9: *
alpar@9: *  SYNOPSIS
alpar@9: *
alpar@9: *  int wclique(int n, const int w[], const unsigned char a[],
alpar@9: *     int ind[]);
alpar@9: *
alpar@9: *  DESCRIPTION
alpar@9: *
alpar@9: *  The routine wclique finds a maximum weight clique in an undirected
alpar@9: *  graph with Ostergard's algorithm.
alpar@9: *
alpar@9: *  INPUT PARAMETERS
alpar@9: *
alpar@9: *  n is the number of vertices, n > 0.
alpar@9: *
alpar@9: *  w[i], i = 1,...,n, is a weight of vertex i.
alpar@9: *
alpar@9: *  a[*] is the strict (without main diagonal) lower triangle of the
alpar@9: *  graph adjacency matrix in packed format.
alpar@9: *
alpar@9: *  OUTPUT PARAMETER
alpar@9: *
alpar@9: *  ind[k], k = 1,...,size, is the number of a vertex included in the
alpar@9: *  clique found, 1 <= ind[k] <= n, where size is the number of vertices
alpar@9: *  in the clique returned on exit.
alpar@9: *
alpar@9: *  RETURNS
alpar@9: *
alpar@9: *  The routine returns the clique size, i.e. the number of vertices in
alpar@9: *  the clique. */
alpar@9: 
alpar@9: struct csa
alpar@9: {     /* common storage area */
alpar@9:       int n;
alpar@9:       /* number of vertices */
alpar@9:       const int *wt; /* int wt[0:n-1]; */
alpar@9:       /* weights */
alpar@9:       const unsigned char *a;
alpar@9:       /* adjacency matrix (packed lower triangle without main diag.) */
alpar@9:       int record;
alpar@9:       /* weight of best clique */
alpar@9:       int rec_level;
alpar@9:       /* number of vertices in best clique */
alpar@9:       int *rec; /* int rec[0:n-1]; */
alpar@9:       /* best clique so far */
alpar@9:       int *clique; /* int clique[0:n-1]; */
alpar@9:       /* table for pruning */
alpar@9:       int *set; /* int set[0:n-1]; */
alpar@9:       /* current clique */
alpar@9: };
alpar@9: 
alpar@9: #define n         (csa->n)
alpar@9: #define wt        (csa->wt)
alpar@9: #define a         (csa->a)
alpar@9: #define record    (csa->record)
alpar@9: #define rec_level (csa->rec_level)
alpar@9: #define rec       (csa->rec)
alpar@9: #define clique    (csa->clique)
alpar@9: #define set       (csa->set)
alpar@9: 
alpar@9: #if 0
alpar@9: static int is_edge(struct csa *csa, int i, int j)
alpar@9: {     /* if there is arc (i,j), the routine returns true; otherwise
alpar@9:          false; 0 <= i, j < n */
alpar@9:       int k;
alpar@9:       xassert(0 <= i && i < n);
alpar@9:       xassert(0 <= j && j < n);
alpar@9:       if (i == j) return 0;
alpar@9:       if (i < j) k = i, i = j, j = k;
alpar@9:       k = (i * (i - 1)) / 2 + j;
alpar@9:       return a[k / CHAR_BIT] &
alpar@9:          (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
alpar@9: }
alpar@9: #else
alpar@9: #define is_edge(csa, i, j) ((i) == (j) ? 0 : \
alpar@9:       (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))
alpar@9: #define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))
alpar@9: #define is_edge2(k) (a[(k) / CHAR_BIT] & \
alpar@9:       (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))
alpar@9: #endif
alpar@9: 
alpar@9: static void sub(struct csa *csa, int ct, int table[], int level,
alpar@9:       int weight, int l_weight)
alpar@9: {     int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
alpar@9:       newtable = xcalloc(n, sizeof(int));
alpar@9:       if (ct <= 0)
alpar@9:       {  /* 0 or 1 elements left; include these */
alpar@9:          if (ct == 0)
alpar@9:          {  set[level++] = table[0];
alpar@9:             weight += l_weight;
alpar@9:          }
alpar@9:          if (weight > record)
alpar@9:          {  record = weight;
alpar@9:             rec_level = level;
alpar@9:             for (i = 0; i < level; i++) rec[i] = set[i];
alpar@9:          }
alpar@9:          goto done;
alpar@9:       }
alpar@9:       for (i = ct; i >= 0; i--)
alpar@9:       {  if ((level == 0) && (i < ct)) goto done;
alpar@9:          k = table[i];
alpar@9:          if ((level > 0) && (clique[k] <= (record - weight)))
alpar@9:             goto done; /* prune */
alpar@9:          set[level] = k;
alpar@9:          curr_weight = weight + wt[k];
alpar@9:          l_weight -= wt[k];
alpar@9:          if (l_weight <= (record - curr_weight))
alpar@9:             goto done; /* prune */
alpar@9:          p1 = newtable;
alpar@9:          p2 = table;
alpar@9:          left_weight = 0;
alpar@9:          while (p2 < table + i)
alpar@9:          {  j = *p2++;
alpar@9:             if (is_edge(csa, j, k))
alpar@9:             {  *p1++ = j;
alpar@9:                left_weight += wt[j];
alpar@9:             }
alpar@9:          }
alpar@9:          if (left_weight <= (record - curr_weight)) continue;
alpar@9:          sub(csa, p1 - newtable - 1, newtable, level + 1, curr_weight,
alpar@9:             left_weight);
alpar@9:       }
alpar@9: done: xfree(newtable);
alpar@9:       return;
alpar@9: }
alpar@9: 
alpar@9: int wclique(int _n, const int w[], const unsigned char _a[], int ind[])
alpar@9: {     struct csa _csa, *csa = &_csa;
alpar@9:       int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
alpar@9:       glp_long timer;
alpar@9:       n = _n;
alpar@9:       xassert(n > 0);
alpar@9:       wt = &w[1];
alpar@9:       a = _a;
alpar@9:       record = 0;
alpar@9:       rec_level = 0;
alpar@9:       rec = &ind[1];
alpar@9:       clique = xcalloc(n, sizeof(int));
alpar@9:       set = xcalloc(n, sizeof(int));
alpar@9:       used = xcalloc(n, sizeof(int));
alpar@9:       nwt = xcalloc(n, sizeof(int));
alpar@9:       pos = xcalloc(n, sizeof(int));
alpar@9:       /* start timer */
alpar@9:       timer = xtime();
alpar@9:       /* order vertices */
alpar@9:       for (i = 0; i < n; i++)
alpar@9:       {  nwt[i] = 0;
alpar@9:          for (j = 0; j < n; j++)
alpar@9:             if (is_edge(csa, i, j)) nwt[i] += wt[j];
alpar@9:       }
alpar@9:       for (i = 0; i < n; i++)
alpar@9:          used[i] = 0;
alpar@9:       for (i = n-1; i >= 0; i--)
alpar@9:       {  max_wt = -1;
alpar@9:          max_nwt = -1;
alpar@9:          for (j = 0; j < n; j++)
alpar@9:          {  if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
alpar@9:                && nwt[j] > max_nwt)))
alpar@9:             {  max_wt = wt[j];
alpar@9:                max_nwt = nwt[j];
alpar@9:                p = j;
alpar@9:             }
alpar@9:          }
alpar@9:          pos[i] = p;
alpar@9:          used[p] = 1;
alpar@9:          for (j = 0; j < n; j++)
alpar@9:             if ((!used[j]) && (j != p) && (is_edge(csa, p, j)))
alpar@9:                nwt[j] -= wt[p];
alpar@9:       }
alpar@9:       /* main routine */
alpar@9:       wth = 0;
alpar@9:       for (i = 0; i < n; i++)
alpar@9:       {  wth += wt[pos[i]];
alpar@9:          sub(csa, i, pos, 0, 0, wth);
alpar@9:          clique[pos[i]] = record;
alpar@9:          if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
alpar@9:          {  /* print current record and reset timer */
alpar@9:             xprintf("level = %d (%d); best = %d\n", i+1, n, record);
alpar@9:             timer = xtime();
alpar@9:          }
alpar@9:       }
alpar@9:       xfree(clique);
alpar@9:       xfree(set);
alpar@9:       xfree(used);
alpar@9:       xfree(nwt);
alpar@9:       xfree(pos);
alpar@9:       /* return the solution found */
alpar@9:       for (i = 1; i <= rec_level; i++) ind[i]++;
alpar@9:       return rec_level;
alpar@9: }
alpar@9: 
alpar@9: #undef n
alpar@9: #undef wt
alpar@9: #undef a
alpar@9: #undef record
alpar@9: #undef rec_level
alpar@9: #undef rec
alpar@9: #undef clique
alpar@9: #undef set
alpar@9: 
alpar@9: /* eof */