You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1448 lines
48 KiB
1448 lines
48 KiB
/* glpios06.c (MIR cut generator) */
|
|
|
|
/***********************************************************************
|
|
* This code is part of GLPK (GNU Linear Programming Kit).
|
|
*
|
|
* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
|
|
* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied
|
|
* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
|
|
* reserved. E-mail: <mao@gnu.org>.
|
|
*
|
|
* GLPK is free software: you can redistribute it and/or modify it
|
|
* under the terms of the GNU General Public License as published by
|
|
* the Free Software Foundation, either version 3 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* GLPK is distributed in the hope that it will be useful, but WITHOUT
|
|
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|
* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
|
|
* License for more details.
|
|
*
|
|
* You should have received a copy of the GNU General Public License
|
|
* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
|
|
***********************************************************************/
|
|
|
|
#include "env.h"
|
|
#include "glpios.h"
|
|
|
|
#define _MIR_DEBUG 0
|
|
|
|
#define MAXAGGR 5
|
|
/* maximal number of rows which can be aggregated */
|
|
|
|
struct MIR
|
|
{ /* MIR cut generator working area */
|
|
/*--------------------------------------------------------------*/
|
|
/* global information valid for the root subproblem */
|
|
int m;
|
|
/* number of rows (in the root subproblem) */
|
|
int n;
|
|
/* number of columns */
|
|
char *skip; /* char skip[1+m]; */
|
|
/* skip[i], 1 <= i <= m, is a flag that means that row i should
|
|
not be used because (1) it is not suitable, or (2) because it
|
|
has been used in the aggregated constraint */
|
|
char *isint; /* char isint[1+m+n]; */
|
|
/* isint[k], 1 <= k <= m+n, is a flag that means that variable
|
|
x[k] is integer (otherwise, continuous) */
|
|
double *lb; /* double lb[1+m+n]; */
|
|
/* lb[k], 1 <= k <= m+n, is lower bound of x[k]; -DBL_MAX means
|
|
that x[k] has no lower bound */
|
|
int *vlb; /* int vlb[1+m+n]; */
|
|
/* vlb[k] = k', 1 <= k <= m+n, is the number of integer variable,
|
|
which defines variable lower bound x[k] >= lb[k] * x[k']; zero
|
|
means that x[k] has simple lower bound */
|
|
double *ub; /* double ub[1+m+n]; */
|
|
/* ub[k], 1 <= k <= m+n, is upper bound of x[k]; +DBL_MAX means
|
|
that x[k] has no upper bound */
|
|
int *vub; /* int vub[1+m+n]; */
|
|
/* vub[k] = k', 1 <= k <= m+n, is the number of integer variable,
|
|
which defines variable upper bound x[k] <= ub[k] * x[k']; zero
|
|
means that x[k] has simple upper bound */
|
|
/*--------------------------------------------------------------*/
|
|
/* current (fractional) point to be separated */
|
|
double *x; /* double x[1+m+n]; */
|
|
/* x[k] is current value of auxiliary (1 <= k <= m) or structural
|
|
(m+1 <= k <= m+n) variable */
|
|
/*--------------------------------------------------------------*/
|
|
/* aggregated constraint sum a[k] * x[k] = b, which is a linear
|
|
combination of original constraints transformed to equalities
|
|
by introducing auxiliary variables */
|
|
int agg_cnt;
|
|
/* number of rows (original constraints) used to build aggregated
|
|
constraint, 1 <= agg_cnt <= MAXAGGR */
|
|
int *agg_row; /* int agg_row[1+MAXAGGR]; */
|
|
/* agg_row[k], 1 <= k <= agg_cnt, is the row number used to build
|
|
aggregated constraint */
|
|
IOSVEC *agg_vec; /* IOSVEC agg_vec[1:m+n]; */
|
|
/* sparse vector of aggregated constraint coefficients, a[k] */
|
|
double agg_rhs;
|
|
/* right-hand side of the aggregated constraint, b */
|
|
/*--------------------------------------------------------------*/
|
|
/* bound substitution flags for modified constraint */
|
|
char *subst; /* char subst[1+m+n]; */
|
|
/* subst[k], 1 <= k <= m+n, is a bound substitution flag used for
|
|
variable x[k]:
|
|
'?' - x[k] is missing in modified constraint
|
|
'L' - x[k] = (lower bound) + x'[k]
|
|
'U' - x[k] = (upper bound) - x'[k] */
|
|
/*--------------------------------------------------------------*/
|
|
/* modified constraint sum a'[k] * x'[k] = b', where x'[k] >= 0,
|
|
derived from aggregated constraint by substituting bounds;
|
|
note that due to substitution of variable bounds there may be
|
|
additional terms in the modified constraint */
|
|
IOSVEC *mod_vec; /* IOSVEC mod_vec[1:m+n]; */
|
|
/* sparse vector of modified constraint coefficients, a'[k] */
|
|
double mod_rhs;
|
|
/* right-hand side of the modified constraint, b' */
|
|
/*--------------------------------------------------------------*/
|
|
/* cutting plane sum alpha[k] * x[k] <= beta */
|
|
IOSVEC *cut_vec; /* IOSVEC cut_vec[1:m+n]; */
|
|
/* sparse vector of cutting plane coefficients, alpha[k] */
|
|
double cut_rhs;
|
|
/* right-hand size of the cutting plane, beta */
|
|
};
|
|
|
|
/***********************************************************************
|
|
* NAME
|
|
*
|
|
* ios_mir_init - initialize MIR cut generator
|
|
*
|
|
* SYNOPSIS
|
|
*
|
|
* #include "glpios.h"
|
|
* void *ios_mir_init(glp_tree *tree);
|
|
*
|
|
* DESCRIPTION
|
|
*
|
|
* The routine ios_mir_init initializes the MIR cut generator assuming
|
|
* that the current subproblem is the root subproblem.
|
|
*
|
|
* RETURNS
|
|
*
|
|
* The routine ios_mir_init returns a pointer to the MIR cut generator
|
|
* working area. */
|
|
|
|
static void set_row_attrib(glp_tree *tree, struct MIR *mir)
|
|
{ /* set global row attributes */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mir->m;
|
|
int k;
|
|
for (k = 1; k <= m; k++)
|
|
{ GLPROW *row = mip->row[k];
|
|
mir->skip[k] = 0;
|
|
mir->isint[k] = 0;
|
|
switch (row->type)
|
|
{ case GLP_FR:
|
|
mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
|
|
case GLP_LO:
|
|
mir->lb[k] = row->lb, mir->ub[k] = +DBL_MAX; break;
|
|
case GLP_UP:
|
|
mir->lb[k] = -DBL_MAX, mir->ub[k] = row->ub; break;
|
|
case GLP_DB:
|
|
mir->lb[k] = row->lb, mir->ub[k] = row->ub; break;
|
|
case GLP_FX:
|
|
mir->lb[k] = mir->ub[k] = row->lb; break;
|
|
default:
|
|
xassert(row != row);
|
|
}
|
|
mir->vlb[k] = mir->vub[k] = 0;
|
|
}
|
|
return;
|
|
}
|
|
|
|
static void set_col_attrib(glp_tree *tree, struct MIR *mir)
|
|
{ /* set global column attributes */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int k;
|
|
for (k = m+1; k <= m+n; k++)
|
|
{ GLPCOL *col = mip->col[k-m];
|
|
switch (col->kind)
|
|
{ case GLP_CV:
|
|
mir->isint[k] = 0; break;
|
|
case GLP_IV:
|
|
mir->isint[k] = 1; break;
|
|
default:
|
|
xassert(col != col);
|
|
}
|
|
switch (col->type)
|
|
{ case GLP_FR:
|
|
mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
|
|
case GLP_LO:
|
|
mir->lb[k] = col->lb, mir->ub[k] = +DBL_MAX; break;
|
|
case GLP_UP:
|
|
mir->lb[k] = -DBL_MAX, mir->ub[k] = col->ub; break;
|
|
case GLP_DB:
|
|
mir->lb[k] = col->lb, mir->ub[k] = col->ub; break;
|
|
case GLP_FX:
|
|
mir->lb[k] = mir->ub[k] = col->lb; break;
|
|
default:
|
|
xassert(col != col);
|
|
}
|
|
mir->vlb[k] = mir->vub[k] = 0;
|
|
}
|
|
return;
|
|
}
|
|
|
|
static void set_var_bounds(glp_tree *tree, struct MIR *mir)
|
|
{ /* set variable bounds */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mir->m;
|
|
GLPAIJ *aij;
|
|
int i, k1, k2;
|
|
double a1, a2;
|
|
for (i = 1; i <= m; i++)
|
|
{ /* we need the row to be '>= 0' or '<= 0' */
|
|
if (!(mir->lb[i] == 0.0 && mir->ub[i] == +DBL_MAX ||
|
|
mir->lb[i] == -DBL_MAX && mir->ub[i] == 0.0)) continue;
|
|
/* take first term */
|
|
aij = mip->row[i]->ptr;
|
|
if (aij == NULL) continue;
|
|
k1 = m + aij->col->j, a1 = aij->val;
|
|
/* take second term */
|
|
aij = aij->r_next;
|
|
if (aij == NULL) continue;
|
|
k2 = m + aij->col->j, a2 = aij->val;
|
|
/* there must be only two terms */
|
|
if (aij->r_next != NULL) continue;
|
|
/* interchange terms, if needed */
|
|
if (!mir->isint[k1] && mir->isint[k2])
|
|
;
|
|
else if (mir->isint[k1] && !mir->isint[k2])
|
|
{ k2 = k1, a2 = a1;
|
|
k1 = m + aij->col->j, a1 = aij->val;
|
|
}
|
|
else
|
|
{ /* both terms are either continuous or integer */
|
|
continue;
|
|
}
|
|
/* x[k2] should be double-bounded */
|
|
if (mir->lb[k2] == -DBL_MAX || mir->ub[k2] == +DBL_MAX ||
|
|
mir->lb[k2] == mir->ub[k2]) continue;
|
|
/* change signs, if necessary */
|
|
if (mir->ub[i] == 0.0) a1 = - a1, a2 = - a2;
|
|
/* now the row has the form a1 * x1 + a2 * x2 >= 0, where x1
|
|
is continuous, x2 is integer */
|
|
if (a1 > 0.0)
|
|
{ /* x1 >= - (a2 / a1) * x2 */
|
|
if (mir->vlb[k1] == 0)
|
|
{ /* set variable lower bound for x1 */
|
|
mir->lb[k1] = - a2 / a1;
|
|
mir->vlb[k1] = k2;
|
|
/* the row should not be used */
|
|
mir->skip[i] = 1;
|
|
}
|
|
}
|
|
else /* a1 < 0.0 */
|
|
{ /* x1 <= - (a2 / a1) * x2 */
|
|
if (mir->vub[k1] == 0)
|
|
{ /* set variable upper bound for x1 */
|
|
mir->ub[k1] = - a2 / a1;
|
|
mir->vub[k1] = k2;
|
|
/* the row should not be used */
|
|
mir->skip[i] = 1;
|
|
}
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
static void mark_useless_rows(glp_tree *tree, struct MIR *mir)
|
|
{ /* mark rows which should not be used */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mir->m;
|
|
GLPAIJ *aij;
|
|
int i, k, nv;
|
|
for (i = 1; i <= m; i++)
|
|
{ /* free rows should not be used */
|
|
if (mir->lb[i] == -DBL_MAX && mir->ub[i] == +DBL_MAX)
|
|
{ mir->skip[i] = 1;
|
|
continue;
|
|
}
|
|
nv = 0;
|
|
for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
|
|
{ k = m + aij->col->j;
|
|
/* rows with free variables should not be used */
|
|
if (mir->lb[k] == -DBL_MAX && mir->ub[k] == +DBL_MAX)
|
|
{ mir->skip[i] = 1;
|
|
break;
|
|
}
|
|
/* rows with integer variables having infinite (lower or
|
|
upper) bound should not be used */
|
|
if (mir->isint[k] && mir->lb[k] == -DBL_MAX ||
|
|
mir->isint[k] && mir->ub[k] == +DBL_MAX)
|
|
{ mir->skip[i] = 1;
|
|
break;
|
|
}
|
|
/* count non-fixed variables */
|
|
if (!(mir->vlb[k] == 0 && mir->vub[k] == 0 &&
|
|
mir->lb[k] == mir->ub[k])) nv++;
|
|
}
|
|
/* rows with all variables fixed should not be used */
|
|
if (nv == 0)
|
|
{ mir->skip[i] = 1;
|
|
continue;
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
void *ios_mir_init(glp_tree *tree)
|
|
{ /* initialize MIR cut generator */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mip->m;
|
|
int n = mip->n;
|
|
struct MIR *mir;
|
|
#if _MIR_DEBUG
|
|
xprintf("ios_mir_init: warning: debug mode enabled\n");
|
|
#endif
|
|
/* allocate working area */
|
|
mir = xmalloc(sizeof(struct MIR));
|
|
mir->m = m;
|
|
mir->n = n;
|
|
mir->skip = xcalloc(1+m, sizeof(char));
|
|
mir->isint = xcalloc(1+m+n, sizeof(char));
|
|
mir->lb = xcalloc(1+m+n, sizeof(double));
|
|
mir->vlb = xcalloc(1+m+n, sizeof(int));
|
|
mir->ub = xcalloc(1+m+n, sizeof(double));
|
|
mir->vub = xcalloc(1+m+n, sizeof(int));
|
|
mir->x = xcalloc(1+m+n, sizeof(double));
|
|
mir->agg_row = xcalloc(1+MAXAGGR, sizeof(int));
|
|
mir->agg_vec = ios_create_vec(m+n);
|
|
mir->subst = xcalloc(1+m+n, sizeof(char));
|
|
mir->mod_vec = ios_create_vec(m+n);
|
|
mir->cut_vec = ios_create_vec(m+n);
|
|
/* set global row attributes */
|
|
set_row_attrib(tree, mir);
|
|
/* set global column attributes */
|
|
set_col_attrib(tree, mir);
|
|
/* set variable bounds */
|
|
set_var_bounds(tree, mir);
|
|
/* mark rows which should not be used */
|
|
mark_useless_rows(tree, mir);
|
|
return mir;
|
|
}
|
|
|
|
/***********************************************************************
|
|
* NAME
|
|
*
|
|
* ios_mir_gen - generate MIR cuts
|
|
*
|
|
* SYNOPSIS
|
|
*
|
|
* #include "glpios.h"
|
|
* void ios_mir_gen(glp_tree *tree, void *gen, IOSPOOL *pool);
|
|
*
|
|
* DESCRIPTION
|
|
*
|
|
* The routine ios_mir_gen generates MIR cuts for the current point and
|
|
* adds them to the cut pool. */
|
|
|
|
static void get_current_point(glp_tree *tree, struct MIR *mir)
|
|
{ /* obtain current point */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int k;
|
|
for (k = 1; k <= m; k++)
|
|
mir->x[k] = mip->row[k]->prim;
|
|
for (k = m+1; k <= m+n; k++)
|
|
mir->x[k] = mip->col[k-m]->prim;
|
|
return;
|
|
}
|
|
|
|
#if _MIR_DEBUG
|
|
static void check_current_point(struct MIR *mir)
|
|
{ /* check current point */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int k, kk;
|
|
double lb, ub, eps;
|
|
for (k = 1; k <= m+n; k++)
|
|
{ /* determine lower bound */
|
|
lb = mir->lb[k];
|
|
kk = mir->vlb[k];
|
|
if (kk != 0)
|
|
{ xassert(lb != -DBL_MAX);
|
|
xassert(!mir->isint[k]);
|
|
xassert(mir->isint[kk]);
|
|
lb *= mir->x[kk];
|
|
}
|
|
/* check lower bound */
|
|
if (lb != -DBL_MAX)
|
|
{ eps = 1e-6 * (1.0 + fabs(lb));
|
|
xassert(mir->x[k] >= lb - eps);
|
|
}
|
|
/* determine upper bound */
|
|
ub = mir->ub[k];
|
|
kk = mir->vub[k];
|
|
if (kk != 0)
|
|
{ xassert(ub != +DBL_MAX);
|
|
xassert(!mir->isint[k]);
|
|
xassert(mir->isint[kk]);
|
|
ub *= mir->x[kk];
|
|
}
|
|
/* check upper bound */
|
|
if (ub != +DBL_MAX)
|
|
{ eps = 1e-6 * (1.0 + fabs(ub));
|
|
xassert(mir->x[k] <= ub + eps);
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
#endif
|
|
|
|
static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i)
|
|
{ /* use original i-th row as initial aggregated constraint */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mir->m;
|
|
GLPAIJ *aij;
|
|
xassert(1 <= i && i <= m);
|
|
xassert(!mir->skip[i]);
|
|
/* mark i-th row in order not to use it in the same aggregated
|
|
constraint */
|
|
mir->skip[i] = 2;
|
|
mir->agg_cnt = 1;
|
|
mir->agg_row[1] = i;
|
|
/* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary
|
|
variable of row i, x[m+j] are structural variables */
|
|
ios_clear_vec(mir->agg_vec);
|
|
ios_set_vj(mir->agg_vec, i, 1.0);
|
|
for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
|
|
ios_set_vj(mir->agg_vec, m + aij->col->j, - aij->val);
|
|
mir->agg_rhs = 0.0;
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->agg_vec);
|
|
#endif
|
|
return;
|
|
}
|
|
|
|
#if _MIR_DEBUG
|
|
static void check_agg_row(struct MIR *mir)
|
|
{ /* check aggregated constraint */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int j, k;
|
|
double r, big;
|
|
/* compute the residual r = sum a[k] * x[k] - b and determine
|
|
big = max(1, |a[k]|, |b|) */
|
|
r = 0.0, big = 1.0;
|
|
for (j = 1; j <= mir->agg_vec->nnz; j++)
|
|
{ k = mir->agg_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
r += mir->agg_vec->val[j] * mir->x[k];
|
|
if (big < fabs(mir->agg_vec->val[j]))
|
|
big = fabs(mir->agg_vec->val[j]);
|
|
}
|
|
r -= mir->agg_rhs;
|
|
if (big < fabs(mir->agg_rhs))
|
|
big = fabs(mir->agg_rhs);
|
|
/* the residual must be close to zero */
|
|
xassert(fabs(r) <= 1e-6 * big);
|
|
return;
|
|
}
|
|
#endif
|
|
|
|
static void subst_fixed_vars(struct MIR *mir)
|
|
{ /* substitute fixed variables into aggregated constraint */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int j, k;
|
|
for (j = 1; j <= mir->agg_vec->nnz; j++)
|
|
{ k = mir->agg_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (mir->vlb[k] == 0 && mir->vub[k] == 0 &&
|
|
mir->lb[k] == mir->ub[k])
|
|
{ /* x[k] is fixed */
|
|
mir->agg_rhs -= mir->agg_vec->val[j] * mir->lb[k];
|
|
mir->agg_vec->val[j] = 0.0;
|
|
}
|
|
}
|
|
/* remove terms corresponding to fixed variables */
|
|
ios_clean_vec(mir->agg_vec, DBL_EPSILON);
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->agg_vec);
|
|
#endif
|
|
return;
|
|
}
|
|
|
|
static void bound_subst_heur(struct MIR *mir)
|
|
{ /* bound substitution heuristic */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int j, k, kk;
|
|
double d1, d2;
|
|
for (j = 1; j <= mir->agg_vec->nnz; j++)
|
|
{ k = mir->agg_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (mir->isint[k]) continue; /* skip integer variable */
|
|
/* compute distance from x[k] to its lower bound */
|
|
kk = mir->vlb[k];
|
|
if (kk == 0)
|
|
{ if (mir->lb[k] == -DBL_MAX)
|
|
d1 = DBL_MAX;
|
|
else
|
|
d1 = mir->x[k] - mir->lb[k];
|
|
}
|
|
else
|
|
{ xassert(1 <= kk && kk <= m+n);
|
|
xassert(mir->isint[kk]);
|
|
xassert(mir->lb[k] != -DBL_MAX);
|
|
d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
|
|
}
|
|
/* compute distance from x[k] to its upper bound */
|
|
kk = mir->vub[k];
|
|
if (kk == 0)
|
|
{ if (mir->vub[k] == +DBL_MAX)
|
|
d2 = DBL_MAX;
|
|
else
|
|
d2 = mir->ub[k] - mir->x[k];
|
|
}
|
|
else
|
|
{ xassert(1 <= kk && kk <= m+n);
|
|
xassert(mir->isint[kk]);
|
|
xassert(mir->ub[k] != +DBL_MAX);
|
|
d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
|
|
}
|
|
/* x[k] cannot be free */
|
|
xassert(d1 != DBL_MAX || d2 != DBL_MAX);
|
|
/* choose the bound which is closer to x[k] */
|
|
xassert(mir->subst[k] == '?');
|
|
if (d1 <= d2)
|
|
mir->subst[k] = 'L';
|
|
else
|
|
mir->subst[k] = 'U';
|
|
}
|
|
return;
|
|
}
|
|
|
|
static void build_mod_row(struct MIR *mir)
|
|
{ /* substitute bounds and build modified constraint */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int j, jj, k, kk;
|
|
/* initially modified constraint is aggregated constraint */
|
|
ios_copy_vec(mir->mod_vec, mir->agg_vec);
|
|
mir->mod_rhs = mir->agg_rhs;
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->mod_vec);
|
|
#endif
|
|
/* substitute bounds for continuous variables; note that due to
|
|
substitution of variable bounds additional terms may appear in
|
|
modified constraint */
|
|
for (j = mir->mod_vec->nnz; j >= 1; j--)
|
|
{ k = mir->mod_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (mir->isint[k]) continue; /* skip integer variable */
|
|
if (mir->subst[k] == 'L')
|
|
{ /* x[k] = (lower bound) + x'[k] */
|
|
xassert(mir->lb[k] != -DBL_MAX);
|
|
kk = mir->vlb[k];
|
|
if (kk == 0)
|
|
{ /* x[k] = lb[k] + x'[k] */
|
|
mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
|
|
}
|
|
else
|
|
{ /* x[k] = lb[k] * x[kk] + x'[k] */
|
|
xassert(mir->isint[kk]);
|
|
jj = mir->mod_vec->pos[kk];
|
|
if (jj == 0)
|
|
{ ios_set_vj(mir->mod_vec, kk, 1.0);
|
|
jj = mir->mod_vec->pos[kk];
|
|
mir->mod_vec->val[jj] = 0.0;
|
|
}
|
|
mir->mod_vec->val[jj] +=
|
|
mir->mod_vec->val[j] * mir->lb[k];
|
|
}
|
|
}
|
|
else if (mir->subst[k] == 'U')
|
|
{ /* x[k] = (upper bound) - x'[k] */
|
|
xassert(mir->ub[k] != +DBL_MAX);
|
|
kk = mir->vub[k];
|
|
if (kk == 0)
|
|
{ /* x[k] = ub[k] - x'[k] */
|
|
mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
|
|
}
|
|
else
|
|
{ /* x[k] = ub[k] * x[kk] - x'[k] */
|
|
xassert(mir->isint[kk]);
|
|
jj = mir->mod_vec->pos[kk];
|
|
if (jj == 0)
|
|
{ ios_set_vj(mir->mod_vec, kk, 1.0);
|
|
jj = mir->mod_vec->pos[kk];
|
|
mir->mod_vec->val[jj] = 0.0;
|
|
}
|
|
mir->mod_vec->val[jj] +=
|
|
mir->mod_vec->val[j] * mir->ub[k];
|
|
}
|
|
mir->mod_vec->val[j] = - mir->mod_vec->val[j];
|
|
}
|
|
else
|
|
xassert(k != k);
|
|
}
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->mod_vec);
|
|
#endif
|
|
/* substitute bounds for integer variables */
|
|
for (j = 1; j <= mir->mod_vec->nnz; j++)
|
|
{ k = mir->mod_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (!mir->isint[k]) continue; /* skip continuous variable */
|
|
xassert(mir->subst[k] == '?');
|
|
xassert(mir->vlb[k] == 0 && mir->vub[k] == 0);
|
|
xassert(mir->lb[k] != -DBL_MAX && mir->ub[k] != +DBL_MAX);
|
|
if (fabs(mir->lb[k]) <= fabs(mir->ub[k]))
|
|
{ /* x[k] = lb[k] + x'[k] */
|
|
mir->subst[k] = 'L';
|
|
mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
|
|
}
|
|
else
|
|
{ /* x[k] = ub[k] - x'[k] */
|
|
mir->subst[k] = 'U';
|
|
mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
|
|
mir->mod_vec->val[j] = - mir->mod_vec->val[j];
|
|
}
|
|
}
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->mod_vec);
|
|
#endif
|
|
return;
|
|
}
|
|
|
|
#if _MIR_DEBUG
|
|
static void check_mod_row(struct MIR *mir)
|
|
{ /* check modified constraint */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int j, k, kk;
|
|
double r, big, x;
|
|
/* compute the residual r = sum a'[k] * x'[k] - b' and determine
|
|
big = max(1, |a[k]|, |b|) */
|
|
r = 0.0, big = 1.0;
|
|
for (j = 1; j <= mir->mod_vec->nnz; j++)
|
|
{ k = mir->mod_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (mir->subst[k] == 'L')
|
|
{ /* x'[k] = x[k] - (lower bound) */
|
|
xassert(mir->lb[k] != -DBL_MAX);
|
|
kk = mir->vlb[k];
|
|
if (kk == 0)
|
|
x = mir->x[k] - mir->lb[k];
|
|
else
|
|
x = mir->x[k] - mir->lb[k] * mir->x[kk];
|
|
}
|
|
else if (mir->subst[k] == 'U')
|
|
{ /* x'[k] = (upper bound) - x[k] */
|
|
xassert(mir->ub[k] != +DBL_MAX);
|
|
kk = mir->vub[k];
|
|
if (kk == 0)
|
|
x = mir->ub[k] - mir->x[k];
|
|
else
|
|
x = mir->ub[k] * mir->x[kk] - mir->x[k];
|
|
}
|
|
else
|
|
xassert(k != k);
|
|
r += mir->mod_vec->val[j] * x;
|
|
if (big < fabs(mir->mod_vec->val[j]))
|
|
big = fabs(mir->mod_vec->val[j]);
|
|
}
|
|
r -= mir->mod_rhs;
|
|
if (big < fabs(mir->mod_rhs))
|
|
big = fabs(mir->mod_rhs);
|
|
/* the residual must be close to zero */
|
|
xassert(fabs(r) <= 1e-6 * big);
|
|
return;
|
|
}
|
|
#endif
|
|
|
|
/***********************************************************************
|
|
* mir_ineq - construct MIR inequality
|
|
*
|
|
* Given the single constraint mixed integer set
|
|
*
|
|
* |N|
|
|
* X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s},
|
|
* + + j in N
|
|
*
|
|
* this routine constructs the mixed integer rounding (MIR) inequality
|
|
*
|
|
* sum alpha[j] * x[j] <= beta + gamma * s,
|
|
* j in N
|
|
*
|
|
* which is valid for X.
|
|
*
|
|
* If the MIR inequality has been successfully constructed, the routine
|
|
* returns zero. Otherwise, if b is close to nearest integer, there may
|
|
* be numeric difficulties due to big coefficients; so in this case the
|
|
* routine returns non-zero. */
|
|
|
|
static int mir_ineq(const int n, const double a[], const double b,
|
|
double alpha[], double *beta, double *gamma)
|
|
{ int j;
|
|
double f, t;
|
|
if (fabs(b - floor(b + .5)) < 0.01)
|
|
return 1;
|
|
f = b - floor(b);
|
|
for (j = 1; j <= n; j++)
|
|
{ t = (a[j] - floor(a[j])) - f;
|
|
if (t <= 0.0)
|
|
alpha[j] = floor(a[j]);
|
|
else
|
|
alpha[j] = floor(a[j]) + t / (1.0 - f);
|
|
}
|
|
*beta = floor(b);
|
|
*gamma = 1.0 / (1.0 - f);
|
|
return 0;
|
|
}
|
|
|
|
/***********************************************************************
|
|
* cmir_ineq - construct c-MIR inequality
|
|
*
|
|
* Given the mixed knapsack set
|
|
*
|
|
* MK |N|
|
|
* X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s,
|
|
* + + j in N
|
|
*
|
|
* x[j] <= u[j]},
|
|
*
|
|
* a subset C of variables to be complemented, and a divisor delta > 0,
|
|
* this routine constructs the complemented MIR (c-MIR) inequality
|
|
*
|
|
* sum alpha[j] * x[j] <= beta + gamma * s,
|
|
* j in N
|
|
* MK
|
|
* which is valid for X .
|
|
*
|
|
* If the c-MIR inequality has been successfully constructed, the
|
|
* routine returns zero. Otherwise, if there is a risk of numerical
|
|
* difficulties due to big coefficients (see comments to the routine
|
|
* mir_ineq), the routine cmir_ineq returns non-zero. */
|
|
|
|
static int cmir_ineq(const int n, const double a[], const double b,
|
|
const double u[], const char cset[], const double delta,
|
|
double alpha[], double *beta, double *gamma)
|
|
{ int j;
|
|
double *aa, bb;
|
|
aa = alpha, bb = b;
|
|
for (j = 1; j <= n; j++)
|
|
{ aa[j] = a[j] / delta;
|
|
if (cset[j])
|
|
aa[j] = - aa[j], bb -= a[j] * u[j];
|
|
}
|
|
bb /= delta;
|
|
if (mir_ineq(n, aa, bb, alpha, beta, gamma)) return 1;
|
|
for (j = 1; j <= n; j++)
|
|
{ if (cset[j])
|
|
alpha[j] = - alpha[j], *beta += alpha[j] * u[j];
|
|
}
|
|
*gamma /= delta;
|
|
return 0;
|
|
}
|
|
|
|
/***********************************************************************
|
|
* cmir_sep - c-MIR separation heuristic
|
|
*
|
|
* Given the mixed knapsack set
|
|
*
|
|
* MK |N|
|
|
* X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s,
|
|
* + + j in N
|
|
*
|
|
* x[j] <= u[j]}
|
|
*
|
|
* * *
|
|
* and a fractional point (x , s ), this routine tries to construct
|
|
* c-MIR inequality
|
|
*
|
|
* sum alpha[j] * x[j] <= beta + gamma * s,
|
|
* j in N
|
|
* MK
|
|
* which is valid for X and has (desirably maximal) violation at the
|
|
* fractional point given. This is attained by choosing an appropriate
|
|
* set C of variables to be complemented and a divisor delta > 0, which
|
|
* together define corresponding c-MIR inequality.
|
|
*
|
|
* If a violated c-MIR inequality has been successfully constructed,
|
|
* the routine returns its violation:
|
|
*
|
|
* * *
|
|
* sum alpha[j] * x [j] - beta - gamma * s ,
|
|
* j in N
|
|
*
|
|
* which is positive. In case of failure the routine returns zero. */
|
|
|
|
struct vset { int j; double v; };
|
|
|
|
static int cmir_cmp(const void *p1, const void *p2)
|
|
{ const struct vset *v1 = p1, *v2 = p2;
|
|
if (v1->v < v2->v) return -1;
|
|
if (v1->v > v2->v) return +1;
|
|
return 0;
|
|
}
|
|
|
|
static double cmir_sep(const int n, const double a[], const double b,
|
|
const double u[], const double x[], const double s,
|
|
double alpha[], double *beta, double *gamma)
|
|
{ int fail, j, k, nv, v;
|
|
double delta, eps, d_try[1+3], r, r_best;
|
|
char *cset;
|
|
struct vset *vset;
|
|
/* allocate working arrays */
|
|
cset = xcalloc(1+n, sizeof(char));
|
|
vset = xcalloc(1+n, sizeof(struct vset));
|
|
/* choose initial C */
|
|
for (j = 1; j <= n; j++)
|
|
cset[j] = (char)(x[j] >= 0.5 * u[j]);
|
|
/* choose initial delta */
|
|
r_best = delta = 0.0;
|
|
for (j = 1; j <= n; j++)
|
|
{ xassert(a[j] != 0.0);
|
|
/* if x[j] is close to its bounds, skip it */
|
|
eps = 1e-9 * (1.0 + fabs(u[j]));
|
|
if (x[j] < eps || x[j] > u[j] - eps) continue;
|
|
/* try delta = |a[j]| to construct c-MIR inequality */
|
|
fail = cmir_ineq(n, a, b, u, cset, fabs(a[j]), alpha, beta,
|
|
gamma);
|
|
if (fail) continue;
|
|
/* compute violation */
|
|
r = - (*beta) - (*gamma) * s;
|
|
for (k = 1; k <= n; k++) r += alpha[k] * x[k];
|
|
if (r_best < r) r_best = r, delta = fabs(a[j]);
|
|
}
|
|
if (r_best < 0.001) r_best = 0.0;
|
|
if (r_best == 0.0) goto done;
|
|
xassert(delta > 0.0);
|
|
/* try to increase violation by dividing delta by 2, 4, and 8,
|
|
respectively */
|
|
d_try[1] = delta / 2.0;
|
|
d_try[2] = delta / 4.0;
|
|
d_try[3] = delta / 8.0;
|
|
for (j = 1; j <= 3; j++)
|
|
{ /* construct c-MIR inequality */
|
|
fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,
|
|
gamma);
|
|
if (fail) continue;
|
|
/* compute violation */
|
|
r = - (*beta) - (*gamma) * s;
|
|
for (k = 1; k <= n; k++) r += alpha[k] * x[k];
|
|
if (r_best < r) r_best = r, delta = d_try[j];
|
|
}
|
|
/* build subset of variables lying strictly between their bounds
|
|
and order it by nondecreasing values of |x[j] - u[j]/2| */
|
|
nv = 0;
|
|
for (j = 1; j <= n; j++)
|
|
{ /* if x[j] is close to its bounds, skip it */
|
|
eps = 1e-9 * (1.0 + fabs(u[j]));
|
|
if (x[j] < eps || x[j] > u[j] - eps) continue;
|
|
/* add x[j] to the subset */
|
|
nv++;
|
|
vset[nv].j = j;
|
|
vset[nv].v = fabs(x[j] - 0.5 * u[j]);
|
|
}
|
|
qsort(&vset[1], nv, sizeof(struct vset), cmir_cmp);
|
|
/* try to increase violation by successively complementing each
|
|
variable in the subset */
|
|
for (v = 1; v <= nv; v++)
|
|
{ j = vset[v].j;
|
|
/* replace x[j] by its complement or vice versa */
|
|
cset[j] = (char)!cset[j];
|
|
/* construct c-MIR inequality */
|
|
fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
|
|
/* restore the variable */
|
|
cset[j] = (char)!cset[j];
|
|
/* do not replace the variable in case of failure */
|
|
if (fail) continue;
|
|
/* compute violation */
|
|
r = - (*beta) - (*gamma) * s;
|
|
for (k = 1; k <= n; k++) r += alpha[k] * x[k];
|
|
if (r_best < r) r_best = r, cset[j] = (char)!cset[j];
|
|
}
|
|
/* construct the best c-MIR inequality chosen */
|
|
fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
|
|
xassert(!fail);
|
|
done: /* free working arrays */
|
|
xfree(cset);
|
|
xfree(vset);
|
|
/* return to the calling routine */
|
|
return r_best;
|
|
}
|
|
|
|
static double generate(struct MIR *mir)
|
|
{ /* try to generate violated c-MIR cut for modified constraint */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int j, k, kk, nint;
|
|
double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma;
|
|
ios_copy_vec(mir->cut_vec, mir->mod_vec);
|
|
mir->cut_rhs = mir->mod_rhs;
|
|
/* remove small terms, which can appear due to substitution of
|
|
variable bounds */
|
|
ios_clean_vec(mir->cut_vec, DBL_EPSILON);
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->cut_vec);
|
|
#endif
|
|
/* remove positive continuous terms to obtain MK relaxation */
|
|
for (j = 1; j <= mir->cut_vec->nnz; j++)
|
|
{ k = mir->cut_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (!mir->isint[k] && mir->cut_vec->val[j] > 0.0)
|
|
mir->cut_vec->val[j] = 0.0;
|
|
}
|
|
ios_clean_vec(mir->cut_vec, 0.0);
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->cut_vec);
|
|
#endif
|
|
/* move integer terms to the beginning of the sparse vector and
|
|
determine the number of integer variables */
|
|
nint = 0;
|
|
for (j = 1; j <= mir->cut_vec->nnz; j++)
|
|
{ k = mir->cut_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (mir->isint[k])
|
|
{ double temp;
|
|
nint++;
|
|
/* interchange elements [nint] and [j] */
|
|
kk = mir->cut_vec->ind[nint];
|
|
mir->cut_vec->pos[k] = nint;
|
|
mir->cut_vec->pos[kk] = j;
|
|
mir->cut_vec->ind[nint] = k;
|
|
mir->cut_vec->ind[j] = kk;
|
|
temp = mir->cut_vec->val[nint];
|
|
mir->cut_vec->val[nint] = mir->cut_vec->val[j];
|
|
mir->cut_vec->val[j] = temp;
|
|
}
|
|
}
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->cut_vec);
|
|
#endif
|
|
/* if there is no integer variable, nothing to generate */
|
|
if (nint == 0) goto done;
|
|
/* allocate working arrays */
|
|
u = xcalloc(1+nint, sizeof(double));
|
|
x = xcalloc(1+nint, sizeof(double));
|
|
alpha = xcalloc(1+nint, sizeof(double));
|
|
/* determine u and x */
|
|
for (j = 1; j <= nint; j++)
|
|
{ k = mir->cut_vec->ind[j];
|
|
xassert(m+1 <= k && k <= m+n);
|
|
xassert(mir->isint[k]);
|
|
u[j] = mir->ub[k] - mir->lb[k];
|
|
xassert(u[j] >= 1.0);
|
|
if (mir->subst[k] == 'L')
|
|
x[j] = mir->x[k] - mir->lb[k];
|
|
else if (mir->subst[k] == 'U')
|
|
x[j] = mir->ub[k] - mir->x[k];
|
|
else
|
|
xassert(k != k);
|
|
xassert(x[j] >= -0.001);
|
|
if (x[j] < 0.0) x[j] = 0.0;
|
|
}
|
|
/* compute s = - sum of continuous terms */
|
|
s = 0.0;
|
|
for (j = nint+1; j <= mir->cut_vec->nnz; j++)
|
|
{ double x;
|
|
k = mir->cut_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
/* must be continuous */
|
|
xassert(!mir->isint[k]);
|
|
if (mir->subst[k] == 'L')
|
|
{ xassert(mir->lb[k] != -DBL_MAX);
|
|
kk = mir->vlb[k];
|
|
if (kk == 0)
|
|
x = mir->x[k] - mir->lb[k];
|
|
else
|
|
x = mir->x[k] - mir->lb[k] * mir->x[kk];
|
|
}
|
|
else if (mir->subst[k] == 'U')
|
|
{ xassert(mir->ub[k] != +DBL_MAX);
|
|
kk = mir->vub[k];
|
|
if (kk == 0)
|
|
x = mir->ub[k] - mir->x[k];
|
|
else
|
|
x = mir->ub[k] * mir->x[kk] - mir->x[k];
|
|
}
|
|
else
|
|
xassert(k != k);
|
|
xassert(x >= -0.001);
|
|
if (x < 0.0) x = 0.0;
|
|
s -= mir->cut_vec->val[j] * x;
|
|
}
|
|
xassert(s >= 0.0);
|
|
/* apply heuristic to obtain most violated c-MIR inequality */
|
|
b = mir->cut_rhs;
|
|
r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,
|
|
&beta, &gamma);
|
|
if (r_best == 0.0) goto skip;
|
|
xassert(r_best > 0.0);
|
|
/* convert to raw cut */
|
|
/* sum alpha[j] * x[j] <= beta + gamma * s */
|
|
for (j = 1; j <= nint; j++)
|
|
mir->cut_vec->val[j] = alpha[j];
|
|
for (j = nint+1; j <= mir->cut_vec->nnz; j++)
|
|
{ k = mir->cut_vec->ind[j];
|
|
if (k <= m+n) mir->cut_vec->val[j] *= gamma;
|
|
}
|
|
mir->cut_rhs = beta;
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->cut_vec);
|
|
#endif
|
|
skip: /* free working arrays */
|
|
xfree(u);
|
|
xfree(x);
|
|
xfree(alpha);
|
|
done: return r_best;
|
|
}
|
|
|
|
#if _MIR_DEBUG
|
|
static void check_raw_cut(struct MIR *mir, double r_best)
|
|
{ /* check raw cut before back bound substitution */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int j, k, kk;
|
|
double r, big, x;
|
|
/* compute the residual r = sum a[k] * x[k] - b and determine
|
|
big = max(1, |a[k]|, |b|) */
|
|
r = 0.0, big = 1.0;
|
|
for (j = 1; j <= mir->cut_vec->nnz; j++)
|
|
{ k = mir->cut_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (mir->subst[k] == 'L')
|
|
{ xassert(mir->lb[k] != -DBL_MAX);
|
|
kk = mir->vlb[k];
|
|
if (kk == 0)
|
|
x = mir->x[k] - mir->lb[k];
|
|
else
|
|
x = mir->x[k] - mir->lb[k] * mir->x[kk];
|
|
}
|
|
else if (mir->subst[k] == 'U')
|
|
{ xassert(mir->ub[k] != +DBL_MAX);
|
|
kk = mir->vub[k];
|
|
if (kk == 0)
|
|
x = mir->ub[k] - mir->x[k];
|
|
else
|
|
x = mir->ub[k] * mir->x[kk] - mir->x[k];
|
|
}
|
|
else
|
|
xassert(k != k);
|
|
r += mir->cut_vec->val[j] * x;
|
|
if (big < fabs(mir->cut_vec->val[j]))
|
|
big = fabs(mir->cut_vec->val[j]);
|
|
}
|
|
r -= mir->cut_rhs;
|
|
if (big < fabs(mir->cut_rhs))
|
|
big = fabs(mir->cut_rhs);
|
|
/* the residual must be close to r_best */
|
|
xassert(fabs(r - r_best) <= 1e-6 * big);
|
|
return;
|
|
}
|
|
#endif
|
|
|
|
static void back_subst(struct MIR *mir)
|
|
{ /* back substitution of original bounds */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int j, jj, k, kk;
|
|
/* at first, restore bounds of integer variables (because on
|
|
restoring variable bounds of continuous variables we need
|
|
original, not shifted, bounds of integer variables) */
|
|
for (j = 1; j <= mir->cut_vec->nnz; j++)
|
|
{ k = mir->cut_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (!mir->isint[k]) continue; /* skip continuous */
|
|
if (mir->subst[k] == 'L')
|
|
{ /* x'[k] = x[k] - lb[k] */
|
|
xassert(mir->lb[k] != -DBL_MAX);
|
|
xassert(mir->vlb[k] == 0);
|
|
mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
|
|
}
|
|
else if (mir->subst[k] == 'U')
|
|
{ /* x'[k] = ub[k] - x[k] */
|
|
xassert(mir->ub[k] != +DBL_MAX);
|
|
xassert(mir->vub[k] == 0);
|
|
mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
|
|
mir->cut_vec->val[j] = - mir->cut_vec->val[j];
|
|
}
|
|
else
|
|
xassert(k != k);
|
|
}
|
|
/* now restore bounds of continuous variables */
|
|
for (j = 1; j <= mir->cut_vec->nnz; j++)
|
|
{ k = mir->cut_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (mir->isint[k]) continue; /* skip integer */
|
|
if (mir->subst[k] == 'L')
|
|
{ /* x'[k] = x[k] - (lower bound) */
|
|
xassert(mir->lb[k] != -DBL_MAX);
|
|
kk = mir->vlb[k];
|
|
if (kk == 0)
|
|
{ /* x'[k] = x[k] - lb[k] */
|
|
mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
|
|
}
|
|
else
|
|
{ /* x'[k] = x[k] - lb[k] * x[kk] */
|
|
jj = mir->cut_vec->pos[kk];
|
|
#if 0
|
|
xassert(jj != 0);
|
|
#else
|
|
if (jj == 0)
|
|
{ ios_set_vj(mir->cut_vec, kk, 1.0);
|
|
jj = mir->cut_vec->pos[kk];
|
|
xassert(jj != 0);
|
|
mir->cut_vec->val[jj] = 0.0;
|
|
}
|
|
#endif
|
|
mir->cut_vec->val[jj] -= mir->cut_vec->val[j] *
|
|
mir->lb[k];
|
|
}
|
|
}
|
|
else if (mir->subst[k] == 'U')
|
|
{ /* x'[k] = (upper bound) - x[k] */
|
|
xassert(mir->ub[k] != +DBL_MAX);
|
|
kk = mir->vub[k];
|
|
if (kk == 0)
|
|
{ /* x'[k] = ub[k] - x[k] */
|
|
mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
|
|
}
|
|
else
|
|
{ /* x'[k] = ub[k] * x[kk] - x[k] */
|
|
jj = mir->cut_vec->pos[kk];
|
|
if (jj == 0)
|
|
{ ios_set_vj(mir->cut_vec, kk, 1.0);
|
|
jj = mir->cut_vec->pos[kk];
|
|
xassert(jj != 0);
|
|
mir->cut_vec->val[jj] = 0.0;
|
|
}
|
|
mir->cut_vec->val[jj] += mir->cut_vec->val[j] *
|
|
mir->ub[k];
|
|
}
|
|
mir->cut_vec->val[j] = - mir->cut_vec->val[j];
|
|
}
|
|
else
|
|
xassert(k != k);
|
|
}
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->cut_vec);
|
|
#endif
|
|
return;
|
|
}
|
|
|
|
#if _MIR_DEBUG
|
|
static void check_cut_row(struct MIR *mir, double r_best)
|
|
{ /* check the cut after back bound substitution or elimination of
|
|
auxiliary variables */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int j, k;
|
|
double r, big;
|
|
/* compute the residual r = sum a[k] * x[k] - b and determine
|
|
big = max(1, |a[k]|, |b|) */
|
|
r = 0.0, big = 1.0;
|
|
for (j = 1; j <= mir->cut_vec->nnz; j++)
|
|
{ k = mir->cut_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
r += mir->cut_vec->val[j] * mir->x[k];
|
|
if (big < fabs(mir->cut_vec->val[j]))
|
|
big = fabs(mir->cut_vec->val[j]);
|
|
}
|
|
r -= mir->cut_rhs;
|
|
if (big < fabs(mir->cut_rhs))
|
|
big = fabs(mir->cut_rhs);
|
|
/* the residual must be close to r_best */
|
|
xassert(fabs(r - r_best) <= 1e-6 * big);
|
|
return;
|
|
}
|
|
#endif
|
|
|
|
static void subst_aux_vars(glp_tree *tree, struct MIR *mir)
|
|
{ /* final substitution to eliminate auxiliary variables */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
GLPAIJ *aij;
|
|
int j, k, kk, jj;
|
|
for (j = mir->cut_vec->nnz; j >= 1; j--)
|
|
{ k = mir->cut_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (k > m) continue; /* skip structurals */
|
|
for (aij = mip->row[k]->ptr; aij != NULL; aij = aij->r_next)
|
|
{ kk = m + aij->col->j; /* structural */
|
|
jj = mir->cut_vec->pos[kk];
|
|
if (jj == 0)
|
|
{ ios_set_vj(mir->cut_vec, kk, 1.0);
|
|
jj = mir->cut_vec->pos[kk];
|
|
mir->cut_vec->val[jj] = 0.0;
|
|
}
|
|
mir->cut_vec->val[jj] += mir->cut_vec->val[j] * aij->val;
|
|
}
|
|
mir->cut_vec->val[j] = 0.0;
|
|
}
|
|
ios_clean_vec(mir->cut_vec, 0.0);
|
|
return;
|
|
}
|
|
|
|
static void add_cut(glp_tree *tree, struct MIR *mir)
|
|
{ /* add constructed cut inequality to the cut pool */
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int j, k, len;
|
|
int *ind = xcalloc(1+n, sizeof(int));
|
|
double *val = xcalloc(1+n, sizeof(double));
|
|
len = 0;
|
|
for (j = mir->cut_vec->nnz; j >= 1; j--)
|
|
{ k = mir->cut_vec->ind[j];
|
|
xassert(m+1 <= k && k <= m+n);
|
|
len++, ind[len] = k - m, val[len] = mir->cut_vec->val[j];
|
|
}
|
|
#if 0
|
|
ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP,
|
|
mir->cut_rhs);
|
|
#else
|
|
glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
|
|
mir->cut_rhs);
|
|
#endif
|
|
xfree(ind);
|
|
xfree(val);
|
|
return;
|
|
}
|
|
|
|
static int aggregate_row(glp_tree *tree, struct MIR *mir)
|
|
{ /* try to aggregate another row */
|
|
glp_prob *mip = tree->mip;
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
GLPAIJ *aij;
|
|
IOSVEC *v;
|
|
int ii, j, jj, k, kk, kappa = 0, ret = 0;
|
|
double d1, d2, d, d_max = 0.0;
|
|
/* choose appropriate structural variable in the aggregated row
|
|
to be substituted */
|
|
for (j = 1; j <= mir->agg_vec->nnz; j++)
|
|
{ k = mir->agg_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
if (k <= m) continue; /* skip auxiliary var */
|
|
if (mir->isint[k]) continue; /* skip integer var */
|
|
if (fabs(mir->agg_vec->val[j]) < 0.001) continue;
|
|
/* compute distance from x[k] to its lower bound */
|
|
kk = mir->vlb[k];
|
|
if (kk == 0)
|
|
{ if (mir->lb[k] == -DBL_MAX)
|
|
d1 = DBL_MAX;
|
|
else
|
|
d1 = mir->x[k] - mir->lb[k];
|
|
}
|
|
else
|
|
{ xassert(1 <= kk && kk <= m+n);
|
|
xassert(mir->isint[kk]);
|
|
xassert(mir->lb[k] != -DBL_MAX);
|
|
d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
|
|
}
|
|
/* compute distance from x[k] to its upper bound */
|
|
kk = mir->vub[k];
|
|
if (kk == 0)
|
|
{ if (mir->vub[k] == +DBL_MAX)
|
|
d2 = DBL_MAX;
|
|
else
|
|
d2 = mir->ub[k] - mir->x[k];
|
|
}
|
|
else
|
|
{ xassert(1 <= kk && kk <= m+n);
|
|
xassert(mir->isint[kk]);
|
|
xassert(mir->ub[k] != +DBL_MAX);
|
|
d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
|
|
}
|
|
/* x[k] cannot be free */
|
|
xassert(d1 != DBL_MAX || d2 != DBL_MAX);
|
|
/* d = min(d1, d2) */
|
|
d = (d1 <= d2 ? d1 : d2);
|
|
xassert(d != DBL_MAX);
|
|
/* should not be close to corresponding bound */
|
|
if (d < 0.001) continue;
|
|
if (d_max < d) d_max = d, kappa = k;
|
|
}
|
|
if (kappa == 0)
|
|
{ /* nothing chosen */
|
|
ret = 1;
|
|
goto done;
|
|
}
|
|
/* x[kappa] has been chosen */
|
|
xassert(m+1 <= kappa && kappa <= m+n);
|
|
xassert(!mir->isint[kappa]);
|
|
/* find another row, which have not been used yet, to eliminate
|
|
x[kappa] from the aggregated row */
|
|
for (ii = 1; ii <= m; ii++)
|
|
{ if (mir->skip[ii]) continue;
|
|
for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
|
|
if (aij->col->j == kappa - m) break;
|
|
if (aij != NULL && fabs(aij->val) >= 0.001) break;
|
|
}
|
|
if (ii > m)
|
|
{ /* nothing found */
|
|
ret = 2;
|
|
goto done;
|
|
}
|
|
/* row ii has been found; include it in the aggregated list */
|
|
mir->agg_cnt++;
|
|
xassert(mir->agg_cnt <= MAXAGGR);
|
|
mir->agg_row[mir->agg_cnt] = ii;
|
|
mir->skip[ii] = 2;
|
|
/* v := new row */
|
|
v = ios_create_vec(m+n);
|
|
ios_set_vj(v, ii, 1.0);
|
|
for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
|
|
ios_set_vj(v, m + aij->col->j, - aij->val);
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(v);
|
|
#endif
|
|
/* perform gaussian elimination to remove x[kappa] */
|
|
j = mir->agg_vec->pos[kappa];
|
|
xassert(j != 0);
|
|
jj = v->pos[kappa];
|
|
xassert(jj != 0);
|
|
ios_linear_comb(mir->agg_vec,
|
|
- mir->agg_vec->val[j] / v->val[jj], v);
|
|
ios_delete_vec(v);
|
|
ios_set_vj(mir->agg_vec, kappa, 0.0);
|
|
#if _MIR_DEBUG
|
|
ios_check_vec(mir->agg_vec);
|
|
#endif
|
|
done: return ret;
|
|
}
|
|
|
|
void ios_mir_gen(glp_tree *tree, void *gen)
|
|
{ /* main routine to generate MIR cuts */
|
|
glp_prob *mip = tree->mip;
|
|
struct MIR *mir = gen;
|
|
int m = mir->m;
|
|
int n = mir->n;
|
|
int i;
|
|
double r_best;
|
|
xassert(mip->m >= m);
|
|
xassert(mip->n == n);
|
|
/* obtain current point */
|
|
get_current_point(tree, mir);
|
|
#if _MIR_DEBUG
|
|
/* check current point */
|
|
check_current_point(mir);
|
|
#endif
|
|
/* reset bound substitution flags */
|
|
memset(&mir->subst[1], '?', m+n);
|
|
/* try to generate a set of violated MIR cuts */
|
|
for (i = 1; i <= m; i++)
|
|
{ if (mir->skip[i]) continue;
|
|
/* use original i-th row as initial aggregated constraint */
|
|
initial_agg_row(tree, mir, i);
|
|
loop: ;
|
|
#if _MIR_DEBUG
|
|
/* check aggregated row */
|
|
check_agg_row(mir);
|
|
#endif
|
|
/* substitute fixed variables into aggregated constraint */
|
|
subst_fixed_vars(mir);
|
|
#if _MIR_DEBUG
|
|
/* check aggregated row */
|
|
check_agg_row(mir);
|
|
#endif
|
|
#if _MIR_DEBUG
|
|
/* check bound substitution flags */
|
|
{ int k;
|
|
for (k = 1; k <= m+n; k++)
|
|
xassert(mir->subst[k] == '?');
|
|
}
|
|
#endif
|
|
/* apply bound substitution heuristic */
|
|
bound_subst_heur(mir);
|
|
/* substitute bounds and build modified constraint */
|
|
build_mod_row(mir);
|
|
#if _MIR_DEBUG
|
|
/* check modified row */
|
|
check_mod_row(mir);
|
|
#endif
|
|
/* try to generate violated c-MIR cut for modified row */
|
|
r_best = generate(mir);
|
|
if (r_best > 0.0)
|
|
{ /* success */
|
|
#if _MIR_DEBUG
|
|
/* check raw cut before back bound substitution */
|
|
check_raw_cut(mir, r_best);
|
|
#endif
|
|
/* back substitution of original bounds */
|
|
back_subst(mir);
|
|
#if _MIR_DEBUG
|
|
/* check the cut after back bound substitution */
|
|
check_cut_row(mir, r_best);
|
|
#endif
|
|
/* final substitution to eliminate auxiliary variables */
|
|
subst_aux_vars(tree, mir);
|
|
#if _MIR_DEBUG
|
|
/* check the cut after elimination of auxiliaries */
|
|
check_cut_row(mir, r_best);
|
|
#endif
|
|
/* add constructed cut inequality to the cut pool */
|
|
add_cut(tree, mir);
|
|
}
|
|
/* reset bound substitution flags */
|
|
{ int j, k;
|
|
for (j = 1; j <= mir->mod_vec->nnz; j++)
|
|
{ k = mir->mod_vec->ind[j];
|
|
xassert(1 <= k && k <= m+n);
|
|
xassert(mir->subst[k] != '?');
|
|
mir->subst[k] = '?';
|
|
}
|
|
}
|
|
if (r_best == 0.0)
|
|
{ /* failure */
|
|
if (mir->agg_cnt < MAXAGGR)
|
|
{ /* try to aggregate another row */
|
|
if (aggregate_row(tree, mir) == 0) goto loop;
|
|
}
|
|
}
|
|
/* unmark rows used in the aggregated constraint */
|
|
{ int k, ii;
|
|
for (k = 1; k <= mir->agg_cnt; k++)
|
|
{ ii = mir->agg_row[k];
|
|
xassert(1 <= ii && ii <= m);
|
|
xassert(mir->skip[ii] == 2);
|
|
mir->skip[ii] = 0;
|
|
}
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
/***********************************************************************
|
|
* NAME
|
|
*
|
|
* ios_mir_term - terminate MIR cut generator
|
|
*
|
|
* SYNOPSIS
|
|
*
|
|
* #include "glpios.h"
|
|
* void ios_mir_term(void *gen);
|
|
*
|
|
* DESCRIPTION
|
|
*
|
|
* The routine ios_mir_term deletes the MIR cut generator working area
|
|
* freeing all the memory allocated to it. */
|
|
|
|
void ios_mir_term(void *gen)
|
|
{ struct MIR *mir = gen;
|
|
xfree(mir->skip);
|
|
xfree(mir->isint);
|
|
xfree(mir->lb);
|
|
xfree(mir->vlb);
|
|
xfree(mir->ub);
|
|
xfree(mir->vub);
|
|
xfree(mir->x);
|
|
xfree(mir->agg_row);
|
|
ios_delete_vec(mir->agg_vec);
|
|
xfree(mir->subst);
|
|
ios_delete_vec(mir->mod_vec);
|
|
ios_delete_vec(mir->cut_vec);
|
|
xfree(mir);
|
|
return;
|
|
}
|
|
|
|
/* eof */
|