/* glpios09.c (branching heuristics) */

/***********************************************************************

*  This code is part of GLPK (GNU Linear Programming Kit).

*

*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,

*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,

*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.

*  E-mail: <mao@gnu.org>.

*

*  GLPK is free software: you can redistribute it and/or modify it

*  under the terms of the GNU General Public License as published by

*  the Free Software Foundation, either version 3 of the License, or

*  (at your option) any later version.

*

*  GLPK is distributed in the hope that it will be useful, but WITHOUT

*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

*  License for more details.

*

*  You should have received a copy of the GNU General Public License

*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.

***********************************************************************/

#include "glpios.h"

/***********************************************************************

*  NAME

*

*  ios_choose_var - select variable to branch on

*

*  SYNOPSIS

*

*  #include "glpios.h"

*  int ios_choose_var(glp_tree *T, int *next);

*

*  The routine ios_choose_var chooses a variable from the candidate

*  list to branch on. Additionally the routine provides a flag stored

*  in the location next to suggests which of the child subproblems

*  should be solved next.

*

*  RETURNS

*

*  The routine ios_choose_var returns the ordinal number of the column

*  choosen. */

static int branch_first(glp_tree *T, int *next);

static int branch_last(glp_tree *T, int *next);

static int branch_mostf(glp_tree *T, int *next);

static int branch_drtom(glp_tree *T, int *next);

int ios_choose_var(glp_tree *T, int *next)

{     int j;

      if (T->parm->br_tech == GLP_BR_FFV)

      {  /* branch on first fractional variable */

         j = branch_first(T, next);

      }

      else if (T->parm->br_tech == GLP_BR_LFV)

      {  /* branch on last fractional variable */

         j = branch_last(T, next);

      }

      else if (T->parm->br_tech == GLP_BR_MFV)

      {  /* branch on most fractional variable */

         j = branch_mostf(T, next);

      }

      else if (T->parm->br_tech == GLP_BR_DTH)

      {  /* branch using the heuristic by Dreebeck and Tomlin */

         j = branch_drtom(T, next);

      }

      else if (T->parm->br_tech == GLP_BR_PCH)

      {  /* hybrid pseudocost heuristic */

         j = ios_pcost_branch(T, next);

      }

      else

         xassert(T != T);

      return j;

}

/***********************************************************************

*  branch_first - choose first branching variable

*

*  This routine looks up the list of structural variables and chooses

*  the first one, which is of integer kind and has fractional value in

*  optimal solution to the current LP relaxation.

*

*  This routine also selects the branch to be solved next where integer

*  infeasibility of the chosen variable is less than in other one. */

static int branch_first(glp_tree *T, int *_next)

{     int j, next;

      double beta;

      /* choose the column to branch on */

      for (j = 1; j <= T->n; j++)

         if (T->non_int[j]) break;

      xassert(1 <= j && j <= T->n);

      /* select the branch to be solved next */

      beta = glp_get_col_prim(T->mip, j);

      if (beta - floor(beta) < ceil(beta) - beta)

         next = GLP_DN_BRNCH;

      else

         next = GLP_UP_BRNCH;

      *_next = next;

      return j;

}

/***********************************************************************

*  branch_last - choose last branching variable

*

*  This routine looks up the list of structural variables and chooses

*  the last one, which is of integer kind and has fractional value in

*  optimal solution to the current LP relaxation.

*

*  This routine also selects the branch to be solved next where integer

*  infeasibility of the chosen variable is less than in other one. */

static int branch_last(glp_tree *T, int *_next)

{     int j, next;

      double beta;

      /* choose the column to branch on */

      for (j = T->n; j >= 1; j--)

         if (T->non_int[j]) break;

      xassert(1 <= j && j <= T->n);

      /* select the branch to be solved next */

      beta = glp_get_col_prim(T->mip, j);

      if (beta - floor(beta) < ceil(beta) - beta)

         next = GLP_DN_BRNCH;

      else

         next = GLP_UP_BRNCH;

      *_next = next;

      return j;

}

