lemon-project-template-glpk: 33de93886c88 deps/glpk/src/glpios10.c

lemon-project-template-glpk

view deps/glpk/src/glpios10.c @ 9:33de93886c88

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

line source

1 /* glpios10.c (feasibility pump heuristic) */

3 /***********************************************************************

4 * This code is part of GLPK (GNU Linear Programming Kit).

5 *

7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,

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

25 #include "glpios.h"

26 #include "glprng.h"

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

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 };

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 }

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