1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpnet08.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,241 @@
1.4 +/* glpnet08.c */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Two subroutines sub() and wclique() below are intended to find a
1.10 +* maximum weight clique in a given undirected graph. These subroutines
1.11 +* are slightly modified version of the program WCLIQUE developed by
1.12 +* Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based
1.13 +* on ideas from the article "P. R. J. Ostergard, A new algorithm for
1.14 +* the maximum-weight clique problem, submitted for publication", which
1.15 +* in turn is a generalization of the algorithm for unweighted graphs
1.16 +* presented in "P. R. J. Ostergard, A fast algorithm for the maximum
1.17 +* clique problem, submitted for publication".
1.18 +*
1.19 +* USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE.
1.20 +*
1.21 +* Changes were made by Andrew Makhorin <mao@gnu.org>.
1.22 +*
1.23 +* GLPK is free software: you can redistribute it and/or modify it
1.24 +* under the terms of the GNU General Public License as published by
1.25 +* the Free Software Foundation, either version 3 of the License, or
1.26 +* (at your option) any later version.
1.27 +*
1.28 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.29 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.30 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.31 +* License for more details.
1.32 +*
1.33 +* You should have received a copy of the GNU General Public License
1.34 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.35 +***********************************************************************/
1.36 +
1.37 +#include "glpenv.h"
1.38 +#include "glpnet.h"
1.39 +
1.40 +/***********************************************************************
1.41 +* NAME
1.42 +*
1.43 +* wclique - find maximum weight clique with Ostergard's algorithm
1.44 +*
1.45 +* SYNOPSIS
1.46 +*
1.47 +* int wclique(int n, const int w[], const unsigned char a[],
1.48 +* int ind[]);
1.49 +*
1.50 +* DESCRIPTION
1.51 +*
1.52 +* The routine wclique finds a maximum weight clique in an undirected
1.53 +* graph with Ostergard's algorithm.
1.54 +*
1.55 +* INPUT PARAMETERS
1.56 +*
1.57 +* n is the number of vertices, n > 0.
1.58 +*
1.59 +* w[i], i = 1,...,n, is a weight of vertex i.
1.60 +*
1.61 +* a[*] is the strict (without main diagonal) lower triangle of the
1.62 +* graph adjacency matrix in packed format.
1.63 +*
1.64 +* OUTPUT PARAMETER
1.65 +*
1.66 +* ind[k], k = 1,...,size, is the number of a vertex included in the
1.67 +* clique found, 1 <= ind[k] <= n, where size is the number of vertices
1.68 +* in the clique returned on exit.
1.69 +*
1.70 +* RETURNS
1.71 +*
1.72 +* The routine returns the clique size, i.e. the number of vertices in
1.73 +* the clique. */
1.74 +
1.75 +struct csa
1.76 +{ /* common storage area */
1.77 + int n;
1.78 + /* number of vertices */
1.79 + const int *wt; /* int wt[0:n-1]; */
1.80 + /* weights */
1.81 + const unsigned char *a;
1.82 + /* adjacency matrix (packed lower triangle without main diag.) */
1.83 + int record;
1.84 + /* weight of best clique */
1.85 + int rec_level;
1.86 + /* number of vertices in best clique */
1.87 + int *rec; /* int rec[0:n-1]; */
1.88 + /* best clique so far */
1.89 + int *clique; /* int clique[0:n-1]; */
1.90 + /* table for pruning */
1.91 + int *set; /* int set[0:n-1]; */
1.92 + /* current clique */
1.93 +};
1.94 +
1.95 +#define n (csa->n)
1.96 +#define wt (csa->wt)
1.97 +#define a (csa->a)
1.98 +#define record (csa->record)
1.99 +#define rec_level (csa->rec_level)
1.100 +#define rec (csa->rec)
1.101 +#define clique (csa->clique)
1.102 +#define set (csa->set)
1.103 +
1.104 +#if 0
1.105 +static int is_edge(struct csa *csa, int i, int j)
1.106 +{ /* if there is arc (i,j), the routine returns true; otherwise
1.107 + false; 0 <= i, j < n */
1.108 + int k;
1.109 + xassert(0 <= i && i < n);
1.110 + xassert(0 <= j && j < n);
1.111 + if (i == j) return 0;
1.112 + if (i < j) k = i, i = j, j = k;
1.113 + k = (i * (i - 1)) / 2 + j;
1.114 + return a[k / CHAR_BIT] &
1.115 + (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
1.116 +}
1.117 +#else
1.118 +#define is_edge(csa, i, j) ((i) == (j) ? 0 : \
1.119 + (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))
1.120 +#define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))
1.121 +#define is_edge2(k) (a[(k) / CHAR_BIT] & \
1.122 + (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))
1.123 +#endif
1.124 +
1.125 +static void sub(struct csa *csa, int ct, int table[], int level,
1.126 + int weight, int l_weight)
1.127 +{ int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
1.128 + newtable = xcalloc(n, sizeof(int));
1.129 + if (ct <= 0)
1.130 + { /* 0 or 1 elements left; include these */
1.131 + if (ct == 0)
1.132 + { set[level++] = table[0];
1.133 + weight += l_weight;
1.134 + }
1.135 + if (weight > record)
1.136 + { record = weight;
1.137 + rec_level = level;
1.138 + for (i = 0; i < level; i++) rec[i] = set[i];
1.139 + }
1.140 + goto done;
1.141 + }
1.142 + for (i = ct; i >= 0; i--)
1.143 + { if ((level == 0) && (i < ct)) goto done;
1.144 + k = table[i];
1.145 + if ((level > 0) && (clique[k] <= (record - weight)))
1.146 + goto done; /* prune */
1.147 + set[level] = k;
1.148 + curr_weight = weight + wt[k];
1.149 + l_weight -= wt[k];
1.150 + if (l_weight <= (record - curr_weight))
1.151 + goto done; /* prune */
1.152 + p1 = newtable;
1.153 + p2 = table;
1.154 + left_weight = 0;
1.155 + while (p2 < table + i)
1.156 + { j = *p2++;
1.157 + if (is_edge(csa, j, k))
1.158 + { *p1++ = j;
1.159 + left_weight += wt[j];
1.160 + }
1.161 + }
1.162 + if (left_weight <= (record - curr_weight)) continue;
1.163 + sub(csa, p1 - newtable - 1, newtable, level + 1, curr_weight,
1.164 + left_weight);
1.165 + }
1.166 +done: xfree(newtable);
1.167 + return;
1.168 +}
1.169 +
1.170 +int wclique(int _n, const int w[], const unsigned char _a[], int ind[])
1.171 +{ struct csa _csa, *csa = &_csa;
1.172 + int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
1.173 + glp_long timer;
1.174 + n = _n;
1.175 + xassert(n > 0);
1.176 + wt = &w[1];
1.177 + a = _a;
1.178 + record = 0;
1.179 + rec_level = 0;
1.180 + rec = &ind[1];
1.181 + clique = xcalloc(n, sizeof(int));
1.182 + set = xcalloc(n, sizeof(int));
1.183 + used = xcalloc(n, sizeof(int));
1.184 + nwt = xcalloc(n, sizeof(int));
1.185 + pos = xcalloc(n, sizeof(int));
1.186 + /* start timer */
1.187 + timer = xtime();
1.188 + /* order vertices */
1.189 + for (i = 0; i < n; i++)
1.190 + { nwt[i] = 0;
1.191 + for (j = 0; j < n; j++)
1.192 + if (is_edge(csa, i, j)) nwt[i] += wt[j];
1.193 + }
1.194 + for (i = 0; i < n; i++)
1.195 + used[i] = 0;
1.196 + for (i = n-1; i >= 0; i--)
1.197 + { max_wt = -1;
1.198 + max_nwt = -1;
1.199 + for (j = 0; j < n; j++)
1.200 + { if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
1.201 + && nwt[j] > max_nwt)))
1.202 + { max_wt = wt[j];
1.203 + max_nwt = nwt[j];
1.204 + p = j;
1.205 + }
1.206 + }
1.207 + pos[i] = p;
1.208 + used[p] = 1;
1.209 + for (j = 0; j < n; j++)
1.210 + if ((!used[j]) && (j != p) && (is_edge(csa, p, j)))
1.211 + nwt[j] -= wt[p];
1.212 + }
1.213 + /* main routine */
1.214 + wth = 0;
1.215 + for (i = 0; i < n; i++)
1.216 + { wth += wt[pos[i]];
1.217 + sub(csa, i, pos, 0, 0, wth);
1.218 + clique[pos[i]] = record;
1.219 + if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
1.220 + { /* print current record and reset timer */
1.221 + xprintf("level = %d (%d); best = %d\n", i+1, n, record);
1.222 + timer = xtime();
1.223 + }
1.224 + }
1.225 + xfree(clique);
1.226 + xfree(set);
1.227 + xfree(used);
1.228 + xfree(nwt);
1.229 + xfree(pos);
1.230 + /* return the solution found */
1.231 + for (i = 1; i <= rec_level; i++) ind[i]++;
1.232 + return rec_level;
1.233 +}
1.234 +
1.235 +#undef n
1.236 +#undef wt
1.237 +#undef a
1.238 +#undef record
1.239 +#undef rec_level
1.240 +#undef rec
1.241 +#undef clique
1.242 +#undef set
1.243 +
1.244 +/* eof */