lemon-project-template-glpk: deps/glpk/src/glpios10.c comparison

lemon-project-template-glpk

comparison deps/glpk/src/glpios10.c @ 9:33de93886c88

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:4d1994fef395
+/* glpios10.c (feasibility pump heuristic) */
+/***********************************************************************
+*  This code is part of GLPK (GNU Linear Programming Kit).
+*
+*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
+*  2009, 2010, 2011 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 "glpios.h"
+#include "glprng.h"
+/***********************************************************************
+*  NAME
+*
+*  ios_feas_pump - feasibility pump heuristic
+*
+*  SYNOPSIS
+*
+*  #include "glpios.h"
+*  void ios_feas_pump(glp_tree *T);
+*
+*  DESCRIPTION
+*
+*  The routine ios_feas_pump is a simple implementation of the Feasi-
+*  bility Pump heuristic.
+*
+*  REFERENCES
+*
+*  M.Fischetti, F.Glover, and A.Lodi. "The feasibility pump." Math.
+*  Program., Ser. A 104, pp. 91-104 (2005). */
+struct VAR
+{     /* binary variable */
+int j;
+/* ordinal number */
+int x;
+/* value in the rounded solution (0 or 1) */
+double d;
+/* sorting key */
+};
+static int fcmp(const void *x, const void *y)
+{     /* comparison routine */
+const struct VAR *vx = x, *vy = y;
+if (vx->d > vy->d)
+return -1;
+else if (vx->d < vy->d)
+return +1;
+else
+return 0;
+}
+void ios_feas_pump(glp_tree *T)
+{     glp_prob *P = T->mip;
+int n = P->n;
+glp_prob *lp = NULL;
+struct VAR *var = NULL;
+RNG *rand = NULL;
+GLPCOL *col;
+glp_smcp parm;
+int j, k, new_x, nfail, npass, nv, ret, stalling;
+double dist, tol;
+xassert(glp_get_status(P) == GLP_OPT);
+/* this heuristic is applied only once on the root level */
+if (!(T->curr->level == 0 && T->curr->solved == 1)) goto done;
+/* determine number of binary variables */
+nv = 0;
+for (j = 1; j <= n; j++)
+{  col = P->col[j];
+/* if x[j] is continuous, skip it */
+if (col->kind == GLP_CV) continue;
+/* if x[j] is fixed, skip it */
+if (col->type == GLP_FX) continue;
+/* x[j] is non-fixed integer */
+xassert(col->kind == GLP_IV);
+if (col->type == GLP_DB && col->lb == 0.0 && col->ub == 1.0)
+{  /* x[j] is binary */
+nv++;
+}
+else
+{  /* x[j] is general integer */
+if (T->parm->msg_lev >= GLP_MSG_ALL)
+xprintf("FPUMP heuristic cannot be applied due to genera"
+"l integer variables\n");
+goto done;
+}
+}
+/* there must be at least one binary variable */
+if (nv == 0) goto done;
+if (T->parm->msg_lev >= GLP_MSG_ALL)
+xprintf("Applying FPUMP heuristic...\n");
+/* build the list of binary variables */
+var = xcalloc(1+nv, sizeof(struct VAR));
+k = 0;
+for (j = 1; j <= n; j++)
+{  col = P->col[j];
+if (col->kind == GLP_IV && col->type == GLP_DB)
+var[++k].j = j;
+}
+xassert(k == nv);
+/* create working problem object */
+lp = glp_create_prob();
+more: /* copy the original problem object to keep it intact */
+glp_copy_prob(lp, P, GLP_OFF);
+/* we are interested to find an integer feasible solution, which
+is better than the best known one */
+if (P->mip_stat == GLP_FEAS)
+{  int *ind;
+double *val, bnd;
+/* add a row and make it identical to the objective row */
+glp_add_rows(lp, 1);
+ind = xcalloc(1+n, sizeof(int));
+val = xcalloc(1+n, sizeof(double));
+for (j = 1; j <= n; j++)
+{  ind[j] = j;
+val[j] = P->col[j]->coef;
+}
+glp_set_mat_row(lp, lp->m, n, ind, val);
+xfree(ind);
+xfree(val);
+/* introduce upper (minimization) or lower (maximization)
+bound to the original objective function; note that this
+additional constraint is not violated at the optimal point
+to LP relaxation */
+#if 0 /* modified by xypron <xypron.glpk@gmx.de> */
+if (P->dir == GLP_MIN)
+{  bnd = P->mip_obj - 0.10 * (1.0 + fabs(P->mip_obj));
+if (bnd < P->obj_val) bnd = P->obj_val;
+glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0);
+}
+else if (P->dir == GLP_MAX)
+{  bnd = P->mip_obj + 0.10 * (1.0 + fabs(P->mip_obj));
+if (bnd > P->obj_val) bnd = P->obj_val;
+glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0);
+}
+else
+xassert(P != P);
+#else
+bnd = 0.1 * P->obj_val + 0.9 * P->mip_obj;
+/* xprintf("bnd = %f\n", bnd); */
+if (P->dir == GLP_MIN)
+glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0);
+else if (P->dir == GLP_MAX)
+glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0);
+else
+xassert(P != P);
+#endif
+}
+/* reset pass count */
+npass = 0;
+/* invalidate the rounded point */
+for (k = 1; k <= nv; k++)
+var[k].x = -1;
+pass: /* next pass starts here */
+npass++;
+if (T->parm->msg_lev >= GLP_MSG_ALL)
+xprintf("Pass %d\n", npass);
+/* initialize minimal distance between the basic point and the
+rounded one obtained during this pass */
+dist = DBL_MAX;
+/* reset failure count (the number of succeeded iterations failed
+to improve the distance) */
+nfail = 0;
+/* if it is not the first pass, perturb the last rounded point
+rather than construct it from the basic solution */
+if (npass > 1)
+{  double rho, temp;
+if (rand == NULL)
+rand = rng_create_rand();
+for (k = 1; k <= nv; k++)
+{  j = var[k].j;
+col = lp->col[j];
+rho = rng_uniform(rand, -0.3, 0.7);
+if (rho < 0.0) rho = 0.0;
+temp = fabs((double)var[k].x - col->prim);
+if (temp + rho > 0.5) var[k].x = 1 - var[k].x;
+}
+goto skip;
+}
+loop: /* innermost loop begins here */
+/* round basic solution (which is assumed primal feasible) */
+stalling = 1;
+for (k = 1; k <= nv; k++)
+{  col = lp->col[var[k].j];
+if (col->prim < 0.5)
+{  /* rounded value is 0 */
+new_x = 0;
+}
+else
+{  /* rounded value is 1 */
+new_x = 1;
+}
+if (var[k].x != new_x)
+{  stalling = 0;
+var[k].x = new_x;
+}
+}
+/* if the rounded point has not changed (stalling), choose and
+flip some its entries heuristically */
+if (stalling)
+{  /* compute d[j] = |x[j] - round(x[j])| */
+for (k = 1; k <= nv; k++)
+{  col = lp->col[var[k].j];
+var[k].d = fabs(col->prim - (double)var[k].x);
+}
+/* sort the list of binary variables by descending d[j] */
+qsort(&var[1], nv, sizeof(struct VAR), fcmp);
+/* choose and flip some rounded components */
+for (k = 1; k <= nv; k++)
+{  if (k >= 5 && var[k].d < 0.35 || k >= 10) break;
+var[k].x = 1 - var[k].x;
+}
+}
+skip: /* check if the time limit has been exhausted */
+if (T->parm->tm_lim < INT_MAX &&
+(double)(T->parm->tm_lim - 1) <=
+1000.0 * xdifftime(xtime(), T->tm_beg)) goto done;
+/* build the objective, which is the distance between the current
+(basic) point and the rounded one */
+lp->dir = GLP_MIN;
+lp->c0 = 0.0;
+for (j = 1; j <= n; j++)
+lp->col[j]->coef = 0.0;
+for (k = 1; k <= nv; k++)
+{  j = var[k].j;
+if (var[k].x == 0)
+lp->col[j]->coef = +1.0;
+else
+{  lp->col[j]->coef = -1.0;
+lp->c0 += 1.0;
+}
+}
+/* minimize the distance with the simplex method */
+glp_init_smcp(&parm);
+if (T->parm->msg_lev <= GLP_MSG_ERR)
+parm.msg_lev = T->parm->msg_lev;
+else if (T->parm->msg_lev <= GLP_MSG_ALL)
+{  parm.msg_lev = GLP_MSG_ON;
+parm.out_dly = 10000;
+}
+ret = glp_simplex(lp, &parm);
+if (ret != 0)
+{  if (T->parm->msg_lev >= GLP_MSG_ERR)
+xprintf("Warning: glp_simplex returned %d\n", ret);
+goto done;
+}
+ret = glp_get_status(lp);
+if (ret != GLP_OPT)
+{  if (T->parm->msg_lev >= GLP_MSG_ERR)
+xprintf("Warning: glp_get_status returned %d\n", ret);
+goto done;
+}
+if (T->parm->msg_lev >= GLP_MSG_DBG)
+xprintf("delta = %g\n", lp->obj_val);
+/* check if the basic solution is integer feasible; note that it
+may be so even if the minimial distance is positive */
+tol = 0.3 * T->parm->tol_int;
+for (k = 1; k <= nv; k++)
+{  col = lp->col[var[k].j];
+if (tol < col->prim && col->prim < 1.0 - tol) break;
+}
+if (k > nv)
+{  /* okay; the basic solution seems to be integer feasible */
+double *x = xcalloc(1+n, sizeof(double));
+for (j = 1; j <= n; j++)
+{  x[j] = lp->col[j]->prim;
+if (P->col[j]->kind == GLP_IV) x[j] = floor(x[j] + 0.5);
+}
+#if 1 /* modified by xypron <xypron.glpk@gmx.de> */
+/* reset direction and right-hand side of objective */
+lp->c0  = P->c0;
+lp->dir = P->dir;
+/* fix integer variables */
+for (k = 1; k <= nv; k++)
+{  lp->col[var[k].j]->lb   = x[var[k].j];
+lp->col[var[k].j]->ub   = x[var[k].j];
+lp->col[var[k].j]->type = GLP_FX;
+}
+/* copy original objective function */
+for (j = 1; j <= n; j++)
+lp->col[j]->coef = P->col[j]->coef;
+/* solve original LP and copy result */
+ret = glp_simplex(lp, &parm);
+if (ret != 0)
+{  if (T->parm->msg_lev >= GLP_MSG_ERR)
+xprintf("Warning: glp_simplex returned %d\n", ret);
+goto done;
+}
+ret = glp_get_status(lp);
+if (ret != GLP_OPT)
+{  if (T->parm->msg_lev >= GLP_MSG_ERR)
+xprintf("Warning: glp_get_status returned %d\n", ret);
+goto done;
+}
+for (j = 1; j <= n; j++)
+if (P->col[j]->kind != GLP_IV) x[j] = lp->col[j]->prim;
+#endif
+ret = glp_ios_heur_sol(T, x);
+xfree(x);
+if (ret == 0)
+{  /* the integer solution is accepted */
+if (ios_is_hopeful(T, T->curr->bound))
+{  /* it is reasonable to apply the heuristic once again */
+goto more;
+}
+else
+{  /* the best known integer feasible solution just found
+is close to optimal solution to LP relaxation */
+goto done;
+}
+}
+}
+/* the basic solution is fractional */
+if (dist == DBL_MAX ||
+lp->obj_val <= dist - 1e-6 * (1.0 + dist))
+{  /* the distance is reducing */
+nfail = 0, dist = lp->obj_val;
+}
+else
+{  /* improving the distance failed */
+nfail++;
+}
+if (nfail < 3) goto loop;
+if (npass < 5) goto pass;
+done: /* delete working objects */
+if (lp != NULL) glp_delete_prob(lp);
+if (var != NULL) xfree(var);
+if (rand != NULL) rng_delete_rand(rand);
+return;
+}
+/* eof */