1 | /* glpios11.c (process cuts stored in the local cut pool) */ |
---|
2 | |
---|
3 | /*********************************************************************** |
---|
4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
5 | * |
---|
6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
7 | * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
---|
8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
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 | ***********************************************************************/ |
---|
24 | |
---|
25 | #include "glpios.h" |
---|
26 | |
---|
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. */ |
---|
53 | |
---|
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 | }; |
---|
66 | |
---|
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 | } |
---|
79 | |
---|
80 | static double parallel(IOSCUT *a, IOSCUT *b, double work[]); |
---|
81 | |
---|
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 */ |
---|