/* 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 */