1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpios08.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,907 @@
1.4 +/* glpios08.c (clique cut generator) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpios.h"
1.29 +
1.30 +static double get_row_lb(LPX *lp, int i)
1.31 +{ /* this routine returns lower bound of row i or -DBL_MAX if the
1.32 + row has no lower bound */
1.33 + double lb;
1.34 + switch (lpx_get_row_type(lp, i))
1.35 + { case LPX_FR:
1.36 + case LPX_UP:
1.37 + lb = -DBL_MAX;
1.38 + break;
1.39 + case LPX_LO:
1.40 + case LPX_DB:
1.41 + case LPX_FX:
1.42 + lb = lpx_get_row_lb(lp, i);
1.43 + break;
1.44 + default:
1.45 + xassert(lp != lp);
1.46 + }
1.47 + return lb;
1.48 +}
1.49 +
1.50 +static double get_row_ub(LPX *lp, int i)
1.51 +{ /* this routine returns upper bound of row i or +DBL_MAX if the
1.52 + row has no upper bound */
1.53 + double ub;
1.54 + switch (lpx_get_row_type(lp, i))
1.55 + { case LPX_FR:
1.56 + case LPX_LO:
1.57 + ub = +DBL_MAX;
1.58 + break;
1.59 + case LPX_UP:
1.60 + case LPX_DB:
1.61 + case LPX_FX:
1.62 + ub = lpx_get_row_ub(lp, i);
1.63 + break;
1.64 + default:
1.65 + xassert(lp != lp);
1.66 + }
1.67 + return ub;
1.68 +}
1.69 +
1.70 +static double get_col_lb(LPX *lp, int j)
1.71 +{ /* this routine returns lower bound of column j or -DBL_MAX if
1.72 + the column has no lower bound */
1.73 + double lb;
1.74 + switch (lpx_get_col_type(lp, j))
1.75 + { case LPX_FR:
1.76 + case LPX_UP:
1.77 + lb = -DBL_MAX;
1.78 + break;
1.79 + case LPX_LO:
1.80 + case LPX_DB:
1.81 + case LPX_FX:
1.82 + lb = lpx_get_col_lb(lp, j);
1.83 + break;
1.84 + default:
1.85 + xassert(lp != lp);
1.86 + }
1.87 + return lb;
1.88 +}
1.89 +
1.90 +static double get_col_ub(LPX *lp, int j)
1.91 +{ /* this routine returns upper bound of column j or +DBL_MAX if
1.92 + the column has no upper bound */
1.93 + double ub;
1.94 + switch (lpx_get_col_type(lp, j))
1.95 + { case LPX_FR:
1.96 + case LPX_LO:
1.97 + ub = +DBL_MAX;
1.98 + break;
1.99 + case LPX_UP:
1.100 + case LPX_DB:
1.101 + case LPX_FX:
1.102 + ub = lpx_get_col_ub(lp, j);
1.103 + break;
1.104 + default:
1.105 + xassert(lp != lp);
1.106 + }
1.107 + return ub;
1.108 +}
1.109 +
1.110 +static int is_binary(LPX *lp, int j)
1.111 +{ /* this routine checks if variable x[j] is binary */
1.112 + return
1.113 + lpx_get_col_kind(lp, j) == LPX_IV &&
1.114 + lpx_get_col_type(lp, j) == LPX_DB &&
1.115 + lpx_get_col_lb(lp, j) == 0.0 && lpx_get_col_ub(lp, j) == 1.0;
1.116 +}
1.117 +
1.118 +static double eval_lf_min(LPX *lp, int len, int ind[], double val[])
1.119 +{ /* this routine computes the minimum of a specified linear form
1.120 +
1.121 + sum a[j]*x[j]
1.122 + j
1.123 +
1.124 + using the formula:
1.125 +
1.126 + min = sum a[j]*lb[j] + sum a[j]*ub[j],
1.127 + j in J+ j in J-
1.128 +
1.129 + where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j]
1.130 + are lower and upper bound of variable x[j], resp. */
1.131 + int j, t;
1.132 + double lb, ub, sum;
1.133 + sum = 0.0;
1.134 + for (t = 1; t <= len; t++)
1.135 + { j = ind[t];
1.136 + if (val[t] > 0.0)
1.137 + { lb = get_col_lb(lp, j);
1.138 + if (lb == -DBL_MAX)
1.139 + { sum = -DBL_MAX;
1.140 + break;
1.141 + }
1.142 + sum += val[t] * lb;
1.143 + }
1.144 + else if (val[t] < 0.0)
1.145 + { ub = get_col_ub(lp, j);
1.146 + if (ub == +DBL_MAX)
1.147 + { sum = -DBL_MAX;
1.148 + break;
1.149 + }
1.150 + sum += val[t] * ub;
1.151 + }
1.152 + else
1.153 + xassert(val != val);
1.154 + }
1.155 + return sum;
1.156 +}
1.157 +
1.158 +static double eval_lf_max(LPX *lp, int len, int ind[], double val[])
1.159 +{ /* this routine computes the maximum of a specified linear form
1.160 +
1.161 + sum a[j]*x[j]
1.162 + j
1.163 +
1.164 + using the formula:
1.165 +
1.166 + max = sum a[j]*ub[j] + sum a[j]*lb[j],
1.167 + j in J+ j in J-
1.168 +
1.169 + where J+ = {j: a[j] > 0}, J- = {j: a[j] < 0}, lb[j] and ub[j]
1.170 + are lower and upper bound of variable x[j], resp. */
1.171 + int j, t;
1.172 + double lb, ub, sum;
1.173 + sum = 0.0;
1.174 + for (t = 1; t <= len; t++)
1.175 + { j = ind[t];
1.176 + if (val[t] > 0.0)
1.177 + { ub = get_col_ub(lp, j);
1.178 + if (ub == +DBL_MAX)
1.179 + { sum = +DBL_MAX;
1.180 + break;
1.181 + }
1.182 + sum += val[t] * ub;
1.183 + }
1.184 + else if (val[t] < 0.0)
1.185 + { lb = get_col_lb(lp, j);
1.186 + if (lb == -DBL_MAX)
1.187 + { sum = +DBL_MAX;
1.188 + break;
1.189 + }
1.190 + sum += val[t] * lb;
1.191 + }
1.192 + else
1.193 + xassert(val != val);
1.194 + }
1.195 + return sum;
1.196 +}
1.197 +
1.198 +/*----------------------------------------------------------------------
1.199 +-- probing - determine logical relation between binary variables.
1.200 +--
1.201 +-- This routine tentatively sets a binary variable to 0 and then to 1
1.202 +-- and examines whether another binary variable is caused to be fixed.
1.203 +--
1.204 +-- The examination is based only on one row (constraint), which is the
1.205 +-- following:
1.206 +--
1.207 +-- L <= sum a[j]*x[j] <= U. (1)
1.208 +-- j
1.209 +--
1.210 +-- Let x[p] be a probing variable, x[q] be an examined variable. Then
1.211 +-- (1) can be written as:
1.212 +--
1.213 +-- L <= sum a[j]*x[j] + a[p]*x[p] + a[q]*x[q] <= U, (2)
1.214 +-- j in J'
1.215 +--
1.216 +-- where J' = {j: j != p and j != q}.
1.217 +--
1.218 +-- Let
1.219 +--
1.220 +-- L' = L - a[p]*x[p], (3)
1.221 +--
1.222 +-- U' = U - a[p]*x[p], (4)
1.223 +--
1.224 +-- where x[p] is assumed to be fixed at 0 or 1. So (2) can be rewritten
1.225 +-- as follows:
1.226 +--
1.227 +-- L' <= sum a[j]*x[j] + a[q]*x[q] <= U', (5)
1.228 +-- j in J'
1.229 +--
1.230 +-- from where we have:
1.231 +--
1.232 +-- L' - sum a[j]*x[j] <= a[q]*x[q] <= U' - sum a[j]*x[j]. (6)
1.233 +-- j in J' j in J'
1.234 +--
1.235 +-- Thus,
1.236 +--
1.237 +-- min a[q]*x[q] = L' - MAX, (7)
1.238 +--
1.239 +-- max a[q]*x[q] = U' - MIN, (8)
1.240 +--
1.241 +-- where
1.242 +--
1.243 +-- MIN = min sum a[j]*x[j], (9)
1.244 +-- j in J'
1.245 +--
1.246 +-- MAX = max sum a[j]*x[j]. (10)
1.247 +-- j in J'
1.248 +--
1.249 +-- Formulae (7) and (8) allows determining implied lower and upper
1.250 +-- bounds of x[q].
1.251 +--
1.252 +-- Parameters len, val, L and U specify the constraint (1).
1.253 +--
1.254 +-- Parameters lf_min and lf_max specify implied lower and upper bounds
1.255 +-- of the linear form (1). It is assumed that these bounds are computed
1.256 +-- with the routines eval_lf_min and eval_lf_max (see above).
1.257 +--
1.258 +-- Parameter p specifies the probing variable x[p], which is set to 0
1.259 +-- (if set is 0) or to 1 (if set is 1).
1.260 +--
1.261 +-- Parameter q specifies the examined variable x[q].
1.262 +--
1.263 +-- On exit the routine returns one of the following codes:
1.264 +--
1.265 +-- 0 - there is no logical relation between x[p] and x[q];
1.266 +-- 1 - x[q] can take only on value 0;
1.267 +-- 2 - x[q] can take only on value 1. */
1.268 +
1.269 +static int probing(int len, double val[], double L, double U,
1.270 + double lf_min, double lf_max, int p, int set, int q)
1.271 +{ double temp;
1.272 + xassert(1 <= p && p < q && q <= len);
1.273 + /* compute L' (3) */
1.274 + if (L != -DBL_MAX && set) L -= val[p];
1.275 + /* compute U' (4) */
1.276 + if (U != +DBL_MAX && set) U -= val[p];
1.277 + /* compute MIN (9) */
1.278 + if (lf_min != -DBL_MAX)
1.279 + { if (val[p] < 0.0) lf_min -= val[p];
1.280 + if (val[q] < 0.0) lf_min -= val[q];
1.281 + }
1.282 + /* compute MAX (10) */
1.283 + if (lf_max != +DBL_MAX)
1.284 + { if (val[p] > 0.0) lf_max -= val[p];
1.285 + if (val[q] > 0.0) lf_max -= val[q];
1.286 + }
1.287 + /* compute implied lower bound of x[q]; see (7), (8) */
1.288 + if (val[q] > 0.0)
1.289 + { if (L == -DBL_MAX || lf_max == +DBL_MAX)
1.290 + temp = -DBL_MAX;
1.291 + else
1.292 + temp = (L - lf_max) / val[q];
1.293 + }
1.294 + else
1.295 + { if (U == +DBL_MAX || lf_min == -DBL_MAX)
1.296 + temp = -DBL_MAX;
1.297 + else
1.298 + temp = (U - lf_min) / val[q];
1.299 + }
1.300 + if (temp > 0.001) return 2;
1.301 + /* compute implied upper bound of x[q]; see (7), (8) */
1.302 + if (val[q] > 0.0)
1.303 + { if (U == +DBL_MAX || lf_min == -DBL_MAX)
1.304 + temp = +DBL_MAX;
1.305 + else
1.306 + temp = (U - lf_min) / val[q];
1.307 + }
1.308 + else
1.309 + { if (L == -DBL_MAX || lf_max == +DBL_MAX)
1.310 + temp = +DBL_MAX;
1.311 + else
1.312 + temp = (L - lf_max) / val[q];
1.313 + }
1.314 + if (temp < 0.999) return 1;
1.315 + /* there is no logical relation between x[p] and x[q] */
1.316 + return 0;
1.317 +}
1.318 +
1.319 +struct COG
1.320 +{ /* conflict graph; it represents logical relations between binary
1.321 + variables and has a vertex for each binary variable and its
1.322 + complement, and an edge between two vertices when at most one
1.323 + of the variables represented by the vertices can equal one in
1.324 + an optimal solution */
1.325 + int n;
1.326 + /* number of variables */
1.327 + int nb;
1.328 + /* number of binary variables represented in the graph (note that
1.329 + not all binary variables can be represented); vertices which
1.330 + correspond to binary variables have numbers 1, ..., nb while
1.331 + vertices which correspond to complements of binary variables
1.332 + have numbers nb+1, ..., nb+nb */
1.333 + int ne;
1.334 + /* number of edges in the graph */
1.335 + int *vert; /* int vert[1+n]; */
1.336 + /* if x[j] is a binary variable represented in the graph, vert[j]
1.337 + is the vertex number corresponding to x[j]; otherwise vert[j]
1.338 + is zero */
1.339 + int *orig; /* int list[1:nb]; */
1.340 + /* if vert[j] = k > 0, then orig[k] = j */
1.341 + unsigned char *a;
1.342 + /* adjacency matrix of the graph having 2*nb rows and columns;
1.343 + only strict lower triangle is stored in dense packed form */
1.344 +};
1.345 +
1.346 +/*----------------------------------------------------------------------
1.347 +-- lpx_create_cog - create the conflict graph.
1.348 +--
1.349 +-- SYNOPSIS
1.350 +--
1.351 +-- #include "glplpx.h"
1.352 +-- void *lpx_create_cog(LPX *lp);
1.353 +--
1.354 +-- DESCRIPTION
1.355 +--
1.356 +-- The routine lpx_create_cog creates the conflict graph for a given
1.357 +-- problem instance.
1.358 +--
1.359 +-- RETURNS
1.360 +--
1.361 +-- If the graph has been created, the routine returns a pointer to it.
1.362 +-- Otherwise the routine returns NULL. */
1.363 +
1.364 +#define MAX_NB 4000
1.365 +#define MAX_ROW_LEN 500
1.366 +
1.367 +static void lpx_add_cog_edge(void *_cog, int i, int j);
1.368 +
1.369 +static void *lpx_create_cog(LPX *lp)
1.370 +{ struct COG *cog = NULL;
1.371 + int m, n, nb, i, j, p, q, len, *ind, *vert, *orig;
1.372 + double L, U, lf_min, lf_max, *val;
1.373 + xprintf("Creating the conflict graph...\n");
1.374 + m = lpx_get_num_rows(lp);
1.375 + n = lpx_get_num_cols(lp);
1.376 + /* determine which binary variables should be included in the
1.377 + conflict graph */
1.378 + nb = 0;
1.379 + vert = xcalloc(1+n, sizeof(int));
1.380 + for (j = 1; j <= n; j++) vert[j] = 0;
1.381 + orig = xcalloc(1+n, sizeof(int));
1.382 + ind = xcalloc(1+n, sizeof(int));
1.383 + val = xcalloc(1+n, sizeof(double));
1.384 + for (i = 1; i <= m; i++)
1.385 + { L = get_row_lb(lp, i);
1.386 + U = get_row_ub(lp, i);
1.387 + if (L == -DBL_MAX && U == +DBL_MAX) continue;
1.388 + len = lpx_get_mat_row(lp, i, ind, val);
1.389 + if (len > MAX_ROW_LEN) continue;
1.390 + lf_min = eval_lf_min(lp, len, ind, val);
1.391 + lf_max = eval_lf_max(lp, len, ind, val);
1.392 + for (p = 1; p <= len; p++)
1.393 + { if (!is_binary(lp, ind[p])) continue;
1.394 + for (q = p+1; q <= len; q++)
1.395 + { if (!is_binary(lp, ind[q])) continue;
1.396 + if (probing(len, val, L, U, lf_min, lf_max, p, 0, q) ||
1.397 + probing(len, val, L, U, lf_min, lf_max, p, 1, q))
1.398 + { /* there is a logical relation */
1.399 + /* include the first variable in the graph */
1.400 + j = ind[p];
1.401 + if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j;
1.402 + /* incude the second variable in the graph */
1.403 + j = ind[q];
1.404 + if (vert[j] == 0) nb++, vert[j] = nb, orig[nb] = j;
1.405 + }
1.406 + }
1.407 + }
1.408 + }
1.409 + /* if the graph is either empty or has too many vertices, do not
1.410 + create it */
1.411 + if (nb == 0 || nb > MAX_NB)
1.412 + { xprintf("The conflict graph is either empty or too big\n");
1.413 + xfree(vert);
1.414 + xfree(orig);
1.415 + goto done;
1.416 + }
1.417 + /* create the conflict graph */
1.418 + cog = xmalloc(sizeof(struct COG));
1.419 + cog->n = n;
1.420 + cog->nb = nb;
1.421 + cog->ne = 0;
1.422 + cog->vert = vert;
1.423 + cog->orig = orig;
1.424 + len = nb + nb; /* number of vertices */
1.425 + len = (len * (len - 1)) / 2; /* number of entries in triangle */
1.426 + len = (len + (CHAR_BIT - 1)) / CHAR_BIT; /* bytes needed */
1.427 + cog->a = xmalloc(len);
1.428 + memset(cog->a, 0, len);
1.429 + for (j = 1; j <= nb; j++)
1.430 + { /* add edge between variable and its complement */
1.431 + lpx_add_cog_edge(cog, +orig[j], -orig[j]);
1.432 + }
1.433 + for (i = 1; i <= m; i++)
1.434 + { L = get_row_lb(lp, i);
1.435 + U = get_row_ub(lp, i);
1.436 + if (L == -DBL_MAX && U == +DBL_MAX) continue;
1.437 + len = lpx_get_mat_row(lp, i, ind, val);
1.438 + if (len > MAX_ROW_LEN) continue;
1.439 + lf_min = eval_lf_min(lp, len, ind, val);
1.440 + lf_max = eval_lf_max(lp, len, ind, val);
1.441 + for (p = 1; p <= len; p++)
1.442 + { if (!is_binary(lp, ind[p])) continue;
1.443 + for (q = p+1; q <= len; q++)
1.444 + { if (!is_binary(lp, ind[q])) continue;
1.445 + /* set x[p] to 0 and examine x[q] */
1.446 + switch (probing(len, val, L, U, lf_min, lf_max, p, 0, q))
1.447 + { case 0:
1.448 + /* no logical relation */
1.449 + break;
1.450 + case 1:
1.451 + /* x[p] = 0 implies x[q] = 0 */
1.452 + lpx_add_cog_edge(cog, -ind[p], +ind[q]);
1.453 + break;
1.454 + case 2:
1.455 + /* x[p] = 0 implies x[q] = 1 */
1.456 + lpx_add_cog_edge(cog, -ind[p], -ind[q]);
1.457 + break;
1.458 + default:
1.459 + xassert(lp != lp);
1.460 + }
1.461 + /* set x[p] to 1 and examine x[q] */
1.462 + switch (probing(len, val, L, U, lf_min, lf_max, p, 1, q))
1.463 + { case 0:
1.464 + /* no logical relation */
1.465 + break;
1.466 + case 1:
1.467 + /* x[p] = 1 implies x[q] = 0 */
1.468 + lpx_add_cog_edge(cog, +ind[p], +ind[q]);
1.469 + break;
1.470 + case 2:
1.471 + /* x[p] = 1 implies x[q] = 1 */
1.472 + lpx_add_cog_edge(cog, +ind[p], -ind[q]);
1.473 + break;
1.474 + default:
1.475 + xassert(lp != lp);
1.476 + }
1.477 + }
1.478 + }
1.479 + }
1.480 + xprintf("The conflict graph has 2*%d vertices and %d edges\n",
1.481 + cog->nb, cog->ne);
1.482 +done: xfree(ind);
1.483 + xfree(val);
1.484 + return cog;
1.485 +}
1.486 +
1.487 +/*----------------------------------------------------------------------
1.488 +-- lpx_add_cog_edge - add edge to the conflict graph.
1.489 +--
1.490 +-- SYNOPSIS
1.491 +--
1.492 +-- #include "glplpx.h"
1.493 +-- void lpx_add_cog_edge(void *cog, int i, int j);
1.494 +--
1.495 +-- DESCRIPTION
1.496 +--
1.497 +-- The routine lpx_add_cog_edge adds an edge to the conflict graph.
1.498 +-- The edge connects x[i] (if i > 0) or its complement (if i < 0) and
1.499 +-- x[j] (if j > 0) or its complement (if j < 0), where i and j are
1.500 +-- original ordinal numbers of corresponding variables. */
1.501 +
1.502 +static void lpx_add_cog_edge(void *_cog, int i, int j)
1.503 +{ struct COG *cog = _cog;
1.504 + int k;
1.505 + xassert(i != j);
1.506 + /* determine indices of corresponding vertices */
1.507 + if (i > 0)
1.508 + { xassert(1 <= i && i <= cog->n);
1.509 + i = cog->vert[i];
1.510 + xassert(i != 0);
1.511 + }
1.512 + else
1.513 + { i = -i;
1.514 + xassert(1 <= i && i <= cog->n);
1.515 + i = cog->vert[i];
1.516 + xassert(i != 0);
1.517 + i += cog->nb;
1.518 + }
1.519 + if (j > 0)
1.520 + { xassert(1 <= j && j <= cog->n);
1.521 + j = cog->vert[j];
1.522 + xassert(j != 0);
1.523 + }
1.524 + else
1.525 + { j = -j;
1.526 + xassert(1 <= j && j <= cog->n);
1.527 + j = cog->vert[j];
1.528 + xassert(j != 0);
1.529 + j += cog->nb;
1.530 + }
1.531 + /* only lower triangle is stored, so we need i > j */
1.532 + if (i < j) k = i, i = j, j = k;
1.533 + k = ((i - 1) * (i - 2)) / 2 + (j - 1);
1.534 + cog->a[k / CHAR_BIT] |=
1.535 + (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
1.536 + cog->ne++;
1.537 + return;
1.538 +}
1.539 +
1.540 +/*----------------------------------------------------------------------
1.541 +-- MAXIMUM WEIGHT CLIQUE
1.542 +--
1.543 +-- Two subroutines sub() and wclique() below are intended to find a
1.544 +-- maximum weight clique in a given undirected graph. These subroutines
1.545 +-- are slightly modified version of the program WCLIQUE developed by
1.546 +-- Patric Ostergard <http://www.tcs.hut.fi/~pat/wclique.html> and based
1.547 +-- on ideas from the article "P. R. J. Ostergard, A new algorithm for
1.548 +-- the maximum-weight clique problem, submitted for publication", which
1.549 +-- in turn is a generalization of the algorithm for unweighted graphs
1.550 +-- presented in "P. R. J. Ostergard, A fast algorithm for the maximum
1.551 +-- clique problem, submitted for publication".
1.552 +--
1.553 +-- USED WITH PERMISSION OF THE AUTHOR OF THE ORIGINAL CODE. */
1.554 +
1.555 +struct dsa
1.556 +{ /* dynamic storage area */
1.557 + int n;
1.558 + /* number of vertices */
1.559 + int *wt; /* int wt[0:n-1]; */
1.560 + /* weights */
1.561 + unsigned char *a;
1.562 + /* adjacency matrix (packed lower triangle without main diag.) */
1.563 + int record;
1.564 + /* weight of best clique */
1.565 + int rec_level;
1.566 + /* number of vertices in best clique */
1.567 + int *rec; /* int rec[0:n-1]; */
1.568 + /* best clique so far */
1.569 + int *clique; /* int clique[0:n-1]; */
1.570 + /* table for pruning */
1.571 + int *set; /* int set[0:n-1]; */
1.572 + /* current clique */
1.573 +};
1.574 +
1.575 +#define n (dsa->n)
1.576 +#define wt (dsa->wt)
1.577 +#define a (dsa->a)
1.578 +#define record (dsa->record)
1.579 +#define rec_level (dsa->rec_level)
1.580 +#define rec (dsa->rec)
1.581 +#define clique (dsa->clique)
1.582 +#define set (dsa->set)
1.583 +
1.584 +#if 0
1.585 +static int is_edge(struct dsa *dsa, int i, int j)
1.586 +{ /* if there is arc (i,j), the routine returns true; otherwise
1.587 + false; 0 <= i, j < n */
1.588 + int k;
1.589 + xassert(0 <= i && i < n);
1.590 + xassert(0 <= j && j < n);
1.591 + if (i == j) return 0;
1.592 + if (i < j) k = i, i = j, j = k;
1.593 + k = (i * (i - 1)) / 2 + j;
1.594 + return a[k / CHAR_BIT] &
1.595 + (unsigned char)(1 << ((CHAR_BIT - 1) - k % CHAR_BIT));
1.596 +}
1.597 +#else
1.598 +#define is_edge(dsa, i, j) ((i) == (j) ? 0 : \
1.599 + (i) > (j) ? is_edge1(i, j) : is_edge1(j, i))
1.600 +#define is_edge1(i, j) is_edge2(((i) * ((i) - 1)) / 2 + (j))
1.601 +#define is_edge2(k) (a[(k) / CHAR_BIT] & \
1.602 + (unsigned char)(1 << ((CHAR_BIT - 1) - (k) % CHAR_BIT)))
1.603 +#endif
1.604 +
1.605 +static void sub(struct dsa *dsa, int ct, int table[], int level,
1.606 + int weight, int l_weight)
1.607 +{ int i, j, k, curr_weight, left_weight, *p1, *p2, *newtable;
1.608 + newtable = xcalloc(n, sizeof(int));
1.609 + if (ct <= 0)
1.610 + { /* 0 or 1 elements left; include these */
1.611 + if (ct == 0)
1.612 + { set[level++] = table[0];
1.613 + weight += l_weight;
1.614 + }
1.615 + if (weight > record)
1.616 + { record = weight;
1.617 + rec_level = level;
1.618 + for (i = 0; i < level; i++) rec[i] = set[i];
1.619 + }
1.620 + goto done;
1.621 + }
1.622 + for (i = ct; i >= 0; i--)
1.623 + { if ((level == 0) && (i < ct)) goto done;
1.624 + k = table[i];
1.625 + if ((level > 0) && (clique[k] <= (record - weight)))
1.626 + goto done; /* prune */
1.627 + set[level] = k;
1.628 + curr_weight = weight + wt[k];
1.629 + l_weight -= wt[k];
1.630 + if (l_weight <= (record - curr_weight))
1.631 + goto done; /* prune */
1.632 + p1 = newtable;
1.633 + p2 = table;
1.634 + left_weight = 0;
1.635 + while (p2 < table + i)
1.636 + { j = *p2++;
1.637 + if (is_edge(dsa, j, k))
1.638 + { *p1++ = j;
1.639 + left_weight += wt[j];
1.640 + }
1.641 + }
1.642 + if (left_weight <= (record - curr_weight)) continue;
1.643 + sub(dsa, p1 - newtable - 1, newtable, level + 1, curr_weight,
1.644 + left_weight);
1.645 + }
1.646 +done: xfree(newtable);
1.647 + return;
1.648 +}
1.649 +
1.650 +static int wclique(int _n, int w[], unsigned char _a[], int sol[])
1.651 +{ struct dsa _dsa, *dsa = &_dsa;
1.652 + int i, j, p, max_wt, max_nwt, wth, *used, *nwt, *pos;
1.653 + glp_long timer;
1.654 + n = _n;
1.655 + wt = &w[1];
1.656 + a = _a;
1.657 + record = 0;
1.658 + rec_level = 0;
1.659 + rec = &sol[1];
1.660 + clique = xcalloc(n, sizeof(int));
1.661 + set = xcalloc(n, sizeof(int));
1.662 + used = xcalloc(n, sizeof(int));
1.663 + nwt = xcalloc(n, sizeof(int));
1.664 + pos = xcalloc(n, sizeof(int));
1.665 + /* start timer */
1.666 + timer = xtime();
1.667 + /* order vertices */
1.668 + for (i = 0; i < n; i++)
1.669 + { nwt[i] = 0;
1.670 + for (j = 0; j < n; j++)
1.671 + if (is_edge(dsa, i, j)) nwt[i] += wt[j];
1.672 + }
1.673 + for (i = 0; i < n; i++)
1.674 + used[i] = 0;
1.675 + for (i = n-1; i >= 0; i--)
1.676 + { max_wt = -1;
1.677 + max_nwt = -1;
1.678 + for (j = 0; j < n; j++)
1.679 + { if ((!used[j]) && ((wt[j] > max_wt) || (wt[j] == max_wt
1.680 + && nwt[j] > max_nwt)))
1.681 + { max_wt = wt[j];
1.682 + max_nwt = nwt[j];
1.683 + p = j;
1.684 + }
1.685 + }
1.686 + pos[i] = p;
1.687 + used[p] = 1;
1.688 + for (j = 0; j < n; j++)
1.689 + if ((!used[j]) && (j != p) && (is_edge(dsa, p, j)))
1.690 + nwt[j] -= wt[p];
1.691 + }
1.692 + /* main routine */
1.693 + wth = 0;
1.694 + for (i = 0; i < n; i++)
1.695 + { wth += wt[pos[i]];
1.696 + sub(dsa, i, pos, 0, 0, wth);
1.697 + clique[pos[i]] = record;
1.698 +#if 0
1.699 + if (utime() >= timer + 5.0)
1.700 +#else
1.701 + if (xdifftime(xtime(), timer) >= 5.0 - 0.001)
1.702 +#endif
1.703 + { /* print current record and reset timer */
1.704 + xprintf("level = %d (%d); best = %d\n", i+1, n, record);
1.705 +#if 0
1.706 + timer = utime();
1.707 +#else
1.708 + timer = xtime();
1.709 +#endif
1.710 + }
1.711 + }
1.712 + xfree(clique);
1.713 + xfree(set);
1.714 + xfree(used);
1.715 + xfree(nwt);
1.716 + xfree(pos);
1.717 + /* return the solution found */
1.718 + for (i = 1; i <= rec_level; i++) sol[i]++;
1.719 + return rec_level;
1.720 +}
1.721 +
1.722 +#undef n
1.723 +#undef wt
1.724 +#undef a
1.725 +#undef record
1.726 +#undef rec_level
1.727 +#undef rec
1.728 +#undef clique
1.729 +#undef set
1.730 +
1.731 +/*----------------------------------------------------------------------
1.732 +-- lpx_clique_cut - generate cluque cut.
1.733 +--
1.734 +-- SYNOPSIS
1.735 +--
1.736 +-- #include "glplpx.h"
1.737 +-- int lpx_clique_cut(LPX *lp, void *cog, int ind[], double val[]);
1.738 +--
1.739 +-- DESCRIPTION
1.740 +--
1.741 +-- The routine lpx_clique_cut generates a clique cut using the conflict
1.742 +-- graph specified by the parameter cog.
1.743 +--
1.744 +-- If a violated clique cut has been found, it has the following form:
1.745 +--
1.746 +-- sum{j in J} a[j]*x[j] <= b.
1.747 +--
1.748 +-- Variable indices j in J are stored in elements ind[1], ..., ind[len]
1.749 +-- while corresponding constraint coefficients are stored in elements
1.750 +-- val[1], ..., val[len], where len is returned on exit. The right-hand
1.751 +-- side b is stored in element val[0].
1.752 +--
1.753 +-- RETURNS
1.754 +--
1.755 +-- If the cutting plane has been successfully generated, the routine
1.756 +-- returns 1 <= len <= n, which is the number of non-zero coefficients
1.757 +-- in the inequality constraint. Otherwise, the routine returns zero. */
1.758 +
1.759 +static int lpx_clique_cut(LPX *lp, void *_cog, int ind[], double val[])
1.760 +{ struct COG *cog = _cog;
1.761 + int n = lpx_get_num_cols(lp);
1.762 + int j, t, v, card, temp, len = 0, *w, *sol;
1.763 + double x, sum, b, *vec;
1.764 + /* allocate working arrays */
1.765 + w = xcalloc(1 + 2 * cog->nb, sizeof(int));
1.766 + sol = xcalloc(1 + 2 * cog->nb, sizeof(int));
1.767 + vec = xcalloc(1+n, sizeof(double));
1.768 + /* assign weights to vertices of the conflict graph */
1.769 + for (t = 1; t <= cog->nb; t++)
1.770 + { j = cog->orig[t];
1.771 + x = lpx_get_col_prim(lp, j);
1.772 + temp = (int)(100.0 * x + 0.5);
1.773 + if (temp < 0) temp = 0;
1.774 + if (temp > 100) temp = 100;
1.775 + w[t] = temp;
1.776 + w[cog->nb + t] = 100 - temp;
1.777 + }
1.778 + /* find a clique of maximum weight */
1.779 + card = wclique(2 * cog->nb, w, cog->a, sol);
1.780 + /* compute the clique weight for unscaled values */
1.781 + sum = 0.0;
1.782 + for ( t = 1; t <= card; t++)
1.783 + { v = sol[t];
1.784 + xassert(1 <= v && v <= 2 * cog->nb);
1.785 + if (v <= cog->nb)
1.786 + { /* vertex v corresponds to binary variable x[j] */
1.787 + j = cog->orig[v];
1.788 + x = lpx_get_col_prim(lp, j);
1.789 + sum += x;
1.790 + }
1.791 + else
1.792 + { /* vertex v corresponds to the complement of x[j] */
1.793 + j = cog->orig[v - cog->nb];
1.794 + x = lpx_get_col_prim(lp, j);
1.795 + sum += 1.0 - x;
1.796 + }
1.797 + }
1.798 + /* if the sum of binary variables and their complements in the
1.799 + clique greater than 1, the clique cut is violated */
1.800 + if (sum >= 1.01)
1.801 + { /* construct the inquality */
1.802 + for (j = 1; j <= n; j++) vec[j] = 0;
1.803 + b = 1.0;
1.804 + for (t = 1; t <= card; t++)
1.805 + { v = sol[t];
1.806 + if (v <= cog->nb)
1.807 + { /* vertex v corresponds to binary variable x[j] */
1.808 + j = cog->orig[v];
1.809 + xassert(1 <= j && j <= n);
1.810 + vec[j] += 1.0;
1.811 + }
1.812 + else
1.813 + { /* vertex v corresponds to the complement of x[j] */
1.814 + j = cog->orig[v - cog->nb];
1.815 + xassert(1 <= j && j <= n);
1.816 + vec[j] -= 1.0;
1.817 + b -= 1.0;
1.818 + }
1.819 + }
1.820 + xassert(len == 0);
1.821 + for (j = 1; j <= n; j++)
1.822 + { if (vec[j] != 0.0)
1.823 + { len++;
1.824 + ind[len] = j, val[len] = vec[j];
1.825 + }
1.826 + }
1.827 + ind[0] = 0, val[0] = b;
1.828 + }
1.829 + /* free working arrays */
1.830 + xfree(w);
1.831 + xfree(sol);
1.832 + xfree(vec);
1.833 + /* return to the calling program */
1.834 + return len;
1.835 +}
1.836 +
1.837 +/*----------------------------------------------------------------------
1.838 +-- lpx_delete_cog - delete the conflict graph.
1.839 +--
1.840 +-- SYNOPSIS
1.841 +--
1.842 +-- #include "glplpx.h"
1.843 +-- void lpx_delete_cog(void *cog);
1.844 +--
1.845 +-- DESCRIPTION
1.846 +--
1.847 +-- The routine lpx_delete_cog deletes the conflict graph, which the
1.848 +-- parameter cog points to, freeing all the memory allocated to this
1.849 +-- object. */
1.850 +
1.851 +static void lpx_delete_cog(void *_cog)
1.852 +{ struct COG *cog = _cog;
1.853 + xfree(cog->vert);
1.854 + xfree(cog->orig);
1.855 + xfree(cog->a);
1.856 + xfree(cog);
1.857 +}
1.858 +
1.859 +/**********************************************************************/
1.860 +
1.861 +void *ios_clq_init(glp_tree *tree)
1.862 +{ /* initialize clique cut generator */
1.863 + glp_prob *mip = tree->mip;
1.864 + xassert(mip != NULL);
1.865 + return lpx_create_cog(mip);
1.866 +}
1.867 +
1.868 +/***********************************************************************
1.869 +* NAME
1.870 +*
1.871 +* ios_clq_gen - generate clique cuts
1.872 +*
1.873 +* SYNOPSIS
1.874 +*
1.875 +* #include "glpios.h"
1.876 +* void ios_clq_gen(glp_tree *tree, void *gen);
1.877 +*
1.878 +* DESCRIPTION
1.879 +*
1.880 +* The routine ios_clq_gen generates clique cuts for the current point
1.881 +* and adds them to the clique pool. */
1.882 +
1.883 +void ios_clq_gen(glp_tree *tree, void *gen)
1.884 +{ int n = lpx_get_num_cols(tree->mip);
1.885 + int len, *ind;
1.886 + double *val;
1.887 + xassert(gen != NULL);
1.888 + ind = xcalloc(1+n, sizeof(int));
1.889 + val = xcalloc(1+n, sizeof(double));
1.890 + len = lpx_clique_cut(tree->mip, gen, ind, val);
1.891 + if (len > 0)
1.892 + { /* xprintf("len = %d\n", len); */
1.893 + glp_ios_add_row(tree, NULL, GLP_RF_CLQ, 0, len, ind, val,
1.894 + GLP_UP, val[0]);
1.895 + }
1.896 + xfree(ind);
1.897 + xfree(val);
1.898 + return;
1.899 +}
1.900 +
1.901 +/**********************************************************************/
1.902 +
1.903 +void ios_clq_term(void *gen)
1.904 +{ /* terminate clique cut generator */
1.905 + xassert(gen != NULL);
1.906 + lpx_delete_cog(gen);
1.907 + return;
1.908 +}
1.909 +
1.910 +/* eof */