1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpios10.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,348 @@
     1.4 +/* glpios10.c (feasibility pump heuristic) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpios.h"
    1.29 +#include "glprng.h"
    1.30 +
    1.31 +/***********************************************************************
    1.32 +*  NAME
    1.33 +*
    1.34 +*  ios_feas_pump - feasibility pump heuristic
    1.35 +*
    1.36 +*  SYNOPSIS
    1.37 +*
    1.38 +*  #include "glpios.h"
    1.39 +*  void ios_feas_pump(glp_tree *T);
    1.40 +*
    1.41 +*  DESCRIPTION
    1.42 +*
    1.43 +*  The routine ios_feas_pump is a simple implementation of the Feasi-
    1.44 +*  bility Pump heuristic.
    1.45 +*
    1.46 +*  REFERENCES
    1.47 +*
    1.48 +*  M.Fischetti, F.Glover, and A.Lodi. "The feasibility pump." Math.
    1.49 +*  Program., Ser. A 104, pp. 91-104 (2005). */
    1.50 +
    1.51 +struct VAR
    1.52 +{     /* binary variable */
    1.53 +      int j;
    1.54 +      /* ordinal number */
    1.55 +      int x;
    1.56 +      /* value in the rounded solution (0 or 1) */
    1.57 +      double d;
    1.58 +      /* sorting key */
    1.59 +};
    1.60 +
    1.61 +static int fcmp(const void *x, const void *y)
    1.62 +{     /* comparison routine */
    1.63 +      const struct VAR *vx = x, *vy = y;
    1.64 +      if (vx->d > vy->d)
    1.65 +         return -1;
    1.66 +      else if (vx->d < vy->d)
    1.67 +         return +1;
    1.68 +      else
    1.69 +         return 0;
    1.70 +}
    1.71 +
    1.72 +void ios_feas_pump(glp_tree *T)
    1.73 +{     glp_prob *P = T->mip;
    1.74 +      int n = P->n;
    1.75 +      glp_prob *lp = NULL;
    1.76 +      struct VAR *var = NULL;
    1.77 +      RNG *rand = NULL;
    1.78 +      GLPCOL *col;
    1.79 +      glp_smcp parm;
    1.80 +      int j, k, new_x, nfail, npass, nv, ret, stalling;
    1.81 +      double dist, tol;
    1.82 +      xassert(glp_get_status(P) == GLP_OPT);
    1.83 +      /* this heuristic is applied only once on the root level */
    1.84 +      if (!(T->curr->level == 0 && T->curr->solved == 1)) goto done;
    1.85 +      /* determine number of binary variables */
    1.86 +      nv = 0;
    1.87 +      for (j = 1; j <= n; j++)
    1.88 +      {  col = P->col[j];
    1.89 +         /* if x[j] is continuous, skip it */
    1.90 +         if (col->kind == GLP_CV) continue;
    1.91 +         /* if x[j] is fixed, skip it */
    1.92 +         if (col->type == GLP_FX) continue;
    1.93 +         /* x[j] is non-fixed integer */
    1.94 +         xassert(col->kind == GLP_IV);
    1.95 +         if (col->type == GLP_DB && col->lb == 0.0 && col->ub == 1.0)
    1.96 +         {  /* x[j] is binary */
    1.97 +            nv++;
    1.98 +         }
    1.99 +         else
   1.100 +         {  /* x[j] is general integer */
   1.101 +            if (T->parm->msg_lev >= GLP_MSG_ALL)
   1.102 +               xprintf("FPUMP heuristic cannot be applied due to genera"
   1.103 +                  "l integer variables\n");
   1.104 +            goto done;
   1.105 +         }
   1.106 +      }
   1.107 +      /* there must be at least one binary variable */
   1.108 +      if (nv == 0) goto done;
   1.109 +      if (T->parm->msg_lev >= GLP_MSG_ALL)
   1.110 +         xprintf("Applying FPUMP heuristic...\n");
   1.111 +      /* build the list of binary variables */
   1.112 +      var = xcalloc(1+nv, sizeof(struct VAR));
   1.113 +      k = 0;
   1.114 +      for (j = 1; j <= n; j++)
   1.115 +      {  col = P->col[j];
   1.116 +         if (col->kind == GLP_IV && col->type == GLP_DB)
   1.117 +            var[++k].j = j;
   1.118 +      }
   1.119 +      xassert(k == nv);
   1.120 +      /* create working problem object */
   1.121 +      lp = glp_create_prob();
   1.122 +more: /* copy the original problem object to keep it intact */
   1.123 +      glp_copy_prob(lp, P, GLP_OFF);
   1.124 +      /* we are interested to find an integer feasible solution, which
   1.125 +         is better than the best known one */
   1.126 +      if (P->mip_stat == GLP_FEAS)
   1.127 +      {  int *ind;
   1.128 +         double *val, bnd;
   1.129 +         /* add a row and make it identical to the objective row */
   1.130 +         glp_add_rows(lp, 1);
   1.131 +         ind = xcalloc(1+n, sizeof(int));
   1.132 +         val = xcalloc(1+n, sizeof(double));
   1.133 +         for (j = 1; j <= n; j++)
   1.134 +         {  ind[j] = j;
   1.135 +            val[j] = P->col[j]->coef;
   1.136 +         }
   1.137 +         glp_set_mat_row(lp, lp->m, n, ind, val);
   1.138 +         xfree(ind);
   1.139 +         xfree(val);
   1.140 +         /* introduce upper (minimization) or lower (maximization)
   1.141 +            bound to the original objective function; note that this
   1.142 +            additional constraint is not violated at the optimal point
   1.143 +            to LP relaxation */
   1.144 +#if 0 /* modified by xypron <xypron.glpk@gmx.de> */
   1.145 +         if (P->dir == GLP_MIN)
   1.146 +         {  bnd = P->mip_obj - 0.10 * (1.0 + fabs(P->mip_obj));
   1.147 +            if (bnd < P->obj_val) bnd = P->obj_val;
   1.148 +            glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0);
   1.149 +         }
   1.150 +         else if (P->dir == GLP_MAX)
   1.151 +         {  bnd = P->mip_obj + 0.10 * (1.0 + fabs(P->mip_obj));
   1.152 +            if (bnd > P->obj_val) bnd = P->obj_val;
   1.153 +            glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0);
   1.154 +         }
   1.155 +         else
   1.156 +            xassert(P != P);
   1.157 +#else
   1.158 +         bnd = 0.1 * P->obj_val + 0.9 * P->mip_obj;
   1.159 +         /* xprintf("bnd = %f\n", bnd); */
   1.160 +         if (P->dir == GLP_MIN)
   1.161 +            glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0);
   1.162 +         else if (P->dir == GLP_MAX)
   1.163 +            glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0);
   1.164 +         else
   1.165 +            xassert(P != P);
   1.166 +#endif
   1.167 +      }
   1.168 +      /* reset pass count */
   1.169 +      npass = 0;
   1.170 +      /* invalidate the rounded point */
   1.171 +      for (k = 1; k <= nv; k++)
   1.172 +         var[k].x = -1;
   1.173 +pass: /* next pass starts here */
   1.174 +      npass++;
   1.175 +      if (T->parm->msg_lev >= GLP_MSG_ALL)
   1.176 +         xprintf("Pass %d\n", npass);
   1.177 +      /* initialize minimal distance between the basic point and the
   1.178 +         rounded one obtained during this pass */
   1.179 +      dist = DBL_MAX;
   1.180 +      /* reset failure count (the number of succeeded iterations failed
   1.181 +         to improve the distance) */
   1.182 +      nfail = 0;
   1.183 +      /* if it is not the first pass, perturb the last rounded point
   1.184 +         rather than construct it from the basic solution */
   1.185 +      if (npass > 1)
   1.186 +      {  double rho, temp;
   1.187 +         if (rand == NULL)
   1.188 +            rand = rng_create_rand();
   1.189 +         for (k = 1; k <= nv; k++)
   1.190 +         {  j = var[k].j;
   1.191 +            col = lp->col[j];
   1.192 +            rho = rng_uniform(rand, -0.3, 0.7);
   1.193 +            if (rho < 0.0) rho = 0.0;
   1.194 +            temp = fabs((double)var[k].x - col->prim);
   1.195 +            if (temp + rho > 0.5) var[k].x = 1 - var[k].x;
   1.196 +         }
   1.197 +         goto skip;
   1.198 +      }
   1.199 +loop: /* innermost loop begins here */
   1.200 +      /* round basic solution (which is assumed primal feasible) */
   1.201 +      stalling = 1;
   1.202 +      for (k = 1; k <= nv; k++)
   1.203 +      {  col = lp->col[var[k].j];
   1.204 +         if (col->prim < 0.5)
   1.205 +         {  /* rounded value is 0 */
   1.206 +            new_x = 0;
   1.207 +         }
   1.208 +         else
   1.209 +         {  /* rounded value is 1 */
   1.210 +            new_x = 1;
   1.211 +         }
   1.212 +         if (var[k].x != new_x)
   1.213 +         {  stalling = 0;
   1.214 +            var[k].x = new_x;
   1.215 +         }
   1.216 +      }
   1.217 +      /* if the rounded point has not changed (stalling), choose and
   1.218 +         flip some its entries heuristically */
   1.219 +      if (stalling)
   1.220 +      {  /* compute d[j] = |x[j] - round(x[j])| */
   1.221 +         for (k = 1; k <= nv; k++)
   1.222 +         {  col = lp->col[var[k].j];
   1.223 +            var[k].d = fabs(col->prim - (double)var[k].x);
   1.224 +         }
   1.225 +         /* sort the list of binary variables by descending d[j] */
   1.226 +         qsort(&var[1], nv, sizeof(struct VAR), fcmp);
   1.227 +         /* choose and flip some rounded components */
   1.228 +         for (k = 1; k <= nv; k++)
   1.229 +         {  if (k >= 5 && var[k].d < 0.35 || k >= 10) break;
   1.230 +            var[k].x = 1 - var[k].x;
   1.231 +         }
   1.232 +      }
   1.233 +skip: /* check if the time limit has been exhausted */
   1.234 +      if (T->parm->tm_lim < INT_MAX &&
   1.235 +         (double)(T->parm->tm_lim - 1) <=
   1.236 +         1000.0 * xdifftime(xtime(), T->tm_beg)) goto done;
   1.237 +      /* build the objective, which is the distance between the current
   1.238 +         (basic) point and the rounded one */
   1.239 +      lp->dir = GLP_MIN;
   1.240 +      lp->c0 = 0.0;
   1.241 +      for (j = 1; j <= n; j++)
   1.242 +         lp->col[j]->coef = 0.0;
   1.243 +      for (k = 1; k <= nv; k++)
   1.244 +      {  j = var[k].j;
   1.245 +         if (var[k].x == 0)
   1.246 +            lp->col[j]->coef = +1.0;
   1.247 +         else
   1.248 +         {  lp->col[j]->coef = -1.0;
   1.249 +            lp->c0 += 1.0;
   1.250 +         }
   1.251 +      }
   1.252 +      /* minimize the distance with the simplex method */
   1.253 +      glp_init_smcp(&parm);
   1.254 +      if (T->parm->msg_lev <= GLP_MSG_ERR)
   1.255 +         parm.msg_lev = T->parm->msg_lev;
   1.256 +      else if (T->parm->msg_lev <= GLP_MSG_ALL)
   1.257 +      {  parm.msg_lev = GLP_MSG_ON;
   1.258 +         parm.out_dly = 10000;
   1.259 +      }
   1.260 +      ret = glp_simplex(lp, &parm);
   1.261 +      if (ret != 0)
   1.262 +      {  if (T->parm->msg_lev >= GLP_MSG_ERR)
   1.263 +            xprintf("Warning: glp_simplex returned %d\n", ret);
   1.264 +         goto done;
   1.265 +      }
   1.266 +      ret = glp_get_status(lp);
   1.267 +      if (ret != GLP_OPT)
   1.268 +      {  if (T->parm->msg_lev >= GLP_MSG_ERR)
   1.269 +            xprintf("Warning: glp_get_status returned %d\n", ret);
   1.270 +         goto done;
   1.271 +      }
   1.272 +      if (T->parm->msg_lev >= GLP_MSG_DBG)
   1.273 +         xprintf("delta = %g\n", lp->obj_val);
   1.274 +      /* check if the basic solution is integer feasible; note that it
   1.275 +         may be so even if the minimial distance is positive */
   1.276 +      tol = 0.3 * T->parm->tol_int;
   1.277 +      for (k = 1; k <= nv; k++)
   1.278 +      {  col = lp->col[var[k].j];
   1.279 +         if (tol < col->prim && col->prim < 1.0 - tol) break;
   1.280 +      }
   1.281 +      if (k > nv)
   1.282 +      {  /* okay; the basic solution seems to be integer feasible */
   1.283 +         double *x = xcalloc(1+n, sizeof(double));
   1.284 +         for (j = 1; j <= n; j++)
   1.285 +         {  x[j] = lp->col[j]->prim;
   1.286 +            if (P->col[j]->kind == GLP_IV) x[j] = floor(x[j] + 0.5);
   1.287 +         }
   1.288 +#if 1 /* modified by xypron <xypron.glpk@gmx.de> */
   1.289 +         /* reset direction and right-hand side of objective */
   1.290 +         lp->c0  = P->c0;
   1.291 +         lp->dir = P->dir;
   1.292 +         /* fix integer variables */
   1.293 +         for (k = 1; k <= nv; k++)
   1.294 +         {  lp->col[var[k].j]->lb   = x[var[k].j];
   1.295 +            lp->col[var[k].j]->ub   = x[var[k].j];
   1.296 +            lp->col[var[k].j]->type = GLP_FX;
   1.297 +         }
   1.298 +         /* copy original objective function */
   1.299 +         for (j = 1; j <= n; j++)
   1.300 +            lp->col[j]->coef = P->col[j]->coef;
   1.301 +         /* solve original LP and copy result */
   1.302 +         ret = glp_simplex(lp, &parm);
   1.303 +         if (ret != 0)
   1.304 +         {  if (T->parm->msg_lev >= GLP_MSG_ERR)
   1.305 +               xprintf("Warning: glp_simplex returned %d\n", ret);
   1.306 +            goto done;
   1.307 +         }
   1.308 +         ret = glp_get_status(lp);
   1.309 +         if (ret != GLP_OPT)
   1.310 +         {  if (T->parm->msg_lev >= GLP_MSG_ERR)
   1.311 +               xprintf("Warning: glp_get_status returned %d\n", ret);
   1.312 +            goto done;
   1.313 +         }
   1.314 +         for (j = 1; j <= n; j++)
   1.315 +            if (P->col[j]->kind != GLP_IV) x[j] = lp->col[j]->prim;
   1.316 +#endif
   1.317 +         ret = glp_ios_heur_sol(T, x);
   1.318 +         xfree(x);
   1.319 +         if (ret == 0)
   1.320 +         {  /* the integer solution is accepted */
   1.321 +            if (ios_is_hopeful(T, T->curr->bound))
   1.322 +            {  /* it is reasonable to apply the heuristic once again */
   1.323 +               goto more;
   1.324 +            }
   1.325 +            else
   1.326 +            {  /* the best known integer feasible solution just found
   1.327 +                  is close to optimal solution to LP relaxation */
   1.328 +               goto done;
   1.329 +            }
   1.330 +         }
   1.331 +      }
   1.332 +      /* the basic solution is fractional */
   1.333 +      if (dist == DBL_MAX ||
   1.334 +          lp->obj_val <= dist - 1e-6 * (1.0 + dist))
   1.335 +      {  /* the distance is reducing */
   1.336 +         nfail = 0, dist = lp->obj_val;
   1.337 +      }
   1.338 +      else
   1.339 +      {  /* improving the distance failed */
   1.340 +         nfail++;
   1.341 +      }
   1.342 +      if (nfail < 3) goto loop;
   1.343 +      if (npass < 5) goto pass;
   1.344 +done: /* delete working objects */
   1.345 +      if (lp != NULL) glp_delete_prob(lp);
   1.346 +      if (var != NULL) xfree(var);
   1.347 +      if (rand != NULL) rng_delete_rand(rand);
   1.348 +      return;
   1.349 +}
   1.350 +
   1.351 +/* eof */