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