1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpios10.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,348 @@
1.4 +/* glpios10.c (feasibility pump heuristic) */
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 +#include "glprng.h"
1.30 +
1.31 +/***********************************************************************
1.32 +* NAME
1.33 +*
1.34 +* ios_feas_pump - feasibility pump heuristic
1.35 +*
1.36 +* SYNOPSIS
1.37 +*
1.38 +* #include "glpios.h"
1.39 +* void ios_feas_pump(glp_tree *T);
1.40 +*
1.41 +* DESCRIPTION
1.42 +*
1.43 +* The routine ios_feas_pump is a simple implementation of the Feasi-
1.44 +* bility Pump heuristic.
1.45 +*
1.46 +* REFERENCES
1.47 +*
1.48 +* M.Fischetti, F.Glover, and A.Lodi. "The feasibility pump." Math.
1.49 +* Program., Ser. A 104, pp. 91-104 (2005). */
1.50 +
1.51 +struct VAR
1.52 +{ /* binary variable */
1.53 + int j;
1.54 + /* ordinal number */
1.55 + int x;
1.56 + /* value in the rounded solution (0 or 1) */
1.57 + double d;
1.58 + /* sorting key */
1.59 +};
1.60 +
1.61 +static int fcmp(const void *x, const void *y)
1.62 +{ /* comparison routine */
1.63 + const struct VAR *vx = x, *vy = y;
1.64 + if (vx->d > vy->d)
1.65 + return -1;
1.66 + else if (vx->d < vy->d)
1.67 + return +1;
1.68 + else
1.69 + return 0;
1.70 +}
1.71 +
1.72 +void ios_feas_pump(glp_tree *T)
1.73 +{ glp_prob *P = T->mip;
1.74 + int n = P->n;
1.75 + glp_prob *lp = NULL;
1.76 + struct VAR *var = NULL;
1.77 + RNG *rand = NULL;
1.78 + GLPCOL *col;
1.79 + glp_smcp parm;
1.80 + int j, k, new_x, nfail, npass, nv, ret, stalling;
1.81 + double dist, tol;
1.82 + xassert(glp_get_status(P) == GLP_OPT);
1.83 + /* this heuristic is applied only once on the root level */
1.84 + if (!(T->curr->level == 0 && T->curr->solved == 1)) goto done;
1.85 + /* determine number of binary variables */
1.86 + nv = 0;
1.87 + for (j = 1; j <= n; j++)
1.88 + { col = P->col[j];
1.89 + /* if x[j] is continuous, skip it */
1.90 + if (col->kind == GLP_CV) continue;
1.91 + /* if x[j] is fixed, skip it */
1.92 + if (col->type == GLP_FX) continue;
1.93 + /* x[j] is non-fixed integer */
1.94 + xassert(col->kind == GLP_IV);
1.95 + if (col->type == GLP_DB && col->lb == 0.0 && col->ub == 1.0)
1.96 + { /* x[j] is binary */
1.97 + nv++;
1.98 + }
1.99 + else
1.100 + { /* x[j] is general integer */
1.101 + if (T->parm->msg_lev >= GLP_MSG_ALL)
1.102 + xprintf("FPUMP heuristic cannot be applied due to genera"
1.103 + "l integer variables\n");
1.104 + goto done;
1.105 + }
1.106 + }
1.107 + /* there must be at least one binary variable */
1.108 + if (nv == 0) goto done;
1.109 + if (T->parm->msg_lev >= GLP_MSG_ALL)
1.110 + xprintf("Applying FPUMP heuristic...\n");
1.111 + /* build the list of binary variables */
1.112 + var = xcalloc(1+nv, sizeof(struct VAR));
1.113 + k = 0;
1.114 + for (j = 1; j <= n; j++)
1.115 + { col = P->col[j];
1.116 + if (col->kind == GLP_IV && col->type == GLP_DB)
1.117 + var[++k].j = j;
1.118 + }
1.119 + xassert(k == nv);
1.120 + /* create working problem object */
1.121 + lp = glp_create_prob();
1.122 +more: /* copy the original problem object to keep it intact */
1.123 + glp_copy_prob(lp, P, GLP_OFF);
1.124 + /* we are interested to find an integer feasible solution, which
1.125 + is better than the best known one */
1.126 + if (P->mip_stat == GLP_FEAS)
1.127 + { int *ind;
1.128 + double *val, bnd;
1.129 + /* add a row and make it identical to the objective row */
1.130 + glp_add_rows(lp, 1);
1.131 + ind = xcalloc(1+n, sizeof(int));
1.132 + val = xcalloc(1+n, sizeof(double));
1.133 + for (j = 1; j <= n; j++)
1.134 + { ind[j] = j;
1.135 + val[j] = P->col[j]->coef;
1.136 + }
1.137 + glp_set_mat_row(lp, lp->m, n, ind, val);
1.138 + xfree(ind);
1.139 + xfree(val);
1.140 + /* introduce upper (minimization) or lower (maximization)
1.141 + bound to the original objective function; note that this
1.142 + additional constraint is not violated at the optimal point
1.143 + to LP relaxation */
1.144 +#if 0 /* modified by xypron <xypron.glpk@gmx.de> */
1.145 + if (P->dir == GLP_MIN)
1.146 + { bnd = P->mip_obj - 0.10 * (1.0 + fabs(P->mip_obj));
1.147 + if (bnd < P->obj_val) bnd = P->obj_val;
1.148 + glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0);
1.149 + }
1.150 + else if (P->dir == GLP_MAX)
1.151 + { bnd = P->mip_obj + 0.10 * (1.0 + fabs(P->mip_obj));
1.152 + if (bnd > P->obj_val) bnd = P->obj_val;
1.153 + glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0);
1.154 + }
1.155 + else
1.156 + xassert(P != P);
1.157 +#else
1.158 + bnd = 0.1 * P->obj_val + 0.9 * P->mip_obj;
1.159 + /* xprintf("bnd = %f\n", bnd); */
1.160 + if (P->dir == GLP_MIN)
1.161 + glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0);
1.162 + else if (P->dir == GLP_MAX)
1.163 + glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0);
1.164 + else
1.165 + xassert(P != P);
1.166 +#endif
1.167 + }
1.168 + /* reset pass count */
1.169 + npass = 0;
1.170 + /* invalidate the rounded point */
1.171 + for (k = 1; k <= nv; k++)
1.172 + var[k].x = -1;
1.173 +pass: /* next pass starts here */
1.174 + npass++;
1.175 + if (T->parm->msg_lev >= GLP_MSG_ALL)
1.176 + xprintf("Pass %d\n", npass);
1.177 + /* initialize minimal distance between the basic point and the
1.178 + rounded one obtained during this pass */
1.179 + dist = DBL_MAX;
1.180 + /* reset failure count (the number of succeeded iterations failed
1.181 + to improve the distance) */
1.182 + nfail = 0;
1.183 + /* if it is not the first pass, perturb the last rounded point
1.184 + rather than construct it from the basic solution */
1.185 + if (npass > 1)
1.186 + { double rho, temp;
1.187 + if (rand == NULL)
1.188 + rand = rng_create_rand();
1.189 + for (k = 1; k <= nv; k++)
1.190 + { j = var[k].j;
1.191 + col = lp->col[j];
1.192 + rho = rng_uniform(rand, -0.3, 0.7);
1.193 + if (rho < 0.0) rho = 0.0;
1.194 + temp = fabs((double)var[k].x - col->prim);
1.195 + if (temp + rho > 0.5) var[k].x = 1 - var[k].x;
1.196 + }
1.197 + goto skip;
1.198 + }
1.199 +loop: /* innermost loop begins here */
1.200 + /* round basic solution (which is assumed primal feasible) */
1.201 + stalling = 1;
1.202 + for (k = 1; k <= nv; k++)
1.203 + { col = lp->col[var[k].j];
1.204 + if (col->prim < 0.5)
1.205 + { /* rounded value is 0 */
1.206 + new_x = 0;
1.207 + }
1.208 + else
1.209 + { /* rounded value is 1 */
1.210 + new_x = 1;
1.211 + }
1.212 + if (var[k].x != new_x)
1.213 + { stalling = 0;
1.214 + var[k].x = new_x;
1.215 + }
1.216 + }
1.217 + /* if the rounded point has not changed (stalling), choose and
1.218 + flip some its entries heuristically */
1.219 + if (stalling)
1.220 + { /* compute d[j] = |x[j] - round(x[j])| */
1.221 + for (k = 1; k <= nv; k++)
1.222 + { col = lp->col[var[k].j];
1.223 + var[k].d = fabs(col->prim - (double)var[k].x);
1.224 + }
1.225 + /* sort the list of binary variables by descending d[j] */
1.226 + qsort(&var[1], nv, sizeof(struct VAR), fcmp);
1.227 + /* choose and flip some rounded components */
1.228 + for (k = 1; k <= nv; k++)
1.229 + { if (k >= 5 && var[k].d < 0.35 || k >= 10) break;
1.230 + var[k].x = 1 - var[k].x;
1.231 + }
1.232 + }
1.233 +skip: /* check if the time limit has been exhausted */
1.234 + if (T->parm->tm_lim < INT_MAX &&
1.235 + (double)(T->parm->tm_lim - 1) <=
1.236 + 1000.0 * xdifftime(xtime(), T->tm_beg)) goto done;
1.237 + /* build the objective, which is the distance between the current
1.238 + (basic) point and the rounded one */
1.239 + lp->dir = GLP_MIN;
1.240 + lp->c0 = 0.0;
1.241 + for (j = 1; j <= n; j++)
1.242 + lp->col[j]->coef = 0.0;
1.243 + for (k = 1; k <= nv; k++)
1.244 + { j = var[k].j;
1.245 + if (var[k].x == 0)
1.246 + lp->col[j]->coef = +1.0;
1.247 + else
1.248 + { lp->col[j]->coef = -1.0;
1.249 + lp->c0 += 1.0;
1.250 + }
1.251 + }
1.252 + /* minimize the distance with the simplex method */
1.253 + glp_init_smcp(&parm);
1.254 + if (T->parm->msg_lev <= GLP_MSG_ERR)
1.255 + parm.msg_lev = T->parm->msg_lev;
1.256 + else if (T->parm->msg_lev <= GLP_MSG_ALL)
1.257 + { parm.msg_lev = GLP_MSG_ON;
1.258 + parm.out_dly = 10000;
1.259 + }
1.260 + ret = glp_simplex(lp, &parm);
1.261 + if (ret != 0)
1.262 + { if (T->parm->msg_lev >= GLP_MSG_ERR)
1.263 + xprintf("Warning: glp_simplex returned %d\n", ret);
1.264 + goto done;
1.265 + }
1.266 + ret = glp_get_status(lp);
1.267 + if (ret != GLP_OPT)
1.268 + { if (T->parm->msg_lev >= GLP_MSG_ERR)
1.269 + xprintf("Warning: glp_get_status returned %d\n", ret);
1.270 + goto done;
1.271 + }
1.272 + if (T->parm->msg_lev >= GLP_MSG_DBG)
1.273 + xprintf("delta = %g\n", lp->obj_val);
1.274 + /* check if the basic solution is integer feasible; note that it
1.275 + may be so even if the minimial distance is positive */
1.276 + tol = 0.3 * T->parm->tol_int;
1.277 + for (k = 1; k <= nv; k++)
1.278 + { col = lp->col[var[k].j];
1.279 + if (tol < col->prim && col->prim < 1.0 - tol) break;
1.280 + }
1.281 + if (k > nv)
1.282 + { /* okay; the basic solution seems to be integer feasible */
1.283 + double *x = xcalloc(1+n, sizeof(double));
1.284 + for (j = 1; j <= n; j++)
1.285 + { x[j] = lp->col[j]->prim;
1.286 + if (P->col[j]->kind == GLP_IV) x[j] = floor(x[j] + 0.5);
1.287 + }
1.288 +#if 1 /* modified by xypron <xypron.glpk@gmx.de> */
1.289 + /* reset direction and right-hand side of objective */
1.290 + lp->c0 = P->c0;
1.291 + lp->dir = P->dir;
1.292 + /* fix integer variables */
1.293 + for (k = 1; k <= nv; k++)
1.294 + { lp->col[var[k].j]->lb = x[var[k].j];
1.295 + lp->col[var[k].j]->ub = x[var[k].j];
1.296 + lp->col[var[k].j]->type = GLP_FX;
1.297 + }
1.298 + /* copy original objective function */
1.299 + for (j = 1; j <= n; j++)
1.300 + lp->col[j]->coef = P->col[j]->coef;
1.301 + /* solve original LP and copy result */
1.302 + ret = glp_simplex(lp, &parm);
1.303 + if (ret != 0)
1.304 + { if (T->parm->msg_lev >= GLP_MSG_ERR)
1.305 + xprintf("Warning: glp_simplex returned %d\n", ret);
1.306 + goto done;
1.307 + }
1.308 + ret = glp_get_status(lp);
1.309 + if (ret != GLP_OPT)
1.310 + { if (T->parm->msg_lev >= GLP_MSG_ERR)
1.311 + xprintf("Warning: glp_get_status returned %d\n", ret);
1.312 + goto done;
1.313 + }
1.314 + for (j = 1; j <= n; j++)
1.315 + if (P->col[j]->kind != GLP_IV) x[j] = lp->col[j]->prim;
1.316 +#endif
1.317 + ret = glp_ios_heur_sol(T, x);
1.318 + xfree(x);
1.319 + if (ret == 0)
1.320 + { /* the integer solution is accepted */
1.321 + if (ios_is_hopeful(T, T->curr->bound))
1.322 + { /* it is reasonable to apply the heuristic once again */
1.323 + goto more;
1.324 + }
1.325 + else
1.326 + { /* the best known integer feasible solution just found
1.327 + is close to optimal solution to LP relaxation */
1.328 + goto done;
1.329 + }
1.330 + }
1.331 + }
1.332 + /* the basic solution is fractional */
1.333 + if (dist == DBL_MAX ||
1.334 + lp->obj_val <= dist - 1e-6 * (1.0 + dist))
1.335 + { /* the distance is reducing */
1.336 + nfail = 0, dist = lp->obj_val;
1.337 + }
1.338 + else
1.339 + { /* improving the distance failed */
1.340 + nfail++;
1.341 + }
1.342 + if (nfail < 3) goto loop;
1.343 + if (npass < 5) goto pass;
1.344 +done: /* delete working objects */
1.345 + if (lp != NULL) glp_delete_prob(lp);
1.346 + if (var != NULL) xfree(var);
1.347 + if (rand != NULL) rng_delete_rand(rand);
1.348 + return;
1.349 +}
1.350 +
1.351 +/* eof */