/***********************************************************************

*  branch_mostf - choose most fractional branching variable

*

*  This routine looks up the list of structural variables and chooses

*  that one, which is of integer kind and has most fractional value in

*  optimal solution to the current LP relaxation.

*

*  This routine also selects the branch to be solved next where integer

*  infeasibility of the chosen variable is less than in other one.

*

*  (Alexander Martin notices that "...most infeasible is as good as

*  random...".) */

static int branch_mostf(glp_tree *T, int *_next)

{     int j, jj, next;

      double beta, most, temp;

      /* choose the column to branch on */

      jj = 0, most = DBL_MAX;

      for (j = 1; j <= T->n; j++)

      {  if (T->non_int[j])

         {  beta = glp_get_col_prim(T->mip, j);

            temp = floor(beta) + 0.5;

            if (most > fabs(beta - temp))

            {  jj = j, most = fabs(beta - temp);

               if (beta < temp)

                  next = GLP_DN_BRNCH;

               else

                  next = GLP_UP_BRNCH;

            }

         }

      }

      *_next = next;

      return jj;

}

/***********************************************************************

*  branch_drtom - choose branching var using Driebeck-Tomlin heuristic

*

*  This routine chooses a structural variable, which is required to be

*  integral and has fractional value in optimal solution of the current

*  LP relaxation, using a heuristic proposed by Driebeck and Tomlin.

*

*  The routine also selects the branch to be solved next, again due to

*  Driebeck and Tomlin.

*

*  This routine is based on the heuristic proposed in:

*

*  Driebeck N.J. An algorithm for the solution of mixed-integer

*  programming problems, Management Science, 12: 576-87 (1966);

*

*  and improved in:

*

*  Tomlin J.A. Branch and bound methods for integer and non-convex

*  programming, in J.Abadie (ed.), Integer and Nonlinear Programming,

*  North-Holland, Amsterdam, pp. 437-50 (1970).

*

*  Must note that this heuristic is time-expensive, because computing

*  one-step degradation (see the routine below) requires one BTRAN for

*  each fractional-valued structural variable. */

static int branch_drtom(glp_tree *T, int *_next)

