source: stable/0.2/Couenne/src/convex/operators/conv-exprPow-getBounds.cpp @ 159

Last change on this file since 159 was 159, checked in by pbelotti, 11 years ago

created new stable branch 0.2 from trunk (rev. 157)

File size: 7.2 KB
Line 
1/* $Id: conv-exprPow-getBounds.cpp 139 2009-06-03 03:56:30Z pbelotti $
2 *
3 * Name:    conv-exprPow-getBounds.cpp
4 * Author:  Pietro Belotti
5 * Purpose: method to get lower and upper bounds of a power x^y
6 *
7 * (C) Carnegie-Mellon University, 2006-09.
8 * This file is licensed under the Common Public License (CPL)
9 */
10
11#include <math.h>
12
13#include "CouenneTypes.hpp"
14#include "exprPow.hpp"
15#include "exprConst.hpp"
16#include "exprClone.hpp"
17#include "exprMax.hpp"
18#include "exprMin.hpp"
19#include "exprOpp.hpp"
20#include "CouennePrecisions.hpp"
21#include "CouenneProblem.hpp"
22
23
24// compute expressions for lower and upper bounds of a power x^y,
25// based on the lower/upper bounds of x and y
26
27void exprPow::getBounds (expression *&lb, expression *&ub) {
28
29  // We have a standardized expression of the form w = x^y, where x or
30  // y could be constant. Let us study each case separately.
31
32  assert (arglist_ [0] -> Type () != CONST);
33
34  // x is not constant, so it has (possibly different) lower and
35  // upper bounds. The expression is x^b, b constant (the case x^y
36  // has been decomposed by simplify() into exp(y log x).
37
38  expression *lbbase, *ubbase;
39  arglist_ [0] -> getBounds (lbbase, ubbase);
40
41  //    printf ("ubbase = "); ubbase -> print (std::cout); printf ("\n");
42
43  if (arglist_ [1] -> Type () == CONST) { 
44
45    // expression = x^b, b!=0. There are four cases:
46    //
47    // 1) b   is integer and odd  (cube, x^5, etc)
48    // 2) b   is integer and even (square, x^8, etc)
49    // 3) 1/b is integer and odd  (cube root, x^(1/7), etc)
50    // 4) 1/b is integer and even (square root, x^(1/4), etc)
51    // 5) none of the above
52    //
53    // For all of these, need to check if the exponent is negative...
54
55    CouNumber expon = arglist_ [1] -> Value ();
56    int rndexp;
57
58    bool isInt    =  fabs (expon - (rndexp = COUENNE_round (expon))) < COUENNE_EPS,
59      isInvInt = !isInt &&
60      ((fabs (expon) > COUENNE_EPS) &&
61       (fabs (1/expon - (rndexp = COUENNE_round (1/expon))) < COUENNE_EPS));
62
63    if ((isInt || isInvInt) && (rndexp % 2) && (rndexp > 0)) { 
64
65      // the exponent is integer (or inverse integer), odd and
66      // positive, hence the function is monotone non decreasing
67
68      lb = new exprPow (lbbase, new exprConst (expon));
69      ub = new exprPow (ubbase, new exprConst (expon));
70    } 
71    else {
72
73      // the exponent is either negative, integer even, or fractional
74
75      expression **all = new expression * [6];
76
77      all [0] = new exprOpp   (lbbase);
78      all [2] = new exprConst (0.);
79      all [4] = ubbase;
80
81      if (expon > 0) 
82        all    [1] = new exprPow (new exprClone (lbbase), new exprConst (expon));
83      else all [1] = new exprPow (new exprClone (ubbase), new exprConst (expon));
84
85      // all [3] is lower bound when lbbase <= 0 <= ubbase
86
87      if (expon > COUENNE_EPS) all [3] = new exprConst (0.);
88      else if (isInt || isInvInt) {
89        if (rndexp % 2) 
90          all [3] = new exprConst (-COUENNE_INFINITY);
91        else all [3] = new exprMin (new exprClone (all [1]),
92                                    new exprPow (new exprClone (lbbase), 
93                                                 new exprConst (expon)));
94      }
95      else all [3] = new exprClone (all [1]);
96
97      // all [5] is the lower bound value when lbbase <= ubbase <= 0
98
99      if (expon > COUENNE_EPS) {
100        if (isInt && !(rndexp % 2))
101          all [5] = new exprPow (new exprClone (ubbase), new exprConst (expon));
102        else all [5] = new exprConst (0.);
103      }
104      else {
105        if (isInt || isInvInt) {
106          if (rndexp % 2)
107            all    [5] = new exprPow (new exprClone (ubbase), new exprConst (expon));
108          else all [5] = new exprPow (new exprClone (lbbase), new exprConst (expon));
109        }
110        else all [5] = new exprConst (0.);
111      }
112
113      lb = new exprMin (all, 6);
114
115      // And now the upper bound ///////////////////////////////////
116
117      if (expon > 0) {
118
119        // special case: upper bound depends to variable bounds only:
120        // $max {lb^k, ub^k}$
121
122        ub = new exprMax (new exprPow (new exprClone (lbbase), new exprConst (expon)),
123                          new exprPow (new exprClone (ubbase), new exprConst (expon)));
124
125      } else { // from this point on, expon < 0
126
127        expression **alu = new expression * [6];
128
129        alu [0] = new exprClone (all [0]);
130        alu [2] = new exprConst (0.);
131        alu [4] = new exprClone (ubbase);
132
133        //if ((isInt || isInvInt) && !(rndexp % 2))
134        //alu    [1] = new exprPow (new exprClone (ubbase), new exprConst (expon));
135        //else
136
137        // if negative exponent and base has nonnegative lower bound,
138        // the upper bound can only be lb^k
139        alu [1] = new exprPow (new exprClone (lbbase), new exprConst (expon));
140
141        // alu [3] is upper bound when lbbase <= 0 <= ubbase
142
143        //if (expon < - COUENNE_EPS)
144        alu [3] = new exprConst (COUENNE_INFINITY);
145        //else if (isInt && !(rndexp % 2))
146        //alu [3] = new exprPow (new exprMax (new exprClone (lbbase), new exprClone (ubbase)),
147        //new exprConst (expon));
148        //else alu [3] = new exprPow (new exprClone (ubbase), new exprConst (expon));
149
150        // alu [5] is the upper bound value when lbbase <= ubbase <= 0
151
152        /*if (expon > COUENNE_EPS) {
153
154          if (isInt && !(rndexp % 2))
155            alu [5] = new exprPow (new exprClone(ubbase), new exprConst(expon));
156          else alu [5] = new exprConst (0.);
157        }
158        else {*/
159        if (isInt || isInvInt) 
160          alu [5] = new exprPow (new exprClone (ubbase), new exprConst (expon));
161        else alu [5] = new exprConst (COUENNE_INFINITY);
162          //}
163
164        ub = new exprMin (alu, 6);
165      }
166    }
167  }
168  else // should NOT happen, exponent is not constant...
169    printf ("exprPow::getBounds(): Warning, exponent not constant\n");
170
171  /*CouNumber l, u;
172  arglist_ [0] -> getBounds (l,u);
173
174  printf ("pow::bound: [");
175  lb -> print (); printf ("=%g, ", (*lb) ());
176  ub -> print (); printf ("=%g [%g,%g]\n", (*ub) (), l, u);*/
177}
178
179
180// get value of lower and upper bound for the expression
181void exprPow::getBounds (CouNumber &lb, CouNumber &ub) {
182
183  CouNumber lba, uba, k = (*(arglist_ [1])) ();
184  arglist_ [0] -> getBounds (lba, uba);
185  int intk;
186
187  bool 
188    isInt    =           fabs (k    - (double) (intk = COUENNE_round (k)))    < COUENNE_EPS,
189    isInvInt = !isInt && fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS;
190
191  if (!isInt && (!isInvInt || !(intk % 2))) {
192
193    // if exponent is fractional or 1/even, the base better be nonnegative
194
195    if (lba < 0.) lba = 0.;
196    if (uba < 0.) uba = 0.;
197  }
198
199  if (isInt && !(intk % 2) && (k > 0)) { // x^{2h}
200
201    if (uba < 0) {
202      lb = safe_pow (-uba, k);
203      ub = safe_pow (-lba, k);
204    } else if (lba > 0) {
205      lb = safe_pow (lba, k);
206      ub = safe_pow (uba, k);
207    } else {
208      lb = 0;
209      ub = safe_pow (CoinMax (-lba, uba), k);
210    }
211
212  } else if (k > 0) { // monotone increasing: x^{2h+1} with h integer, or x^h with h real
213
214    lb = safe_pow (lba, k);
215    ub = safe_pow (uba, k);
216
217  } else if (isInt && !(intk % 2)) { // x^{-2h} or x^{-1/2h} with h integer
218
219    if (uba < 0) {
220      lb = safe_pow (-lba, k);
221      ub = safe_pow (-uba, k);
222    } else if (lba > 0) {
223      lb = safe_pow (uba, k);
224      ub = safe_pow (lba, k);
225    } else {
226      lb = safe_pow (CoinMax (-lba, uba), k);
227      ub = COUENNE_INFINITY;
228    }
229
230  } else { // x^k, k<0
231
232    if (uba < 0) {
233      lb = safe_pow (uba, k);
234      ub = safe_pow (lba, k);
235    } else if (lba > 0) {
236      lb = safe_pow (uba, k);
237      ub = safe_pow (lba, k);
238    } else {
239      lb = -COIN_DBL_MAX; // !!! may not be reached
240      ub =  COIN_DBL_MAX;
241    }
242  }
243}
Note: See TracBrowser for help on using the repository browser.