1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpios09.c	Sun Nov 06 20:59:10 2011 +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, 2011 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 */