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