1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpios09.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,659 @@
1.4 +/* glpios09.c (branching heuristics) */
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 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 +
1.30 +/***********************************************************************
1.31 +* NAME
1.32 +*
1.33 +* ios_choose_var - select variable to branch on
1.34 +*
1.35 +* SYNOPSIS
1.36 +*
1.37 +* #include "glpios.h"
1.38 +* int ios_choose_var(glp_tree *T, int *next);
1.39 +*
1.40 +* The routine ios_choose_var chooses a variable from the candidate
1.41 +* list to branch on. Additionally the routine provides a flag stored
1.42 +* in the location next to suggests which of the child subproblems
1.43 +* should be solved next.
1.44 +*
1.45 +* RETURNS
1.46 +*
1.47 +* The routine ios_choose_var returns the ordinal number of the column
1.48 +* choosen. */
1.49 +
1.50 +static int branch_first(glp_tree *T, int *next);
1.51 +static int branch_last(glp_tree *T, int *next);
1.52 +static int branch_mostf(glp_tree *T, int *next);
1.53 +static int branch_drtom(glp_tree *T, int *next);
1.54 +
1.55 +int ios_choose_var(glp_tree *T, int *next)
1.56 +{ int j;
1.57 + if (T->parm->br_tech == GLP_BR_FFV)
1.58 + { /* branch on first fractional variable */
1.59 + j = branch_first(T, next);
1.60 + }
1.61 + else if (T->parm->br_tech == GLP_BR_LFV)
1.62 + { /* branch on last fractional variable */
1.63 + j = branch_last(T, next);
1.64 + }
1.65 + else if (T->parm->br_tech == GLP_BR_MFV)
1.66 + { /* branch on most fractional variable */
1.67 + j = branch_mostf(T, next);
1.68 + }
1.69 + else if (T->parm->br_tech == GLP_BR_DTH)
1.70 + { /* branch using the heuristic by Dreebeck and Tomlin */
1.71 + j = branch_drtom(T, next);
1.72 + }
1.73 + else if (T->parm->br_tech == GLP_BR_PCH)
1.74 + { /* hybrid pseudocost heuristic */
1.75 + j = ios_pcost_branch(T, next);
1.76 + }
1.77 + else
1.78 + xassert(T != T);
1.79 + return j;
1.80 +}
1.81 +
1.82 +/***********************************************************************
1.83 +* branch_first - choose first branching variable
1.84 +*
1.85 +* This routine looks up the list of structural variables and chooses
1.86 +* the first one, which is of integer kind and has fractional value in
1.87 +* optimal solution to the current LP relaxation.
1.88 +*
1.89 +* This routine also selects the branch to be solved next where integer
1.90 +* infeasibility of the chosen variable is less than in other one. */
1.91 +
1.92 +static int branch_first(glp_tree *T, int *_next)
1.93 +{ int j, next;
1.94 + double beta;
1.95 + /* choose the column to branch on */
1.96 + for (j = 1; j <= T->n; j++)
1.97 + if (T->non_int[j]) break;
1.98 + xassert(1 <= j && j <= T->n);
1.99 + /* select the branch to be solved next */
1.100 + beta = glp_get_col_prim(T->mip, j);
1.101 + if (beta - floor(beta) < ceil(beta) - beta)
1.102 + next = GLP_DN_BRNCH;
1.103 + else
1.104 + next = GLP_UP_BRNCH;
1.105 + *_next = next;
1.106 + return j;
1.107 +}
1.108 +
1.109 +/***********************************************************************
1.110 +* branch_last - choose last branching variable
1.111 +*
1.112 +* This routine looks up the list of structural variables and chooses
1.113 +* the last one, which is of integer kind and has fractional value in
1.114 +* optimal solution to the current LP relaxation.
1.115 +*
1.116 +* This routine also selects the branch to be solved next where integer
1.117 +* infeasibility of the chosen variable is less than in other one. */
1.118 +
1.119 +static int branch_last(glp_tree *T, int *_next)
1.120 +{ int j, next;
1.121 + double beta;
1.122 + /* choose the column to branch on */
1.123 + for (j = T->n; j >= 1; j--)
1.124 + if (T->non_int[j]) break;
1.125 + xassert(1 <= j && j <= T->n);
1.126 + /* select the branch to be solved next */
1.127 + beta = glp_get_col_prim(T->mip, j);
1.128 + if (beta - floor(beta) < ceil(beta) - beta)
1.129 + next = GLP_DN_BRNCH;
1.130 + else
1.131 + next = GLP_UP_BRNCH;
1.132 + *_next = next;
1.133 + return j;
1.134 +}
1.135 +
1.136 +/***********************************************************************
1.137 +* branch_mostf - choose most fractional branching variable
1.138 +*
1.139 +* This routine looks up the list of structural variables and chooses
1.140 +* that one, which is of integer kind and has most fractional value in
1.141 +* optimal solution to the current LP relaxation.
1.142 +*
1.143 +* This routine also selects the branch to be solved next where integer
1.144 +* infeasibility of the chosen variable is less than in other one.
1.145 +*
1.146 +* (Alexander Martin notices that "...most infeasible is as good as
1.147 +* random...".) */
1.148 +
1.149 +static int branch_mostf(glp_tree *T, int *_next)
1.150 +{ int j, jj, next;
1.151 + double beta, most, temp;
1.152 + /* choose the column to branch on */
1.153 + jj = 0, most = DBL_MAX;
1.154 + for (j = 1; j <= T->n; j++)
1.155 + { if (T->non_int[j])
1.156 + { beta = glp_get_col_prim(T->mip, j);
1.157 + temp = floor(beta) + 0.5;
1.158 + if (most > fabs(beta - temp))
1.159 + { jj = j, most = fabs(beta - temp);
1.160 + if (beta < temp)
1.161 + next = GLP_DN_BRNCH;
1.162 + else
1.163 + next = GLP_UP_BRNCH;
1.164 + }
1.165 + }
1.166 + }
1.167 + *_next = next;
1.168 + return jj;
1.169 +}
1.170 +
1.171 +/***********************************************************************
1.172 +* branch_drtom - choose branching var using Driebeck-Tomlin heuristic
1.173 +*
1.174 +* This routine chooses a structural variable, which is required to be
1.175 +* integral and has fractional value in optimal solution of the current
1.176 +* LP relaxation, using a heuristic proposed by Driebeck and Tomlin.
1.177 +*
1.178 +* The routine also selects the branch to be solved next, again due to
1.179 +* Driebeck and Tomlin.
1.180 +*
1.181 +* This routine is based on the heuristic proposed in:
1.182 +*
1.183 +* Driebeck N.J. An algorithm for the solution of mixed-integer
1.184 +* programming problems, Management Science, 12: 576-87 (1966);
1.185 +*
1.186 +* and improved in:
1.187 +*
1.188 +* Tomlin J.A. Branch and bound methods for integer and non-convex
1.189 +* programming, in J.Abadie (ed.), Integer and Nonlinear Programming,
1.190 +* North-Holland, Amsterdam, pp. 437-50 (1970).
1.191 +*
1.192 +* Must note that this heuristic is time-expensive, because computing
1.193 +* one-step degradation (see the routine below) requires one BTRAN for
1.194 +* each fractional-valued structural variable. */
1.195 +
1.196 +static int branch_drtom(glp_tree *T, int *_next)
1.197 +{ glp_prob *mip = T->mip;
1.198 + int m = mip->m;
1.199 + int n = mip->n;
1.200 + char *non_int = T->non_int;
1.201 + int j, jj, k, t, next, kase, len, stat, *ind;
1.202 + double x, dk, alfa, delta_j, delta_k, delta_z, dz_dn, dz_up,
1.203 + dd_dn, dd_up, degrad, *val;
1.204 + /* basic solution of LP relaxation must be optimal */
1.205 + xassert(glp_get_status(mip) == GLP_OPT);
1.206 + /* allocate working arrays */
1.207 + ind = xcalloc(1+n, sizeof(int));
1.208 + val = xcalloc(1+n, sizeof(double));
1.209 + /* nothing has been chosen so far */
1.210 + jj = 0, degrad = -1.0;
1.211 + /* walk through the list of columns (structural variables) */
1.212 + for (j = 1; j <= n; j++)
1.213 + { /* if j-th column is not marked as fractional, skip it */
1.214 + if (!non_int[j]) continue;
1.215 + /* obtain (fractional) value of j-th column in basic solution
1.216 + of LP relaxation */
1.217 + x = glp_get_col_prim(mip, j);
1.218 + /* since the value of j-th column is fractional, the column is
1.219 + basic; compute corresponding row of the simplex table */
1.220 + len = glp_eval_tab_row(mip, m+j, ind, val);
1.221 + /* the following fragment computes a change in the objective
1.222 + function: delta Z = new Z - old Z, where old Z is the
1.223 + objective value in the current optimal basis, and new Z is
1.224 + the objective value in the adjacent basis, for two cases:
1.225 + 1) if new upper bound ub' = floor(x[j]) is introduced for
1.226 + j-th column (down branch);
1.227 + 2) if new lower bound lb' = ceil(x[j]) is introduced for
1.228 + j-th column (up branch);
1.229 + since in both cases the solution remaining dual feasible
1.230 + becomes primal infeasible, one implicit simplex iteration
1.231 + is performed to determine the change delta Z;
1.232 + it is obvious that new Z, which is never better than old Z,
1.233 + is a lower (minimization) or upper (maximization) bound of
1.234 + the objective function for down- and up-branches. */
1.235 + for (kase = -1; kase <= +1; kase += 2)
1.236 + { /* if kase < 0, the new upper bound of x[j] is introduced;
1.237 + in this case x[j] should decrease in order to leave the
1.238 + basis and go to its new upper bound */
1.239 + /* if kase > 0, the new lower bound of x[j] is introduced;
1.240 + in this case x[j] should increase in order to leave the
1.241 + basis and go to its new lower bound */
1.242 + /* apply the dual ratio test in order to determine which
1.243 + auxiliary or structural variable should enter the basis
1.244 + to keep dual feasibility */
1.245 + k = glp_dual_rtest(mip, len, ind, val, kase, 1e-9);
1.246 + if (k != 0) k = ind[k];
1.247 + /* if no non-basic variable has been chosen, LP relaxation
1.248 + of corresponding branch being primal infeasible and dual
1.249 + unbounded has no primal feasible solution; in this case
1.250 + the change delta Z is formally set to infinity */
1.251 + if (k == 0)
1.252 + { delta_z =
1.253 + (T->mip->dir == GLP_MIN ? +DBL_MAX : -DBL_MAX);
1.254 + goto skip;
1.255 + }
1.256 + /* row of the simplex table that corresponds to non-basic
1.257 + variable x[k] choosen by the dual ratio test is:
1.258 + x[j] = ... + alfa * x[k] + ...
1.259 + where alfa is the influence coefficient (an element of
1.260 + the simplex table row) */
1.261 + /* determine the coefficient alfa */
1.262 + for (t = 1; t <= len; t++) if (ind[t] == k) break;
1.263 + xassert(1 <= t && t <= len);
1.264 + alfa = val[t];
1.265 + /* since in the adjacent basis the variable x[j] becomes
1.266 + non-basic, knowing its value in the current basis we can
1.267 + determine its change delta x[j] = new x[j] - old x[j] */
1.268 + delta_j = (kase < 0 ? floor(x) : ceil(x)) - x;
1.269 + /* and knowing the coefficient alfa we can determine the
1.270 + corresponding change delta x[k] = new x[k] - old x[k],
1.271 + where old x[k] is a value of x[k] in the current basis,
1.272 + and new x[k] is a value of x[k] in the adjacent basis */
1.273 + delta_k = delta_j / alfa;
1.274 + /* Tomlin noticed that if the variable x[k] is of integer
1.275 + kind, its change cannot be less (eventually) than one in
1.276 + the magnitude */
1.277 + if (k > m && glp_get_col_kind(mip, k-m) != GLP_CV)
1.278 + { /* x[k] is structural integer variable */
1.279 + if (fabs(delta_k - floor(delta_k + 0.5)) > 1e-3)
1.280 + { if (delta_k > 0.0)
1.281 + delta_k = ceil(delta_k); /* +3.14 -> +4 */
1.282 + else
1.283 + delta_k = floor(delta_k); /* -3.14 -> -4 */
1.284 + }
1.285 + }
1.286 + /* now determine the status and reduced cost of x[k] in the
1.287 + current basis */
1.288 + if (k <= m)
1.289 + { stat = glp_get_row_stat(mip, k);
1.290 + dk = glp_get_row_dual(mip, k);
1.291 + }
1.292 + else
1.293 + { stat = glp_get_col_stat(mip, k-m);
1.294 + dk = glp_get_col_dual(mip, k-m);
1.295 + }
1.296 + /* if the current basis is dual degenerate, some reduced
1.297 + costs which are close to zero may have wrong sign due to
1.298 + round-off errors, so correct the sign of d[k] */
1.299 + switch (T->mip->dir)
1.300 + { case GLP_MIN:
1.301 + if (stat == GLP_NL && dk < 0.0 ||
1.302 + stat == GLP_NU && dk > 0.0 ||
1.303 + stat == GLP_NF) dk = 0.0;
1.304 + break;
1.305 + case GLP_MAX:
1.306 + if (stat == GLP_NL && dk > 0.0 ||
1.307 + stat == GLP_NU && dk < 0.0 ||
1.308 + stat == GLP_NF) dk = 0.0;
1.309 + break;
1.310 + default:
1.311 + xassert(T != T);
1.312 + }
1.313 + /* now knowing the change of x[k] and its reduced cost d[k]
1.314 + we can compute the corresponding change in the objective
1.315 + function delta Z = new Z - old Z = d[k] * delta x[k];
1.316 + note that due to Tomlin's modification new Z can be even
1.317 + worse than in the adjacent basis */
1.318 + delta_z = dk * delta_k;
1.319 +skip: /* new Z is never better than old Z, therefore the change
1.320 + delta Z is always non-negative (in case of minimization)
1.321 + or non-positive (in case of maximization) */
1.322 + switch (T->mip->dir)
1.323 + { case GLP_MIN: xassert(delta_z >= 0.0); break;
1.324 + case GLP_MAX: xassert(delta_z <= 0.0); break;
1.325 + default: xassert(T != T);
1.326 + }
1.327 + /* save the change in the objective fnction for down- and
1.328 + up-branches, respectively */
1.329 + if (kase < 0) dz_dn = delta_z; else dz_up = delta_z;
1.330 + }
1.331 + /* thus, in down-branch no integer feasible solution can be
1.332 + better than Z + dz_dn, and in up-branch no integer feasible
1.333 + solution can be better than Z + dz_up, where Z is value of
1.334 + the objective function in the current basis */
1.335 + /* following the heuristic by Driebeck and Tomlin we choose a
1.336 + column (i.e. structural variable) which provides largest
1.337 + degradation of the objective function in some of branches;
1.338 + besides, we select the branch with smaller degradation to
1.339 + be solved next and keep other branch with larger degradation
1.340 + in the active list hoping to minimize the number of further
1.341 + backtrackings */
1.342 + if (degrad < fabs(dz_dn) || degrad < fabs(dz_up))
1.343 + { jj = j;
1.344 + if (fabs(dz_dn) < fabs(dz_up))
1.345 + { /* select down branch to be solved next */
1.346 + next = GLP_DN_BRNCH;
1.347 + degrad = fabs(dz_up);
1.348 + }
1.349 + else
1.350 + { /* select up branch to be solved next */
1.351 + next = GLP_UP_BRNCH;
1.352 + degrad = fabs(dz_dn);
1.353 + }
1.354 + /* save the objective changes for printing */
1.355 + dd_dn = dz_dn, dd_up = dz_up;
1.356 + /* if down- or up-branch has no feasible solution, we does
1.357 + not need to consider other candidates (in principle, the
1.358 + corresponding branch could be pruned right now) */
1.359 + if (degrad == DBL_MAX) break;
1.360 + }
1.361 + }
1.362 + /* free working arrays */
1.363 + xfree(ind);
1.364 + xfree(val);
1.365 + /* something must be chosen */
1.366 + xassert(1 <= jj && jj <= n);
1.367 +#if 1 /* 02/XI-2009 */
1.368 + if (degrad < 1e-6 * (1.0 + 0.001 * fabs(mip->obj_val)))
1.369 + { jj = branch_mostf(T, &next);
1.370 + goto done;
1.371 + }
1.372 +#endif
1.373 + if (T->parm->msg_lev >= GLP_MSG_DBG)
1.374 + { xprintf("branch_drtom: column %d chosen to branch on\n", jj);
1.375 + if (fabs(dd_dn) == DBL_MAX)
1.376 + xprintf("branch_drtom: down-branch is infeasible\n");
1.377 + else
1.378 + xprintf("branch_drtom: down-branch bound is %.9e\n",
1.379 + lpx_get_obj_val(mip) + dd_dn);
1.380 + if (fabs(dd_up) == DBL_MAX)
1.381 + xprintf("branch_drtom: up-branch is infeasible\n");
1.382 + else
1.383 + xprintf("branch_drtom: up-branch bound is %.9e\n",
1.384 + lpx_get_obj_val(mip) + dd_up);
1.385 + }
1.386 +done: *_next = next;
1.387 + return jj;
1.388 +}
1.389 +
1.390 +/**********************************************************************/
1.391 +
1.392 +struct csa
1.393 +{ /* common storage area */
1.394 + int *dn_cnt; /* int dn_cnt[1+n]; */
1.395 + /* dn_cnt[j] is the number of subproblems, whose LP relaxations
1.396 + have been solved and which are down-branches for variable x[j];
1.397 + dn_cnt[j] = 0 means the down pseudocost is uninitialized */
1.398 + double *dn_sum; /* double dn_sum[1+n]; */
1.399 + /* dn_sum[j] is the sum of per unit degradations of the objective
1.400 + over all dn_cnt[j] subproblems */
1.401 + int *up_cnt; /* int up_cnt[1+n]; */
1.402 + /* up_cnt[j] is the number of subproblems, whose LP relaxations
1.403 + have been solved and which are up-branches for variable x[j];
1.404 + up_cnt[j] = 0 means the up pseudocost is uninitialized */
1.405 + double *up_sum; /* double up_sum[1+n]; */
1.406 + /* up_sum[j] is the sum of per unit degradations of the objective
1.407 + over all up_cnt[j] subproblems */
1.408 +};
1.409 +
1.410 +void *ios_pcost_init(glp_tree *tree)
1.411 +{ /* initialize working data used on pseudocost branching */
1.412 + struct csa *csa;
1.413 + int n = tree->n, j;
1.414 + csa = xmalloc(sizeof(struct csa));
1.415 + csa->dn_cnt = xcalloc(1+n, sizeof(int));
1.416 + csa->dn_sum = xcalloc(1+n, sizeof(double));
1.417 + csa->up_cnt = xcalloc(1+n, sizeof(int));
1.418 + csa->up_sum = xcalloc(1+n, sizeof(double));
1.419 + for (j = 1; j <= n; j++)
1.420 + { csa->dn_cnt[j] = csa->up_cnt[j] = 0;
1.421 + csa->dn_sum[j] = csa->up_sum[j] = 0.0;
1.422 + }
1.423 + return csa;
1.424 +}
1.425 +
1.426 +static double eval_degrad(glp_prob *P, int j, double bnd)
1.427 +{ /* compute degradation of the objective on fixing x[j] at given
1.428 + value with a limited number of dual simplex iterations */
1.429 + /* this routine fixes column x[j] at specified value bnd,
1.430 + solves resulting LP, and returns a lower bound to degradation
1.431 + of the objective, degrad >= 0 */
1.432 + glp_prob *lp;
1.433 + glp_smcp parm;
1.434 + int ret;
1.435 + double degrad;
1.436 + /* the current basis must be optimal */
1.437 + xassert(glp_get_status(P) == GLP_OPT);
1.438 + /* create a copy of P */
1.439 + lp = glp_create_prob();
1.440 + glp_copy_prob(lp, P, 0);
1.441 + /* fix column x[j] at specified value */
1.442 + glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd);
1.443 + /* try to solve resulting LP */
1.444 + glp_init_smcp(&parm);
1.445 + parm.msg_lev = GLP_MSG_OFF;
1.446 + parm.meth = GLP_DUAL;
1.447 + parm.it_lim = 30;
1.448 + parm.out_dly = 1000;
1.449 + parm.meth = GLP_DUAL;
1.450 + ret = glp_simplex(lp, &parm);
1.451 + if (ret == 0 || ret == GLP_EITLIM)
1.452 + { if (glp_get_prim_stat(lp) == GLP_NOFEAS)
1.453 + { /* resulting LP has no primal feasible solution */
1.454 + degrad = DBL_MAX;
1.455 + }
1.456 + else if (glp_get_dual_stat(lp) == GLP_FEAS)
1.457 + { /* resulting basis is optimal or at least dual feasible,
1.458 + so we have the correct lower bound to degradation */
1.459 + if (P->dir == GLP_MIN)
1.460 + degrad = lp->obj_val - P->obj_val;
1.461 + else if (P->dir == GLP_MAX)
1.462 + degrad = P->obj_val - lp->obj_val;
1.463 + else
1.464 + xassert(P != P);
1.465 + /* degradation cannot be negative by definition */
1.466 + /* note that the lower bound to degradation may be close
1.467 + to zero even if its exact value is zero due to round-off
1.468 + errors on computing the objective value */
1.469 + if (degrad < 1e-6 * (1.0 + 0.001 * fabs(P->obj_val)))
1.470 + degrad = 0.0;
1.471 + }
1.472 + else
1.473 + { /* the final basis reported by the simplex solver is dual
1.474 + infeasible, so we cannot determine a non-trivial lower
1.475 + bound to degradation */
1.476 + degrad = 0.0;
1.477 + }
1.478 + }
1.479 + else
1.480 + { /* the simplex solver failed */
1.481 + degrad = 0.0;
1.482 + }
1.483 + /* delete the copy of P */
1.484 + glp_delete_prob(lp);
1.485 + return degrad;
1.486 +}
1.487 +
1.488 +void ios_pcost_update(glp_tree *tree)
1.489 +{ /* update history information for pseudocost branching */
1.490 + /* this routine is called every time when LP relaxation of the
1.491 + current subproblem has been solved to optimality with all lazy
1.492 + and cutting plane constraints included */
1.493 + int j;
1.494 + double dx, dz, psi;
1.495 + struct csa *csa = tree->pcost;
1.496 + xassert(csa != NULL);
1.497 + xassert(tree->curr != NULL);
1.498 + /* if the current subproblem is the root, skip updating */
1.499 + if (tree->curr->up == NULL) goto skip;
1.500 + /* determine branching variable x[j], which was used in the
1.501 + parent subproblem to create the current subproblem */
1.502 + j = tree->curr->up->br_var;
1.503 + xassert(1 <= j && j <= tree->n);
1.504 + /* determine the change dx[j] = new x[j] - old x[j],
1.505 + where new x[j] is a value of x[j] in optimal solution to LP
1.506 + relaxation of the current subproblem, old x[j] is a value of
1.507 + x[j] in optimal solution to LP relaxation of the parent
1.508 + subproblem */
1.509 + dx = tree->mip->col[j]->prim - tree->curr->up->br_val;
1.510 + xassert(dx != 0.0);
1.511 + /* determine corresponding change dz = new dz - old dz in the
1.512 + objective function value */
1.513 + dz = tree->mip->obj_val - tree->curr->up->lp_obj;
1.514 + /* determine per unit degradation of the objective function */
1.515 + psi = fabs(dz / dx);
1.516 + /* update history information */
1.517 + if (dx < 0.0)
1.518 + { /* the current subproblem is down-branch */
1.519 + csa->dn_cnt[j]++;
1.520 + csa->dn_sum[j] += psi;
1.521 + }
1.522 + else /* dx > 0.0 */
1.523 + { /* the current subproblem is up-branch */
1.524 + csa->up_cnt[j]++;
1.525 + csa->up_sum[j] += psi;
1.526 + }
1.527 +skip: return;
1.528 +}
1.529 +
1.530 +void ios_pcost_free(glp_tree *tree)
1.531 +{ /* free working area used on pseudocost branching */
1.532 + struct csa *csa = tree->pcost;
1.533 + xassert(csa != NULL);
1.534 + xfree(csa->dn_cnt);
1.535 + xfree(csa->dn_sum);
1.536 + xfree(csa->up_cnt);
1.537 + xfree(csa->up_sum);
1.538 + xfree(csa);
1.539 + tree->pcost = NULL;
1.540 + return;
1.541 +}
1.542 +
1.543 +static double eval_psi(glp_tree *T, int j, int brnch)
1.544 +{ /* compute estimation of pseudocost of variable x[j] for down-
1.545 + or up-branch */
1.546 + struct csa *csa = T->pcost;
1.547 + double beta, degrad, psi;
1.548 + xassert(csa != NULL);
1.549 + xassert(1 <= j && j <= T->n);
1.550 + if (brnch == GLP_DN_BRNCH)
1.551 + { /* down-branch */
1.552 + if (csa->dn_cnt[j] == 0)
1.553 + { /* initialize down pseudocost */
1.554 + beta = T->mip->col[j]->prim;
1.555 + degrad = eval_degrad(T->mip, j, floor(beta));
1.556 + if (degrad == DBL_MAX)
1.557 + { psi = DBL_MAX;
1.558 + goto done;
1.559 + }
1.560 + csa->dn_cnt[j] = 1;
1.561 + csa->dn_sum[j] = degrad / (beta - floor(beta));
1.562 + }
1.563 + psi = csa->dn_sum[j] / (double)csa->dn_cnt[j];
1.564 + }
1.565 + else if (brnch == GLP_UP_BRNCH)
1.566 + { /* up-branch */
1.567 + if (csa->up_cnt[j] == 0)
1.568 + { /* initialize up pseudocost */
1.569 + beta = T->mip->col[j]->prim;
1.570 + degrad = eval_degrad(T->mip, j, ceil(beta));
1.571 + if (degrad == DBL_MAX)
1.572 + { psi = DBL_MAX;
1.573 + goto done;
1.574 + }
1.575 + csa->up_cnt[j] = 1;
1.576 + csa->up_sum[j] = degrad / (ceil(beta) - beta);
1.577 + }
1.578 + psi = csa->up_sum[j] / (double)csa->up_cnt[j];
1.579 + }
1.580 + else
1.581 + xassert(brnch != brnch);
1.582 +done: return psi;
1.583 +}
1.584 +
1.585 +static void progress(glp_tree *T)
1.586 +{ /* display progress of pseudocost initialization */
1.587 + struct csa *csa = T->pcost;
1.588 + int j, nv = 0, ni = 0;
1.589 + for (j = 1; j <= T->n; j++)
1.590 + { if (glp_ios_can_branch(T, j))
1.591 + { nv++;
1.592 + if (csa->dn_cnt[j] > 0 && csa->up_cnt[j] > 0) ni++;
1.593 + }
1.594 + }
1.595 + xprintf("Pseudocosts initialized for %d of %d variables\n",
1.596 + ni, nv);
1.597 + return;
1.598 +}
1.599 +
1.600 +int ios_pcost_branch(glp_tree *T, int *_next)
1.601 +{ /* choose branching variable with pseudocost branching */
1.602 + glp_long t = xtime();
1.603 + int j, jjj, sel;
1.604 + double beta, psi, d1, d2, d, dmax;
1.605 + /* initialize the working arrays */
1.606 + if (T->pcost == NULL)
1.607 + T->pcost = ios_pcost_init(T);
1.608 + /* nothing has been chosen so far */
1.609 + jjj = 0, dmax = -1.0;
1.610 + /* go through the list of branching candidates */
1.611 + for (j = 1; j <= T->n; j++)
1.612 + { if (!glp_ios_can_branch(T, j)) continue;
1.613 + /* determine primal value of x[j] in optimal solution to LP
1.614 + relaxation of the current subproblem */
1.615 + beta = T->mip->col[j]->prim;
1.616 + /* estimate pseudocost of x[j] for down-branch */
1.617 + psi = eval_psi(T, j, GLP_DN_BRNCH);
1.618 + if (psi == DBL_MAX)
1.619 + { /* down-branch has no primal feasible solution */
1.620 + jjj = j, sel = GLP_DN_BRNCH;
1.621 + goto done;
1.622 + }
1.623 + /* estimate degradation of the objective for down-branch */
1.624 + d1 = psi * (beta - floor(beta));
1.625 + /* estimate pseudocost of x[j] for up-branch */
1.626 + psi = eval_psi(T, j, GLP_UP_BRNCH);
1.627 + if (psi == DBL_MAX)
1.628 + { /* up-branch has no primal feasible solution */
1.629 + jjj = j, sel = GLP_UP_BRNCH;
1.630 + goto done;
1.631 + }
1.632 + /* estimate degradation of the objective for up-branch */
1.633 + d2 = psi * (ceil(beta) - beta);
1.634 + /* determine d = max(d1, d2) */
1.635 + d = (d1 > d2 ? d1 : d2);
1.636 + /* choose x[j] which provides maximal estimated degradation of
1.637 + the objective either in down- or up-branch */
1.638 + if (dmax < d)
1.639 + { dmax = d;
1.640 + jjj = j;
1.641 + /* continue the search from a subproblem, where degradation
1.642 + is less than in other one */
1.643 + sel = (d1 <= d2 ? GLP_DN_BRNCH : GLP_UP_BRNCH);
1.644 + }
1.645 + /* display progress of pseudocost initialization */
1.646 + if (T->parm->msg_lev >= GLP_ON)
1.647 + { if (xdifftime(xtime(), t) >= 10.0)
1.648 + { progress(T);
1.649 + t = xtime();
1.650 + }
1.651 + }
1.652 + }
1.653 + if (dmax == 0.0)
1.654 + { /* no degradation is indicated; choose a variable having most
1.655 + fractional value */
1.656 + jjj = branch_mostf(T, &sel);
1.657 + }
1.658 +done: *_next = sel;
1.659 + return jjj;
1.660 +}
1.661 +
1.662 +/* eof */