lemon-project-template-glpk

annotate deps/glpk/src/glpios09.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
rev   line source
alpar@9 1 /* glpios09.c (branching heuristics) */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpios.h"
alpar@9 26
alpar@9 27 /***********************************************************************
alpar@9 28 * NAME
alpar@9 29 *
alpar@9 30 * ios_choose_var - select variable to branch on
alpar@9 31 *
alpar@9 32 * SYNOPSIS
alpar@9 33 *
alpar@9 34 * #include "glpios.h"
alpar@9 35 * int ios_choose_var(glp_tree *T, int *next);
alpar@9 36 *
alpar@9 37 * The routine ios_choose_var chooses a variable from the candidate
alpar@9 38 * list to branch on. Additionally the routine provides a flag stored
alpar@9 39 * in the location next to suggests which of the child subproblems
alpar@9 40 * should be solved next.
alpar@9 41 *
alpar@9 42 * RETURNS
alpar@9 43 *
alpar@9 44 * The routine ios_choose_var returns the ordinal number of the column
alpar@9 45 * choosen. */
alpar@9 46
alpar@9 47 static int branch_first(glp_tree *T, int *next);
alpar@9 48 static int branch_last(glp_tree *T, int *next);
alpar@9 49 static int branch_mostf(glp_tree *T, int *next);
alpar@9 50 static int branch_drtom(glp_tree *T, int *next);
alpar@9 51
alpar@9 52 int ios_choose_var(glp_tree *T, int *next)
alpar@9 53 { int j;
alpar@9 54 if (T->parm->br_tech == GLP_BR_FFV)
alpar@9 55 { /* branch on first fractional variable */
alpar@9 56 j = branch_first(T, next);
alpar@9 57 }
alpar@9 58 else if (T->parm->br_tech == GLP_BR_LFV)
alpar@9 59 { /* branch on last fractional variable */
alpar@9 60 j = branch_last(T, next);
alpar@9 61 }
alpar@9 62 else if (T->parm->br_tech == GLP_BR_MFV)
alpar@9 63 { /* branch on most fractional variable */
alpar@9 64 j = branch_mostf(T, next);
alpar@9 65 }
alpar@9 66 else if (T->parm->br_tech == GLP_BR_DTH)
alpar@9 67 { /* branch using the heuristic by Dreebeck and Tomlin */
alpar@9 68 j = branch_drtom(T, next);
alpar@9 69 }
alpar@9 70 else if (T->parm->br_tech == GLP_BR_PCH)
alpar@9 71 { /* hybrid pseudocost heuristic */
alpar@9 72 j = ios_pcost_branch(T, next);
alpar@9 73 }
alpar@9 74 else
alpar@9 75 xassert(T != T);
alpar@9 76 return j;
alpar@9 77 }
alpar@9 78
alpar@9 79 /***********************************************************************
alpar@9 80 * branch_first - choose first branching variable
alpar@9 81 *
alpar@9 82 * This routine looks up the list of structural variables and chooses
alpar@9 83 * the first one, which is of integer kind and has fractional value in
alpar@9 84 * optimal solution to the current LP relaxation.
alpar@9 85 *
alpar@9 86 * This routine also selects the branch to be solved next where integer
alpar@9 87 * infeasibility of the chosen variable is less than in other one. */
alpar@9 88
alpar@9 89 static int branch_first(glp_tree *T, int *_next)
alpar@9 90 { int j, next;
alpar@9 91 double beta;
alpar@9 92 /* choose the column to branch on */
alpar@9 93 for (j = 1; j <= T->n; j++)
alpar@9 94 if (T->non_int[j]) break;
alpar@9 95 xassert(1 <= j && j <= T->n);
alpar@9 96 /* select the branch to be solved next */
alpar@9 97 beta = glp_get_col_prim(T->mip, j);
alpar@9 98 if (beta - floor(beta) < ceil(beta) - beta)
alpar@9 99 next = GLP_DN_BRNCH;
alpar@9 100 else
alpar@9 101 next = GLP_UP_BRNCH;
alpar@9 102 *_next = next;
alpar@9 103 return j;
alpar@9 104 }
alpar@9 105
alpar@9 106 /***********************************************************************
alpar@9 107 * branch_last - choose last branching variable
alpar@9 108 *
alpar@9 109 * This routine looks up the list of structural variables and chooses
alpar@9 110 * the last one, which is of integer kind and has fractional value in
alpar@9 111 * optimal solution to the current LP relaxation.
alpar@9 112 *
alpar@9 113 * This routine also selects the branch to be solved next where integer
alpar@9 114 * infeasibility of the chosen variable is less than in other one. */
alpar@9 115
alpar@9 116 static int branch_last(glp_tree *T, int *_next)
alpar@9 117 { int j, next;
alpar@9 118 double beta;
alpar@9 119 /* choose the column to branch on */
alpar@9 120 for (j = T->n; j >= 1; j--)
alpar@9 121 if (T->non_int[j]) break;
alpar@9 122 xassert(1 <= j && j <= T->n);
alpar@9 123 /* select the branch to be solved next */
alpar@9 124 beta = glp_get_col_prim(T->mip, j);
alpar@9 125 if (beta - floor(beta) < ceil(beta) - beta)
alpar@9 126 next = GLP_DN_BRNCH;
alpar@9 127 else
alpar@9 128 next = GLP_UP_BRNCH;
alpar@9 129 *_next = next;
alpar@9 130 return j;
alpar@9 131 }
alpar@9 132
alpar@9 133 /***********************************************************************
alpar@9 134 * branch_mostf - choose most fractional branching variable
alpar@9 135 *
alpar@9 136 * This routine looks up the list of structural variables and chooses
alpar@9 137 * that one, which is of integer kind and has most fractional value in
alpar@9 138 * optimal solution to the current LP relaxation.
alpar@9 139 *
alpar@9 140 * This routine also selects the branch to be solved next where integer
alpar@9 141 * infeasibility of the chosen variable is less than in other one.
alpar@9 142 *
alpar@9 143 * (Alexander Martin notices that "...most infeasible is as good as
alpar@9 144 * random...".) */
alpar@9 145
alpar@9 146 static int branch_mostf(glp_tree *T, int *_next)
alpar@9 147 { int j, jj, next;
alpar@9 148 double beta, most, temp;
alpar@9 149 /* choose the column to branch on */
alpar@9 150 jj = 0, most = DBL_MAX;
alpar@9 151 for (j = 1; j <= T->n; j++)
alpar@9 152 { if (T->non_int[j])
alpar@9 153 { beta = glp_get_col_prim(T->mip, j);
alpar@9 154 temp = floor(beta) + 0.5;
alpar@9 155 if (most > fabs(beta - temp))
alpar@9 156 { jj = j, most = fabs(beta - temp);
alpar@9 157 if (beta < temp)
alpar@9 158 next = GLP_DN_BRNCH;
alpar@9 159 else
alpar@9 160 next = GLP_UP_BRNCH;
alpar@9 161 }
alpar@9 162 }
alpar@9 163 }
alpar@9 164 *_next = next;
alpar@9 165 return jj;
alpar@9 166 }
alpar@9 167
alpar@9 168 /***********************************************************************
alpar@9 169 * branch_drtom - choose branching var using Driebeck-Tomlin heuristic
alpar@9 170 *
alpar@9 171 * This routine chooses a structural variable, which is required to be
alpar@9 172 * integral and has fractional value in optimal solution of the current
alpar@9 173 * LP relaxation, using a heuristic proposed by Driebeck and Tomlin.
alpar@9 174 *
alpar@9 175 * The routine also selects the branch to be solved next, again due to
alpar@9 176 * Driebeck and Tomlin.
alpar@9 177 *
alpar@9 178 * This routine is based on the heuristic proposed in:
alpar@9 179 *
alpar@9 180 * Driebeck N.J. An algorithm for the solution of mixed-integer
alpar@9 181 * programming problems, Management Science, 12: 576-87 (1966);
alpar@9 182 *
alpar@9 183 * and improved in:
alpar@9 184 *
alpar@9 185 * Tomlin J.A. Branch and bound methods for integer and non-convex
alpar@9 186 * programming, in J.Abadie (ed.), Integer and Nonlinear Programming,
alpar@9 187 * North-Holland, Amsterdam, pp. 437-50 (1970).
alpar@9 188 *
alpar@9 189 * Must note that this heuristic is time-expensive, because computing
alpar@9 190 * one-step degradation (see the routine below) requires one BTRAN for
alpar@9 191 * each fractional-valued structural variable. */
alpar@9 192
alpar@9 193 static int branch_drtom(glp_tree *T, int *_next)
alpar@9 194 { glp_prob *mip = T->mip;
alpar@9 195 int m = mip->m;
alpar@9 196 int n = mip->n;
alpar@9 197 char *non_int = T->non_int;
alpar@9 198 int j, jj, k, t, next, kase, len, stat, *ind;
alpar@9 199 double x, dk, alfa, delta_j, delta_k, delta_z, dz_dn, dz_up,
alpar@9 200 dd_dn, dd_up, degrad, *val;
alpar@9 201 /* basic solution of LP relaxation must be optimal */
alpar@9 202 xassert(glp_get_status(mip) == GLP_OPT);
alpar@9 203 /* allocate working arrays */
alpar@9 204 ind = xcalloc(1+n, sizeof(int));
alpar@9 205 val = xcalloc(1+n, sizeof(double));
alpar@9 206 /* nothing has been chosen so far */
alpar@9 207 jj = 0, degrad = -1.0;
alpar@9 208 /* walk through the list of columns (structural variables) */
alpar@9 209 for (j = 1; j <= n; j++)
alpar@9 210 { /* if j-th column is not marked as fractional, skip it */
alpar@9 211 if (!non_int[j]) continue;
alpar@9 212 /* obtain (fractional) value of j-th column in basic solution
alpar@9 213 of LP relaxation */
alpar@9 214 x = glp_get_col_prim(mip, j);
alpar@9 215 /* since the value of j-th column is fractional, the column is
alpar@9 216 basic; compute corresponding row of the simplex table */
alpar@9 217 len = glp_eval_tab_row(mip, m+j, ind, val);
alpar@9 218 /* the following fragment computes a change in the objective
alpar@9 219 function: delta Z = new Z - old Z, where old Z is the
alpar@9 220 objective value in the current optimal basis, and new Z is
alpar@9 221 the objective value in the adjacent basis, for two cases:
alpar@9 222 1) if new upper bound ub' = floor(x[j]) is introduced for
alpar@9 223 j-th column (down branch);
alpar@9 224 2) if new lower bound lb' = ceil(x[j]) is introduced for
alpar@9 225 j-th column (up branch);
alpar@9 226 since in both cases the solution remaining dual feasible
alpar@9 227 becomes primal infeasible, one implicit simplex iteration
alpar@9 228 is performed to determine the change delta Z;
alpar@9 229 it is obvious that new Z, which is never better than old Z,
alpar@9 230 is a lower (minimization) or upper (maximization) bound of
alpar@9 231 the objective function for down- and up-branches. */
alpar@9 232 for (kase = -1; kase <= +1; kase += 2)
alpar@9 233 { /* if kase < 0, the new upper bound of x[j] is introduced;
alpar@9 234 in this case x[j] should decrease in order to leave the
alpar@9 235 basis and go to its new upper bound */
alpar@9 236 /* if kase > 0, the new lower bound of x[j] is introduced;
alpar@9 237 in this case x[j] should increase in order to leave the
alpar@9 238 basis and go to its new lower bound */
alpar@9 239 /* apply the dual ratio test in order to determine which
alpar@9 240 auxiliary or structural variable should enter the basis
alpar@9 241 to keep dual feasibility */
alpar@9 242 k = glp_dual_rtest(mip, len, ind, val, kase, 1e-9);
alpar@9 243 if (k != 0) k = ind[k];
alpar@9 244 /* if no non-basic variable has been chosen, LP relaxation
alpar@9 245 of corresponding branch being primal infeasible and dual
alpar@9 246 unbounded has no primal feasible solution; in this case
alpar@9 247 the change delta Z is formally set to infinity */
alpar@9 248 if (k == 0)
alpar@9 249 { delta_z =
alpar@9 250 (T->mip->dir == GLP_MIN ? +DBL_MAX : -DBL_MAX);
alpar@9 251 goto skip;
alpar@9 252 }
alpar@9 253 /* row of the simplex table that corresponds to non-basic
alpar@9 254 variable x[k] choosen by the dual ratio test is:
alpar@9 255 x[j] = ... + alfa * x[k] + ...
alpar@9 256 where alfa is the influence coefficient (an element of
alpar@9 257 the simplex table row) */
alpar@9 258 /* determine the coefficient alfa */
alpar@9 259 for (t = 1; t <= len; t++) if (ind[t] == k) break;
alpar@9 260 xassert(1 <= t && t <= len);
alpar@9 261 alfa = val[t];
alpar@9 262 /* since in the adjacent basis the variable x[j] becomes
alpar@9 263 non-basic, knowing its value in the current basis we can
alpar@9 264 determine its change delta x[j] = new x[j] - old x[j] */
alpar@9 265 delta_j = (kase < 0 ? floor(x) : ceil(x)) - x;
alpar@9 266 /* and knowing the coefficient alfa we can determine the
alpar@9 267 corresponding change delta x[k] = new x[k] - old x[k],
alpar@9 268 where old x[k] is a value of x[k] in the current basis,
alpar@9 269 and new x[k] is a value of x[k] in the adjacent basis */
alpar@9 270 delta_k = delta_j / alfa;
alpar@9 271 /* Tomlin noticed that if the variable x[k] is of integer
alpar@9 272 kind, its change cannot be less (eventually) than one in
alpar@9 273 the magnitude */
alpar@9 274 if (k > m && glp_get_col_kind(mip, k-m) != GLP_CV)
alpar@9 275 { /* x[k] is structural integer variable */
alpar@9 276 if (fabs(delta_k - floor(delta_k + 0.5)) > 1e-3)
alpar@9 277 { if (delta_k > 0.0)
alpar@9 278 delta_k = ceil(delta_k); /* +3.14 -> +4 */
alpar@9 279 else
alpar@9 280 delta_k = floor(delta_k); /* -3.14 -> -4 */
alpar@9 281 }
alpar@9 282 }
alpar@9 283 /* now determine the status and reduced cost of x[k] in the
alpar@9 284 current basis */
alpar@9 285 if (k <= m)
alpar@9 286 { stat = glp_get_row_stat(mip, k);
alpar@9 287 dk = glp_get_row_dual(mip, k);
alpar@9 288 }
alpar@9 289 else
alpar@9 290 { stat = glp_get_col_stat(mip, k-m);
alpar@9 291 dk = glp_get_col_dual(mip, k-m);
alpar@9 292 }
alpar@9 293 /* if the current basis is dual degenerate, some reduced
alpar@9 294 costs which are close to zero may have wrong sign due to
alpar@9 295 round-off errors, so correct the sign of d[k] */
alpar@9 296 switch (T->mip->dir)
alpar@9 297 { case GLP_MIN:
alpar@9 298 if (stat == GLP_NL && dk < 0.0 ||
alpar@9 299 stat == GLP_NU && dk > 0.0 ||
alpar@9 300 stat == GLP_NF) dk = 0.0;
alpar@9 301 break;
alpar@9 302 case GLP_MAX:
alpar@9 303 if (stat == GLP_NL && dk > 0.0 ||
alpar@9 304 stat == GLP_NU && dk < 0.0 ||
alpar@9 305 stat == GLP_NF) dk = 0.0;
alpar@9 306 break;
alpar@9 307 default:
alpar@9 308 xassert(T != T);
alpar@9 309 }
alpar@9 310 /* now knowing the change of x[k] and its reduced cost d[k]
alpar@9 311 we can compute the corresponding change in the objective
alpar@9 312 function delta Z = new Z - old Z = d[k] * delta x[k];
alpar@9 313 note that due to Tomlin's modification new Z can be even
alpar@9 314 worse than in the adjacent basis */
alpar@9 315 delta_z = dk * delta_k;
alpar@9 316 skip: /* new Z is never better than old Z, therefore the change
alpar@9 317 delta Z is always non-negative (in case of minimization)
alpar@9 318 or non-positive (in case of maximization) */
alpar@9 319 switch (T->mip->dir)
alpar@9 320 { case GLP_MIN: xassert(delta_z >= 0.0); break;
alpar@9 321 case GLP_MAX: xassert(delta_z <= 0.0); break;
alpar@9 322 default: xassert(T != T);
alpar@9 323 }
alpar@9 324 /* save the change in the objective fnction for down- and
alpar@9 325 up-branches, respectively */
alpar@9 326 if (kase < 0) dz_dn = delta_z; else dz_up = delta_z;
alpar@9 327 }
alpar@9 328 /* thus, in down-branch no integer feasible solution can be
alpar@9 329 better than Z + dz_dn, and in up-branch no integer feasible
alpar@9 330 solution can be better than Z + dz_up, where Z is value of
alpar@9 331 the objective function in the current basis */
alpar@9 332 /* following the heuristic by Driebeck and Tomlin we choose a
alpar@9 333 column (i.e. structural variable) which provides largest
alpar@9 334 degradation of the objective function in some of branches;
alpar@9 335 besides, we select the branch with smaller degradation to
alpar@9 336 be solved next and keep other branch with larger degradation
alpar@9 337 in the active list hoping to minimize the number of further
alpar@9 338 backtrackings */
alpar@9 339 if (degrad < fabs(dz_dn) || degrad < fabs(dz_up))
alpar@9 340 { jj = j;
alpar@9 341 if (fabs(dz_dn) < fabs(dz_up))
alpar@9 342 { /* select down branch to be solved next */
alpar@9 343 next = GLP_DN_BRNCH;
alpar@9 344 degrad = fabs(dz_up);
alpar@9 345 }
alpar@9 346 else
alpar@9 347 { /* select up branch to be solved next */
alpar@9 348 next = GLP_UP_BRNCH;
alpar@9 349 degrad = fabs(dz_dn);
alpar@9 350 }
alpar@9 351 /* save the objective changes for printing */
alpar@9 352 dd_dn = dz_dn, dd_up = dz_up;
alpar@9 353 /* if down- or up-branch has no feasible solution, we does
alpar@9 354 not need to consider other candidates (in principle, the
alpar@9 355 corresponding branch could be pruned right now) */
alpar@9 356 if (degrad == DBL_MAX) break;
alpar@9 357 }
alpar@9 358 }
alpar@9 359 /* free working arrays */
alpar@9 360 xfree(ind);
alpar@9 361 xfree(val);
alpar@9 362 /* something must be chosen */
alpar@9 363 xassert(1 <= jj && jj <= n);
alpar@9 364 #if 1 /* 02/XI-2009 */
alpar@9 365 if (degrad < 1e-6 * (1.0 + 0.001 * fabs(mip->obj_val)))
alpar@9 366 { jj = branch_mostf(T, &next);
alpar@9 367 goto done;
alpar@9 368 }
alpar@9 369 #endif
alpar@9 370 if (T->parm->msg_lev >= GLP_MSG_DBG)
alpar@9 371 { xprintf("branch_drtom: column %d chosen to branch on\n", jj);
alpar@9 372 if (fabs(dd_dn) == DBL_MAX)
alpar@9 373 xprintf("branch_drtom: down-branch is infeasible\n");
alpar@9 374 else
alpar@9 375 xprintf("branch_drtom: down-branch bound is %.9e\n",
alpar@9 376 lpx_get_obj_val(mip) + dd_dn);
alpar@9 377 if (fabs(dd_up) == DBL_MAX)
alpar@9 378 xprintf("branch_drtom: up-branch is infeasible\n");
alpar@9 379 else
alpar@9 380 xprintf("branch_drtom: up-branch bound is %.9e\n",
alpar@9 381 lpx_get_obj_val(mip) + dd_up);
alpar@9 382 }
alpar@9 383 done: *_next = next;
alpar@9 384 return jj;
alpar@9 385 }
alpar@9 386
alpar@9 387 /**********************************************************************/
alpar@9 388
alpar@9 389 struct csa
alpar@9 390 { /* common storage area */
alpar@9 391 int *dn_cnt; /* int dn_cnt[1+n]; */
alpar@9 392 /* dn_cnt[j] is the number of subproblems, whose LP relaxations
alpar@9 393 have been solved and which are down-branches for variable x[j];
alpar@9 394 dn_cnt[j] = 0 means the down pseudocost is uninitialized */
alpar@9 395 double *dn_sum; /* double dn_sum[1+n]; */
alpar@9 396 /* dn_sum[j] is the sum of per unit degradations of the objective
alpar@9 397 over all dn_cnt[j] subproblems */
alpar@9 398 int *up_cnt; /* int up_cnt[1+n]; */
alpar@9 399 /* up_cnt[j] is the number of subproblems, whose LP relaxations
alpar@9 400 have been solved and which are up-branches for variable x[j];
alpar@9 401 up_cnt[j] = 0 means the up pseudocost is uninitialized */
alpar@9 402 double *up_sum; /* double up_sum[1+n]; */
alpar@9 403 /* up_sum[j] is the sum of per unit degradations of the objective
alpar@9 404 over all up_cnt[j] subproblems */
alpar@9 405 };
alpar@9 406
alpar@9 407 void *ios_pcost_init(glp_tree *tree)
alpar@9 408 { /* initialize working data used on pseudocost branching */
alpar@9 409 struct csa *csa;
alpar@9 410 int n = tree->n, j;
alpar@9 411 csa = xmalloc(sizeof(struct csa));
alpar@9 412 csa->dn_cnt = xcalloc(1+n, sizeof(int));
alpar@9 413 csa->dn_sum = xcalloc(1+n, sizeof(double));
alpar@9 414 csa->up_cnt = xcalloc(1+n, sizeof(int));
alpar@9 415 csa->up_sum = xcalloc(1+n, sizeof(double));
alpar@9 416 for (j = 1; j <= n; j++)
alpar@9 417 { csa->dn_cnt[j] = csa->up_cnt[j] = 0;
alpar@9 418 csa->dn_sum[j] = csa->up_sum[j] = 0.0;
alpar@9 419 }
alpar@9 420 return csa;
alpar@9 421 }
alpar@9 422
alpar@9 423 static double eval_degrad(glp_prob *P, int j, double bnd)
alpar@9 424 { /* compute degradation of the objective on fixing x[j] at given
alpar@9 425 value with a limited number of dual simplex iterations */
alpar@9 426 /* this routine fixes column x[j] at specified value bnd,
alpar@9 427 solves resulting LP, and returns a lower bound to degradation
alpar@9 428 of the objective, degrad >= 0 */
alpar@9 429 glp_prob *lp;
alpar@9 430 glp_smcp parm;
alpar@9 431 int ret;
alpar@9 432 double degrad;
alpar@9 433 /* the current basis must be optimal */
alpar@9 434 xassert(glp_get_status(P) == GLP_OPT);
alpar@9 435 /* create a copy of P */
alpar@9 436 lp = glp_create_prob();
alpar@9 437 glp_copy_prob(lp, P, 0);
alpar@9 438 /* fix column x[j] at specified value */
alpar@9 439 glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd);
alpar@9 440 /* try to solve resulting LP */
alpar@9 441 glp_init_smcp(&parm);
alpar@9 442 parm.msg_lev = GLP_MSG_OFF;
alpar@9 443 parm.meth = GLP_DUAL;
alpar@9 444 parm.it_lim = 30;
alpar@9 445 parm.out_dly = 1000;
alpar@9 446 parm.meth = GLP_DUAL;
alpar@9 447 ret = glp_simplex(lp, &parm);
alpar@9 448 if (ret == 0 || ret == GLP_EITLIM)
alpar@9 449 { if (glp_get_prim_stat(lp) == GLP_NOFEAS)
alpar@9 450 { /* resulting LP has no primal feasible solution */
alpar@9 451 degrad = DBL_MAX;
alpar@9 452 }
alpar@9 453 else if (glp_get_dual_stat(lp) == GLP_FEAS)
alpar@9 454 { /* resulting basis is optimal or at least dual feasible,
alpar@9 455 so we have the correct lower bound to degradation */
alpar@9 456 if (P->dir == GLP_MIN)
alpar@9 457 degrad = lp->obj_val - P->obj_val;
alpar@9 458 else if (P->dir == GLP_MAX)
alpar@9 459 degrad = P->obj_val - lp->obj_val;
alpar@9 460 else
alpar@9 461 xassert(P != P);
alpar@9 462 /* degradation cannot be negative by definition */
alpar@9 463 /* note that the lower bound to degradation may be close
alpar@9 464 to zero even if its exact value is zero due to round-off
alpar@9 465 errors on computing the objective value */
alpar@9 466 if (degrad < 1e-6 * (1.0 + 0.001 * fabs(P->obj_val)))
alpar@9 467 degrad = 0.0;
alpar@9 468 }
alpar@9 469 else
alpar@9 470 { /* the final basis reported by the simplex solver is dual
alpar@9 471 infeasible, so we cannot determine a non-trivial lower
alpar@9 472 bound to degradation */
alpar@9 473 degrad = 0.0;
alpar@9 474 }
alpar@9 475 }
alpar@9 476 else
alpar@9 477 { /* the simplex solver failed */
alpar@9 478 degrad = 0.0;
alpar@9 479 }
alpar@9 480 /* delete the copy of P */
alpar@9 481 glp_delete_prob(lp);
alpar@9 482 return degrad;
alpar@9 483 }
alpar@9 484
alpar@9 485 void ios_pcost_update(glp_tree *tree)
alpar@9 486 { /* update history information for pseudocost branching */
alpar@9 487 /* this routine is called every time when LP relaxation of the
alpar@9 488 current subproblem has been solved to optimality with all lazy
alpar@9 489 and cutting plane constraints included */
alpar@9 490 int j;
alpar@9 491 double dx, dz, psi;
alpar@9 492 struct csa *csa = tree->pcost;
alpar@9 493 xassert(csa != NULL);
alpar@9 494 xassert(tree->curr != NULL);
alpar@9 495 /* if the current subproblem is the root, skip updating */
alpar@9 496 if (tree->curr->up == NULL) goto skip;
alpar@9 497 /* determine branching variable x[j], which was used in the
alpar@9 498 parent subproblem to create the current subproblem */
alpar@9 499 j = tree->curr->up->br_var;
alpar@9 500 xassert(1 <= j && j <= tree->n);
alpar@9 501 /* determine the change dx[j] = new x[j] - old x[j],
alpar@9 502 where new x[j] is a value of x[j] in optimal solution to LP
alpar@9 503 relaxation of the current subproblem, old x[j] is a value of
alpar@9 504 x[j] in optimal solution to LP relaxation of the parent
alpar@9 505 subproblem */
alpar@9 506 dx = tree->mip->col[j]->prim - tree->curr->up->br_val;
alpar@9 507 xassert(dx != 0.0);
alpar@9 508 /* determine corresponding change dz = new dz - old dz in the
alpar@9 509 objective function value */
alpar@9 510 dz = tree->mip->obj_val - tree->curr->up->lp_obj;
alpar@9 511 /* determine per unit degradation of the objective function */
alpar@9 512 psi = fabs(dz / dx);
alpar@9 513 /* update history information */
alpar@9 514 if (dx < 0.0)
alpar@9 515 { /* the current subproblem is down-branch */
alpar@9 516 csa->dn_cnt[j]++;
alpar@9 517 csa->dn_sum[j] += psi;
alpar@9 518 }
alpar@9 519 else /* dx > 0.0 */
alpar@9 520 { /* the current subproblem is up-branch */
alpar@9 521 csa->up_cnt[j]++;
alpar@9 522 csa->up_sum[j] += psi;
alpar@9 523 }
alpar@9 524 skip: return;
alpar@9 525 }
alpar@9 526
alpar@9 527 void ios_pcost_free(glp_tree *tree)
alpar@9 528 { /* free working area used on pseudocost branching */
alpar@9 529 struct csa *csa = tree->pcost;
alpar@9 530 xassert(csa != NULL);
alpar@9 531 xfree(csa->dn_cnt);
alpar@9 532 xfree(csa->dn_sum);
alpar@9 533 xfree(csa->up_cnt);
alpar@9 534 xfree(csa->up_sum);
alpar@9 535 xfree(csa);
alpar@9 536 tree->pcost = NULL;
alpar@9 537 return;
alpar@9 538 }
alpar@9 539
alpar@9 540 static double eval_psi(glp_tree *T, int j, int brnch)
alpar@9 541 { /* compute estimation of pseudocost of variable x[j] for down-
alpar@9 542 or up-branch */
alpar@9 543 struct csa *csa = T->pcost;
alpar@9 544 double beta, degrad, psi;
alpar@9 545 xassert(csa != NULL);
alpar@9 546 xassert(1 <= j && j <= T->n);
alpar@9 547 if (brnch == GLP_DN_BRNCH)
alpar@9 548 { /* down-branch */
alpar@9 549 if (csa->dn_cnt[j] == 0)
alpar@9 550 { /* initialize down pseudocost */
alpar@9 551 beta = T->mip->col[j]->prim;
alpar@9 552 degrad = eval_degrad(T->mip, j, floor(beta));
alpar@9 553 if (degrad == DBL_MAX)
alpar@9 554 { psi = DBL_MAX;
alpar@9 555 goto done;
alpar@9 556 }
alpar@9 557 csa->dn_cnt[j] = 1;
alpar@9 558 csa->dn_sum[j] = degrad / (beta - floor(beta));
alpar@9 559 }
alpar@9 560 psi = csa->dn_sum[j] / (double)csa->dn_cnt[j];
alpar@9 561 }
alpar@9 562 else if (brnch == GLP_UP_BRNCH)
alpar@9 563 { /* up-branch */
alpar@9 564 if (csa->up_cnt[j] == 0)
alpar@9 565 { /* initialize up pseudocost */
alpar@9 566 beta = T->mip->col[j]->prim;
alpar@9 567 degrad = eval_degrad(T->mip, j, ceil(beta));
alpar@9 568 if (degrad == DBL_MAX)
alpar@9 569 { psi = DBL_MAX;
alpar@9 570 goto done;
alpar@9 571 }
alpar@9 572 csa->up_cnt[j] = 1;
alpar@9 573 csa->up_sum[j] = degrad / (ceil(beta) - beta);
alpar@9 574 }
alpar@9 575 psi = csa->up_sum[j] / (double)csa->up_cnt[j];
alpar@9 576 }
alpar@9 577 else
alpar@9 578 xassert(brnch != brnch);
alpar@9 579 done: return psi;
alpar@9 580 }
alpar@9 581
alpar@9 582 static void progress(glp_tree *T)
alpar@9 583 { /* display progress of pseudocost initialization */
alpar@9 584 struct csa *csa = T->pcost;
alpar@9 585 int j, nv = 0, ni = 0;
alpar@9 586 for (j = 1; j <= T->n; j++)
alpar@9 587 { if (glp_ios_can_branch(T, j))
alpar@9 588 { nv++;
alpar@9 589 if (csa->dn_cnt[j] > 0 && csa->up_cnt[j] > 0) ni++;
alpar@9 590 }
alpar@9 591 }
alpar@9 592 xprintf("Pseudocosts initialized for %d of %d variables\n",
alpar@9 593 ni, nv);
alpar@9 594 return;
alpar@9 595 }
alpar@9 596
alpar@9 597 int ios_pcost_branch(glp_tree *T, int *_next)
alpar@9 598 { /* choose branching variable with pseudocost branching */
alpar@9 599 glp_long t = xtime();
alpar@9 600 int j, jjj, sel;
alpar@9 601 double beta, psi, d1, d2, d, dmax;
alpar@9 602 /* initialize the working arrays */
alpar@9 603 if (T->pcost == NULL)
alpar@9 604 T->pcost = ios_pcost_init(T);
alpar@9 605 /* nothing has been chosen so far */
alpar@9 606 jjj = 0, dmax = -1.0;
alpar@9 607 /* go through the list of branching candidates */
alpar@9 608 for (j = 1; j <= T->n; j++)
alpar@9 609 { if (!glp_ios_can_branch(T, j)) continue;
alpar@9 610 /* determine primal value of x[j] in optimal solution to LP
alpar@9 611 relaxation of the current subproblem */
alpar@9 612 beta = T->mip->col[j]->prim;
alpar@9 613 /* estimate pseudocost of x[j] for down-branch */
alpar@9 614 psi = eval_psi(T, j, GLP_DN_BRNCH);
alpar@9 615 if (psi == DBL_MAX)
alpar@9 616 { /* down-branch has no primal feasible solution */
alpar@9 617 jjj = j, sel = GLP_DN_BRNCH;
alpar@9 618 goto done;
alpar@9 619 }
alpar@9 620 /* estimate degradation of the objective for down-branch */
alpar@9 621 d1 = psi * (beta - floor(beta));
alpar@9 622 /* estimate pseudocost of x[j] for up-branch */
alpar@9 623 psi = eval_psi(T, j, GLP_UP_BRNCH);
alpar@9 624 if (psi == DBL_MAX)
alpar@9 625 { /* up-branch has no primal feasible solution */
alpar@9 626 jjj = j, sel = GLP_UP_BRNCH;
alpar@9 627 goto done;
alpar@9 628 }
alpar@9 629 /* estimate degradation of the objective for up-branch */
alpar@9 630 d2 = psi * (ceil(beta) - beta);
alpar@9 631 /* determine d = max(d1, d2) */
alpar@9 632 d = (d1 > d2 ? d1 : d2);
alpar@9 633 /* choose x[j] which provides maximal estimated degradation of
alpar@9 634 the objective either in down- or up-branch */
alpar@9 635 if (dmax < d)
alpar@9 636 { dmax = d;
alpar@9 637 jjj = j;
alpar@9 638 /* continue the search from a subproblem, where degradation
alpar@9 639 is less than in other one */
alpar@9 640 sel = (d1 <= d2 ? GLP_DN_BRNCH : GLP_UP_BRNCH);
alpar@9 641 }
alpar@9 642 /* display progress of pseudocost initialization */
alpar@9 643 if (T->parm->msg_lev >= GLP_ON)
alpar@9 644 { if (xdifftime(xtime(), t) >= 10.0)
alpar@9 645 { progress(T);
alpar@9 646 t = xtime();
alpar@9 647 }
alpar@9 648 }
alpar@9 649 }
alpar@9 650 if (dmax == 0.0)
alpar@9 651 { /* no degradation is indicated; choose a variable having most
alpar@9 652 fractional value */
alpar@9 653 jjj = branch_mostf(T, &sel);
alpar@9 654 }
alpar@9 655 done: *_next = sel;
alpar@9 656 return jjj;
alpar@9 657 }
alpar@9 658
alpar@9 659 /* eof */