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 */ |
---|