lemon-project-template-glpk: 33de93886c88 deps/glpk/src/glpios11.c

lemon-project-template-glpk

view deps/glpk/src/glpios11.c @ 9:33de93886c88

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

line source

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

3 /***********************************************************************

4 * This code is part of GLPK (GNU Linear Programming Kit).

5 *

7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,

9 * E-mail: <mao@gnu.org>.

10 *

11 * GLPK is free software: you can redistribute it and/or modify it

12 * under the terms of the GNU General Public License as published by

13 * the Free Software Foundation, either version 3 of the License, or

14 * (at your option) any later version.

15 *

16 * GLPK is distributed in the hope that it will be useful, but WITHOUT

17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

19 * License for more details.

20 *

21 * You should have received a copy of the GNU General Public License

22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.

23 ***********************************************************************/

25 #include "glpios.h"

27 /***********************************************************************

28 * NAME

29 *

30 * ios_process_cuts - process cuts stored in the local cut pool

31 *

32 * SYNOPSIS

33 *

34 * #include "glpios.h"

35 * void ios_process_cuts(glp_tree *T);

36 *

37 * DESCRIPTION

38 *

39 * The routine ios_process_cuts analyzes each cut currently stored in

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

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

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

43 *

44 * REFERENCES

45 *

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

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

48 * 42 (1996) 1229-1246.

49 *

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

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

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

54 struct info

55 { /* estimated cut efficiency */

56 IOSCUT *cut;

57 /* pointer to cut in the cut pool */

58 char flag;

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

60 subproblem */

61 double eff;

62 /* cut efficacy (normalized residual) */

63 double deg;

64 /* lower bound to objective degradation */

65 };

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

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

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

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

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

72 }

73 else

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

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

76 }

77 return 0;

78 }

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

82 void ios_process_cuts(glp_tree *T)

83 { IOSPOOL *pool;

84 IOSCUT *cut;

85 IOSAIJ *aij;

86 struct info *info;

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

88 double *val, *work;

89 /* the current subproblem must exist */

90 xassert(T->curr != NULL);

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

92 pool = T->local;

93 xassert(pool != NULL);

94 xassert(pool->size > 0);

95 /* allocate working arrays */

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

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

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

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

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

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

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

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

104 xassert(k == pool->size);

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

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

107 { double temp, dy, dz;

108 cut = info[k].cut;

109 /* build the vector of cut coefficients and compute its

110 Euclidean norm */

111 len = 0; temp = 0.0;

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

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

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

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

116 }

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

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

119 (auxiliary and structural) variables */

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

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

122 value for the adjacent basis by simulating one step of the

123 dual simplex */

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

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

126 /* determine normalized residual and lower bound to objective

127 degradation */

128 if (ret == 0)

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

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

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

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

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

134 { if (dz < 0.0) dz = 0.0;

135 info[k].deg = + dz;

136 }

137 else /* GLP_MAX */

138 { if (dz > 0.0) dz = 0.0;

139 info[k].deg = - dz;

140 }

141 }

142 else if (ret == 1)

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

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

145 }

146 else if (ret == 2)

147 { /* no dual feasible adjacent basis exists */

148 info[k].eff = 1.0;

149 info[k].deg = DBL_MAX;

150 }

151 else

152 xassert(ret != ret);

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

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

155 }

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

157 then by decreasing efficacy */

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

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

160 as candidates to be added to the current subproblem */

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

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

163 /* add cuts to the current subproblem */

164 #if 0

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

166 #endif

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

168 { int i, len;

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

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

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

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

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

174 { if (info[kk].flag)

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

176 break;

177 }

178 }

179 if (kk < k) continue;

180 /* add this cut to the current subproblem */

181 #if 0

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

183 #endif

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

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

186 if (cut->name != NULL)

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

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

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

190 len = 0;

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

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

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

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

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

196 }

197 /* free working arrays */

198 xfree(info);

199 xfree(ind);

200 xfree(val);

201 xfree(work);

202 return;

203 }

204

205 #if 0

206 /***********************************************************************

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

208 * efficacy as follows:

209 *

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

211 *

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

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

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

215 * of the vector of cut coefficients.

216 *

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

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

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

220 *

221 * Following geometrical intuition, it is quite natural to consider

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

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

224

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

226 { glp_prob *mip = T->mip;

227 IOSAIJ *aij;

228 double s = 0.0, t = 0.0, temp;

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

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

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

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

233 }

234 temp = sqrt(t);

235 if (temp < DBL_EPSILON) temp = DBL_EPSILON;

236 if (cut->type == GLP_LO)

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

238 else if (cut->type == GLP_UP)

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

240 else

241 xassert(cut != cut);

242 return temp;

243 }

244 #endif

245

246 /***********************************************************************

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

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

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

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

251 * follows:

252 *

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

254 *

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

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

257 *

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

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

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

261

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

263 { IOSAIJ *aij;

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

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

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

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

268 }

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

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

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

272 }

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

274 work[aij->j] = 0.0;

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

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

277 return s / temp;

278 }

279

280 /* eof */