/* 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 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 */