lemon-project-template-glpk

annotate deps/glpk/src/glpios11.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
rev   line source
alpar@9 1 /* glpios11.c (process cuts stored in the local cut pool) */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpios.h"
alpar@9 26
alpar@9 27 /***********************************************************************
alpar@9 28 * NAME
alpar@9 29 *
alpar@9 30 * ios_process_cuts - process cuts stored in the local cut pool
alpar@9 31 *
alpar@9 32 * SYNOPSIS
alpar@9 33 *
alpar@9 34 * #include "glpios.h"
alpar@9 35 * void ios_process_cuts(glp_tree *T);
alpar@9 36 *
alpar@9 37 * DESCRIPTION
alpar@9 38 *
alpar@9 39 * The routine ios_process_cuts analyzes each cut currently stored in
alpar@9 40 * the local cut pool, which must be non-empty, and either adds the cut
alpar@9 41 * to the current subproblem or just discards it. All cuts are assumed
alpar@9 42 * to be locally valid. On exit the local cut pool remains unchanged.
alpar@9 43 *
alpar@9 44 * REFERENCES
alpar@9 45 *
alpar@9 46 * 1. E.Balas, S.Ceria, G.Cornuejols, "Mixed 0-1 Programming by
alpar@9 47 * Lift-and-Project in a Branch-and-Cut Framework", Management Sc.,
alpar@9 48 * 42 (1996) 1229-1246.
alpar@9 49 *
alpar@9 50 * 2. G.Andreello, A.Caprara, and M.Fischetti, "Embedding Cuts in
alpar@9 51 * a Branch&Cut Framework: a Computational Study with {0,1/2}-Cuts",
alpar@9 52 * Preliminary Draft, October 28, 2003, pp.6-8. */
alpar@9 53
alpar@9 54 struct info
alpar@9 55 { /* estimated cut efficiency */
alpar@9 56 IOSCUT *cut;
alpar@9 57 /* pointer to cut in the cut pool */
alpar@9 58 char flag;
alpar@9 59 /* if this flag is set, the cut is included into the current
alpar@9 60 subproblem */
alpar@9 61 double eff;
alpar@9 62 /* cut efficacy (normalized residual) */
alpar@9 63 double deg;
alpar@9 64 /* lower bound to objective degradation */
alpar@9 65 };
alpar@9 66
alpar@9 67 static int fcmp(const void *arg1, const void *arg2)
alpar@9 68 { const struct info *info1 = arg1, *info2 = arg2;
alpar@9 69 if (info1->deg == 0.0 && info2->deg == 0.0)
alpar@9 70 { if (info1->eff > info2->eff) return -1;
alpar@9 71 if (info1->eff < info2->eff) return +1;
alpar@9 72 }
alpar@9 73 else
alpar@9 74 { if (info1->deg > info2->deg) return -1;
alpar@9 75 if (info1->deg < info2->deg) return +1;
alpar@9 76 }
alpar@9 77 return 0;
alpar@9 78 }
alpar@9 79
alpar@9 80 static double parallel(IOSCUT *a, IOSCUT *b, double work[]);
alpar@9 81
alpar@9 82 void ios_process_cuts(glp_tree *T)
alpar@9 83 { IOSPOOL *pool;
alpar@9 84 IOSCUT *cut;
alpar@9 85 IOSAIJ *aij;
alpar@9 86 struct info *info;
alpar@9 87 int k, kk, max_cuts, len, ret, *ind;
alpar@9 88 double *val, *work;
alpar@9 89 /* the current subproblem must exist */
alpar@9 90 xassert(T->curr != NULL);
alpar@9 91 /* the pool must exist and be non-empty */
alpar@9 92 pool = T->local;
alpar@9 93 xassert(pool != NULL);
alpar@9 94 xassert(pool->size > 0);
alpar@9 95 /* allocate working arrays */
alpar@9 96 info = xcalloc(1+pool->size, sizeof(struct info));
alpar@9 97 ind = xcalloc(1+T->n, sizeof(int));
alpar@9 98 val = xcalloc(1+T->n, sizeof(double));
alpar@9 99 work = xcalloc(1+T->n, sizeof(double));
alpar@9 100 for (k = 1; k <= T->n; k++) work[k] = 0.0;
alpar@9 101 /* build the list of cuts stored in the cut pool */
alpar@9 102 for (k = 0, cut = pool->head; cut != NULL; cut = cut->next)
alpar@9 103 k++, info[k].cut = cut, info[k].flag = 0;
alpar@9 104 xassert(k == pool->size);
alpar@9 105 /* estimate efficiency of all cuts in the cut pool */
alpar@9 106 for (k = 1; k <= pool->size; k++)
alpar@9 107 { double temp, dy, dz;
alpar@9 108 cut = info[k].cut;
alpar@9 109 /* build the vector of cut coefficients and compute its
alpar@9 110 Euclidean norm */
alpar@9 111 len = 0; temp = 0.0;
alpar@9 112 for (aij = cut->ptr; aij != NULL; aij = aij->next)
alpar@9 113 { xassert(1 <= aij->j && aij->j <= T->n);
alpar@9 114 len++, ind[len] = aij->j, val[len] = aij->val;
alpar@9 115 temp += aij->val * aij->val;
alpar@9 116 }
alpar@9 117 if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON;
alpar@9 118 /* transform the cut to express it only through non-basic
alpar@9 119 (auxiliary and structural) variables */
alpar@9 120 len = glp_transform_row(T->mip, len, ind, val);
alpar@9 121 /* determine change in the cut value and in the objective
alpar@9 122 value for the adjacent basis by simulating one step of the
alpar@9 123 dual simplex */
alpar@9 124 ret = _glp_analyze_row(T->mip, len, ind, val, cut->type,
alpar@9 125 cut->rhs, 1e-9, NULL, NULL, NULL, NULL, &dy, &dz);
alpar@9 126 /* determine normalized residual and lower bound to objective
alpar@9 127 degradation */
alpar@9 128 if (ret == 0)
alpar@9 129 { info[k].eff = fabs(dy) / sqrt(temp);
alpar@9 130 /* if some reduced costs violates (slightly) their zero
alpar@9 131 bounds (i.e. have wrong signs) due to round-off errors,
alpar@9 132 dz also may have wrong sign being close to zero */
alpar@9 133 if (T->mip->dir == GLP_MIN)
alpar@9 134 { if (dz < 0.0) dz = 0.0;
alpar@9 135 info[k].deg = + dz;
alpar@9 136 }
alpar@9 137 else /* GLP_MAX */
alpar@9 138 { if (dz > 0.0) dz = 0.0;
alpar@9 139 info[k].deg = - dz;
alpar@9 140 }
alpar@9 141 }
alpar@9 142 else if (ret == 1)
alpar@9 143 { /* the constraint is not violated at the current point */
alpar@9 144 info[k].eff = info[k].deg = 0.0;
alpar@9 145 }
alpar@9 146 else if (ret == 2)
alpar@9 147 { /* no dual feasible adjacent basis exists */
alpar@9 148 info[k].eff = 1.0;
alpar@9 149 info[k].deg = DBL_MAX;
alpar@9 150 }
alpar@9 151 else
alpar@9 152 xassert(ret != ret);
alpar@9 153 /* if the degradation is too small, just ignore it */
alpar@9 154 if (info[k].deg < 0.01) info[k].deg = 0.0;
alpar@9 155 }
alpar@9 156 /* sort the list of cuts by decreasing objective degradation and
alpar@9 157 then by decreasing efficacy */
alpar@9 158 qsort(&info[1], pool->size, sizeof(struct info), fcmp);
alpar@9 159 /* only first (most efficient) max_cuts in the list are qualified
alpar@9 160 as candidates to be added to the current subproblem */
alpar@9 161 max_cuts = (T->curr->level == 0 ? 90 : 10);
alpar@9 162 if (max_cuts > pool->size) max_cuts = pool->size;
alpar@9 163 /* add cuts to the current subproblem */
alpar@9 164 #if 0
alpar@9 165 xprintf("*** adding cuts ***\n");
alpar@9 166 #endif
alpar@9 167 for (k = 1; k <= max_cuts; k++)
alpar@9 168 { int i, len;
alpar@9 169 /* if this cut seems to be inefficient, skip it */
alpar@9 170 if (info[k].deg < 0.01 && info[k].eff < 0.01) continue;
alpar@9 171 /* if the angle between this cut and every other cut included
alpar@9 172 in the current subproblem is small, skip this cut */
alpar@9 173 for (kk = 1; kk < k; kk++)
alpar@9 174 { if (info[kk].flag)
alpar@9 175 { if (parallel(info[k].cut, info[kk].cut, work) > 0.90)
alpar@9 176 break;
alpar@9 177 }
alpar@9 178 }
alpar@9 179 if (kk < k) continue;
alpar@9 180 /* add this cut to the current subproblem */
alpar@9 181 #if 0
alpar@9 182 xprintf("eff = %g; deg = %g\n", info[k].eff, info[k].deg);
alpar@9 183 #endif
alpar@9 184 cut = info[k].cut, info[k].flag = 1;
alpar@9 185 i = glp_add_rows(T->mip, 1);
alpar@9 186 if (cut->name != NULL)
alpar@9 187 glp_set_row_name(T->mip, i, cut->name);
alpar@9 188 xassert(T->mip->row[i]->origin == GLP_RF_CUT);
alpar@9 189 T->mip->row[i]->klass = cut->klass;
alpar@9 190 len = 0;
alpar@9 191 for (aij = cut->ptr; aij != NULL; aij = aij->next)
alpar@9 192 len++, ind[len] = aij->j, val[len] = aij->val;
alpar@9 193 glp_set_mat_row(T->mip, i, len, ind, val);
alpar@9 194 xassert(cut->type == GLP_LO || cut->type == GLP_UP);
alpar@9 195 glp_set_row_bnds(T->mip, i, cut->type, cut->rhs, cut->rhs);
alpar@9 196 }
alpar@9 197 /* free working arrays */
alpar@9 198 xfree(info);
alpar@9 199 xfree(ind);
alpar@9 200 xfree(val);
alpar@9 201 xfree(work);
alpar@9 202 return;
alpar@9 203 }
alpar@9 204
alpar@9 205 #if 0
alpar@9 206 /***********************************************************************
alpar@9 207 * Given a cut a * x >= b (<= b) the routine efficacy computes the cut
alpar@9 208 * efficacy as follows:
alpar@9 209 *
alpar@9 210 * eff = d * (a * x~ - b) / ||a||,
alpar@9 211 *
alpar@9 212 * where d is -1 (in case of '>= b') or +1 (in case of '<= b'), x~ is
alpar@9 213 * the vector of values of structural variables in optimal solution to
alpar@9 214 * LP relaxation of the current subproblem, ||a|| is the Euclidean norm
alpar@9 215 * of the vector of cut coefficients.
alpar@9 216 *
alpar@9 217 * If the cut is violated at point x~, the efficacy eff is positive,
alpar@9 218 * and its value is the Euclidean distance between x~ and the cut plane
alpar@9 219 * a * x = b in the space of structural variables.
alpar@9 220 *
alpar@9 221 * Following geometrical intuition, it is quite natural to consider
alpar@9 222 * this distance as a first-order measure of the expected efficacy of
alpar@9 223 * the cut: the larger the distance the better the cut [1]. */
alpar@9 224
alpar@9 225 static double efficacy(glp_tree *T, IOSCUT *cut)
alpar@9 226 { glp_prob *mip = T->mip;
alpar@9 227 IOSAIJ *aij;
alpar@9 228 double s = 0.0, t = 0.0, temp;
alpar@9 229 for (aij = cut->ptr; aij != NULL; aij = aij->next)
alpar@9 230 { xassert(1 <= aij->j && aij->j <= mip->n);
alpar@9 231 s += aij->val * mip->col[aij->j]->prim;
alpar@9 232 t += aij->val * aij->val;
alpar@9 233 }
alpar@9 234 temp = sqrt(t);
alpar@9 235 if (temp < DBL_EPSILON) temp = DBL_EPSILON;
alpar@9 236 if (cut->type == GLP_LO)
alpar@9 237 temp = (s >= cut->rhs ? 0.0 : (cut->rhs - s) / temp);
alpar@9 238 else if (cut->type == GLP_UP)
alpar@9 239 temp = (s <= cut->rhs ? 0.0 : (s - cut->rhs) / temp);
alpar@9 240 else
alpar@9 241 xassert(cut != cut);
alpar@9 242 return temp;
alpar@9 243 }
alpar@9 244 #endif
alpar@9 245
alpar@9 246 /***********************************************************************
alpar@9 247 * Given two cuts a1 * x >= b1 (<= b1) and a2 * x >= b2 (<= b2) the
alpar@9 248 * routine parallel computes the cosine of angle between the cut planes
alpar@9 249 * a1 * x = b1 and a2 * x = b2 (which is the acute angle between two
alpar@9 250 * normals to these planes) in the space of structural variables as
alpar@9 251 * follows:
alpar@9 252 *
alpar@9 253 * cos phi = (a1' * a2) / (||a1|| * ||a2||),
alpar@9 254 *
alpar@9 255 * where (a1' * a2) is a dot product of vectors of cut coefficients,
alpar@9 256 * ||a1|| and ||a2|| are Euclidean norms of vectors a1 and a2.
alpar@9 257 *
alpar@9 258 * Note that requirement cos phi = 0 forces the cuts to be orthogonal,
alpar@9 259 * i.e. with disjoint support, while requirement cos phi <= 0.999 means
alpar@9 260 * only avoiding duplicate (parallel) cuts [1]. */
alpar@9 261
alpar@9 262 static double parallel(IOSCUT *a, IOSCUT *b, double work[])
alpar@9 263 { IOSAIJ *aij;
alpar@9 264 double s = 0.0, sa = 0.0, sb = 0.0, temp;
alpar@9 265 for (aij = a->ptr; aij != NULL; aij = aij->next)
alpar@9 266 { work[aij->j] = aij->val;
alpar@9 267 sa += aij->val * aij->val;
alpar@9 268 }
alpar@9 269 for (aij = b->ptr; aij != NULL; aij = aij->next)
alpar@9 270 { s += work[aij->j] * aij->val;
alpar@9 271 sb += aij->val * aij->val;
alpar@9 272 }
alpar@9 273 for (aij = a->ptr; aij != NULL; aij = aij->next)
alpar@9 274 work[aij->j] = 0.0;
alpar@9 275 temp = sqrt(sa) * sqrt(sb);
alpar@9 276 if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON;
alpar@9 277 return s / temp;
alpar@9 278 }
alpar@9 279
alpar@9 280 /* eof */