lemon-project-template-glpk

annotate deps/glpk/src/glpios10.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
rev   line source
alpar@9 1 /* glpios10.c (feasibility pump heuristic) */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpios.h"
alpar@9 26 #include "glprng.h"
alpar@9 27
alpar@9 28 /***********************************************************************
alpar@9 29 * NAME
alpar@9 30 *
alpar@9 31 * ios_feas_pump - feasibility pump heuristic
alpar@9 32 *
alpar@9 33 * SYNOPSIS
alpar@9 34 *
alpar@9 35 * #include "glpios.h"
alpar@9 36 * void ios_feas_pump(glp_tree *T);
alpar@9 37 *
alpar@9 38 * DESCRIPTION
alpar@9 39 *
alpar@9 40 * The routine ios_feas_pump is a simple implementation of the Feasi-
alpar@9 41 * bility Pump heuristic.
alpar@9 42 *
alpar@9 43 * REFERENCES
alpar@9 44 *
alpar@9 45 * M.Fischetti, F.Glover, and A.Lodi. "The feasibility pump." Math.
alpar@9 46 * Program., Ser. A 104, pp. 91-104 (2005). */
alpar@9 47
alpar@9 48 struct VAR
alpar@9 49 { /* binary variable */
alpar@9 50 int j;
alpar@9 51 /* ordinal number */
alpar@9 52 int x;
alpar@9 53 /* value in the rounded solution (0 or 1) */
alpar@9 54 double d;
alpar@9 55 /* sorting key */
alpar@9 56 };
alpar@9 57
alpar@9 58 static int fcmp(const void *x, const void *y)
alpar@9 59 { /* comparison routine */
alpar@9 60 const struct VAR *vx = x, *vy = y;
alpar@9 61 if (vx->d > vy->d)
alpar@9 62 return -1;
alpar@9 63 else if (vx->d < vy->d)
alpar@9 64 return +1;
alpar@9 65 else
alpar@9 66 return 0;
alpar@9 67 }
alpar@9 68
alpar@9 69 void ios_feas_pump(glp_tree *T)
alpar@9 70 { glp_prob *P = T->mip;
alpar@9 71 int n = P->n;
alpar@9 72 glp_prob *lp = NULL;
alpar@9 73 struct VAR *var = NULL;
alpar@9 74 RNG *rand = NULL;
alpar@9 75 GLPCOL *col;
alpar@9 76 glp_smcp parm;
alpar@9 77 int j, k, new_x, nfail, npass, nv, ret, stalling;
alpar@9 78 double dist, tol;
alpar@9 79 xassert(glp_get_status(P) == GLP_OPT);
alpar@9 80 /* this heuristic is applied only once on the root level */
alpar@9 81 if (!(T->curr->level == 0 && T->curr->solved == 1)) goto done;
alpar@9 82 /* determine number of binary variables */
alpar@9 83 nv = 0;
alpar@9 84 for (j = 1; j <= n; j++)
alpar@9 85 { col = P->col[j];
alpar@9 86 /* if x[j] is continuous, skip it */
alpar@9 87 if (col->kind == GLP_CV) continue;
alpar@9 88 /* if x[j] is fixed, skip it */
alpar@9 89 if (col->type == GLP_FX) continue;
alpar@9 90 /* x[j] is non-fixed integer */
alpar@9 91 xassert(col->kind == GLP_IV);
alpar@9 92 if (col->type == GLP_DB && col->lb == 0.0 && col->ub == 1.0)
alpar@9 93 { /* x[j] is binary */
alpar@9 94 nv++;
alpar@9 95 }
alpar@9 96 else
alpar@9 97 { /* x[j] is general integer */
alpar@9 98 if (T->parm->msg_lev >= GLP_MSG_ALL)
alpar@9 99 xprintf("FPUMP heuristic cannot be applied due to genera"
alpar@9 100 "l integer variables\n");
alpar@9 101 goto done;
alpar@9 102 }
alpar@9 103 }
alpar@9 104 /* there must be at least one binary variable */
alpar@9 105 if (nv == 0) goto done;
alpar@9 106 if (T->parm->msg_lev >= GLP_MSG_ALL)
alpar@9 107 xprintf("Applying FPUMP heuristic...\n");
alpar@9 108 /* build the list of binary variables */
alpar@9 109 var = xcalloc(1+nv, sizeof(struct VAR));
alpar@9 110 k = 0;
alpar@9 111 for (j = 1; j <= n; j++)
alpar@9 112 { col = P->col[j];
alpar@9 113 if (col->kind == GLP_IV && col->type == GLP_DB)
alpar@9 114 var[++k].j = j;
alpar@9 115 }
alpar@9 116 xassert(k == nv);
alpar@9 117 /* create working problem object */
alpar@9 118 lp = glp_create_prob();
alpar@9 119 more: /* copy the original problem object to keep it intact */
alpar@9 120 glp_copy_prob(lp, P, GLP_OFF);
alpar@9 121 /* we are interested to find an integer feasible solution, which
alpar@9 122 is better than the best known one */
alpar@9 123 if (P->mip_stat == GLP_FEAS)
alpar@9 124 { int *ind;
alpar@9 125 double *val, bnd;
alpar@9 126 /* add a row and make it identical to the objective row */
alpar@9 127 glp_add_rows(lp, 1);
alpar@9 128 ind = xcalloc(1+n, sizeof(int));
alpar@9 129 val = xcalloc(1+n, sizeof(double));
alpar@9 130 for (j = 1; j <= n; j++)
alpar@9 131 { ind[j] = j;
alpar@9 132 val[j] = P->col[j]->coef;
alpar@9 133 }
alpar@9 134 glp_set_mat_row(lp, lp->m, n, ind, val);
alpar@9 135 xfree(ind);
alpar@9 136 xfree(val);
alpar@9 137 /* introduce upper (minimization) or lower (maximization)
alpar@9 138 bound to the original objective function; note that this
alpar@9 139 additional constraint is not violated at the optimal point
alpar@9 140 to LP relaxation */
alpar@9 141 #if 0 /* modified by xypron <xypron.glpk@gmx.de> */
alpar@9 142 if (P->dir == GLP_MIN)
alpar@9 143 { bnd = P->mip_obj - 0.10 * (1.0 + fabs(P->mip_obj));
alpar@9 144 if (bnd < P->obj_val) bnd = P->obj_val;
alpar@9 145 glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0);
alpar@9 146 }
alpar@9 147 else if (P->dir == GLP_MAX)
alpar@9 148 { bnd = P->mip_obj + 0.10 * (1.0 + fabs(P->mip_obj));
alpar@9 149 if (bnd > P->obj_val) bnd = P->obj_val;
alpar@9 150 glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0);
alpar@9 151 }
alpar@9 152 else
alpar@9 153 xassert(P != P);
alpar@9 154 #else
alpar@9 155 bnd = 0.1 * P->obj_val + 0.9 * P->mip_obj;
alpar@9 156 /* xprintf("bnd = %f\n", bnd); */
alpar@9 157 if (P->dir == GLP_MIN)
alpar@9 158 glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0);
alpar@9 159 else if (P->dir == GLP_MAX)
alpar@9 160 glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0);
alpar@9 161 else
alpar@9 162 xassert(P != P);
alpar@9 163 #endif
alpar@9 164 }
alpar@9 165 /* reset pass count */
alpar@9 166 npass = 0;
alpar@9 167 /* invalidate the rounded point */
alpar@9 168 for (k = 1; k <= nv; k++)
alpar@9 169 var[k].x = -1;
alpar@9 170 pass: /* next pass starts here */
alpar@9 171 npass++;
alpar@9 172 if (T->parm->msg_lev >= GLP_MSG_ALL)
alpar@9 173 xprintf("Pass %d\n", npass);
alpar@9 174 /* initialize minimal distance between the basic point and the
alpar@9 175 rounded one obtained during this pass */
alpar@9 176 dist = DBL_MAX;
alpar@9 177 /* reset failure count (the number of succeeded iterations failed
alpar@9 178 to improve the distance) */
alpar@9 179 nfail = 0;
alpar@9 180 /* if it is not the first pass, perturb the last rounded point
alpar@9 181 rather than construct it from the basic solution */
alpar@9 182 if (npass > 1)
alpar@9 183 { double rho, temp;
alpar@9 184 if (rand == NULL)
alpar@9 185 rand = rng_create_rand();
alpar@9 186 for (k = 1; k <= nv; k++)
alpar@9 187 { j = var[k].j;
alpar@9 188 col = lp->col[j];
alpar@9 189 rho = rng_uniform(rand, -0.3, 0.7);
alpar@9 190 if (rho < 0.0) rho = 0.0;
alpar@9 191 temp = fabs((double)var[k].x - col->prim);
alpar@9 192 if (temp + rho > 0.5) var[k].x = 1 - var[k].x;
alpar@9 193 }
alpar@9 194 goto skip;
alpar@9 195 }
alpar@9 196 loop: /* innermost loop begins here */
alpar@9 197 /* round basic solution (which is assumed primal feasible) */
alpar@9 198 stalling = 1;
alpar@9 199 for (k = 1; k <= nv; k++)
alpar@9 200 { col = lp->col[var[k].j];
alpar@9 201 if (col->prim < 0.5)
alpar@9 202 { /* rounded value is 0 */
alpar@9 203 new_x = 0;
alpar@9 204 }
alpar@9 205 else
alpar@9 206 { /* rounded value is 1 */
alpar@9 207 new_x = 1;
alpar@9 208 }
alpar@9 209 if (var[k].x != new_x)
alpar@9 210 { stalling = 0;
alpar@9 211 var[k].x = new_x;
alpar@9 212 }
alpar@9 213 }
alpar@9 214 /* if the rounded point has not changed (stalling), choose and
alpar@9 215 flip some its entries heuristically */
alpar@9 216 if (stalling)
alpar@9 217 { /* compute d[j] = |x[j] - round(x[j])| */
alpar@9 218 for (k = 1; k <= nv; k++)
alpar@9 219 { col = lp->col[var[k].j];
alpar@9 220 var[k].d = fabs(col->prim - (double)var[k].x);
alpar@9 221 }
alpar@9 222 /* sort the list of binary variables by descending d[j] */
alpar@9 223 qsort(&var[1], nv, sizeof(struct VAR), fcmp);
alpar@9 224 /* choose and flip some rounded components */
alpar@9 225 for (k = 1; k <= nv; k++)
alpar@9 226 { if (k >= 5 && var[k].d < 0.35 || k >= 10) break;
alpar@9 227 var[k].x = 1 - var[k].x;
alpar@9 228 }
alpar@9 229 }
alpar@9 230 skip: /* check if the time limit has been exhausted */
alpar@9 231 if (T->parm->tm_lim < INT_MAX &&
alpar@9 232 (double)(T->parm->tm_lim - 1) <=
alpar@9 233 1000.0 * xdifftime(xtime(), T->tm_beg)) goto done;
alpar@9 234 /* build the objective, which is the distance between the current
alpar@9 235 (basic) point and the rounded one */
alpar@9 236 lp->dir = GLP_MIN;
alpar@9 237 lp->c0 = 0.0;
alpar@9 238 for (j = 1; j <= n; j++)
alpar@9 239 lp->col[j]->coef = 0.0;
alpar@9 240 for (k = 1; k <= nv; k++)
alpar@9 241 { j = var[k].j;
alpar@9 242 if (var[k].x == 0)
alpar@9 243 lp->col[j]->coef = +1.0;
alpar@9 244 else
alpar@9 245 { lp->col[j]->coef = -1.0;
alpar@9 246 lp->c0 += 1.0;
alpar@9 247 }
alpar@9 248 }
alpar@9 249 /* minimize the distance with the simplex method */
alpar@9 250 glp_init_smcp(&parm);
alpar@9 251 if (T->parm->msg_lev <= GLP_MSG_ERR)
alpar@9 252 parm.msg_lev = T->parm->msg_lev;
alpar@9 253 else if (T->parm->msg_lev <= GLP_MSG_ALL)
alpar@9 254 { parm.msg_lev = GLP_MSG_ON;
alpar@9 255 parm.out_dly = 10000;
alpar@9 256 }
alpar@9 257 ret = glp_simplex(lp, &parm);
alpar@9 258 if (ret != 0)
alpar@9 259 { if (T->parm->msg_lev >= GLP_MSG_ERR)
alpar@9 260 xprintf("Warning: glp_simplex returned %d\n", ret);
alpar@9 261 goto done;
alpar@9 262 }
alpar@9 263 ret = glp_get_status(lp);
alpar@9 264 if (ret != GLP_OPT)
alpar@9 265 { if (T->parm->msg_lev >= GLP_MSG_ERR)
alpar@9 266 xprintf("Warning: glp_get_status returned %d\n", ret);
alpar@9 267 goto done;
alpar@9 268 }
alpar@9 269 if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@9 270 xprintf("delta = %g\n", lp->obj_val);
alpar@9 271 /* check if the basic solution is integer feasible; note that it
alpar@9 272 may be so even if the minimial distance is positive */
alpar@9 273 tol = 0.3 * T->parm->tol_int;
alpar@9 274 for (k = 1; k <= nv; k++)
alpar@9 275 { col = lp->col[var[k].j];
alpar@9 276 if (tol < col->prim && col->prim < 1.0 - tol) break;
alpar@9 277 }
alpar@9 278 if (k > nv)
alpar@9 279 { /* okay; the basic solution seems to be integer feasible */
alpar@9 280 double *x = xcalloc(1+n, sizeof(double));
alpar@9 281 for (j = 1; j <= n; j++)
alpar@9 282 { x[j] = lp->col[j]->prim;
alpar@9 283 if (P->col[j]->kind == GLP_IV) x[j] = floor(x[j] + 0.5);
alpar@9 284 }
alpar@9 285 #if 1 /* modified by xypron <xypron.glpk@gmx.de> */
alpar@9 286 /* reset direction and right-hand side of objective */
alpar@9 287 lp->c0 = P->c0;
alpar@9 288 lp->dir = P->dir;
alpar@9 289 /* fix integer variables */
alpar@9 290 for (k = 1; k <= nv; k++)
alpar@9 291 { lp->col[var[k].j]->lb = x[var[k].j];
alpar@9 292 lp->col[var[k].j]->ub = x[var[k].j];
alpar@9 293 lp->col[var[k].j]->type = GLP_FX;
alpar@9 294 }
alpar@9 295 /* copy original objective function */
alpar@9 296 for (j = 1; j <= n; j++)
alpar@9 297 lp->col[j]->coef = P->col[j]->coef;
alpar@9 298 /* solve original LP and copy result */
alpar@9 299 ret = glp_simplex(lp, &parm);
alpar@9 300 if (ret != 0)
alpar@9 301 { if (T->parm->msg_lev >= GLP_MSG_ERR)
alpar@9 302 xprintf("Warning: glp_simplex returned %d\n", ret);
alpar@9 303 goto done;
alpar@9 304 }
alpar@9 305 ret = glp_get_status(lp);
alpar@9 306 if (ret != GLP_OPT)
alpar@9 307 { if (T->parm->msg_lev >= GLP_MSG_ERR)
alpar@9 308 xprintf("Warning: glp_get_status returned %d\n", ret);
alpar@9 309 goto done;
alpar@9 310 }
alpar@9 311 for (j = 1; j <= n; j++)
alpar@9 312 if (P->col[j]->kind != GLP_IV) x[j] = lp->col[j]->prim;
alpar@9 313 #endif
alpar@9 314 ret = glp_ios_heur_sol(T, x);
alpar@9 315 xfree(x);
alpar@9 316 if (ret == 0)
alpar@9 317 { /* the integer solution is accepted */
alpar@9 318 if (ios_is_hopeful(T, T->curr->bound))
alpar@9 319 { /* it is reasonable to apply the heuristic once again */
alpar@9 320 goto more;
alpar@9 321 }
alpar@9 322 else
alpar@9 323 { /* the best known integer feasible solution just found
alpar@9 324 is close to optimal solution to LP relaxation */
alpar@9 325 goto done;
alpar@9 326 }
alpar@9 327 }
alpar@9 328 }
alpar@9 329 /* the basic solution is fractional */
alpar@9 330 if (dist == DBL_MAX ||
alpar@9 331 lp->obj_val <= dist - 1e-6 * (1.0 + dist))
alpar@9 332 { /* the distance is reducing */
alpar@9 333 nfail = 0, dist = lp->obj_val;
alpar@9 334 }
alpar@9 335 else
alpar@9 336 { /* improving the distance failed */
alpar@9 337 nfail++;
alpar@9 338 }
alpar@9 339 if (nfail < 3) goto loop;
alpar@9 340 if (npass < 5) goto pass;
alpar@9 341 done: /* delete working objects */
alpar@9 342 if (lp != NULL) glp_delete_prob(lp);
alpar@9 343 if (var != NULL) xfree(var);
alpar@9 344 if (rand != NULL) rng_delete_rand(rand);
alpar@9 345 return;
alpar@9 346 }
alpar@9 347
alpar@9 348 /* eof */