COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpios10.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 11.8 KB
Line 
1/* glpios10.c (feasibility pump heuristic) */
2
3/***********************************************************************
4*  This code is part of GLPK (GNU Linear Programming Kit).
5*
6*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
8*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9*  E-mail: <mao@gnu.org>.
10*
11*  GLPK is free software: you can redistribute it and/or modify it
12*  under the terms of the GNU General Public License as published by
13*  the Free Software Foundation, either version 3 of the License, or
14*  (at your option) any later version.
15*
16*  GLPK is distributed in the hope that it will be useful, but WITHOUT
17*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19*  License for more details.
20*
21*  You should have received a copy of the GNU General Public License
22*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23***********************************************************************/
24
25#include "glpios.h"
26#include "glprng.h"
27
28/***********************************************************************
29*  NAME
30*
31*  ios_feas_pump - feasibility pump heuristic
32*
33*  SYNOPSIS
34*
35*  #include "glpios.h"
36*  void ios_feas_pump(glp_tree *T);
37*
38*  DESCRIPTION
39*
40*  The routine ios_feas_pump is a simple implementation of the Feasi-
41*  bility Pump heuristic.
42*
43*  REFERENCES
44*
45*  M.Fischetti, F.Glover, and A.Lodi. "The feasibility pump." Math.
46*  Program., Ser. A 104, pp. 91-104 (2005). */
47
48struct VAR
49{     /* binary variable */
50      int j;
51      /* ordinal number */
52      int x;
53      /* value in the rounded solution (0 or 1) */
54      double d;
55      /* sorting key */
56};
57
58static int fcmp(const void *x, const void *y)
59{     /* comparison routine */
60      const struct VAR *vx = x, *vy = y;
61      if (vx->d > vy->d)
62         return -1;
63      else if (vx->d < vy->d)
64         return +1;
65      else
66         return 0;
67}
68
69void ios_feas_pump(glp_tree *T)
70{     glp_prob *P = T->mip;
71      int n = P->n;
72      glp_prob *lp = NULL;
73      struct VAR *var = NULL;
74      RNG *rand = NULL;
75      GLPCOL *col;
76      glp_smcp parm;
77      int j, k, new_x, nfail, npass, nv, ret, stalling;
78      double dist, tol;
79      xassert(glp_get_status(P) == GLP_OPT);
80      /* this heuristic is applied only once on the root level */
81      if (!(T->curr->level == 0 && T->curr->solved == 1)) goto done;
82      /* determine number of binary variables */
83      nv = 0;
84      for (j = 1; j <= n; j++)
85      {  col = P->col[j];
86         /* if x[j] is continuous, skip it */
87         if (col->kind == GLP_CV) continue;
88         /* if x[j] is fixed, skip it */
89         if (col->type == GLP_FX) continue;
90         /* x[j] is non-fixed integer */
91         xassert(col->kind == GLP_IV);
92         if (col->type == GLP_DB && col->lb == 0.0 && col->ub == 1.0)
93         {  /* x[j] is binary */
94            nv++;
95         }
96         else
97         {  /* x[j] is general integer */
98            if (T->parm->msg_lev >= GLP_MSG_ALL)
99               xprintf("FPUMP heuristic cannot be applied due to genera"
100                  "l integer variables\n");
101            goto done;
102         }
103      }
104      /* there must be at least one binary variable */
105      if (nv == 0) goto done;
106      if (T->parm->msg_lev >= GLP_MSG_ALL)
107         xprintf("Applying FPUMP heuristic...\n");
108      /* build the list of binary variables */
109      var = xcalloc(1+nv, sizeof(struct VAR));
110      k = 0;
111      for (j = 1; j <= n; j++)
112      {  col = P->col[j];
113         if (col->kind == GLP_IV && col->type == GLP_DB)
114            var[++k].j = j;
115      }
116      xassert(k == nv);
117      /* create working problem object */
118      lp = glp_create_prob();
119more: /* copy the original problem object to keep it intact */
120      glp_copy_prob(lp, P, GLP_OFF);
121      /* we are interested to find an integer feasible solution, which
122         is better than the best known one */
123      if (P->mip_stat == GLP_FEAS)
124      {  int *ind;
125         double *val, bnd;
126         /* add a row and make it identical to the objective row */
127         glp_add_rows(lp, 1);
128         ind = xcalloc(1+n, sizeof(int));
129         val = xcalloc(1+n, sizeof(double));
130         for (j = 1; j <= n; j++)
131         {  ind[j] = j;
132            val[j] = P->col[j]->coef;
133         }
134         glp_set_mat_row(lp, lp->m, n, ind, val);
135         xfree(ind);
136         xfree(val);
137         /* introduce upper (minimization) or lower (maximization)
138            bound to the original objective function; note that this
139            additional constraint is not violated at the optimal point
140            to LP relaxation */
141#if 0 /* modified by xypron <xypron.glpk@gmx.de> */
142         if (P->dir == GLP_MIN)
143         {  bnd = P->mip_obj - 0.10 * (1.0 + fabs(P->mip_obj));
144            if (bnd < P->obj_val) bnd = P->obj_val;
145            glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0);
146         }
147         else if (P->dir == GLP_MAX)
148         {  bnd = P->mip_obj + 0.10 * (1.0 + fabs(P->mip_obj));
149            if (bnd > P->obj_val) bnd = P->obj_val;
150            glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0);
151         }
152         else
153            xassert(P != P);
154#else
155         bnd = 0.1 * P->obj_val + 0.9 * P->mip_obj;
156         /* xprintf("bnd = %f\n", bnd); */
157         if (P->dir == GLP_MIN)
158            glp_set_row_bnds(lp, lp->m, GLP_UP, 0.0, bnd - P->c0);
159         else if (P->dir == GLP_MAX)
160            glp_set_row_bnds(lp, lp->m, GLP_LO, bnd - P->c0, 0.0);
161         else
162            xassert(P != P);
163#endif
164      }
165      /* reset pass count */
166      npass = 0;
167      /* invalidate the rounded point */
168      for (k = 1; k <= nv; k++)
169         var[k].x = -1;
170pass: /* next pass starts here */
171      npass++;
172      if (T->parm->msg_lev >= GLP_MSG_ALL)
173         xprintf("Pass %d\n", npass);
174      /* initialize minimal distance between the basic point and the
175         rounded one obtained during this pass */
176      dist = DBL_MAX;
177      /* reset failure count (the number of succeeded iterations failed
178         to improve the distance) */
179      nfail = 0;
180      /* if it is not the first pass, perturb the last rounded point
181         rather than construct it from the basic solution */
182      if (npass > 1)
183      {  double rho, temp;
184         if (rand == NULL)
185            rand = rng_create_rand();
186         for (k = 1; k <= nv; k++)
187         {  j = var[k].j;
188            col = lp->col[j];
189            rho = rng_uniform(rand, -0.3, 0.7);
190            if (rho < 0.0) rho = 0.0;
191            temp = fabs((double)var[k].x - col->prim);
192            if (temp + rho > 0.5) var[k].x = 1 - var[k].x;
193         }
194         goto skip;
195      }
196loop: /* innermost loop begins here */
197      /* round basic solution (which is assumed primal feasible) */
198      stalling = 1;
199      for (k = 1; k <= nv; k++)
200      {  col = lp->col[var[k].j];
201         if (col->prim < 0.5)
202         {  /* rounded value is 0 */
203            new_x = 0;
204         }
205         else
206         {  /* rounded value is 1 */
207            new_x = 1;
208         }
209         if (var[k].x != new_x)
210         {  stalling = 0;
211            var[k].x = new_x;
212         }
213      }
214      /* if the rounded point has not changed (stalling), choose and
215         flip some its entries heuristically */
216      if (stalling)
217      {  /* compute d[j] = |x[j] - round(x[j])| */
218         for (k = 1; k <= nv; k++)
219         {  col = lp->col[var[k].j];
220            var[k].d = fabs(col->prim - (double)var[k].x);
221         }
222         /* sort the list of binary variables by descending d[j] */
223         qsort(&var[1], nv, sizeof(struct VAR), fcmp);
224         /* choose and flip some rounded components */
225         for (k = 1; k <= nv; k++)
226         {  if (k >= 5 && var[k].d < 0.35 || k >= 10) break;
227            var[k].x = 1 - var[k].x;
228         }
229      }
230skip: /* check if the time limit has been exhausted */
231      if (T->parm->tm_lim < INT_MAX &&
232         (double)(T->parm->tm_lim - 1) <=
233         1000.0 * xdifftime(xtime(), T->tm_beg)) goto done;
234      /* build the objective, which is the distance between the current
235         (basic) point and the rounded one */
236      lp->dir = GLP_MIN;
237      lp->c0 = 0.0;
238      for (j = 1; j <= n; j++)
239         lp->col[j]->coef = 0.0;
240      for (k = 1; k <= nv; k++)
241      {  j = var[k].j;
242         if (var[k].x == 0)
243            lp->col[j]->coef = +1.0;
244         else
245         {  lp->col[j]->coef = -1.0;
246            lp->c0 += 1.0;
247         }
248      }
249      /* minimize the distance with the simplex method */
250      glp_init_smcp(&parm);
251      if (T->parm->msg_lev <= GLP_MSG_ERR)
252         parm.msg_lev = T->parm->msg_lev;
253      else if (T->parm->msg_lev <= GLP_MSG_ALL)
254      {  parm.msg_lev = GLP_MSG_ON;
255         parm.out_dly = 10000;
256      }
257      ret = glp_simplex(lp, &parm);
258      if (ret != 0)
259      {  if (T->parm->msg_lev >= GLP_MSG_ERR)
260            xprintf("Warning: glp_simplex returned %d\n", ret);
261         goto done;
262      }
263      ret = glp_get_status(lp);
264      if (ret != GLP_OPT)
265      {  if (T->parm->msg_lev >= GLP_MSG_ERR)
266            xprintf("Warning: glp_get_status returned %d\n", ret);
267         goto done;
268      }
269      if (T->parm->msg_lev >= GLP_MSG_DBG)
270         xprintf("delta = %g\n", lp->obj_val);
271      /* check if the basic solution is integer feasible; note that it
272         may be so even if the minimial distance is positive */
273      tol = 0.3 * T->parm->tol_int;
274      for (k = 1; k <= nv; k++)
275      {  col = lp->col[var[k].j];
276         if (tol < col->prim && col->prim < 1.0 - tol) break;
277      }
278      if (k > nv)
279      {  /* okay; the basic solution seems to be integer feasible */
280         double *x = xcalloc(1+n, sizeof(double));
281         for (j = 1; j <= n; j++)
282         {  x[j] = lp->col[j]->prim;
283            if (P->col[j]->kind == GLP_IV) x[j] = floor(x[j] + 0.5);
284         }
285#if 1 /* modified by xypron <xypron.glpk@gmx.de> */
286         /* reset direction and right-hand side of objective */
287         lp->c0  = P->c0;
288         lp->dir = P->dir;
289         /* fix integer variables */
290         for (k = 1; k <= nv; k++)
291         {  lp->col[var[k].j]->lb   = x[var[k].j];
292            lp->col[var[k].j]->ub   = x[var[k].j];
293            lp->col[var[k].j]->type = GLP_FX;
294         }
295         /* copy original objective function */
296         for (j = 1; j <= n; j++)
297            lp->col[j]->coef = P->col[j]->coef;
298         /* solve original LP and copy result */
299         ret = glp_simplex(lp, &parm);
300         if (ret != 0)
301         {  if (T->parm->msg_lev >= GLP_MSG_ERR)
302               xprintf("Warning: glp_simplex returned %d\n", ret);
303            goto done;
304         }
305         ret = glp_get_status(lp);
306         if (ret != GLP_OPT)
307         {  if (T->parm->msg_lev >= GLP_MSG_ERR)
308               xprintf("Warning: glp_get_status returned %d\n", ret);
309            goto done;
310         }
311         for (j = 1; j <= n; j++)
312            if (P->col[j]->kind != GLP_IV) x[j] = lp->col[j]->prim;
313#endif
314         ret = glp_ios_heur_sol(T, x);
315         xfree(x);
316         if (ret == 0)
317         {  /* the integer solution is accepted */
318            if (ios_is_hopeful(T, T->curr->bound))
319            {  /* it is reasonable to apply the heuristic once again */
320               goto more;
321            }
322            else
323            {  /* the best known integer feasible solution just found
324                  is close to optimal solution to LP relaxation */
325               goto done;
326            }
327         }
328      }
329      /* the basic solution is fractional */
330      if (dist == DBL_MAX ||
331          lp->obj_val <= dist - 1e-6 * (1.0 + dist))
332      {  /* the distance is reducing */
333         nfail = 0, dist = lp->obj_val;
334      }
335      else
336      {  /* improving the distance failed */
337         nfail++;
338      }
339      if (nfail < 3) goto loop;
340      if (npass < 5) goto pass;
341done: /* delete working objects */
342      if (lp != NULL) glp_delete_prob(lp);
343      if (var != NULL) xfree(var);
344      if (rand != NULL) rng_delete_rand(rand);
345      return;
346}
347
348/* eof */
Note: See TracBrowser for help on using the repository browser.