source: stable/0.1/Couenne/src/bound_tightening/operators/impliedBounds-exprSum.cpp @ 113

Last change on this file since 113 was 113, checked in by pbelotti, 13 years ago

skip restore unused if there are none (see bug pointed out by C. D'Ambrosio). Restored precedence in if+assignment+OR in impliedBounds-exprSum while (hopefully) avoiding warnings Stefan mentioned.

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