/* glpnet08.c */

/***********************************************************************

*  This code is part of GLPK (GNU Linear Programming Kit).

*

*  Two subroutines sub() and wclique() below are intended to find a

*  maximum weight clique in a given undirected graph. These subroutines

*  are slightly modified version of the program WCLIQUE developed by

*  Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based

*  on ideas from the article "P. R. J. Ostergard, A new algorithm for

*  the maximum-weight clique problem, submitted for publication", which

*  in turn is a generalization of the algorithm for unweighted graphs

*  presented in "P. R. J. Ostergard, A fast algorithm for the maximum

*  clique problem, submitted for publication".

*

*  USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE.

*

*  Changes were made by Andrew Makhorin <mao@gnu.org>.

*

*  GLPK is free software: you can redistribute it and/or modify it

*  under the terms of the GNU General Public License as published by

*  the Free Software Foundation, either version 3 of the License, or

*  (at your option) any later version.

*

*  GLPK is distributed in the hope that it will be useful, but WITHOUT

*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

*  License for more details.

*

*  You should have received a copy of the GNU General Public License

*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.

***********************************************************************/

#include "glpenv.h"

#include "glpnet.h"

/***********************************************************************

*  NAME

*

*  wclique - find maximum weight clique with Ostergard's algorithm

*

*  SYNOPSIS

*

*  int wclique(int n, const int w[], const unsigned char a[],

*     int ind[]);

*

*  DESCRIPTION

*

*  The routine wclique finds a maximum weight clique in an undirected

*  graph with Ostergard's algorithm.

*

*  INPUT PARAMETERS

*

*  n is the number of vertices, n > 0.

*

*  w[i], i = 1,...,n, is a weight of vertex i.

*

*  a[*] is the strict (without main diagonal) lower triangle of the

*  graph adjacency matrix in packed format.

*

*  OUTPUT PARAMETER

*

*  ind[k], k = 1,...,size, is the number of a vertex included in the

*  clique found, 1 <= ind[k] <= n, where size is the number of vertices

*  in the clique returned on exit.

*

*  RETURNS

*

*  The routine returns the clique size, i.e. the number of vertices in

*  the clique. */

struct csa

{     /* common storage area */

      int n;

      /* number of vertices */

      const int *wt; /* int wt[0:n-1]; */

      /* weights */

      const unsigned char *a;

      /* adjacency matrix (packed lower triangle without main diag.) */

      int record;

      /* weight of best clique */

      int rec_level;

      /* number of vertices in best clique */

      int *rec; /* int rec[0:n-1]; */

      /* best clique so far */

      int *clique; /* int clique[0:n-1]; */

      /* table for pruning */

      int *set; /* int set[0:n-1]; */

      /* current clique */

};

#define n         (csa->n)

#define wt        (csa->wt)

#define a         (csa->a)

#define record    (csa->record)

#define rec_level (csa->rec_level)

#define rec       (csa->rec)

#define clique    (csa->clique)

#define set       (csa->set)

#if 0

static int is_edge(struct csa *csa, int i, int j)

{     /* if there is arc (i,j), the routine returns true; otherwise

         false; 0 <= i, j < n */

      int k;

      xassert(0 <= i && i < n);

      xassert(0 <= j && j < n);

      if (i == j) return 0;

      if (i < j) k = i, i = j, j = k;

      k = (i * (i - 1)) / 2 + j;

      return a[k / CHAR_BIT] &

         (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));

}

#else

#define is_edge(csa, i, j) ((i) == (j) ? 0 : \

      (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))

#define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))

#define is_edge2(k) (a[(k) / CHAR_BIT] & \

      (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))

#endif

static void sub(struct csa *csa, int ct, int table[], int level,

      int weight, int l_weight)

{     int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;

      newtable = xcalloc(n, sizeof(int));

      if (ct <= 0)

      {  /* 0 or 1 elements left; include these */

         if (ct == 0)

         {  set[level++] = table[0];

            weight += l_weight;

         }

         if (weight > record)

         {  record = weight;

            rec_level = level;

            for (i = 0; i < level; i++) rec[i] = set[i];

         }

         goto done;

      }

      for (i = ct; i >= 0; i--)

      {  if ((level == 0) && (i < ct)) goto done;

         k = table[i];

         if ((level > 0) && (clique[k] <= (record - weight)))

            goto done; /* prune */

         set[level] = k;

         curr_weight = weight + wt[k];

         l_weight -= wt[k];

         if (l_weight <= (record - curr_weight))

            goto done; /* prune */

         p1 = newtable;

         p2 = table;

         left_weight = 0;

         while (p2 < table + i)

         {  j = *p2++;

            if (is_edge(csa, j, k))

            {  *p1++ = j;

               left_weight += wt[j];

            }

         }

         if (left_weight <= (record - curr_weight)) continue;

         sub(csa, p1 - newtable - 1, newtable, level + 1, curr_weight,

            left_weight);

      }

done: xfree(newtable);

      return;

}

int wclique(int _n, const int w[], const unsigned char _a[], int ind[])

{     struct csa _csa, *csa = &_csa;

      int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;

      glp_long timer;

      n = _n;

      xassert(n > 0);

      wt = &w[1];

      a = _a;

      record = 0;

      rec_level = 0;

      rec = &ind[1];

      clique = xcalloc(n, sizeof(int));

      set = xcalloc(n, sizeof(int));

      used = xcalloc(n, sizeof(int));

      nwt = xcalloc(n, sizeof(int));

      pos = xcalloc(n, sizeof(int));

      /* start timer */

      timer = xtime();

      /* order vertices */

      for (i = 0; i < n; i++)

      {  nwt[i] = 0;

         for (j = 0; j < n; j++)

            if (is_edge(csa, i, j)) nwt[i] += wt[j];

      }

      for (i = 0; i < n; i++)

         used[i] = 0;

      for (i = n-1; i >= 0; i--)

      {  max_wt = -1;

         max_nwt = -1;

         for (j = 0; j < n; j++)

         {  if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt

               && nwt[j] > max_nwt)))

            {  max_wt = wt[j];

               max_nwt = nwt[j];

               p = j;

            }

         }

         pos[i] = p;

         used[p] = 1;

         for (j = 0; j < n; j++)

            if ((!used[j]) && (j != p) && (is_edge(csa, p, j)))

               nwt[j] -= wt[p];

      }

      /* main routine */

      wth = 0;

      for (i = 0; i < n; i++)

      {  wth += wt[pos[i]];

         sub(csa, i, pos, 0, 0, wth);

         clique[pos[i]] = record;

         if (xdifftime(xtime(), timer) >= 5.0 - 0.001)

         {  /* print current record and reset timer */

            xprintf("level = %d (%d); best = %d\n", i+1, n, record);

            timer = xtime();

         }

      }

      xfree(clique);

      xfree(set);

      xfree(used);

      xfree(nwt);

      xfree(pos);

      /* return the solution found */

      for (i = 1; i <= rec_level; i++) ind[i]++;

      return rec_level;

}

#undef n

#undef wt

#undef a

#undef record

#undef rec_level

#undef rec

#undef clique

#undef set

/* eof */