1 | /* glpnet08.c */ |
---|
2 | |
---|
3 | /*********************************************************************** |
---|
4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
5 | * |
---|
6 | * Two subroutines sub() and wclique() below are intended to find a |
---|
7 | * maximum weight clique in a given undirected graph. These subroutines |
---|
8 | * are slightly modified version of the program WCLIQUE developed by |
---|
9 | * Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based |
---|
10 | * on ideas from the article "P. R. J. Ostergard, A new algorithm for |
---|
11 | * the maximum-weight clique problem, submitted for publication", which |
---|
12 | * in turn is a generalization of the algorithm for unweighted graphs |
---|
13 | * presented in "P. R. J. Ostergard, A fast algorithm for the maximum |
---|
14 | * clique problem, submitted for publication". |
---|
15 | * |
---|
16 | * USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE. |
---|
17 | * |
---|
18 | * Changes were made by Andrew Makhorin <mao@gnu.org>. |
---|
19 | * |
---|
20 | * GLPK is free software: you can redistribute it and/or modify it |
---|
21 | * under the terms of the GNU General Public License as published by |
---|
22 | * the Free Software Foundation, either version 3 of the License, or |
---|
23 | * (at your option) any later version. |
---|
24 | * |
---|
25 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
26 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
27 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
28 | * License for more details. |
---|
29 | * |
---|
30 | * You should have received a copy of the GNU General Public License |
---|
31 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
32 | ***********************************************************************/ |
---|
33 | |
---|
34 | #include "glpenv.h" |
---|
35 | #include "glpnet.h" |
---|
36 | |
---|
37 | /*********************************************************************** |
---|
38 | * NAME |
---|
39 | * |
---|
40 | * wclique - find maximum weight clique with Ostergard's algorithm |
---|
41 | * |
---|
42 | * SYNOPSIS |
---|
43 | * |
---|
44 | * int wclique(int n, const int w[], const unsigned char a[], |
---|
45 | * int ind[]); |
---|
46 | * |
---|
47 | * DESCRIPTION |
---|
48 | * |
---|
49 | * The routine wclique finds a maximum weight clique in an undirected |
---|
50 | * graph with Ostergard's algorithm. |
---|
51 | * |
---|
52 | * INPUT PARAMETERS |
---|
53 | * |
---|
54 | * n is the number of vertices, n > 0. |
---|
55 | * |
---|
56 | * w[i], i = 1,...,n, is a weight of vertex i. |
---|
57 | * |
---|
58 | * a[*] is the strict (without main diagonal) lower triangle of the |
---|
59 | * graph adjacency matrix in packed format. |
---|
60 | * |
---|
61 | * OUTPUT PARAMETER |
---|
62 | * |
---|
63 | * ind[k], k = 1,...,size, is the number of a vertex included in the |
---|
64 | * clique found, 1 <= ind[k] <= n, where size is the number of vertices |
---|
65 | * in the clique returned on exit. |
---|
66 | * |
---|
67 | * RETURNS |
---|
68 | * |
---|
69 | * The routine returns the clique size, i.e. the number of vertices in |
---|
70 | * the clique. */ |
---|
71 | |
---|
72 | struct csa |
---|
73 | { /* common storage area */ |
---|
74 | int n; |
---|
75 | /* number of vertices */ |
---|
76 | const int *wt; /* int wt[0:n-1]; */ |
---|
77 | /* weights */ |
---|
78 | const unsigned char *a; |
---|
79 | /* adjacency matrix (packed lower triangle without main diag.) */ |
---|
80 | int record; |
---|
81 | /* weight of best clique */ |
---|
82 | int rec_level; |
---|
83 | /* number of vertices in best clique */ |
---|
84 | int *rec; /* int rec[0:n-1]; */ |
---|
85 | /* best clique so far */ |
---|
86 | int *clique; /* int clique[0:n-1]; */ |
---|
87 | /* table for pruning */ |
---|
88 | int *set; /* int set[0:n-1]; */ |
---|
89 | /* current clique */ |
---|
90 | }; |
---|
91 | |
---|
92 | #define n (csa->n) |
---|
93 | #define wt (csa->wt) |
---|
94 | #define a (csa->a) |
---|
95 | #define record (csa->record) |
---|
96 | #define rec_level (csa->rec_level) |
---|
97 | #define rec (csa->rec) |
---|
98 | #define clique (csa->clique) |
---|
99 | #define set (csa->set) |
---|
100 | |
---|
101 | #if 0 |
---|
102 | static int is_edge(struct csa *csa, int i, int j) |
---|
103 | { /* if there is arc (i,j), the routine returns true; otherwise |
---|
104 | false; 0 <= i, j < n */ |
---|
105 | int k; |
---|
106 | xassert(0 <= i && i < n); |
---|
107 | xassert(0 <= j && j < n); |
---|
108 | if (i == j) return 0; |
---|
109 | if (i < j) k = i, i = j, j = k; |
---|
110 | k = (i * (i - 1)) / 2 + j; |
---|
111 | return a[k / CHAR_BIT] & |
---|
112 | (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT)); |
---|
113 | } |
---|
114 | #else |
---|
115 | #define is_edge(csa, i, j) ((i) == (j) ? 0 : \ |
---|
116 | (i) > (j) ? is_edge1(i, j) : is_edge1(j, i)) |
---|
117 | #define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j)) |
---|
118 | #define is_edge2(k) (a[(k) / CHAR_BIT] & \ |
---|
119 | (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT))) |
---|
120 | #endif |
---|
121 | |
---|
122 | static void sub(struct csa *csa, int ct, int table[], int level, |
---|
123 | int weight, int l_weight) |
---|
124 | { int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable; |
---|
125 | newtable = xcalloc(n, sizeof(int)); |
---|
126 | if (ct <= 0) |
---|
127 | { /* 0 or 1 elements left; include these */ |
---|
128 | if (ct == 0) |
---|
129 | { set[level++] = table[0]; |
---|
130 | weight += l_weight; |
---|
131 | } |
---|
132 | if (weight > record) |
---|
133 | { record = weight; |
---|
134 | rec_level = level; |
---|
135 | for (i = 0; i < level; i++) rec[i] = set[i]; |
---|
136 | } |
---|
137 | goto done; |
---|
138 | } |
---|
139 | for (i = ct; i >= 0; i--) |
---|
140 | { if ((level == 0) && (i < ct)) goto done; |
---|
141 | k = table[i]; |
---|
142 | if ((level > 0) && (clique[k] <= (record - weight))) |
---|
143 | goto done; /* prune */ |
---|
144 | set[level] = k; |
---|
145 | curr_weight = weight + wt[k]; |
---|
146 | l_weight -= wt[k]; |
---|
147 | if (l_weight <= (record - curr_weight)) |
---|
148 | goto done; /* prune */ |
---|
149 | p1 = newtable; |
---|
150 | p2 = table; |
---|
151 | left_weight = 0; |
---|
152 | while (p2 < table + i) |
---|
153 | { j = *p2++; |
---|
154 | if (is_edge(csa, j, k)) |
---|
155 | { *p1++ = j; |
---|
156 | left_weight += wt[j]; |
---|
157 | } |
---|
158 | } |
---|
159 | if (left_weight <= (record - curr_weight)) continue; |
---|
160 | sub(csa, p1 - newtable - 1, newtable, level + 1, curr_weight, |
---|
161 | left_weight); |
---|
162 | } |
---|
163 | done: xfree(newtable); |
---|
164 | return; |
---|
165 | } |
---|
166 | |
---|
167 | int wclique(int _n, const int w[], const unsigned char _a[], int ind[]) |
---|
168 | { struct csa _csa, *csa = &_csa; |
---|
169 | int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos; |
---|
170 | glp_long timer; |
---|
171 | n = _n; |
---|
172 | xassert(n > 0); |
---|
173 | wt = &w[1]; |
---|
174 | a = _a; |
---|
175 | record = 0; |
---|
176 | rec_level = 0; |
---|
177 | rec = &ind[1]; |
---|
178 | clique = xcalloc(n, sizeof(int)); |
---|
179 | set = xcalloc(n, sizeof(int)); |
---|
180 | used = xcalloc(n, sizeof(int)); |
---|
181 | nwt = xcalloc(n, sizeof(int)); |
---|
182 | pos = xcalloc(n, sizeof(int)); |
---|
183 | /* start timer */ |
---|
184 | timer = xtime(); |
---|
185 | /* order vertices */ |
---|
186 | for (i = 0; i < n; i++) |
---|
187 | { nwt[i] = 0; |
---|
188 | for (j = 0; j < n; j++) |
---|
189 | if (is_edge(csa, i, j)) nwt[i] += wt[j]; |
---|
190 | } |
---|
191 | for (i = 0; i < n; i++) |
---|
192 | used[i] = 0; |
---|
193 | for (i = n-1; i >= 0; i--) |
---|
194 | { max_wt = -1; |
---|
195 | max_nwt = -1; |
---|
196 | for (j = 0; j < n; j++) |
---|
197 | { if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt |
---|
198 | && nwt[j] > max_nwt))) |
---|
199 | { max_wt = wt[j]; |
---|
200 | max_nwt = nwt[j]; |
---|
201 | p = j; |
---|
202 | } |
---|
203 | } |
---|
204 | pos[i] = p; |
---|
205 | used[p] = 1; |
---|
206 | for (j = 0; j < n; j++) |
---|
207 | if ((!used[j]) && (j != p) && (is_edge(csa, p, j))) |
---|
208 | nwt[j] -= wt[p]; |
---|
209 | } |
---|
210 | /* main routine */ |
---|
211 | wth = 0; |
---|
212 | for (i = 0; i < n; i++) |
---|
213 | { wth += wt[pos[i]]; |
---|
214 | sub(csa, i, pos, 0, 0, wth); |
---|
215 | clique[pos[i]] = record; |
---|
216 | if (xdifftime(xtime(), timer) >= 5.0 - 0.001) |
---|
217 | { /* print current record and reset timer */ |
---|
218 | xprintf("level = %d (%d); best = %d\n", i+1, n, record); |
---|
219 | timer = xtime(); |
---|
220 | } |
---|
221 | } |
---|
222 | xfree(clique); |
---|
223 | xfree(set); |
---|
224 | xfree(used); |
---|
225 | xfree(nwt); |
---|
226 | xfree(pos); |
---|
227 | /* return the solution found */ |
---|
228 | for (i = 1; i <= rec_level; i++) ind[i]++; |
---|
229 | return rec_level; |
---|
230 | } |
---|
231 | |
---|
232 | #undef n |
---|
233 | #undef wt |
---|
234 | #undef a |
---|
235 | #undef record |
---|
236 | #undef rec_level |
---|
237 | #undef rec |
---|
238 | #undef clique |
---|
239 | #undef set |
---|
240 | |
---|
241 | /* eof */ |
---|