1 /* $Id$
2 *
3 * Name: impliedBounds-exprSum.cpp
4 * Author: Pietro Belotti
5 * Purpose: implied bound enforcing for exprSum and exprGroup
6 *
7 * (C) Carnegie-Mellon University, 2006-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include "CouenneExprSum.hpp"
12 #include "CouenneExprGroup.hpp"
13 #include "CouenneConfig.h"
14 #include "CoinHelperFunctions.hpp"
15 #include "CoinFinite.hpp"
16
17 using namespace Couenne;
18
19 /// vector operation to find bound to variable in a sum
20 static CouNumber scanBounds (int, int, int *, CouNumber *, CouNumber *, int *);
21
22
23 /// implied bound processing for expression w = x+y+t+..., upon change
24 /// in lower- and/or upper bound of w, whose index is wind
25
impliedBound(int wind,CouNumber * l,CouNumber * u,t_chg_bounds * chg,enum auxSign sign)26 bool exprSum::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg, enum auxSign sign) {
27
28 /**
29 * An expression
30 *
31 * \f$w = a0 + \sum_{i\in I1} a_i x_i + \sum_{i\in I2} a_i x_i\f$
32 *
33 * is given such that all $a_i$ are positive for $i \in I1$ and
34 * negative for $i in I2$. If the bounds on $w \in [l,b]$, implied
35 * bounds on all $x_i, i\in I1 \cup I2$ are as follows:
36 *
37 * \f$\forall i\in I1\f$
38 *
39 * \f$x_i \ge (l - a0 - sum_{j\in I1 | j\neq i} a_j u_j - sum_{j\in I2} a_j l_j) / a_i\f$
40 * \f$x_i \le (u - a0 - sum_{j\in I1 | j\neq i} a_j l_j - sum_{j\in I2} a_j u_j) / a_i\f$
41 *
42 * \f$\forall i\in I2\f$
43 *
44 * \f$x_i \ge (u - a0 - sum_{j\in I1} a_j u_j - sum_{j\in I2 | j\neq i} a_j l_j) / a_i\f$
45 * \f$x_i \le (l - a0 - sum_{j\in I1} a_j l_j - sum_{j\in I2 | j\neq i} a_j u_j) / a_i\f$,
46 *
47 * where \f$l_i\f$ and \f$u_i\f$ are lower and upper bound,
48 * respectively, of \f$x_i\f$. We also have to check if some of
49 * these bounds are infinite.
50 */
51
52 // compute number of pos/neg coefficients and sum up constants
53
54 CouNumber
55 a0 = 0., // constant term in the sum
56 wl = sign == expression::AUX_GEQ ? -COIN_DBL_MAX : l [wind],
57 wu = sign == expression::AUX_LEQ ? COIN_DBL_MAX : u [wind];
58
59 // quick check: if both are infinite, nothing is implied...
60
61 if ((wl < -COUENNE_INFINITY) &&
62 (wu > COUENNE_INFINITY))
63 return false;
64
65 int
66 nterms = nargs_, // # nonconstant terms
67 nlin = 0; // # linear terms
68
69 exprGroup *eg = NULL;
70
71 bool isExprGroup = (code () == COU_EXPRGROUP);
72
73 if (isExprGroup) { // if exprGroup, compute no. linear terms
74
75 eg = dynamic_cast <exprGroup *> (this);
76
77 a0 += eg -> getc0 ();
78 nlin += eg -> lcoeff (). size ();
79 }
80
81 nterms += nlin;
82
83 // Coefficients and indices of the positive and the negative
84 // non-constant parts of the sum (at most nlin are negative, as all
85 // "nonlinear" terms have coefficient 1)
86
87 CouNumber *C1 = (CouNumber *) malloc (nterms * sizeof (CouNumber)),
88 *C2 = (CouNumber *) malloc (nlin * sizeof (CouNumber));
89 int *I1 = (int *) malloc (nterms * sizeof (int)),
90 *I2 = (int *) malloc (nlin * sizeof (int));
91
92 int ipos, ineg = ipos = 0; // #pos and #neg terms
93
94 std::set <int> intSet; // contains indices of integer variables
95
96 // classify terms as positive or constant for the exprSum
97
98 for (int i = nargs_; i--;) {
99
100 int index = arglist_ [i] -> Index ();
101
102 if (index < 0) // assume a non variable is a constant
103 a0 += arglist_ [i] -> Value ();
104 else {
105 I1 [ipos] = index;
106 C1 [ipos++] = 1.;
107
108 // add entry to integer variable
109 if (arglist_ [i] -> isDefinedInteger ())
110 intSet.insert (index);
111 }
112 }
113
114 // classify terms as pos/neg or constant for the remaining of exprGroup
115
116 if (isExprGroup) {
117
118 exprGroup::lincoeff &lcoe = eg -> lcoeff ();
119
120 for (register exprGroup::lincoeff::iterator el = lcoe.begin ();
121 el != lcoe.end (); ++el) {
122
123 register CouNumber coe = el -> second;
124 register int ind = el -> first -> Index ();
125
126 if (coe > 0.) {I1 [ipos] = ind; C1 [ipos++] = coe;}
127 else if (coe < -0.) {I2 [ineg] = ind; C2 [ineg++] = coe;}
128
129 if (el -> first -> isDefinedInteger ())
130 intSet.insert (ind);
131 }
132 }
133
134 // realloc to save some memory
135
136 /*printf (" a0 = %g, w in [%g, %g]\n", a0, wl, wu);
137 printf (" pos: "); for (int i=0; i<ipos; i++) printf ("%g x%d [%g,%g] ", C1 [i], I1 [i], l [I1 [i]], u [I1 [i]]); printf ("\n");
138 printf (" neg: "); for (int i=0; i<ineg; i++) printf ("%g x%d [%g,%g] ", C2 [i], I2 [i], l [I2 [i]], u [I2 [i]]); printf ("\n");*/
139
140 // now we have correct I1, I2, C1, C2, ipos, ineg, and a0
141
142 // indices of the variable in I1 or I2 with infinite lower or upper
143 // bound. If more than one is found, it is set to -2
144
145 int infLo1 = -1, infLo2 = -1,
146 infUp1 = -1, infUp2 = -1;
147
148 // Upper bound of the sum, considering lower/upper bounds of the
149 // variables but negliging the infinite ones:
150 //
151 // lower = $a0 + \sum_{i in I_1: a_i < \infinity} a_i l_i
152 // + \sum_{i in I_2: a_i > -\infinity} a_i u_i$
153 //
154 // upper = $a0 + \sum_{i in I_1: a_i < \infinity} a_i u_i
155 // + \sum_{i in I_2: a_i > -\infinity} a_i l_i$
156
157 CouNumber
158
159 lower = a0 + scanBounds (ipos, -1, I1, C1, l, &infLo1)
160 + scanBounds (ineg, +1, I2, C2, u, &infUp2),
161
162 upper = a0 + scanBounds (ipos, +1, I1, C1, u, &infUp1)
163 + scanBounds (ineg, -1, I2, C2, l, &infLo2);
164
165 // Compute lower bound for all or for some of the variables: There
166 // is a bound for all variables only if both infUp1 and infLo2 are
167 // -1, otherwise there is a bound only for one variable if one is -1
168 // and the other is nonnegative.
169
170 bool tighter = false;
171
172 // check if there is room for improvement
173
174 {
175 CouNumber
176 slackL = lower - wl, // each is positive if no implication
177 slackU = wu - upper; //
178
179 // if lower < wl or upper > wu, some bounds can be tightened.
180 //
181 // otherwise, there is no implication, but if lower
182 //
183 // steal some work to bound propagation...
184
185 if ((slackL > 0.) &&
186 (infLo1 == -1) &&
187 (infUp2 == -1)) { // no implication on lower
188
189 // propagate lower bound to w
190 if ((sign != expression::AUX_LEQ) && (updateBound (-1, l + wind, lower))) {
191 chg [wind].setLower(t_chg_bounds::CHANGED);
192 tighter = true;
193 if (intSet.find (wind)!= intSet.end ())
194 l [wind] = ceil (l [wind] - COUENNE_EPS);
195 }
196
197 if ((slackU > 0.) &&
198 (infLo2 == -1) &&
199 (infUp1 == -1)) { // no implication on upper
200
201 // propagate upper bound to w
202 if ((sign != expression::AUX_GEQ) && (updateBound (+1, u + wind, upper))) {
203 chg [wind].setUpper(t_chg_bounds::CHANGED);
204 tighter = true;
205 if (intSet.find (wind)!= intSet.end ())
206 u [wind] = floor (u [wind] + COUENNE_EPS);
207 }
208
209 free (I1); free (I2);
210 free (C1); free (C2);
211
212 return false; // both bounds were weak, no implication possible
213 }
214 }
215 else if ((slackU > 0.) &&
216 (infLo2 == -1) &&
217 (infUp1 == -1)) {
218
219 // propagate upper bound to w
220 if ((sign != expression::AUX_GEQ) && (updateBound (+1, u + wind, upper))) {
221 tighter = true;
222 chg [wind].setUpper(t_chg_bounds::CHANGED);
223 if (intSet.find (wind)!= intSet.end ())
224 u [wind] = floor (u [wind] + COUENNE_EPS);
225 }
226 }
227 }
228
229 // Subtle... make two copies of lower and upper to avoid updating
230 // bounds with some previously updated bounds
231
232 // first of all, find maximum index in I1 and I2
233
234 int maxind = -1;
235
236 for (register int i=ipos; i--; I1++) if (*I1 > maxind) maxind = *I1;
237 for (register int i=ineg; i--; I2++) if (*I2 > maxind) maxind = *I2;
238
239 I1 -= ipos;
240 I2 -= ineg;
241
242 CouNumber *lc = (CouNumber *) malloc (++maxind * sizeof (CouNumber));
243 CouNumber *uc = (CouNumber *) malloc (maxind * sizeof (CouNumber));
244
245 CoinCopyN (l, maxind, lc);
246 CoinCopyN (u, maxind, uc);
247
248 // Update lowers in I1 and uppers in I2
249
250 if ((infLo1 == -1) && (infUp2 == -1) && (wu < COUENNE_INFINITY / 1e10)) {
251 // All finite bounds. All var. bounds can be tightened.
252
253 // tighten upper bound of variables in I1
254 for (register int i=ipos; i--;) {
255 int ind = I1 [i];
256 if ((tighter = (updateBound (+1, u + ind, (wu - lower) / C1 [i] + lc [ind]) || tighter))) {
257 chg [ind].setUpper(t_chg_bounds::CHANGED);
258 if (intSet.find (ind)!= intSet.end ())
259 u [ind] = floor (u [ind] + COUENNE_EPS);
260 }
261 }
262
263 // tighten lower bound of variables in I2
264 for (register int i=ineg; i--;) {
265 int ind = I2 [i];
266 if ((tighter = (updateBound (-1, l + ind, (wu - lower) / C2 [i] + uc [ind]) || tighter))) {
267 chg [ind].setLower(t_chg_bounds::CHANGED);
268 if (intSet.find (ind)!= intSet.end ())
269 l [ind] = ceil (l [ind] - COUENNE_EPS);
270 }
271 }
272 } else
273
274 if ((infLo1 >= 0) && (infUp2 == -1)) { // There is one infinite lower bound in I1
275 int ind = I1 [infLo1];
276 if ((tighter = (updateBound (+1, u + ind, (wu - lower) / C1 [infLo1]) || tighter))) {
277 chg [ind].setUpper(t_chg_bounds::CHANGED);
278 if (intSet.find (ind)!= intSet.end ())
279 u [ind] = floor (u [ind] + COUENNE_EPS);
280 }
281 }
282 else
283 if ((infLo1 == -1) && (infUp2 >= 0)) { // There is one infinite upper bound in I2
284 int ind = I2 [infUp2];
285 if ((tighter = (updateBound (-1, l + ind, (wu - lower) / C2 [infUp2]) || tighter))) {
286 chg [ind].setLower(t_chg_bounds::CHANGED);
287 if (intSet.find (ind)!= intSet.end ())
288 l [ind] = ceil (l [ind] - COUENNE_EPS);
289 }
290 }
291
292 // Update uppers in I1 and lowers in I2
293
294 if ((infUp1 == -1) && (infLo2 == -1) && (wl > -COUENNE_INFINITY / 1e10)) {
295 // All finite bounds. All var. bounds can be tightened.
296
297 for (register int i=ipos; i--;) {
298 int ind = I1 [i];
299 if ((tighter = (updateBound (-1, l + ind, (wl - upper) / C1 [i] + uc [ind]) || tighter))) {
300 chg [ind].setLower(t_chg_bounds::CHANGED);
301 if (intSet.find (ind) != intSet.end ())
302 l [ind] = ceil (l [ind] - COUENNE_EPS);
303 }
304 }
305
306 for (register int i=ineg; i--;) {
307 int ind = I2 [i];
308 if ((tighter = (updateBound (+1, u + ind, (wl - upper) / C2 [i] + lc [ind]) || tighter))) {
309 chg [ind].setUpper(t_chg_bounds::CHANGED);
310 if (intSet.find (ind) != intSet.end ())
311 u [ind] = floor (u [ind] + COUENNE_EPS);
312 }
313 }
314 } else
315
316 if ((infUp1 >= 0) && (infLo2 == -1)) { // There is one infinite lower bound in I2
317 int ind = I1 [infUp1];
318 if ((tighter = (updateBound (-1, l + ind, (wl - upper) / C1 [infUp1]) || tighter))) {
319 chg [ind].setLower(t_chg_bounds::CHANGED);
320 if (intSet.find (ind) != intSet.end ())
321 l [ind] = ceil (l [ind] - COUENNE_EPS);
322 }
323 }
324 else
325 if ((infUp1 == -1) && (infLo2 >= 0)) { // There is one infinite upper bound in I1
326 int ind = I2 [infLo2];
327 if ((tighter = (updateBound (+1, u + ind, (wl - upper) / C2 [infLo2]) || tighter))) {
328 chg [ind].setUpper(t_chg_bounds::CHANGED);
329 if (intSet.find (ind) != intSet.end ())
330 u [ind] = floor (u [ind] + COUENNE_EPS);
331 }
332 }
333
334 // ...phew!
335
336 free (I1); free (I2);
337 free (C1); free (C2);
338 free (lc); free (uc);
339
340 return tighter;
341 }
342
343
344 /// sum bounds depending on coefficients
345
scanBounds(int num,int sign,int * indices,CouNumber * coeff,CouNumber * bounds,int * infnum)346 static CouNumber scanBounds (int num, /// cardinality of the set (I1 or I2)
347 int sign, /// +1: check for +inf, -1: check for -inf
348 int *indices, /// indices in the set, $i\in I_k$
349 CouNumber *coeff, /// coefficients in the set
350 CouNumber *bounds, /// variable bounds to check
351 int *infnum) { /// return value: index of variable with infinite
352 /// bound, or -2 if more than one exist
353 CouNumber bound = 0.;
354
355 for (register int i = num; i--;) {
356
357 CouNumber bd = bounds [indices [i]];
358
359 // be sensitive here, check for bounds a little within the finite realm
360
361 if (((sign > 0) ? bd : -bd) > COUENNE_INFINITY / 1e10) {
362
363 bounds [indices [i]] = (sign > 0) ? COIN_DBL_MAX : -COIN_DBL_MAX;
364 //bounds [indices [i]] = (sign > 0) ? 1e300 : -1e300;
365
366 // this variable has an infinite bound, mark it
367 if (*infnum == -1) *infnum = i; // first variable with infinite bound, so far
368 else if (*infnum >= 0) *infnum = -2; // oops... more than one found, no finite bound
369 }
370 else bound += coeff [i] * bd; // no infinity detected, sum a weighted, finite bound
371 }
372
373 return bound;
374 }
375