lemon-project-template-glpk

comparison deps/glpk/src/glpios10.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:4d1994fef395
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, 2011 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
48 struct 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
58 static 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
69 void 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();
119 more: /* 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;
170 pass: /* 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 }
196 loop: /* 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 }
230 skip: /* 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;
341 done: /* 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 */