|
1 /* glpios06.c (MIR cut generator) */ |
|
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 #define _MIR_DEBUG 0 |
|
28 |
|
29 #define MAXAGGR 5 |
|
30 /* maximal number of rows which can be aggregated */ |
|
31 |
|
32 struct MIR |
|
33 { /* MIR cut generator working area */ |
|
34 /*--------------------------------------------------------------*/ |
|
35 /* global information valid for the root subproblem */ |
|
36 int m; |
|
37 /* number of rows (in the root subproblem) */ |
|
38 int n; |
|
39 /* number of columns */ |
|
40 char *skip; /* char skip[1+m]; */ |
|
41 /* skip[i], 1 <= i <= m, is a flag that means that row i should |
|
42 not be used because (1) it is not suitable, or (2) because it |
|
43 has been used in the aggregated constraint */ |
|
44 char *isint; /* char isint[1+m+n]; */ |
|
45 /* isint[k], 1 <= k <= m+n, is a flag that means that variable |
|
46 x[k] is integer (otherwise, continuous) */ |
|
47 double *lb; /* double lb[1+m+n]; */ |
|
48 /* lb[k], 1 <= k <= m+n, is lower bound of x[k]; -DBL_MAX means |
|
49 that x[k] has no lower bound */ |
|
50 int *vlb; /* int vlb[1+m+n]; */ |
|
51 /* vlb[k] = k', 1 <= k <= m+n, is the number of integer variable, |
|
52 which defines variable lower bound x[k] >= lb[k] * x[k']; zero |
|
53 means that x[k] has simple lower bound */ |
|
54 double *ub; /* double ub[1+m+n]; */ |
|
55 /* ub[k], 1 <= k <= m+n, is upper bound of x[k]; +DBL_MAX means |
|
56 that x[k] has no upper bound */ |
|
57 int *vub; /* int vub[1+m+n]; */ |
|
58 /* vub[k] = k', 1 <= k <= m+n, is the number of integer variable, |
|
59 which defines variable upper bound x[k] <= ub[k] * x[k']; zero |
|
60 means that x[k] has simple upper bound */ |
|
61 /*--------------------------------------------------------------*/ |
|
62 /* current (fractional) point to be separated */ |
|
63 double *x; /* double x[1+m+n]; */ |
|
64 /* x[k] is current value of auxiliary (1 <= k <= m) or structural |
|
65 (m+1 <= k <= m+n) variable */ |
|
66 /*--------------------------------------------------------------*/ |
|
67 /* aggregated constraint sum a[k] * x[k] = b, which is a linear |
|
68 combination of original constraints transformed to equalities |
|
69 by introducing auxiliary variables */ |
|
70 int agg_cnt; |
|
71 /* number of rows (original constraints) used to build aggregated |
|
72 constraint, 1 <= agg_cnt <= MAXAGGR */ |
|
73 int *agg_row; /* int agg_row[1+MAXAGGR]; */ |
|
74 /* agg_row[k], 1 <= k <= agg_cnt, is the row number used to build |
|
75 aggregated constraint */ |
|
76 IOSVEC *agg_vec; /* IOSVEC agg_vec[1:m+n]; */ |
|
77 /* sparse vector of aggregated constraint coefficients, a[k] */ |
|
78 double agg_rhs; |
|
79 /* right-hand side of the aggregated constraint, b */ |
|
80 /*--------------------------------------------------------------*/ |
|
81 /* bound substitution flags for modified constraint */ |
|
82 char *subst; /* char subst[1+m+n]; */ |
|
83 /* subst[k], 1 <= k <= m+n, is a bound substitution flag used for |
|
84 variable x[k]: |
|
85 '?' - x[k] is missing in modified constraint |
|
86 'L' - x[k] = (lower bound) + x'[k] |
|
87 'U' - x[k] = (upper bound) - x'[k] */ |
|
88 /*--------------------------------------------------------------*/ |
|
89 /* modified constraint sum a'[k] * x'[k] = b', where x'[k] >= 0, |
|
90 derived from aggregated constraint by substituting bounds; |
|
91 note that due to substitution of variable bounds there may be |
|
92 additional terms in the modified constraint */ |
|
93 IOSVEC *mod_vec; /* IOSVEC mod_vec[1:m+n]; */ |
|
94 /* sparse vector of modified constraint coefficients, a'[k] */ |
|
95 double mod_rhs; |
|
96 /* right-hand side of the modified constraint, b' */ |
|
97 /*--------------------------------------------------------------*/ |
|
98 /* cutting plane sum alpha[k] * x[k] <= beta */ |
|
99 IOSVEC *cut_vec; /* IOSVEC cut_vec[1:m+n]; */ |
|
100 /* sparse vector of cutting plane coefficients, alpha[k] */ |
|
101 double cut_rhs; |
|
102 /* right-hand size of the cutting plane, beta */ |
|
103 }; |
|
104 |
|
105 /*********************************************************************** |
|
106 * NAME |
|
107 * |
|
108 * ios_mir_init - initialize MIR cut generator |
|
109 * |
|
110 * SYNOPSIS |
|
111 * |
|
112 * #include "glpios.h" |
|
113 * void *ios_mir_init(glp_tree *tree); |
|
114 * |
|
115 * DESCRIPTION |
|
116 * |
|
117 * The routine ios_mir_init initializes the MIR cut generator assuming |
|
118 * that the current subproblem is the root subproblem. |
|
119 * |
|
120 * RETURNS |
|
121 * |
|
122 * The routine ios_mir_init returns a pointer to the MIR cut generator |
|
123 * working area. */ |
|
124 |
|
125 static void set_row_attrib(glp_tree *tree, struct MIR *mir) |
|
126 { /* set global row attributes */ |
|
127 glp_prob *mip = tree->mip; |
|
128 int m = mir->m; |
|
129 int k; |
|
130 for (k = 1; k <= m; k++) |
|
131 { GLPROW *row = mip->row[k]; |
|
132 mir->skip[k] = 0; |
|
133 mir->isint[k] = 0; |
|
134 switch (row->type) |
|
135 { case GLP_FR: |
|
136 mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break; |
|
137 case GLP_LO: |
|
138 mir->lb[k] = row->lb, mir->ub[k] = +DBL_MAX; break; |
|
139 case GLP_UP: |
|
140 mir->lb[k] = -DBL_MAX, mir->ub[k] = row->ub; break; |
|
141 case GLP_DB: |
|
142 mir->lb[k] = row->lb, mir->ub[k] = row->ub; break; |
|
143 case GLP_FX: |
|
144 mir->lb[k] = mir->ub[k] = row->lb; break; |
|
145 default: |
|
146 xassert(row != row); |
|
147 } |
|
148 mir->vlb[k] = mir->vub[k] = 0; |
|
149 } |
|
150 return; |
|
151 } |
|
152 |
|
153 static void set_col_attrib(glp_tree *tree, struct MIR *mir) |
|
154 { /* set global column attributes */ |
|
155 glp_prob *mip = tree->mip; |
|
156 int m = mir->m; |
|
157 int n = mir->n; |
|
158 int k; |
|
159 for (k = m+1; k <= m+n; k++) |
|
160 { GLPCOL *col = mip->col[k-m]; |
|
161 switch (col->kind) |
|
162 { case GLP_CV: |
|
163 mir->isint[k] = 0; break; |
|
164 case GLP_IV: |
|
165 mir->isint[k] = 1; break; |
|
166 default: |
|
167 xassert(col != col); |
|
168 } |
|
169 switch (col->type) |
|
170 { case GLP_FR: |
|
171 mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break; |
|
172 case GLP_LO: |
|
173 mir->lb[k] = col->lb, mir->ub[k] = +DBL_MAX; break; |
|
174 case GLP_UP: |
|
175 mir->lb[k] = -DBL_MAX, mir->ub[k] = col->ub; break; |
|
176 case GLP_DB: |
|
177 mir->lb[k] = col->lb, mir->ub[k] = col->ub; break; |
|
178 case GLP_FX: |
|
179 mir->lb[k] = mir->ub[k] = col->lb; break; |
|
180 default: |
|
181 xassert(col != col); |
|
182 } |
|
183 mir->vlb[k] = mir->vub[k] = 0; |
|
184 } |
|
185 return; |
|
186 } |
|
187 |
|
188 static void set_var_bounds(glp_tree *tree, struct MIR *mir) |
|
189 { /* set variable bounds */ |
|
190 glp_prob *mip = tree->mip; |
|
191 int m = mir->m; |
|
192 GLPAIJ *aij; |
|
193 int i, k1, k2; |
|
194 double a1, a2; |
|
195 for (i = 1; i <= m; i++) |
|
196 { /* we need the row to be '>= 0' or '<= 0' */ |
|
197 if (!(mir->lb[i] == 0.0 && mir->ub[i] == +DBL_MAX || |
|
198 mir->lb[i] == -DBL_MAX && mir->ub[i] == 0.0)) continue; |
|
199 /* take first term */ |
|
200 aij = mip->row[i]->ptr; |
|
201 if (aij == NULL) continue; |
|
202 k1 = m + aij->col->j, a1 = aij->val; |
|
203 /* take second term */ |
|
204 aij = aij->r_next; |
|
205 if (aij == NULL) continue; |
|
206 k2 = m + aij->col->j, a2 = aij->val; |
|
207 /* there must be only two terms */ |
|
208 if (aij->r_next != NULL) continue; |
|
209 /* interchange terms, if needed */ |
|
210 if (!mir->isint[k1] && mir->isint[k2]) |
|
211 ; |
|
212 else if (mir->isint[k1] && !mir->isint[k2]) |
|
213 { k2 = k1, a2 = a1; |
|
214 k1 = m + aij->col->j, a1 = aij->val; |
|
215 } |
|
216 else |
|
217 { /* both terms are either continuous or integer */ |
|
218 continue; |
|
219 } |
|
220 /* x[k2] should be double-bounded */ |
|
221 if (mir->lb[k2] == -DBL_MAX || mir->ub[k2] == +DBL_MAX || |
|
222 mir->lb[k2] == mir->ub[k2]) continue; |
|
223 /* change signs, if necessary */ |
|
224 if (mir->ub[i] == 0.0) a1 = - a1, a2 = - a2; |
|
225 /* now the row has the form a1 * x1 + a2 * x2 >= 0, where x1 |
|
226 is continuous, x2 is integer */ |
|
227 if (a1 > 0.0) |
|
228 { /* x1 >= - (a2 / a1) * x2 */ |
|
229 if (mir->vlb[k1] == 0) |
|
230 { /* set variable lower bound for x1 */ |
|
231 mir->lb[k1] = - a2 / a1; |
|
232 mir->vlb[k1] = k2; |
|
233 /* the row should not be used */ |
|
234 mir->skip[i] = 1; |
|
235 } |
|
236 } |
|
237 else /* a1 < 0.0 */ |
|
238 { /* x1 <= - (a2 / a1) * x2 */ |
|
239 if (mir->vub[k1] == 0) |
|
240 { /* set variable upper bound for x1 */ |
|
241 mir->ub[k1] = - a2 / a1; |
|
242 mir->vub[k1] = k2; |
|
243 /* the row should not be used */ |
|
244 mir->skip[i] = 1; |
|
245 } |
|
246 } |
|
247 } |
|
248 return; |
|
249 } |
|
250 |
|
251 static void mark_useless_rows(glp_tree *tree, struct MIR *mir) |
|
252 { /* mark rows which should not be used */ |
|
253 glp_prob *mip = tree->mip; |
|
254 int m = mir->m; |
|
255 GLPAIJ *aij; |
|
256 int i, k, nv; |
|
257 for (i = 1; i <= m; i++) |
|
258 { /* free rows should not be used */ |
|
259 if (mir->lb[i] == -DBL_MAX && mir->ub[i] == +DBL_MAX) |
|
260 { mir->skip[i] = 1; |
|
261 continue; |
|
262 } |
|
263 nv = 0; |
|
264 for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next) |
|
265 { k = m + aij->col->j; |
|
266 /* rows with free variables should not be used */ |
|
267 if (mir->lb[k] == -DBL_MAX && mir->ub[k] == +DBL_MAX) |
|
268 { mir->skip[i] = 1; |
|
269 break; |
|
270 } |
|
271 /* rows with integer variables having infinite (lower or |
|
272 upper) bound should not be used */ |
|
273 if (mir->isint[k] && mir->lb[k] == -DBL_MAX || |
|
274 mir->isint[k] && mir->ub[k] == +DBL_MAX) |
|
275 { mir->skip[i] = 1; |
|
276 break; |
|
277 } |
|
278 /* count non-fixed variables */ |
|
279 if (!(mir->vlb[k] == 0 && mir->vub[k] == 0 && |
|
280 mir->lb[k] == mir->ub[k])) nv++; |
|
281 } |
|
282 /* rows with all variables fixed should not be used */ |
|
283 if (nv == 0) |
|
284 { mir->skip[i] = 1; |
|
285 continue; |
|
286 } |
|
287 } |
|
288 return; |
|
289 } |
|
290 |
|
291 void *ios_mir_init(glp_tree *tree) |
|
292 { /* initialize MIR cut generator */ |
|
293 glp_prob *mip = tree->mip; |
|
294 int m = mip->m; |
|
295 int n = mip->n; |
|
296 struct MIR *mir; |
|
297 #if _MIR_DEBUG |
|
298 xprintf("ios_mir_init: warning: debug mode enabled\n"); |
|
299 #endif |
|
300 /* allocate working area */ |
|
301 mir = xmalloc(sizeof(struct MIR)); |
|
302 mir->m = m; |
|
303 mir->n = n; |
|
304 mir->skip = xcalloc(1+m, sizeof(char)); |
|
305 mir->isint = xcalloc(1+m+n, sizeof(char)); |
|
306 mir->lb = xcalloc(1+m+n, sizeof(double)); |
|
307 mir->vlb = xcalloc(1+m+n, sizeof(int)); |
|
308 mir->ub = xcalloc(1+m+n, sizeof(double)); |
|
309 mir->vub = xcalloc(1+m+n, sizeof(int)); |
|
310 mir->x = xcalloc(1+m+n, sizeof(double)); |
|
311 mir->agg_row = xcalloc(1+MAXAGGR, sizeof(int)); |
|
312 mir->agg_vec = ios_create_vec(m+n); |
|
313 mir->subst = xcalloc(1+m+n, sizeof(char)); |
|
314 mir->mod_vec = ios_create_vec(m+n); |
|
315 mir->cut_vec = ios_create_vec(m+n); |
|
316 /* set global row attributes */ |
|
317 set_row_attrib(tree, mir); |
|
318 /* set global column attributes */ |
|
319 set_col_attrib(tree, mir); |
|
320 /* set variable bounds */ |
|
321 set_var_bounds(tree, mir); |
|
322 /* mark rows which should not be used */ |
|
323 mark_useless_rows(tree, mir); |
|
324 return mir; |
|
325 } |
|
326 |
|
327 /*********************************************************************** |
|
328 * NAME |
|
329 * |
|
330 * ios_mir_gen - generate MIR cuts |
|
331 * |
|
332 * SYNOPSIS |
|
333 * |
|
334 * #include "glpios.h" |
|
335 * void ios_mir_gen(glp_tree *tree, void *gen, IOSPOOL *pool); |
|
336 * |
|
337 * DESCRIPTION |
|
338 * |
|
339 * The routine ios_mir_gen generates MIR cuts for the current point and |
|
340 * adds them to the cut pool. */ |
|
341 |
|
342 static void get_current_point(glp_tree *tree, struct MIR *mir) |
|
343 { /* obtain current point */ |
|
344 glp_prob *mip = tree->mip; |
|
345 int m = mir->m; |
|
346 int n = mir->n; |
|
347 int k; |
|
348 for (k = 1; k <= m; k++) |
|
349 mir->x[k] = mip->row[k]->prim; |
|
350 for (k = m+1; k <= m+n; k++) |
|
351 mir->x[k] = mip->col[k-m]->prim; |
|
352 return; |
|
353 } |
|
354 |
|
355 #if _MIR_DEBUG |
|
356 static void check_current_point(struct MIR *mir) |
|
357 { /* check current point */ |
|
358 int m = mir->m; |
|
359 int n = mir->n; |
|
360 int k, kk; |
|
361 double lb, ub, eps; |
|
362 for (k = 1; k <= m+n; k++) |
|
363 { /* determine lower bound */ |
|
364 lb = mir->lb[k]; |
|
365 kk = mir->vlb[k]; |
|
366 if (kk != 0) |
|
367 { xassert(lb != -DBL_MAX); |
|
368 xassert(!mir->isint[k]); |
|
369 xassert(mir->isint[kk]); |
|
370 lb *= mir->x[kk]; |
|
371 } |
|
372 /* check lower bound */ |
|
373 if (lb != -DBL_MAX) |
|
374 { eps = 1e-6 * (1.0 + fabs(lb)); |
|
375 xassert(mir->x[k] >= lb - eps); |
|
376 } |
|
377 /* determine upper bound */ |
|
378 ub = mir->ub[k]; |
|
379 kk = mir->vub[k]; |
|
380 if (kk != 0) |
|
381 { xassert(ub != +DBL_MAX); |
|
382 xassert(!mir->isint[k]); |
|
383 xassert(mir->isint[kk]); |
|
384 ub *= mir->x[kk]; |
|
385 } |
|
386 /* check upper bound */ |
|
387 if (ub != +DBL_MAX) |
|
388 { eps = 1e-6 * (1.0 + fabs(ub)); |
|
389 xassert(mir->x[k] <= ub + eps); |
|
390 } |
|
391 } |
|
392 return; |
|
393 } |
|
394 #endif |
|
395 |
|
396 static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i) |
|
397 { /* use original i-th row as initial aggregated constraint */ |
|
398 glp_prob *mip = tree->mip; |
|
399 int m = mir->m; |
|
400 GLPAIJ *aij; |
|
401 xassert(1 <= i && i <= m); |
|
402 xassert(!mir->skip[i]); |
|
403 /* mark i-th row in order not to use it in the same aggregated |
|
404 constraint */ |
|
405 mir->skip[i] = 2; |
|
406 mir->agg_cnt = 1; |
|
407 mir->agg_row[1] = i; |
|
408 /* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary |
|
409 variable of row i, x[m+j] are structural variables */ |
|
410 ios_clear_vec(mir->agg_vec); |
|
411 ios_set_vj(mir->agg_vec, i, 1.0); |
|
412 for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next) |
|
413 ios_set_vj(mir->agg_vec, m + aij->col->j, - aij->val); |
|
414 mir->agg_rhs = 0.0; |
|
415 #if _MIR_DEBUG |
|
416 ios_check_vec(mir->agg_vec); |
|
417 #endif |
|
418 return; |
|
419 } |
|
420 |
|
421 #if _MIR_DEBUG |
|
422 static void check_agg_row(struct MIR *mir) |
|
423 { /* check aggregated constraint */ |
|
424 int m = mir->m; |
|
425 int n = mir->n; |
|
426 int j, k; |
|
427 double r, big; |
|
428 /* compute the residual r = sum a[k] * x[k] - b and determine |
|
429 big = max(1, |a[k]|, |b|) */ |
|
430 r = 0.0, big = 1.0; |
|
431 for (j = 1; j <= mir->agg_vec->nnz; j++) |
|
432 { k = mir->agg_vec->ind[j]; |
|
433 xassert(1 <= k && k <= m+n); |
|
434 r += mir->agg_vec->val[j] * mir->x[k]; |
|
435 if (big < fabs(mir->agg_vec->val[j])) |
|
436 big = fabs(mir->agg_vec->val[j]); |
|
437 } |
|
438 r -= mir->agg_rhs; |
|
439 if (big < fabs(mir->agg_rhs)) |
|
440 big = fabs(mir->agg_rhs); |
|
441 /* the residual must be close to zero */ |
|
442 xassert(fabs(r) <= 1e-6 * big); |
|
443 return; |
|
444 } |
|
445 #endif |
|
446 |
|
447 static void subst_fixed_vars(struct MIR *mir) |
|
448 { /* substitute fixed variables into aggregated constraint */ |
|
449 int m = mir->m; |
|
450 int n = mir->n; |
|
451 int j, k; |
|
452 for (j = 1; j <= mir->agg_vec->nnz; j++) |
|
453 { k = mir->agg_vec->ind[j]; |
|
454 xassert(1 <= k && k <= m+n); |
|
455 if (mir->vlb[k] == 0 && mir->vub[k] == 0 && |
|
456 mir->lb[k] == mir->ub[k]) |
|
457 { /* x[k] is fixed */ |
|
458 mir->agg_rhs -= mir->agg_vec->val[j] * mir->lb[k]; |
|
459 mir->agg_vec->val[j] = 0.0; |
|
460 } |
|
461 } |
|
462 /* remove terms corresponding to fixed variables */ |
|
463 ios_clean_vec(mir->agg_vec, DBL_EPSILON); |
|
464 #if _MIR_DEBUG |
|
465 ios_check_vec(mir->agg_vec); |
|
466 #endif |
|
467 return; |
|
468 } |
|
469 |
|
470 static void bound_subst_heur(struct MIR *mir) |
|
471 { /* bound substitution heuristic */ |
|
472 int m = mir->m; |
|
473 int n = mir->n; |
|
474 int j, k, kk; |
|
475 double d1, d2; |
|
476 for (j = 1; j <= mir->agg_vec->nnz; j++) |
|
477 { k = mir->agg_vec->ind[j]; |
|
478 xassert(1 <= k && k <= m+n); |
|
479 if (mir->isint[k]) continue; /* skip integer variable */ |
|
480 /* compute distance from x[k] to its lower bound */ |
|
481 kk = mir->vlb[k]; |
|
482 if (kk == 0) |
|
483 { if (mir->lb[k] == -DBL_MAX) |
|
484 d1 = DBL_MAX; |
|
485 else |
|
486 d1 = mir->x[k] - mir->lb[k]; |
|
487 } |
|
488 else |
|
489 { xassert(1 <= kk && kk <= m+n); |
|
490 xassert(mir->isint[kk]); |
|
491 xassert(mir->lb[k] != -DBL_MAX); |
|
492 d1 = mir->x[k] - mir->lb[k] * mir->x[kk]; |
|
493 } |
|
494 /* compute distance from x[k] to its upper bound */ |
|
495 kk = mir->vub[k]; |
|
496 if (kk == 0) |
|
497 { if (mir->vub[k] == +DBL_MAX) |
|
498 d2 = DBL_MAX; |
|
499 else |
|
500 d2 = mir->ub[k] - mir->x[k]; |
|
501 } |
|
502 else |
|
503 { xassert(1 <= kk && kk <= m+n); |
|
504 xassert(mir->isint[kk]); |
|
505 xassert(mir->ub[k] != +DBL_MAX); |
|
506 d2 = mir->ub[k] * mir->x[kk] - mir->x[k]; |
|
507 } |
|
508 /* x[k] cannot be free */ |
|
509 xassert(d1 != DBL_MAX || d2 != DBL_MAX); |
|
510 /* choose the bound which is closer to x[k] */ |
|
511 xassert(mir->subst[k] == '?'); |
|
512 if (d1 <= d2) |
|
513 mir->subst[k] = 'L'; |
|
514 else |
|
515 mir->subst[k] = 'U'; |
|
516 } |
|
517 return; |
|
518 } |
|
519 |
|
520 static void build_mod_row(struct MIR *mir) |
|
521 { /* substitute bounds and build modified constraint */ |
|
522 int m = mir->m; |
|
523 int n = mir->n; |
|
524 int j, jj, k, kk; |
|
525 /* initially modified constraint is aggregated constraint */ |
|
526 ios_copy_vec(mir->mod_vec, mir->agg_vec); |
|
527 mir->mod_rhs = mir->agg_rhs; |
|
528 #if _MIR_DEBUG |
|
529 ios_check_vec(mir->mod_vec); |
|
530 #endif |
|
531 /* substitute bounds for continuous variables; note that due to |
|
532 substitution of variable bounds additional terms may appear in |
|
533 modified constraint */ |
|
534 for (j = mir->mod_vec->nnz; j >= 1; j--) |
|
535 { k = mir->mod_vec->ind[j]; |
|
536 xassert(1 <= k && k <= m+n); |
|
537 if (mir->isint[k]) continue; /* skip integer variable */ |
|
538 if (mir->subst[k] == 'L') |
|
539 { /* x[k] = (lower bound) + x'[k] */ |
|
540 xassert(mir->lb[k] != -DBL_MAX); |
|
541 kk = mir->vlb[k]; |
|
542 if (kk == 0) |
|
543 { /* x[k] = lb[k] + x'[k] */ |
|
544 mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k]; |
|
545 } |
|
546 else |
|
547 { /* x[k] = lb[k] * x[kk] + x'[k] */ |
|
548 xassert(mir->isint[kk]); |
|
549 jj = mir->mod_vec->pos[kk]; |
|
550 if (jj == 0) |
|
551 { ios_set_vj(mir->mod_vec, kk, 1.0); |
|
552 jj = mir->mod_vec->pos[kk]; |
|
553 mir->mod_vec->val[jj] = 0.0; |
|
554 } |
|
555 mir->mod_vec->val[jj] += |
|
556 mir->mod_vec->val[j] * mir->lb[k]; |
|
557 } |
|
558 } |
|
559 else if (mir->subst[k] == 'U') |
|
560 { /* x[k] = (upper bound) - x'[k] */ |
|
561 xassert(mir->ub[k] != +DBL_MAX); |
|
562 kk = mir->vub[k]; |
|
563 if (kk == 0) |
|
564 { /* x[k] = ub[k] - x'[k] */ |
|
565 mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k]; |
|
566 } |
|
567 else |
|
568 { /* x[k] = ub[k] * x[kk] - x'[k] */ |
|
569 xassert(mir->isint[kk]); |
|
570 jj = mir->mod_vec->pos[kk]; |
|
571 if (jj == 0) |
|
572 { ios_set_vj(mir->mod_vec, kk, 1.0); |
|
573 jj = mir->mod_vec->pos[kk]; |
|
574 mir->mod_vec->val[jj] = 0.0; |
|
575 } |
|
576 mir->mod_vec->val[jj] += |
|
577 mir->mod_vec->val[j] * mir->ub[k]; |
|
578 } |
|
579 mir->mod_vec->val[j] = - mir->mod_vec->val[j]; |
|
580 } |
|
581 else |
|
582 xassert(k != k); |
|
583 } |
|
584 #if _MIR_DEBUG |
|
585 ios_check_vec(mir->mod_vec); |
|
586 #endif |
|
587 /* substitute bounds for integer variables */ |
|
588 for (j = 1; j <= mir->mod_vec->nnz; j++) |
|
589 { k = mir->mod_vec->ind[j]; |
|
590 xassert(1 <= k && k <= m+n); |
|
591 if (!mir->isint[k]) continue; /* skip continuous variable */ |
|
592 xassert(mir->subst[k] == '?'); |
|
593 xassert(mir->vlb[k] == 0 && mir->vub[k] == 0); |
|
594 xassert(mir->lb[k] != -DBL_MAX && mir->ub[k] != +DBL_MAX); |
|
595 if (fabs(mir->lb[k]) <= fabs(mir->ub[k])) |
|
596 { /* x[k] = lb[k] + x'[k] */ |
|
597 mir->subst[k] = 'L'; |
|
598 mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k]; |
|
599 } |
|
600 else |
|
601 { /* x[k] = ub[k] - x'[k] */ |
|
602 mir->subst[k] = 'U'; |
|
603 mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k]; |
|
604 mir->mod_vec->val[j] = - mir->mod_vec->val[j]; |
|
605 } |
|
606 } |
|
607 #if _MIR_DEBUG |
|
608 ios_check_vec(mir->mod_vec); |
|
609 #endif |
|
610 return; |
|
611 } |
|
612 |
|
613 #if _MIR_DEBUG |
|
614 static void check_mod_row(struct MIR *mir) |
|
615 { /* check modified constraint */ |
|
616 int m = mir->m; |
|
617 int n = mir->n; |
|
618 int j, k, kk; |
|
619 double r, big, x; |
|
620 /* compute the residual r = sum a'[k] * x'[k] - b' and determine |
|
621 big = max(1, |a[k]|, |b|) */ |
|
622 r = 0.0, big = 1.0; |
|
623 for (j = 1; j <= mir->mod_vec->nnz; j++) |
|
624 { k = mir->mod_vec->ind[j]; |
|
625 xassert(1 <= k && k <= m+n); |
|
626 if (mir->subst[k] == 'L') |
|
627 { /* x'[k] = x[k] - (lower bound) */ |
|
628 xassert(mir->lb[k] != -DBL_MAX); |
|
629 kk = mir->vlb[k]; |
|
630 if (kk == 0) |
|
631 x = mir->x[k] - mir->lb[k]; |
|
632 else |
|
633 x = mir->x[k] - mir->lb[k] * mir->x[kk]; |
|
634 } |
|
635 else if (mir->subst[k] == 'U') |
|
636 { /* x'[k] = (upper bound) - x[k] */ |
|
637 xassert(mir->ub[k] != +DBL_MAX); |
|
638 kk = mir->vub[k]; |
|
639 if (kk == 0) |
|
640 x = mir->ub[k] - mir->x[k]; |
|
641 else |
|
642 x = mir->ub[k] * mir->x[kk] - mir->x[k]; |
|
643 } |
|
644 else |
|
645 xassert(k != k); |
|
646 r += mir->mod_vec->val[j] * x; |
|
647 if (big < fabs(mir->mod_vec->val[j])) |
|
648 big = fabs(mir->mod_vec->val[j]); |
|
649 } |
|
650 r -= mir->mod_rhs; |
|
651 if (big < fabs(mir->mod_rhs)) |
|
652 big = fabs(mir->mod_rhs); |
|
653 /* the residual must be close to zero */ |
|
654 xassert(fabs(r) <= 1e-6 * big); |
|
655 return; |
|
656 } |
|
657 #endif |
|
658 |
|
659 /*********************************************************************** |
|
660 * mir_ineq - construct MIR inequality |
|
661 * |
|
662 * Given the single constraint mixed integer set |
|
663 * |
|
664 * |N| |
|
665 * X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s}, |
|
666 * + + j in N |
|
667 * |
|
668 * this routine constructs the mixed integer rounding (MIR) inequality |
|
669 * |
|
670 * sum alpha[j] * x[j] <= beta + gamma * s, |
|
671 * j in N |
|
672 * |
|
673 * which is valid for X. |
|
674 * |
|
675 * If the MIR inequality has been successfully constructed, the routine |
|
676 * returns zero. Otherwise, if b is close to nearest integer, there may |
|
677 * be numeric difficulties due to big coefficients; so in this case the |
|
678 * routine returns non-zero. */ |
|
679 |
|
680 static int mir_ineq(const int n, const double a[], const double b, |
|
681 double alpha[], double *beta, double *gamma) |
|
682 { int j; |
|
683 double f, t; |
|
684 if (fabs(b - floor(b + .5)) < 0.01) |
|
685 return 1; |
|
686 f = b - floor(b); |
|
687 for (j = 1; j <= n; j++) |
|
688 { t = (a[j] - floor(a[j])) - f; |
|
689 if (t <= 0.0) |
|
690 alpha[j] = floor(a[j]); |
|
691 else |
|
692 alpha[j] = floor(a[j]) + t / (1.0 - f); |
|
693 } |
|
694 *beta = floor(b); |
|
695 *gamma = 1.0 / (1.0 - f); |
|
696 return 0; |
|
697 } |
|
698 |
|
699 /*********************************************************************** |
|
700 * cmir_ineq - construct c-MIR inequality |
|
701 * |
|
702 * Given the mixed knapsack set |
|
703 * |
|
704 * MK |N| |
|
705 * X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s, |
|
706 * + + j in N |
|
707 * |
|
708 * x[j] <= u[j]}, |
|
709 * |
|
710 * a subset C of variables to be complemented, and a divisor delta > 0, |
|
711 * this routine constructs the complemented MIR (c-MIR) inequality |
|
712 * |
|
713 * sum alpha[j] * x[j] <= beta + gamma * s, |
|
714 * j in N |
|
715 * MK |
|
716 * which is valid for X . |
|
717 * |
|
718 * If the c-MIR inequality has been successfully constructed, the |
|
719 * routine returns zero. Otherwise, if there is a risk of numerical |
|
720 * difficulties due to big coefficients (see comments to the routine |
|
721 * mir_ineq), the routine cmir_ineq returns non-zero. */ |
|
722 |
|
723 static int cmir_ineq(const int n, const double a[], const double b, |
|
724 const double u[], const char cset[], const double delta, |
|
725 double alpha[], double *beta, double *gamma) |
|
726 { int j; |
|
727 double *aa, bb; |
|
728 aa = alpha, bb = b; |
|
729 for (j = 1; j <= n; j++) |
|
730 { aa[j] = a[j] / delta; |
|
731 if (cset[j]) |
|
732 aa[j] = - aa[j], bb -= a[j] * u[j]; |
|
733 } |
|
734 bb /= delta; |
|
735 if (mir_ineq(n, aa, bb, alpha, beta, gamma)) return 1; |
|
736 for (j = 1; j <= n; j++) |
|
737 { if (cset[j]) |
|
738 alpha[j] = - alpha[j], *beta += alpha[j] * u[j]; |
|
739 } |
|
740 *gamma /= delta; |
|
741 return 0; |
|
742 } |
|
743 |
|
744 /*********************************************************************** |
|
745 * cmir_sep - c-MIR separation heuristic |
|
746 * |
|
747 * Given the mixed knapsack set |
|
748 * |
|
749 * MK |N| |
|
750 * X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s, |
|
751 * + + j in N |
|
752 * |
|
753 * x[j] <= u[j]} |
|
754 * |
|
755 * * * |
|
756 * and a fractional point (x , s ), this routine tries to construct |
|
757 * c-MIR inequality |
|
758 * |
|
759 * sum alpha[j] * x[j] <= beta + gamma * s, |
|
760 * j in N |
|
761 * MK |
|
762 * which is valid for X and has (desirably maximal) violation at the |
|
763 * fractional point given. This is attained by choosing an appropriate |
|
764 * set C of variables to be complemented and a divisor delta > 0, which |
|
765 * together define corresponding c-MIR inequality. |
|
766 * |
|
767 * If a violated c-MIR inequality has been successfully constructed, |
|
768 * the routine returns its violation: |
|
769 * |
|
770 * * * |
|
771 * sum alpha[j] * x [j] - beta - gamma * s , |
|
772 * j in N |
|
773 * |
|
774 * which is positive. In case of failure the routine returns zero. */ |
|
775 |
|
776 struct vset { int j; double v; }; |
|
777 |
|
778 static int cmir_cmp(const void *p1, const void *p2) |
|
779 { const struct vset *v1 = p1, *v2 = p2; |
|
780 if (v1->v < v2->v) return -1; |
|
781 if (v1->v > v2->v) return +1; |
|
782 return 0; |
|
783 } |
|
784 |
|
785 static double cmir_sep(const int n, const double a[], const double b, |
|
786 const double u[], const double x[], const double s, |
|
787 double alpha[], double *beta, double *gamma) |
|
788 { int fail, j, k, nv, v; |
|
789 double delta, eps, d_try[1+3], r, r_best; |
|
790 char *cset; |
|
791 struct vset *vset; |
|
792 /* allocate working arrays */ |
|
793 cset = xcalloc(1+n, sizeof(char)); |
|
794 vset = xcalloc(1+n, sizeof(struct vset)); |
|
795 /* choose initial C */ |
|
796 for (j = 1; j <= n; j++) |
|
797 cset[j] = (char)(x[j] >= 0.5 * u[j]); |
|
798 /* choose initial delta */ |
|
799 r_best = delta = 0.0; |
|
800 for (j = 1; j <= n; j++) |
|
801 { xassert(a[j] != 0.0); |
|
802 /* if x[j] is close to its bounds, skip it */ |
|
803 eps = 1e-9 * (1.0 + fabs(u[j])); |
|
804 if (x[j] < eps || x[j] > u[j] - eps) continue; |
|
805 /* try delta = |a[j]| to construct c-MIR inequality */ |
|
806 fail = cmir_ineq(n, a, b, u, cset, fabs(a[j]), alpha, beta, |
|
807 gamma); |
|
808 if (fail) continue; |
|
809 /* compute violation */ |
|
810 r = - (*beta) - (*gamma) * s; |
|
811 for (k = 1; k <= n; k++) r += alpha[k] * x[k]; |
|
812 if (r_best < r) r_best = r, delta = fabs(a[j]); |
|
813 } |
|
814 if (r_best < 0.001) r_best = 0.0; |
|
815 if (r_best == 0.0) goto done; |
|
816 xassert(delta > 0.0); |
|
817 /* try to increase violation by dividing delta by 2, 4, and 8, |
|
818 respectively */ |
|
819 d_try[1] = delta / 2.0; |
|
820 d_try[2] = delta / 4.0; |
|
821 d_try[3] = delta / 8.0; |
|
822 for (j = 1; j <= 3; j++) |
|
823 { /* construct c-MIR inequality */ |
|
824 fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta, |
|
825 gamma); |
|
826 if (fail) continue; |
|
827 /* compute violation */ |
|
828 r = - (*beta) - (*gamma) * s; |
|
829 for (k = 1; k <= n; k++) r += alpha[k] * x[k]; |
|
830 if (r_best < r) r_best = r, delta = d_try[j]; |
|
831 } |
|
832 /* build subset of variables lying strictly between their bounds |
|
833 and order it by nondecreasing values of |x[j] - u[j]/2| */ |
|
834 nv = 0; |
|
835 for (j = 1; j <= n; j++) |
|
836 { /* if x[j] is close to its bounds, skip it */ |
|
837 eps = 1e-9 * (1.0 + fabs(u[j])); |
|
838 if (x[j] < eps || x[j] > u[j] - eps) continue; |
|
839 /* add x[j] to the subset */ |
|
840 nv++; |
|
841 vset[nv].j = j; |
|
842 vset[nv].v = fabs(x[j] - 0.5 * u[j]); |
|
843 } |
|
844 qsort(&vset[1], nv, sizeof(struct vset), cmir_cmp); |
|
845 /* try to increase violation by successively complementing each |
|
846 variable in the subset */ |
|
847 for (v = 1; v <= nv; v++) |
|
848 { j = vset[v].j; |
|
849 /* replace x[j] by its complement or vice versa */ |
|
850 cset[j] = (char)!cset[j]; |
|
851 /* construct c-MIR inequality */ |
|
852 fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma); |
|
853 /* restore the variable */ |
|
854 cset[j] = (char)!cset[j]; |
|
855 /* do not replace the variable in case of failure */ |
|
856 if (fail) continue; |
|
857 /* compute violation */ |
|
858 r = - (*beta) - (*gamma) * s; |
|
859 for (k = 1; k <= n; k++) r += alpha[k] * x[k]; |
|
860 if (r_best < r) r_best = r, cset[j] = (char)!cset[j]; |
|
861 } |
|
862 /* construct the best c-MIR inequality chosen */ |
|
863 fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma); |
|
864 xassert(!fail); |
|
865 done: /* free working arrays */ |
|
866 xfree(cset); |
|
867 xfree(vset); |
|
868 /* return to the calling routine */ |
|
869 return r_best; |
|
870 } |
|
871 |
|
872 static double generate(struct MIR *mir) |
|
873 { /* try to generate violated c-MIR cut for modified constraint */ |
|
874 int m = mir->m; |
|
875 int n = mir->n; |
|
876 int j, k, kk, nint; |
|
877 double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma; |
|
878 ios_copy_vec(mir->cut_vec, mir->mod_vec); |
|
879 mir->cut_rhs = mir->mod_rhs; |
|
880 /* remove small terms, which can appear due to substitution of |
|
881 variable bounds */ |
|
882 ios_clean_vec(mir->cut_vec, DBL_EPSILON); |
|
883 #if _MIR_DEBUG |
|
884 ios_check_vec(mir->cut_vec); |
|
885 #endif |
|
886 /* remove positive continuous terms to obtain MK relaxation */ |
|
887 for (j = 1; j <= mir->cut_vec->nnz; j++) |
|
888 { k = mir->cut_vec->ind[j]; |
|
889 xassert(1 <= k && k <= m+n); |
|
890 if (!mir->isint[k] && mir->cut_vec->val[j] > 0.0) |
|
891 mir->cut_vec->val[j] = 0.0; |
|
892 } |
|
893 ios_clean_vec(mir->cut_vec, 0.0); |
|
894 #if _MIR_DEBUG |
|
895 ios_check_vec(mir->cut_vec); |
|
896 #endif |
|
897 /* move integer terms to the beginning of the sparse vector and |
|
898 determine the number of integer variables */ |
|
899 nint = 0; |
|
900 for (j = 1; j <= mir->cut_vec->nnz; j++) |
|
901 { k = mir->cut_vec->ind[j]; |
|
902 xassert(1 <= k && k <= m+n); |
|
903 if (mir->isint[k]) |
|
904 { double temp; |
|
905 nint++; |
|
906 /* interchange elements [nint] and [j] */ |
|
907 kk = mir->cut_vec->ind[nint]; |
|
908 mir->cut_vec->pos[k] = nint; |
|
909 mir->cut_vec->pos[kk] = j; |
|
910 mir->cut_vec->ind[nint] = k; |
|
911 mir->cut_vec->ind[j] = kk; |
|
912 temp = mir->cut_vec->val[nint]; |
|
913 mir->cut_vec->val[nint] = mir->cut_vec->val[j]; |
|
914 mir->cut_vec->val[j] = temp; |
|
915 } |
|
916 } |
|
917 #if _MIR_DEBUG |
|
918 ios_check_vec(mir->cut_vec); |
|
919 #endif |
|
920 /* if there is no integer variable, nothing to generate */ |
|
921 if (nint == 0) goto done; |
|
922 /* allocate working arrays */ |
|
923 u = xcalloc(1+nint, sizeof(double)); |
|
924 x = xcalloc(1+nint, sizeof(double)); |
|
925 alpha = xcalloc(1+nint, sizeof(double)); |
|
926 /* determine u and x */ |
|
927 for (j = 1; j <= nint; j++) |
|
928 { k = mir->cut_vec->ind[j]; |
|
929 xassert(m+1 <= k && k <= m+n); |
|
930 xassert(mir->isint[k]); |
|
931 u[j] = mir->ub[k] - mir->lb[k]; |
|
932 xassert(u[j] >= 1.0); |
|
933 if (mir->subst[k] == 'L') |
|
934 x[j] = mir->x[k] - mir->lb[k]; |
|
935 else if (mir->subst[k] == 'U') |
|
936 x[j] = mir->ub[k] - mir->x[k]; |
|
937 else |
|
938 xassert(k != k); |
|
939 xassert(x[j] >= -0.001); |
|
940 if (x[j] < 0.0) x[j] = 0.0; |
|
941 } |
|
942 /* compute s = - sum of continuous terms */ |
|
943 s = 0.0; |
|
944 for (j = nint+1; j <= mir->cut_vec->nnz; j++) |
|
945 { double x; |
|
946 k = mir->cut_vec->ind[j]; |
|
947 xassert(1 <= k && k <= m+n); |
|
948 /* must be continuous */ |
|
949 xassert(!mir->isint[k]); |
|
950 if (mir->subst[k] == 'L') |
|
951 { xassert(mir->lb[k] != -DBL_MAX); |
|
952 kk = mir->vlb[k]; |
|
953 if (kk == 0) |
|
954 x = mir->x[k] - mir->lb[k]; |
|
955 else |
|
956 x = mir->x[k] - mir->lb[k] * mir->x[kk]; |
|
957 } |
|
958 else if (mir->subst[k] == 'U') |
|
959 { xassert(mir->ub[k] != +DBL_MAX); |
|
960 kk = mir->vub[k]; |
|
961 if (kk == 0) |
|
962 x = mir->ub[k] - mir->x[k]; |
|
963 else |
|
964 x = mir->ub[k] * mir->x[kk] - mir->x[k]; |
|
965 } |
|
966 else |
|
967 xassert(k != k); |
|
968 xassert(x >= -0.001); |
|
969 if (x < 0.0) x = 0.0; |
|
970 s -= mir->cut_vec->val[j] * x; |
|
971 } |
|
972 xassert(s >= 0.0); |
|
973 /* apply heuristic to obtain most violated c-MIR inequality */ |
|
974 b = mir->cut_rhs; |
|
975 r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha, |
|
976 &beta, &gamma); |
|
977 if (r_best == 0.0) goto skip; |
|
978 xassert(r_best > 0.0); |
|
979 /* convert to raw cut */ |
|
980 /* sum alpha[j] * x[j] <= beta + gamma * s */ |
|
981 for (j = 1; j <= nint; j++) |
|
982 mir->cut_vec->val[j] = alpha[j]; |
|
983 for (j = nint+1; j <= mir->cut_vec->nnz; j++) |
|
984 { k = mir->cut_vec->ind[j]; |
|
985 if (k <= m+n) mir->cut_vec->val[j] *= gamma; |
|
986 } |
|
987 mir->cut_rhs = beta; |
|
988 #if _MIR_DEBUG |
|
989 ios_check_vec(mir->cut_vec); |
|
990 #endif |
|
991 skip: /* free working arrays */ |
|
992 xfree(u); |
|
993 xfree(x); |
|
994 xfree(alpha); |
|
995 done: return r_best; |
|
996 } |
|
997 |
|
998 #if _MIR_DEBUG |
|
999 static void check_raw_cut(struct MIR *mir, double r_best) |
|
1000 { /* check raw cut before back bound substitution */ |
|
1001 int m = mir->m; |
|
1002 int n = mir->n; |
|
1003 int j, k, kk; |
|
1004 double r, big, x; |
|
1005 /* compute the residual r = sum a[k] * x[k] - b and determine |
|
1006 big = max(1, |a[k]|, |b|) */ |
|
1007 r = 0.0, big = 1.0; |
|
1008 for (j = 1; j <= mir->cut_vec->nnz; j++) |
|
1009 { k = mir->cut_vec->ind[j]; |
|
1010 xassert(1 <= k && k <= m+n); |
|
1011 if (mir->subst[k] == 'L') |
|
1012 { xassert(mir->lb[k] != -DBL_MAX); |
|
1013 kk = mir->vlb[k]; |
|
1014 if (kk == 0) |
|
1015 x = mir->x[k] - mir->lb[k]; |
|
1016 else |
|
1017 x = mir->x[k] - mir->lb[k] * mir->x[kk]; |
|
1018 } |
|
1019 else if (mir->subst[k] == 'U') |
|
1020 { xassert(mir->ub[k] != +DBL_MAX); |
|
1021 kk = mir->vub[k]; |
|
1022 if (kk == 0) |
|
1023 x = mir->ub[k] - mir->x[k]; |
|
1024 else |
|
1025 x = mir->ub[k] * mir->x[kk] - mir->x[k]; |
|
1026 } |
|
1027 else |
|
1028 xassert(k != k); |
|
1029 r += mir->cut_vec->val[j] * x; |
|
1030 if (big < fabs(mir->cut_vec->val[j])) |
|
1031 big = fabs(mir->cut_vec->val[j]); |
|
1032 } |
|
1033 r -= mir->cut_rhs; |
|
1034 if (big < fabs(mir->cut_rhs)) |
|
1035 big = fabs(mir->cut_rhs); |
|
1036 /* the residual must be close to r_best */ |
|
1037 xassert(fabs(r - r_best) <= 1e-6 * big); |
|
1038 return; |
|
1039 } |
|
1040 #endif |
|
1041 |
|
1042 static void back_subst(struct MIR *mir) |
|
1043 { /* back substitution of original bounds */ |
|
1044 int m = mir->m; |
|
1045 int n = mir->n; |
|
1046 int j, jj, k, kk; |
|
1047 /* at first, restore bounds of integer variables (because on |
|
1048 restoring variable bounds of continuous variables we need |
|
1049 original, not shifted, bounds of integer variables) */ |
|
1050 for (j = 1; j <= mir->cut_vec->nnz; j++) |
|
1051 { k = mir->cut_vec->ind[j]; |
|
1052 xassert(1 <= k && k <= m+n); |
|
1053 if (!mir->isint[k]) continue; /* skip continuous */ |
|
1054 if (mir->subst[k] == 'L') |
|
1055 { /* x'[k] = x[k] - lb[k] */ |
|
1056 xassert(mir->lb[k] != -DBL_MAX); |
|
1057 xassert(mir->vlb[k] == 0); |
|
1058 mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k]; |
|
1059 } |
|
1060 else if (mir->subst[k] == 'U') |
|
1061 { /* x'[k] = ub[k] - x[k] */ |
|
1062 xassert(mir->ub[k] != +DBL_MAX); |
|
1063 xassert(mir->vub[k] == 0); |
|
1064 mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k]; |
|
1065 mir->cut_vec->val[j] = - mir->cut_vec->val[j]; |
|
1066 } |
|
1067 else |
|
1068 xassert(k != k); |
|
1069 } |
|
1070 /* now restore bounds of continuous variables */ |
|
1071 for (j = 1; j <= mir->cut_vec->nnz; j++) |
|
1072 { k = mir->cut_vec->ind[j]; |
|
1073 xassert(1 <= k && k <= m+n); |
|
1074 if (mir->isint[k]) continue; /* skip integer */ |
|
1075 if (mir->subst[k] == 'L') |
|
1076 { /* x'[k] = x[k] - (lower bound) */ |
|
1077 xassert(mir->lb[k] != -DBL_MAX); |
|
1078 kk = mir->vlb[k]; |
|
1079 if (kk == 0) |
|
1080 { /* x'[k] = x[k] - lb[k] */ |
|
1081 mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k]; |
|
1082 } |
|
1083 else |
|
1084 { /* x'[k] = x[k] - lb[k] * x[kk] */ |
|
1085 jj = mir->cut_vec->pos[kk]; |
|
1086 #if 0 |
|
1087 xassert(jj != 0); |
|
1088 #else |
|
1089 if (jj == 0) |
|
1090 { ios_set_vj(mir->cut_vec, kk, 1.0); |
|
1091 jj = mir->cut_vec->pos[kk]; |
|
1092 xassert(jj != 0); |
|
1093 mir->cut_vec->val[jj] = 0.0; |
|
1094 } |
|
1095 #endif |
|
1096 mir->cut_vec->val[jj] -= mir->cut_vec->val[j] * |
|
1097 mir->lb[k]; |
|
1098 } |
|
1099 } |
|
1100 else if (mir->subst[k] == 'U') |
|
1101 { /* x'[k] = (upper bound) - x[k] */ |
|
1102 xassert(mir->ub[k] != +DBL_MAX); |
|
1103 kk = mir->vub[k]; |
|
1104 if (kk == 0) |
|
1105 { /* x'[k] = ub[k] - x[k] */ |
|
1106 mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k]; |
|
1107 } |
|
1108 else |
|
1109 { /* x'[k] = ub[k] * x[kk] - x[k] */ |
|
1110 jj = mir->cut_vec->pos[kk]; |
|
1111 if (jj == 0) |
|
1112 { ios_set_vj(mir->cut_vec, kk, 1.0); |
|
1113 jj = mir->cut_vec->pos[kk]; |
|
1114 xassert(jj != 0); |
|
1115 mir->cut_vec->val[jj] = 0.0; |
|
1116 } |
|
1117 mir->cut_vec->val[jj] += mir->cut_vec->val[j] * |
|
1118 mir->ub[k]; |
|
1119 } |
|
1120 mir->cut_vec->val[j] = - mir->cut_vec->val[j]; |
|
1121 } |
|
1122 else |
|
1123 xassert(k != k); |
|
1124 } |
|
1125 #if _MIR_DEBUG |
|
1126 ios_check_vec(mir->cut_vec); |
|
1127 #endif |
|
1128 return; |
|
1129 } |
|
1130 |
|
1131 #if _MIR_DEBUG |
|
1132 static void check_cut_row(struct MIR *mir, double r_best) |
|
1133 { /* check the cut after back bound substitution or elimination of |
|
1134 auxiliary variables */ |
|
1135 int m = mir->m; |
|
1136 int n = mir->n; |
|
1137 int j, k; |
|
1138 double r, big; |
|
1139 /* compute the residual r = sum a[k] * x[k] - b and determine |
|
1140 big = max(1, |a[k]|, |b|) */ |
|
1141 r = 0.0, big = 1.0; |
|
1142 for (j = 1; j <= mir->cut_vec->nnz; j++) |
|
1143 { k = mir->cut_vec->ind[j]; |
|
1144 xassert(1 <= k && k <= m+n); |
|
1145 r += mir->cut_vec->val[j] * mir->x[k]; |
|
1146 if (big < fabs(mir->cut_vec->val[j])) |
|
1147 big = fabs(mir->cut_vec->val[j]); |
|
1148 } |
|
1149 r -= mir->cut_rhs; |
|
1150 if (big < fabs(mir->cut_rhs)) |
|
1151 big = fabs(mir->cut_rhs); |
|
1152 /* the residual must be close to r_best */ |
|
1153 xassert(fabs(r - r_best) <= 1e-6 * big); |
|
1154 return; |
|
1155 } |
|
1156 #endif |
|
1157 |
|
1158 static void subst_aux_vars(glp_tree *tree, struct MIR *mir) |
|
1159 { /* final substitution to eliminate auxiliary variables */ |
|
1160 glp_prob *mip = tree->mip; |
|
1161 int m = mir->m; |
|
1162 int n = mir->n; |
|
1163 GLPAIJ *aij; |
|
1164 int j, k, kk, jj; |
|
1165 for (j = mir->cut_vec->nnz; j >= 1; j--) |
|
1166 { k = mir->cut_vec->ind[j]; |
|
1167 xassert(1 <= k && k <= m+n); |
|
1168 if (k > m) continue; /* skip structurals */ |
|
1169 for (aij = mip->row[k]->ptr; aij != NULL; aij = aij->r_next) |
|
1170 { kk = m + aij->col->j; /* structural */ |
|
1171 jj = mir->cut_vec->pos[kk]; |
|
1172 if (jj == 0) |
|
1173 { ios_set_vj(mir->cut_vec, kk, 1.0); |
|
1174 jj = mir->cut_vec->pos[kk]; |
|
1175 mir->cut_vec->val[jj] = 0.0; |
|
1176 } |
|
1177 mir->cut_vec->val[jj] += mir->cut_vec->val[j] * aij->val; |
|
1178 } |
|
1179 mir->cut_vec->val[j] = 0.0; |
|
1180 } |
|
1181 ios_clean_vec(mir->cut_vec, 0.0); |
|
1182 return; |
|
1183 } |
|
1184 |
|
1185 static void add_cut(glp_tree *tree, struct MIR *mir) |
|
1186 { /* add constructed cut inequality to the cut pool */ |
|
1187 int m = mir->m; |
|
1188 int n = mir->n; |
|
1189 int j, k, len; |
|
1190 int *ind = xcalloc(1+n, sizeof(int)); |
|
1191 double *val = xcalloc(1+n, sizeof(double)); |
|
1192 len = 0; |
|
1193 for (j = mir->cut_vec->nnz; j >= 1; j--) |
|
1194 { k = mir->cut_vec->ind[j]; |
|
1195 xassert(m+1 <= k && k <= m+n); |
|
1196 len++, ind[len] = k - m, val[len] = mir->cut_vec->val[j]; |
|
1197 } |
|
1198 #if 0 |
|
1199 ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP, |
|
1200 mir->cut_rhs); |
|
1201 #else |
|
1202 glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP, |
|
1203 mir->cut_rhs); |
|
1204 #endif |
|
1205 xfree(ind); |
|
1206 xfree(val); |
|
1207 return; |
|
1208 } |
|
1209 |
|
1210 static int aggregate_row(glp_tree *tree, struct MIR *mir) |
|
1211 { /* try to aggregate another row */ |
|
1212 glp_prob *mip = tree->mip; |
|
1213 int m = mir->m; |
|
1214 int n = mir->n; |
|
1215 GLPAIJ *aij; |
|
1216 IOSVEC *v; |
|
1217 int ii, j, jj, k, kk, kappa = 0, ret = 0; |
|
1218 double d1, d2, d, d_max = 0.0; |
|
1219 /* choose appropriate structural variable in the aggregated row |
|
1220 to be substituted */ |
|
1221 for (j = 1; j <= mir->agg_vec->nnz; j++) |
|
1222 { k = mir->agg_vec->ind[j]; |
|
1223 xassert(1 <= k && k <= m+n); |
|
1224 if (k <= m) continue; /* skip auxiliary var */ |
|
1225 if (mir->isint[k]) continue; /* skip integer var */ |
|
1226 if (fabs(mir->agg_vec->val[j]) < 0.001) continue; |
|
1227 /* compute distance from x[k] to its lower bound */ |
|
1228 kk = mir->vlb[k]; |
|
1229 if (kk == 0) |
|
1230 { if (mir->lb[k] == -DBL_MAX) |
|
1231 d1 = DBL_MAX; |
|
1232 else |
|
1233 d1 = mir->x[k] - mir->lb[k]; |
|
1234 } |
|
1235 else |
|
1236 { xassert(1 <= kk && kk <= m+n); |
|
1237 xassert(mir->isint[kk]); |
|
1238 xassert(mir->lb[k] != -DBL_MAX); |
|
1239 d1 = mir->x[k] - mir->lb[k] * mir->x[kk]; |
|
1240 } |
|
1241 /* compute distance from x[k] to its upper bound */ |
|
1242 kk = mir->vub[k]; |
|
1243 if (kk == 0) |
|
1244 { if (mir->vub[k] == +DBL_MAX) |
|
1245 d2 = DBL_MAX; |
|
1246 else |
|
1247 d2 = mir->ub[k] - mir->x[k]; |
|
1248 } |
|
1249 else |
|
1250 { xassert(1 <= kk && kk <= m+n); |
|
1251 xassert(mir->isint[kk]); |
|
1252 xassert(mir->ub[k] != +DBL_MAX); |
|
1253 d2 = mir->ub[k] * mir->x[kk] - mir->x[k]; |
|
1254 } |
|
1255 /* x[k] cannot be free */ |
|
1256 xassert(d1 != DBL_MAX || d2 != DBL_MAX); |
|
1257 /* d = min(d1, d2) */ |
|
1258 d = (d1 <= d2 ? d1 : d2); |
|
1259 xassert(d != DBL_MAX); |
|
1260 /* should not be close to corresponding bound */ |
|
1261 if (d < 0.001) continue; |
|
1262 if (d_max < d) d_max = d, kappa = k; |
|
1263 } |
|
1264 if (kappa == 0) |
|
1265 { /* nothing chosen */ |
|
1266 ret = 1; |
|
1267 goto done; |
|
1268 } |
|
1269 /* x[kappa] has been chosen */ |
|
1270 xassert(m+1 <= kappa && kappa <= m+n); |
|
1271 xassert(!mir->isint[kappa]); |
|
1272 /* find another row, which have not been used yet, to eliminate |
|
1273 x[kappa] from the aggregated row */ |
|
1274 for (ii = 1; ii <= m; ii++) |
|
1275 { if (mir->skip[ii]) continue; |
|
1276 for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next) |
|
1277 if (aij->col->j == kappa - m) break; |
|
1278 if (aij != NULL && fabs(aij->val) >= 0.001) break; |
|
1279 } |
|
1280 if (ii > m) |
|
1281 { /* nothing found */ |
|
1282 ret = 2; |
|
1283 goto done; |
|
1284 } |
|
1285 /* row ii has been found; include it in the aggregated list */ |
|
1286 mir->agg_cnt++; |
|
1287 xassert(mir->agg_cnt <= MAXAGGR); |
|
1288 mir->agg_row[mir->agg_cnt] = ii; |
|
1289 mir->skip[ii] = 2; |
|
1290 /* v := new row */ |
|
1291 v = ios_create_vec(m+n); |
|
1292 ios_set_vj(v, ii, 1.0); |
|
1293 for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next) |
|
1294 ios_set_vj(v, m + aij->col->j, - aij->val); |
|
1295 #if _MIR_DEBUG |
|
1296 ios_check_vec(v); |
|
1297 #endif |
|
1298 /* perform gaussian elimination to remove x[kappa] */ |
|
1299 j = mir->agg_vec->pos[kappa]; |
|
1300 xassert(j != 0); |
|
1301 jj = v->pos[kappa]; |
|
1302 xassert(jj != 0); |
|
1303 ios_linear_comb(mir->agg_vec, |
|
1304 - mir->agg_vec->val[j] / v->val[jj], v); |
|
1305 ios_delete_vec(v); |
|
1306 ios_set_vj(mir->agg_vec, kappa, 0.0); |
|
1307 #if _MIR_DEBUG |
|
1308 ios_check_vec(mir->agg_vec); |
|
1309 #endif |
|
1310 done: return ret; |
|
1311 } |
|
1312 |
|
1313 void ios_mir_gen(glp_tree *tree, void *gen) |
|
1314 { /* main routine to generate MIR cuts */ |
|
1315 glp_prob *mip = tree->mip; |
|
1316 struct MIR *mir = gen; |
|
1317 int m = mir->m; |
|
1318 int n = mir->n; |
|
1319 int i; |
|
1320 double r_best; |
|
1321 xassert(mip->m >= m); |
|
1322 xassert(mip->n == n); |
|
1323 /* obtain current point */ |
|
1324 get_current_point(tree, mir); |
|
1325 #if _MIR_DEBUG |
|
1326 /* check current point */ |
|
1327 check_current_point(mir); |
|
1328 #endif |
|
1329 /* reset bound substitution flags */ |
|
1330 memset(&mir->subst[1], '?', m+n); |
|
1331 /* try to generate a set of violated MIR cuts */ |
|
1332 for (i = 1; i <= m; i++) |
|
1333 { if (mir->skip[i]) continue; |
|
1334 /* use original i-th row as initial aggregated constraint */ |
|
1335 initial_agg_row(tree, mir, i); |
|
1336 loop: ; |
|
1337 #if _MIR_DEBUG |
|
1338 /* check aggregated row */ |
|
1339 check_agg_row(mir); |
|
1340 #endif |
|
1341 /* substitute fixed variables into aggregated constraint */ |
|
1342 subst_fixed_vars(mir); |
|
1343 #if _MIR_DEBUG |
|
1344 /* check aggregated row */ |
|
1345 check_agg_row(mir); |
|
1346 #endif |
|
1347 #if _MIR_DEBUG |
|
1348 /* check bound substitution flags */ |
|
1349 { int k; |
|
1350 for (k = 1; k <= m+n; k++) |
|
1351 xassert(mir->subst[k] == '?'); |
|
1352 } |
|
1353 #endif |
|
1354 /* apply bound substitution heuristic */ |
|
1355 bound_subst_heur(mir); |
|
1356 /* substitute bounds and build modified constraint */ |
|
1357 build_mod_row(mir); |
|
1358 #if _MIR_DEBUG |
|
1359 /* check modified row */ |
|
1360 check_mod_row(mir); |
|
1361 #endif |
|
1362 /* try to generate violated c-MIR cut for modified row */ |
|
1363 r_best = generate(mir); |
|
1364 if (r_best > 0.0) |
|
1365 { /* success */ |
|
1366 #if _MIR_DEBUG |
|
1367 /* check raw cut before back bound substitution */ |
|
1368 check_raw_cut(mir, r_best); |
|
1369 #endif |
|
1370 /* back substitution of original bounds */ |
|
1371 back_subst(mir); |
|
1372 #if _MIR_DEBUG |
|
1373 /* check the cut after back bound substitution */ |
|
1374 check_cut_row(mir, r_best); |
|
1375 #endif |
|
1376 /* final substitution to eliminate auxiliary variables */ |
|
1377 subst_aux_vars(tree, mir); |
|
1378 #if _MIR_DEBUG |
|
1379 /* check the cut after elimination of auxiliaries */ |
|
1380 check_cut_row(mir, r_best); |
|
1381 #endif |
|
1382 /* add constructed cut inequality to the cut pool */ |
|
1383 add_cut(tree, mir); |
|
1384 } |
|
1385 /* reset bound substitution flags */ |
|
1386 { int j, k; |
|
1387 for (j = 1; j <= mir->mod_vec->nnz; j++) |
|
1388 { k = mir->mod_vec->ind[j]; |
|
1389 xassert(1 <= k && k <= m+n); |
|
1390 xassert(mir->subst[k] != '?'); |
|
1391 mir->subst[k] = '?'; |
|
1392 } |
|
1393 } |
|
1394 if (r_best == 0.0) |
|
1395 { /* failure */ |
|
1396 if (mir->agg_cnt < MAXAGGR) |
|
1397 { /* try to aggregate another row */ |
|
1398 if (aggregate_row(tree, mir) == 0) goto loop; |
|
1399 } |
|
1400 } |
|
1401 /* unmark rows used in the aggregated constraint */ |
|
1402 { int k, ii; |
|
1403 for (k = 1; k <= mir->agg_cnt; k++) |
|
1404 { ii = mir->agg_row[k]; |
|
1405 xassert(1 <= ii && ii <= m); |
|
1406 xassert(mir->skip[ii] == 2); |
|
1407 mir->skip[ii] = 0; |
|
1408 } |
|
1409 } |
|
1410 } |
|
1411 return; |
|
1412 } |
|
1413 |
|
1414 /*********************************************************************** |
|
1415 * NAME |
|
1416 * |
|
1417 * ios_mir_term - terminate MIR cut generator |
|
1418 * |
|
1419 * SYNOPSIS |
|
1420 * |
|
1421 * #include "glpios.h" |
|
1422 * void ios_mir_term(void *gen); |
|
1423 * |
|
1424 * DESCRIPTION |
|
1425 * |
|
1426 * The routine ios_mir_term deletes the MIR cut generator working area |
|
1427 * freeing all the memory allocated to it. */ |
|
1428 |
|
1429 void ios_mir_term(void *gen) |
|
1430 { struct MIR *mir = gen; |
|
1431 xfree(mir->skip); |
|
1432 xfree(mir->isint); |
|
1433 xfree(mir->lb); |
|
1434 xfree(mir->vlb); |
|
1435 xfree(mir->ub); |
|
1436 xfree(mir->vub); |
|
1437 xfree(mir->x); |
|
1438 xfree(mir->agg_row); |
|
1439 ios_delete_vec(mir->agg_vec); |
|
1440 xfree(mir->subst); |
|
1441 ios_delete_vec(mir->mod_vec); |
|
1442 ios_delete_vec(mir->cut_vec); |
|
1443 xfree(mir); |
|
1444 return; |
|
1445 } |
|
1446 |
|
1447 /* eof */ |