/* glpios11.c (process cuts stored in the local cut pool) */

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

*  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_process_cuts - process cuts stored in the local cut pool

*

*  SYNOPSIS

*

*  #include "glpios.h"

*  void ios_process_cuts(glp_tree *T);

*

*  DESCRIPTION

*

*  The routine ios_process_cuts analyzes each cut currently stored in

*  the local cut pool, which must be non-empty, and either adds the cut

*  to the current subproblem or just discards it. All cuts are assumed

*  to be locally valid. On exit the local cut pool remains unchanged.

*

*  REFERENCES

*

*  1. E.Balas, S.Ceria, G.Cornuejols, "Mixed 0-1 Programming by

*     Lift-and-Project in a Branch-and-Cut Framework", Management Sc.,

*     42 (1996) 1229-1246.

*

*  2. G.Andreello, A.Caprara, and M.Fischetti, "Embedding Cuts in

*     a Branch&Cut Framework: a Computational Study with {0,1/2}-Cuts",

*     Preliminary Draft, October 28, 2003, pp.6-8. */

struct info

{     /* estimated cut efficiency */

      IOSCUT *cut;

      /* pointer to cut in the cut pool */

      char flag;

      /* if this flag is set, the cut is included into the current

         subproblem */

      double eff;

      /* cut efficacy (normalized residual) */

      double deg;

      /* lower bound to objective degradation */

};

static int fcmp(const void *arg1, const void *arg2)

{     const struct info *info1 = arg1, *info2 = arg2;

      if (info1->deg == 0.0 && info2->deg == 0.0)

      {  if (info1->eff > info2->eff) return -1;

         if (info1->eff < info2->eff) return +1;

      }

      else

      {  if (info1->deg > info2->deg) return -1;

         if (info1->deg < info2->deg) return +1;

      }

      return 0;

}

static double parallel(IOSCUT *a, IOSCUT *b, double work[]);

void ios_process_cuts(glp_tree *T)

{     IOSPOOL *pool;

      IOSCUT *cut;

      IOSAIJ *aij;

      struct info *info;

      int k, kk, max_cuts, len, ret, *ind;

      double *val, *work;

      /* the current subproblem must exist */

      xassert(T->curr != NULL);

      /* the pool must exist and be non-empty */

      pool = T->local;

      xassert(pool != NULL);

      xassert(pool->size > 0);

      /* allocate working arrays */

      info = xcalloc(1+pool->size, sizeof(struct info));

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

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

      work = xcalloc(1+T->n, sizeof(double));

      for (k = 1; k <= T->n; k++) work[k] = 0.0;

      /* build the list of cuts stored in the cut pool */

      for (k = 0, cut = pool->head; cut != NULL; cut = cut->next)

         k++, info[k].cut = cut, info[k].flag = 0;

      xassert(k == pool->size);

      /* estimate efficiency of all cuts in the cut pool */

      for (k = 1; k <= pool->size; k++)

      {  double temp, dy, dz;

         cut = info[k].cut;

         /* build the vector of cut coefficients and compute its

            Euclidean norm */

         len = 0; temp = 0.0;

         for (aij = cut->ptr; aij != NULL; aij = aij->next)

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

            len++, ind[len] = aij->j, val[len] = aij->val;

            temp += aij->val * aij->val;

         }

         if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON;

         /* transform the cut to express it only through non-basic

            (auxiliary and structural) variables */

         len = glp_transform_row(T->mip, len, ind, val);

         /* determine change in the cut value and in the objective

            value for the adjacent basis by simulating one step of the

            dual simplex */

         ret = _glp_analyze_row(T->mip, len, ind, val, cut->type,

            cut->rhs, 1e-9, NULL, NULL, NULL, NULL, &dy, &dz);

         /* determine normalized residual and lower bound to objective

            degradation */

         if (ret == 0)

         {  info[k].eff = fabs(dy) / sqrt(temp);

            /* if some reduced costs violates (slightly) their zero

               bounds (i.e. have wrong signs) due to round-off errors,

               dz also may have wrong sign being close to zero */

            if (T->mip->dir == GLP_MIN)

            {  if (dz < 0.0) dz = 0.0;

               info[k].deg = + dz;

            }

            else /* GLP_MAX */

            {  if (dz > 0.0) dz = 0.0;

               info[k].deg = - dz;

            }

         }

         else if (ret == 1)

         {  /* the constraint is not violated at the current point */

            info[k].eff = info[k].deg = 0.0;

         }

         else if (ret == 2)

         {  /* no dual feasible adjacent basis exists */

            info[k].eff = 1.0;

            info[k].deg = DBL_MAX;

         }

         else

            xassert(ret != ret);

         /* if the degradation is too small, just ignore it */

         if (info[k].deg < 0.01) info[k].deg = 0.0;

      }

      /* sort the list of cuts by decreasing objective degradation and

         then by decreasing efficacy */

      qsort(&info[1], pool->size, sizeof(struct info), fcmp);

      /* only first (most efficient) max_cuts in the list are qualified

         as candidates to be added to the current subproblem */

      max_cuts = (T->curr->level == 0 ? 90 : 10);

      if (max_cuts > pool->size) max_cuts = pool->size;

      /* add cuts to the current subproblem */

