/* fpump.c (feasibility pump heuristic) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*  Copyright (C) 2009-2018 Free Software Foundation, Inc.
*  Written by Andrew Makhorin <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 "env.h"
#include "ios.h"
#include "rng.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 CDECL 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++)
#if 0 /* 18/VI-2013; fixed by mao
       * this bug causes numerical instability, because column statuses
       * are not changed appropriately */
         {  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;
         }
#else
            glp_set_col_bnds(lp, var[k].j, GLP_FX, x[var[k].j], 0.);
#endif
         /* 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);
#if 1 /* 17/III-2016: fix memory leak */
            xfree(x);
#endif
            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);
#if 1 /* 17/III-2016: fix memory leak */
            xfree(x);
#endif
            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 */