{     glp_prob *mip = T->mip;

      int m = mip->m;

      int n = mip->n;

      char *non_int = T->non_int;

      int j, jj, k, t, next, kase, len, stat, *ind;

      double x, dk, alfa, delta_j, delta_k, delta_z, dz_dn, dz_up,

         dd_dn, dd_up, degrad, *val;

      /* basic solution of LP relaxation must be optimal */

      xassert(glp_get_status(mip) == GLP_OPT);

      /* allocate working arrays */

      ind = xcalloc(1+n, sizeof(int));

      val = xcalloc(1+n, sizeof(double));

      /* nothing has been chosen so far */

      jj = 0, degrad = -1.0;

      /* walk through the list of columns (structural variables) */

      for (j = 1; j <= n; j++)

      {  /* if j-th column is not marked as fractional, skip it */

         if (!non_int[j]) continue;

         /* obtain (fractional) value of j-th column in basic solution

            of LP relaxation */

         x = glp_get_col_prim(mip, j);

         /* since the value of j-th column is fractional, the column is

            basic; compute corresponding row of the simplex table */

         len = glp_eval_tab_row(mip, m+j, ind, val);

         /* the following fragment computes a change in the objective

            function: delta Z = new Z - old Z, where old Z is the

            objective value in the current optimal basis, and new Z is

            the objective value in the adjacent basis, for two cases:

            1) if new upper bound ub' = floor(x[j]) is introduced for

               j-th column (down branch);

            2) if new lower bound lb' = ceil(x[j]) is introduced for

               j-th column (up branch);

            since in both cases the solution remaining dual feasible

            becomes primal infeasible, one implicit simplex iteration

            is performed to determine the change delta Z;

            it is obvious that new Z, which is never better than old Z,

            is a lower (minimization) or upper (maximization) bound of

            the objective function for down- and up-branches. */

         for (kase = -1; kase <= +1; kase += 2)

         {  /* if kase < 0, the new upper bound of x[j] is introduced;

               in this case x[j] should decrease in order to leave the

               basis and go to its new upper bound */

            /* if kase > 0, the new lower bound of x[j] is introduced;

               in this case x[j] should increase in order to leave the

               basis and go to its new lower bound */

            /* apply the dual ratio test in order to determine which

               auxiliary or structural variable should enter the basis

               to keep dual feasibility */

            k = glp_dual_rtest(mip, len, ind, val, kase, 1e-9);

            if (k != 0) k = ind[k];

            /* if no non-basic variable has been chosen, LP relaxation

               of corresponding branch being primal infeasible and dual

               unbounded has no primal feasible solution; in this case

               the change delta Z is formally set to infinity */

            if (k == 0)

            {  delta_z =

                  (T->mip->dir == GLP_MIN ? +DBL_MAX : -DBL_MAX);

               goto skip;

            }

            /* row of the simplex table that corresponds to non-basic

               variable x[k] choosen by the dual ratio test is:

                  x[j] = ... + alfa * x[k] + ...

               where alfa is the influence coefficient (an element of

               the simplex table row) */

            /* determine the coefficient alfa */

            for (t = 1; t <= len; t++) if (ind[t] == k) break;

            xassert(1 <= t && t <= len);

            alfa = val[t];

            /* since in the adjacent basis the variable x[j] becomes

               non-basic, knowing its value in the current basis we can

               determine its change delta x[j] = new x[j] - old x[j] */

            delta_j = (kase < 0 ? floor(x) : ceil(x)) - x;

            /* and knowing the coefficient alfa we can determine the

               corresponding change delta x[k] = new x[k] - old x[k],

               where old x[k] is a value of x[k] in the current basis,

               and new x[k] is a value of x[k] in the adjacent basis */

            delta_k = delta_j / alfa;

            /* Tomlin noticed that if the variable x[k] is of integer

               kind, its change cannot be less (eventually) than one in

               the magnitude */

            if (k > m && glp_get_col_kind(mip, k-m) != GLP_CV)

            {  /* x[k] is structural integer variable */

               if (fabs(delta_k - floor(delta_k + 0.5)) > 1e-3)

               {  if (delta_k > 0.0)

                     delta_k = ceil(delta_k);  /* +3.14 -> +4 */

                  else

                     delta_k = floor(delta_k); /* -3.14 -> -4 */

               }

            }

            /* now determine the status and reduced cost of x[k] in the

               current basis */

            if (k <= m)

            {  stat = glp_get_row_stat(mip, k);

               dk = glp_get_row_dual(mip, k);

            }

            else

            {  stat = glp_get_col_stat(mip, k-m);

               dk = glp_get_col_dual(mip, k-m);

            }

            /* if the current basis is dual degenerate, some reduced

               costs which are close to zero may have wrong sign due to

               round-off errors, so correct the sign of d[k] */

            switch (T->mip->dir)

            {  case GLP_MIN:

                  if (stat == GLP_NL && dk < 0.0 ||

                      stat == GLP_NU && dk > 0.0 ||

                      stat == GLP_NF) dk = 0.0;

                  break;

               case GLP_MAX:

                  if (stat == GLP_NL && dk > 0.0 ||

                      stat == GLP_NU && dk < 0.0 ||

                      stat == GLP_NF) dk = 0.0;

                  break;

               default:

                  xassert(T != T);

            }

            /* now knowing the change of x[k] and its reduced cost d[k]

               we can compute the corresponding change in the objective

               function delta Z = new Z - old Z = d[k] * delta x[k];

               note that due to Tomlin's modification new Z can be even

               worse than in the adjacent basis */

            delta_z = dk * delta_k;

skip:       /* new Z is never better than old Z, therefore the change

               delta Z is always non-negative (in case of minimization)

               or non-positive (in case of maximization) */

            switch (T->mip->dir)

            {  case GLP_MIN: xassert(delta_z >= 0.0); break;

               case GLP_MAX: xassert(delta_z <= 0.0); break;

               default: xassert(T != T);

            }

            /* save the change in the objective fnction for down- and

               up-branches, respectively */

            if (kase < 0) dz_dn = delta_z; else dz_up = delta_z;

         }

         /* thus, in down-branch no integer feasible solution can be

            better than Z + dz_dn, and in up-branch no integer feasible

            solution can be better than Z + dz_up, where Z is value of

            the objective function in the current basis */

         /* following the heuristic by Driebeck and Tomlin we choose a

            column (i.e. structural variable) which provides largest

            degradation of the objective function in some of branches;

            besides, we select the branch with smaller degradation to

            be solved next and keep other branch with larger degradation

            in the active list hoping to minimize the number of further

            backtrackings */

         if (degrad < fabs(dz_dn) || degrad < fabs(dz_up))

         {  jj = j;

            if (fabs(dz_dn) < fabs(dz_up))

            {  /* select down branch to be solved next */

               next = GLP_DN_BRNCH;

               degrad = fabs(dz_up);

            }

            else

            {  /* select up branch to be solved next */

               next = GLP_UP_BRNCH;

               degrad = fabs(dz_dn);

            }

            /* save the objective changes for printing */

            dd_dn = dz_dn, dd_up = dz_up;

            /* if down- or up-branch has no feasible solution, we does

               not need to consider other candidates (in principle, the

               corresponding branch could be pruned right now) */

            if (degrad == DBL_MAX) break;

         }

      }

      /* free working arrays */

      xfree(ind);

      xfree(val);

      /* something must be chosen */

      xassert(1 <= jj && jj <= n);

#if 1 /* 02/XI-2009 */

      if (degrad < 1e-6 * (1.0 + 0.001 * fabs(mip->obj_val)))

      {  jj = branch_mostf(T, &next);

         goto done;

      }

#endif

      if (T->parm->msg_lev >= GLP_MSG_DBG)

      {  xprintf("branch_drtom: column %d chosen to branch on\n", jj);

         if (fabs(dd_dn) == DBL_MAX)

            xprintf("branch_drtom: down-branch is infeasible\n");

         else

            xprintf("branch_drtom: down-branch bound is %.9e\n",

               lpx_get_obj_val(mip) + dd_dn);

         if (fabs(dd_up) == DBL_MAX)

            xprintf("branch_drtom: up-branch   is infeasible\n");

         else

            xprintf("branch_drtom: up-branch   bound is %.9e\n",

               lpx_get_obj_val(mip) + dd_up);

      }

done: *_next = next;

      return jj;

}

