lemon-project-template-glpk
diff deps/glpk/src/glpnpp04.c @ 9:33de93886c88
Import GLPK 4.47
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 20:59:10 +0100 |
| parents | |
| children |
line diff
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000 1.2 +++ b/deps/glpk/src/glpnpp04.c Sun Nov 06 20:59:10 2011 +0100 1.3 @@ -0,0 +1,1414 @@ 1.4 +/* glpnpp04.c */ 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, 2011 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 "glpnpp.h" 1.29 + 1.30 +/*********************************************************************** 1.31 +* NAME 1.32 +* 1.33 +* npp_binarize_prob - binarize MIP problem 1.34 +* 1.35 +* SYNOPSIS 1.36 +* 1.37 +* #include "glpnpp.h" 1.38 +* int npp_binarize_prob(NPP *npp); 1.39 +* 1.40 +* DESCRIPTION 1.41 +* 1.42 +* The routine npp_binarize_prob replaces in the original MIP problem 1.43 +* every integer variable: 1.44 +* 1.45 +* l[q] <= x[q] <= u[q], (1) 1.46 +* 1.47 +* where l[q] < u[q], by an equivalent sum of binary variables. 1.48 +* 1.49 +* RETURNS 1.50 +* 1.51 +* The routine returns the number of integer variables for which the 1.52 +* transformation failed, because u[q] - l[q] > d_max. 1.53 +* 1.54 +* PROBLEM TRANSFORMATION 1.55 +* 1.56 +* If variable x[q] has non-zero lower bound, it is first processed 1.57 +* with the routine npp_lbnd_col. Thus, we can assume that: 1.58 +* 1.59 +* 0 <= x[q] <= u[q]. (2) 1.60 +* 1.61 +* If u[q] = 1, variable x[q] is already binary, so further processing 1.62 +* is not needed. Let, therefore, that 2 <= u[q] <= d_max, and n be a 1.63 +* smallest integer such that u[q] <= 2^n - 1 (n >= 2, since u[q] >= 2). 1.64 +* Then variable x[q] can be replaced by the following sum: 1.65 +* 1.66 +* n-1 1.67 +* x[q] = sum 2^k x[k], (3) 1.68 +* k=0 1.69 +* 1.70 +* where x[k] are binary columns (variables). If u[q] < 2^n - 1, the 1.71 +* following additional inequality constraint must be also included in 1.72 +* the transformed problem: 1.73 +* 1.74 +* n-1 1.75 +* sum 2^k x[k] <= u[q]. (4) 1.76 +* k=0 1.77 +* 1.78 +* Note: Assuming that in the transformed problem x[q] becomes binary 1.79 +* variable x[0], this transformation causes new n-1 binary variables 1.80 +* to appear. 1.81 +* 1.82 +* Substituting x[q] from (3) to the objective row gives: 1.83 +* 1.84 +* z = sum c[j] x[j] + c[0] = 1.85 +* j 1.86 +* 1.87 +* = sum c[j] x[j] + c[q] x[q] + c[0] = 1.88 +* j!=q 1.89 +* n-1 1.90 +* = sum c[j] x[j] + c[q] sum 2^k x[k] + c[0] = 1.91 +* j!=q k=0 1.92 +* n-1 1.93 +* = sum c[j] x[j] + sum c[k] x[k] + c[0], 1.94 +* j!=q k=0 1.95 +* 1.96 +* where: 1.97 +* 1.98 +* c[k] = 2^k c[q], k = 0, ..., n-1. (5) 1.99 +* 1.100 +* And substituting x[q] from (3) to i-th constraint row i gives: 1.101 +* 1.102 +* L[i] <= sum a[i,j] x[j] <= U[i] ==> 1.103 +* j 1.104 +* 1.105 +* L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==> 1.106 +* j!=q 1.107 +* n-1 1.108 +* L[i] <= sum a[i,j] x[j] + a[i,q] sum 2^k x[k] <= U[i] ==> 1.109 +* j!=q k=0 1.110 +* n-1 1.111 +* L[i] <= sum a[i,j] x[j] + sum a[i,k] x[k] <= U[i], 1.112 +* j!=q k=0 1.113 +* 1.114 +* where: 1.115 +* 1.116 +* a[i,k] = 2^k a[i,q], k = 0, ..., n-1. (6) 1.117 +* 1.118 +* RECOVERING SOLUTION 1.119 +* 1.120 +* Value of variable x[q] is computed with formula (3). */ 1.121 + 1.122 +struct binarize 1.123 +{ int q; 1.124 + /* column reference number for x[q] = x[0] */ 1.125 + int j; 1.126 + /* column reference number for x[1]; x[2] has reference number 1.127 + j+1, x[3] - j+2, etc. */ 1.128 + int n; 1.129 + /* total number of binary variables, n >= 2 */ 1.130 +}; 1.131 + 1.132 +static int rcv_binarize_prob(NPP *npp, void *info); 1.133 + 1.134 +int npp_binarize_prob(NPP *npp) 1.135 +{ /* binarize MIP problem */ 1.136 + struct binarize *info; 1.137 + NPPROW *row; 1.138 + NPPCOL *col, *bin; 1.139 + NPPAIJ *aij; 1.140 + int u, n, k, temp, nfails, nvars, nbins, nrows; 1.141 + /* new variables will be added to the end of the column list, so 1.142 + we go from the end to beginning of the column list */ 1.143 + nfails = nvars = nbins = nrows = 0; 1.144 + for (col = npp->c_tail; col != NULL; col = col->prev) 1.145 + { /* skip continuous variable */ 1.146 + if (!col->is_int) continue; 1.147 + /* skip fixed variable */ 1.148 + if (col->lb == col->ub) continue; 1.149 + /* skip binary variable */ 1.150 + if (col->lb == 0.0 && col->ub == 1.0) continue; 1.151 + /* check if the transformation is applicable */ 1.152 + if (col->lb < -1e6 || col->ub > +1e6 || 1.153 + col->ub - col->lb > 4095.0) 1.154 + { /* unfortunately, not */ 1.155 + nfails++; 1.156 + continue; 1.157 + } 1.158 + /* process integer non-binary variable x[q] */ 1.159 + nvars++; 1.160 + /* make x[q] non-negative, if its lower bound is non-zero */ 1.161 + if (col->lb != 0.0) 1.162 + npp_lbnd_col(npp, col); 1.163 + /* now 0 <= x[q] <= u[q] */ 1.164 + xassert(col->lb == 0.0); 1.165 + u = (int)col->ub; 1.166 + xassert(col->ub == (double)u); 1.167 + /* if x[q] is binary, further processing is not needed */ 1.168 + if (u == 1) continue; 1.169 + /* determine smallest n such that u <= 2^n - 1 (thus, n is the 1.170 + number of binary variables needed) */ 1.171 + n = 2, temp = 4; 1.172 + while (u >= temp) 1.173 + n++, temp += temp; 1.174 + nbins += n; 1.175 + /* create transformation stack entry */ 1.176 + info = npp_push_tse(npp, 1.177 + rcv_binarize_prob, sizeof(struct binarize)); 1.178 + info->q = col->j; 1.179 + info->j = 0; /* will be set below */ 1.180 + info->n = n; 1.181 + /* if u < 2^n - 1, we need one additional row for (4) */ 1.182 + if (u < temp - 1) 1.183 + { row = npp_add_row(npp), nrows++; 1.184 + row->lb = -DBL_MAX, row->ub = u; 1.185 + } 1.186 + else 1.187 + row = NULL; 1.188 + /* in the transformed problem variable x[q] becomes binary 1.189 + variable x[0], so its objective and constraint coefficients 1.190 + are not changed */ 1.191 + col->ub = 1.0; 1.192 + /* include x[0] into constraint (4) */ 1.193 + if (row != NULL) 1.194 + npp_add_aij(npp, row, col, 1.0); 1.195 + /* add other binary variables x[1], ..., x[n-1] */ 1.196 + for (k = 1, temp = 2; k < n; k++, temp += temp) 1.197 + { /* add new binary variable x[k] */ 1.198 + bin = npp_add_col(npp); 1.199 + bin->is_int = 1; 1.200 + bin->lb = 0.0, bin->ub = 1.0; 1.201 + bin->coef = (double)temp * col->coef; 1.202 + /* store column reference number for x[1] */ 1.203 + if (info->j == 0) 1.204 + info->j = bin->j; 1.205 + else 1.206 + xassert(info->j + (k-1) == bin->j); 1.207 + /* duplicate constraint coefficients for x[k]; this also 1.208 + automatically includes x[k] into constraint (4) */ 1.209 + for (aij = col->ptr; aij != NULL; aij = aij->c_next) 1.210 + npp_add_aij(npp, aij->row, bin, (double)temp * aij->val); 1.211 + } 1.212 + } 1.213 + if (nvars > 0) 1.214 + xprintf("%d integer variable(s) were replaced by %d binary one" 1.215 + "s\n", nvars, nbins); 1.216 + if (nrows > 0) 1.217 + xprintf("%d row(s) were added due to binarization\n", nrows); 1.218 + if (nfails > 0) 1.219 + xprintf("Binarization failed for %d integer variable(s)\n", 1.220 + nfails); 1.221 + return nfails; 1.222 +} 1.223 + 1.224 +static int rcv_binarize_prob(NPP *npp, void *_info) 1.225 +{ /* recovery binarized variable */ 1.226 + struct binarize *info = _info; 1.227 + int k, temp; 1.228 + double sum; 1.229 + /* compute value of x[q]; see formula (3) */ 1.230 + sum = npp->c_value[info->q]; 1.231 + for (k = 1, temp = 2; k < info->n; k++, temp += temp) 1.232 + sum += (double)temp * npp->c_value[info->j + (k-1)]; 1.233 + npp->c_value[info->q] = sum; 1.234 + return 0; 1.235 +} 1.236 + 1.237 +/**********************************************************************/ 1.238 + 1.239 +struct elem 1.240 +{ /* linear form element a[j] x[j] */ 1.241 + double aj; 1.242 + /* non-zero coefficient value */ 1.243 + NPPCOL *xj; 1.244 + /* pointer to variable (column) */ 1.245 + struct elem *next; 1.246 + /* pointer to another term */ 1.247 +}; 1.248 + 1.249 +static struct elem *copy_form(NPP *npp, NPPROW *row, double s) 1.250 +{ /* copy linear form */ 1.251 + NPPAIJ *aij; 1.252 + struct elem *ptr, *e; 1.253 + ptr = NULL; 1.254 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.255 + { e = dmp_get_atom(npp->pool, sizeof(struct elem)); 1.256 + e->aj = s * aij->val; 1.257 + e->xj = aij->col; 1.258 + e->next = ptr; 1.259 + ptr = e; 1.260 + } 1.261 + return ptr; 1.262 +} 1.263 + 1.264 +static void drop_form(NPP *npp, struct elem *ptr) 1.265 +{ /* drop linear form */ 1.266 + struct elem *e; 1.267 + while (ptr != NULL) 1.268 + { e = ptr; 1.269 + ptr = e->next; 1.270 + dmp_free_atom(npp->pool, e, sizeof(struct elem)); 1.271 + } 1.272 + return; 1.273 +} 1.274 + 1.275 +/*********************************************************************** 1.276 +* NAME 1.277 +* 1.278 +* npp_is_packing - test if constraint is packing inequality 1.279 +* 1.280 +* SYNOPSIS 1.281 +* 1.282 +* #include "glpnpp.h" 1.283 +* int npp_is_packing(NPP *npp, NPPROW *row); 1.284 +* 1.285 +* RETURNS 1.286 +* 1.287 +* If the specified row (constraint) is packing inequality (see below), 1.288 +* the routine npp_is_packing returns non-zero. Otherwise, it returns 1.289 +* zero. 1.290 +* 1.291 +* PACKING INEQUALITIES 1.292 +* 1.293 +* In canonical format the packing inequality is the following: 1.294 +* 1.295 +* sum x[j] <= 1, (1) 1.296 +* j in J 1.297 +* 1.298 +* where all variables x[j] are binary. This inequality expresses the 1.299 +* condition that in any integer feasible solution at most one variable 1.300 +* from set J can take non-zero (unity) value while other variables 1.301 +* must be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because 1.302 +* if J is empty or |J| = 1, the inequality (1) is redundant. 1.303 +* 1.304 +* In general case the packing inequality may include original variables 1.305 +* x[j] as well as their complements x~[j]: 1.306 +* 1.307 +* sum x[j] + sum x~[j] <= 1, (2) 1.308 +* j in Jp j in Jn 1.309 +* 1.310 +* where Jp and Jn are not intersected. Therefore, using substitution 1.311 +* x~[j] = 1 - x[j] gives the packing inequality in generalized format: 1.312 +* 1.313 +* sum x[j] - sum x[j] <= 1 - |Jn|. (3) 1.314 +* j in Jp j in Jn */ 1.315 + 1.316 +int npp_is_packing(NPP *npp, NPPROW *row) 1.317 +{ /* test if constraint is packing inequality */ 1.318 + NPPCOL *col; 1.319 + NPPAIJ *aij; 1.320 + int b; 1.321 + xassert(npp == npp); 1.322 + if (!(row->lb == -DBL_MAX && row->ub != +DBL_MAX)) 1.323 + return 0; 1.324 + b = 1; 1.325 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.326 + { col = aij->col; 1.327 + if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0)) 1.328 + return 0; 1.329 + if (aij->val == +1.0) 1.330 + ; 1.331 + else if (aij->val == -1.0) 1.332 + b--; 1.333 + else 1.334 + return 0; 1.335 + } 1.336 + if (row->ub != (double)b) return 0; 1.337 + return 1; 1.338 +} 1.339 + 1.340 +/*********************************************************************** 1.341 +* NAME 1.342 +* 1.343 +* npp_hidden_packing - identify hidden packing inequality 1.344 +* 1.345 +* SYNOPSIS 1.346 +* 1.347 +* #include "glpnpp.h" 1.348 +* int npp_hidden_packing(NPP *npp, NPPROW *row); 1.349 +* 1.350 +* DESCRIPTION 1.351 +* 1.352 +* The routine npp_hidden_packing processes specified inequality 1.353 +* constraint, which includes only binary variables, and the number of 1.354 +* the variables is not less than two. If the original inequality is 1.355 +* equivalent to a packing inequality, the routine replaces it by this 1.356 +* equivalent inequality. If the original constraint is double-sided 1.357 +* inequality, it is replaced by a pair of single-sided inequalities, 1.358 +* if necessary. 1.359 +* 1.360 +* RETURNS 1.361 +* 1.362 +* If the original inequality constraint was replaced by equivalent 1.363 +* packing inequality, the routine npp_hidden_packing returns non-zero. 1.364 +* Otherwise, it returns zero. 1.365 +* 1.366 +* PROBLEM TRANSFORMATION 1.367 +* 1.368 +* Consider an inequality constraint: 1.369 +* 1.370 +* sum a[j] x[j] <= b, (1) 1.371 +* j in J 1.372 +* 1.373 +* where all variables x[j] are binary, and |J| >= 2. (In case of '>=' 1.374 +* inequality it can be transformed to '<=' format by multiplying both 1.375 +* its sides by -1.) 1.376 +* 1.377 +* Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution 1.378 +* x[j] = 1 - x~[j] for all j in Jn, we have: 1.379 +* 1.380 +* sum a[j] x[j] <= b ==> 1.381 +* j in J 1.382 +* 1.383 +* sum a[j] x[j] + sum a[j] x[j] <= b ==> 1.384 +* j in Jp j in Jn 1.385 +* 1.386 +* sum a[j] x[j] + sum a[j] (1 - x~[j]) <= b ==> 1.387 +* j in Jp j in Jn 1.388 +* 1.389 +* sum a[j] x[j] - sum a[j] x~[j] <= b - sum a[j]. 1.390 +* j in Jp j in Jn j in Jn 1.391 +* 1.392 +* Thus, meaning the transformation above, we can assume that in 1.393 +* inequality (1) all coefficients a[j] are positive. Moreover, we can 1.394 +* assume that a[j] <= b. In fact, let a[j] > b; then the following 1.395 +* three cases are possible: 1.396 +* 1.397 +* 1) b < 0. In this case inequality (1) is infeasible, so the problem 1.398 +* has no feasible solution (see the routine npp_analyze_row); 1.399 +* 1.400 +* 2) b = 0. In this case inequality (1) is a forcing inequality on its 1.401 +* upper bound (see the routine npp_forcing row), from which it 1.402 +* follows that all variables x[j] should be fixed at zero; 1.403 +* 1.404 +* 3) b > 0. In this case inequality (1) defines an implied zero upper 1.405 +* bound for variable x[j] (see the routine npp_implied_bounds), from 1.406 +* which it follows that x[j] should be fixed at zero. 1.407 +* 1.408 +* It is assumed that all three cases listed above have been recognized 1.409 +* by the routine npp_process_prob, which performs basic MIP processing 1.410 +* prior to a call the routine npp_hidden_packing. So, if one of these 1.411 +* cases occurs, we should just skip processing such constraint. 1.412 +* 1.413 +* Thus, let 0 < a[j] <= b. Then it is obvious that constraint (1) is 1.414 +* equivalent to packing inquality only if: 1.415 +* 1.416 +* a[j] + a[k] > b + eps (2) 1.417 +* 1.418 +* for all j, k in J, j != k, where eps is an absolute tolerance for 1.419 +* row (linear form) value. Checking the condition (2) for all j and k, 1.420 +* j != k, requires time O(|J|^2). However, this time can be reduced to 1.421 +* O(|J|), if use minimal a[j] and a[k], in which case it is sufficient 1.422 +* to check the condition (2) only once. 1.423 +* 1.424 +* Once the original inequality (1) is replaced by equivalent packing 1.425 +* inequality, we need to perform back substitution x~[j] = 1 - x[j] for 1.426 +* all j in Jn (see above). 1.427 +* 1.428 +* RECOVERING SOLUTION 1.429 +* 1.430 +* None needed. */ 1.431 + 1.432 +static int hidden_packing(NPP *npp, struct elem *ptr, double *_b) 1.433 +{ /* process inequality constraint: sum a[j] x[j] <= b; 1.434 + 0 - specified row is NOT hidden packing inequality; 1.435 + 1 - specified row is packing inequality; 1.436 + 2 - specified row is hidden packing inequality. */ 1.437 + struct elem *e, *ej, *ek; 1.438 + int neg; 1.439 + double b = *_b, eps; 1.440 + xassert(npp == npp); 1.441 + /* a[j] must be non-zero, x[j] must be binary, for all j in J */ 1.442 + for (e = ptr; e != NULL; e = e->next) 1.443 + { xassert(e->aj != 0.0); 1.444 + xassert(e->xj->is_int); 1.445 + xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0); 1.446 + } 1.447 + /* check if the specified inequality constraint already has the 1.448 + form of packing inequality */ 1.449 + neg = 0; /* neg is |Jn| */ 1.450 + for (e = ptr; e != NULL; e = e->next) 1.451 + { if (e->aj == +1.0) 1.452 + ; 1.453 + else if (e->aj == -1.0) 1.454 + neg++; 1.455 + else 1.456 + break; 1.457 + } 1.458 + if (e == NULL) 1.459 + { /* all coefficients a[j] are +1 or -1; check rhs b */ 1.460 + if (b == (double)(1 - neg)) 1.461 + { /* it is packing inequality; no processing is needed */ 1.462 + return 1; 1.463 + } 1.464 + } 1.465 + /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j] 1.466 + positive; the result is a~[j] = |a[j]| and new rhs b */ 1.467 + for (e = ptr; e != NULL; e = e->next) 1.468 + if (e->aj < 0) b -= e->aj; 1.469 + /* now a[j] > 0 for all j in J (actually |a[j]| are used) */ 1.470 + /* if a[j] > b, skip processing--this case must not appear */ 1.471 + for (e = ptr; e != NULL; e = e->next) 1.472 + if (fabs(e->aj) > b) return 0; 1.473 + /* now 0 < a[j] <= b for all j in J */ 1.474 + /* find two minimal coefficients a[j] and a[k], j != k */ 1.475 + ej = NULL; 1.476 + for (e = ptr; e != NULL; e = e->next) 1.477 + if (ej == NULL || fabs(ej->aj) > fabs(e->aj)) ej = e; 1.478 + xassert(ej != NULL); 1.479 + ek = NULL; 1.480 + for (e = ptr; e != NULL; e = e->next) 1.481 + if (e != ej) 1.482 + if (ek == NULL || fabs(ek->aj) > fabs(e->aj)) ek = e; 1.483 + xassert(ek != NULL); 1.484 + /* the specified constraint is equivalent to packing inequality 1.485 + iff a[j] + a[k] > b + eps */ 1.486 + eps = 1e-3 + 1e-6 * fabs(b); 1.487 + if (fabs(ej->aj) + fabs(ek->aj) <= b + eps) return 0; 1.488 + /* perform back substitution x~[j] = 1 - x[j] and construct the 1.489 + final equivalent packing inequality in generalized format */ 1.490 + b = 1.0; 1.491 + for (e = ptr; e != NULL; e = e->next) 1.492 + { if (e->aj > 0.0) 1.493 + e->aj = +1.0; 1.494 + else /* e->aj < 0.0 */ 1.495 + e->aj = -1.0, b -= 1.0; 1.496 + } 1.497 + *_b = b; 1.498 + return 2; 1.499 +} 1.500 + 1.501 +int npp_hidden_packing(NPP *npp, NPPROW *row) 1.502 +{ /* identify hidden packing inequality */ 1.503 + NPPROW *copy; 1.504 + NPPAIJ *aij; 1.505 + struct elem *ptr, *e; 1.506 + int kase, ret, count = 0; 1.507 + double b; 1.508 + /* the row must be inequality constraint */ 1.509 + xassert(row->lb < row->ub); 1.510 + for (kase = 0; kase <= 1; kase++) 1.511 + { if (kase == 0) 1.512 + { /* process row upper bound */ 1.513 + if (row->ub == +DBL_MAX) continue; 1.514 + ptr = copy_form(npp, row, +1.0); 1.515 + b = + row->ub; 1.516 + } 1.517 + else 1.518 + { /* process row lower bound */ 1.519 + if (row->lb == -DBL_MAX) continue; 1.520 + ptr = copy_form(npp, row, -1.0); 1.521 + b = - row->lb; 1.522 + } 1.523 + /* now the inequality has the form "sum a[j] x[j] <= b" */ 1.524 + ret = hidden_packing(npp, ptr, &b); 1.525 + xassert(0 <= ret && ret <= 2); 1.526 + if (kase == 1 && ret == 1 || ret == 2) 1.527 + { /* the original inequality has been identified as hidden 1.528 + packing inequality */ 1.529 + count++; 1.530 +#ifdef GLP_DEBUG 1.531 + xprintf("Original constraint:\n"); 1.532 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.533 + xprintf(" %+g x%d", aij->val, aij->col->j); 1.534 + if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb); 1.535 + if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub); 1.536 + xprintf("\n"); 1.537 + xprintf("Equivalent packing inequality:\n"); 1.538 + for (e = ptr; e != NULL; e = e->next) 1.539 + xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j); 1.540 + xprintf(", <= %g\n", b); 1.541 +#endif 1.542 + if (row->lb == -DBL_MAX || row->ub == +DBL_MAX) 1.543 + { /* the original row is single-sided inequality; no copy 1.544 + is needed */ 1.545 + copy = NULL; 1.546 + } 1.547 + else 1.548 + { /* the original row is double-sided inequality; we need 1.549 + to create its copy for other bound before replacing it 1.550 + with the equivalent inequality */ 1.551 + copy = npp_add_row(npp); 1.552 + if (kase == 0) 1.553 + { /* the copy is for lower bound */ 1.554 + copy->lb = row->lb, copy->ub = +DBL_MAX; 1.555 + } 1.556 + else 1.557 + { /* the copy is for upper bound */ 1.558 + copy->lb = -DBL_MAX, copy->ub = row->ub; 1.559 + } 1.560 + /* copy original row coefficients */ 1.561 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.562 + npp_add_aij(npp, copy, aij->col, aij->val); 1.563 + } 1.564 + /* replace the original inequality by equivalent one */ 1.565 + npp_erase_row(npp, row); 1.566 + row->lb = -DBL_MAX, row->ub = b; 1.567 + for (e = ptr; e != NULL; e = e->next) 1.568 + npp_add_aij(npp, row, e->xj, e->aj); 1.569 + /* continue processing lower bound for the copy */ 1.570 + if (copy != NULL) row = copy; 1.571 + } 1.572 + drop_form(npp, ptr); 1.573 + } 1.574 + return count; 1.575 +} 1.576 + 1.577 +/*********************************************************************** 1.578 +* NAME 1.579 +* 1.580 +* npp_implied_packing - identify implied packing inequality 1.581 +* 1.582 +* SYNOPSIS 1.583 +* 1.584 +* #include "glpnpp.h" 1.585 +* int npp_implied_packing(NPP *npp, NPPROW *row, int which, 1.586 +* NPPCOL *var[], char set[]); 1.587 +* 1.588 +* DESCRIPTION 1.589 +* 1.590 +* The routine npp_implied_packing processes specified row (constraint) 1.591 +* of general format: 1.592 +* 1.593 +* L <= sum a[j] x[j] <= U. (1) 1.594 +* j 1.595 +* 1.596 +* If which = 0, only lower bound L, which must exist, is considered, 1.597 +* while upper bound U is ignored. Similarly, if which = 1, only upper 1.598 +* bound U, which must exist, is considered, while lower bound L is 1.599 +* ignored. Thus, if the specified row is a double-sided inequality or 1.600 +* equality constraint, this routine should be called twice for both 1.601 +* lower and upper bounds. 1.602 +* 1.603 +* The routine npp_implied_packing attempts to find a non-trivial (i.e. 1.604 +* having not less than two binary variables) packing inequality: 1.605 +* 1.606 +* sum x[j] - sum x[j] <= 1 - |Jn|, (2) 1.607 +* j in Jp j in Jn 1.608 +* 1.609 +* which is relaxation of the constraint (1) in the sense that any 1.610 +* solution satisfying to that constraint also satisfies to the packing 1.611 +* inequality (2). If such relaxation exists, the routine stores 1.612 +* pointers to descriptors of corresponding binary variables and their 1.613 +* flags, resp., to locations var[1], var[2], ..., var[len] and set[1], 1.614 +* set[2], ..., set[len], where set[j] = 0 means that j in Jp and 1.615 +* set[j] = 1 means that j in Jn. 1.616 +* 1.617 +* RETURNS 1.618 +* 1.619 +* The routine npp_implied_packing returns len, which is the total 1.620 +* number of binary variables in the packing inequality found, len >= 2. 1.621 +* However, if the relaxation does not exist, the routine returns zero. 1.622 +* 1.623 +* ALGORITHM 1.624 +* 1.625 +* If which = 0, the constraint coefficients (1) are multiplied by -1 1.626 +* and b is assigned -L; if which = 1, the constraint coefficients (1) 1.627 +* are not changed and b is assigned +U. In both cases the specified 1.628 +* constraint gets the following format: 1.629 +* 1.630 +* sum a[j] x[j] <= b. (3) 1.631 +* j 1.632 +* 1.633 +* (Note that (3) is a relaxation of (1), because one of bounds L or U 1.634 +* is ignored.) 1.635 +* 1.636 +* Let J be set of binary variables, Kp be set of non-binary (integer 1.637 +* or continuous) variables with a[j] > 0, and Kn be set of non-binary 1.638 +* variables with a[j] < 0. Then the inequality (3) can be written as 1.639 +* follows: 1.640 +* 1.641 +* sum a[j] x[j] <= b - sum a[j] x[j] - sum a[j] x[j]. (4) 1.642 +* j in J j in Kp j in Kn 1.643 +* 1.644 +* To get rid of non-binary variables we can replace the inequality (4) 1.645 +* by the following relaxed inequality: 1.646 +* 1.647 +* sum a[j] x[j] <= b~, (5) 1.648 +* j in J 1.649 +* 1.650 +* where: 1.651 +* 1.652 +* b~ = sup(b - sum a[j] x[j] - sum a[j] x[j]) = 1.653 +* j in Kp j in Kn 1.654 +* 1.655 +* = b - inf sum a[j] x[j] - inf sum a[j] x[j] = (6) 1.656 +* j in Kp j in Kn 1.657 +* 1.658 +* = b - sum a[j] l[j] - sum a[j] u[j]. 1.659 +* j in Kp j in Kn 1.660 +* 1.661 +* Note that if lower bound l[j] (if j in Kp) or upper bound u[j] 1.662 +* (if j in Kn) of some non-binary variable x[j] does not exist, then 1.663 +* formally b = +oo, in which case further analysis is not performed. 1.664 +* 1.665 +* Let Bp = {j in J: a[j] > 0}, Bn = {j in J: a[j] < 0}. To make all 1.666 +* the inequality coefficients in (5) positive, we replace all x[j] in 1.667 +* Bn by their complementaries, substituting x[j] = 1 - x~[j] for all 1.668 +* j in Bn, that gives: 1.669 +* 1.670 +* sum a[j] x[j] - sum a[j] x~[j] <= b~ - sum a[j]. (7) 1.671 +* j in Bp j in Bn j in Bn 1.672 +* 1.673 +* This inequality is a relaxation of the original constraint (1), and 1.674 +* it is a binary knapsack inequality. Writing it in the standard format 1.675 +* we have: 1.676 +* 1.677 +* sum alfa[j] z[j] <= beta, (8) 1.678 +* j in J 1.679 +* 1.680 +* where: 1.681 +* ( + a[j], if j in Bp, 1.682 +* alfa[j] = < (9) 1.683 +* ( - a[j], if j in Bn, 1.684 +* 1.685 +* ( x[j], if j in Bp, 1.686 +* z[j] = < (10) 1.687 +* ( 1 - x[j], if j in Bn, 1.688 +* 1.689 +* beta = b~ - sum a[j]. (11) 1.690 +* j in Bn 1.691 +* 1.692 +* In the inequality (8) all coefficients are positive, therefore, the 1.693 +* packing relaxation to be found for this inequality is the following: 1.694 +* 1.695 +* sum z[j] <= 1. (12) 1.696 +* j in P 1.697 +* 1.698 +* It is obvious that set P within J, which we would like to find, must 1.699 +* satisfy to the following condition: 1.700 +* 1.701 +* alfa[j] + alfa[k] > beta + eps for all j, k in P, j != k, (13) 1.702 +* 1.703 +* where eps is an absolute tolerance for value of the linear form. 1.704 +* Thus, it is natural to take P = {j: alpha[j] > (beta + eps) / 2}. 1.705 +* Moreover, if in the equality (8) there exist coefficients alfa[k], 1.706 +* for which alfa[k] <= (beta + eps) / 2, but which, nevertheless, 1.707 +* satisfies to the condition (13) for all j in P, *one* corresponding 1.708 +* variable z[k] (having, for example, maximal coefficient alfa[k]) can 1.709 +* be included in set P, that allows increasing the number of binary 1.710 +* variables in (12) by one. 1.711 +* 1.712 +* Once the set P has been built, for the inequality (12) we need to 1.713 +* perform back substitution according to (10) in order to express it 1.714 +* through the original binary variables. As the result of such back 1.715 +* substitution the relaxed packing inequality get its final format (2), 1.716 +* where Jp = J intersect Bp, and Jn = J intersect Bn. */ 1.717 + 1.718 +int npp_implied_packing(NPP *npp, NPPROW *row, int which, 1.719 + NPPCOL *var[], char set[]) 1.720 +{ struct elem *ptr, *e, *i, *k; 1.721 + int len = 0; 1.722 + double b, eps; 1.723 + /* build inequality (3) */ 1.724 + if (which == 0) 1.725 + { ptr = copy_form(npp, row, -1.0); 1.726 + xassert(row->lb != -DBL_MAX); 1.727 + b = - row->lb; 1.728 + } 1.729 + else if (which == 1) 1.730 + { ptr = copy_form(npp, row, +1.0); 1.731 + xassert(row->ub != +DBL_MAX); 1.732 + b = + row->ub; 1.733 + } 1.734 + /* remove non-binary variables to build relaxed inequality (5); 1.735 + compute its right-hand side b~ with formula (6) */ 1.736 + for (e = ptr; e != NULL; e = e->next) 1.737 + { if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0)) 1.738 + { /* x[j] is non-binary variable */ 1.739 + if (e->aj > 0.0) 1.740 + { if (e->xj->lb == -DBL_MAX) goto done; 1.741 + b -= e->aj * e->xj->lb; 1.742 + } 1.743 + else /* e->aj < 0.0 */ 1.744 + { if (e->xj->ub == +DBL_MAX) goto done; 1.745 + b -= e->aj * e->xj->ub; 1.746 + } 1.747 + /* a[j] = 0 means that variable x[j] is removed */ 1.748 + e->aj = 0.0; 1.749 + } 1.750 + } 1.751 + /* substitute x[j] = 1 - x~[j] to build knapsack inequality (8); 1.752 + compute its right-hand side beta with formula (11) */ 1.753 + for (e = ptr; e != NULL; e = e->next) 1.754 + if (e->aj < 0.0) b -= e->aj; 1.755 + /* if beta is close to zero, the knapsack inequality is either 1.756 + infeasible or forcing inequality; this must never happen, so 1.757 + we skip further analysis */ 1.758 + if (b < 1e-3) goto done; 1.759 + /* build set P as well as sets Jp and Jn, and determine x[k] as 1.760 + explained above in comments to the routine */ 1.761 + eps = 1e-3 + 1e-6 * b; 1.762 + i = k = NULL; 1.763 + for (e = ptr; e != NULL; e = e->next) 1.764 + { /* note that alfa[j] = |a[j]| */ 1.765 + if (fabs(e->aj) > 0.5 * (b + eps)) 1.766 + { /* alfa[j] > (b + eps) / 2; include x[j] in set P, i.e. in 1.767 + set Jp or Jn */ 1.768 + var[++len] = e->xj; 1.769 + set[len] = (char)(e->aj > 0.0 ? 0 : 1); 1.770 + /* alfa[i] = min alfa[j] over all j included in set P */ 1.771 + if (i == NULL || fabs(i->aj) > fabs(e->aj)) i = e; 1.772 + } 1.773 + else if (fabs(e->aj) >= 1e-3) 1.774 + { /* alfa[k] = max alfa[j] over all j not included in set P; 1.775 + we skip coefficient a[j] if it is close to zero to avoid 1.776 + numerically unreliable results */ 1.777 + if (k == NULL || fabs(k->aj) < fabs(e->aj)) k = e; 1.778 + } 1.779 + } 1.780 + /* if alfa[k] satisfies to condition (13) for all j in P, include 1.781 + x[k] in P */ 1.782 + if (i != NULL && k != NULL && fabs(i->aj) + fabs(k->aj) > b + eps) 1.783 + { var[++len] = k->xj; 1.784 + set[len] = (char)(k->aj > 0.0 ? 0 : 1); 1.785 + } 1.786 + /* trivial packing inequality being redundant must never appear, 1.787 + so we just ignore it */ 1.788 + if (len < 2) len = 0; 1.789 +done: drop_form(npp, ptr); 1.790 + return len; 1.791 +} 1.792 + 1.793 +/*********************************************************************** 1.794 +* NAME 1.795 +* 1.796 +* npp_is_covering - test if constraint is covering inequality 1.797 +* 1.798 +* SYNOPSIS 1.799 +* 1.800 +* #include "glpnpp.h" 1.801 +* int npp_is_covering(NPP *npp, NPPROW *row); 1.802 +* 1.803 +* RETURNS 1.804 +* 1.805 +* If the specified row (constraint) is covering inequality (see below), 1.806 +* the routine npp_is_covering returns non-zero. Otherwise, it returns 1.807 +* zero. 1.808 +* 1.809 +* COVERING INEQUALITIES 1.810 +* 1.811 +* In canonical format the covering inequality is the following: 1.812 +* 1.813 +* sum x[j] >= 1, (1) 1.814 +* j in J 1.815 +* 1.816 +* where all variables x[j] are binary. This inequality expresses the 1.817 +* condition that in any integer feasible solution variables in set J 1.818 +* cannot be all equal to zero at the same time, i.e. at least one 1.819 +* variable must take non-zero (unity) value. W.l.o.g. it is assumed 1.820 +* that |J| >= 2, because if J is empty, the inequality (1) is 1.821 +* infeasible, and if |J| = 1, the inequality (1) is a forcing row. 1.822 +* 1.823 +* In general case the covering inequality may include original 1.824 +* variables x[j] as well as their complements x~[j]: 1.825 +* 1.826 +* sum x[j] + sum x~[j] >= 1, (2) 1.827 +* j in Jp j in Jn 1.828 +* 1.829 +* where Jp and Jn are not intersected. Therefore, using substitution 1.830 +* x~[j] = 1 - x[j] gives the packing inequality in generalized format: 1.831 +* 1.832 +* sum x[j] - sum x[j] >= 1 - |Jn|. (3) 1.833 +* j in Jp j in Jn 1.834 +* 1.835 +* (May note that the inequality (3) cuts off infeasible solutions, 1.836 +* where x[j] = 0 for all j in Jp and x[j] = 1 for all j in Jn.) 1.837 +* 1.838 +* NOTE: If |J| = 2, the inequality (3) is equivalent to packing 1.839 +* inequality (see the routine npp_is_packing). */ 1.840 + 1.841 +int npp_is_covering(NPP *npp, NPPROW *row) 1.842 +{ /* test if constraint is covering inequality */ 1.843 + NPPCOL *col; 1.844 + NPPAIJ *aij; 1.845 + int b; 1.846 + xassert(npp == npp); 1.847 + if (!(row->lb != -DBL_MAX && row->ub == +DBL_MAX)) 1.848 + return 0; 1.849 + b = 1; 1.850 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.851 + { col = aij->col; 1.852 + if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0)) 1.853 + return 0; 1.854 + if (aij->val == +1.0) 1.855 + ; 1.856 + else if (aij->val == -1.0) 1.857 + b--; 1.858 + else 1.859 + return 0; 1.860 + } 1.861 + if (row->lb != (double)b) return 0; 1.862 + return 1; 1.863 +} 1.864 + 1.865 +/*********************************************************************** 1.866 +* NAME 1.867 +* 1.868 +* npp_hidden_covering - identify hidden covering inequality 1.869 +* 1.870 +* SYNOPSIS 1.871 +* 1.872 +* #include "glpnpp.h" 1.873 +* int npp_hidden_covering(NPP *npp, NPPROW *row); 1.874 +* 1.875 +* DESCRIPTION 1.876 +* 1.877 +* The routine npp_hidden_covering processes specified inequality 1.878 +* constraint, which includes only binary variables, and the number of 1.879 +* the variables is not less than three. If the original inequality is 1.880 +* equivalent to a covering inequality (see below), the routine 1.881 +* replaces it by the equivalent inequality. If the original constraint 1.882 +* is double-sided inequality, it is replaced by a pair of single-sided 1.883 +* inequalities, if necessary. 1.884 +* 1.885 +* RETURNS 1.886 +* 1.887 +* If the original inequality constraint was replaced by equivalent 1.888 +* covering inequality, the routine npp_hidden_covering returns 1.889 +* non-zero. Otherwise, it returns zero. 1.890 +* 1.891 +* PROBLEM TRANSFORMATION 1.892 +* 1.893 +* Consider an inequality constraint: 1.894 +* 1.895 +* sum a[j] x[j] >= b, (1) 1.896 +* j in J 1.897 +* 1.898 +* where all variables x[j] are binary, and |J| >= 3. (In case of '<=' 1.899 +* inequality it can be transformed to '>=' format by multiplying both 1.900 +* its sides by -1.) 1.901 +* 1.902 +* Let Jp = {j: a[j] > 0}, Jn = {j: a[j] < 0}. Performing substitution 1.903 +* x[j] = 1 - x~[j] for all j in Jn, we have: 1.904 +* 1.905 +* sum a[j] x[j] >= b ==> 1.906 +* j in J 1.907 +* 1.908 +* sum a[j] x[j] + sum a[j] x[j] >= b ==> 1.909 +* j in Jp j in Jn 1.910 +* 1.911 +* sum a[j] x[j] + sum a[j] (1 - x~[j]) >= b ==> 1.912 +* j in Jp j in Jn 1.913 +* 1.914 +* sum m a[j] x[j] - sum a[j] x~[j] >= b - sum a[j]. 1.915 +* j in Jp j in Jn j in Jn 1.916 +* 1.917 +* Thus, meaning the transformation above, we can assume that in 1.918 +* inequality (1) all coefficients a[j] are positive. Moreover, we can 1.919 +* assume that b > 0, because otherwise the inequality (1) would be 1.920 +* redundant (see the routine npp_analyze_row). It is then obvious that 1.921 +* constraint (1) is equivalent to covering inequality only if: 1.922 +* 1.923 +* a[j] >= b, (2) 1.924 +* 1.925 +* for all j in J. 1.926 +* 1.927 +* Once the original inequality (1) is replaced by equivalent covering 1.928 +* inequality, we need to perform back substitution x~[j] = 1 - x[j] for 1.929 +* all j in Jn (see above). 1.930 +* 1.931 +* RECOVERING SOLUTION 1.932 +* 1.933 +* None needed. */ 1.934 + 1.935 +static int hidden_covering(NPP *npp, struct elem *ptr, double *_b) 1.936 +{ /* process inequality constraint: sum a[j] x[j] >= b; 1.937 + 0 - specified row is NOT hidden covering inequality; 1.938 + 1 - specified row is covering inequality; 1.939 + 2 - specified row is hidden covering inequality. */ 1.940 + struct elem *e; 1.941 + int neg; 1.942 + double b = *_b, eps; 1.943 + xassert(npp == npp); 1.944 + /* a[j] must be non-zero, x[j] must be binary, for all j in J */ 1.945 + for (e = ptr; e != NULL; e = e->next) 1.946 + { xassert(e->aj != 0.0); 1.947 + xassert(e->xj->is_int); 1.948 + xassert(e->xj->lb == 0.0 && e->xj->ub == 1.0); 1.949 + } 1.950 + /* check if the specified inequality constraint already has the 1.951 + form of covering inequality */ 1.952 + neg = 0; /* neg is |Jn| */ 1.953 + for (e = ptr; e != NULL; e = e->next) 1.954 + { if (e->aj == +1.0) 1.955 + ; 1.956 + else if (e->aj == -1.0) 1.957 + neg++; 1.958 + else 1.959 + break; 1.960 + } 1.961 + if (e == NULL) 1.962 + { /* all coefficients a[j] are +1 or -1; check rhs b */ 1.963 + if (b == (double)(1 - neg)) 1.964 + { /* it is covering inequality; no processing is needed */ 1.965 + return 1; 1.966 + } 1.967 + } 1.968 + /* substitute x[j] = 1 - x~[j] for all j in Jn to make all a[j] 1.969 + positive; the result is a~[j] = |a[j]| and new rhs b */ 1.970 + for (e = ptr; e != NULL; e = e->next) 1.971 + if (e->aj < 0) b -= e->aj; 1.972 + /* now a[j] > 0 for all j in J (actually |a[j]| are used) */ 1.973 + /* if b <= 0, skip processing--this case must not appear */ 1.974 + if (b < 1e-3) return 0; 1.975 + /* now a[j] > 0 for all j in J, and b > 0 */ 1.976 + /* the specified constraint is equivalent to covering inequality 1.977 + iff a[j] >= b for all j in J */ 1.978 + eps = 1e-9 + 1e-12 * fabs(b); 1.979 + for (e = ptr; e != NULL; e = e->next) 1.980 + if (fabs(e->aj) < b - eps) return 0; 1.981 + /* perform back substitution x~[j] = 1 - x[j] and construct the 1.982 + final equivalent covering inequality in generalized format */ 1.983 + b = 1.0; 1.984 + for (e = ptr; e != NULL; e = e->next) 1.985 + { if (e->aj > 0.0) 1.986 + e->aj = +1.0; 1.987 + else /* e->aj < 0.0 */ 1.988 + e->aj = -1.0, b -= 1.0; 1.989 + } 1.990 + *_b = b; 1.991 + return 2; 1.992 +} 1.993 + 1.994 +int npp_hidden_covering(NPP *npp, NPPROW *row) 1.995 +{ /* identify hidden covering inequality */ 1.996 + NPPROW *copy; 1.997 + NPPAIJ *aij; 1.998 + struct elem *ptr, *e; 1.999 + int kase, ret, count = 0; 1.1000 + double b; 1.1001 + /* the row must be inequality constraint */ 1.1002 + xassert(row->lb < row->ub); 1.1003 + for (kase = 0; kase <= 1; kase++) 1.1004 + { if (kase == 0) 1.1005 + { /* process row lower bound */ 1.1006 + if (row->lb == -DBL_MAX) continue; 1.1007 + ptr = copy_form(npp, row, +1.0); 1.1008 + b = + row->lb; 1.1009 + } 1.1010 + else 1.1011 + { /* process row upper bound */ 1.1012 + if (row->ub == +DBL_MAX) continue; 1.1013 + ptr = copy_form(npp, row, -1.0); 1.1014 + b = - row->ub; 1.1015 + } 1.1016 + /* now the inequality has the form "sum a[j] x[j] >= b" */ 1.1017 + ret = hidden_covering(npp, ptr, &b); 1.1018 + xassert(0 <= ret && ret <= 2); 1.1019 + if (kase == 1 && ret == 1 || ret == 2) 1.1020 + { /* the original inequality has been identified as hidden 1.1021 + covering inequality */ 1.1022 + count++; 1.1023 +#ifdef GLP_DEBUG 1.1024 + xprintf("Original constraint:\n"); 1.1025 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.1026 + xprintf(" %+g x%d", aij->val, aij->col->j); 1.1027 + if (row->lb != -DBL_MAX) xprintf(", >= %g", row->lb); 1.1028 + if (row->ub != +DBL_MAX) xprintf(", <= %g", row->ub); 1.1029 + xprintf("\n"); 1.1030 + xprintf("Equivalent covering inequality:\n"); 1.1031 + for (e = ptr; e != NULL; e = e->next) 1.1032 + xprintf(" %sx%d", e->aj > 0.0 ? "+" : "-", e->xj->j); 1.1033 + xprintf(", >= %g\n", b); 1.1034 +#endif 1.1035 + if (row->lb == -DBL_MAX || row->ub == +DBL_MAX) 1.1036 + { /* the original row is single-sided inequality; no copy 1.1037 + is needed */ 1.1038 + copy = NULL; 1.1039 + } 1.1040 + else 1.1041 + { /* the original row is double-sided inequality; we need 1.1042 + to create its copy for other bound before replacing it 1.1043 + with the equivalent inequality */ 1.1044 + copy = npp_add_row(npp); 1.1045 + if (kase == 0) 1.1046 + { /* the copy is for upper bound */ 1.1047 + copy->lb = -DBL_MAX, copy->ub = row->ub; 1.1048 + } 1.1049 + else 1.1050 + { /* the copy is for lower bound */ 1.1051 + copy->lb = row->lb, copy->ub = +DBL_MAX; 1.1052 + } 1.1053 + /* copy original row coefficients */ 1.1054 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.1055 + npp_add_aij(npp, copy, aij->col, aij->val); 1.1056 + } 1.1057 + /* replace the original inequality by equivalent one */ 1.1058 + npp_erase_row(npp, row); 1.1059 + row->lb = b, row->ub = +DBL_MAX; 1.1060 + for (e = ptr; e != NULL; e = e->next) 1.1061 + npp_add_aij(npp, row, e->xj, e->aj); 1.1062 + /* continue processing upper bound for the copy */ 1.1063 + if (copy != NULL) row = copy; 1.1064 + } 1.1065 + drop_form(npp, ptr); 1.1066 + } 1.1067 + return count; 1.1068 +} 1.1069 + 1.1070 +/*********************************************************************** 1.1071 +* NAME 1.1072 +* 1.1073 +* npp_is_partitioning - test if constraint is partitioning equality 1.1074 +* 1.1075 +* SYNOPSIS 1.1076 +* 1.1077 +* #include "glpnpp.h" 1.1078 +* int npp_is_partitioning(NPP *npp, NPPROW *row); 1.1079 +* 1.1080 +* RETURNS 1.1081 +* 1.1082 +* If the specified row (constraint) is partitioning equality (see 1.1083 +* below), the routine npp_is_partitioning returns non-zero. Otherwise, 1.1084 +* it returns zero. 1.1085 +* 1.1086 +* PARTITIONING EQUALITIES 1.1087 +* 1.1088 +* In canonical format the partitioning equality is the following: 1.1089 +* 1.1090 +* sum x[j] = 1, (1) 1.1091 +* j in J 1.1092 +* 1.1093 +* where all variables x[j] are binary. This equality expresses the 1.1094 +* condition that in any integer feasible solution exactly one variable 1.1095 +* in set J must take non-zero (unity) value while other variables must 1.1096 +* be equal to zero. W.l.o.g. it is assumed that |J| >= 2, because if 1.1097 +* J is empty, the inequality (1) is infeasible, and if |J| = 1, the 1.1098 +* inequality (1) is a fixing row. 1.1099 +* 1.1100 +* In general case the partitioning equality may include original 1.1101 +* variables x[j] as well as their complements x~[j]: 1.1102 +* 1.1103 +* sum x[j] + sum x~[j] = 1, (2) 1.1104 +* j in Jp j in Jn 1.1105 +* 1.1106 +* where Jp and Jn are not intersected. Therefore, using substitution 1.1107 +* x~[j] = 1 - x[j] leads to the partitioning equality in generalized 1.1108 +* format: 1.1109 +* 1.1110 +* sum x[j] - sum x[j] = 1 - |Jn|. (3) 1.1111 +* j in Jp j in Jn */ 1.1112 + 1.1113 +int npp_is_partitioning(NPP *npp, NPPROW *row) 1.1114 +{ /* test if constraint is partitioning equality */ 1.1115 + NPPCOL *col; 1.1116 + NPPAIJ *aij; 1.1117 + int b; 1.1118 + xassert(npp == npp); 1.1119 + if (row->lb != row->ub) return 0; 1.1120 + b = 1; 1.1121 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.1122 + { col = aij->col; 1.1123 + if (!(col->is_int && col->lb == 0.0 && col->ub == 1.0)) 1.1124 + return 0; 1.1125 + if (aij->val == +1.0) 1.1126 + ; 1.1127 + else if (aij->val == -1.0) 1.1128 + b--; 1.1129 + else 1.1130 + return 0; 1.1131 + } 1.1132 + if (row->lb != (double)b) return 0; 1.1133 + return 1; 1.1134 +} 1.1135 + 1.1136 +/*********************************************************************** 1.1137 +* NAME 1.1138 +* 1.1139 +* npp_reduce_ineq_coef - reduce inequality constraint coefficients 1.1140 +* 1.1141 +* SYNOPSIS 1.1142 +* 1.1143 +* #include "glpnpp.h" 1.1144 +* int npp_reduce_ineq_coef(NPP *npp, NPPROW *row); 1.1145 +* 1.1146 +* DESCRIPTION 1.1147 +* 1.1148 +* The routine npp_reduce_ineq_coef processes specified inequality 1.1149 +* constraint attempting to replace it by an equivalent constraint, 1.1150 +* where magnitude of coefficients at binary variables is smaller than 1.1151 +* in the original constraint. If the inequality is double-sided, it is 1.1152 +* replaced by a pair of single-sided inequalities, if necessary. 1.1153 +* 1.1154 +* RETURNS 1.1155 +* 1.1156 +* The routine npp_reduce_ineq_coef returns the number of coefficients 1.1157 +* reduced. 1.1158 +* 1.1159 +* BACKGROUND 1.1160 +* 1.1161 +* Consider an inequality constraint: 1.1162 +* 1.1163 +* sum a[j] x[j] >= b. (1) 1.1164 +* j in J 1.1165 +* 1.1166 +* (In case of '<=' inequality it can be transformed to '>=' format by 1.1167 +* multiplying both its sides by -1.) Let x[k] be a binary variable; 1.1168 +* other variables can be integer as well as continuous. We can write 1.1169 +* constraint (1) as follows: 1.1170 +* 1.1171 +* a[k] x[k] + t[k] >= b, (2) 1.1172 +* 1.1173 +* where: 1.1174 +* 1.1175 +* t[k] = sum a[j] x[j]. (3) 1.1176 +* j in J\{k} 1.1177 +* 1.1178 +* Since x[k] is binary, constraint (2) is equivalent to disjunction of 1.1179 +* the following two constraints: 1.1180 +* 1.1181 +* x[k] = 0, t[k] >= b (4) 1.1182 +* 1.1183 +* OR 1.1184 +* 1.1185 +* x[k] = 1, t[k] >= b - a[k]. (5) 1.1186 +* 1.1187 +* Let also that for the partial sum t[k] be known some its implied 1.1188 +* lower bound inf t[k]. 1.1189 +* 1.1190 +* Case a[k] > 0. Let inf t[k] < b, since otherwise both constraints 1.1191 +* (4) and (5) and therefore constraint (2) are redundant. 1.1192 +* If inf t[k] > b - a[k], only constraint (5) is redundant, in which 1.1193 +* case it can be replaced with the following redundant and therefore 1.1194 +* equivalent constraint: 1.1195 +* 1.1196 +* t[k] >= b - a'[k] = inf t[k], (6) 1.1197 +* 1.1198 +* where: 1.1199 +* 1.1200 +* a'[k] = b - inf t[k]. (7) 1.1201 +* 1.1202 +* Thus, the original constraint (2) is equivalent to the following 1.1203 +* constraint with coefficient at variable x[k] changed: 1.1204 +* 1.1205 +* a'[k] x[k] + t[k] >= b. (8) 1.1206 +* 1.1207 +* From inf t[k] < b it follows that a'[k] > 0, i.e. the coefficient 1.1208 +* at x[k] keeps its sign. And from inf t[k] > b - a[k] it follows that 1.1209 +* a'[k] < a[k], i.e. the coefficient reduces in magnitude. 1.1210 +* 1.1211 +* Case a[k] < 0. Let inf t[k] < b - a[k], since otherwise both 1.1212 +* constraints (4) and (5) and therefore constraint (2) are redundant. 1.1213 +* If inf t[k] > b, only constraint (4) is redundant, in which case it 1.1214 +* can be replaced with the following redundant and therefore equivalent 1.1215 +* constraint: 1.1216 +* 1.1217 +* t[k] >= b' = inf t[k]. (9) 1.1218 +* 1.1219 +* Rewriting constraint (5) as follows: 1.1220 +* 1.1221 +* t[k] >= b - a[k] = b' - a'[k], (10) 1.1222 +* 1.1223 +* where: 1.1224 +* 1.1225 +* a'[k] = a[k] + b' - b = a[k] + inf t[k] - b, (11) 1.1226 +* 1.1227 +* we can see that disjunction of constraint (9) and (10) is equivalent 1.1228 +* to disjunction of constraint (4) and (5), from which it follows that 1.1229 +* the original constraint (2) is equivalent to the following constraint 1.1230 +* with both coefficient at variable x[k] and right-hand side changed: 1.1231 +* 1.1232 +* a'[k] x[k] + t[k] >= b'. (12) 1.1233 +* 1.1234 +* From inf t[k] < b - a[k] it follows that a'[k] < 0, i.e. the 1.1235 +* coefficient at x[k] keeps its sign. And from inf t[k] > b it follows 1.1236 +* that a'[k] > a[k], i.e. the coefficient reduces in magnitude. 1.1237 +* 1.1238 +* PROBLEM TRANSFORMATION 1.1239 +* 1.1240 +* In the routine npp_reduce_ineq_coef the following implied lower 1.1241 +* bound of the partial sum (3) is used: 1.1242 +* 1.1243 +* inf t[k] = sum a[j] l[j] + sum a[j] u[j], (13) 1.1244 +* j in Jp\{k} k in Jn\{k} 1.1245 +* 1.1246 +* where Jp = {j : a[j] > 0}, Jn = {j : a[j] < 0}, l[j] and u[j] are 1.1247 +* lower and upper bounds, resp., of variable x[j]. 1.1248 +* 1.1249 +* In order to compute inf t[k] more efficiently, the following formula, 1.1250 +* which is equivalent to (13), is actually used: 1.1251 +* 1.1252 +* ( h - a[k] l[k] = h, if a[k] > 0, 1.1253 +* inf t[k] = < (14) 1.1254 +* ( h - a[k] u[k] = h - a[k], if a[k] < 0, 1.1255 +* 1.1256 +* where: 1.1257 +* 1.1258 +* h = sum a[j] l[j] + sum a[j] u[j] (15) 1.1259 +* j in Jp j in Jn 1.1260 +* 1.1261 +* is the implied lower bound of row (1). 1.1262 +* 1.1263 +* Reduction of positive coefficient (a[k] > 0) does not change value 1.1264 +* of h, since l[k] = 0. In case of reduction of negative coefficient 1.1265 +* (a[k] < 0) from (11) it follows that: 1.1266 +* 1.1267 +* delta a[k] = a'[k] - a[k] = inf t[k] - b (> 0), (16) 1.1268 +* 1.1269 +* so new value of h (accounting that u[k] = 1) can be computed as 1.1270 +* follows: 1.1271 +* 1.1272 +* h := h + delta a[k] = h + (inf t[k] - b). (17) 1.1273 +* 1.1274 +* RECOVERING SOLUTION 1.1275 +* 1.1276 +* None needed. */ 1.1277 + 1.1278 +static int reduce_ineq_coef(NPP *npp, struct elem *ptr, double *_b) 1.1279 +{ /* process inequality constraint: sum a[j] x[j] >= b */ 1.1280 + /* returns: the number of coefficients reduced */ 1.1281 + struct elem *e; 1.1282 + int count = 0; 1.1283 + double h, inf_t, new_a, b = *_b; 1.1284 + xassert(npp == npp); 1.1285 + /* compute h; see (15) */ 1.1286 + h = 0.0; 1.1287 + for (e = ptr; e != NULL; e = e->next) 1.1288 + { if (e->aj > 0.0) 1.1289 + { if (e->xj->lb == -DBL_MAX) goto done; 1.1290 + h += e->aj * e->xj->lb; 1.1291 + } 1.1292 + else /* e->aj < 0.0 */ 1.1293 + { if (e->xj->ub == +DBL_MAX) goto done; 1.1294 + h += e->aj * e->xj->ub; 1.1295 + } 1.1296 + } 1.1297 + /* perform reduction of coefficients at binary variables */ 1.1298 + for (e = ptr; e != NULL; e = e->next) 1.1299 + { /* skip non-binary variable */ 1.1300 + if (!(e->xj->is_int && e->xj->lb == 0.0 && e->xj->ub == 1.0)) 1.1301 + continue; 1.1302 + if (e->aj > 0.0) 1.1303 + { /* compute inf t[k]; see (14) */ 1.1304 + inf_t = h; 1.1305 + if (b - e->aj < inf_t && inf_t < b) 1.1306 + { /* compute reduced coefficient a'[k]; see (7) */ 1.1307 + new_a = b - inf_t; 1.1308 + if (new_a >= +1e-3 && 1.1309 + e->aj - new_a >= 0.01 * (1.0 + e->aj)) 1.1310 + { /* accept a'[k] */ 1.1311 +#ifdef GLP_DEBUG 1.1312 + xprintf("+"); 1.1313 +#endif 1.1314 + e->aj = new_a; 1.1315 + count++; 1.1316 + } 1.1317 + } 1.1318 + } 1.1319 + else /* e->aj < 0.0 */ 1.1320 + { /* compute inf t[k]; see (14) */ 1.1321 + inf_t = h - e->aj; 1.1322 + if (b < inf_t && inf_t < b - e->aj) 1.1323 + { /* compute reduced coefficient a'[k]; see (11) */ 1.1324 + new_a = e->aj + (inf_t - b); 1.1325 + if (new_a <= -1e-3 && 1.1326 + new_a - e->aj >= 0.01 * (1.0 - e->aj)) 1.1327 + { /* accept a'[k] */ 1.1328 +#ifdef GLP_DEBUG 1.1329 + xprintf("-"); 1.1330 +#endif 1.1331 + e->aj = new_a; 1.1332 + /* update h; see (17) */ 1.1333 + h += (inf_t - b); 1.1334 + /* compute b'; see (9) */ 1.1335 + b = inf_t; 1.1336 + count++; 1.1337 + } 1.1338 + } 1.1339 + } 1.1340 + } 1.1341 + *_b = b; 1.1342 +done: return count; 1.1343 +} 1.1344 + 1.1345 +int npp_reduce_ineq_coef(NPP *npp, NPPROW *row) 1.1346 +{ /* reduce inequality constraint coefficients */ 1.1347 + NPPROW *copy; 1.1348 + NPPAIJ *aij; 1.1349 + struct elem *ptr, *e; 1.1350 + int kase, count[2]; 1.1351 + double b; 1.1352 + /* the row must be inequality constraint */ 1.1353 + xassert(row->lb < row->ub); 1.1354 + count[0] = count[1] = 0; 1.1355 + for (kase = 0; kase <= 1; kase++) 1.1356 + { if (kase == 0) 1.1357 + { /* process row lower bound */ 1.1358 + if (row->lb == -DBL_MAX) continue; 1.1359 +#ifdef GLP_DEBUG 1.1360 + xprintf("L"); 1.1361 +#endif 1.1362 + ptr = copy_form(npp, row, +1.0); 1.1363 + b = + row->lb; 1.1364 + } 1.1365 + else 1.1366 + { /* process row upper bound */ 1.1367 + if (row->ub == +DBL_MAX) continue; 1.1368 +#ifdef GLP_DEBUG 1.1369 + xprintf("U"); 1.1370 +#endif 1.1371 + ptr = copy_form(npp, row, -1.0); 1.1372 + b = - row->ub; 1.1373 + } 1.1374 + /* now the inequality has the form "sum a[j] x[j] >= b" */ 1.1375 + count[kase] = reduce_ineq_coef(npp, ptr, &b); 1.1376 + if (count[kase] > 0) 1.1377 + { /* the original inequality has been replaced by equivalent 1.1378 + one with coefficients reduced */ 1.1379 + if (row->lb == -DBL_MAX || row->ub == +DBL_MAX) 1.1380 + { /* the original row is single-sided inequality; no copy 1.1381 + is needed */ 1.1382 + copy = NULL; 1.1383 + } 1.1384 + else 1.1385 + { /* the original row is double-sided inequality; we need 1.1386 + to create its copy for other bound before replacing it 1.1387 + with the equivalent inequality */ 1.1388 +#ifdef GLP_DEBUG 1.1389 + xprintf("*"); 1.1390 +#endif 1.1391 + copy = npp_add_row(npp); 1.1392 + if (kase == 0) 1.1393 + { /* the copy is for upper bound */ 1.1394 + copy->lb = -DBL_MAX, copy->ub = row->ub; 1.1395 + } 1.1396 + else 1.1397 + { /* the copy is for lower bound */ 1.1398 + copy->lb = row->lb, copy->ub = +DBL_MAX; 1.1399 + } 1.1400 + /* copy original row coefficients */ 1.1401 + for (aij = row->ptr; aij != NULL; aij = aij->r_next) 1.1402 + npp_add_aij(npp, copy, aij->col, aij->val); 1.1403 + } 1.1404 + /* replace the original inequality by equivalent one */ 1.1405 + npp_erase_row(npp, row); 1.1406 + row->lb = b, row->ub = +DBL_MAX; 1.1407 + for (e = ptr; e != NULL; e = e->next) 1.1408 + npp_add_aij(npp, row, e->xj, e->aj); 1.1409 + /* continue processing upper bound for the copy */ 1.1410 + if (copy != NULL) row = copy; 1.1411 + } 1.1412 + drop_form(npp, ptr); 1.1413 + } 1.1414 + return count[0] + count[1]; 1.1415 +} 1.1416 + 1.1417 +/* eof */