1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpios11.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,280 @@
     1.4 +/* glpios11.c (process cuts stored in the local cut pool) */
     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_process_cuts - process cuts stored in the local cut pool
    1.34 +*
    1.35 +*  SYNOPSIS
    1.36 +*
    1.37 +*  #include "glpios.h"
    1.38 +*  void ios_process_cuts(glp_tree *T);
    1.39 +*
    1.40 +*  DESCRIPTION
    1.41 +*
    1.42 +*  The routine ios_process_cuts analyzes each cut currently stored in
    1.43 +*  the local cut pool, which must be non-empty, and either adds the cut
    1.44 +*  to the current subproblem or just discards it. All cuts are assumed
    1.45 +*  to be locally valid. On exit the local cut pool remains unchanged.
    1.46 +*
    1.47 +*  REFERENCES
    1.48 +*
    1.49 +*  1. E.Balas, S.Ceria, G.Cornuejols, "Mixed 0-1 Programming by
    1.50 +*     Lift-and-Project in a Branch-and-Cut Framework", Management Sc.,
    1.51 +*     42 (1996) 1229-1246.
    1.52 +*
    1.53 +*  2. G.Andreello, A.Caprara, and M.Fischetti, "Embedding Cuts in
    1.54 +*     a Branch&Cut Framework: a Computational Study with {0,1/2}-Cuts",
    1.55 +*     Preliminary Draft, October 28, 2003, pp.6-8. */
    1.56 +
    1.57 +struct info
    1.58 +{     /* estimated cut efficiency */
    1.59 +      IOSCUT *cut;
    1.60 +      /* pointer to cut in the cut pool */
    1.61 +      char flag;
    1.62 +      /* if this flag is set, the cut is included into the current
    1.63 +         subproblem */
    1.64 +      double eff;
    1.65 +      /* cut efficacy (normalized residual) */
    1.66 +      double deg;
    1.67 +      /* lower bound to objective degradation */
    1.68 +};
    1.69 +
    1.70 +static int fcmp(const void *arg1, const void *arg2)
    1.71 +{     const struct info *info1 = arg1, *info2 = arg2;
    1.72 +      if (info1->deg == 0.0 && info2->deg == 0.0)
    1.73 +      {  if (info1->eff > info2->eff) return -1;
    1.74 +         if (info1->eff < info2->eff) return +1;
    1.75 +      }
    1.76 +      else
    1.77 +      {  if (info1->deg > info2->deg) return -1;
    1.78 +         if (info1->deg < info2->deg) return +1;
    1.79 +      }
    1.80 +      return 0;
    1.81 +}
    1.82 +
    1.83 +static double parallel(IOSCUT *a, IOSCUT *b, double work[]);
    1.84 +
    1.85 +void ios_process_cuts(glp_tree *T)
    1.86 +{     IOSPOOL *pool;
    1.87 +      IOSCUT *cut;
    1.88 +      IOSAIJ *aij;
    1.89 +      struct info *info;
    1.90 +      int k, kk, max_cuts, len, ret, *ind;
    1.91 +      double *val, *work;
    1.92 +      /* the current subproblem must exist */
    1.93 +      xassert(T->curr != NULL);
    1.94 +      /* the pool must exist and be non-empty */
    1.95 +      pool = T->local;
    1.96 +      xassert(pool != NULL);
    1.97 +      xassert(pool->size > 0);
    1.98 +      /* allocate working arrays */
    1.99 +      info = xcalloc(1+pool->size, sizeof(struct info));
   1.100 +      ind = xcalloc(1+T->n, sizeof(int));
   1.101 +      val = xcalloc(1+T->n, sizeof(double));
   1.102 +      work = xcalloc(1+T->n, sizeof(double));
   1.103 +      for (k = 1; k <= T->n; k++) work[k] = 0.0;
   1.104 +      /* build the list of cuts stored in the cut pool */
   1.105 +      for (k = 0, cut = pool->head; cut != NULL; cut = cut->next)
   1.106 +         k++, info[k].cut = cut, info[k].flag = 0;
   1.107 +      xassert(k == pool->size);
   1.108 +      /* estimate efficiency of all cuts in the cut pool */
   1.109 +      for (k = 1; k <= pool->size; k++)
   1.110 +      {  double temp, dy, dz;
   1.111 +         cut = info[k].cut;
   1.112 +         /* build the vector of cut coefficients and compute its
   1.113 +            Euclidean norm */
   1.114 +         len = 0; temp = 0.0;
   1.115 +         for (aij = cut->ptr; aij != NULL; aij = aij->next)
   1.116 +         {  xassert(1 <= aij->j && aij->j <= T->n);
   1.117 +            len++, ind[len] = aij->j, val[len] = aij->val;
   1.118 +            temp += aij->val * aij->val;
   1.119 +         }
   1.120 +         if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON;
   1.121 +         /* transform the cut to express it only through non-basic
   1.122 +            (auxiliary and structural) variables */
   1.123 +         len = glp_transform_row(T->mip, len, ind, val);
   1.124 +         /* determine change in the cut value and in the objective
   1.125 +            value for the adjacent basis by simulating one step of the
   1.126 +            dual simplex */
   1.127 +         ret = _glp_analyze_row(T->mip, len, ind, val, cut->type,
   1.128 +            cut->rhs, 1e-9, NULL, NULL, NULL, NULL, &dy, &dz);
   1.129 +         /* determine normalized residual and lower bound to objective
   1.130 +            degradation */
   1.131 +         if (ret == 0)
   1.132 +         {  info[k].eff = fabs(dy) / sqrt(temp);
   1.133 +            /* if some reduced costs violates (slightly) their zero
   1.134 +               bounds (i.e. have wrong signs) due to round-off errors,
   1.135 +               dz also may have wrong sign being close to zero */
   1.136 +            if (T->mip->dir == GLP_MIN)
   1.137 +            {  if (dz < 0.0) dz = 0.0;
   1.138 +               info[k].deg = + dz;
   1.139 +            }
   1.140 +            else /* GLP_MAX */
   1.141 +            {  if (dz > 0.0) dz = 0.0;
   1.142 +               info[k].deg = - dz;
   1.143 +            }
   1.144 +         }
   1.145 +         else if (ret == 1)
   1.146 +         {  /* the constraint is not violated at the current point */
   1.147 +            info[k].eff = info[k].deg = 0.0;
   1.148 +         }
   1.149 +         else if (ret == 2)
   1.150 +         {  /* no dual feasible adjacent basis exists */
   1.151 +            info[k].eff = 1.0;
   1.152 +            info[k].deg = DBL_MAX;
   1.153 +         }
   1.154 +         else
   1.155 +            xassert(ret != ret);
   1.156 +         /* if the degradation is too small, just ignore it */
   1.157 +         if (info[k].deg < 0.01) info[k].deg = 0.0;
   1.158 +      }
   1.159 +      /* sort the list of cuts by decreasing objective degradation and
   1.160 +         then by decreasing efficacy */
   1.161 +      qsort(&info[1], pool->size, sizeof(struct info), fcmp);
   1.162 +      /* only first (most efficient) max_cuts in the list are qualified
   1.163 +         as candidates to be added to the current subproblem */
   1.164 +      max_cuts = (T->curr->level == 0 ? 90 : 10);
   1.165 +      if (max_cuts > pool->size) max_cuts = pool->size;
   1.166 +      /* add cuts to the current subproblem */
   1.167 +#if 0
   1.168 +      xprintf("*** adding cuts ***\n");
   1.169 +#endif
   1.170 +      for (k = 1; k <= max_cuts; k++)
   1.171 +      {  int i, len;
   1.172 +         /* if this cut seems to be inefficient, skip it */
   1.173 +         if (info[k].deg < 0.01 && info[k].eff < 0.01) continue;
   1.174 +         /* if the angle between this cut and every other cut included
   1.175 +            in the current subproblem is small, skip this cut */
   1.176 +         for (kk = 1; kk < k; kk++)
   1.177 +         {  if (info[kk].flag)
   1.178 +            {  if (parallel(info[k].cut, info[kk].cut, work) > 0.90)
   1.179 +                  break;
   1.180 +            }
   1.181 +         }
   1.182 +         if (kk < k) continue;
   1.183 +         /* add this cut to the current subproblem */
   1.184 +#if 0
   1.185 +         xprintf("eff = %g; deg = %g\n", info[k].eff, info[k].deg);
   1.186 +#endif
   1.187 +         cut = info[k].cut, info[k].flag = 1;
   1.188 +         i = glp_add_rows(T->mip, 1);
   1.189 +         if (cut->name != NULL)
   1.190 +            glp_set_row_name(T->mip, i, cut->name);
   1.191 +         xassert(T->mip->row[i]->origin == GLP_RF_CUT);
   1.192 +         T->mip->row[i]->klass = cut->klass;
   1.193 +         len = 0;
   1.194 +         for (aij = cut->ptr; aij != NULL; aij = aij->next)
   1.195 +            len++, ind[len] = aij->j, val[len] = aij->val;
   1.196 +         glp_set_mat_row(T->mip, i, len, ind, val);
   1.197 +         xassert(cut->type == GLP_LO || cut->type == GLP_UP);
   1.198 +         glp_set_row_bnds(T->mip, i, cut->type, cut->rhs, cut->rhs);
   1.199 +      }
   1.200 +      /* free working arrays */
   1.201 +      xfree(info);
   1.202 +      xfree(ind);
   1.203 +      xfree(val);
   1.204 +      xfree(work);
   1.205 +      return;
   1.206 +}
   1.207 +
   1.208 +#if 0
   1.209 +/***********************************************************************
   1.210 +*  Given a cut a * x >= b (<= b) the routine efficacy computes the cut
   1.211 +*  efficacy as follows:
   1.212 +*
   1.213 +*     eff = d * (a * x~ - b) / ||a||,
   1.214 +*
   1.215 +*  where d is -1 (in case of '>= b') or +1 (in case of '<= b'), x~ is
   1.216 +*  the vector of values of structural variables in optimal solution to
   1.217 +*  LP relaxation of the current subproblem, ||a|| is the Euclidean norm
   1.218 +*  of the vector of cut coefficients.
   1.219 +*
   1.220 +*  If the cut is violated at point x~, the efficacy eff is positive,
   1.221 +*  and its value is the Euclidean distance between x~ and the cut plane
   1.222 +*  a * x = b in the space of structural variables.
   1.223 +*
   1.224 +*  Following geometrical intuition, it is quite natural to consider
   1.225 +*  this distance as a first-order measure of the expected efficacy of
   1.226 +*  the cut: the larger the distance the better the cut [1]. */
   1.227 +
   1.228 +static double efficacy(glp_tree *T, IOSCUT *cut)
   1.229 +{     glp_prob *mip = T->mip;
   1.230 +      IOSAIJ *aij;
   1.231 +      double s = 0.0, t = 0.0, temp;
   1.232 +      for (aij = cut->ptr; aij != NULL; aij = aij->next)
   1.233 +      {  xassert(1 <= aij->j && aij->j <= mip->n);
   1.234 +         s += aij->val * mip->col[aij->j]->prim;
   1.235 +         t += aij->val * aij->val;
   1.236 +      }
   1.237 +      temp = sqrt(t);
   1.238 +      if (temp < DBL_EPSILON) temp = DBL_EPSILON;
   1.239 +      if (cut->type == GLP_LO)
   1.240 +         temp = (s >= cut->rhs ? 0.0 : (cut->rhs - s) / temp);
   1.241 +      else if (cut->type == GLP_UP)
   1.242 +         temp = (s <= cut->rhs ? 0.0 : (s - cut->rhs) / temp);
   1.243 +      else
   1.244 +         xassert(cut != cut);
   1.245 +      return temp;
   1.246 +}
   1.247 +#endif
   1.248 +
   1.249 +/***********************************************************************
   1.250 +*  Given two cuts a1 * x >= b1 (<= b1) and a2 * x >= b2 (<= b2) the
   1.251 +*  routine parallel computes the cosine of angle between the cut planes
   1.252 +*  a1 * x = b1 and a2 * x = b2 (which is the acute angle between two
   1.253 +*  normals to these planes) in the space of structural variables as
   1.254 +*  follows:
   1.255 +*
   1.256 +*     cos phi = (a1' * a2) / (||a1|| * ||a2||),
   1.257 +*
   1.258 +*  where (a1' * a2) is a dot product of vectors of cut coefficients,
   1.259 +*  ||a1|| and ||a2|| are Euclidean norms of vectors a1 and a2.
   1.260 +*
   1.261 +*  Note that requirement cos phi = 0 forces the cuts to be orthogonal,
   1.262 +*  i.e. with disjoint support, while requirement cos phi <= 0.999 means
   1.263 +*  only avoiding duplicate (parallel) cuts [1]. */
   1.264 +
   1.265 +static double parallel(IOSCUT *a, IOSCUT *b, double work[])
   1.266 +{     IOSAIJ *aij;
   1.267 +      double s = 0.0, sa = 0.0, sb = 0.0, temp;
   1.268 +      for (aij = a->ptr; aij != NULL; aij = aij->next)
   1.269 +      {  work[aij->j] = aij->val;
   1.270 +         sa += aij->val * aij->val;
   1.271 +      }
   1.272 +      for (aij = b->ptr; aij != NULL; aij = aij->next)
   1.273 +      {  s += work[aij->j] * aij->val;
   1.274 +         sb += aij->val * aij->val;
   1.275 +      }
   1.276 +      for (aij = a->ptr; aij != NULL; aij = aij->next)
   1.277 +         work[aij->j] = 0.0;
   1.278 +      temp = sqrt(sa) * sqrt(sb);
   1.279 +      if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON;
   1.280 +      return s / temp;
   1.281 +}
   1.282 +
   1.283 +/* eof */