/**********************************************************************/

struct csa

{     /* common storage area */

      int *dn_cnt; /* int dn_cnt[1+n]; */

      /* dn_cnt[j] is the number of subproblems, whose LP relaxations

         have been solved and which are down-branches for variable x[j];

         dn_cnt[j] = 0 means the down pseudocost is uninitialized */

      double *dn_sum; /* double dn_sum[1+n]; */

      /* dn_sum[j] is the sum of per unit degradations of the objective

         over all dn_cnt[j] subproblems */

      int *up_cnt; /* int up_cnt[1+n]; */

      /* up_cnt[j] is the number of subproblems, whose LP relaxations

         have been solved and which are up-branches for variable x[j];

         up_cnt[j] = 0 means the up pseudocost is uninitialized */

      double *up_sum; /* double up_sum[1+n]; */

      /* up_sum[j] is the sum of per unit degradations of the objective

         over all up_cnt[j] subproblems */

};

void *ios_pcost_init(glp_tree *tree)

{     /* initialize working data used on pseudocost branching */

      struct csa *csa;

      int n = tree->n, j;

      csa = xmalloc(sizeof(struct csa));

      csa->dn_cnt = xcalloc(1+n, sizeof(int));

      csa->dn_sum = xcalloc(1+n, sizeof(double));

      csa->up_cnt = xcalloc(1+n, sizeof(int));

      csa->up_sum = xcalloc(1+n, sizeof(double));

      for (j = 1; j <= n; j++)

      {  csa->dn_cnt[j] = csa->up_cnt[j] = 0;

         csa->dn_sum[j] = csa->up_sum[j] = 0.0;

      }

      return csa;

}

static double eval_degrad(glp_prob *P, int j, double bnd)

