impliedBounds-exprMul.cpp
Go to the documentation of this file.
1 /* $Id: impliedBounds-exprMul.cpp 746 2011-07-22 21:56:44Z pbelotti $
2  *
3  * Name: impliedBounds-exprMul.cpp
4  * Author: Pietro Belotti
5  * Purpose: implied bounds for multiplications
6  *
7  * (C) Carnegie-Mellon University, 2006-10.
8  * This file is licensed under the Eclipse Public License (EPL)
9  */
10 
11 #include <stdio.h>
12 #include <stdlib.h>
13 
14 #include "CouenneExprMul.hpp"
15 #include "CouennePrecisions.hpp"
16 #include "CouenneConfig.h"
17 #include "CoinFinite.hpp"
18 
19 //#define FM_CHECK
20 //#define FM_MOD
21 
22 using namespace Couenne;
23 
24 
27 
28 bool exprMul::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg, enum auxSign sign) {
29 
30  bool resL, resU = resL = false;
31  int ind;
32 
33  CouNumber
34  wl = sign == expression::AUX_GEQ ? -COIN_DBL_MAX : l [wind],
35  wu = sign == expression::AUX_LEQ ? COIN_DBL_MAX : u [wind];
36 
37  if ((arglist_ [ind=0] -> Type () <= CONST) ||
38  (arglist_ [ind=1] -> Type () <= CONST)) {
39 
40  // at least one constant in product w=cx:
41  //
42  // wl/c <= x <= wu/c, if c is positive
43  // wu/c <= x <= wl/c, if c is negative
44 
45  CouNumber c = arglist_ [ind] -> Value ();
46 
47  bool argInt = arglist_ [1-ind] -> isInteger ();
48 
49  // get the index of the nonconstant part
50  ind = arglist_ [1-ind] -> Index ();
51 
52  if (ind==-1) // should not happen, it is a product of constants
53  return false;
54 
55  if (c > COUENNE_EPS) {
56 
57  resL = (wl > - COUENNE_INFINITY) && updateBound (-1, l + ind, argInt ? ceil (wl / c - COUENNE_EPS) : (wl / c));
58  resU = (wu < COUENNE_INFINITY) && updateBound ( 1, u + ind, argInt ? floor (wu / c + COUENNE_EPS) : (wu / c));
59 
60  } else if (c < - COUENNE_EPS) {
61 
62  resL = (wu < COUENNE_INFINITY) && updateBound (-1, l + ind, argInt ? ceil (wu / c - COUENNE_EPS) : (wu / c));
63  resU = (wl > - COUENNE_INFINITY) && updateBound ( 1, u + ind, argInt ? floor (wl / c + COUENNE_EPS) : (wl / c));
64  }
65 
66  if (resL) chg [ind].setLower(t_chg_bounds::CHANGED);
67  if (resU) chg [ind].setUpper(t_chg_bounds::CHANGED);
68 
69  /*printf ("w_%d [%g,%g] -------> x_%d in [%g,%g] ",
70  wind, l [wind], u [wind],
71  ind, l [ind], u [ind]);*/
72  } else {
73 
74  // these bounds would be implied by McCormick's convexification,
75  // however we write them explicitly for internal use within bound
76  // tightening, as otherwise they would be known to Clp only.
77 
78  int xi = arglist_ [0] -> Index (),
79  yi = arglist_ [1] -> Index ();
80 
81  CouNumber
82  *xl = l + xi, *yl = l + yi,
83  *xu = u + xi, *yu = u + yi;
84 
85  // w's lower bound ///////////////////////////////////////////
86 
87  bool resxL, resxU, resyL, resyU =
88  resxL = resxU = resyL = false;
89 
90  bool xlIsZero = (fabs(*xl) < COUENNE_EPS);
91  bool xuIsZero = (fabs(*xu) < COUENNE_EPS);
92  bool ylIsZero = (fabs(*yl) < COUENNE_EPS);
93  bool yuIsZero = (fabs(*yu) < COUENNE_EPS);
94  bool wlIsZero = (fabs(wl) < COUENNE_EPS);
95  bool wuIsZero = (fabs(wu) < COUENNE_EPS);
96 
97 #ifdef FM_MOD
98 
99  if(wlIsZero) {
100  if ((!xuIsZero) && (!yuIsZero) && (*xu * *yu < wl)) {
101  if (!ylIsZero) {
102  resxU = (*xu * *yl < wl) && updateBound (+1, xu, wl / *yl);
103  xuIsZero = (fabs(*xu) < COUENNE_EPS ? true : false);
104  }
105  if (!xlIsZero) {
106  resyU = (*xl * *yu < wl) && updateBound (+1, yu, wl / *xl);
107  yuIsZero = (fabs(*yu) < COUENNE_EPS ? true : false);
108  }
109  }
110 
111  // point C in central infeasible area
112 
113  if ((!xlIsZero) && (!ylIsZero) && (*xl * *yl < wl)) {
114  if (!yuIsZero) {
115  resxL = (*xl * *yu < wl) && updateBound (-1, xl, wl / *yu);
116  xlIsZero = (fabs(*xl) < COUENNE_EPS ? true : false);
117  }
118  if (!xuIsZero) {
119  resyL = (*xu * *yl < wl) && updateBound (-1, yl, wl / *xu);
120  ylIsZero = (fabs(*yl) < COUENNE_EPS ? true : false);
121  }
122  }
123  }
124  else // wl is not zero
125 
126 #endif
127  {
128  if (wl >= 0.) {
129  // point B in central infeasible area
130  if (*xu * *yu < wl) {
131 
132 #ifdef FM_CHECK
133  if((*xu * *yl < wl) && (ylIsZero)) {
134  printf("ylIsZero (A): %g\n", *yl);
135  exit(1);
136  }
137  if((*xl * *yu < wl) && (xlIsZero)) {
138  printf("xlIsZero (B): %g\n", *xl);
139  exit(1);
140  }
141 #endif
142 
143  resxU = (*xu * *yl < wl) && updateBound (+1, xu, wl / *yl);
144  xuIsZero = (fabs(*xu) < COUENNE_EPS ? true : false);
145  resyU = (*xl * *yu < wl) && updateBound (+1, yu, wl / *xl);
146  yuIsZero = (fabs(*yu) < COUENNE_EPS ? true : false);
147  }
148 
149  // point C in central infeasible area
150 
151  if (*xl * *yl < wl) {
152 
153 #ifdef FM_CHECK
154  if((*xl * *yu < wl) && (yuIsZero)) {
155  printf("yuIsZero (C): %g\n", *yu);
156  exit(1);
157  }
158  if((*xu * *yl < wl) && (xuIsZero)) {
159  printf("xuIsZero (D): %g\n", *xu);
160  exit(1);
161  }
162 #endif
163  resxL = (*xl * *yu < wl) && (!yuIsZero) && updateBound (-1, xl, wl / *yu);
164  xlIsZero = (fabs(*xl) < COUENNE_EPS ? true : false);
165  resyL = (*xu * *yl < wl) && (!xuIsZero) && updateBound (-1, yl, wl / *xu);
166  ylIsZero = (fabs(*yl) < COUENNE_EPS ? true : false);
167  }
168 
169  } else if (wl > -COUENNE_INFINITY) {
170 
171  // the infeasible set is a hyperbola with two branches
172 
173  // upper left
174 
175 #ifdef FM_CHECK
176  if(((*xl * *yl < wl) && (*yl > 0)) && (ylIsZero)) {
177  printf("ylIsZero (E): %g\n", *yl);
178  }
179  if(((*xu * *yu < wl) && (*yu > 0)) && (xuIsZero)) {
180  printf("xuIsZero (F): %g\n", *xu);
181  exit(1);
182  }
183  if(((*xl * *yl < wl) && (*yl < 0)) && (xlIsZero)) {
184  printf("xlIsZero (G): %g\n", *xl);
185  exit(1);
186  }
187  if(( (*xu * *yu < wl) && (*yu < 0)) && (yuIsZero)) {
188  printf("yuIsZero (H): %g\n", *yu);
189  }
190 #endif
191 
192  resxL = (*xl * *yl < wl) && (*yl > 0.) && (!ylIsZero) && updateBound (-1, xl, wl / *yl); // point C
193  xlIsZero = (fabs(*xl) < COUENNE_EPS ? true : false);
194  resyU = (*xu * *yu < wl) && (*yu > 0.) && (!xuIsZero) && updateBound (+1, yu, wl / *xu); // point B
195  yuIsZero = (fabs(*yu) < COUENNE_EPS ? true : false);
196 
197  // lower right
198  resyL = (*xl * *yl < wl) && (*yl < 0.) && (!xlIsZero) && updateBound (-1, yl, wl / *xl); // point C
199  ylIsZero = (fabs(*yl) < COUENNE_EPS ? true : false);
200  resxU = (*xu * *yu < wl) && (*yu < 0.) && (!yuIsZero) && updateBound (+1, xu, wl / *yu); // point B
201  xuIsZero = (fabs(*xu) < COUENNE_EPS ? true : false);
202  }
203  }
204 
205  bool
206  xInt = arglist_ [0] -> isInteger (),
207  yInt = arglist_ [1] -> isInteger ();
208 
209  if (resxL && xInt) *xl = ceil (*xl - COUENNE_EPS);
210  if (resxU && xInt) *xu = floor (*xu + COUENNE_EPS);
211  if (resyL && yInt) *yl = ceil (*yl - COUENNE_EPS);
212  if (resyU && yInt) *yu = floor (*yu + COUENNE_EPS);
213 
214  // w's upper bound ///////////////////////////////////////////
215 #ifdef FM_MOD
216 
217  if(wuIsZero) {
218 
219  if((!xuIsZero) && (!ylIsZero) && (*xu * *yl > wu)) {
220  if(!yuIsZero) {
221  resxU = (*xu * *yu > wu) && updateBound (+1, xu, wu / *yu) || resxU;
222  xuIsZero = (fabs(*xu) < COUENNE_EPS ? true : false);
223  }
224  if(!xlIsZero) {
225  resyL = (*xl * *yl > wu) && updateBound (-1, yl, wu / *xl) || resyL;
226  ylIsZero = (fabs(*yl) < COUENNE_EPS ? true : false);
227  }
228  }
229 
230  if((!xlIsZero) && (!yuIsZero) && (*xl * *yu > wu)) {
231  if(!ylIsZero) {
232  resxL = (*xl * *yl > wu) && updateBound (-1, xl, wu / *yl) || resxL;
233  xlIsZero = (fabs(*xl) < COUENNE_EPS ? true : false);
234  }
235  if(!xuIsZero) {
236  resyU = (*xu * *yu > wu) && updateBound (+1, yu, wu / *xu) || resyU;
237  yuIsZero = (fabs(*yu) < COUENNE_EPS ? true : false);
238  }
239  }
240  }
241 
242  else // wu is not zero
243 
244 #endif /* FM_MOD */
245 
246  {
247  if (wu >= 0.) {
248  if (wu < COUENNE_INFINITY) {
249  // the infeasible set is a hyperbola with two branches
250 
251 #ifdef FM_CHECK
252  if(((*xu * *yl > wu) && (*yl > 0)) && (ylIsZero)) {
253  printf("ylIsZero (A2): yl: %g xu: %g wu: %g\n", *yl, *xu, wu);
254  exit(1);
255  }
256  if(((*xl * *yu > wu) && (*yu > 0)) && (xlIsZero)) {
257  printf("xlIsZero (B2): %g\n", *xl);
258  exit(1);
259  }
260 #endif
261 
262  // upper right
263  resxU = ((*xu * *yl > wu) && (*yl > 0.) && (!ylIsZero) && updateBound (+1, xu, wu / *yl)) || resxU; // point D
264  xuIsZero = (fabs(*xu) < COUENNE_EPS ? true : false);
265  resyU = ((*xl * *yu > wu) && (*yu > 0.) && (!xlIsZero) && updateBound (+1, yu, wu / *xl)) || resyU; // point A
266  yuIsZero = (fabs(*yu) < COUENNE_EPS ? true : false);
267 
268 #ifdef FM_CHECK
269  if(((*xl * *yu > wu) && (*yu < 0)) && (yuIsZero)) {
270  printf("yuIsZero (C2): %g\n", *yu);
271  }
272  if(((*xu * *yl > wu) && (*yl < 0)) && (xuIsZero)) {
273  printf("xuIsZero (D2): %g\n", *xu);
274  exit(1);
275  }
276 #endif
277 
278  // lower left
279  resxL = ((*xl * *yu > wu) && (*yu < 0.) && (!yuIsZero) && updateBound (-1, xl, wu / *yu)) || resxL; // point A
280  xlIsZero = (fabs(*xl) < COUENNE_EPS ? true : false);
281  resyL = ((*xu * *yl > wu) && (*yl < 0.) && (!xuIsZero) && updateBound (-1, yl, wu / *xu)) || resyL; // point D
282  ylIsZero = (fabs(*yl) < COUENNE_EPS ? true : false);
283  }
284 
285  } else {
286 
287  // point D in central infeasible area
288 
289  if (*xu * *yl > wu) {
290 
291 #ifdef FM_CHECK
292  if((*xu * *yu > wu) && (yuIsZero)) {
293  printf("yuIsZero (E2): %g\n", *yu);
294  exit(1);
295  }
296  if((*xl * *yl > wu) && (xlIsZero)) {
297  printf("xlIsZero (F2): %g\n", *xl);
298  exit(1);
299  }
300 #endif
301 
302  resxU = ((*xu * *yu > wu) && (!yuIsZero) && updateBound (+1, xu, wu / *yu)) || resxU;
303  xuIsZero = (fabs(*xu) < COUENNE_EPS ? true : false);
304  resyL = ((*xl * *yl > wu) && (!xlIsZero) && updateBound (-1, yl, wu / *xl)) || resyL;
305  ylIsZero = (fabs(*yl) < COUENNE_EPS ? true : false);
306  }
307 
308  // point A in central infeasible area
309 
310  if (*xl * *yu > wu) {
311 
312 #ifdef FM_CHECK
313  if((*xl * *yl > wu) && (ylIsZero)) {
314  printf("ylIsZero (G2): %g\n", *yl);
315  exit(1);
316  }
317  if((*xu * *yu > wu) && (xuIsZero)) {
318  printf("xuIsZero (H2): %g\n", *xu);
319  exit(1);
320  }
321 #endif
322 
323  resxL = ((*xl * *yl > wu) && (!ylIsZero) && updateBound (-1, xl, wu / *yl)) || resxL;
324  xlIsZero = (fabs(*xl) < COUENNE_EPS ? true : false);
325  resyU = ((*xu * *yu > wu) && (!xuIsZero) && updateBound (+1, yu, wu / *xu)) || resyU;
326  yuIsZero = (fabs(*yu) < COUENNE_EPS ? true : false);
327  }
328  }
329  }
330 
331  // extra integrality check
332 
333  if (resxL) {chg [xi].setLower(t_chg_bounds::CHANGED); if (xInt) *xl = ceil (*xl - COUENNE_EPS);}
334  if (resxU) {chg [xi].setUpper(t_chg_bounds::CHANGED); if (xInt) *xu = floor (*xu + COUENNE_EPS);}
335  if (resyL) {chg [yi].setLower(t_chg_bounds::CHANGED); if (yInt) *yl = ceil (*yl - COUENNE_EPS);}
336  if (resyU) {chg [yi].setUpper(t_chg_bounds::CHANGED); if (yInt) *yu = floor (*yu + COUENNE_EPS);}
337 
338  resL = resxL || resyL;
339  resU = resxU || resyU;
340  }
341 
342  return (resL || resU);
343 }
Couenne::updateBound
bool updateBound(register int sign, register CouNumber *dst, register CouNumber src)
updates maximum violation.
Definition: CouenneExpression.hpp:279
CouenneExprMul.hpp
Couenne::t_chg_bounds
status of lower/upper bound of a variable, to be checked/modified in bound tightening ...
Definition: CouenneTypes.hpp:66
Couenne::exprOp::isInteger
virtual bool isInteger()
is this expression integer?
Definition: exprOp.cpp:188
Couenne::t_chg_bounds::setLower
void setLower(ChangeStatus lower)
Definition: CouenneTypes.hpp:85
Couenne::t_chg_bounds::setUpper
void setUpper(ChangeStatus upper)
Definition: CouenneTypes.hpp:86
Couenne::exprOp::Type
virtual enum nodeType Type() const
Node type.
Definition: CouenneExprOp.hpp:41
Couenne::expression::Index
virtual int Index() const
Return index of variable (only valid for exprVar and exprAux)
Definition: CouenneExpression.hpp:73
Couenne::expression::auxSign
auxSign
&quot;sign&quot; of the constraint defining an auxiliary.
Definition: CouenneExpression.hpp:55
Couenne::t_chg_bounds::CHANGED
Definition: CouenneTypes.hpp:71
CouenneConfig.h
Couenne::exprMul::impliedBound
virtual bool impliedBound(int, CouNumber *, CouNumber *, t_chg_bounds *, enum Couenne::expression::auxSign=Couenne::expression::AUX_EQ)
implied bound processing
Definition: impliedBounds-exprMul.cpp:28
COUENNE_EPS
#define COUENNE_EPS
Definition: CouennePrecisions.hpp:19
Couenne::exprOp::arglist_
expression ** arglist_
argument list is an array of pointers to other expressions
Definition: CouenneExprOp.hpp:35
Couenne::CouNumber
double CouNumber
main number type in Couenne
Definition: CouenneTypes.hpp:100
COUENNE_INFINITY
#define COUENNE_INFINITY
Definition: CouennePrecisions.hpp:32
CouennePrecisions.hpp
Couenne::expression::Value
virtual CouNumber Value() const
value (empty)
Definition: CouenneExpression.hpp:108
c
real c
Definition: BonBqpdSolver.cpp:98
Generated by   doxygen 1.8.5