alpar@1: /* glpios09.c (branching heuristics) */ alpar@1: alpar@1: /*********************************************************************** alpar@1: * This code is part of GLPK (GNU Linear Programming Kit). alpar@1: * alpar@1: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@1: * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, alpar@1: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@1: * E-mail: . alpar@1: * alpar@1: * GLPK is free software: you can redistribute it and/or modify it alpar@1: * under the terms of the GNU General Public License as published by alpar@1: * the Free Software Foundation, either version 3 of the License, or alpar@1: * (at your option) any later version. alpar@1: * alpar@1: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@1: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@1: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@1: * License for more details. alpar@1: * alpar@1: * You should have received a copy of the GNU General Public License alpar@1: * along with GLPK. If not, see . alpar@1: ***********************************************************************/ alpar@1: alpar@1: #include "glpios.h" alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * ios_choose_var - select variable to branch on alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpios.h" alpar@1: * int ios_choose_var(glp_tree *T, int *next); alpar@1: * alpar@1: * The routine ios_choose_var chooses a variable from the candidate alpar@1: * list to branch on. Additionally the routine provides a flag stored alpar@1: * in the location next to suggests which of the child subproblems alpar@1: * should be solved next. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine ios_choose_var returns the ordinal number of the column alpar@1: * choosen. */ alpar@1: alpar@1: static int branch_first(glp_tree *T, int *next); alpar@1: static int branch_last(glp_tree *T, int *next); alpar@1: static int branch_mostf(glp_tree *T, int *next); alpar@1: static int branch_drtom(glp_tree *T, int *next); alpar@1: alpar@1: int ios_choose_var(glp_tree *T, int *next) alpar@1: { int j; alpar@1: if (T->parm->br_tech == GLP_BR_FFV) alpar@1: { /* branch on first fractional variable */ alpar@1: j = branch_first(T, next); alpar@1: } alpar@1: else if (T->parm->br_tech == GLP_BR_LFV) alpar@1: { /* branch on last fractional variable */ alpar@1: j = branch_last(T, next); alpar@1: } alpar@1: else if (T->parm->br_tech == GLP_BR_MFV) alpar@1: { /* branch on most fractional variable */ alpar@1: j = branch_mostf(T, next); alpar@1: } alpar@1: else if (T->parm->br_tech == GLP_BR_DTH) alpar@1: { /* branch using the heuristic by Dreebeck and Tomlin */ alpar@1: j = branch_drtom(T, next); alpar@1: } alpar@1: else if (T->parm->br_tech == GLP_BR_PCH) alpar@1: { /* hybrid pseudocost heuristic */ alpar@1: j = ios_pcost_branch(T, next); alpar@1: } alpar@1: else alpar@1: xassert(T != T); alpar@1: return j; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * branch_first - choose first branching variable alpar@1: * alpar@1: * This routine looks up the list of structural variables and chooses alpar@1: * the first one, which is of integer kind and has fractional value in alpar@1: * optimal solution to the current LP relaxation. alpar@1: * alpar@1: * This routine also selects the branch to be solved next where integer alpar@1: * infeasibility of the chosen variable is less than in other one. */ alpar@1: alpar@1: static int branch_first(glp_tree *T, int *_next) alpar@1: { int j, next; alpar@1: double beta; alpar@1: /* choose the column to branch on */ alpar@1: for (j = 1; j <= T->n; j++) alpar@1: if (T->non_int[j]) break; alpar@1: xassert(1 <= j && j <= T->n); alpar@1: /* select the branch to be solved next */ alpar@1: beta = glp_get_col_prim(T->mip, j); alpar@1: if (beta - floor(beta) < ceil(beta) - beta) alpar@1: next = GLP_DN_BRNCH; alpar@1: else alpar@1: next = GLP_UP_BRNCH; alpar@1: *_next = next; alpar@1: return j; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * branch_last - choose last branching variable alpar@1: * alpar@1: * This routine looks up the list of structural variables and chooses alpar@1: * the last one, which is of integer kind and has fractional value in alpar@1: * optimal solution to the current LP relaxation. alpar@1: * alpar@1: * This routine also selects the branch to be solved next where integer alpar@1: * infeasibility of the chosen variable is less than in other one. */ alpar@1: alpar@1: static int branch_last(glp_tree *T, int *_next) alpar@1: { int j, next; alpar@1: double beta; alpar@1: /* choose the column to branch on */ alpar@1: for (j = T->n; j >= 1; j--) alpar@1: if (T->non_int[j]) break; alpar@1: xassert(1 <= j && j <= T->n); alpar@1: /* select the branch to be solved next */ alpar@1: beta = glp_get_col_prim(T->mip, j); alpar@1: if (beta - floor(beta) < ceil(beta) - beta) alpar@1: next = GLP_DN_BRNCH; alpar@1: else alpar@1: next = GLP_UP_BRNCH; alpar@1: *_next = next; alpar@1: return j; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * branch_mostf - choose most fractional branching variable alpar@1: * alpar@1: * This routine looks up the list of structural variables and chooses alpar@1: * that one, which is of integer kind and has most fractional value in alpar@1: * optimal solution to the current LP relaxation. alpar@1: * alpar@1: * This routine also selects the branch to be solved next where integer alpar@1: * infeasibility of the chosen variable is less than in other one. alpar@1: * alpar@1: * (Alexander Martin notices that "...most infeasible is as good as alpar@1: * random...".) */ alpar@1: alpar@1: static int branch_mostf(glp_tree *T, int *_next) alpar@1: { int j, jj, next; alpar@1: double beta, most, temp; alpar@1: /* choose the column to branch on */ alpar@1: jj = 0, most = DBL_MAX; alpar@1: for (j = 1; j <= T->n; j++) alpar@1: { if (T->non_int[j]) alpar@1: { beta = glp_get_col_prim(T->mip, j); alpar@1: temp = floor(beta) + 0.5; alpar@1: if (most > fabs(beta - temp)) alpar@1: { jj = j, most = fabs(beta - temp); alpar@1: if (beta < temp) alpar@1: next = GLP_DN_BRNCH; alpar@1: else alpar@1: next = GLP_UP_BRNCH; alpar@1: } alpar@1: } alpar@1: } alpar@1: *_next = next; alpar@1: return jj; alpar@1: } alpar@1: alpar@1: /*********************************************************************** alpar@1: * branch_drtom - choose branching var using Driebeck-Tomlin heuristic alpar@1: * alpar@1: * This routine chooses a structural variable, which is required to be alpar@1: * integral and has fractional value in optimal solution of the current alpar@1: * LP relaxation, using a heuristic proposed by Driebeck and Tomlin. alpar@1: * alpar@1: * The routine also selects the branch to be solved next, again due to alpar@1: * Driebeck and Tomlin. alpar@1: * alpar@1: * This routine is based on the heuristic proposed in: alpar@1: * alpar@1: * Driebeck N.J. An algorithm for the solution of mixed-integer alpar@1: * programming problems, Management Science, 12: 576-87 (1966); alpar@1: * alpar@1: * and improved in: alpar@1: * alpar@1: * Tomlin J.A. Branch and bound methods for integer and non-convex alpar@1: * programming, in J.Abadie (ed.), Integer and Nonlinear Programming, alpar@1: * North-Holland, Amsterdam, pp. 437-50 (1970). alpar@1: * alpar@1: * Must note that this heuristic is time-expensive, because computing alpar@1: * one-step degradation (see the routine below) requires one BTRAN for alpar@1: * each fractional-valued structural variable. */ alpar@1: alpar@1: static int branch_drtom(glp_tree *T, int *_next) alpar@1: { glp_prob *mip = T->mip; alpar@1: int m = mip->m; alpar@1: int n = mip->n; alpar@1: char *non_int = T->non_int; alpar@1: int j, jj, k, t, next, kase, len, stat, *ind; alpar@1: double x, dk, alfa, delta_j, delta_k, delta_z, dz_dn, dz_up, alpar@1: dd_dn, dd_up, degrad, *val; alpar@1: /* basic solution of LP relaxation must be optimal */ alpar@1: xassert(glp_get_status(mip) == GLP_OPT); alpar@1: /* allocate working arrays */ alpar@1: ind = xcalloc(1+n, sizeof(int)); alpar@1: val = xcalloc(1+n, sizeof(double)); alpar@1: /* nothing has been chosen so far */ alpar@1: jj = 0, degrad = -1.0; alpar@1: /* walk through the list of columns (structural variables) */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { /* if j-th column is not marked as fractional, skip it */ alpar@1: if (!non_int[j]) continue; alpar@1: /* obtain (fractional) value of j-th column in basic solution alpar@1: of LP relaxation */ alpar@1: x = glp_get_col_prim(mip, j); alpar@1: /* since the value of j-th column is fractional, the column is alpar@1: basic; compute corresponding row of the simplex table */ alpar@1: len = glp_eval_tab_row(mip, m+j, ind, val); alpar@1: /* the following fragment computes a change in the objective alpar@1: function: delta Z = new Z - old Z, where old Z is the alpar@1: objective value in the current optimal basis, and new Z is alpar@1: the objective value in the adjacent basis, for two cases: alpar@1: 1) if new upper bound ub' = floor(x[j]) is introduced for alpar@1: j-th column (down branch); alpar@1: 2) if new lower bound lb' = ceil(x[j]) is introduced for alpar@1: j-th column (up branch); alpar@1: since in both cases the solution remaining dual feasible alpar@1: becomes primal infeasible, one implicit simplex iteration alpar@1: is performed to determine the change delta Z; alpar@1: it is obvious that new Z, which is never better than old Z, alpar@1: is a lower (minimization) or upper (maximization) bound of alpar@1: the objective function for down- and up-branches. */ alpar@1: for (kase = -1; kase <= +1; kase += 2) alpar@1: { /* if kase < 0, the new upper bound of x[j] is introduced; alpar@1: in this case x[j] should decrease in order to leave the alpar@1: basis and go to its new upper bound */ alpar@1: /* if kase > 0, the new lower bound of x[j] is introduced; alpar@1: in this case x[j] should increase in order to leave the alpar@1: basis and go to its new lower bound */ alpar@1: /* apply the dual ratio test in order to determine which alpar@1: auxiliary or structural variable should enter the basis alpar@1: to keep dual feasibility */ alpar@1: k = glp_dual_rtest(mip, len, ind, val, kase, 1e-9); alpar@1: if (k != 0) k = ind[k]; alpar@1: /* if no non-basic variable has been chosen, LP relaxation alpar@1: of corresponding branch being primal infeasible and dual alpar@1: unbounded has no primal feasible solution; in this case alpar@1: the change delta Z is formally set to infinity */ alpar@1: if (k == 0) alpar@1: { delta_z = alpar@1: (T->mip->dir == GLP_MIN ? +DBL_MAX : -DBL_MAX); alpar@1: goto skip; alpar@1: } alpar@1: /* row of the simplex table that corresponds to non-basic alpar@1: variable x[k] choosen by the dual ratio test is: alpar@1: x[j] = ... + alfa * x[k] + ... alpar@1: where alfa is the influence coefficient (an element of alpar@1: the simplex table row) */ alpar@1: /* determine the coefficient alfa */ alpar@1: for (t = 1; t <= len; t++) if (ind[t] == k) break; alpar@1: xassert(1 <= t && t <= len); alpar@1: alfa = val[t]; alpar@1: /* since in the adjacent basis the variable x[j] becomes alpar@1: non-basic, knowing its value in the current basis we can alpar@1: determine its change delta x[j] = new x[j] - old x[j] */ alpar@1: delta_j = (kase < 0 ? floor(x) : ceil(x)) - x; alpar@1: /* and knowing the coefficient alfa we can determine the alpar@1: corresponding change delta x[k] = new x[k] - old x[k], alpar@1: where old x[k] is a value of x[k] in the current basis, alpar@1: and new x[k] is a value of x[k] in the adjacent basis */ alpar@1: delta_k = delta_j / alfa; alpar@1: /* Tomlin noticed that if the variable x[k] is of integer alpar@1: kind, its change cannot be less (eventually) than one in alpar@1: the magnitude */ alpar@1: if (k > m && glp_get_col_kind(mip, k-m) != GLP_CV) alpar@1: { /* x[k] is structural integer variable */ alpar@1: if (fabs(delta_k - floor(delta_k + 0.5)) > 1e-3) alpar@1: { if (delta_k > 0.0) alpar@1: delta_k = ceil(delta_k); /* +3.14 -> +4 */ alpar@1: else alpar@1: delta_k = floor(delta_k); /* -3.14 -> -4 */ alpar@1: } alpar@1: } alpar@1: /* now determine the status and reduced cost of x[k] in the alpar@1: current basis */ alpar@1: if (k <= m) alpar@1: { stat = glp_get_row_stat(mip, k); alpar@1: dk = glp_get_row_dual(mip, k); alpar@1: } alpar@1: else alpar@1: { stat = glp_get_col_stat(mip, k-m); alpar@1: dk = glp_get_col_dual(mip, k-m); alpar@1: } alpar@1: /* if the current basis is dual degenerate, some reduced alpar@1: costs which are close to zero may have wrong sign due to alpar@1: round-off errors, so correct the sign of d[k] */ alpar@1: switch (T->mip->dir) alpar@1: { case GLP_MIN: alpar@1: if (stat == GLP_NL && dk < 0.0 || alpar@1: stat == GLP_NU && dk > 0.0 || alpar@1: stat == GLP_NF) dk = 0.0; alpar@1: break; alpar@1: case GLP_MAX: alpar@1: if (stat == GLP_NL && dk > 0.0 || alpar@1: stat == GLP_NU && dk < 0.0 || alpar@1: stat == GLP_NF) dk = 0.0; alpar@1: break; alpar@1: default: alpar@1: xassert(T != T); alpar@1: } alpar@1: /* now knowing the change of x[k] and its reduced cost d[k] alpar@1: we can compute the corresponding change in the objective alpar@1: function delta Z = new Z - old Z = d[k] * delta x[k]; alpar@1: note that due to Tomlin's modification new Z can be even alpar@1: worse than in the adjacent basis */ alpar@1: delta_z = dk * delta_k; alpar@1: skip: /* new Z is never better than old Z, therefore the change alpar@1: delta Z is always non-negative (in case of minimization) alpar@1: or non-positive (in case of maximization) */ alpar@1: switch (T->mip->dir) alpar@1: { case GLP_MIN: xassert(delta_z >= 0.0); break; alpar@1: case GLP_MAX: xassert(delta_z <= 0.0); break; alpar@1: default: xassert(T != T); alpar@1: } alpar@1: /* save the change in the objective fnction for down- and alpar@1: up-branches, respectively */ alpar@1: if (kase < 0) dz_dn = delta_z; else dz_up = delta_z; alpar@1: } alpar@1: /* thus, in down-branch no integer feasible solution can be alpar@1: better than Z + dz_dn, and in up-branch no integer feasible alpar@1: solution can be better than Z + dz_up, where Z is value of alpar@1: the objective function in the current basis */ alpar@1: /* following the heuristic by Driebeck and Tomlin we choose a alpar@1: column (i.e. structural variable) which provides largest alpar@1: degradation of the objective function in some of branches; alpar@1: besides, we select the branch with smaller degradation to alpar@1: be solved next and keep other branch with larger degradation alpar@1: in the active list hoping to minimize the number of further alpar@1: backtrackings */ alpar@1: if (degrad < fabs(dz_dn) || degrad < fabs(dz_up)) alpar@1: { jj = j; alpar@1: if (fabs(dz_dn) < fabs(dz_up)) alpar@1: { /* select down branch to be solved next */ alpar@1: next = GLP_DN_BRNCH; alpar@1: degrad = fabs(dz_up); alpar@1: } alpar@1: else alpar@1: { /* select up branch to be solved next */ alpar@1: next = GLP_UP_BRNCH; alpar@1: degrad = fabs(dz_dn); alpar@1: } alpar@1: /* save the objective changes for printing */ alpar@1: dd_dn = dz_dn, dd_up = dz_up; alpar@1: /* if down- or up-branch has no feasible solution, we does alpar@1: not need to consider other candidates (in principle, the alpar@1: corresponding branch could be pruned right now) */ alpar@1: if (degrad == DBL_MAX) break; alpar@1: } alpar@1: } alpar@1: /* free working arrays */ alpar@1: xfree(ind); alpar@1: xfree(val); alpar@1: /* something must be chosen */ alpar@1: xassert(1 <= jj && jj <= n); alpar@1: #if 1 /* 02/XI-2009 */ alpar@1: if (degrad < 1e-6 * (1.0 + 0.001 * fabs(mip->obj_val))) alpar@1: { jj = branch_mostf(T, &next); alpar@1: goto done; alpar@1: } alpar@1: #endif alpar@1: if (T->parm->msg_lev >= GLP_MSG_DBG) alpar@1: { xprintf("branch_drtom: column %d chosen to branch on\n", jj); alpar@1: if (fabs(dd_dn) == DBL_MAX) alpar@1: xprintf("branch_drtom: down-branch is infeasible\n"); alpar@1: else alpar@1: xprintf("branch_drtom: down-branch bound is %.9e\n", alpar@1: lpx_get_obj_val(mip) + dd_dn); alpar@1: if (fabs(dd_up) == DBL_MAX) alpar@1: xprintf("branch_drtom: up-branch is infeasible\n"); alpar@1: else alpar@1: xprintf("branch_drtom: up-branch bound is %.9e\n", alpar@1: lpx_get_obj_val(mip) + dd_up); alpar@1: } alpar@1: done: *_next = next; alpar@1: return jj; alpar@1: } alpar@1: alpar@1: /**********************************************************************/ alpar@1: alpar@1: struct csa alpar@1: { /* common storage area */ alpar@1: int *dn_cnt; /* int dn_cnt[1+n]; */ alpar@1: /* dn_cnt[j] is the number of subproblems, whose LP relaxations alpar@1: have been solved and which are down-branches for variable x[j]; alpar@1: dn_cnt[j] = 0 means the down pseudocost is uninitialized */ alpar@1: double *dn_sum; /* double dn_sum[1+n]; */ alpar@1: /* dn_sum[j] is the sum of per unit degradations of the objective alpar@1: over all dn_cnt[j] subproblems */ alpar@1: int *up_cnt; /* int up_cnt[1+n]; */ alpar@1: /* up_cnt[j] is the number of subproblems, whose LP relaxations alpar@1: have been solved and which are up-branches for variable x[j]; alpar@1: up_cnt[j] = 0 means the up pseudocost is uninitialized */ alpar@1: double *up_sum; /* double up_sum[1+n]; */ alpar@1: /* up_sum[j] is the sum of per unit degradations of the objective alpar@1: over all up_cnt[j] subproblems */ alpar@1: }; alpar@1: alpar@1: void *ios_pcost_init(glp_tree *tree) alpar@1: { /* initialize working data used on pseudocost branching */ alpar@1: struct csa *csa; alpar@1: int n = tree->n, j; alpar@1: csa = xmalloc(sizeof(struct csa)); alpar@1: csa->dn_cnt = xcalloc(1+n, sizeof(int)); alpar@1: csa->dn_sum = xcalloc(1+n, sizeof(double)); alpar@1: csa->up_cnt = xcalloc(1+n, sizeof(int)); alpar@1: csa->up_sum = xcalloc(1+n, sizeof(double)); alpar@1: for (j = 1; j <= n; j++) alpar@1: { csa->dn_cnt[j] = csa->up_cnt[j] = 0; alpar@1: csa->dn_sum[j] = csa->up_sum[j] = 0.0; alpar@1: } alpar@1: return csa; alpar@1: } alpar@1: alpar@1: static double eval_degrad(glp_prob *P, int j, double bnd) alpar@1: { /* compute degradation of the objective on fixing x[j] at given alpar@1: value with a limited number of dual simplex iterations */ alpar@1: /* this routine fixes column x[j] at specified value bnd, alpar@1: solves resulting LP, and returns a lower bound to degradation alpar@1: of the objective, degrad >= 0 */ alpar@1: glp_prob *lp; alpar@1: glp_smcp parm; alpar@1: int ret; alpar@1: double degrad; alpar@1: /* the current basis must be optimal */ alpar@1: xassert(glp_get_status(P) == GLP_OPT); alpar@1: /* create a copy of P */ alpar@1: lp = glp_create_prob(); alpar@1: glp_copy_prob(lp, P, 0); alpar@1: /* fix column x[j] at specified value */ alpar@1: glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd); alpar@1: /* try to solve resulting LP */ alpar@1: glp_init_smcp(&parm); alpar@1: parm.msg_lev = GLP_MSG_OFF; alpar@1: parm.meth = GLP_DUAL; alpar@1: parm.it_lim = 30; alpar@1: parm.out_dly = 1000; alpar@1: parm.meth = GLP_DUAL; alpar@1: ret = glp_simplex(lp, &parm); alpar@1: if (ret == 0 || ret == GLP_EITLIM) alpar@1: { if (glp_get_prim_stat(lp) == GLP_NOFEAS) alpar@1: { /* resulting LP has no primal feasible solution */ alpar@1: degrad = DBL_MAX; alpar@1: } alpar@1: else if (glp_get_dual_stat(lp) == GLP_FEAS) alpar@1: { /* resulting basis is optimal or at least dual feasible, alpar@1: so we have the correct lower bound to degradation */ alpar@1: if (P->dir == GLP_MIN) alpar@1: degrad = lp->obj_val - P->obj_val; alpar@1: else if (P->dir == GLP_MAX) alpar@1: degrad = P->obj_val - lp->obj_val; alpar@1: else alpar@1: xassert(P != P); alpar@1: /* degradation cannot be negative by definition */ alpar@1: /* note that the lower bound to degradation may be close alpar@1: to zero even if its exact value is zero due to round-off alpar@1: errors on computing the objective value */ alpar@1: if (degrad < 1e-6 * (1.0 + 0.001 * fabs(P->obj_val))) alpar@1: degrad = 0.0; alpar@1: } alpar@1: else alpar@1: { /* the final basis reported by the simplex solver is dual alpar@1: infeasible, so we cannot determine a non-trivial lower alpar@1: bound to degradation */ alpar@1: degrad = 0.0; alpar@1: } alpar@1: } alpar@1: else alpar@1: { /* the simplex solver failed */ alpar@1: degrad = 0.0; alpar@1: } alpar@1: /* delete the copy of P */ alpar@1: glp_delete_prob(lp); alpar@1: return degrad; alpar@1: } alpar@1: alpar@1: void ios_pcost_update(glp_tree *tree) alpar@1: { /* update history information for pseudocost branching */ alpar@1: /* this routine is called every time when LP relaxation of the alpar@1: current subproblem has been solved to optimality with all lazy alpar@1: and cutting plane constraints included */ alpar@1: int j; alpar@1: double dx, dz, psi; alpar@1: struct csa *csa = tree->pcost; alpar@1: xassert(csa != NULL); alpar@1: xassert(tree->curr != NULL); alpar@1: /* if the current subproblem is the root, skip updating */ alpar@1: if (tree->curr->up == NULL) goto skip; alpar@1: /* determine branching variable x[j], which was used in the alpar@1: parent subproblem to create the current subproblem */ alpar@1: j = tree->curr->up->br_var; alpar@1: xassert(1 <= j && j <= tree->n); alpar@1: /* determine the change dx[j] = new x[j] - old x[j], alpar@1: where new x[j] is a value of x[j] in optimal solution to LP alpar@1: relaxation of the current subproblem, old x[j] is a value of alpar@1: x[j] in optimal solution to LP relaxation of the parent alpar@1: subproblem */ alpar@1: dx = tree->mip->col[j]->prim - tree->curr->up->br_val; alpar@1: xassert(dx != 0.0); alpar@1: /* determine corresponding change dz = new dz - old dz in the alpar@1: objective function value */ alpar@1: dz = tree->mip->obj_val - tree->curr->up->lp_obj; alpar@1: /* determine per unit degradation of the objective function */ alpar@1: psi = fabs(dz / dx); alpar@1: /* update history information */ alpar@1: if (dx < 0.0) alpar@1: { /* the current subproblem is down-branch */ alpar@1: csa->dn_cnt[j]++; alpar@1: csa->dn_sum[j] += psi; alpar@1: } alpar@1: else /* dx > 0.0 */ alpar@1: { /* the current subproblem is up-branch */ alpar@1: csa->up_cnt[j]++; alpar@1: csa->up_sum[j] += psi; alpar@1: } alpar@1: skip: return; alpar@1: } alpar@1: alpar@1: void ios_pcost_free(glp_tree *tree) alpar@1: { /* free working area used on pseudocost branching */ alpar@1: struct csa *csa = tree->pcost; alpar@1: xassert(csa != NULL); alpar@1: xfree(csa->dn_cnt); alpar@1: xfree(csa->dn_sum); alpar@1: xfree(csa->up_cnt); alpar@1: xfree(csa->up_sum); alpar@1: xfree(csa); alpar@1: tree->pcost = NULL; alpar@1: return; alpar@1: } alpar@1: alpar@1: static double eval_psi(glp_tree *T, int j, int brnch) alpar@1: { /* compute estimation of pseudocost of variable x[j] for down- alpar@1: or up-branch */ alpar@1: struct csa *csa = T->pcost; alpar@1: double beta, degrad, psi; alpar@1: xassert(csa != NULL); alpar@1: xassert(1 <= j && j <= T->n); alpar@1: if (brnch == GLP_DN_BRNCH) alpar@1: { /* down-branch */ alpar@1: if (csa->dn_cnt[j] == 0) alpar@1: { /* initialize down pseudocost */ alpar@1: beta = T->mip->col[j]->prim; alpar@1: degrad = eval_degrad(T->mip, j, floor(beta)); alpar@1: if (degrad == DBL_MAX) alpar@1: { psi = DBL_MAX; alpar@1: goto done; alpar@1: } alpar@1: csa->dn_cnt[j] = 1; alpar@1: csa->dn_sum[j] = degrad / (beta - floor(beta)); alpar@1: } alpar@1: psi = csa->dn_sum[j] / (double)csa->dn_cnt[j]; alpar@1: } alpar@1: else if (brnch == GLP_UP_BRNCH) alpar@1: { /* up-branch */ alpar@1: if (csa->up_cnt[j] == 0) alpar@1: { /* initialize up pseudocost */ alpar@1: beta = T->mip->col[j]->prim; alpar@1: degrad = eval_degrad(T->mip, j, ceil(beta)); alpar@1: if (degrad == DBL_MAX) alpar@1: { psi = DBL_MAX; alpar@1: goto done; alpar@1: } alpar@1: csa->up_cnt[j] = 1; alpar@1: csa->up_sum[j] = degrad / (ceil(beta) - beta); alpar@1: } alpar@1: psi = csa->up_sum[j] / (double)csa->up_cnt[j]; alpar@1: } alpar@1: else alpar@1: xassert(brnch != brnch); alpar@1: done: return psi; alpar@1: } alpar@1: alpar@1: static void progress(glp_tree *T) alpar@1: { /* display progress of pseudocost initialization */ alpar@1: struct csa *csa = T->pcost; alpar@1: int j, nv = 0, ni = 0; alpar@1: for (j = 1; j <= T->n; j++) alpar@1: { if (glp_ios_can_branch(T, j)) alpar@1: { nv++; alpar@1: if (csa->dn_cnt[j] > 0 && csa->up_cnt[j] > 0) ni++; alpar@1: } alpar@1: } alpar@1: xprintf("Pseudocosts initialized for %d of %d variables\n", alpar@1: ni, nv); alpar@1: return; alpar@1: } alpar@1: alpar@1: int ios_pcost_branch(glp_tree *T, int *_next) alpar@1: { /* choose branching variable with pseudocost branching */ alpar@1: glp_long t = xtime(); alpar@1: int j, jjj, sel; alpar@1: double beta, psi, d1, d2, d, dmax; alpar@1: /* initialize the working arrays */ alpar@1: if (T->pcost == NULL) alpar@1: T->pcost = ios_pcost_init(T); alpar@1: /* nothing has been chosen so far */ alpar@1: jjj = 0, dmax = -1.0; alpar@1: /* go through the list of branching candidates */ alpar@1: for (j = 1; j <= T->n; j++) alpar@1: { if (!glp_ios_can_branch(T, j)) continue; alpar@1: /* determine primal value of x[j] in optimal solution to LP alpar@1: relaxation of the current subproblem */ alpar@1: beta = T->mip->col[j]->prim; alpar@1: /* estimate pseudocost of x[j] for down-branch */ alpar@1: psi = eval_psi(T, j, GLP_DN_BRNCH); alpar@1: if (psi == DBL_MAX) alpar@1: { /* down-branch has no primal feasible solution */ alpar@1: jjj = j, sel = GLP_DN_BRNCH; alpar@1: goto done; alpar@1: } alpar@1: /* estimate degradation of the objective for down-branch */ alpar@1: d1 = psi * (beta - floor(beta)); alpar@1: /* estimate pseudocost of x[j] for up-branch */ alpar@1: psi = eval_psi(T, j, GLP_UP_BRNCH); alpar@1: if (psi == DBL_MAX) alpar@1: { /* up-branch has no primal feasible solution */ alpar@1: jjj = j, sel = GLP_UP_BRNCH; alpar@1: goto done; alpar@1: } alpar@1: /* estimate degradation of the objective for up-branch */ alpar@1: d2 = psi * (ceil(beta) - beta); alpar@1: /* determine d = max(d1, d2) */ alpar@1: d = (d1 > d2 ? d1 : d2); alpar@1: /* choose x[j] which provides maximal estimated degradation of alpar@1: the objective either in down- or up-branch */ alpar@1: if (dmax < d) alpar@1: { dmax = d; alpar@1: jjj = j; alpar@1: /* continue the search from a subproblem, where degradation alpar@1: is less than in other one */ alpar@1: sel = (d1 <= d2 ? GLP_DN_BRNCH : GLP_UP_BRNCH); alpar@1: } alpar@1: /* display progress of pseudocost initialization */ alpar@1: if (T->parm->msg_lev >= GLP_ON) alpar@1: { if (xdifftime(xtime(), t) >= 10.0) alpar@1: { progress(T); alpar@1: t = xtime(); alpar@1: } alpar@1: } alpar@1: } alpar@1: if (dmax == 0.0) alpar@1: { /* no degradation is indicated; choose a variable having most alpar@1: fractional value */ alpar@1: jjj = branch_mostf(T, &sel); alpar@1: } alpar@1: done: *_next = sel; alpar@1: return jjj; alpar@1: } alpar@1: alpar@1: /* eof */