{     /* compute degradation of the objective on fixing x[j] at given

         value with a limited number of dual simplex iterations */

      /* this routine fixes column x[j] at specified value bnd,

         solves resulting LP, and returns a lower bound to degradation

         of the objective, degrad >= 0 */

      glp_prob *lp;

      glp_smcp parm;

      int ret;

      double degrad;

      /* the current basis must be optimal */

      xassert(glp_get_status(P) == GLP_OPT);

      /* create a copy of P */

      lp = glp_create_prob();

      glp_copy_prob(lp, P, 0);

      /* fix column x[j] at specified value */

      glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd);

      /* try to solve resulting LP */

      glp_init_smcp(&parm);

      parm.msg_lev = GLP_MSG_OFF;

      parm.meth = GLP_DUAL;

      parm.it_lim = 30;

      parm.out_dly = 1000;

      parm.meth = GLP_DUAL;

      ret = glp_simplex(lp, &parm);

      if (ret == 0 || ret == GLP_EITLIM)

      {  if (glp_get_prim_stat(lp) == GLP_NOFEAS)

         {  /* resulting LP has no primal feasible solution */

            degrad = DBL_MAX;

         }

         else if (glp_get_dual_stat(lp) == GLP_FEAS)

         {  /* resulting basis is optimal or at least dual feasible,

               so we have the correct lower bound to degradation */

            if (P->dir == GLP_MIN)

               degrad = lp->obj_val - P->obj_val;

            else if (P->dir == GLP_MAX)

               degrad = P->obj_val - lp->obj_val;

            else

               xassert(P != P);

            /* degradation cannot be negative by definition */

            /* note that the lower bound to degradation may be close

               to zero even if its exact value is zero due to round-off

               errors on computing the objective value */

            if (degrad < 1e-6 * (1.0 + 0.001 * fabs(P->obj_val)))

               degrad = 0.0;

         }

         else

         {  /* the final basis reported by the simplex solver is dual

               infeasible, so we cannot determine a non-trivial lower

               bound to degradation */

            degrad = 0.0;

         }

      }

      else

      {  /* the simplex solver failed */

         degrad = 0.0;

      }

      /* delete the copy of P */

      glp_delete_prob(lp);

      return degrad;

}

void ios_pcost_update(glp_tree *tree)

{     /* update history information for pseudocost branching */

      /* this routine is called every time when LP relaxation of the

         current subproblem has been solved to optimality with all lazy

         and cutting plane constraints included */

      int j;

      double dx, dz, psi;

      struct csa *csa = tree->pcost;

      xassert(csa != NULL);

      xassert(tree->curr != NULL);

      /* if the current subproblem is the root, skip updating */

      if (tree->curr->up == NULL) goto skip;

      /* determine branching variable x[j], which was used in the

         parent subproblem to create the current subproblem */

      j = tree->curr->up->br_var;

      xassert(1 <= j && j <= tree->n);

      /* determine the change dx[j] = new x[j] - old x[j],

         where new x[j] is a value of x[j] in optimal solution to LP

         relaxation of the current subproblem, old x[j] is a value of

         x[j] in optimal solution to LP relaxation of the parent

         subproblem */

      dx = tree->mip->col[j]->prim - tree->curr->up->br_val;

      xassert(dx != 0.0);

      /* determine corresponding change dz = new dz - old dz in the

         objective function value */

      dz = tree->mip->obj_val - tree->curr->up->lp_obj;

      /* determine per unit degradation of the objective function */

      psi = fabs(dz / dx);

      /* update history information */

      if (dx < 0.0)

      {  /* the current subproblem is down-branch */

         csa->dn_cnt[j]++;

         csa->dn_sum[j] += psi;

      }

      else /* dx > 0.0 */

      {  /* the current subproblem is up-branch */

         csa->up_cnt[j]++;

         csa->up_sum[j] += psi;

      }

skip: return;

}

void ios_pcost_free(glp_tree *tree)

{     /* free working area used on pseudocost branching */

      struct csa *csa = tree->pcost;

      xassert(csa != NULL);

      xfree(csa->dn_cnt);

      xfree(csa->dn_sum);

      xfree(csa->up_cnt);

      xfree(csa->up_sum);

      xfree(csa);

      tree->pcost = NULL;

      return;

}

static double eval_psi(glp_tree *T, int j, int brnch)

