1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpios06.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,1447 @@
1.4 +/* glpios06.c (MIR cut generator) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpios.h"
1.29 +
1.30 +#define _MIR_DEBUG 0
1.31 +
1.32 +#define MAXAGGR 5
1.33 +/* maximal number of rows which can be aggregated */
1.34 +
1.35 +struct MIR
1.36 +{ /* MIR cut generator working area */
1.37 + /*--------------------------------------------------------------*/
1.38 + /* global information valid for the root subproblem */
1.39 + int m;
1.40 + /* number of rows (in the root subproblem) */
1.41 + int n;
1.42 + /* number of columns */
1.43 + char *skip; /* char skip[1+m]; */
1.44 + /* skip[i], 1 <= i <= m, is a flag that means that row i should
1.45 + not be used because (1) it is not suitable, or (2) because it
1.46 + has been used in the aggregated constraint */
1.47 + char *isint; /* char isint[1+m+n]; */
1.48 + /* isint[k], 1 <= k <= m+n, is a flag that means that variable
1.49 + x[k] is integer (otherwise, continuous) */
1.50 + double *lb; /* double lb[1+m+n]; */
1.51 + /* lb[k], 1 <= k <= m+n, is lower bound of x[k]; -DBL_MAX means
1.52 + that x[k] has no lower bound */
1.53 + int *vlb; /* int vlb[1+m+n]; */
1.54 + /* vlb[k] = k', 1 <= k <= m+n, is the number of integer variable,
1.55 + which defines variable lower bound x[k] >= lb[k] * x[k']; zero
1.56 + means that x[k] has simple lower bound */
1.57 + double *ub; /* double ub[1+m+n]; */
1.58 + /* ub[k], 1 <= k <= m+n, is upper bound of x[k]; +DBL_MAX means
1.59 + that x[k] has no upper bound */
1.60 + int *vub; /* int vub[1+m+n]; */
1.61 + /* vub[k] = k', 1 <= k <= m+n, is the number of integer variable,
1.62 + which defines variable upper bound x[k] <= ub[k] * x[k']; zero
1.63 + means that x[k] has simple upper bound */
1.64 + /*--------------------------------------------------------------*/
1.65 + /* current (fractional) point to be separated */
1.66 + double *x; /* double x[1+m+n]; */
1.67 + /* x[k] is current value of auxiliary (1 <= k <= m) or structural
1.68 + (m+1 <= k <= m+n) variable */
1.69 + /*--------------------------------------------------------------*/
1.70 + /* aggregated constraint sum a[k] * x[k] = b, which is a linear
1.71 + combination of original constraints transformed to equalities
1.72 + by introducing auxiliary variables */
1.73 + int agg_cnt;
1.74 + /* number of rows (original constraints) used to build aggregated
1.75 + constraint, 1 <= agg_cnt <= MAXAGGR */
1.76 + int *agg_row; /* int agg_row[1+MAXAGGR]; */
1.77 + /* agg_row[k], 1 <= k <= agg_cnt, is the row number used to build
1.78 + aggregated constraint */
1.79 + IOSVEC *agg_vec; /* IOSVEC agg_vec[1:m+n]; */
1.80 + /* sparse vector of aggregated constraint coefficients, a[k] */
1.81 + double agg_rhs;
1.82 + /* right-hand side of the aggregated constraint, b */
1.83 + /*--------------------------------------------------------------*/
1.84 + /* bound substitution flags for modified constraint */
1.85 + char *subst; /* char subst[1+m+n]; */
1.86 + /* subst[k], 1 <= k <= m+n, is a bound substitution flag used for
1.87 + variable x[k]:
1.88 + '?' - x[k] is missing in modified constraint
1.89 + 'L' - x[k] = (lower bound) + x'[k]
1.90 + 'U' - x[k] = (upper bound) - x'[k] */
1.91 + /*--------------------------------------------------------------*/
1.92 + /* modified constraint sum a'[k] * x'[k] = b', where x'[k] >= 0,
1.93 + derived from aggregated constraint by substituting bounds;
1.94 + note that due to substitution of variable bounds there may be
1.95 + additional terms in the modified constraint */
1.96 + IOSVEC *mod_vec; /* IOSVEC mod_vec[1:m+n]; */
1.97 + /* sparse vector of modified constraint coefficients, a'[k] */
1.98 + double mod_rhs;
1.99 + /* right-hand side of the modified constraint, b' */
1.100 + /*--------------------------------------------------------------*/
1.101 + /* cutting plane sum alpha[k] * x[k] <= beta */
1.102 + IOSVEC *cut_vec; /* IOSVEC cut_vec[1:m+n]; */
1.103 + /* sparse vector of cutting plane coefficients, alpha[k] */
1.104 + double cut_rhs;
1.105 + /* right-hand size of the cutting plane, beta */
1.106 +};
1.107 +
1.108 +/***********************************************************************
1.109 +* NAME
1.110 +*
1.111 +* ios_mir_init - initialize MIR cut generator
1.112 +*
1.113 +* SYNOPSIS
1.114 +*
1.115 +* #include "glpios.h"
1.116 +* void *ios_mir_init(glp_tree *tree);
1.117 +*
1.118 +* DESCRIPTION
1.119 +*
1.120 +* The routine ios_mir_init initializes the MIR cut generator assuming
1.121 +* that the current subproblem is the root subproblem.
1.122 +*
1.123 +* RETURNS
1.124 +*
1.125 +* The routine ios_mir_init returns a pointer to the MIR cut generator
1.126 +* working area. */
1.127 +
1.128 +static void set_row_attrib(glp_tree *tree, struct MIR *mir)
1.129 +{ /* set global row attributes */
1.130 + glp_prob *mip = tree->mip;
1.131 + int m = mir->m;
1.132 + int k;
1.133 + for (k = 1; k <= m; k++)
1.134 + { GLPROW *row = mip->row[k];
1.135 + mir->skip[k] = 0;
1.136 + mir->isint[k] = 0;
1.137 + switch (row->type)
1.138 + { case GLP_FR:
1.139 + mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
1.140 + case GLP_LO:
1.141 + mir->lb[k] = row->lb, mir->ub[k] = +DBL_MAX; break;
1.142 + case GLP_UP:
1.143 + mir->lb[k] = -DBL_MAX, mir->ub[k] = row->ub; break;
1.144 + case GLP_DB:
1.145 + mir->lb[k] = row->lb, mir->ub[k] = row->ub; break;
1.146 + case GLP_FX:
1.147 + mir->lb[k] = mir->ub[k] = row->lb; break;
1.148 + default:
1.149 + xassert(row != row);
1.150 + }
1.151 + mir->vlb[k] = mir->vub[k] = 0;
1.152 + }
1.153 + return;
1.154 +}
1.155 +
1.156 +static void set_col_attrib(glp_tree *tree, struct MIR *mir)
1.157 +{ /* set global column attributes */
1.158 + glp_prob *mip = tree->mip;
1.159 + int m = mir->m;
1.160 + int n = mir->n;
1.161 + int k;
1.162 + for (k = m+1; k <= m+n; k++)
1.163 + { GLPCOL *col = mip->col[k-m];
1.164 + switch (col->kind)
1.165 + { case GLP_CV:
1.166 + mir->isint[k] = 0; break;
1.167 + case GLP_IV:
1.168 + mir->isint[k] = 1; break;
1.169 + default:
1.170 + xassert(col != col);
1.171 + }
1.172 + switch (col->type)
1.173 + { case GLP_FR:
1.174 + mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
1.175 + case GLP_LO:
1.176 + mir->lb[k] = col->lb, mir->ub[k] = +DBL_MAX; break;
1.177 + case GLP_UP:
1.178 + mir->lb[k] = -DBL_MAX, mir->ub[k] = col->ub; break;
1.179 + case GLP_DB:
1.180 + mir->lb[k] = col->lb, mir->ub[k] = col->ub; break;
1.181 + case GLP_FX:
1.182 + mir->lb[k] = mir->ub[k] = col->lb; break;
1.183 + default:
1.184 + xassert(col != col);
1.185 + }
1.186 + mir->vlb[k] = mir->vub[k] = 0;
1.187 + }
1.188 + return;
1.189 +}
1.190 +
1.191 +static void set_var_bounds(glp_tree *tree, struct MIR *mir)
1.192 +{ /* set variable bounds */
1.193 + glp_prob *mip = tree->mip;
1.194 + int m = mir->m;
1.195 + GLPAIJ *aij;
1.196 + int i, k1, k2;
1.197 + double a1, a2;
1.198 + for (i = 1; i <= m; i++)
1.199 + { /* we need the row to be '>= 0' or '<= 0' */
1.200 + if (!(mir->lb[i] == 0.0 && mir->ub[i] == +DBL_MAX ||
1.201 + mir->lb[i] == -DBL_MAX && mir->ub[i] == 0.0)) continue;
1.202 + /* take first term */
1.203 + aij = mip->row[i]->ptr;
1.204 + if (aij == NULL) continue;
1.205 + k1 = m + aij->col->j, a1 = aij->val;
1.206 + /* take second term */
1.207 + aij = aij->r_next;
1.208 + if (aij == NULL) continue;
1.209 + k2 = m + aij->col->j, a2 = aij->val;
1.210 + /* there must be only two terms */
1.211 + if (aij->r_next != NULL) continue;
1.212 + /* interchange terms, if needed */
1.213 + if (!mir->isint[k1] && mir->isint[k2])
1.214 + ;
1.215 + else if (mir->isint[k1] && !mir->isint[k2])
1.216 + { k2 = k1, a2 = a1;
1.217 + k1 = m + aij->col->j, a1 = aij->val;
1.218 + }
1.219 + else
1.220 + { /* both terms are either continuous or integer */
1.221 + continue;
1.222 + }
1.223 + /* x[k2] should be double-bounded */
1.224 + if (mir->lb[k2] == -DBL_MAX || mir->ub[k2] == +DBL_MAX ||
1.225 + mir->lb[k2] == mir->ub[k2]) continue;
1.226 + /* change signs, if necessary */
1.227 + if (mir->ub[i] == 0.0) a1 = - a1, a2 = - a2;
1.228 + /* now the row has the form a1 * x1 + a2 * x2 >= 0, where x1
1.229 + is continuous, x2 is integer */
1.230 + if (a1 > 0.0)
1.231 + { /* x1 >= - (a2 / a1) * x2 */
1.232 + if (mir->vlb[k1] == 0)
1.233 + { /* set variable lower bound for x1 */
1.234 + mir->lb[k1] = - a2 / a1;
1.235 + mir->vlb[k1] = k2;
1.236 + /* the row should not be used */
1.237 + mir->skip[i] = 1;
1.238 + }
1.239 + }
1.240 + else /* a1 < 0.0 */
1.241 + { /* x1 <= - (a2 / a1) * x2 */
1.242 + if (mir->vub[k1] == 0)
1.243 + { /* set variable upper bound for x1 */
1.244 + mir->ub[k1] = - a2 / a1;
1.245 + mir->vub[k1] = k2;
1.246 + /* the row should not be used */
1.247 + mir->skip[i] = 1;
1.248 + }
1.249 + }
1.250 + }
1.251 + return;
1.252 +}
1.253 +
1.254 +static void mark_useless_rows(glp_tree *tree, struct MIR *mir)
1.255 +{ /* mark rows which should not be used */
1.256 + glp_prob *mip = tree->mip;
1.257 + int m = mir->m;
1.258 + GLPAIJ *aij;
1.259 + int i, k, nv;
1.260 + for (i = 1; i <= m; i++)
1.261 + { /* free rows should not be used */
1.262 + if (mir->lb[i] == -DBL_MAX && mir->ub[i] == +DBL_MAX)
1.263 + { mir->skip[i] = 1;
1.264 + continue;
1.265 + }
1.266 + nv = 0;
1.267 + for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
1.268 + { k = m + aij->col->j;
1.269 + /* rows with free variables should not be used */
1.270 + if (mir->lb[k] == -DBL_MAX && mir->ub[k] == +DBL_MAX)
1.271 + { mir->skip[i] = 1;
1.272 + break;
1.273 + }
1.274 + /* rows with integer variables having infinite (lower or
1.275 + upper) bound should not be used */
1.276 + if (mir->isint[k] && mir->lb[k] == -DBL_MAX ||
1.277 + mir->isint[k] && mir->ub[k] == +DBL_MAX)
1.278 + { mir->skip[i] = 1;
1.279 + break;
1.280 + }
1.281 + /* count non-fixed variables */
1.282 + if (!(mir->vlb[k] == 0 && mir->vub[k] == 0 &&
1.283 + mir->lb[k] == mir->ub[k])) nv++;
1.284 + }
1.285 + /* rows with all variables fixed should not be used */
1.286 + if (nv == 0)
1.287 + { mir->skip[i] = 1;
1.288 + continue;
1.289 + }
1.290 + }
1.291 + return;
1.292 +}
1.293 +
1.294 +void *ios_mir_init(glp_tree *tree)
1.295 +{ /* initialize MIR cut generator */
1.296 + glp_prob *mip = tree->mip;
1.297 + int m = mip->m;
1.298 + int n = mip->n;
1.299 + struct MIR *mir;
1.300 +#if _MIR_DEBUG
1.301 + xprintf("ios_mir_init: warning: debug mode enabled\n");
1.302 +#endif
1.303 + /* allocate working area */
1.304 + mir = xmalloc(sizeof(struct MIR));
1.305 + mir->m = m;
1.306 + mir->n = n;
1.307 + mir->skip = xcalloc(1+m, sizeof(char));
1.308 + mir->isint = xcalloc(1+m+n, sizeof(char));
1.309 + mir->lb = xcalloc(1+m+n, sizeof(double));
1.310 + mir->vlb = xcalloc(1+m+n, sizeof(int));
1.311 + mir->ub = xcalloc(1+m+n, sizeof(double));
1.312 + mir->vub = xcalloc(1+m+n, sizeof(int));
1.313 + mir->x = xcalloc(1+m+n, sizeof(double));
1.314 + mir->agg_row = xcalloc(1+MAXAGGR, sizeof(int));
1.315 + mir->agg_vec = ios_create_vec(m+n);
1.316 + mir->subst = xcalloc(1+m+n, sizeof(char));
1.317 + mir->mod_vec = ios_create_vec(m+n);
1.318 + mir->cut_vec = ios_create_vec(m+n);
1.319 + /* set global row attributes */
1.320 + set_row_attrib(tree, mir);
1.321 + /* set global column attributes */
1.322 + set_col_attrib(tree, mir);
1.323 + /* set variable bounds */
1.324 + set_var_bounds(tree, mir);
1.325 + /* mark rows which should not be used */
1.326 + mark_useless_rows(tree, mir);
1.327 + return mir;
1.328 +}
1.329 +
1.330 +/***********************************************************************
1.331 +* NAME
1.332 +*
1.333 +* ios_mir_gen - generate MIR cuts
1.334 +*
1.335 +* SYNOPSIS
1.336 +*
1.337 +* #include "glpios.h"
1.338 +* void ios_mir_gen(glp_tree *tree, void *gen, IOSPOOL *pool);
1.339 +*
1.340 +* DESCRIPTION
1.341 +*
1.342 +* The routine ios_mir_gen generates MIR cuts for the current point and
1.343 +* adds them to the cut pool. */
1.344 +
1.345 +static void get_current_point(glp_tree *tree, struct MIR *mir)
1.346 +{ /* obtain current point */
1.347 + glp_prob *mip = tree->mip;
1.348 + int m = mir->m;
1.349 + int n = mir->n;
1.350 + int k;
1.351 + for (k = 1; k <= m; k++)
1.352 + mir->x[k] = mip->row[k]->prim;
1.353 + for (k = m+1; k <= m+n; k++)
1.354 + mir->x[k] = mip->col[k-m]->prim;
1.355 + return;
1.356 +}
1.357 +
1.358 +#if _MIR_DEBUG
1.359 +static void check_current_point(struct MIR *mir)
1.360 +{ /* check current point */
1.361 + int m = mir->m;
1.362 + int n = mir->n;
1.363 + int k, kk;
1.364 + double lb, ub, eps;
1.365 + for (k = 1; k <= m+n; k++)
1.366 + { /* determine lower bound */
1.367 + lb = mir->lb[k];
1.368 + kk = mir->vlb[k];
1.369 + if (kk != 0)
1.370 + { xassert(lb != -DBL_MAX);
1.371 + xassert(!mir->isint[k]);
1.372 + xassert(mir->isint[kk]);
1.373 + lb *= mir->x[kk];
1.374 + }
1.375 + /* check lower bound */
1.376 + if (lb != -DBL_MAX)
1.377 + { eps = 1e-6 * (1.0 + fabs(lb));
1.378 + xassert(mir->x[k] >= lb - eps);
1.379 + }
1.380 + /* determine upper bound */
1.381 + ub = mir->ub[k];
1.382 + kk = mir->vub[k];
1.383 + if (kk != 0)
1.384 + { xassert(ub != +DBL_MAX);
1.385 + xassert(!mir->isint[k]);
1.386 + xassert(mir->isint[kk]);
1.387 + ub *= mir->x[kk];
1.388 + }
1.389 + /* check upper bound */
1.390 + if (ub != +DBL_MAX)
1.391 + { eps = 1e-6 * (1.0 + fabs(ub));
1.392 + xassert(mir->x[k] <= ub + eps);
1.393 + }
1.394 + }
1.395 + return;
1.396 +}
1.397 +#endif
1.398 +
1.399 +static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i)
1.400 +{ /* use original i-th row as initial aggregated constraint */
1.401 + glp_prob *mip = tree->mip;
1.402 + int m = mir->m;
1.403 + GLPAIJ *aij;
1.404 + xassert(1 <= i && i <= m);
1.405 + xassert(!mir->skip[i]);
1.406 + /* mark i-th row in order not to use it in the same aggregated
1.407 + constraint */
1.408 + mir->skip[i] = 2;
1.409 + mir->agg_cnt = 1;
1.410 + mir->agg_row[1] = i;
1.411 + /* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary
1.412 + variable of row i, x[m+j] are structural variables */
1.413 + ios_clear_vec(mir->agg_vec);
1.414 + ios_set_vj(mir->agg_vec, i, 1.0);
1.415 + for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
1.416 + ios_set_vj(mir->agg_vec, m + aij->col->j, - aij->val);
1.417 + mir->agg_rhs = 0.0;
1.418 +#if _MIR_DEBUG
1.419 + ios_check_vec(mir->agg_vec);
1.420 +#endif
1.421 + return;
1.422 +}
1.423 +
1.424 +#if _MIR_DEBUG
1.425 +static void check_agg_row(struct MIR *mir)
1.426 +{ /* check aggregated constraint */
1.427 + int m = mir->m;
1.428 + int n = mir->n;
1.429 + int j, k;
1.430 + double r, big;
1.431 + /* compute the residual r = sum a[k] * x[k] - b and determine
1.432 + big = max(1, |a[k]|, |b|) */
1.433 + r = 0.0, big = 1.0;
1.434 + for (j = 1; j <= mir->agg_vec->nnz; j++)
1.435 + { k = mir->agg_vec->ind[j];
1.436 + xassert(1 <= k && k <= m+n);
1.437 + r += mir->agg_vec->val[j] * mir->x[k];
1.438 + if (big < fabs(mir->agg_vec->val[j]))
1.439 + big = fabs(mir->agg_vec->val[j]);
1.440 + }
1.441 + r -= mir->agg_rhs;
1.442 + if (big < fabs(mir->agg_rhs))
1.443 + big = fabs(mir->agg_rhs);
1.444 + /* the residual must be close to zero */
1.445 + xassert(fabs(r) <= 1e-6 * big);
1.446 + return;
1.447 +}
1.448 +#endif
1.449 +
1.450 +static void subst_fixed_vars(struct MIR *mir)
1.451 +{ /* substitute fixed variables into aggregated constraint */
1.452 + int m = mir->m;
1.453 + int n = mir->n;
1.454 + int j, k;
1.455 + for (j = 1; j <= mir->agg_vec->nnz; j++)
1.456 + { k = mir->agg_vec->ind[j];
1.457 + xassert(1 <= k && k <= m+n);
1.458 + if (mir->vlb[k] == 0 && mir->vub[k] == 0 &&
1.459 + mir->lb[k] == mir->ub[k])
1.460 + { /* x[k] is fixed */
1.461 + mir->agg_rhs -= mir->agg_vec->val[j] * mir->lb[k];
1.462 + mir->agg_vec->val[j] = 0.0;
1.463 + }
1.464 + }
1.465 + /* remove terms corresponding to fixed variables */
1.466 + ios_clean_vec(mir->agg_vec, DBL_EPSILON);
1.467 +#if _MIR_DEBUG
1.468 + ios_check_vec(mir->agg_vec);
1.469 +#endif
1.470 + return;
1.471 +}
1.472 +
1.473 +static void bound_subst_heur(struct MIR *mir)
1.474 +{ /* bound substitution heuristic */
1.475 + int m = mir->m;
1.476 + int n = mir->n;
1.477 + int j, k, kk;
1.478 + double d1, d2;
1.479 + for (j = 1; j <= mir->agg_vec->nnz; j++)
1.480 + { k = mir->agg_vec->ind[j];
1.481 + xassert(1 <= k && k <= m+n);
1.482 + if (mir->isint[k]) continue; /* skip integer variable */
1.483 + /* compute distance from x[k] to its lower bound */
1.484 + kk = mir->vlb[k];
1.485 + if (kk == 0)
1.486 + { if (mir->lb[k] == -DBL_MAX)
1.487 + d1 = DBL_MAX;
1.488 + else
1.489 + d1 = mir->x[k] - mir->lb[k];
1.490 + }
1.491 + else
1.492 + { xassert(1 <= kk && kk <= m+n);
1.493 + xassert(mir->isint[kk]);
1.494 + xassert(mir->lb[k] != -DBL_MAX);
1.495 + d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
1.496 + }
1.497 + /* compute distance from x[k] to its upper bound */
1.498 + kk = mir->vub[k];
1.499 + if (kk == 0)
1.500 + { if (mir->vub[k] == +DBL_MAX)
1.501 + d2 = DBL_MAX;
1.502 + else
1.503 + d2 = mir->ub[k] - mir->x[k];
1.504 + }
1.505 + else
1.506 + { xassert(1 <= kk && kk <= m+n);
1.507 + xassert(mir->isint[kk]);
1.508 + xassert(mir->ub[k] != +DBL_MAX);
1.509 + d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
1.510 + }
1.511 + /* x[k] cannot be free */
1.512 + xassert(d1 != DBL_MAX || d2 != DBL_MAX);
1.513 + /* choose the bound which is closer to x[k] */
1.514 + xassert(mir->subst[k] == '?');
1.515 + if (d1 <= d2)
1.516 + mir->subst[k] = 'L';
1.517 + else
1.518 + mir->subst[k] = 'U';
1.519 + }
1.520 + return;
1.521 +}
1.522 +
1.523 +static void build_mod_row(struct MIR *mir)
1.524 +{ /* substitute bounds and build modified constraint */
1.525 + int m = mir->m;
1.526 + int n = mir->n;
1.527 + int j, jj, k, kk;
1.528 + /* initially modified constraint is aggregated constraint */
1.529 + ios_copy_vec(mir->mod_vec, mir->agg_vec);
1.530 + mir->mod_rhs = mir->agg_rhs;
1.531 +#if _MIR_DEBUG
1.532 + ios_check_vec(mir->mod_vec);
1.533 +#endif
1.534 + /* substitute bounds for continuous variables; note that due to
1.535 + substitution of variable bounds additional terms may appear in
1.536 + modified constraint */
1.537 + for (j = mir->mod_vec->nnz; j >= 1; j--)
1.538 + { k = mir->mod_vec->ind[j];
1.539 + xassert(1 <= k && k <= m+n);
1.540 + if (mir->isint[k]) continue; /* skip integer variable */
1.541 + if (mir->subst[k] == 'L')
1.542 + { /* x[k] = (lower bound) + x'[k] */
1.543 + xassert(mir->lb[k] != -DBL_MAX);
1.544 + kk = mir->vlb[k];
1.545 + if (kk == 0)
1.546 + { /* x[k] = lb[k] + x'[k] */
1.547 + mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
1.548 + }
1.549 + else
1.550 + { /* x[k] = lb[k] * x[kk] + x'[k] */
1.551 + xassert(mir->isint[kk]);
1.552 + jj = mir->mod_vec->pos[kk];
1.553 + if (jj == 0)
1.554 + { ios_set_vj(mir->mod_vec, kk, 1.0);
1.555 + jj = mir->mod_vec->pos[kk];
1.556 + mir->mod_vec->val[jj] = 0.0;
1.557 + }
1.558 + mir->mod_vec->val[jj] +=
1.559 + mir->mod_vec->val[j] * mir->lb[k];
1.560 + }
1.561 + }
1.562 + else if (mir->subst[k] == 'U')
1.563 + { /* x[k] = (upper bound) - x'[k] */
1.564 + xassert(mir->ub[k] != +DBL_MAX);
1.565 + kk = mir->vub[k];
1.566 + if (kk == 0)
1.567 + { /* x[k] = ub[k] - x'[k] */
1.568 + mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
1.569 + }
1.570 + else
1.571 + { /* x[k] = ub[k] * x[kk] - x'[k] */
1.572 + xassert(mir->isint[kk]);
1.573 + jj = mir->mod_vec->pos[kk];
1.574 + if (jj == 0)
1.575 + { ios_set_vj(mir->mod_vec, kk, 1.0);
1.576 + jj = mir->mod_vec->pos[kk];
1.577 + mir->mod_vec->val[jj] = 0.0;
1.578 + }
1.579 + mir->mod_vec->val[jj] +=
1.580 + mir->mod_vec->val[j] * mir->ub[k];
1.581 + }
1.582 + mir->mod_vec->val[j] = - mir->mod_vec->val[j];
1.583 + }
1.584 + else
1.585 + xassert(k != k);
1.586 + }
1.587 +#if _MIR_DEBUG
1.588 + ios_check_vec(mir->mod_vec);
1.589 +#endif
1.590 + /* substitute bounds for integer variables */
1.591 + for (j = 1; j <= mir->mod_vec->nnz; j++)
1.592 + { k = mir->mod_vec->ind[j];
1.593 + xassert(1 <= k && k <= m+n);
1.594 + if (!mir->isint[k]) continue; /* skip continuous variable */
1.595 + xassert(mir->subst[k] == '?');
1.596 + xassert(mir->vlb[k] == 0 && mir->vub[k] == 0);
1.597 + xassert(mir->lb[k] != -DBL_MAX && mir->ub[k] != +DBL_MAX);
1.598 + if (fabs(mir->lb[k]) <= fabs(mir->ub[k]))
1.599 + { /* x[k] = lb[k] + x'[k] */
1.600 + mir->subst[k] = 'L';
1.601 + mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
1.602 + }
1.603 + else
1.604 + { /* x[k] = ub[k] - x'[k] */
1.605 + mir->subst[k] = 'U';
1.606 + mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
1.607 + mir->mod_vec->val[j] = - mir->mod_vec->val[j];
1.608 + }
1.609 + }
1.610 +#if _MIR_DEBUG
1.611 + ios_check_vec(mir->mod_vec);
1.612 +#endif
1.613 + return;
1.614 +}
1.615 +
1.616 +#if _MIR_DEBUG
1.617 +static void check_mod_row(struct MIR *mir)
1.618 +{ /* check modified constraint */
1.619 + int m = mir->m;
1.620 + int n = mir->n;
1.621 + int j, k, kk;
1.622 + double r, big, x;
1.623 + /* compute the residual r = sum a'[k] * x'[k] - b' and determine
1.624 + big = max(1, |a[k]|, |b|) */
1.625 + r = 0.0, big = 1.0;
1.626 + for (j = 1; j <= mir->mod_vec->nnz; j++)
1.627 + { k = mir->mod_vec->ind[j];
1.628 + xassert(1 <= k && k <= m+n);
1.629 + if (mir->subst[k] == 'L')
1.630 + { /* x'[k] = x[k] - (lower bound) */
1.631 + xassert(mir->lb[k] != -DBL_MAX);
1.632 + kk = mir->vlb[k];
1.633 + if (kk == 0)
1.634 + x = mir->x[k] - mir->lb[k];
1.635 + else
1.636 + x = mir->x[k] - mir->lb[k] * mir->x[kk];
1.637 + }
1.638 + else if (mir->subst[k] == 'U')
1.639 + { /* x'[k] = (upper bound) - x[k] */
1.640 + xassert(mir->ub[k] != +DBL_MAX);
1.641 + kk = mir->vub[k];
1.642 + if (kk == 0)
1.643 + x = mir->ub[k] - mir->x[k];
1.644 + else
1.645 + x = mir->ub[k] * mir->x[kk] - mir->x[k];
1.646 + }
1.647 + else
1.648 + xassert(k != k);
1.649 + r += mir->mod_vec->val[j] * x;
1.650 + if (big < fabs(mir->mod_vec->val[j]))
1.651 + big = fabs(mir->mod_vec->val[j]);
1.652 + }
1.653 + r -= mir->mod_rhs;
1.654 + if (big < fabs(mir->mod_rhs))
1.655 + big = fabs(mir->mod_rhs);
1.656 + /* the residual must be close to zero */
1.657 + xassert(fabs(r) <= 1e-6 * big);
1.658 + return;
1.659 +}
1.660 +#endif
1.661 +
1.662 +/***********************************************************************
1.663 +* mir_ineq - construct MIR inequality
1.664 +*
1.665 +* Given the single constraint mixed integer set
1.666 +*
1.667 +* |N|
1.668 +* X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s},
1.669 +* + + j in N
1.670 +*
1.671 +* this routine constructs the mixed integer rounding (MIR) inequality
1.672 +*
1.673 +* sum alpha[j] * x[j] <= beta + gamma * s,
1.674 +* j in N
1.675 +*
1.676 +* which is valid for X.
1.677 +*
1.678 +* If the MIR inequality has been successfully constructed, the routine
1.679 +* returns zero. Otherwise, if b is close to nearest integer, there may
1.680 +* be numeric difficulties due to big coefficients; so in this case the
1.681 +* routine returns non-zero. */
1.682 +
1.683 +static int mir_ineq(const int n, const double a[], const double b,
1.684 + double alpha[], double *beta, double *gamma)
1.685 +{ int j;
1.686 + double f, t;
1.687 + if (fabs(b - floor(b + .5)) < 0.01)
1.688 + return 1;
1.689 + f = b - floor(b);
1.690 + for (j = 1; j <= n; j++)
1.691 + { t = (a[j] - floor(a[j])) - f;
1.692 + if (t <= 0.0)
1.693 + alpha[j] = floor(a[j]);
1.694 + else
1.695 + alpha[j] = floor(a[j]) + t / (1.0 - f);
1.696 + }
1.697 + *beta = floor(b);
1.698 + *gamma = 1.0 / (1.0 - f);
1.699 + return 0;
1.700 +}
1.701 +
1.702 +/***********************************************************************
1.703 +* cmir_ineq - construct c-MIR inequality
1.704 +*
1.705 +* Given the mixed knapsack set
1.706 +*
1.707 +* MK |N|
1.708 +* X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s,
1.709 +* + + j in N
1.710 +*
1.711 +* x[j] <= u[j]},
1.712 +*
1.713 +* a subset C of variables to be complemented, and a divisor delta > 0,
1.714 +* this routine constructs the complemented MIR (c-MIR) inequality
1.715 +*
1.716 +* sum alpha[j] * x[j] <= beta + gamma * s,
1.717 +* j in N
1.718 +* MK
1.719 +* which is valid for X .
1.720 +*
1.721 +* If the c-MIR inequality has been successfully constructed, the
1.722 +* routine returns zero. Otherwise, if there is a risk of numerical
1.723 +* difficulties due to big coefficients (see comments to the routine
1.724 +* mir_ineq), the routine cmir_ineq returns non-zero. */
1.725 +
1.726 +static int cmir_ineq(const int n, const double a[], const double b,
1.727 + const double u[], const char cset[], const double delta,
1.728 + double alpha[], double *beta, double *gamma)
1.729 +{ int j;
1.730 + double *aa, bb;
1.731 + aa = alpha, bb = b;
1.732 + for (j = 1; j <= n; j++)
1.733 + { aa[j] = a[j] / delta;
1.734 + if (cset[j])
1.735 + aa[j] = - aa[j], bb -= a[j] * u[j];
1.736 + }
1.737 + bb /= delta;
1.738 + if (mir_ineq(n, aa, bb, alpha, beta, gamma)) return 1;
1.739 + for (j = 1; j <= n; j++)
1.740 + { if (cset[j])
1.741 + alpha[j] = - alpha[j], *beta += alpha[j] * u[j];
1.742 + }
1.743 + *gamma /= delta;
1.744 + return 0;
1.745 +}
1.746 +
1.747 +/***********************************************************************
1.748 +* cmir_sep - c-MIR separation heuristic
1.749 +*
1.750 +* Given the mixed knapsack set
1.751 +*
1.752 +* MK |N|
1.753 +* X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s,
1.754 +* + + j in N
1.755 +*
1.756 +* x[j] <= u[j]}
1.757 +*
1.758 +* * *
1.759 +* and a fractional point (x , s ), this routine tries to construct
1.760 +* c-MIR inequality
1.761 +*
1.762 +* sum alpha[j] * x[j] <= beta + gamma * s,
1.763 +* j in N
1.764 +* MK
1.765 +* which is valid for X and has (desirably maximal) violation at the
1.766 +* fractional point given. This is attained by choosing an appropriate
1.767 +* set C of variables to be complemented and a divisor delta > 0, which
1.768 +* together define corresponding c-MIR inequality.
1.769 +*
1.770 +* If a violated c-MIR inequality has been successfully constructed,
1.771 +* the routine returns its violation:
1.772 +*
1.773 +* * *
1.774 +* sum alpha[j] * x [j] - beta - gamma * s ,
1.775 +* j in N
1.776 +*
1.777 +* which is positive. In case of failure the routine returns zero. */
1.778 +
1.779 +struct vset { int j; double v; };
1.780 +
1.781 +static int cmir_cmp(const void *p1, const void *p2)
1.782 +{ const struct vset *v1 = p1, *v2 = p2;
1.783 + if (v1->v < v2->v) return -1;
1.784 + if (v1->v > v2->v) return +1;
1.785 + return 0;
1.786 +}
1.787 +
1.788 +static double cmir_sep(const int n, const double a[], const double b,
1.789 + const double u[], const double x[], const double s,
1.790 + double alpha[], double *beta, double *gamma)
1.791 +{ int fail, j, k, nv, v;
1.792 + double delta, eps, d_try[1+3], r, r_best;
1.793 + char *cset;
1.794 + struct vset *vset;
1.795 + /* allocate working arrays */
1.796 + cset = xcalloc(1+n, sizeof(char));
1.797 + vset = xcalloc(1+n, sizeof(struct vset));
1.798 + /* choose initial C */
1.799 + for (j = 1; j <= n; j++)
1.800 + cset[j] = (char)(x[j] >= 0.5 * u[j]);
1.801 + /* choose initial delta */
1.802 + r_best = delta = 0.0;
1.803 + for (j = 1; j <= n; j++)
1.804 + { xassert(a[j] != 0.0);
1.805 + /* if x[j] is close to its bounds, skip it */
1.806 + eps = 1e-9 * (1.0 + fabs(u[j]));
1.807 + if (x[j] < eps || x[j] > u[j] - eps) continue;
1.808 + /* try delta = |a[j]| to construct c-MIR inequality */
1.809 + fail = cmir_ineq(n, a, b, u, cset, fabs(a[j]), alpha, beta,
1.810 + gamma);
1.811 + if (fail) continue;
1.812 + /* compute violation */
1.813 + r = - (*beta) - (*gamma) * s;
1.814 + for (k = 1; k <= n; k++) r += alpha[k] * x[k];
1.815 + if (r_best < r) r_best = r, delta = fabs(a[j]);
1.816 + }
1.817 + if (r_best < 0.001) r_best = 0.0;
1.818 + if (r_best == 0.0) goto done;
1.819 + xassert(delta > 0.0);
1.820 + /* try to increase violation by dividing delta by 2, 4, and 8,
1.821 + respectively */
1.822 + d_try[1] = delta / 2.0;
1.823 + d_try[2] = delta / 4.0;
1.824 + d_try[3] = delta / 8.0;
1.825 + for (j = 1; j <= 3; j++)
1.826 + { /* construct c-MIR inequality */
1.827 + fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,
1.828 + gamma);
1.829 + if (fail) continue;
1.830 + /* compute violation */
1.831 + r = - (*beta) - (*gamma) * s;
1.832 + for (k = 1; k <= n; k++) r += alpha[k] * x[k];
1.833 + if (r_best < r) r_best = r, delta = d_try[j];
1.834 + }
1.835 + /* build subset of variables lying strictly between their bounds
1.836 + and order it by nondecreasing values of |x[j] - u[j]/2| */
1.837 + nv = 0;
1.838 + for (j = 1; j <= n; j++)
1.839 + { /* if x[j] is close to its bounds, skip it */
1.840 + eps = 1e-9 * (1.0 + fabs(u[j]));
1.841 + if (x[j] < eps || x[j] > u[j] - eps) continue;
1.842 + /* add x[j] to the subset */
1.843 + nv++;
1.844 + vset[nv].j = j;
1.845 + vset[nv].v = fabs(x[j] - 0.5 * u[j]);
1.846 + }
1.847 + qsort(&vset[1], nv, sizeof(struct vset), cmir_cmp);
1.848 + /* try to increase violation by successively complementing each
1.849 + variable in the subset */
1.850 + for (v = 1; v <= nv; v++)
1.851 + { j = vset[v].j;
1.852 + /* replace x[j] by its complement or vice versa */
1.853 + cset[j] = (char)!cset[j];
1.854 + /* construct c-MIR inequality */
1.855 + fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
1.856 + /* restore the variable */
1.857 + cset[j] = (char)!cset[j];
1.858 + /* do not replace the variable in case of failure */
1.859 + if (fail) continue;
1.860 + /* compute violation */
1.861 + r = - (*beta) - (*gamma) * s;
1.862 + for (k = 1; k <= n; k++) r += alpha[k] * x[k];
1.863 + if (r_best < r) r_best = r, cset[j] = (char)!cset[j];
1.864 + }
1.865 + /* construct the best c-MIR inequality chosen */
1.866 + fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
1.867 + xassert(!fail);
1.868 +done: /* free working arrays */
1.869 + xfree(cset);
1.870 + xfree(vset);
1.871 + /* return to the calling routine */
1.872 + return r_best;
1.873 +}
1.874 +
1.875 +static double generate(struct MIR *mir)
1.876 +{ /* try to generate violated c-MIR cut for modified constraint */
1.877 + int m = mir->m;
1.878 + int n = mir->n;
1.879 + int j, k, kk, nint;
1.880 + double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma;
1.881 + ios_copy_vec(mir->cut_vec, mir->mod_vec);
1.882 + mir->cut_rhs = mir->mod_rhs;
1.883 + /* remove small terms, which can appear due to substitution of
1.884 + variable bounds */
1.885 + ios_clean_vec(mir->cut_vec, DBL_EPSILON);
1.886 +#if _MIR_DEBUG
1.887 + ios_check_vec(mir->cut_vec);
1.888 +#endif
1.889 + /* remove positive continuous terms to obtain MK relaxation */
1.890 + for (j = 1; j <= mir->cut_vec->nnz; j++)
1.891 + { k = mir->cut_vec->ind[j];
1.892 + xassert(1 <= k && k <= m+n);
1.893 + if (!mir->isint[k] && mir->cut_vec->val[j] > 0.0)
1.894 + mir->cut_vec->val[j] = 0.0;
1.895 + }
1.896 + ios_clean_vec(mir->cut_vec, 0.0);
1.897 +#if _MIR_DEBUG
1.898 + ios_check_vec(mir->cut_vec);
1.899 +#endif
1.900 + /* move integer terms to the beginning of the sparse vector and
1.901 + determine the number of integer variables */
1.902 + nint = 0;
1.903 + for (j = 1; j <= mir->cut_vec->nnz; j++)
1.904 + { k = mir->cut_vec->ind[j];
1.905 + xassert(1 <= k && k <= m+n);
1.906 + if (mir->isint[k])
1.907 + { double temp;
1.908 + nint++;
1.909 + /* interchange elements [nint] and [j] */
1.910 + kk = mir->cut_vec->ind[nint];
1.911 + mir->cut_vec->pos[k] = nint;
1.912 + mir->cut_vec->pos[kk] = j;
1.913 + mir->cut_vec->ind[nint] = k;
1.914 + mir->cut_vec->ind[j] = kk;
1.915 + temp = mir->cut_vec->val[nint];
1.916 + mir->cut_vec->val[nint] = mir->cut_vec->val[j];
1.917 + mir->cut_vec->val[j] = temp;
1.918 + }
1.919 + }
1.920 +#if _MIR_DEBUG
1.921 + ios_check_vec(mir->cut_vec);
1.922 +#endif
1.923 + /* if there is no integer variable, nothing to generate */
1.924 + if (nint == 0) goto done;
1.925 + /* allocate working arrays */
1.926 + u = xcalloc(1+nint, sizeof(double));
1.927 + x = xcalloc(1+nint, sizeof(double));
1.928 + alpha = xcalloc(1+nint, sizeof(double));
1.929 + /* determine u and x */
1.930 + for (j = 1; j <= nint; j++)
1.931 + { k = mir->cut_vec->ind[j];
1.932 + xassert(m+1 <= k && k <= m+n);
1.933 + xassert(mir->isint[k]);
1.934 + u[j] = mir->ub[k] - mir->lb[k];
1.935 + xassert(u[j] >= 1.0);
1.936 + if (mir->subst[k] == 'L')
1.937 + x[j] = mir->x[k] - mir->lb[k];
1.938 + else if (mir->subst[k] == 'U')
1.939 + x[j] = mir->ub[k] - mir->x[k];
1.940 + else
1.941 + xassert(k != k);
1.942 + xassert(x[j] >= -0.001);
1.943 + if (x[j] < 0.0) x[j] = 0.0;
1.944 + }
1.945 + /* compute s = - sum of continuous terms */
1.946 + s = 0.0;
1.947 + for (j = nint+1; j <= mir->cut_vec->nnz; j++)
1.948 + { double x;
1.949 + k = mir->cut_vec->ind[j];
1.950 + xassert(1 <= k && k <= m+n);
1.951 + /* must be continuous */
1.952 + xassert(!mir->isint[k]);
1.953 + if (mir->subst[k] == 'L')
1.954 + { xassert(mir->lb[k] != -DBL_MAX);
1.955 + kk = mir->vlb[k];
1.956 + if (kk == 0)
1.957 + x = mir->x[k] - mir->lb[k];
1.958 + else
1.959 + x = mir->x[k] - mir->lb[k] * mir->x[kk];
1.960 + }
1.961 + else if (mir->subst[k] == 'U')
1.962 + { xassert(mir->ub[k] != +DBL_MAX);
1.963 + kk = mir->vub[k];
1.964 + if (kk == 0)
1.965 + x = mir->ub[k] - mir->x[k];
1.966 + else
1.967 + x = mir->ub[k] * mir->x[kk] - mir->x[k];
1.968 + }
1.969 + else
1.970 + xassert(k != k);
1.971 + xassert(x >= -0.001);
1.972 + if (x < 0.0) x = 0.0;
1.973 + s -= mir->cut_vec->val[j] * x;
1.974 + }
1.975 + xassert(s >= 0.0);
1.976 + /* apply heuristic to obtain most violated c-MIR inequality */
1.977 + b = mir->cut_rhs;
1.978 + r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,
1.979 + &beta, &gamma);
1.980 + if (r_best == 0.0) goto skip;
1.981 + xassert(r_best > 0.0);
1.982 + /* convert to raw cut */
1.983 + /* sum alpha[j] * x[j] <= beta + gamma * s */
1.984 + for (j = 1; j <= nint; j++)
1.985 + mir->cut_vec->val[j] = alpha[j];
1.986 + for (j = nint+1; j <= mir->cut_vec->nnz; j++)
1.987 + { k = mir->cut_vec->ind[j];
1.988 + if (k <= m+n) mir->cut_vec->val[j] *= gamma;
1.989 + }
1.990 + mir->cut_rhs = beta;
1.991 +#if _MIR_DEBUG
1.992 + ios_check_vec(mir->cut_vec);
1.993 +#endif
1.994 +skip: /* free working arrays */
1.995 + xfree(u);
1.996 + xfree(x);
1.997 + xfree(alpha);
1.998 +done: return r_best;
1.999 +}
1.1000 +
1.1001 +#if _MIR_DEBUG
1.1002 +static void check_raw_cut(struct MIR *mir, double r_best)
1.1003 +{ /* check raw cut before back bound substitution */
1.1004 + int m = mir->m;
1.1005 + int n = mir->n;
1.1006 + int j, k, kk;
1.1007 + double r, big, x;
1.1008 + /* compute the residual r = sum a[k] * x[k] - b and determine
1.1009 + big = max(1, |a[k]|, |b|) */
1.1010 + r = 0.0, big = 1.0;
1.1011 + for (j = 1; j <= mir->cut_vec->nnz; j++)
1.1012 + { k = mir->cut_vec->ind[j];
1.1013 + xassert(1 <= k && k <= m+n);
1.1014 + if (mir->subst[k] == 'L')
1.1015 + { xassert(mir->lb[k] != -DBL_MAX);
1.1016 + kk = mir->vlb[k];
1.1017 + if (kk == 0)
1.1018 + x = mir->x[k] - mir->lb[k];
1.1019 + else
1.1020 + x = mir->x[k] - mir->lb[k] * mir->x[kk];
1.1021 + }
1.1022 + else if (mir->subst[k] == 'U')
1.1023 + { xassert(mir->ub[k] != +DBL_MAX);
1.1024 + kk = mir->vub[k];
1.1025 + if (kk == 0)
1.1026 + x = mir->ub[k] - mir->x[k];
1.1027 + else
1.1028 + x = mir->ub[k] * mir->x[kk] - mir->x[k];
1.1029 + }
1.1030 + else
1.1031 + xassert(k != k);
1.1032 + r += mir->cut_vec->val[j] * x;
1.1033 + if (big < fabs(mir->cut_vec->val[j]))
1.1034 + big = fabs(mir->cut_vec->val[j]);
1.1035 + }
1.1036 + r -= mir->cut_rhs;
1.1037 + if (big < fabs(mir->cut_rhs))
1.1038 + big = fabs(mir->cut_rhs);
1.1039 + /* the residual must be close to r_best */
1.1040 + xassert(fabs(r - r_best) <= 1e-6 * big);
1.1041 + return;
1.1042 +}
1.1043 +#endif
1.1044 +
1.1045 +static void back_subst(struct MIR *mir)
1.1046 +{ /* back substitution of original bounds */
1.1047 + int m = mir->m;
1.1048 + int n = mir->n;
1.1049 + int j, jj, k, kk;
1.1050 + /* at first, restore bounds of integer variables (because on
1.1051 + restoring variable bounds of continuous variables we need
1.1052 + original, not shifted, bounds of integer variables) */
1.1053 + for (j = 1; j <= mir->cut_vec->nnz; j++)
1.1054 + { k = mir->cut_vec->ind[j];
1.1055 + xassert(1 <= k && k <= m+n);
1.1056 + if (!mir->isint[k]) continue; /* skip continuous */
1.1057 + if (mir->subst[k] == 'L')
1.1058 + { /* x'[k] = x[k] - lb[k] */
1.1059 + xassert(mir->lb[k] != -DBL_MAX);
1.1060 + xassert(mir->vlb[k] == 0);
1.1061 + mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
1.1062 + }
1.1063 + else if (mir->subst[k] == 'U')
1.1064 + { /* x'[k] = ub[k] - x[k] */
1.1065 + xassert(mir->ub[k] != +DBL_MAX);
1.1066 + xassert(mir->vub[k] == 0);
1.1067 + mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
1.1068 + mir->cut_vec->val[j] = - mir->cut_vec->val[j];
1.1069 + }
1.1070 + else
1.1071 + xassert(k != k);
1.1072 + }
1.1073 + /* now restore bounds of continuous variables */
1.1074 + for (j = 1; j <= mir->cut_vec->nnz; j++)
1.1075 + { k = mir->cut_vec->ind[j];
1.1076 + xassert(1 <= k && k <= m+n);
1.1077 + if (mir->isint[k]) continue; /* skip integer */
1.1078 + if (mir->subst[k] == 'L')
1.1079 + { /* x'[k] = x[k] - (lower bound) */
1.1080 + xassert(mir->lb[k] != -DBL_MAX);
1.1081 + kk = mir->vlb[k];
1.1082 + if (kk == 0)
1.1083 + { /* x'[k] = x[k] - lb[k] */
1.1084 + mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
1.1085 + }
1.1086 + else
1.1087 + { /* x'[k] = x[k] - lb[k] * x[kk] */
1.1088 + jj = mir->cut_vec->pos[kk];
1.1089 +#if 0
1.1090 + xassert(jj != 0);
1.1091 +#else
1.1092 + if (jj == 0)
1.1093 + { ios_set_vj(mir->cut_vec, kk, 1.0);
1.1094 + jj = mir->cut_vec->pos[kk];
1.1095 + xassert(jj != 0);
1.1096 + mir->cut_vec->val[jj] = 0.0;
1.1097 + }
1.1098 +#endif
1.1099 + mir->cut_vec->val[jj] -= mir->cut_vec->val[j] *
1.1100 + mir->lb[k];
1.1101 + }
1.1102 + }
1.1103 + else if (mir->subst[k] == 'U')
1.1104 + { /* x'[k] = (upper bound) - x[k] */
1.1105 + xassert(mir->ub[k] != +DBL_MAX);
1.1106 + kk = mir->vub[k];
1.1107 + if (kk == 0)
1.1108 + { /* x'[k] = ub[k] - x[k] */
1.1109 + mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
1.1110 + }
1.1111 + else
1.1112 + { /* x'[k] = ub[k] * x[kk] - x[k] */
1.1113 + jj = mir->cut_vec->pos[kk];
1.1114 + if (jj == 0)
1.1115 + { ios_set_vj(mir->cut_vec, kk, 1.0);
1.1116 + jj = mir->cut_vec->pos[kk];
1.1117 + xassert(jj != 0);
1.1118 + mir->cut_vec->val[jj] = 0.0;
1.1119 + }
1.1120 + mir->cut_vec->val[jj] += mir->cut_vec->val[j] *
1.1121 + mir->ub[k];
1.1122 + }
1.1123 + mir->cut_vec->val[j] = - mir->cut_vec->val[j];
1.1124 + }
1.1125 + else
1.1126 + xassert(k != k);
1.1127 + }
1.1128 +#if _MIR_DEBUG
1.1129 + ios_check_vec(mir->cut_vec);
1.1130 +#endif
1.1131 + return;
1.1132 +}
1.1133 +
1.1134 +#if _MIR_DEBUG
1.1135 +static void check_cut_row(struct MIR *mir, double r_best)
1.1136 +{ /* check the cut after back bound substitution or elimination of
1.1137 + auxiliary variables */
1.1138 + int m = mir->m;
1.1139 + int n = mir->n;
1.1140 + int j, k;
1.1141 + double r, big;
1.1142 + /* compute the residual r = sum a[k] * x[k] - b and determine
1.1143 + big = max(1, |a[k]|, |b|) */
1.1144 + r = 0.0, big = 1.0;
1.1145 + for (j = 1; j <= mir->cut_vec->nnz; j++)
1.1146 + { k = mir->cut_vec->ind[j];
1.1147 + xassert(1 <= k && k <= m+n);
1.1148 + r += mir->cut_vec->val[j] * mir->x[k];
1.1149 + if (big < fabs(mir->cut_vec->val[j]))
1.1150 + big = fabs(mir->cut_vec->val[j]);
1.1151 + }
1.1152 + r -= mir->cut_rhs;
1.1153 + if (big < fabs(mir->cut_rhs))
1.1154 + big = fabs(mir->cut_rhs);
1.1155 + /* the residual must be close to r_best */
1.1156 + xassert(fabs(r - r_best) <= 1e-6 * big);
1.1157 + return;
1.1158 +}
1.1159 +#endif
1.1160 +
1.1161 +static void subst_aux_vars(glp_tree *tree, struct MIR *mir)
1.1162 +{ /* final substitution to eliminate auxiliary variables */
1.1163 + glp_prob *mip = tree->mip;
1.1164 + int m = mir->m;
1.1165 + int n = mir->n;
1.1166 + GLPAIJ *aij;
1.1167 + int j, k, kk, jj;
1.1168 + for (j = mir->cut_vec->nnz; j >= 1; j--)
1.1169 + { k = mir->cut_vec->ind[j];
1.1170 + xassert(1 <= k && k <= m+n);
1.1171 + if (k > m) continue; /* skip structurals */
1.1172 + for (aij = mip->row[k]->ptr; aij != NULL; aij = aij->r_next)
1.1173 + { kk = m + aij->col->j; /* structural */
1.1174 + jj = mir->cut_vec->pos[kk];
1.1175 + if (jj == 0)
1.1176 + { ios_set_vj(mir->cut_vec, kk, 1.0);
1.1177 + jj = mir->cut_vec->pos[kk];
1.1178 + mir->cut_vec->val[jj] = 0.0;
1.1179 + }
1.1180 + mir->cut_vec->val[jj] += mir->cut_vec->val[j] * aij->val;
1.1181 + }
1.1182 + mir->cut_vec->val[j] = 0.0;
1.1183 + }
1.1184 + ios_clean_vec(mir->cut_vec, 0.0);
1.1185 + return;
1.1186 +}
1.1187 +
1.1188 +static void add_cut(glp_tree *tree, struct MIR *mir)
1.1189 +{ /* add constructed cut inequality to the cut pool */
1.1190 + int m = mir->m;
1.1191 + int n = mir->n;
1.1192 + int j, k, len;
1.1193 + int *ind = xcalloc(1+n, sizeof(int));
1.1194 + double *val = xcalloc(1+n, sizeof(double));
1.1195 + len = 0;
1.1196 + for (j = mir->cut_vec->nnz; j >= 1; j--)
1.1197 + { k = mir->cut_vec->ind[j];
1.1198 + xassert(m+1 <= k && k <= m+n);
1.1199 + len++, ind[len] = k - m, val[len] = mir->cut_vec->val[j];
1.1200 + }
1.1201 +#if 0
1.1202 + ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP,
1.1203 + mir->cut_rhs);
1.1204 +#else
1.1205 + glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
1.1206 + mir->cut_rhs);
1.1207 +#endif
1.1208 + xfree(ind);
1.1209 + xfree(val);
1.1210 + return;
1.1211 +}
1.1212 +
1.1213 +static int aggregate_row(glp_tree *tree, struct MIR *mir)
1.1214 +{ /* try to aggregate another row */
1.1215 + glp_prob *mip = tree->mip;
1.1216 + int m = mir->m;
1.1217 + int n = mir->n;
1.1218 + GLPAIJ *aij;
1.1219 + IOSVEC *v;
1.1220 + int ii, j, jj, k, kk, kappa = 0, ret = 0;
1.1221 + double d1, d2, d, d_max = 0.0;
1.1222 + /* choose appropriate structural variable in the aggregated row
1.1223 + to be substituted */
1.1224 + for (j = 1; j <= mir->agg_vec->nnz; j++)
1.1225 + { k = mir->agg_vec->ind[j];
1.1226 + xassert(1 <= k && k <= m+n);
1.1227 + if (k <= m) continue; /* skip auxiliary var */
1.1228 + if (mir->isint[k]) continue; /* skip integer var */
1.1229 + if (fabs(mir->agg_vec->val[j]) < 0.001) continue;
1.1230 + /* compute distance from x[k] to its lower bound */
1.1231 + kk = mir->vlb[k];
1.1232 + if (kk == 0)
1.1233 + { if (mir->lb[k] == -DBL_MAX)
1.1234 + d1 = DBL_MAX;
1.1235 + else
1.1236 + d1 = mir->x[k] - mir->lb[k];
1.1237 + }
1.1238 + else
1.1239 + { xassert(1 <= kk && kk <= m+n);
1.1240 + xassert(mir->isint[kk]);
1.1241 + xassert(mir->lb[k] != -DBL_MAX);
1.1242 + d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
1.1243 + }
1.1244 + /* compute distance from x[k] to its upper bound */
1.1245 + kk = mir->vub[k];
1.1246 + if (kk == 0)
1.1247 + { if (mir->vub[k] == +DBL_MAX)
1.1248 + d2 = DBL_MAX;
1.1249 + else
1.1250 + d2 = mir->ub[k] - mir->x[k];
1.1251 + }
1.1252 + else
1.1253 + { xassert(1 <= kk && kk <= m+n);
1.1254 + xassert(mir->isint[kk]);
1.1255 + xassert(mir->ub[k] != +DBL_MAX);
1.1256 + d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
1.1257 + }
1.1258 + /* x[k] cannot be free */
1.1259 + xassert(d1 != DBL_MAX || d2 != DBL_MAX);
1.1260 + /* d = min(d1, d2) */
1.1261 + d = (d1 <= d2 ? d1 : d2);
1.1262 + xassert(d != DBL_MAX);
1.1263 + /* should not be close to corresponding bound */
1.1264 + if (d < 0.001) continue;
1.1265 + if (d_max < d) d_max = d, kappa = k;
1.1266 + }
1.1267 + if (kappa == 0)
1.1268 + { /* nothing chosen */
1.1269 + ret = 1;
1.1270 + goto done;
1.1271 + }
1.1272 + /* x[kappa] has been chosen */
1.1273 + xassert(m+1 <= kappa && kappa <= m+n);
1.1274 + xassert(!mir->isint[kappa]);
1.1275 + /* find another row, which have not been used yet, to eliminate
1.1276 + x[kappa] from the aggregated row */
1.1277 + for (ii = 1; ii <= m; ii++)
1.1278 + { if (mir->skip[ii]) continue;
1.1279 + for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
1.1280 + if (aij->col->j == kappa - m) break;
1.1281 + if (aij != NULL && fabs(aij->val) >= 0.001) break;
1.1282 + }
1.1283 + if (ii > m)
1.1284 + { /* nothing found */
1.1285 + ret = 2;
1.1286 + goto done;
1.1287 + }
1.1288 + /* row ii has been found; include it in the aggregated list */
1.1289 + mir->agg_cnt++;
1.1290 + xassert(mir->agg_cnt <= MAXAGGR);
1.1291 + mir->agg_row[mir->agg_cnt] = ii;
1.1292 + mir->skip[ii] = 2;
1.1293 + /* v := new row */
1.1294 + v = ios_create_vec(m+n);
1.1295 + ios_set_vj(v, ii, 1.0);
1.1296 + for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
1.1297 + ios_set_vj(v, m + aij->col->j, - aij->val);
1.1298 +#if _MIR_DEBUG
1.1299 + ios_check_vec(v);
1.1300 +#endif
1.1301 + /* perform gaussian elimination to remove x[kappa] */
1.1302 + j = mir->agg_vec->pos[kappa];
1.1303 + xassert(j != 0);
1.1304 + jj = v->pos[kappa];
1.1305 + xassert(jj != 0);
1.1306 + ios_linear_comb(mir->agg_vec,
1.1307 + - mir->agg_vec->val[j] / v->val[jj], v);
1.1308 + ios_delete_vec(v);
1.1309 + ios_set_vj(mir->agg_vec, kappa, 0.0);
1.1310 +#if _MIR_DEBUG
1.1311 + ios_check_vec(mir->agg_vec);
1.1312 +#endif
1.1313 +done: return ret;
1.1314 +}
1.1315 +
1.1316 +void ios_mir_gen(glp_tree *tree, void *gen)
1.1317 +{ /* main routine to generate MIR cuts */
1.1318 + glp_prob *mip = tree->mip;
1.1319 + struct MIR *mir = gen;
1.1320 + int m = mir->m;
1.1321 + int n = mir->n;
1.1322 + int i;
1.1323 + double r_best;
1.1324 + xassert(mip->m >= m);
1.1325 + xassert(mip->n == n);
1.1326 + /* obtain current point */
1.1327 + get_current_point(tree, mir);
1.1328 +#if _MIR_DEBUG
1.1329 + /* check current point */
1.1330 + check_current_point(mir);
1.1331 +#endif
1.1332 + /* reset bound substitution flags */
1.1333 + memset(&mir->subst[1], '?', m+n);
1.1334 + /* try to generate a set of violated MIR cuts */
1.1335 + for (i = 1; i <= m; i++)
1.1336 + { if (mir->skip[i]) continue;
1.1337 + /* use original i-th row as initial aggregated constraint */
1.1338 + initial_agg_row(tree, mir, i);
1.1339 +loop: ;
1.1340 +#if _MIR_DEBUG
1.1341 + /* check aggregated row */
1.1342 + check_agg_row(mir);
1.1343 +#endif
1.1344 + /* substitute fixed variables into aggregated constraint */
1.1345 + subst_fixed_vars(mir);
1.1346 +#if _MIR_DEBUG
1.1347 + /* check aggregated row */
1.1348 + check_agg_row(mir);
1.1349 +#endif
1.1350 +#if _MIR_DEBUG
1.1351 + /* check bound substitution flags */
1.1352 + { int k;
1.1353 + for (k = 1; k <= m+n; k++)
1.1354 + xassert(mir->subst[k] == '?');
1.1355 + }
1.1356 +#endif
1.1357 + /* apply bound substitution heuristic */
1.1358 + bound_subst_heur(mir);
1.1359 + /* substitute bounds and build modified constraint */
1.1360 + build_mod_row(mir);
1.1361 +#if _MIR_DEBUG
1.1362 + /* check modified row */
1.1363 + check_mod_row(mir);
1.1364 +#endif
1.1365 + /* try to generate violated c-MIR cut for modified row */
1.1366 + r_best = generate(mir);
1.1367 + if (r_best > 0.0)
1.1368 + { /* success */
1.1369 +#if _MIR_DEBUG
1.1370 + /* check raw cut before back bound substitution */
1.1371 + check_raw_cut(mir, r_best);
1.1372 +#endif
1.1373 + /* back substitution of original bounds */
1.1374 + back_subst(mir);
1.1375 +#if _MIR_DEBUG
1.1376 + /* check the cut after back bound substitution */
1.1377 + check_cut_row(mir, r_best);
1.1378 +#endif
1.1379 + /* final substitution to eliminate auxiliary variables */
1.1380 + subst_aux_vars(tree, mir);
1.1381 +#if _MIR_DEBUG
1.1382 + /* check the cut after elimination of auxiliaries */
1.1383 + check_cut_row(mir, r_best);
1.1384 +#endif
1.1385 + /* add constructed cut inequality to the cut pool */
1.1386 + add_cut(tree, mir);
1.1387 + }
1.1388 + /* reset bound substitution flags */
1.1389 + { int j, k;
1.1390 + for (j = 1; j <= mir->mod_vec->nnz; j++)
1.1391 + { k = mir->mod_vec->ind[j];
1.1392 + xassert(1 <= k && k <= m+n);
1.1393 + xassert(mir->subst[k] != '?');
1.1394 + mir->subst[k] = '?';
1.1395 + }
1.1396 + }
1.1397 + if (r_best == 0.0)
1.1398 + { /* failure */
1.1399 + if (mir->agg_cnt < MAXAGGR)
1.1400 + { /* try to aggregate another row */
1.1401 + if (aggregate_row(tree, mir) == 0) goto loop;
1.1402 + }
1.1403 + }
1.1404 + /* unmark rows used in the aggregated constraint */
1.1405 + { int k, ii;
1.1406 + for (k = 1; k <= mir->agg_cnt; k++)
1.1407 + { ii = mir->agg_row[k];
1.1408 + xassert(1 <= ii && ii <= m);
1.1409 + xassert(mir->skip[ii] == 2);
1.1410 + mir->skip[ii] = 0;
1.1411 + }
1.1412 + }
1.1413 + }
1.1414 + return;
1.1415 +}
1.1416 +
1.1417 +/***********************************************************************
1.1418 +* NAME
1.1419 +*
1.1420 +* ios_mir_term - terminate MIR cut generator
1.1421 +*
1.1422 +* SYNOPSIS
1.1423 +*
1.1424 +* #include "glpios.h"
1.1425 +* void ios_mir_term(void *gen);
1.1426 +*
1.1427 +* DESCRIPTION
1.1428 +*
1.1429 +* The routine ios_mir_term deletes the MIR cut generator working area
1.1430 +* freeing all the memory allocated to it. */
1.1431 +
1.1432 +void ios_mir_term(void *gen)
1.1433 +{ struct MIR *mir = gen;
1.1434 + xfree(mir->skip);
1.1435 + xfree(mir->isint);
1.1436 + xfree(mir->lb);
1.1437 + xfree(mir->vlb);
1.1438 + xfree(mir->ub);
1.1439 + xfree(mir->vub);
1.1440 + xfree(mir->x);
1.1441 + xfree(mir->agg_row);
1.1442 + ios_delete_vec(mir->agg_vec);
1.1443 + xfree(mir->subst);
1.1444 + ios_delete_vec(mir->mod_vec);
1.1445 + ios_delete_vec(mir->cut_vec);
1.1446 + xfree(mir);
1.1447 + return;
1.1448 +}
1.1449 +
1.1450 +/* eof */