source: trunk/Couenne/src/bound_tightening/operators/impliedBounds-exprSum.cpp @ 115

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

fix right value of objective computed in checkNLP -- numerics can make a bound interval infeasible. Some cleanup

File size: 11.8 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-09.
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 (register 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  CoinCopyN (l, maxind, lc);
241  CoinCopyN (u, maxind, uc);
242
243  /*for (register int i = maxind; i--;) {
244    *lc++ = *l++;
245    *uc++ = *u++;
246    //lc [i] = l [i];
247    //uc [i] = u [i];
248  }
249
250  lc -= maxind; uc -= maxind;
251  l  -= maxind; u  -= maxind;*/
252
253  // Update lowers in I1 and uppers in I2
254
255  if ((infLo1 == -1) && (infUp2 == -1) && (wu < COUENNE_INFINITY / 1e10)) { 
256    // All finite bounds. All var. bounds can be tightened.
257
258    // tighten upper bound of variables in I1
259    for (register int i=ipos; i--;) {
260      int ind = I1 [i];
261      if (tighter = updateBound (+1, u + ind, (wu - lower) / C1 [i] + lc [ind]) || tighter) {
262        chg [ind].setUpper(t_chg_bounds::CHANGED);
263        if (intSet.find (ind)!= intSet.end ()) 
264          u [ind] = floor (u [ind] + COUENNE_EPS);
265      }
266    }
267
268    // tighten lower bound of variables in I2
269    for (register int i=ineg; i--;) {
270      int ind = I2 [i];
271      if (tighter = updateBound (-1, l + ind, (wu - lower) / C2 [i] + uc [ind]) || tighter) {
272        chg [ind].setLower(t_chg_bounds::CHANGED);
273        if (intSet.find (ind)!= intSet.end ()) 
274          l [ind] = ceil (l [ind] - COUENNE_EPS);
275      }
276    }
277  } else
278
279    if ((infLo1 >= 0) && (infUp2 == -1)) {    // There is one infinite lower bound in I1
280      int ind = I1 [infLo1];
281      if (tighter = updateBound (+1, u + ind, (wu - lower) / C1 [infLo1]) || tighter) {
282        chg [ind].setUpper(t_chg_bounds::CHANGED);
283        if (intSet.find (ind)!= intSet.end ()) 
284          u [ind] = floor (u [ind] + COUENNE_EPS);
285      }
286    }
287    else 
288      if ((infLo1 == -1) && (infUp2 >= 0)) {  // There is one infinite upper bound in I2
289        int ind = I2 [infUp2];
290        if (tighter = updateBound (-1, l + ind, (wu - lower) / C2 [infUp2]) || tighter) {
291          chg [ind].setLower(t_chg_bounds::CHANGED);
292          if (intSet.find (ind)!= intSet.end ()) 
293            l [ind] = ceil (l [ind] - COUENNE_EPS);
294        }
295      }
296
297  // Update uppers in I1 and lowers in I2
298
299  if ((infUp1 == -1) && (infLo2 == -1) && (wl > -COUENNE_INFINITY / 1e10)) { 
300    // All finite bounds. All var. bounds can be tightened.
301
302    for (register int i=ipos; i--;) {
303      int ind = I1 [i];
304      if (tighter = updateBound (-1, l + ind, (wl - upper) / C1 [i] + uc [ind]) || tighter) {
305        chg [ind].setLower(t_chg_bounds::CHANGED); 
306        if (intSet.find (ind) != intSet.end ()) 
307          l [ind] = ceil (l [ind] - COUENNE_EPS);
308      }
309    }
310
311    for (register int i=ineg; i--;) {
312      int ind = I2 [i];
313      if (tighter = updateBound (+1, u + ind, (wl - upper) / C2 [i] + lc [ind]) || tighter) {
314        chg [ind].setUpper(t_chg_bounds::CHANGED);
315        if (intSet.find (ind) != intSet.end ()) 
316          u [ind] = floor (u [ind] + COUENNE_EPS);
317      }
318    }
319  } else 
320
321    if ((infUp1 >= 0) && (infLo2 == -1)) { // There is one infinite lower bound in I2
322      int ind = I1 [infUp1];
323      if (tighter = updateBound (-1, l + ind, (wl - upper) / C1 [infUp1]) || tighter) {
324        chg [ind].setLower(t_chg_bounds::CHANGED);
325        if (intSet.find (ind) != intSet.end ()) 
326          l [ind] = ceil (l [ind] - COUENNE_EPS);
327      }
328    }
329    else 
330      if ((infUp1 == -1) && (infLo2 >= 0)) {  // There is one infinite upper bound in I1
331        int ind = I2 [infLo2];
332        if (tighter = updateBound (+1, u + ind, (wl - upper) / C2 [infLo2]) || tighter) {
333          chg [ind].setUpper(t_chg_bounds::CHANGED);
334          if (intSet.find (ind) != intSet.end ()) 
335            u [ind] = floor (u [ind] + COUENNE_EPS);
336        }
337      }
338
339  // ...phew!
340
341  free (I1); free (I2);
342  free (C1); free (C2);
343  free (lc); free (uc);
344
345  return tighter;
346}
347
348
349/// sum bounds depending on coefficients
350
351static CouNumber scanBounds (int        num,      /// cardinality of the set (I1 or I2)
352                             int        sign,     /// +1: check for +inf, -1: check for -inf
353                             int       *indices,  /// indices in the set, $i\in I_k$
354                             CouNumber *coeff,    /// coefficients in the set
355                             CouNumber *bounds,   /// variable bounds to check
356                             int       *infnum) { /// return value: index of variable with infinite
357                                                  /// bound, or -2 if more than one exist
358  CouNumber bound = 0.;
359
360  for (register int i = num; i--;) {
361
362    CouNumber bd = bounds [indices [i]];
363
364    // be sensitive here, check for bounds a little within the finite realm
365
366    if (((sign > 0) ? bd : -bd) > COUENNE_INFINITY / 1e10) {
367
368      bounds [indices [i]] = (sign > 0) ? COIN_DBL_MAX : -COIN_DBL_MAX;
369      //bounds [indices [i]] = (sign > 0) ? 1e300 : -1e300;
370
371      // this variable has an infinite bound, mark it
372      if      (*infnum == -1) *infnum =  i; // first variable with infinite bound, so far
373      else if (*infnum >=  0) *infnum = -2; // oops... more than one found, no finite bound
374    }
375    else bound += coeff [i] * bd; // no infinity detected, sum a weighted, finite bound
376  }
377
378  return bound;
379}
Note: See TracBrowser for help on using the repository browser.