1 /* $Id$
2 *
3 * Name: exprInv.cpp
4 * Author: Pietro Belotti
5 * Purpose: definition of inverse of a function (1/f(x))
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> // ! must go
12
13 #include "CouenneExprInv.hpp"
14 #include "CouenneExprClone.hpp"
15 #include "CouenneExprMul.hpp"
16 #include "CouenneExprOpp.hpp"
17 #include "CouenneExprDiv.hpp"
18 #include "CouenneExprPow.hpp"
19 #include "CouenneProblem.hpp"
20 #include "CouenneExpression.hpp"
21
22 #include "CoinFinite.hpp"
23
24 using namespace Couenne;
25
26 // differentiation
differentiate(int index)27 expression *exprInv::differentiate (int index) {
28
29 return new exprOpp (new exprDiv (argument_ -> differentiate (index),
30 new exprPow (new exprClone (argument_),
31 new exprConst (2.))));
32 }
33
34
35 // printing
print(std::ostream & out,bool descend) const36 void exprInv::print (std::ostream &out,
37 bool descend) const {
38 out << "(1/";
39 argument_ -> print (out, descend);
40 out << ")";
41 }
42 // exprUnary::print (out, "1/", PRE);}
43
44
45 /// general function to tighten implied bounds of a function w = x^k,
46 /// k negative, integer or inverse integer, and odd
invPowImplBounds(int wind,int index,CouNumber * l,CouNumber * u,CouNumber k,bool & resL,bool & resU,enum expression::auxSign sign)47 void invPowImplBounds (int wind, int index,
48 CouNumber *l, CouNumber *u, CouNumber k,
49 bool &resL, bool &resU,
50 enum expression::auxSign sign) {
51
52 CouNumber wl = sign == expression::AUX_GEQ ? -COIN_DBL_MAX : l [wind],
53 wu = sign == expression::AUX_LEQ ? COIN_DBL_MAX : u [wind];
54
55 // 0 <= l <= w <= u
56
57 if (wl >= 0.) {
58 if (wu > COUENNE_EPS) {
59 if (wu < COUENNE_INFINITY) resL = updateBound (-1, l + index, pow (wu, k));
60 else resL = updateBound (-1, l + index, 0.);
61 }
62 if (wl > COUENNE_EPS) resU = updateBound (+1, u + index, pow (wl, k));
63 }
64
65 // l <= w <= u <= 0
66
67 if (wu <= -0.) {
68 if (wl < - COUENNE_EPS) {
69 if (wl > - COUENNE_INFINITY) resU = updateBound (+1, u + index, pow (wl, k)) || resU;
70 else resU = updateBound (+1, u + index, 0.) || resU;
71 }
72 if (wu < - COUENNE_EPS) resL = updateBound (-1, l + index, pow (wu, k)) || resL;
73 }
74 }
75
76
77 /// implied bound processing for expression w = 1/x, upon change in
78 /// lower- and/or upper bound of w, whose index is wind
impliedBound(int wind,CouNumber * l,CouNumber * u,t_chg_bounds * chg,enum auxSign sign)79 bool exprInv::impliedBound (int wind, CouNumber *l, CouNumber *u, t_chg_bounds *chg, enum auxSign sign) {
80
81 // Expression w = 1/x: we can only improve the bounds if
82 //
83 // 0 <= l <= w <= u or
84 // l <= w <= u <= 0.
85 //
86 // Then 1/u <= x <= 1/l (given l, u finite and nonzero)
87
88 int index = argument_ -> Index ();
89
90 bool resL, resU = resL = false;
91
92 invPowImplBounds (wind, index, l, u, -1., resL, resU, sign);
93
94 bool argInt = argument_ -> isInteger ();
95
96 if (resL) {
97 chg [index].setLower(t_chg_bounds::CHANGED);
98 if (argInt) l [index] = ceil (l [index] - COUENNE_EPS);
99 }
100
101 if (resU) {
102 chg [index].setUpper(t_chg_bounds::CHANGED);
103 if (argInt) u [index] = floor (u [index] + COUENNE_EPS);
104 }
105
106 return (resL || resU);
107 }
108
109
110 /// return l-2 norm of gradient at given point
gradientNorm(const double * x)111 CouNumber exprInv::gradientNorm (const double *x) {
112 int ind = argument_ -> Index ();
113 CouNumber xx;
114 if (ind < 0) return 0.;
115 else {
116 xx = x [argument_ -> Index ()];
117 return 1. / (xx*xx);
118 }
119 }
120
121
122 /// can this expression be further linearized or are we on its
123 /// concave ("bad") side
isCuttable(CouenneProblem * problem,int index) const124 bool exprInv::isCuttable (CouenneProblem *problem, int index) const {
125
126 int xind = argument_ -> Index ();
127
128 double
129 x = problem -> X (xind),
130 y = problem -> X (index);
131
132 return (((problem -> Lb (xind) >= 0) && (x > 0) && (y*x <= 1)) ||
133 ((problem -> Ub (xind) <= 0) && (x < 0) && (y*x <= 1)));
134 }
135