#if 0

      xprintf("*** adding cuts ***\n");

#endif

      for (k = 1; k <= max_cuts; k++)

      {  int i, len;

         /* if this cut seems to be inefficient, skip it */

         if (info[k].deg < 0.01 && info[k].eff < 0.01) continue;

         /* if the angle between this cut and every other cut included

            in the current subproblem is small, skip this cut */

         for (kk = 1; kk < k; kk++)

         {  if (info[kk].flag)

            {  if (parallel(info[k].cut, info[kk].cut, work) > 0.90)

                  break;

            }

         }

         if (kk < k) continue;

         /* add this cut to the current subproblem */

#if 0

         xprintf("eff = %g; deg = %g\n", info[k].eff, info[k].deg);

#endif

         cut = info[k].cut, info[k].flag = 1;

         i = glp_add_rows(T->mip, 1);

         if (cut->name != NULL)

            glp_set_row_name(T->mip, i, cut->name);

         xassert(T->mip->row[i]->origin == GLP_RF_CUT);

         T->mip->row[i]->klass = cut->klass;

         len = 0;

         for (aij = cut->ptr; aij != NULL; aij = aij->next)

            len++, ind[len] = aij->j, val[len] = aij->val;

         glp_set_mat_row(T->mip, i, len, ind, val);

         xassert(cut->type == GLP_LO || cut->type == GLP_UP);

         glp_set_row_bnds(T->mip, i, cut->type, cut->rhs, cut->rhs);

      }

      /* free working arrays */

      xfree(info);

      xfree(ind);

      xfree(val);

      xfree(work);

      return;

}

#if 0

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

*  Given a cut a * x >= b (<= b) the routine efficacy computes the cut

*  efficacy as follows:

*

*     eff = d * (a * x~ - b) / ||a||,

*

*  where d is -1 (in case of '>= b') or +1 (in case of '<= b'), x~ is

*  the vector of values of structural variables in optimal solution to

*  LP relaxation of the current subproblem, ||a|| is the Euclidean norm

*  of the vector of cut coefficients.

*

*  If the cut is violated at point x~, the efficacy eff is positive,

*  and its value is the Euclidean distance between x~ and the cut plane

*  a * x = b in the space of structural variables.

*

*  Following geometrical intuition, it is quite natural to consider

*  this distance as a first-order measure of the expected efficacy of

*  the cut: the larger the distance the better the cut [1]. */

static double efficacy(glp_tree *T, IOSCUT *cut)

{     glp_prob *mip = T->mip;

      IOSAIJ *aij;

      double s = 0.0, t = 0.0, temp;

      for (aij = cut->ptr; aij != NULL; aij = aij->next)

      {  xassert(1 <= aij->j && aij->j <= mip->n);

         s += aij->val * mip->col[aij->j]->prim;

         t += aij->val * aij->val;

      }

      temp = sqrt(t);

      if (temp < DBL_EPSILON) temp = DBL_EPSILON;

      if (cut->type == GLP_LO)

         temp = (s >= cut->rhs ? 0.0 : (cut->rhs - s) / temp);

      else if (cut->type == GLP_UP)

         temp = (s <= cut->rhs ? 0.0 : (s - cut->rhs) / temp);

      else

         xassert(cut != cut);

      return temp;

}

#endif

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

*  Given two cuts a1 * x >= b1 (<= b1) and a2 * x >= b2 (<= b2) the

*  routine parallel computes the cosine of angle between the cut planes

*  a1 * x = b1 and a2 * x = b2 (which is the acute angle between two

*  normals to these planes) in the space of structural variables as

*  follows:

*

*     cos phi = (a1' * a2) / (||a1|| * ||a2||),

*

*  where (a1' * a2) is a dot product of vectors of cut coefficients,

*  ||a1|| and ||a2|| are Euclidean norms of vectors a1 and a2.

*

*  Note that requirement cos phi = 0 forces the cuts to be orthogonal,

*  i.e. with disjoint support, while requirement cos phi <= 0.999 means

*  only avoiding duplicate (parallel) cuts [1]. */

static double parallel(IOSCUT *a, IOSCUT *b, double work[])

{     IOSAIJ *aij;

      double s = 0.0, sa = 0.0, sb = 0.0, temp;

      for (aij = a->ptr; aij != NULL; aij = aij->next)

      {  work[aij->j] = aij->val;

         sa += aij->val * aij->val;

      }

      for (aij = b->ptr; aij != NULL; aij = aij->next)

      {  s += work[aij->j] * aij->val;

         sb += aij->val * aij->val;

      }

      for (aij = a->ptr; aij != NULL; aij = aij->next)

         work[aij->j] = 0.0;

      temp = sqrt(sa) * sqrt(sb);

      if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON;

      return s / temp;

}

/* eof */