{     /* compute estimation of pseudocost of variable x[j] for down-

         or up-branch */

      struct csa *csa = T->pcost;

      double beta, degrad, psi;

      xassert(csa != NULL);

      xassert(1 <= j && j <= T->n);

      if (brnch == GLP_DN_BRNCH)

      {  /* down-branch */

         if (csa->dn_cnt[j] == 0)

         {  /* initialize down pseudocost */

            beta = T->mip->col[j]->prim;

            degrad = eval_degrad(T->mip, j, floor(beta));

            if (degrad == DBL_MAX)

            {  psi = DBL_MAX;

               goto done;

            }

            csa->dn_cnt[j] = 1;

            csa->dn_sum[j] = degrad / (beta - floor(beta));

         }

         psi = csa->dn_sum[j] / (double)csa->dn_cnt[j];

      }

      else if (brnch == GLP_UP_BRNCH)

      {  /* up-branch */

         if (csa->up_cnt[j] == 0)

         {  /* initialize up pseudocost */

            beta = T->mip->col[j]->prim;

            degrad = eval_degrad(T->mip, j, ceil(beta));

            if (degrad == DBL_MAX)

            {  psi = DBL_MAX;

               goto done;

            }

            csa->up_cnt[j] = 1;

            csa->up_sum[j] = degrad / (ceil(beta) - beta);

         }

         psi = csa->up_sum[j] / (double)csa->up_cnt[j];

      }

      else

         xassert(brnch != brnch);

done: return psi;

}

static void progress(glp_tree *T)

{     /* display progress of pseudocost initialization */

      struct csa *csa = T->pcost;

      int j, nv = 0, ni = 0;

      for (j = 1; j <= T->n; j++)

      {  if (glp_ios_can_branch(T, j))

         {  nv++;

            if (csa->dn_cnt[j] > 0 && csa->up_cnt[j] > 0) ni++;

         }

      }

      xprintf("Pseudocosts initialized for %d of %d variables\n",

         ni, nv);

      return;

}

int ios_pcost_branch(glp_tree *T, int *_next)

{     /* choose branching variable with pseudocost branching */

      glp_long t = xtime();

      int j, jjj, sel;

      double beta, psi, d1, d2, d, dmax;

      /* initialize the working arrays */

      if (T->pcost == NULL)

         T->pcost = ios_pcost_init(T);

      /* nothing has been chosen so far */

      jjj = 0, dmax = -1.0;

      /* go through the list of branching candidates */

      for (j = 1; j <= T->n; j++)

      {  if (!glp_ios_can_branch(T, j)) continue;

         /* determine primal value of x[j] in optimal solution to LP

            relaxation of the current subproblem */

         beta = T->mip->col[j]->prim;

         /* estimate pseudocost of x[j] for down-branch */

         psi = eval_psi(T, j, GLP_DN_BRNCH);

         if (psi == DBL_MAX)

         {  /* down-branch has no primal feasible solution */

            jjj = j, sel = GLP_DN_BRNCH;

            goto done;

         }

         /* estimate degradation of the objective for down-branch */

         d1 = psi * (beta - floor(beta));

         /* estimate pseudocost of x[j] for up-branch */

         psi = eval_psi(T, j, GLP_UP_BRNCH);

         if (psi == DBL_MAX)

         {  /* up-branch has no primal feasible solution */

            jjj = j, sel = GLP_UP_BRNCH;

            goto done;

         }

         /* estimate degradation of the objective for up-branch */

         d2 = psi * (ceil(beta) - beta);

         /* determine d = max(d1, d2) */

         d = (d1 > d2 ? d1 : d2);

         /* choose x[j] which provides maximal estimated degradation of

            the objective either in down- or up-branch */

         if (dmax < d)

         {  dmax = d;

            jjj = j;

            /* continue the search from a subproblem, where degradation

               is less than in other one */

            sel = (d1 <= d2 ? GLP_DN_BRNCH : GLP_UP_BRNCH);

         }

         /* display progress of pseudocost initialization */

         if (T->parm->msg_lev >= GLP_ON)

         {  if (xdifftime(xtime(), t) >= 10.0)

            {  progress(T);

               t = xtime();

            }

         }

      }

      if (dmax == 0.0)

      {  /* no degradation is indicated; choose a variable having most

            fractional value */

         jjj = branch_mostf(T, &sel);

      }

done: *_next = sel;

      return jjj;

}

/* eof */