1 /* $Id$
2 *
3 * Name: branchExprExp.cpp
4 * Author: Pietro Belotti
5 * Purpose: return branch gain and branch object for exponentials
6 *
7 * (C) Carnegie-Mellon University, 2006-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include "CoinHelperFunctions.hpp"
12
13 #include "CouenneExprExp.hpp"
14 #include "CouenneObject.hpp"
15 #include "CouenneBranchingObject.hpp"
16 #include "CouenneProjections.hpp"
17 #include "CouenneFunTriplets.hpp"
18
19 using namespace Couenne;
20
21 /// set up branching object by evaluating many branching points for
22 /// each expression's arguments
selectBranch(const CouenneObject * obj,const OsiBranchingInformation * info,expression * & var,double * & brpts,double * & brDist,int & way)23 CouNumber exprExp::selectBranch (const CouenneObject *obj,
24 const OsiBranchingInformation *info,
25 expression *&var,
26 double * &brpts,
27 double * &brDist, // distance of current LP
28 // point to new convexifications
29 int &way) {
30
31 // two cases: inside and outside the curve.
32 //
33 // Inside: the distance depends on the projection of the current
34 // point onto the would-be upper envelopes, which forces us to look
35 // at it numerically. If both bounds are infinite, create a ThreeWay
36 // branch.
37 //
38 // Outside: it suffices to project the current point on the line
39 // (i.e. to get the closest point on the line) to get the maxi-min
40 // displacement.
41 //
42 // As for all monotone functions, after choosing *brpts it is
43 // equivalent to choose w's or x's index as ind, as the implied- and
44 // propagated bounds will do the rest.
45
46 var = argument_;
47
48 brDist = (double *) realloc (brDist, 2 * sizeof (double));
49 brpts = (double *) realloc (brpts, sizeof (double));
50
51 int
52 ind = var -> Index (),
53 wi = obj -> Reference () -> Index ();
54
55 assert ((ind >= 0) && (wi >= 0));
56
57 CouNumber y0 = info -> solution_ [wi],
58 x0 = info -> solution_ [ind],
59 l = info -> lower_ [ind],
60 u = info -> upper_ [ind];
61
62 // Outside //////////////////////////////////////////////////////////////////
63
64 if (y0 < exp (x0)) {
65
66 // Look for point (x1,y1) on curve y=exp(x) closest to current
67 // (x0,y0), branch at x1
68
69 *brpts = obj -> midInterval (powNewton (x0, y0, exp, exp, exp), l, u, info);
70
71 way = TWO_RAND;
72
73 y0 -= exp (*brpts);
74 x0 -= *brpts;
75
76 return sqrt (brDist [0] = brDist [1] = sqrt (x0*x0 + y0*y0)); // exact distance
77 }
78
79 // Inside. Four cases: ///////////////////////////////////////////////////////
80
81 if ((l < -COUENNE_INFINITY) &&
82 (u > COUENNE_INFINITY)) { // unbounded in both directions
83
84 /* // TODO: restore when we can do three-way branching
85 #if 0
86 brpts = (double *) realloc (brpts, 2 * sizeof (double));
87 way = THREE_CENTER; // focus on central convexification first
88 brpts [0] = x0; // draw vertical from (x0,y0) south to curve y=exp(x)
89 brpts [1] = log (y0); // horizontal east
90 CouNumber
91 a = y0 - exp (x0), // sides of a triangle with (x0,y0)
92 b = log (y0) - x0; // as one of the vertices
93 // exact distance from central interval, from others it's a and b
94 return (a * cos (atan (a/b)));
95 #endif */
96
97 // follow South-East diagonal to find point on curve
98 // so that current point is surely cut
99 *brpts = 0.5 * (x0 + log (y0));
100 way = TWO_RAND;
101
102 return CoinMin (brDist [0] = log (y0) - x0,
103 brDist [1] = y0 - exp (x0));
104 }
105
106 // 2,3) at least one of them is finite
107
108 if (l < - COUENNE_INFINITY) { // 2) unbounded from below --> break vertically
109
110 *brpts = obj -> midInterval (x0, l, u, info);
111
112 way = TWO_RIGHT;
113 return CoinMin (brDist [0] = y0 - exp (x0),
114 brDist [1] = projectSeg (x0, y0, *brpts, exp (*brpts), u, exp (u), -1));
115 }
116
117 if (u > COUENNE_INFINITY) { // 3) unbounded from above -- break horizontally
118
119 *brpts = obj -> midInterval (log (y0), l, u, info);
120
121 way = TWO_LEFT;
122 return CoinMin (brDist [0] = projectSeg (x0, y0, l, exp (l), *brpts, exp (*brpts), -1),
123 brDist [1] = log (y0) - x0);
124
125 }
126
127 // 4) both are finite
128
129 simpletriplet ft (exp, exp, exp, log);
130 *brpts = obj -> getBrPoint (&ft, x0, l, u, info); // select based on strategy
131
132 way = TWO_RAND;
133
134 // exact distance
135 return CoinMin (brDist [0] = projectSeg (x0, y0, l, exp (l), *brpts, exp (*brpts), -1),
136 brDist [1] = projectSeg (x0, y0, *brpts, exp (*brpts), u, exp (u), -1));
137 }
138