1 /* $Id$
2 *
3 * Name: branchExprInv.cpp
4 * Author: Pietro Belotti
5 * Purpose: return branch selection for 1/x
6 *
7 * (C) Carnegie-Mellon University, 2006-07.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include "CoinHelperFunctions.hpp"
12
13 #include "CouenneExprInv.hpp"
14 #include "CouenneObject.hpp"
15 #include "CouenneBranchingObject.hpp"
16 #include "CouenneProjections.hpp"
17 #include "CouenneFunTriplets.hpp"
18
19 using namespace Couenne;
20
21 /// generic approach for negative powers (commom with exprInv::selectBranch())
negPowSelectBranch(const CouenneObject * obj,const OsiBranchingInformation * info,double * & brpts,double * & brDist,int & way,CouNumber k,CouNumber x0,CouNumber y0,CouNumber l,CouNumber u)22 CouNumber negPowSelectBranch (const CouenneObject *obj,
23 const OsiBranchingInformation *info,
24 double * &brpts,
25 double * &brDist, // distance of current LP
26 // point to new convexifications
27 int &way,
28 CouNumber k,
29 CouNumber x0, CouNumber y0,
30 CouNumber l, CouNumber u) {
31
32 brDist = (double *) realloc (brDist, 2 * sizeof (double));
33 brpts = (double *) realloc (brpts, sizeof (CouNumber));
34
35 // two cases: inside or outside the curves (there are two branches
36 // of the hyperbola).
37 //
38 // Inside: the distance depends on the projection of the current
39 // point onto the would-be upper envelopes, which forces us to look
40 // at it numerically. If both bounds are infinite, create a ThreeWay
41 // branch.
42 //
43 // Outside: it suffices to project the current point on the line
44 // (i.e. to get the closest point on the line) to get the maxi-min
45 // displacement.
46 //
47 // As for all monotonous functions, after choosing *brpts it is
48 // equivalent to choose w's or x's index as ind, as the implied- and
49 // propagated bounds will do the rest.
50
51 if ((l < -COUENNE_EPS) && (u > COUENNE_EPS)) { // handle discontinuity
52
53 // no matter if the (negative) exponent is odd or even, we better
54 // branch on 0 lest we have no good convexification (especially
55 // with odd exponent)
56
57 *brpts = 0.;
58 way = TWO_RAND;
59
60 // Closest branch of the hyperbola is on the same side of y+x=0 as
61 // (x0,y0) => need only one powNewton
62
63 if (fabs (x0) < COUENNE_EPS)
64 x0 = (x0 <= -0.) ? -COUENNE_EPS : COUENNE_EPS;
65
66 CouNumber xp, xx0 = x0, yy0 = y0, exponent = k;
67
68 // invert dependent and independent if
69 if (((x0+y0 < 0.) && (x0 > 0.)) || // in lower half of fourth orthant, or
70 ((x0+y0 > 0.) && (x0 < 0.))) { // in upper half of second orthant
71
72 exponent = 1. / k;
73 xx0 = y0;
74 yy0 = x0;
75 }
76
77 powertriplet pt (exponent);
78
79 xp = (xx0 >= 0) ?
80 powNewton (xx0, yy0, &pt) :
81 -powNewton (-xx0, -yy0, &pt);
82
83 CouNumber diff = x0 - xp;
84 y0 -= safe_pow (xp, 1. / k);
85
86 brDist [0] = sqrt (diff*diff + y0*y0); // exact distance
87 brDist [1] = CoinMax (fabs (x0), 1.);
88
89 if (x0 > 0.) {
90 double swap = brDist [0];
91 brDist [0] = brDist [1];
92 brDist [1] = swap;
93 }
94
95 return CoinMin (brDist [0], brDist [1]);
96 }
97
98 int intk = 0;
99
100 bool
101 isInt = fabs (k - (double) (intk = COUENNE_round (k))) < COUENNE_EPS,
102 isInvInt = !isInt && (fabs (1./k - (double) (intk = COUENNE_round (1./k))) < COUENNE_EPS);
103
104 // case 2: bound interval does not contain zero. Look if inside or
105 // outside of belly (refer to branchExprExp.cpp)
106
107 if (((x0 >= 0.) && (y0 < safe_pow (x0,k))) || // x0>0, or
108 ((x0 <= -0.) && // x0<0, and
109 (((isInt && !(intk % 2) && (y0 < safe_pow (x0,k)))) || // integer, even
110 (((isInt || isInvInt) && (intk % 2) && (y0 > -safe_pow (-x0,k))))))) { // (inv)integer, odd
111
112 // Inside. Branch on closest point on curve, computed with a
113 // Newton method
114
115 way = (u < -0.) ? TWO_RIGHT : TWO_LEFT; // explore finite interval first
116
117 powertriplet pt (k);
118
119 *brpts = obj -> midInterval ((x0 >= 0.) ?
120 powNewton ( x0, y0, &pt) :
121 -powNewton (-x0, -y0, &pt), l, u, info);
122
123 CouNumber dy = y0 - safe_pow (*brpts >= 0 ? *brpts : - *brpts, 1. / k);
124 x0 -= *brpts;
125 return (brDist [0] = brDist [1] = sqrt (x0*x0 + dy*dy)); // distance is exact
126 }
127
128 // Inside, (x0^k) * y0 >= 1. Two cases: /////////////////////////////////////////////////
129
130 // 1) bounds are infinite both horizontally and vertically
131 // (i.e. [-inf,0] or [0,inf]) --> as for exprExp, pick point on
132 // diagonal from current to curve, to be sure current will be cut by
133 // branching rule
134
135 if (((l < COUENNE_EPS) && (u > COUENNE_INFINITY)) ||
136 ((u > - COUENNE_EPS) && (l < - COUENNE_INFINITY))) {
137
138 /* brpts = (double *) realloc (brpts, 2 * sizeof (double));
139 way = THREE_CENTER; // focus on central convexification first
140 brpts [0] = x0; // draw vertical from (x0,y0) south (north) to curve y=1/x
141 brpts [1] = 1. / y0; // horizontal west (east)
142 CouNumber a = fabs (y0 - 1 / x0), // sides of a triangle with (x0,y0)
143 b = fabs (x0 - 1 / y0), // as one of the vertices
144 c = a * cos (atan (a/b)); */
145
146 //brpts = (double *) realloc (brpts, sizeof (double));
147
148 //if (x0 > COUENNE_EPS)
149 *brpts = 0.5 * (fabs (x0) + pow (fabs (y0), 1./k));
150
151 if (x0 < 0.) {
152 *brpts = - *brpts;
153 brDist [0] = fabs (fabs (y0) - safe_pow (fabs (x0), k));
154 brDist [1] = *brpts - x0;
155 } else {
156 brDist [0] = x0 - *brpts;
157 brDist [1] = fabs (y0 - safe_pow (x0, k));
158 }
159
160 //else
161 //*brpts = 0.5 * (x0 + pow (y0, 1./k));
162
163 // follow South-East diagonal to find point on curve
164 // so that current point is surely cut
165 //*brpts = 0.5 * (x0 + log (y0));
166 //way = TWO_RAND;
167 way = (x0 > *brpts) ? TWO_RIGHT : TWO_LEFT;
168
169 return CoinMin (brDist [0], brDist [1]);
170 //x0 - pow (fabs (y0), 1./k), y0 - pow (x0,k));
171 //return CoinMin (a, CoinMin (b, c)); // distance is exact
172 }
173
174 // 2) at least one of them is finite
175
176 if (l < - COUENNE_INFINITY) { // u << -0
177
178 way = TWO_RIGHT;
179 *brpts = obj -> midInterval (x0, l, u, info);
180
181 return CoinMin (brDist [0] = y0 - safe_pow (*brpts, 1. / k),
182 brDist [1] = projectSeg (x0, y0, l, safe_pow (l, k),
183 *brpts, safe_pow (*brpts, k), -1)); // distance is exact
184 }
185
186 if (u > COUENNE_INFINITY) { // l >> +0
187
188 way = TWO_LEFT;
189 *brpts = obj -> midInterval (x0, l, u, info);
190
191 return CoinMin (brDist [1] = y0 - safe_pow (*brpts, 1. / k),
192 brDist [0] = projectSeg (x0, y0, l, safe_pow (l, k),
193 *brpts, safe_pow (*brpts, k), +1)); // distance is exact
194 }
195
196 // last case: nice finite interval and limited curve
197
198 powertriplet ft (k);
199 *brpts = obj -> getBrPoint (&ft, x0, l, u, info);
200
201 /* // TODO: check if it works with all exponents
202 if (u > l + COUENNE_EPS) {
203
204 powertriplet ft (k);
205 *brpts = maxHeight (&ft, l, u); // min area
206
207 // *brpts = safe_pow ((safe_pow (u,k) - safe_pow (l,k)) / (k * (u-l)), 1/(k-1));
208 // if (u < 0)
209 // *brpts = - *brpts;
210 }
211 else *brpts = midInterval (x0, l, u, info);*/
212
213 way = TWO_RAND;
214
215 x0 -= *brpts;
216 y0 -= safe_pow (*brpts, k);
217
218 brDist [0] = projectSeg (x0,y0, l, safe_pow (l, k), *brpts, safe_pow (*brpts, k), 0);
219 brDist [1] = projectSeg (x0,y0, *brpts, safe_pow (*brpts, k), u, safe_pow (u, k), 0);
220
221 return CoinMin (brDist [0], brDist [1]);//sqrt (x0*x0 + y0*y0); // distance is exact
222 }
223
224
225
226 /// set up branching object by evaluating many branching points for
227 /// each expression's arguments
228
selectBranch(const CouenneObject * obj,const OsiBranchingInformation * info,expression * & var,double * & brpts,double * & brDist,int & way)229 CouNumber exprInv::selectBranch (const CouenneObject *obj,
230 const OsiBranchingInformation *info,
231 expression *&var,
232 double * &brpts,
233 double * &brDist, // distance of current LP
234 // point to new convexifications
235 int &way) {
236
237 var = argument_;
238
239 int
240 ind = argument_ -> Index (),
241 wi = obj -> Reference () -> Index ();
242
243 assert ((ind >= 0) && (wi >= 0));
244
245 CouNumber y0 = info -> solution_ [wi],
246 x0 = info -> solution_ [ind],
247 l = info -> lower_ [ind],
248 u = info -> upper_ [ind];
249
250 return negPowSelectBranch (obj, info, brpts, brDist, way, -1, x0, y0, l, u);
251 }
252