1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpios11.c Mon Dec 06 13:09:21 2010 +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 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 */