[9] | 1 | /* glpnpp04.c */ |
---|
| 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, 2011 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 "glpnpp.h" |
---|
| 26 | |
---|
| 27 | /*********************************************************************** |
---|
| 28 | * NAME |
---|
| 29 | * |
---|
| 30 | * npp_binarize_prob - binarize MIP problem |
---|
| 31 | * |
---|
| 32 | * SYNOPSIS |
---|
| 33 | * |
---|
| 34 | * #include "glpnpp.h" |
---|
| 35 | * int npp_binarize_prob(NPP *npp); |
---|
| 36 | * |
---|
| 37 | * DESCRIPTION |
---|
| 38 | * |
---|
| 39 | * The routine npp_binarize_prob replaces in the original MIP problem |
---|
| 40 | * every integer variable: |
---|
| 41 | * |
---|
| 42 | * l[q] <= x[q] <= u[q], (1) |
---|
| 43 | * |
---|
| 44 | * where l[q] < u[q], by an equivalent sum of binary variables. |
---|
| 45 | * |
---|
| 46 | * RETURNS |
---|
| 47 | * |
---|
| 48 | * The routine returns the number of integer variables for which the |
---|
| 49 | * transformation failed, because u[q] - l[q] > d_max. |
---|
| 50 | * |
---|
| 51 | * PROBLEM TRANSFORMATION |
---|
| 52 | * |
---|
| 53 | * If variable x[q] has non-zero lower bound, it is first processed |
---|
| 54 | * with the routine npp_lbnd_col. Thus, we can assume that: |
---|
| 55 | * |
---|
| 56 | * 0 <= x[q] <= u[q]. (2) |
---|
| 57 | * |
---|
| 58 | * If u[q] = 1, variable x[q] is already binary, so further processing |
---|
| 59 | * is not needed. Let, therefore, that 2 <= u[q] <= d_max, and n be a |
---|
| 60 | * smallest integer such that u[q] <= 2^n - 1 (n >= 2, since u[q] >= 2). |
---|
| 61 | * Then variable x[q] can be replaced by the following sum: |
---|
| 62 | * |
---|
| 63 | * n-1 |
---|
| 64 | * x[q] = sum 2^k x[k], (3) |
---|
| 65 | * k=0 |
---|
| 66 | * |
---|
| 67 | * where x[k] are binary columns (variables). If u[q] < 2^n - 1, the |
---|
| 68 | * following additional inequality constraint must be also included in |
---|
| 69 | * the transformed problem: |
---|
| 70 | * |
---|
| 71 | * n-1 |
---|
| 72 | * sum 2^k x[k] <= u[q]. (4) |
---|
| 73 | * k=0 |
---|
| 74 | * |
---|
| 75 | * Note: Assuming that in the transformed problem x[q] becomes binary |
---|
| 76 | * variable x[0], this transformation causes new n-1 binary variables |
---|
| 77 | * to appear. |
---|
| 78 | * |
---|
| 79 | * Substituting x[q] from (3) to the objective row gives: |
---|
| 80 | * |
---|
| 81 | * z = sum c[j] x[j] + c[0] = |
---|
| 82 | * j |
---|
| 83 | * |
---|
| 84 | * = sum c[j] x[j] + c[q] x[q] + c[0] = |
---|
| 85 | * j!=q |
---|
| 86 | * n-1 |
---|
| 87 | * = sum c[j] x[j] + c[q] sum 2^k x[k] + c[0] = |
---|
| 88 | * j!=q k=0 |
---|
| 89 | * n-1 |
---|
| 90 | * = sum c[j] x[j] + sum c[k] x[k] + c[0], |
---|
| 91 | * j!=q k=0 |
---|
| 92 | * |
---|
| 93 | * where: |
---|
| 94 | * |
---|
| 95 | * c[k] = 2^k c[q], k = 0, ..., n-1. (5) |
---|
| 96 | * |
---|
| 97 | * And substituting x[q] from (3) to i-th constraint row i gives: |
---|
| 98 | * |
---|
| 99 | * L[i] <= sum a[i,j] x[j] <= U[i] ==> |
---|
| 100 | * j |
---|
| 101 | * |
---|
| 102 | * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==> |
---|
| 103 | * j!=q |
---|
| 104 | * n-1 |
---|
| 105 | * L[i] <= sum a[i,j] x[j] + a[i,q] sum 2^k x[k] <= U[i] ==> |
---|
| 106 | * j!=q k=0 |
---|
| 107 | * n-1 |
---|
| 108 | * L[i] <= sum a[i,j] x[j] + sum a[i,k] x[k] <= U[i], |
---|
| 109 | * j!=q k=0 |
---|
| 110 | * |
---|
| 111 | * where: |
---|
| 112 | * |
---|
| 113 | * a[i,k] = 2^k a[i,q], k = 0, ..., n-1. (6) |
---|
| 114 | * |
---|
| 115 | * RECOVERING SOLUTION |
---|
| 116 | * |
---|
| 117 | * Value of variable x[q] is computed with formula (3). */ |
---|
| 118 | |
---|
| 119 | struct binarize |
---|
| 120 | { int q; |
---|
| 121 | /* column reference number for x[q] = x[0] */ |
---|
| 122 | int j; |
---|
| 123 | /* column reference number for x[1]; x[2] has reference number |
---|
| 124 | j+1, x[3] - j+2, etc. */ |
---|
| 125 | int n; |
---|
| 126 | /* total number of binary variables, n >= 2 */ |
---|
| 127 | }; |
---|
| 128 | |
---|
| 129 | static int rcv_binarize_prob(NPP *npp, void *info); |
---|
| 130 | |
---|
| 131 | int npp_binarize_prob(NPP *npp) |
---|
| 132 | { /* binarize MIP problem */ |
---|
| 133 | struct binarize *info; |
---|
| 134 | NPPROW *row; |
---|
| 135 | NPPCOL *col, *bin; |
---|
| 136 | NPPAIJ *aij; |
---|
| 137 | int u, n, k, temp, nfails, nvars, nbins, nrows; |
---|
| 138 | /* new variables will be added to the end of the column list, so |
---|
| 139 | we go from the end to beginning of the column list */ |
---|
| 140 | nfails = nvars = nbins = nrows = 0; |
---|
| 141 | for (col = npp->c_tail; col != NULL; col = col->prev) |
---|
| 142 | { /* skip continuous variable */ |
---|
| 143 | if (!col->is_int) continue; |
---|
| 144 | /* skip fixed variable */ |
---|
| 145 | if (col->lb == col->ub) continue; |
---|
| 146 | /* skip binary variable */ |
---|
| 147 | if (col->lb == 0.0 && col->ub == 1.0) continue; |
---|
| 148 | /* check if the transformation is applicable */ |
---|
| 149 | if (col->lb < -1e6 || col->ub > +1e6 || |
---|
| 150 | col->ub - col->lb > 4095.0) |
---|
| 151 | { /* unfortunately, not */ |
---|
| 152 | nfails++; |
---|
| 153 | continue; |
---|
| 154 | } |
---|
| 155 | /* process integer non-binary variable x[q] */ |
---|
| 156 | nvars++; |
---|
| 157 | /* make x[q] non-negative, if its lower bound is non-zero */ |
---|
| 158 | if (col->lb != 0.0) |
---|
| 159 | npp_lbnd_col(npp, col); |
---|
| 160 | /* now 0 <= x[q] <= u[q] */ |
---|
| 161 | xassert(col->lb == 0.0); |
---|
| 162 | u = (int)col->ub; |
---|
| 163 | xassert(col->ub == (double)u); |
---|
| 164 | /* if x[q] is binary, further processing is not needed */ |
---|
| 165 | if (u == 1) continue; |
---|
| 166 | /* determine smallest n such that u <= 2^n - 1 (thus, n is the |
---|
| 167 | number of binary variables needed) */ |
---|
| 168 | n = 2, temp = 4; |
---|
| 169 | while (u >= temp) |
---|
| 170 | n++, temp += temp; |
---|
| 171 | nbins += n; |
---|
| 172 | /* create transformation stack entry */ |
---|
| 173 | info = npp_push_tse(npp, |
---|
| 174 | rcv_binarize_prob, sizeof(struct binarize)); |
---|
| 175 | info->q = col->j; |
---|
| 176 | info->j = 0; /* will be set below */ |
---|
| 177 | info->n = n; |
---|
| 178 | /* if u < 2^n - 1, we need one additional row for (4) */ |
---|
| 179 | if (u < temp - 1) |
---|
| 180 | { row = npp_add_row(npp), nrows++; |
---|
| 181 | row->lb = -DBL_MAX, row->ub = u; |
---|
| 182 | } |
---|
| 183 | else |
---|
| 184 | row = NULL; |
---|
| 185 | /* in the transformed problem variable x[q] becomes binary |
---|
| 186 | variable x[0], so its objective and constraint coefficients |
---|
| 187 | are not changed */ |
---|
| 188 | col->ub = 1.0; |
---|
| 189 | /* include x[0] into constraint (4) */ |
---|
| 190 | if (row != NULL) |
---|
| 191 | npp_add_aij(npp, row, col, 1.0); |
---|
| 192 | /* add other binary variables x[1], ..., x[n-1] */ |
---|
| 193 | for (k = 1, temp = 2; k < n; k++, temp += temp) |
---|
| 194 | { /* add new binary variable x[k] */ |
---|
| 195 | bin = npp_add_col(npp); |
---|
| 196 | bin->is_int = 1; |
---|
| 197 | bin->lb = 0.0, bin->ub = 1.0; |
---|
| 198 | bin->coef = (double)temp * col->coef; |
---|
| 199 | /* store column reference number for x[1] */ |
---|
| 200 | if (info->j == 0) |
---|
| 201 | info->j = bin->j; |
---|
| 202 | else |
---|
| 203 | xassert(info->j + (k-1) == bin->j); |
---|
| 204 | /* duplicate constraint coefficients for x[k]; this also |
---|
| 205 | automatically includes x[k] into constraint (4) */ |
---|
| 206 | for (aij = col->ptr; aij != NULL; aij = aij->c_next) |
---|
| 207 | npp_add_aij(npp, aij->row, bin, (double)temp * aij->val); |
---|
| 208 | } |
---|
| 209 | } |
---|
| 210 | if (nvars > 0) |
---|
| 211 | xprintf("%d integer variable(s) were replaced by %d binary one" |
---|
| 212 | "s\n", nvars, nbins); |
---|
| 213 | if (nrows > 0) |
---|
| 214 | xprintf("%d row(s) were added due to binarization\n", nrows); |
---|
| 215 | if (nfails > 0) |
---|
| 216 | xprintf("Binarization failed for %d integer variable(s)\n", |
---|
| 217 | nfails); |
---|
| 218 | return nfails; |
---|
| 219 | } |
---|
| 220 | |
---|
| 221 | static int rcv_binarize_prob(NPP *npp, void *_info) |
---|
| 222 | { /* recovery binarized variable */ |
---|
| 223 | struct binarize *info = _info; |
---|
| 224 | int k, temp; |
---|
| 225 | double sum; |
---|
| 226 | /* compute value of x[q]; see formula (3) */ |
---|
| 227 | sum = npp->c_value[info->q]; |
---|
| 228 | for (k = 1, temp = 2; k < info->n; k++, temp += temp) |
---|
| 229 | sum += (double)temp * npp->c_value[info->j + (k-1)]; |
---|
| 230 | npp->c_value[info->q] = sum; |
---|
| 231 | return 0; |
---|
| 232 | } |
---|
| 233 | |
---|
| 234 | /**********************************************************************/ |
---|
| 235 | |
---|
| 236 | struct elem |
---|
| 237 | { /* linear form element a[j] x[j] */ |
---|
| 238 | double aj; |
---|
| 239 | /* non-zero coefficient value */ |
---|
| 240 | NPPCOL *xj; |
---|
| 241 | /* pointer to variable (column) */ |
---|
| 242 | struct elem *next; |
---|
| 243 | /* pointer to another term */ |
---|
| 244 | }; |
---|
| 245 | |
---|
| 246 | static struct elem *copy_form(NPP *npp, NPPROW *row, double s) |
---|
| 247 | { /* copy linear form */ |
---|
| 248 | NPPAIJ *aij; |
---|
| 249 | struct elem *ptr, *e; |
---|
| 250 | ptr = NULL; |
---|
| 251 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
| 252 | { e = dmp_get_atom(npp->pool, sizeof(struct elem)); |
---|
| 253 | e->aj = s * aij->val; |
---|
| 254 | e->xj = aij->col; |
---|
| 255 | e->next = ptr; |
---|
| 256 | ptr = e; |
---|
| 257 | } |
---|
| 258 | return ptr; |
---|
| 259 | } |
---|
| 260 | |
---|
| 261 | static void drop_form(NPP *npp, struct elem *ptr) |
---|
| 262 | { /* drop linear form */ |
---|
| 263 | struct elem *e; |
---|
| 264 | while (ptr != NULL) |
---|
| 265 | { e = ptr; |
---|
| 266 | ptr = e->next; |
---|
| 267 | dmp_free_atom(npp->pool, e, sizeof(struct elem)); |
---|
| 268 | } |
---|
| 269 | return; |
---|
| 270 | } |
---|
| 271 | |
---|
| 272 | /*********************************************************************** |
---|
| 273 | * NAME |
---|
| 274 | * |
---|
| 275 | * npp_is_packing - test if constraint is packing inequality |
---|
| 276 | * |
---|
| 277 | * SYNOPSIS |
---|
| 278 | * |
---|
| 279 | * #include "glpnpp.h" |
---|
| 280 | * int npp_is_packing(NPP *npp, NPPROW *row); |
---|
| 281 | * |
---|
| 282 | * RETURNS |
---|
| 283 | * |
---|
| 284 | * If the specified row (constraint) is packing inequality (see below), |
---|
| 285 | * the routine npp_is_packing returns non-zero. Otherwise, it returns |
---|
| 286 | * zero. |
---|
| 287 | * |
---|
| 288 | * PACKING INEQUALITIES |
---|
| 289 | * |
---|
| 290 | * In canonical format the packing inequality is the following: |
---|
| 291 | * |
---|
| 292 | * sum x[j] <= 1, (1) |
---|
| 293 | * j in J |
---|
| 294 | * |
---|
| 295 | * where all variables x[j] are binary. This inequality expresses the |
---|
| 296 | * condition that in any integer feasible solution at most one variable |
---|
| 297 | * from set J can take non-zero (unity) value while other variables |
---|
| 298 | * must be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because |
---|
| 299 | * if J is empty or |J| = 1, the inequality (1) is redundant. |
---|
| 300 | * |
---|
| 301 | * In general case the packing inequality may include original variables |
---|
| 302 | * x[j] as well as their complements x~[j]: |
---|
| 303 | * |
---|
| 304 | * sum x[j] + sum x~[j] <= 1, (2) |
---|
| 305 | * j in Jp j in Jn |
---|
| 306 | * |
---|
| 307 | * where Jp and Jn are not intersected. Therefore, using substitution |
---|
| 308 | * x~[j] = 1 - x[j] gives the packing inequality in generalized format: |
---|
| 309 | * |
---|
| 310 | * sum x[j] - sum x[j] <= 1 - |Jn|. (3) |
---|
| 311 | * j in Jp j in Jn */ |
---|
| 312 | |
---|
| 313 | int npp_is_packing(NPP *npp, NPPROW *row) |
---|
| 314 | { /* test if constraint is packing inequality */ |
---|
| 315 | NPPCOL *col; |
---|
| 316 | NPPAIJ *aij; |
---|
| 317 | int b; |
---|
| 318 | xassert(npp == npp); |
---|
| 319 | if (!(row->lb == -DBL_MAX && row->ub != +DBL_MAX)) |
---|
| 320 | return 0; |
---|
| 321 | b = 1; |
---|
| 322 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
| 323 | { col = aij->col; |
---|
| 324 | if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0)) |
---|
| 325 | return 0; |
---|
| 326 | if (aij->val == +1.0) |
---|
| 327 | ; |
---|
| 328 | else if (aij->val == -1.0) |
---|
| 329 | b--; |
---|
| 330 | else |
---|
| 331 | return 0; |
---|
| 332 | } |
---|
| 333 | if (row->ub != (double)b) return 0; |
---|
| 334 | return 1; |
---|
| 335 | } |
---|
| 336 | |
---|
| 337 | /*********************************************************************** |
---|
| 338 | * NAME |
---|
| 339 | * |
---|
| 340 | * npp_hidden_packing - identify hidden packing inequality |
---|
| 341 | * |
---|
| 342 | * SYNOPSIS |
---|
| 343 | * |
---|
| 344 | * #include "glpnpp.h" |
---|
| 345 | * int npp_hidden_packing(NPP *npp, NPPROW *row); |
---|
| 346 | * |
---|
| 347 | * DESCRIPTION |
---|
| 348 | * |
---|
| 349 | * The routine npp_hidden_packing processes specified inequality |
---|
| 350 | * constraint, which includes only binary variables, and the number of |
---|
| 351 | * the variables is not less than two. If the original inequality is |
---|
| 352 | * equivalent to a packing inequality, the routine replaces it by this |
---|
| 353 | * equivalent inequality. If the original constraint is double-sided |
---|
| 354 | * inequality, it is replaced by a pair of single-sided inequalities, |
---|
| 355 | * if necessary. |
---|
| 356 | * |
---|
| 357 | * RETURNS |
---|
| 358 | * |
---|
| 359 | * If the original inequality constraint was replaced by equivalent |
---|
| 360 | * packing inequality, the routine npp_hidden_packing returns non-zero. |
---|
| 361 | * Otherwise, it returns zero. |
---|
| 362 | * |
---|
| 363 | * PROBLEM TRANSFORMATION |
---|
| 364 | * |
---|
| 365 | * Consider an inequality constraint: |
---|
| 366 | * |
---|
| 367 | * sum a[j] x[j] <= b, (1) |
---|
| 368 | * j in J |
---|
| 369 | * |
---|
| 370 | * where all variables x[j] are binary, and |J| >= 2. (In case of '>=' |
---|
| 371 | * inequality it can be transformed to '<=' format by multiplying both |
---|
| 372 | * its sides by -1.) |
---|
| 373 | * |
---|
| 374 | * Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution |
---|
| 375 | * x[j] = 1 - x~[j] for all j in Jn, we have: |
---|
| 376 | * |
---|
| 377 | * sum a[j] x[j] <= b ==> |
---|
| 378 | * j in J |
---|
| 379 | * |
---|
| 380 | * sum a[j] x[j] + sum a[j] x[j] <= b ==> |
---|
| 381 | * j in Jp j in Jn |
---|
| 382 | * |
---|
| 383 | * sum a[j] x[j] + sum a[j] (1 - x~[j]) <= b ==> |
---|
| 384 | * j in Jp j in Jn |
---|
| 385 | * |
---|
| 386 | * sum a[j] x[j] - sum a[j] x~[j] <= b - sum a[j]. |
---|
| 387 | * j in Jp j in Jn j in Jn |
---|
| 388 | * |
---|
| 389 | * Thus, meaning the transformation above, we can assume that in |
---|
| 390 | * inequality (1) all coefficients a[j] are positive. Moreover, we can |
---|
| 391 | * assume that a[j] <= b. In fact, let a[j] > b; then the following |
---|
| 392 | * three cases are possible: |
---|
| 393 | * |
---|
| 394 | * 1) b < 0. In this case inequality (1) is infeasible, so the problem |
---|
| 395 | * has no feasible solution (see the routine npp_analyze_row); |
---|
| 396 | * |
---|
| 397 | * 2) b = 0. In this case inequality (1) is a forcing inequality on its |
---|
| 398 | * upper bound (see the routine npp_forcing row), from which it |
---|
| 399 | * follows that all variables x[j] should be fixed at zero; |
---|
| 400 | * |
---|
| 401 | * 3) b > 0. In this case inequality (1) defines an implied zero upper |
---|
| 402 | * bound for variable x[j] (see the routine npp_implied_bounds), from |
---|
| 403 | * which it follows that x[j] should be fixed at zero. |
---|
| 404 | * |
---|
| 405 | * It is assumed that all three cases listed above have been recognized |
---|
| 406 | * by the routine npp_process_prob, which performs basic MIP processing |
---|
| 407 | * prior to a call the routine npp_hidden_packing. So, if one of these |
---|
| 408 | * cases occurs, we should just skip processing such constraint. |
---|
| 409 | * |
---|
| 410 | * Thus, let 0 < a[j] <= b. Then it is obvious that constraint (1) is |
---|
| 411 | * equivalent to packing inquality only if: |
---|
| 412 | * |
---|
| 413 | * a[j] + a[k] > b + eps (2) |
---|
| 414 | * |
---|
| 415 | * for all j, k in J, j != k, where eps is an absolute tolerance for |
---|
| 416 | * row (linear form) value. Checking the condition (2) for all j and k, |
---|
| 417 | * j != k, requires time O(|J|^2). However, this time can be reduced to |
---|
| 418 | * O(|J|), if use minimal a[j] and a[k], in which case it is sufficient |
---|
| 419 | * to check the condition (2) only once. |
---|
| 420 | * |
---|
| 421 | * Once the original inequality (1) is replaced by equivalent packing |
---|
| 422 | * inequality, we need to perform back substitution x~[j] = 1 - x[j] for |
---|
| 423 | * all j in Jn (see above). |
---|
| 424 | * |
---|
| 425 | * RECOVERING SOLUTION |
---|
| 426 | * |
---|
| 427 | * None needed. */ |
---|
| 428 | |
---|
| 429 | static int hidden_packing(NPP *npp, struct elem *ptr, double *_b) |
---|
| 430 | { /* process inequality constraint: sum a[j] x[j] <= b; |
---|
| 431 | 0 - specified row is NOT hidden packing inequality; |
---|
| 432 | 1 - specified row is packing inequality; |
---|
| 433 | 2 - specified row is hidden packing inequality. */ |
---|
| 434 | struct elem *e, *ej, *ek; |
---|
| 435 | int neg; |
---|
| 436 | double b = *_b, eps; |
---|
| 437 | xassert(npp == npp); |
---|
| 438 | /* a[j] must be non-zero, x[j] must be binary, for all j in J */ |
---|
| 439 | for (e = ptr; e != NULL; e = e->next) |
---|
| 440 | { xassert(e->aj != 0.0); |
---|
| 441 | xassert(e->xj->is_int); |
---|
| 442 | xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0); |
---|
| 443 | } |
---|
| 444 | /* check if the specified inequality constraint already has the |
---|
| 445 | form of packing inequality */ |
---|
| 446 | neg = 0; /* neg is |Jn| */ |
---|
| 447 | for (e = ptr; e != NULL; e = e->next) |
---|
| 448 | { if (e->aj == +1.0) |
---|
| 449 | ; |
---|
| 450 | else if (e->aj == -1.0) |
---|
| 451 | neg++; |
---|
| 452 | else |
---|
| 453 | break; |
---|
| 454 | } |
---|
| 455 | if (e == NULL) |
---|
| 456 | { /* all coefficients a[j] are +1 or -1; check rhs b */ |
---|
| 457 | if (b == (double)(1 - neg)) |
---|
| 458 | { /* it is packing inequality; no processing is needed */ |
---|
| 459 | return 1; |
---|
| 460 | } |
---|
| 461 | } |
---|
| 462 | /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j] |
---|
| 463 | positive; the result is a~[j] = |a[j]| and new rhs b */ |
---|
| 464 | for (e = ptr; e != NULL; e = e->next) |
---|
| 465 | if (e->aj < 0) b -= e->aj; |
---|
| 466 | /* now a[j] > 0 for all j in J (actually |a[j]| are used) */ |
---|
| 467 | /* if a[j] > b, skip processing--this case must not appear */ |
---|
| 468 | for (e = ptr; e != NULL; e = e->next) |
---|
| 469 | if (fabs(e->aj) > b) return 0; |
---|
| 470 | /* now 0 < a[j] <= b for all j in J */ |
---|
| 471 | /* find two minimal coefficients a[j] and a[k], j != k */ |
---|
| 472 | ej = NULL; |
---|
| 473 | for (e = ptr; e != NULL; e = e->next) |
---|
| 474 | if (ej == NULL || fabs(ej->aj) > fabs(e->aj)) ej = e; |
---|
| 475 | xassert(ej != NULL); |
---|
| 476 | ek = NULL; |
---|
| 477 | for (e = ptr; e != NULL; e = e->next) |
---|
| 478 | if (e != ej) |
---|
| 479 | if (ek == NULL || fabs(ek->aj) > fabs(e->aj)) ek = e; |
---|
| 480 | xassert(ek != NULL); |
---|
| 481 | /* the specified constraint is equivalent to packing inequality |
---|
| 482 | iff a[j] + a[k] > b + eps */ |
---|
| 483 | eps = 1e-3 + 1e-6 * fabs(b); |
---|
| 484 | if (fabs(ej->aj) + fabs(ek->aj) <= b + eps) return 0; |
---|
| 485 | /* perform back substitution x~[j] = 1 - x[j] and construct the |
---|
| 486 | final equivalent packing inequality in generalized format */ |
---|
| 487 | b = 1.0; |
---|
| 488 | for (e = ptr; e != NULL; e = e->next) |
---|
| 489 | { if (e->aj > 0.0) |
---|
| 490 | e->aj = +1.0; |
---|
| 491 | else /* e->aj < 0.0 */ |
---|
| 492 | e->aj = -1.0, b -= 1.0; |
---|
| 493 | } |
---|
| 494 | *_b = b; |
---|
| 495 | return 2; |
---|
| 496 | } |
---|
| 497 | |
---|
| 498 | int npp_hidden_packing(NPP *npp, NPPROW *row) |
---|
| 499 | { /* identify hidden packing inequality */ |
---|
| 500 | NPPROW *copy; |
---|
| 501 | NPPAIJ *aij; |
---|
| 502 | struct elem *ptr, *e; |
---|
| 503 | int kase, ret, count = 0; |
---|
| 504 | double b; |
---|
| 505 | /* the row must be inequality constraint */ |
---|
| 506 | xassert(row->lb < row->ub); |
---|
| 507 | for (kase = 0; kase <= 1; kase++) |
---|
| 508 | { if (kase == 0) |
---|
| 509 | { /* process row upper bound */ |
---|
| 510 | if (row->ub == +DBL_MAX) continue; |
---|
| 511 | ptr = copy_form(npp, row, +1.0); |
---|
| 512 | b = + row->ub; |
---|
| 513 | } |
---|
| 514 | else |
---|
| 515 | { /* process row lower bound */ |
---|
| 516 | if (row->lb == -DBL_MAX) continue; |
---|
| 517 | ptr = copy_form(npp, row, -1.0); |
---|
| 518 | b = - row->lb; |
---|
| 519 | } |
---|
| 520 | /* now the inequality has the form "sum a[j] x[j] <= b" */ |
---|
| 521 | ret = hidden_packing(npp, ptr, &b); |
---|
| 522 | xassert(0 <= ret && ret <= 2); |
---|
| 523 | if (kase == 1 && ret == 1 || ret == 2) |
---|
| 524 | { /* the original inequality has been identified as hidden |
---|
| 525 | packing inequality */ |
---|
| 526 | count++; |
---|
| 527 | #ifdef GLP_DEBUG |
---|
| 528 | xprintf("Original constraint:\n"); |
---|
| 529 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
| 530 | xprintf(" %+g x%d", aij->val, aij->col->j); |
---|
| 531 | if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb); |
---|
| 532 | if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub); |
---|
| 533 | xprintf("\n"); |
---|
| 534 | xprintf("Equivalent packing inequality:\n"); |
---|
| 535 | for (e = ptr; e != NULL; e = e->next) |
---|
| 536 | xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j); |
---|
| 537 | xprintf(", <= %g\n", b); |
---|
| 538 | #endif |
---|
| 539 | if (row->lb == -DBL_MAX || row->ub == +DBL_MAX) |
---|
| 540 | { /* the original row is single-sided inequality; no copy |
---|
| 541 | is needed */ |
---|
| 542 | copy = NULL; |
---|
| 543 | } |
---|
| 544 | else |
---|
| 545 | { /* the original row is double-sided inequality; we need |
---|
| 546 | to create its copy for other bound before replacing it |
---|
| 547 | with the equivalent inequality */ |
---|
| 548 | copy = npp_add_row(npp); |
---|
| 549 | if (kase == 0) |
---|
| 550 | { /* the copy is for lower bound */ |
---|
| 551 | copy->lb = row->lb, copy->ub = +DBL_MAX; |
---|
| 552 | } |
---|
| 553 | else |
---|
| 554 | { /* the copy is for upper bound */ |
---|
| 555 | copy->lb = -DBL_MAX, copy->ub = row->ub; |
---|
| 556 | } |
---|
| 557 | /* copy original row coefficients */ |
---|
| 558 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
| 559 | npp_add_aij(npp, copy, aij->col, aij->val); |
---|
| 560 | } |
---|
| 561 | /* replace the original inequality by equivalent one */ |
---|
| 562 | npp_erase_row(npp, row); |
---|
| 563 | row->lb = -DBL_MAX, row->ub = b; |
---|
| 564 | for (e = ptr; e != NULL; e = e->next) |
---|
| 565 | npp_add_aij(npp, row, e->xj, e->aj); |
---|
| 566 | /* continue processing lower bound for the copy */ |
---|
| 567 | if (copy != NULL) row = copy; |
---|
| 568 | } |
---|
| 569 | drop_form(npp, ptr); |
---|
| 570 | } |
---|
| 571 | return count; |
---|
| 572 | } |
---|
| 573 | |
---|
| 574 | /*********************************************************************** |
---|
| 575 | * NAME |
---|
| 576 | * |
---|
| 577 | * npp_implied_packing - identify implied packing inequality |
---|
| 578 | * |
---|
| 579 | * SYNOPSIS |
---|
| 580 | * |
---|
| 581 | * #include "glpnpp.h" |
---|
| 582 | * int npp_implied_packing(NPP *npp, NPPROW *row, int which, |
---|
| 583 | * NPPCOL *var[], char set[]); |
---|
| 584 | * |
---|
| 585 | * DESCRIPTION |
---|
| 586 | * |
---|
| 587 | * The routine npp_implied_packing processes specified row (constraint) |
---|
| 588 | * of general format: |
---|
| 589 | * |
---|
| 590 | * L <= sum a[j] x[j] <= U. (1) |
---|
| 591 | * j |
---|
| 592 | * |
---|
| 593 | * If which = 0, only lower bound L, which must exist, is considered, |
---|
| 594 | * while upper bound U is ignored. Similarly, if which = 1, only upper |
---|
| 595 | * bound U, which must exist, is considered, while lower bound L is |
---|
| 596 | * ignored. Thus, if the specified row is a double-sided inequality or |
---|
| 597 | * equality constraint, this routine should be called twice for both |
---|
| 598 | * lower and upper bounds. |
---|
| 599 | * |
---|
| 600 | * The routine npp_implied_packing attempts to find a non-trivial (i.e. |
---|
| 601 | * having not less than two binary variables) packing inequality: |
---|
| 602 | * |
---|
| 603 | * sum x[j] - sum x[j] <= 1 - |Jn|, (2) |
---|
| 604 | * j in Jp j in Jn |
---|
| 605 | * |
---|
| 606 | * which is relaxation of the constraint (1) in the sense that any |
---|
| 607 | * solution satisfying to that constraint also satisfies to the packing |
---|
| 608 | * inequality (2). If such relaxation exists, the routine stores |
---|
| 609 | * pointers to descriptors of corresponding binary variables and their |
---|
| 610 | * flags, resp., to locations var[1], var[2], ..., var[len] and set[1], |
---|
| 611 | * set[2], ..., set[len], where set[j] = 0 means that j in Jp and |
---|
| 612 | * set[j] = 1 means that j in Jn. |
---|
| 613 | * |
---|
| 614 | * RETURNS |
---|
| 615 | * |
---|
| 616 | * The routine npp_implied_packing returns len, which is the total |
---|
| 617 | * number of binary variables in the packing inequality found, len >= 2. |
---|
| 618 | * However, if the relaxation does not exist, the routine returns zero. |
---|
| 619 | * |
---|
| 620 | * ALGORITHM |
---|
| 621 | * |
---|
| 622 | * If which = 0, the constraint coefficients (1) are multiplied by -1 |
---|
| 623 | * and b is assigned -L; if which = 1, the constraint coefficients (1) |
---|
| 624 | * are not changed and b is assigned +U. In both cases the specified |
---|
| 625 | * constraint gets the following format: |
---|
| 626 | * |
---|
| 627 | * sum a[j] x[j] <= b. (3) |
---|
| 628 | * j |
---|
| 629 | * |
---|
| 630 | * (Note that (3) is a relaxation of (1), because one of bounds L or U |
---|
| 631 | * is ignored.) |
---|
| 632 | * |
---|
| 633 | * Let J be set of binary variables, Kp be set of non-binary (integer |
---|
| 634 | * or continuous) variables with a[j] > 0, and Kn be set of non-binary |
---|
| 635 | * variables with a[j] < 0. Then the inequality (3) can be written as |
---|
| 636 | * follows: |
---|
| 637 | * |
---|
| 638 | * sum a[j] x[j] <= b - sum a[j] x[j] - sum a[j] x[j]. (4) |
---|
| 639 | * j in J j in Kp j in Kn |
---|
| 640 | * |
---|
| 641 | * To get rid of non-binary variables we can replace the inequality (4) |
---|
| 642 | * by the following relaxed inequality: |
---|
| 643 | * |
---|
| 644 | * sum a[j] x[j] <= b~, (5) |
---|
| 645 | * j in J |
---|
| 646 | * |
---|
| 647 | * where: |
---|
| 648 | * |
---|
| 649 | * b~ = sup(b - sum a[j] x[j] - sum a[j] x[j]) = |
---|
| 650 | * j in Kp j in Kn |
---|
| 651 | * |
---|
| 652 | * = b - inf sum a[j] x[j] - inf sum a[j] x[j] = (6) |
---|
| 653 | * j in Kp j in Kn |
---|
| 654 | * |
---|
| 655 | * = b - sum a[j] l[j] - sum a[j] u[j]. |
---|
| 656 | * j in Kp j in Kn |
---|
| 657 | * |
---|
| 658 | * Note that if lower bound l[j] (if j in Kp) or upper bound u[j] |
---|
| 659 | * (if j in Kn) of some non-binary variable x[j] does not exist, then |
---|
| 660 | * formally b = +oo, in which case further analysis is not performed. |
---|
| 661 | * |
---|
| 662 | * Let Bp = {j in J: a[j] > 0}, Bn = {j in J: a[j] < 0}. To make all |
---|
| 663 | * the inequality coefficients in (5) positive, we replace all x[j] in |
---|
| 664 | * Bn by their complementaries, substituting x[j] = 1 - x~[j] for all |
---|
| 665 | * j in Bn, that gives: |
---|
| 666 | * |
---|
| 667 | * sum a[j] x[j] - sum a[j] x~[j] <= b~ - sum a[j]. (7) |
---|
| 668 | * j in Bp j in Bn j in Bn |
---|
| 669 | * |
---|
| 670 | * This inequality is a relaxation of the original constraint (1), and |
---|
| 671 | * it is a binary knapsack inequality. Writing it in the standard format |
---|
| 672 | * we have: |
---|
| 673 | * |
---|
| 674 | * sum alfa[j] z[j] <= beta, (8) |
---|
| 675 | * j in J |
---|
| 676 | * |
---|
| 677 | * where: |
---|
| 678 | * ( + a[j], if j in Bp, |
---|
| 679 | * alfa[j] = < (9) |
---|
| 680 | * ( - a[j], if j in Bn, |
---|
| 681 | * |
---|
| 682 | * ( x[j], if j in Bp, |
---|
| 683 | * z[j] = < (10) |
---|
| 684 | * ( 1 - x[j], if j in Bn, |
---|
| 685 | * |
---|
| 686 | * beta = b~ - sum a[j]. (11) |
---|
| 687 | * j in Bn |
---|
| 688 | * |
---|
| 689 | * In the inequality (8) all coefficients are positive, therefore, the |
---|
| 690 | * packing relaxation to be found for this inequality is the following: |
---|
| 691 | * |
---|
| 692 | * sum z[j] <= 1. (12) |
---|
| 693 | * j in P |
---|
| 694 | * |
---|
| 695 | * It is obvious that set P within J, which we would like to find, must |
---|
| 696 | * satisfy to the following condition: |
---|
| 697 | * |
---|
| 698 | * alfa[j] + alfa[k] > beta + eps for all j, k in P, j != k, (13) |
---|
| 699 | * |
---|
| 700 | * where eps is an absolute tolerance for value of the linear form. |
---|
| 701 | * Thus, it is natural to take P = {j: alpha[j] > (beta + eps) / 2}. |
---|
| 702 | * Moreover, if in the equality (8) there exist coefficients alfa[k], |
---|
| 703 | * for which alfa[k] <= (beta + eps) / 2, but which, nevertheless, |
---|
| 704 | * satisfies to the condition (13) for all j in P, *one* corresponding |
---|
| 705 | * variable z[k] (having, for example, maximal coefficient alfa[k]) can |
---|
| 706 | * be included in set P, that allows increasing the number of binary |
---|
| 707 | * variables in (12) by one. |
---|
| 708 | * |
---|
| 709 | * Once the set P has been built, for the inequality (12) we need to |
---|
| 710 | * perform back substitution according to (10) in order to express it |
---|
| 711 | * through the original binary variables. As the result of such back |
---|
| 712 | * substitution the relaxed packing inequality get its final format (2), |
---|
| 713 | * where Jp = J intersect Bp, and Jn = J intersect Bn. */ |
---|
| 714 | |
---|
| 715 | int npp_implied_packing(NPP *npp, NPPROW *row, int which, |
---|
| 716 | NPPCOL *var[], char set[]) |
---|
| 717 | { struct elem *ptr, *e, *i, *k; |
---|
| 718 | int len = 0; |
---|
| 719 | double b, eps; |
---|
| 720 | /* build inequality (3) */ |
---|
| 721 | if (which == 0) |
---|
| 722 | { ptr = copy_form(npp, row, -1.0); |
---|
| 723 | xassert(row->lb != -DBL_MAX); |
---|
| 724 | b = - row->lb; |
---|
| 725 | } |
---|
| 726 | else if (which == 1) |
---|
| 727 | { ptr = copy_form(npp, row, +1.0); |
---|
| 728 | xassert(row->ub != +DBL_MAX); |
---|
| 729 | b = + row->ub; |
---|
| 730 | } |
---|
| 731 | /* remove non-binary variables to build relaxed inequality (5); |
---|
| 732 | compute its right-hand side b~ with formula (6) */ |
---|
| 733 | for (e = ptr; e != NULL; e = e->next) |
---|
| 734 | { if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0)) |
---|
| 735 | { /* x[j] is non-binary variable */ |
---|
| 736 | if (e->aj > 0.0) |
---|
| 737 | { if (e->xj->lb == -DBL_MAX) goto done; |
---|
| 738 | b -= e->aj * e->xj->lb; |
---|
| 739 | } |
---|
| 740 | else /* e->aj < 0.0 */ |
---|
| 741 | { if (e->xj->ub == +DBL_MAX) goto done; |
---|
| 742 | b -= e->aj * e->xj->ub; |
---|
| 743 | } |
---|
| 744 | /* a[j] = 0 means that variable x[j] is removed */ |
---|
| 745 | e->aj = 0.0; |
---|
| 746 | } |
---|
| 747 | } |
---|
| 748 | /* substitute x[j] = 1 - x~[j] to build knapsack inequality (8); |
---|
| 749 | compute its right-hand side beta with formula (11) */ |
---|
| 750 | for (e = ptr; e != NULL; e = e->next) |
---|
| 751 | if (e->aj < 0.0) b -= e->aj; |
---|
| 752 | /* if beta is close to zero, the knapsack inequality is either |
---|
| 753 | infeasible or forcing inequality; this must never happen, so |
---|
| 754 | we skip further analysis */ |
---|
| 755 | if (b < 1e-3) goto done; |
---|
| 756 | /* build set P as well as sets Jp and Jn, and determine x[k] as |
---|
| 757 | explained above in comments to the routine */ |
---|
| 758 | eps = 1e-3 + 1e-6 * b; |
---|
| 759 | i = k = NULL; |
---|
| 760 | for (e = ptr; e != NULL; e = e->next) |
---|
| 761 | { /* note that alfa[j] = |a[j]| */ |
---|
| 762 | if (fabs(e->aj) > 0.5 * (b + eps)) |
---|
| 763 | { /* alfa[j] > (b + eps) / 2; include x[j] in set P, i.e. in |
---|
| 764 | set Jp or Jn */ |
---|
| 765 | var[++len] = e->xj; |
---|
| 766 | set[len] = (char)(e->aj > 0.0 ? 0 : 1); |
---|
| 767 | /* alfa[i] = min alfa[j] over all j included in set P */ |
---|
| 768 | if (i == NULL || fabs(i->aj) > fabs(e->aj)) i = e; |
---|
| 769 | } |
---|
| 770 | else if (fabs(e->aj) >= 1e-3) |
---|
| 771 | { /* alfa[k] = max alfa[j] over all j not included in set P; |
---|
| 772 | we skip coefficient a[j] if it is close to zero to avoid |
---|
| 773 | numerically unreliable results */ |
---|
| 774 | if (k == NULL || fabs(k->aj) < fabs(e->aj)) k = e; |
---|
| 775 | } |
---|
| 776 | } |
---|
| 777 | /* if alfa[k] satisfies to condition (13) for all j in P, include |
---|
| 778 | x[k] in P */ |
---|
| 779 | if (i != NULL && k != NULL && fabs(i->aj) + fabs(k->aj) > b + eps) |
---|
| 780 | { var[++len] = k->xj; |
---|
| 781 | set[len] = (char)(k->aj > 0.0 ? 0 : 1); |
---|
| 782 | } |
---|
| 783 | /* trivial packing inequality being redundant must never appear, |
---|
| 784 | so we just ignore it */ |
---|
| 785 | if (len < 2) len = 0; |
---|
| 786 | done: drop_form(npp, ptr); |
---|
| 787 | return len; |
---|
| 788 | } |
---|
| 789 | |
---|
| 790 | /*********************************************************************** |
---|
| 791 | * NAME |
---|
| 792 | * |
---|
| 793 | * npp_is_covering - test if constraint is covering inequality |
---|
| 794 | * |
---|
| 795 | * SYNOPSIS |
---|
| 796 | * |
---|
| 797 | * #include "glpnpp.h" |
---|
| 798 | * int npp_is_covering(NPP *npp, NPPROW *row); |
---|
| 799 | * |
---|
| 800 | * RETURNS |
---|
| 801 | * |
---|
| 802 | * If the specified row (constraint) is covering inequality (see below), |
---|
| 803 | * the routine npp_is_covering returns non-zero. Otherwise, it returns |
---|
| 804 | * zero. |
---|
| 805 | * |
---|
| 806 | * COVERING INEQUALITIES |
---|
| 807 | * |
---|
| 808 | * In canonical format the covering inequality is the following: |
---|
| 809 | * |
---|
| 810 | * sum x[j] >= 1, (1) |
---|
| 811 | * j in J |
---|
| 812 | * |
---|
| 813 | * where all variables x[j] are binary. This inequality expresses the |
---|
| 814 | * condition that in any integer feasible solution variables in set J |
---|
| 815 | * cannot be all equal to zero at the same time, i.e. at least one |
---|
| 816 | * variable must take non-zero (unity) value. W.l.o.g. it is assumed |
---|
| 817 | * that |J| >= 2, because if J is empty, the inequality (1) is |
---|
| 818 | * infeasible, and if |J| = 1, the inequality (1) is a forcing row. |
---|
| 819 | * |
---|
| 820 | * In general case the covering inequality may include original |
---|
| 821 | * variables x[j] as well as their complements x~[j]: |
---|
| 822 | * |
---|
| 823 | * sum x[j] + sum x~[j] >= 1, (2) |
---|
| 824 | * j in Jp j in Jn |
---|
| 825 | * |
---|
| 826 | * where Jp and Jn are not intersected. Therefore, using substitution |
---|
| 827 | * x~[j] = 1 - x[j] gives the packing inequality in generalized format: |
---|
| 828 | * |
---|
| 829 | * sum x[j] - sum x[j] >= 1 - |Jn|. (3) |
---|
| 830 | * j in Jp j in Jn |
---|
| 831 | * |
---|
| 832 | * (May note that the inequality (3) cuts off infeasible solutions, |
---|
| 833 | * where x[j] = 0 for all j in Jp and x[j] = 1 for all j in Jn.) |
---|
| 834 | * |
---|
| 835 | * NOTE: If |J| = 2, the inequality (3) is equivalent to packing |
---|
| 836 | * inequality (see the routine npp_is_packing). */ |
---|
| 837 | |
---|
| 838 | int npp_is_covering(NPP *npp, NPPROW *row) |
---|
| 839 | { /* test if constraint is covering inequality */ |
---|
| 840 | NPPCOL *col; |
---|
| 841 | NPPAIJ *aij; |
---|
| 842 | int b; |
---|
| 843 | xassert(npp == npp); |
---|
| 844 | if (!(row->lb != -DBL_MAX && row->ub == +DBL_MAX)) |
---|
| 845 | return 0; |
---|
| 846 | b = 1; |
---|
| 847 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
| 848 | { col = aij->col; |
---|
| 849 | if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0)) |
---|
| 850 | return 0; |
---|
| 851 | if (aij->val == +1.0) |
---|
| 852 | ; |
---|
| 853 | else if (aij->val == -1.0) |
---|
| 854 | b--; |
---|
| 855 | else |
---|
| 856 | return 0; |
---|
| 857 | } |
---|
| 858 | if (row->lb != (double)b) return 0; |
---|
| 859 | return 1; |
---|
| 860 | } |
---|
| 861 | |
---|
| 862 | /*********************************************************************** |
---|
| 863 | * NAME |
---|
| 864 | * |
---|
| 865 | * npp_hidden_covering - identify hidden covering inequality |
---|
| 866 | * |
---|
| 867 | * SYNOPSIS |
---|
| 868 | * |
---|
| 869 | * #include "glpnpp.h" |
---|
| 870 | * int npp_hidden_covering(NPP *npp, NPPROW *row); |
---|
| 871 | * |
---|
| 872 | * DESCRIPTION |
---|
| 873 | * |
---|
| 874 | * The routine npp_hidden_covering processes specified inequality |
---|
| 875 | * constraint, which includes only binary variables, and the number of |
---|
| 876 | * the variables is not less than three. If the original inequality is |
---|
| 877 | * equivalent to a covering inequality (see below), the routine |
---|
| 878 | * replaces it by the equivalent inequality. If the original constraint |
---|
| 879 | * is double-sided inequality, it is replaced by a pair of single-sided |
---|
| 880 | * inequalities, if necessary. |
---|
| 881 | * |
---|
| 882 | * RETURNS |
---|
| 883 | * |
---|
| 884 | * If the original inequality constraint was replaced by equivalent |
---|
| 885 | * covering inequality, the routine npp_hidden_covering returns |
---|
| 886 | * non-zero. Otherwise, it returns zero. |
---|
| 887 | * |
---|
| 888 | * PROBLEM TRANSFORMATION |
---|
| 889 | * |
---|
| 890 | * Consider an inequality constraint: |
---|
| 891 | * |
---|
| 892 | * sum a[j] x[j] >= b, (1) |
---|
| 893 | * j in J |
---|
| 894 | * |
---|
| 895 | * where all variables x[j] are binary, and |J| >= 3. (In case of '<=' |
---|
| 896 | * inequality it can be transformed to '>=' format by multiplying both |
---|
| 897 | * its sides by -1.) |
---|
| 898 | * |
---|
| 899 | * Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution |
---|
| 900 | * x[j] = 1 - x~[j] for all j in Jn, we have: |
---|
| 901 | * |
---|
| 902 | * sum a[j] x[j] >= b ==> |
---|
| 903 | * j in J |
---|
| 904 | * |
---|
| 905 | * sum a[j] x[j] + sum a[j] x[j] >= b ==> |
---|
| 906 | * j in Jp j in Jn |
---|
| 907 | * |
---|
| 908 | * sum a[j] x[j] + sum a[j] (1 - x~[j]) >= b ==> |
---|
| 909 | * j in Jp j in Jn |
---|
| 910 | * |
---|
| 911 | * sum m a[j] x[j] - sum a[j] x~[j] >= b - sum a[j]. |
---|
| 912 | * j in Jp j in Jn j in Jn |
---|
| 913 | * |
---|
| 914 | * Thus, meaning the transformation above, we can assume that in |
---|
| 915 | * inequality (1) all coefficients a[j] are positive. Moreover, we can |
---|
| 916 | * assume that b > 0, because otherwise the inequality (1) would be |
---|
| 917 | * redundant (see the routine npp_analyze_row). It is then obvious that |
---|
| 918 | * constraint (1) is equivalent to covering inequality only if: |
---|
| 919 | * |
---|
| 920 | * a[j] >= b, (2) |
---|
| 921 | * |
---|
| 922 | * for all j in J. |
---|
| 923 | * |
---|
| 924 | * Once the original inequality (1) is replaced by equivalent covering |
---|
| 925 | * inequality, we need to perform back substitution x~[j] = 1 - x[j] for |
---|
| 926 | * all j in Jn (see above). |
---|
| 927 | * |
---|
| 928 | * RECOVERING SOLUTION |
---|
| 929 | * |
---|
| 930 | * None needed. */ |
---|
| 931 | |
---|
| 932 | static int hidden_covering(NPP *npp, struct elem *ptr, double *_b) |
---|
| 933 | { /* process inequality constraint: sum a[j] x[j] >= b; |
---|
| 934 | 0 - specified row is NOT hidden covering inequality; |
---|
| 935 | 1 - specified row is covering inequality; |
---|
| 936 | 2 - specified row is hidden covering inequality. */ |
---|
| 937 | struct elem *e; |
---|
| 938 | int neg; |
---|
| 939 | double b = *_b, eps; |
---|
| 940 | xassert(npp == npp); |
---|
| 941 | /* a[j] must be non-zero, x[j] must be binary, for all j in J */ |
---|
| 942 | for (e = ptr; e != NULL; e = e->next) |
---|
| 943 | { xassert(e->aj != 0.0); |
---|
| 944 | xassert(e->xj->is_int); |
---|
| 945 | xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0); |
---|
| 946 | } |
---|
| 947 | /* check if the specified inequality constraint already has the |
---|
| 948 | form of covering inequality */ |
---|
| 949 | neg = 0; /* neg is |Jn| */ |
---|
| 950 | for (e = ptr; e != NULL; e = e->next) |
---|
| 951 | { if (e->aj == +1.0) |
---|
| 952 | ; |
---|
| 953 | else if (e->aj == -1.0) |
---|
| 954 | neg++; |
---|
| 955 | else |
---|
| 956 | break; |
---|
| 957 | } |
---|
| 958 | if (e == NULL) |
---|
| 959 | { /* all coefficients a[j] are +1 or -1; check rhs b */ |
---|
| 960 | if (b == (double)(1 - neg)) |
---|
| 961 | { /* it is covering inequality; no processing is needed */ |
---|
| 962 | return 1; |
---|
| 963 | } |
---|
| 964 | } |
---|
| 965 | /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j] |
---|
| 966 | positive; the result is a~[j] = |a[j]| and new rhs b */ |
---|
| 967 | for (e = ptr; e != NULL; e = e->next) |
---|
| 968 | if (e->aj < 0) b -= e->aj; |
---|
| 969 | /* now a[j] > 0 for all j in J (actually |a[j]| are used) */ |
---|
| 970 | /* if b <= 0, skip processing--this case must not appear */ |
---|
| 971 | if (b < 1e-3) return 0; |
---|
| 972 | /* now a[j] > 0 for all j in J, and b > 0 */ |
---|
| 973 | /* the specified constraint is equivalent to covering inequality |
---|
| 974 | iff a[j] >= b for all j in J */ |
---|
| 975 | eps = 1e-9 + 1e-12 * fabs(b); |
---|
| 976 | for (e = ptr; e != NULL; e = e->next) |
---|
| 977 | if (fabs(e->aj) < b - eps) return 0; |
---|
| 978 | /* perform back substitution x~[j] = 1 - x[j] and construct the |
---|
| 979 | final equivalent covering inequality in generalized format */ |
---|
| 980 | b = 1.0; |
---|
| 981 | for (e = ptr; e != NULL; e = e->next) |
---|
| 982 | { if (e->aj > 0.0) |
---|
| 983 | e->aj = +1.0; |
---|
| 984 | else /* e->aj < 0.0 */ |
---|
| 985 | e->aj = -1.0, b -= 1.0; |
---|
| 986 | } |
---|
| 987 | *_b = b; |
---|
| 988 | return 2; |
---|
| 989 | } |
---|
| 990 | |
---|
| 991 | int npp_hidden_covering(NPP *npp, NPPROW *row) |
---|
| 992 | { /* identify hidden covering inequality */ |
---|
| 993 | NPPROW *copy; |
---|
| 994 | NPPAIJ *aij; |
---|
| 995 | struct elem *ptr, *e; |
---|
| 996 | int kase, ret, count = 0; |
---|
| 997 | double b; |
---|
| 998 | /* the row must be inequality constraint */ |
---|
| 999 | xassert(row->lb < row->ub); |
---|
| 1000 | for (kase = 0; kase <= 1; kase++) |
---|
| 1001 | { if (kase == 0) |
---|
| 1002 | { /* process row lower bound */ |
---|
| 1003 | if (row->lb == -DBL_MAX) continue; |
---|
| 1004 | ptr = copy_form(npp, row, +1.0); |
---|
| 1005 | b = + row->lb; |
---|
| 1006 | } |
---|
| 1007 | else |
---|
| 1008 | { /* process row upper bound */ |
---|
| 1009 | if (row->ub == +DBL_MAX) continue; |
---|
| 1010 | ptr = copy_form(npp, row, -1.0); |
---|
| 1011 | b = - row->ub; |
---|
| 1012 | } |
---|
| 1013 | /* now the inequality has the form "sum a[j] x[j] >= b" */ |
---|
| 1014 | ret = hidden_covering(npp, ptr, &b); |
---|
| 1015 | xassert(0 <= ret && ret <= 2); |
---|
| 1016 | if (kase == 1 && ret == 1 || ret == 2) |
---|
| 1017 | { /* the original inequality has been identified as hidden |
---|
| 1018 | covering inequality */ |
---|
| 1019 | count++; |
---|
| 1020 | #ifdef GLP_DEBUG |
---|
| 1021 | xprintf("Original constraint:\n"); |
---|
| 1022 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
| 1023 | xprintf(" %+g x%d", aij->val, aij->col->j); |
---|
| 1024 | if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb); |
---|
| 1025 | if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub); |
---|
| 1026 | xprintf("\n"); |
---|
| 1027 | xprintf("Equivalent covering inequality:\n"); |
---|
| 1028 | for (e = ptr; e != NULL; e = e->next) |
---|
| 1029 | xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j); |
---|
| 1030 | xprintf(", >= %g\n", b); |
---|
| 1031 | #endif |
---|
| 1032 | if (row->lb == -DBL_MAX || row->ub == +DBL_MAX) |
---|
| 1033 | { /* the original row is single-sided inequality; no copy |
---|
| 1034 | is needed */ |
---|
| 1035 | copy = NULL; |
---|
| 1036 | } |
---|
| 1037 | else |
---|
| 1038 | { /* the original row is double-sided inequality; we need |
---|
| 1039 | to create its copy for other bound before replacing it |
---|
| 1040 | with the equivalent inequality */ |
---|
| 1041 | copy = npp_add_row(npp); |
---|
| 1042 | if (kase == 0) |
---|
| 1043 | { /* the copy is for upper bound */ |
---|
| 1044 | copy->lb = -DBL_MAX, copy->ub = row->ub; |
---|
| 1045 | } |
---|
| 1046 | else |
---|
| 1047 | { /* the copy is for lower bound */ |
---|
| 1048 | copy->lb = row->lb, copy->ub = +DBL_MAX; |
---|
| 1049 | } |
---|
| 1050 | /* copy original row coefficients */ |
---|
| 1051 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
| 1052 | npp_add_aij(npp, copy, aij->col, aij->val); |
---|
| 1053 | } |
---|
| 1054 | /* replace the original inequality by equivalent one */ |
---|
| 1055 | npp_erase_row(npp, row); |
---|
| 1056 | row->lb = b, row->ub = +DBL_MAX; |
---|
| 1057 | for (e = ptr; e != NULL; e = e->next) |
---|
| 1058 | npp_add_aij(npp, row, e->xj, e->aj); |
---|
| 1059 | /* continue processing upper bound for the copy */ |
---|
| 1060 | if (copy != NULL) row = copy; |
---|
| 1061 | } |
---|
| 1062 | drop_form(npp, ptr); |
---|
| 1063 | } |
---|
| 1064 | return count; |
---|
| 1065 | } |
---|
| 1066 | |
---|
| 1067 | /*********************************************************************** |
---|
| 1068 | * NAME |
---|
| 1069 | * |
---|
| 1070 | * npp_is_partitioning - test if constraint is partitioning equality |
---|
| 1071 | * |
---|
| 1072 | * SYNOPSIS |
---|
| 1073 | * |
---|
| 1074 | * #include "glpnpp.h" |
---|
| 1075 | * int npp_is_partitioning(NPP *npp, NPPROW *row); |
---|
| 1076 | * |
---|
| 1077 | * RETURNS |
---|
| 1078 | * |
---|
| 1079 | * If the specified row (constraint) is partitioning equality (see |
---|
| 1080 | * below), the routine npp_is_partitioning returns non-zero. Otherwise, |
---|
| 1081 | * it returns zero. |
---|
| 1082 | * |
---|
| 1083 | * PARTITIONING EQUALITIES |
---|
| 1084 | * |
---|
| 1085 | * In canonical format the partitioning equality is the following: |
---|
| 1086 | * |
---|
| 1087 | * sum x[j] = 1, (1) |
---|
| 1088 | * j in J |
---|
| 1089 | * |
---|
| 1090 | * where all variables x[j] are binary. This equality expresses the |
---|
| 1091 | * condition that in any integer feasible solution exactly one variable |
---|
| 1092 | * in set J must take non-zero (unity) value while other variables must |
---|
| 1093 | * be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because if |
---|
| 1094 | * J is empty, the inequality (1) is infeasible, and if |J| = 1, the |
---|
| 1095 | * inequality (1) is a fixing row. |
---|
| 1096 | * |
---|
| 1097 | * In general case the partitioning equality may include original |
---|
| 1098 | * variables x[j] as well as their complements x~[j]: |
---|
| 1099 | * |
---|
| 1100 | * sum x[j] + sum x~[j] = 1, (2) |
---|
| 1101 | * j in Jp j in Jn |
---|
| 1102 | * |
---|
| 1103 | * where Jp and Jn are not intersected. Therefore, using substitution |
---|
| 1104 | * x~[j] = 1 - x[j] leads to the partitioning equality in generalized |
---|
| 1105 | * format: |
---|
| 1106 | * |
---|
| 1107 | * sum x[j] - sum x[j] = 1 - |Jn|. (3) |
---|
| 1108 | * j in Jp j in Jn */ |
---|
| 1109 | |
---|
| 1110 | int npp_is_partitioning(NPP *npp, NPPROW *row) |
---|
| 1111 | { /* test if constraint is partitioning equality */ |
---|
| 1112 | NPPCOL *col; |
---|
| 1113 | NPPAIJ *aij; |
---|
| 1114 | int b; |
---|
| 1115 | xassert(npp == npp); |
---|
| 1116 | if (row->lb != row->ub) return 0; |
---|
| 1117 | b = 1; |
---|
| 1118 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
| 1119 | { col = aij->col; |
---|
| 1120 | if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0)) |
---|
| 1121 | return 0; |
---|
| 1122 | if (aij->val == +1.0) |
---|
| 1123 | ; |
---|
| 1124 | else if (aij->val == -1.0) |
---|
| 1125 | b--; |
---|
| 1126 | else |
---|
| 1127 | return 0; |
---|
| 1128 | } |
---|
| 1129 | if (row->lb != (double)b) return 0; |
---|
| 1130 | return 1; |
---|
| 1131 | } |
---|
| 1132 | |
---|
| 1133 | /*********************************************************************** |
---|
| 1134 | * NAME |
---|
| 1135 | * |
---|
| 1136 | * npp_reduce_ineq_coef - reduce inequality constraint coefficients |
---|
| 1137 | * |
---|
| 1138 | * SYNOPSIS |
---|
| 1139 | * |
---|
| 1140 | * #include "glpnpp.h" |
---|
| 1141 | * int npp_reduce_ineq_coef(NPP *npp, NPPROW *row); |
---|
| 1142 | * |
---|
| 1143 | * DESCRIPTION |
---|
| 1144 | * |
---|
| 1145 | * The routine npp_reduce_ineq_coef processes specified inequality |
---|
| 1146 | * constraint attempting to replace it by an equivalent constraint, |
---|
| 1147 | * where magnitude of coefficients at binary variables is smaller than |
---|
| 1148 | * in the original constraint. If the inequality is double-sided, it is |
---|
| 1149 | * replaced by a pair of single-sided inequalities, if necessary. |
---|
| 1150 | * |
---|
| 1151 | * RETURNS |
---|
| 1152 | * |
---|
| 1153 | * The routine npp_reduce_ineq_coef returns the number of coefficients |
---|
| 1154 | * reduced. |
---|
| 1155 | * |
---|
| 1156 | * BACKGROUND |
---|
| 1157 | * |
---|
| 1158 | * Consider an inequality constraint: |
---|
| 1159 | * |
---|
| 1160 | * sum a[j] x[j] >= b. (1) |
---|
| 1161 | * j in J |
---|
| 1162 | * |
---|
| 1163 | * (In case of '<=' inequality it can be transformed to '>=' format by |
---|
| 1164 | * multiplying both its sides by -1.) Let x[k] be a binary variable; |
---|
| 1165 | * other variables can be integer as well as continuous. We can write |
---|
| 1166 | * constraint (1) as follows: |
---|
| 1167 | * |
---|
| 1168 | * a[k] x[k] + t[k] >= b, (2) |
---|
| 1169 | * |
---|
| 1170 | * where: |
---|
| 1171 | * |
---|
| 1172 | * t[k] = sum a[j] x[j]. (3) |
---|
| 1173 | * j in J\{k} |
---|
| 1174 | * |
---|
| 1175 | * Since x[k] is binary, constraint (2) is equivalent to disjunction of |
---|
| 1176 | * the following two constraints: |
---|
| 1177 | * |
---|
| 1178 | * x[k] = 0, t[k] >= b (4) |
---|
| 1179 | * |
---|
| 1180 | * OR |
---|
| 1181 | * |
---|
| 1182 | * x[k] = 1, t[k] >= b - a[k]. (5) |
---|
| 1183 | * |
---|
| 1184 | * Let also that for the partial sum t[k] be known some its implied |
---|
| 1185 | * lower bound inf t[k]. |
---|
| 1186 | * |
---|
| 1187 | * Case a[k] > 0. Let inf t[k] < b, since otherwise both constraints |
---|
| 1188 | * (4) and (5) and therefore constraint (2) are redundant. |
---|
| 1189 | * If inf t[k] > b - a[k], only constraint (5) is redundant, in which |
---|
| 1190 | * case it can be replaced with the following redundant and therefore |
---|
| 1191 | * equivalent constraint: |
---|
| 1192 | * |
---|
| 1193 | * t[k] >= b - a'[k] = inf t[k], (6) |
---|
| 1194 | * |
---|
| 1195 | * where: |
---|
| 1196 | * |
---|
| 1197 | * a'[k] = b - inf t[k]. (7) |
---|
| 1198 | * |
---|
| 1199 | * Thus, the original constraint (2) is equivalent to the following |
---|
| 1200 | * constraint with coefficient at variable x[k] changed: |
---|
| 1201 | * |
---|
| 1202 | * a'[k] x[k] + t[k] >= b. (8) |
---|
| 1203 | * |
---|
| 1204 | * From inf t[k] < b it follows that a'[k] > 0, i.e. the coefficient |
---|
| 1205 | * at x[k] keeps its sign. And from inf t[k] > b - a[k] it follows that |
---|
| 1206 | * a'[k] < a[k], i.e. the coefficient reduces in magnitude. |
---|
| 1207 | * |
---|
| 1208 | * Case a[k] < 0. Let inf t[k] < b - a[k], since otherwise both |
---|
| 1209 | * constraints (4) and (5) and therefore constraint (2) are redundant. |
---|
| 1210 | * If inf t[k] > b, only constraint (4) is redundant, in which case it |
---|
| 1211 | * can be replaced with the following redundant and therefore equivalent |
---|
| 1212 | * constraint: |
---|
| 1213 | * |
---|
| 1214 | * t[k] >= b' = inf t[k]. (9) |
---|
| 1215 | * |
---|
| 1216 | * Rewriting constraint (5) as follows: |
---|
| 1217 | * |
---|
| 1218 | * t[k] >= b - a[k] = b' - a'[k], (10) |
---|
| 1219 | * |
---|
| 1220 | * where: |
---|
| 1221 | * |
---|
| 1222 | * a'[k] = a[k] + b' - b = a[k] + inf t[k] - b, (11) |
---|
| 1223 | * |
---|
| 1224 | * we can see that disjunction of constraint (9) and (10) is equivalent |
---|
| 1225 | * to disjunction of constraint (4) and (5), from which it follows that |
---|
| 1226 | * the original constraint (2) is equivalent to the following constraint |
---|
| 1227 | * with both coefficient at variable x[k] and right-hand side changed: |
---|
| 1228 | * |
---|
| 1229 | * a'[k] x[k] + t[k] >= b'. (12) |
---|
| 1230 | * |
---|
| 1231 | * From inf t[k] < b - a[k] it follows that a'[k] < 0, i.e. the |
---|
| 1232 | * coefficient at x[k] keeps its sign. And from inf t[k] > b it follows |
---|
| 1233 | * that a'[k] > a[k], i.e. the coefficient reduces in magnitude. |
---|
| 1234 | * |
---|
| 1235 | * PROBLEM TRANSFORMATION |
---|
| 1236 | * |
---|
| 1237 | * In the routine npp_reduce_ineq_coef the following implied lower |
---|
| 1238 | * bound of the partial sum (3) is used: |
---|
| 1239 | * |
---|
| 1240 | * inf t[k] = sum a[j] l[j] + sum a[j] u[j], (13) |
---|
| 1241 | * j in Jp\{k} k in Jn\{k} |
---|
| 1242 | * |
---|
| 1243 | * where Jp = {j : a[j] > 0}, Jn = {j : a[j] < 0}, l[j] and u[j] are |
---|
| 1244 | * lower and upper bounds, resp., of variable x[j]. |
---|
| 1245 | * |
---|
| 1246 | * In order to compute inf t[k] more efficiently, the following formula, |
---|
| 1247 | * which is equivalent to (13), is actually used: |
---|
| 1248 | * |
---|
| 1249 | * ( h - a[k] l[k] = h, if a[k] > 0, |
---|
| 1250 | * inf t[k] = < (14) |
---|
| 1251 | * ( h - a[k] u[k] = h - a[k], if a[k] < 0, |
---|
| 1252 | * |
---|
| 1253 | * where: |
---|
| 1254 | * |
---|
| 1255 | * h = sum a[j] l[j] + sum a[j] u[j] (15) |
---|
| 1256 | * j in Jp j in Jn |
---|
| 1257 | * |
---|
| 1258 | * is the implied lower bound of row (1). |
---|
| 1259 | * |
---|
| 1260 | * Reduction of positive coefficient (a[k] > 0) does not change value |
---|
| 1261 | * of h, since l[k] = 0. In case of reduction of negative coefficient |
---|
| 1262 | * (a[k] < 0) from (11) it follows that: |
---|
| 1263 | * |
---|
| 1264 | * delta a[k] = a'[k] - a[k] = inf t[k] - b (> 0), (16) |
---|
| 1265 | * |
---|
| 1266 | * so new value of h (accounting that u[k] = 1) can be computed as |
---|
| 1267 | * follows: |
---|
| 1268 | * |
---|
| 1269 | * h := h + delta a[k] = h + (inf t[k] - b). (17) |
---|
| 1270 | * |
---|
| 1271 | * RECOVERING SOLUTION |
---|
| 1272 | * |
---|
| 1273 | * None needed. */ |
---|
| 1274 | |
---|
| 1275 | static int reduce_ineq_coef(NPP *npp, struct elem *ptr, double *_b) |
---|
| 1276 | { /* process inequality constraint: sum a[j] x[j] >= b */ |
---|
| 1277 | /* returns: the number of coefficients reduced */ |
---|
| 1278 | struct elem *e; |
---|
| 1279 | int count = 0; |
---|
| 1280 | double h, inf_t, new_a, b = *_b; |
---|
| 1281 | xassert(npp == npp); |
---|
| 1282 | /* compute h; see (15) */ |
---|
| 1283 | h = 0.0; |
---|
| 1284 | for (e = ptr; e != NULL; e = e->next) |
---|
| 1285 | { if (e->aj > 0.0) |
---|
| 1286 | { if (e->xj->lb == -DBL_MAX) goto done; |
---|
| 1287 | h += e->aj * e->xj->lb; |
---|
| 1288 | } |
---|
| 1289 | else /* e->aj < 0.0 */ |
---|
| 1290 | { if (e->xj->ub == +DBL_MAX) goto done; |
---|
| 1291 | h += e->aj * e->xj->ub; |
---|
| 1292 | } |
---|
| 1293 | } |
---|
| 1294 | /* perform reduction of coefficients at binary variables */ |
---|
| 1295 | for (e = ptr; e != NULL; e = e->next) |
---|
| 1296 | { /* skip non-binary variable */ |
---|
| 1297 | if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0)) |
---|
| 1298 | continue; |
---|
| 1299 | if (e->aj > 0.0) |
---|
| 1300 | { /* compute inf t[k]; see (14) */ |
---|
| 1301 | inf_t = h; |
---|
| 1302 | if (b - e->aj < inf_t && inf_t < b) |
---|
| 1303 | { /* compute reduced coefficient a'[k]; see (7) */ |
---|
| 1304 | new_a = b - inf_t; |
---|
| 1305 | if (new_a >= +1e-3 && |
---|
| 1306 | e->aj - new_a >= 0.01 * (1.0 + e->aj)) |
---|
| 1307 | { /* accept a'[k] */ |
---|
| 1308 | #ifdef GLP_DEBUG |
---|
| 1309 | xprintf("+"); |
---|
| 1310 | #endif |
---|
| 1311 | e->aj = new_a; |
---|
| 1312 | count++; |
---|
| 1313 | } |
---|
| 1314 | } |
---|
| 1315 | } |
---|
| 1316 | else /* e->aj < 0.0 */ |
---|
| 1317 | { /* compute inf t[k]; see (14) */ |
---|
| 1318 | inf_t = h - e->aj; |
---|
| 1319 | if (b < inf_t && inf_t < b - e->aj) |
---|
| 1320 | { /* compute reduced coefficient a'[k]; see (11) */ |
---|
| 1321 | new_a = e->aj + (inf_t - b); |
---|
| 1322 | if (new_a <= -1e-3 && |
---|
| 1323 | new_a - e->aj >= 0.01 * (1.0 - e->aj)) |
---|
| 1324 | { /* accept a'[k] */ |
---|
| 1325 | #ifdef GLP_DEBUG |
---|
| 1326 | xprintf("-"); |
---|
| 1327 | #endif |
---|
| 1328 | e->aj = new_a; |
---|
| 1329 | /* update h; see (17) */ |
---|
| 1330 | h += (inf_t - b); |
---|
| 1331 | /* compute b'; see (9) */ |
---|
| 1332 | b = inf_t; |
---|
| 1333 | count++; |
---|
| 1334 | } |
---|
| 1335 | } |
---|
| 1336 | } |
---|
| 1337 | } |
---|
| 1338 | *_b = b; |
---|
| 1339 | done: return count; |
---|
| 1340 | } |
---|
| 1341 | |
---|
| 1342 | int npp_reduce_ineq_coef(NPP *npp, NPPROW *row) |
---|
| 1343 | { /* reduce inequality constraint coefficients */ |
---|
| 1344 | NPPROW *copy; |
---|
| 1345 | NPPAIJ *aij; |
---|
| 1346 | struct elem *ptr, *e; |
---|
| 1347 | int kase, count[2]; |
---|
| 1348 | double b; |
---|
| 1349 | /* the row must be inequality constraint */ |
---|
| 1350 | xassert(row->lb < row->ub); |
---|
| 1351 | count[0] = count[1] = 0; |
---|
| 1352 | for (kase = 0; kase <= 1; kase++) |
---|
| 1353 | { if (kase == 0) |
---|
| 1354 | { /* process row lower bound */ |
---|
| 1355 | if (row->lb == -DBL_MAX) continue; |
---|
| 1356 | #ifdef GLP_DEBUG |
---|
| 1357 | xprintf("L"); |
---|
| 1358 | #endif |
---|
| 1359 | ptr = copy_form(npp, row, +1.0); |
---|
| 1360 | b = + row->lb; |
---|
| 1361 | } |
---|
| 1362 | else |
---|
| 1363 | { /* process row upper bound */ |
---|
| 1364 | if (row->ub == +DBL_MAX) continue; |
---|
| 1365 | #ifdef GLP_DEBUG |
---|
| 1366 | xprintf("U"); |
---|
| 1367 | #endif |
---|
| 1368 | ptr = copy_form(npp, row, -1.0); |
---|
| 1369 | b = - row->ub; |
---|
| 1370 | } |
---|
| 1371 | /* now the inequality has the form "sum a[j] x[j] >= b" */ |
---|
| 1372 | count[kase] = reduce_ineq_coef(npp, ptr, &b); |
---|
| 1373 | if (count[kase] > 0) |
---|
| 1374 | { /* the original inequality has been replaced by equivalent |
---|
| 1375 | one with coefficients reduced */ |
---|
| 1376 | if (row->lb == -DBL_MAX || row->ub == +DBL_MAX) |
---|
| 1377 | { /* the original row is single-sided inequality; no copy |
---|
| 1378 | is needed */ |
---|
| 1379 | copy = NULL; |
---|
| 1380 | } |
---|
| 1381 | else |
---|
| 1382 | { /* the original row is double-sided inequality; we need |
---|
| 1383 | to create its copy for other bound before replacing it |
---|
| 1384 | with the equivalent inequality */ |
---|
| 1385 | #ifdef GLP_DEBUG |
---|
| 1386 | xprintf("*"); |
---|
| 1387 | #endif |
---|
| 1388 | copy = npp_add_row(npp); |
---|
| 1389 | if (kase == 0) |
---|
| 1390 | { /* the copy is for upper bound */ |
---|
| 1391 | copy->lb = -DBL_MAX, copy->ub = row->ub; |
---|
| 1392 | } |
---|
| 1393 | else |
---|
| 1394 | { /* the copy is for lower bound */ |
---|
| 1395 | copy->lb = row->lb, copy->ub = +DBL_MAX; |
---|
| 1396 | } |
---|
| 1397 | /* copy original row coefficients */ |
---|
| 1398 | for (aij = row->ptr; aij != NULL; aij = aij->r_next) |
---|
| 1399 | npp_add_aij(npp, copy, aij->col, aij->val); |
---|
| 1400 | } |
---|
| 1401 | /* replace the original inequality by equivalent one */ |
---|
| 1402 | npp_erase_row(npp, row); |
---|
| 1403 | row->lb = b, row->ub = +DBL_MAX; |
---|
| 1404 | for (e = ptr; e != NULL; e = e->next) |
---|
| 1405 | npp_add_aij(npp, row, e->xj, e->aj); |
---|
| 1406 | /* continue processing upper bound for the copy */ |
---|
| 1407 | if (copy != NULL) row = copy; |
---|
| 1408 | } |
---|
| 1409 | drop_form(npp, ptr); |
---|
| 1410 | } |
---|
| 1411 | return count[0] + count[1]; |
---|
| 1412 | } |
---|
| 1413 | |
---|
| 1414 | /* eof */ |
---|