|
1 /* glpios08.c (clique cut generator) */ |
|
2 |
|
3 /*********************************************************************** |
|
4 * This code is part of GLPK (GNU Linear Programming Kit). |
|
5 * |
|
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
|
7 * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
|
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
|
9 * E-mail: <mao@gnu.org>. |
|
10 * |
|
11 * GLPK is free software: you can redistribute it and/or modify it |
|
12 * under the terms of the GNU General Public License as published by |
|
13 * the Free Software Foundation, either version 3 of the License, or |
|
14 * (at your option) any later version. |
|
15 * |
|
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT |
|
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
|
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
|
19 * License for more details. |
|
20 * |
|
21 * You should have received a copy of the GNU General Public License |
|
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
|
23 ***********************************************************************/ |
|
24 |
|
25 #include "glpios.h" |
|
26 |
|
27 static double get_row_lb(LPX *lp, int i) |
|
28 { /* this routine returns lower bound of row i or -DBL_MAX if the |
|
29 row has no lower bound */ |
|
30 double lb; |
|
31 switch (lpx_get_row_type(lp, i)) |
|
32 { case LPX_FR: |
|
33 case LPX_UP: |
|
34 lb = -DBL_MAX; |
|
35 break; |
|
36 case LPX_LO: |
|
37 case LPX_DB: |
|
38 case LPX_FX: |
|
39 lb = lpx_get_row_lb(lp, i); |
|
40 break; |
|
41 default: |
|
42 xassert(lp != lp); |
|
43 } |
|
44 return lb; |
|
45 } |
|
46 |
|
47 static double get_row_ub(LPX *lp, int i) |
|
48 { /* this routine returns upper bound of row i or +DBL_MAX if the |
|
49 row has no upper bound */ |
|
50 double ub; |
|
51 switch (lpx_get_row_type(lp, i)) |
|
52 { case LPX_FR: |
|
53 case LPX_LO: |
|
54 ub = +DBL_MAX; |
|
55 break; |
|
56 case LPX_UP: |
|
57 case LPX_DB: |
|
58 case LPX_FX: |
|
59 ub = lpx_get_row_ub(lp, i); |
|
60 break; |
|
61 default: |
|
62 xassert(lp != lp); |
|
63 } |
|
64 return ub; |
|
65 } |
|
66 |
|
67 static double get_col_lb(LPX *lp, int j) |
|
68 { /* this routine returns lower bound of column j or -DBL_MAX if |
|
69 the column has no lower bound */ |
|
70 double lb; |
|
71 switch (lpx_get_col_type(lp, j)) |
|
72 { case LPX_FR: |
|
73 case LPX_UP: |
|
74 lb = -DBL_MAX; |
|
75 break; |
|
76 case LPX_LO: |
|
77 case LPX_DB: |
|
78 case LPX_FX: |
|
79 lb = lpx_get_col_lb(lp, j); |
|
80 break; |
|
81 default: |
|
82 xassert(lp != lp); |
|
83 } |
|
84 return lb; |
|
85 } |
|
86 |
|
87 static double get_col_ub(LPX *lp, int j) |
|
88 { /* this routine returns upper bound of column j or +DBL_MAX if |
|
89 the column has no upper bound */ |
|
90 double ub; |
|
91 switch (lpx_get_col_type(lp, j)) |
|
92 { case LPX_FR: |
|
93 case LPX_LO: |
|
94 ub = +DBL_MAX; |
|
95 break; |
|
96 case LPX_UP: |
|
97 case LPX_DB: |
|
98 case LPX_FX: |
|
99 ub = lpx_get_col_ub(lp, j); |
|
100 break; |
|
101 default: |
|
102 xassert(lp != lp); |
|
103 } |
|
104 return ub; |
|
105 } |
|
106 |
|
107 static int is_binary(LPX *lp, int j) |
|
108 { /* this routine checks if variable x[j] is binary */ |
|
109 return |
|
110 lpx_get_col_kind(lp, j) == LPX_IV && |
|
111 lpx_get_col_type(lp, j) == LPX_DB && |
|
112 lpx_get_col_lb(lp, j) == 0.0 && lpx_get_col_ub(lp, j) == 1.0; |
|
113 } |
|
114 |
|
115 static double eval_lf_min(LPX *lp, int len, int ind[], double val[]) |
|
116 { /* this routine computes the minimum of a specified linear form |
|
117 |
|
118 sum a[j]*x[j] |
|
119 j |
|
120 |
|
121 using the formula: |
|
122 |
|
123 min = sum a[j]*lb[j] + sum a[j]*ub[j], |
|
124 j in J+ j in J- |
|
125 |
|
126 where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j] |
|
127 are lower and upper bound of variable x[j], resp. */ |
|
128 int j, t; |
|
129 double lb, ub, sum; |
|
130 sum = 0.0; |
|
131 for (t = 1; t <= len; t++) |
|
132 { j = ind[t]; |
|
133 if (val[t] > 0.0) |
|
134 { lb = get_col_lb(lp, j); |
|
135 if (lb == -DBL_MAX) |
|
136 { sum = -DBL_MAX; |
|
137 break; |
|
138 } |
|
139 sum += val[t] * lb; |
|
140 } |
|
141 else if (val[t] < 0.0) |
|
142 { ub = get_col_ub(lp, j); |
|
143 if (ub == +DBL_MAX) |
|
144 { sum = -DBL_MAX; |
|
145 break; |
|
146 } |
|
147 sum += val[t] * ub; |
|
148 } |
|
149 else |
|
150 xassert(val != val); |
|
151 } |
|
152 return sum; |
|
153 } |
|
154 |
|
155 static double eval_lf_max(LPX *lp, int len, int ind[], double val[]) |
|
156 { /* this routine computes the maximum of a specified linear form |
|
157 |
|
158 sum a[j]*x[j] |
|
159 j |
|
160 |
|
161 using the formula: |
|
162 |
|
163 max = sum a[j]*ub[j] + sum a[j]*lb[j], |
|
164 j in J+ j in J- |
|
165 |
|
166 where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j] |
|
167 are lower and upper bound of variable x[j], resp. */ |
|
168 int j, t; |
|
169 double lb, ub, sum; |
|
170 sum = 0.0; |
|
171 for (t = 1; t <= len; t++) |
|
172 { j = ind[t]; |
|
173 if (val[t] > 0.0) |
|
174 { ub = get_col_ub(lp, j); |
|
175 if (ub == +DBL_MAX) |
|
176 { sum = +DBL_MAX; |
|
177 break; |
|
178 } |
|
179 sum += val[t] * ub; |
|
180 } |
|
181 else if (val[t] < 0.0) |
|
182 { lb = get_col_lb(lp, j); |
|
183 if (lb == -DBL_MAX) |
|
184 { sum = +DBL_MAX; |
|
185 break; |
|
186 } |
|
187 sum += val[t] * lb; |
|
188 } |
|
189 else |
|
190 xassert(val != val); |
|
191 } |
|
192 return sum; |
|
193 } |
|
194 |
|
195 /*---------------------------------------------------------------------- |
|
196 -- probing - determine logical relation between binary variables. |
|
197 -- |
|
198 -- This routine tentatively sets a binary variable to 0 and then to 1 |
|
199 -- and examines whether another binary variable is caused to be fixed. |
|
200 -- |
|
201 -- The examination is based only on one row (constraint), which is the |
|
202 -- following: |
|
203 -- |
|
204 -- L <= sum a[j]*x[j] <= U. (1) |
|
205 -- j |
|
206 -- |
|
207 -- Let x[p] be a probing variable, x[q] be an examined variable. Then |
|
208 -- (1) can be written as: |
|
209 -- |
|
210 -- L <= sum a[j]*x[j] + a[p]*x[p] + a[q]*x[q] <= U, (2) |
|
211 -- j in J' |
|
212 -- |
|
213 -- where J' = {j: j != p and j != q}. |
|
214 -- |
|
215 -- Let |
|
216 -- |
|
217 -- L' = L - a[p]*x[p], (3) |
|
218 -- |
|
219 -- U' = U - a[p]*x[p], (4) |
|
220 -- |
|
221 -- where x[p] is assumed to be fixed at 0 or 1. So (2) can be rewritten |
|
222 -- as follows: |
|
223 -- |
|
224 -- L' <= sum a[j]*x[j] + a[q]*x[q] <= U', (5) |
|
225 -- j in J' |
|
226 -- |
|
227 -- from where we have: |
|
228 -- |
|
229 -- L' - sum a[j]*x[j] <= a[q]*x[q] <= U' - sum a[j]*x[j]. (6) |
|
230 -- j in J' j in J' |
|
231 -- |
|
232 -- Thus, |
|
233 -- |
|
234 -- min a[q]*x[q] = L' - MAX, (7) |
|
235 -- |
|
236 -- max a[q]*x[q] = U' - MIN, (8) |
|
237 -- |
|
238 -- where |
|
239 -- |
|
240 -- MIN = min sum a[j]*x[j], (9) |
|
241 -- j in J' |
|
242 -- |
|
243 -- MAX = max sum a[j]*x[j]. (10) |
|
244 -- j in J' |
|
245 -- |
|
246 -- Formulae (7) and (8) allows determining implied lower and upper |
|
247 -- bounds of x[q]. |
|
248 -- |
|
249 -- Parameters len, val, L and U specify the constraint (1). |
|
250 -- |
|
251 -- Parameters lf_min and lf_max specify implied lower and upper bounds |
|
252 -- of the linear form (1). It is assumed that these bounds are computed |
|
253 -- with the routines eval_lf_min and eval_lf_max (see above). |
|
254 -- |
|
255 -- Parameter p specifies the probing variable x[p], which is set to 0 |
|
256 -- (if set is 0) or to 1 (if set is 1). |
|
257 -- |
|
258 -- Parameter q specifies the examined variable x[q]. |
|
259 -- |
|
260 -- On exit the routine returns one of the following codes: |
|
261 -- |
|
262 -- 0 - there is no logical relation between x[p] and x[q]; |
|
263 -- 1 - x[q] can take only on value 0; |
|
264 -- 2 - x[q] can take only on value 1. */ |
|
265 |
|
266 static int probing(int len, double val[], double L, double U, |
|
267 double lf_min, double lf_max, int p, int set, int q) |
|
268 { double temp; |
|
269 xassert(1 <= p && p < q && q <= len); |
|
270 /* compute L' (3) */ |
|
271 if (L != -DBL_MAX && set) L -= val[p]; |
|
272 /* compute U' (4) */ |
|
273 if (U != +DBL_MAX && set) U -= val[p]; |
|
274 /* compute MIN (9) */ |
|
275 if (lf_min != -DBL_MAX) |
|
276 { if (val[p] < 0.0) lf_min -= val[p]; |
|
277 if (val[q] < 0.0) lf_min -= val[q]; |
|
278 } |
|
279 /* compute MAX (10) */ |
|
280 if (lf_max != +DBL_MAX) |
|
281 { if (val[p] > 0.0) lf_max -= val[p]; |
|
282 if (val[q] > 0.0) lf_max -= val[q]; |
|
283 } |
|
284 /* compute implied lower bound of x[q]; see (7), (8) */ |
|
285 if (val[q] > 0.0) |
|
286 { if (L == -DBL_MAX || lf_max == +DBL_MAX) |
|
287 temp = -DBL_MAX; |
|
288 else |
|
289 temp = (L - lf_max) / val[q]; |
|
290 } |
|
291 else |
|
292 { if (U == +DBL_MAX || lf_min == -DBL_MAX) |
|
293 temp = -DBL_MAX; |
|
294 else |
|
295 temp = (U - lf_min) / val[q]; |
|
296 } |
|
297 if (temp > 0.001) return 2; |
|
298 /* compute implied upper bound of x[q]; see (7), (8) */ |
|
299 if (val[q] > 0.0) |
|
300 { if (U == +DBL_MAX || lf_min == -DBL_MAX) |
|
301 temp = +DBL_MAX; |
|
302 else |
|
303 temp = (U - lf_min) / val[q]; |
|
304 } |
|
305 else |
|
306 { if (L == -DBL_MAX || lf_max == +DBL_MAX) |
|
307 temp = +DBL_MAX; |
|
308 else |
|
309 temp = (L - lf_max) / val[q]; |
|
310 } |
|
311 if (temp < 0.999) return 1; |
|
312 /* there is no logical relation between x[p] and x[q] */ |
|
313 return 0; |
|
314 } |
|
315 |
|
316 struct COG |
|
317 { /* conflict graph; it represents logical relations between binary |
|
318 variables and has a vertex for each binary variable and its |
|
319 complement, and an edge between two vertices when at most one |
|
320 of the variables represented by the vertices can equal one in |
|
321 an optimal solution */ |
|
322 int n; |
|
323 /* number of variables */ |
|
324 int nb; |
|
325 /* number of binary variables represented in the graph (note that |
|
326 not all binary variables can be represented); vertices which |
|
327 correspond to binary variables have numbers 1, ..., nb while |
|
328 vertices which correspond to complements of binary variables |
|
329 have numbers nb+1, ..., nb+nb */ |
|
330 int ne; |
|
331 /* number of edges in the graph */ |
|
332 int *vert; /* int vert[1+n]; */ |
|
333 /* if x[j] is a binary variable represented in the graph, vert[j] |
|
334 is the vertex number corresponding to x[j]; otherwise vert[j] |
|
335 is zero */ |
|
336 int *orig; /* int list[1:nb]; */ |
|
337 /* if vert[j] = k > 0, then orig[k] = j */ |
|
338 unsigned char *a; |
|
339 /* adjacency matrix of the graph having 2*nb rows and columns; |
|
340 only strict lower triangle is stored in dense packed form */ |
|
341 }; |
|
342 |
|
343 /*---------------------------------------------------------------------- |
|
344 -- lpx_create_cog - create the conflict graph. |
|
345 -- |
|
346 -- SYNOPSIS |
|
347 -- |
|
348 -- #include "glplpx.h" |
|
349 -- void *lpx_create_cog(LPX *lp); |
|
350 -- |
|
351 -- DESCRIPTION |
|
352 -- |
|
353 -- The routine lpx_create_cog creates the conflict graph for a given |
|
354 -- problem instance. |
|
355 -- |
|
356 -- RETURNS |
|
357 -- |
|
358 -- If the graph has been created, the routine returns a pointer to it. |
|
359 -- Otherwise the routine returns NULL. */ |
|
360 |
|
361 #define MAX_NB 4000 |
|
362 #define MAX_ROW_LEN 500 |
|
363 |
|
364 static void lpx_add_cog_edge(void *_cog, int i, int j); |
|
365 |
|
366 static void *lpx_create_cog(LPX *lp) |
|
367 { struct COG *cog = NULL; |
|
368 int m, n, nb, i, j, p, q, len, *ind, *vert, *orig; |
|
369 double L, U, lf_min, lf_max, *val; |
|
370 xprintf("Creating the conflict graph...\n"); |
|
371 m = lpx_get_num_rows(lp); |
|
372 n = lpx_get_num_cols(lp); |
|
373 /* determine which binary variables should be included in the |
|
374 conflict graph */ |
|
375 nb = 0; |
|
376 vert = xcalloc(1+n, sizeof(int)); |
|
377 for (j = 1; j <= n; j++) vert[j] = 0; |
|
378 orig = xcalloc(1+n, sizeof(int)); |
|
379 ind = xcalloc(1+n, sizeof(int)); |
|
380 val = xcalloc(1+n, sizeof(double)); |
|
381 for (i = 1; i <= m; i++) |
|
382 { L = get_row_lb(lp, i); |
|
383 U = get_row_ub(lp, i); |
|
384 if (L == -DBL_MAX && U == +DBL_MAX) continue; |
|
385 len = lpx_get_mat_row(lp, i, ind, val); |
|
386 if (len > MAX_ROW_LEN) continue; |
|
387 lf_min = eval_lf_min(lp, len, ind, val); |
|
388 lf_max = eval_lf_max(lp, len, ind, val); |
|
389 for (p = 1; p <= len; p++) |
|
390 { if (!is_binary(lp, ind[p])) continue; |
|
391 for (q = p+1; q <= len; q++) |
|
392 { if (!is_binary(lp, ind[q])) continue; |
|
393 if (probing(len, val, L, U, lf_min, lf_max, p, 0, q) || |
|
394 probing(len, val, L, U, lf_min, lf_max, p, 1, q)) |
|
395 { /* there is a logical relation */ |
|
396 /* include the first variable in the graph */ |
|
397 j = ind[p]; |
|
398 if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j; |
|
399 /* incude the second variable in the graph */ |
|
400 j = ind[q]; |
|
401 if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j; |
|
402 } |
|
403 } |
|
404 } |
|
405 } |
|
406 /* if the graph is either empty or has too many vertices, do not |
|
407 create it */ |
|
408 if (nb == 0 || nb > MAX_NB) |
|
409 { xprintf("The conflict graph is either empty or too big\n"); |
|
410 xfree(vert); |
|
411 xfree(orig); |
|
412 goto done; |
|
413 } |
|
414 /* create the conflict graph */ |
|
415 cog = xmalloc(sizeof(struct COG)); |
|
416 cog->n = n; |
|
417 cog->nb = nb; |
|
418 cog->ne = 0; |
|
419 cog->vert = vert; |
|
420 cog->orig = orig; |
|
421 len = nb + nb; /* number of vertices */ |
|
422 len = (len * (len - 1)) / 2; /* number of entries in triangle */ |
|
423 len = (len + (CHAR_BIT - 1)) / CHAR_BIT; /* bytes needed */ |
|
424 cog->a = xmalloc(len); |
|
425 memset(cog->a, 0, len); |
|
426 for (j = 1; j <= nb; j++) |
|
427 { /* add edge between variable and its complement */ |
|
428 lpx_add_cog_edge(cog, +orig[j], -orig[j]); |
|
429 } |
|
430 for (i = 1; i <= m; i++) |
|
431 { L = get_row_lb(lp, i); |
|
432 U = get_row_ub(lp, i); |
|
433 if (L == -DBL_MAX && U == +DBL_MAX) continue; |
|
434 len = lpx_get_mat_row(lp, i, ind, val); |
|
435 if (len > MAX_ROW_LEN) continue; |
|
436 lf_min = eval_lf_min(lp, len, ind, val); |
|
437 lf_max = eval_lf_max(lp, len, ind, val); |
|
438 for (p = 1; p <= len; p++) |
|
439 { if (!is_binary(lp, ind[p])) continue; |
|
440 for (q = p+1; q <= len; q++) |
|
441 { if (!is_binary(lp, ind[q])) continue; |
|
442 /* set x[p] to 0 and examine x[q] */ |
|
443 switch (probing(len, val, L, U, lf_min, lf_max, p, 0, q)) |
|
444 { case 0: |
|
445 /* no logical relation */ |
|
446 break; |
|
447 case 1: |
|
448 /* x[p] = 0 implies x[q] = 0 */ |
|
449 lpx_add_cog_edge(cog, -ind[p], +ind[q]); |
|
450 break; |
|
451 case 2: |
|
452 /* x[p] = 0 implies x[q] = 1 */ |
|
453 lpx_add_cog_edge(cog, -ind[p], -ind[q]); |
|
454 break; |
|
455 default: |
|
456 xassert(lp != lp); |
|
457 } |
|
458 /* set x[p] to 1 and examine x[q] */ |
|
459 switch (probing(len, val, L, U, lf_min, lf_max, p, 1, q)) |
|
460 { case 0: |
|
461 /* no logical relation */ |
|
462 break; |
|
463 case 1: |
|
464 /* x[p] = 1 implies x[q] = 0 */ |
|
465 lpx_add_cog_edge(cog, +ind[p], +ind[q]); |
|
466 break; |
|
467 case 2: |
|
468 /* x[p] = 1 implies x[q] = 1 */ |
|
469 lpx_add_cog_edge(cog, +ind[p], -ind[q]); |
|
470 break; |
|
471 default: |
|
472 xassert(lp != lp); |
|
473 } |
|
474 } |
|
475 } |
|
476 } |
|
477 xprintf("The conflict graph has 2*%d vertices and %d edges\n", |
|
478 cog->nb, cog->ne); |
|
479 done: xfree(ind); |
|
480 xfree(val); |
|
481 return cog; |
|
482 } |
|
483 |
|
484 /*---------------------------------------------------------------------- |
|
485 -- lpx_add_cog_edge - add edge to the conflict graph. |
|
486 -- |
|
487 -- SYNOPSIS |
|
488 -- |
|
489 -- #include "glplpx.h" |
|
490 -- void lpx_add_cog_edge(void *cog, int i, int j); |
|
491 -- |
|
492 -- DESCRIPTION |
|
493 -- |
|
494 -- The routine lpx_add_cog_edge adds an edge to the conflict graph. |
|
495 -- The edge connects x[i] (if i > 0) or its complement (if i < 0) and |
|
496 -- x[j] (if j > 0) or its complement (if j < 0), where i and j are |
|
497 -- original ordinal numbers of corresponding variables. */ |
|
498 |
|
499 static void lpx_add_cog_edge(void *_cog, int i, int j) |
|
500 { struct COG *cog = _cog; |
|
501 int k; |
|
502 xassert(i != j); |
|
503 /* determine indices of corresponding vertices */ |
|
504 if (i > 0) |
|
505 { xassert(1 <= i && i <= cog->n); |
|
506 i = cog->vert[i]; |
|
507 xassert(i != 0); |
|
508 } |
|
509 else |
|
510 { i = -i; |
|
511 xassert(1 <= i && i <= cog->n); |
|
512 i = cog->vert[i]; |
|
513 xassert(i != 0); |
|
514 i += cog->nb; |
|
515 } |
|
516 if (j > 0) |
|
517 { xassert(1 <= j && j <= cog->n); |
|
518 j = cog->vert[j]; |
|
519 xassert(j != 0); |
|
520 } |
|
521 else |
|
522 { j = -j; |
|
523 xassert(1 <= j && j <= cog->n); |
|
524 j = cog->vert[j]; |
|
525 xassert(j != 0); |
|
526 j += cog->nb; |
|
527 } |
|
528 /* only lower triangle is stored, so we need i > j */ |
|
529 if (i < j) k = i, i = j, j = k; |
|
530 k = ((i - 1) * (i - 2)) / 2 + (j - 1); |
|
531 cog->a[k / CHAR_BIT] |= |
|
532 (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT)); |
|
533 cog->ne++; |
|
534 return; |
|
535 } |
|
536 |
|
537 /*---------------------------------------------------------------------- |
|
538 -- MAXIMUM WEIGHT CLIQUE |
|
539 -- |
|
540 -- Two subroutines sub() and wclique() below are intended to find a |
|
541 -- maximum weight clique in a given undirected graph. These subroutines |
|
542 -- are slightly modified version of the program WCLIQUE developed by |
|
543 -- Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based |
|
544 -- on ideas from the article "P. R. J. Ostergard, A new algorithm for |
|
545 -- the maximum-weight clique problem, submitted for publication", which |
|
546 -- in turn is a generalization of the algorithm for unweighted graphs |
|
547 -- presented in "P. R. J. Ostergard, A fast algorithm for the maximum |
|
548 -- clique problem, submitted for publication". |
|
549 -- |
|
550 -- USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE. */ |
|
551 |
|
552 struct dsa |
|
553 { /* dynamic storage area */ |
|
554 int n; |
|
555 /* number of vertices */ |
|
556 int *wt; /* int wt[0:n-1]; */ |
|
557 /* weights */ |
|
558 unsigned char *a; |
|
559 /* adjacency matrix (packed lower triangle without main diag.) */ |
|
560 int record; |
|
561 /* weight of best clique */ |
|
562 int rec_level; |
|
563 /* number of vertices in best clique */ |
|
564 int *rec; /* int rec[0:n-1]; */ |
|
565 /* best clique so far */ |
|
566 int *clique; /* int clique[0:n-1]; */ |
|
567 /* table for pruning */ |
|
568 int *set; /* int set[0:n-1]; */ |
|
569 /* current clique */ |
|
570 }; |
|
571 |
|
572 #define n (dsa->n) |
|
573 #define wt (dsa->wt) |
|
574 #define a (dsa->a) |
|
575 #define record (dsa->record) |
|
576 #define rec_level (dsa->rec_level) |
|
577 #define rec (dsa->rec) |
|
578 #define clique (dsa->clique) |
|
579 #define set (dsa->set) |
|
580 |
|
581 #if 0 |
|
582 static int is_edge(struct dsa *dsa, int i, int j) |
|
583 { /* if there is arc (i,j), the routine returns true; otherwise |
|
584 false; 0 <= i, j < n */ |
|
585 int k; |
|
586 xassert(0 <= i && i < n); |
|
587 xassert(0 <= j && j < n); |
|
588 if (i == j) return 0; |
|
589 if (i < j) k = i, i = j, j = k; |
|
590 k = (i * (i - 1)) / 2 + j; |
|
591 return a[k / CHAR_BIT] & |
|
592 (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT)); |
|
593 } |
|
594 #else |
|
595 #define is_edge(dsa, i, j) ((i) == (j) ? 0 : \ |
|
596 (i) > (j) ? is_edge1(i, j) : is_edge1(j, i)) |
|
597 #define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j)) |
|
598 #define is_edge2(k) (a[(k) / CHAR_BIT] & \ |
|
599 (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT))) |
|
600 #endif |
|
601 |
|
602 static void sub(struct dsa *dsa, int ct, int table[], int level, |
|
603 int weight, int l_weight) |
|
604 { int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable; |
|
605 newtable = xcalloc(n, sizeof(int)); |
|
606 if (ct <= 0) |
|
607 { /* 0 or 1 elements left; include these */ |
|
608 if (ct == 0) |
|
609 { set[level++] = table[0]; |
|
610 weight += l_weight; |
|
611 } |
|
612 if (weight > record) |
|
613 { record = weight; |
|
614 rec_level = level; |
|
615 for (i = 0; i < level; i++) rec[i] = set[i]; |
|
616 } |
|
617 goto done; |
|
618 } |
|
619 for (i = ct; i >= 0; i--) |
|
620 { if ((level == 0) && (i < ct)) goto done; |
|
621 k = table[i]; |
|
622 if ((level > 0) && (clique[k] <= (record - weight))) |
|
623 goto done; /* prune */ |
|
624 set[level] = k; |
|
625 curr_weight = weight + wt[k]; |
|
626 l_weight -= wt[k]; |
|
627 if (l_weight <= (record - curr_weight)) |
|
628 goto done; /* prune */ |
|
629 p1 = newtable; |
|
630 p2 = table; |
|
631 left_weight = 0; |
|
632 while (p2 < table + i) |
|
633 { j = *p2++; |
|
634 if (is_edge(dsa, j, k)) |
|
635 { *p1++ = j; |
|
636 left_weight += wt[j]; |
|
637 } |
|
638 } |
|
639 if (left_weight <= (record - curr_weight)) continue; |
|
640 sub(dsa, p1 - newtable - 1, newtable, level + 1, curr_weight, |
|
641 left_weight); |
|
642 } |
|
643 done: xfree(newtable); |
|
644 return; |
|
645 } |
|
646 |
|
647 static int wclique(int _n, int w[], unsigned char _a[], int sol[]) |
|
648 { struct dsa _dsa, *dsa = &_dsa; |
|
649 int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos; |
|
650 glp_long timer; |
|
651 n = _n; |
|
652 wt = &w[1]; |
|
653 a = _a; |
|
654 record = 0; |
|
655 rec_level = 0; |
|
656 rec = &sol[1]; |
|
657 clique = xcalloc(n, sizeof(int)); |
|
658 set = xcalloc(n, sizeof(int)); |
|
659 used = xcalloc(n, sizeof(int)); |
|
660 nwt = xcalloc(n, sizeof(int)); |
|
661 pos = xcalloc(n, sizeof(int)); |
|
662 /* start timer */ |
|
663 timer = xtime(); |
|
664 /* order vertices */ |
|
665 for (i = 0; i < n; i++) |
|
666 { nwt[i] = 0; |
|
667 for (j = 0; j < n; j++) |
|
668 if (is_edge(dsa, i, j)) nwt[i] += wt[j]; |
|
669 } |
|
670 for (i = 0; i < n; i++) |
|
671 used[i] = 0; |
|
672 for (i = n-1; i >= 0; i--) |
|
673 { max_wt = -1; |
|
674 max_nwt = -1; |
|
675 for (j = 0; j < n; j++) |
|
676 { if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt |
|
677 && nwt[j] > max_nwt))) |
|
678 { max_wt = wt[j]; |
|
679 max_nwt = nwt[j]; |
|
680 p = j; |
|
681 } |
|
682 } |
|
683 pos[i] = p; |
|
684 used[p] = 1; |
|
685 for (j = 0; j < n; j++) |
|
686 if ((!used[j]) && (j != p) && (is_edge(dsa, p, j))) |
|
687 nwt[j] -= wt[p]; |
|
688 } |
|
689 /* main routine */ |
|
690 wth = 0; |
|
691 for (i = 0; i < n; i++) |
|
692 { wth += wt[pos[i]]; |
|
693 sub(dsa, i, pos, 0, 0, wth); |
|
694 clique[pos[i]] = record; |
|
695 #if 0 |
|
696 if (utime() >= timer + 5.0) |
|
697 #else |
|
698 if (xdifftime(xtime(), timer) >= 5.0 - 0.001) |
|
699 #endif |
|
700 { /* print current record and reset timer */ |
|
701 xprintf("level = %d (%d); best = %d\n", i+1, n, record); |
|
702 #if 0 |
|
703 timer = utime(); |
|
704 #else |
|
705 timer = xtime(); |
|
706 #endif |
|
707 } |
|
708 } |
|
709 xfree(clique); |
|
710 xfree(set); |
|
711 xfree(used); |
|
712 xfree(nwt); |
|
713 xfree(pos); |
|
714 /* return the solution found */ |
|
715 for (i = 1; i <= rec_level; i++) sol[i]++; |
|
716 return rec_level; |
|
717 } |
|
718 |
|
719 #undef n |
|
720 #undef wt |
|
721 #undef a |
|
722 #undef record |
|
723 #undef rec_level |
|
724 #undef rec |
|
725 #undef clique |
|
726 #undef set |
|
727 |
|
728 /*---------------------------------------------------------------------- |
|
729 -- lpx_clique_cut - generate cluque cut. |
|
730 -- |
|
731 -- SYNOPSIS |
|
732 -- |
|
733 -- #include "glplpx.h" |
|
734 -- int lpx_clique_cut(LPX *lp, void *cog, int ind[], double val[]); |
|
735 -- |
|
736 -- DESCRIPTION |
|
737 -- |
|
738 -- The routine lpx_clique_cut generates a clique cut using the conflict |
|
739 -- graph specified by the parameter cog. |
|
740 -- |
|
741 -- If a violated clique cut has been found, it has the following form: |
|
742 -- |
|
743 -- sum{j in J} a[j]*x[j] <= b. |
|
744 -- |
|
745 -- Variable indices j in J are stored in elements ind[1], ..., ind[len] |
|
746 -- while corresponding constraint coefficients are stored in elements |
|
747 -- val[1], ..., val[len], where len is returned on exit. The right-hand |
|
748 -- side b is stored in element val[0]. |
|
749 -- |
|
750 -- RETURNS |
|
751 -- |
|
752 -- If the cutting plane has been successfully generated, the routine |
|
753 -- returns 1 <= len <= n, which is the number of non-zero coefficients |
|
754 -- in the inequality constraint. Otherwise, the routine returns zero. */ |
|
755 |
|
756 static int lpx_clique_cut(LPX *lp, void *_cog, int ind[], double val[]) |
|
757 { struct COG *cog = _cog; |
|
758 int n = lpx_get_num_cols(lp); |
|
759 int j, t, v, card, temp, len = 0, *w, *sol; |
|
760 double x, sum, b, *vec; |
|
761 /* allocate working arrays */ |
|
762 w = xcalloc(1 + 2 * cog->nb, sizeof(int)); |
|
763 sol = xcalloc(1 + 2 * cog->nb, sizeof(int)); |
|
764 vec = xcalloc(1+n, sizeof(double)); |
|
765 /* assign weights to vertices of the conflict graph */ |
|
766 for (t = 1; t <= cog->nb; t++) |
|
767 { j = cog->orig[t]; |
|
768 x = lpx_get_col_prim(lp, j); |
|
769 temp = (int)(100.0 * x + 0.5); |
|
770 if (temp < 0) temp = 0; |
|
771 if (temp > 100) temp = 100; |
|
772 w[t] = temp; |
|
773 w[cog->nb + t] = 100 - temp; |
|
774 } |
|
775 /* find a clique of maximum weight */ |
|
776 card = wclique(2 * cog->nb, w, cog->a, sol); |
|
777 /* compute the clique weight for unscaled values */ |
|
778 sum = 0.0; |
|
779 for ( t = 1; t <= card; t++) |
|
780 { v = sol[t]; |
|
781 xassert(1 <= v && v <= 2 * cog->nb); |
|
782 if (v <= cog->nb) |
|
783 { /* vertex v corresponds to binary variable x[j] */ |
|
784 j = cog->orig[v]; |
|
785 x = lpx_get_col_prim(lp, j); |
|
786 sum += x; |
|
787 } |
|
788 else |
|
789 { /* vertex v corresponds to the complement of x[j] */ |
|
790 j = cog->orig[v - cog->nb]; |
|
791 x = lpx_get_col_prim(lp, j); |
|
792 sum += 1.0 - x; |
|
793 } |
|
794 } |
|
795 /* if the sum of binary variables and their complements in the |
|
796 clique greater than 1, the clique cut is violated */ |
|
797 if (sum >= 1.01) |
|
798 { /* construct the inquality */ |
|
799 for (j = 1; j <= n; j++) vec[j] = 0; |
|
800 b = 1.0; |
|
801 for (t = 1; t <= card; t++) |
|
802 { v = sol[t]; |
|
803 if (v <= cog->nb) |
|
804 { /* vertex v corresponds to binary variable x[j] */ |
|
805 j = cog->orig[v]; |
|
806 xassert(1 <= j && j <= n); |
|
807 vec[j] += 1.0; |
|
808 } |
|
809 else |
|
810 { /* vertex v corresponds to the complement of x[j] */ |
|
811 j = cog->orig[v - cog->nb]; |
|
812 xassert(1 <= j && j <= n); |
|
813 vec[j] -= 1.0; |
|
814 b -= 1.0; |
|
815 } |
|
816 } |
|
817 xassert(len == 0); |
|
818 for (j = 1; j <= n; j++) |
|
819 { if (vec[j] != 0.0) |
|
820 { len++; |
|
821 ind[len] = j, val[len] = vec[j]; |
|
822 } |
|
823 } |
|
824 ind[0] = 0, val[0] = b; |
|
825 } |
|
826 /* free working arrays */ |
|
827 xfree(w); |
|
828 xfree(sol); |
|
829 xfree(vec); |
|
830 /* return to the calling program */ |
|
831 return len; |
|
832 } |
|
833 |
|
834 /*---------------------------------------------------------------------- |
|
835 -- lpx_delete_cog - delete the conflict graph. |
|
836 -- |
|
837 -- SYNOPSIS |
|
838 -- |
|
839 -- #include "glplpx.h" |
|
840 -- void lpx_delete_cog(void *cog); |
|
841 -- |
|
842 -- DESCRIPTION |
|
843 -- |
|
844 -- The routine lpx_delete_cog deletes the conflict graph, which the |
|
845 -- parameter cog points to, freeing all the memory allocated to this |
|
846 -- object. */ |
|
847 |
|
848 static void lpx_delete_cog(void *_cog) |
|
849 { struct COG *cog = _cog; |
|
850 xfree(cog->vert); |
|
851 xfree(cog->orig); |
|
852 xfree(cog->a); |
|
853 xfree(cog); |
|
854 } |
|
855 |
|
856 /**********************************************************************/ |
|
857 |
|
858 void *ios_clq_init(glp_tree *tree) |
|
859 { /* initialize clique cut generator */ |
|
860 glp_prob *mip = tree->mip; |
|
861 xassert(mip != NULL); |
|
862 return lpx_create_cog(mip); |
|
863 } |
|
864 |
|
865 /*********************************************************************** |
|
866 * NAME |
|
867 * |
|
868 * ios_clq_gen - generate clique cuts |
|
869 * |
|
870 * SYNOPSIS |
|
871 * |
|
872 * #include "glpios.h" |
|
873 * void ios_clq_gen(glp_tree *tree, void *gen); |
|
874 * |
|
875 * DESCRIPTION |
|
876 * |
|
877 * The routine ios_clq_gen generates clique cuts for the current point |
|
878 * and adds them to the clique pool. */ |
|
879 |
|
880 void ios_clq_gen(glp_tree *tree, void *gen) |
|
881 { int n = lpx_get_num_cols(tree->mip); |
|
882 int len, *ind; |
|
883 double *val; |
|
884 xassert(gen != NULL); |
|
885 ind = xcalloc(1+n, sizeof(int)); |
|
886 val = xcalloc(1+n, sizeof(double)); |
|
887 len = lpx_clique_cut(tree->mip, gen, ind, val); |
|
888 if (len > 0) |
|
889 { /* xprintf("len = %d\n", len); */ |
|
890 glp_ios_add_row(tree, NULL, GLP_RF_CLQ, 0, len, ind, val, |
|
891 GLP_UP, val[0]); |
|
892 } |
|
893 xfree(ind); |
|
894 xfree(val); |
|
895 return; |
|
896 } |
|
897 |
|
898 /**********************************************************************/ |
|
899 |
|
900 void ios_clq_term(void *gen) |
|
901 { /* terminate clique cut generator */ |
|
902 xassert(gen != NULL); |
|
903 lpx_delete_cog(gen); |
|
904 return; |
|
905 } |
|
906 |
|
907 /* eof */ |