branchExprLog.cpp
Go to the documentation of this file.
1 /* $Id: branchExprLog.cpp 708 2011-06-23 14:04:59Z pbelotti $
2  *
3  * Name: branchExprLog.cpp
4  * Author: Pietro Belotti
5  * Purpose: return branch gain and branch object for logarithms
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 "CouenneExprLog.hpp"
14 #include "CouenneObject.hpp"
16 #include "CouenneProjections.hpp"
17 #include "CouenneFunTriplets.hpp"
18 
19 using namespace Couenne;
20 
21 #define SQ_COUENNE_EPS COUENNE_EPS * COUENNE_EPS
22 
26  const OsiBranchingInformation *info,
27  expression *&var,
28  double * &brpts,
29  double * &brDist, // distance of current LP
30  // point to new convexifications
31  int &way) {
32 
33  brpts = (double *) realloc (brpts, sizeof (double));
34  brDist = (double *) realloc (brDist, 2*sizeof (double));
35 
36  // quite similar to exprExp::selectBranch() (see branchExprExp.cpp)
37  //
38  // two cases: inside or outside the belly.
39  //
40  // Inside: the distance depends on the projection of the current
41  // point onto the would-be upper envelopes, which forces us to look
42  // at it numerically. If both bounds are infinite, create a ThreeWay
43  // branch.
44  //
45  // Outside: it suffices to project the current point on the line
46  // (i.e. to get the closest point on the line) to get the maxi-min
47  // displacement.
48  //
49  // As for all monotonous functions, after choosing *brpts it is
50  // equivalent to choose w's or x's index as ind, as the implied- and
51  // propagated bounds will do the rest.
52 
53  var = argument_;
54 
55  int
56  ind = var -> Index (),
57  wi = obj -> Reference () -> Index ();
58 
59  assert ((ind >= 0) && (wi >= 0));
60 
61  CouNumber
62  y0 = info -> solution_ [wi],
63  x0 = info -> solution_ [ind],
64  l = info -> lower_ [ind],
65  u = info -> upper_ [ind];
66 
67  if (u < COUENNE_EPS) { // strange case, return default branching rule
68  var = NULL;
69  return 0.;
70  }
71 
72  if (x0 < SQ_COUENNE_EPS) // very unlikely...
73  x0 = SQ_COUENNE_EPS;
74 
75  if (y0 > log (x0)) {
76 
77  // Outside -> branch on closest point on curve
78 
79  *brpts = obj -> midInterval (powNewton (x0, y0, log, inv, oppInvSqr), l, u, info);
80 
81  way = TWO_LEFT;
82  CouNumber dy = y0 - log (*brpts);
83  x0 -= *brpts;
84 
85  return (brDist [0] = brDist [1] = sqrt (x0*x0 + dy*dy)); // exact distance
86  }
87 
88  // Inside. Two cases: ////////////////////////////////////////////////////
89 
90  if ((l <= SQ_COUENNE_EPS) &&
91  (u > COUENNE_INFINITY)) {
92 
93  // 1) curve is unlimited in both senses
94 
95  /*// restore when we can do three-way branching
96  #if 0
97  brpts = (double *) realloc (brpts, 2 * sizeof (double));
98  way = THREE_CENTER; // focus on central convexification first
99  brpts [0] = exp (y0); // draw horizontal from (x0,y0) east to curve y=log(x)
100  brpts [1] = x0; // vertical north
101  CouNumber a = x0 - exp (y0), // sides of a triangle with (x0,y0)
102  b = log (x0) - y0; // as one of the vertices
103 
104  return a * cos (atan (a/b)); // exact distance
105  #endif*/
106 
107  // follow South-East diagonal to find point on curve
108  // so that current point is surely cut
109  brpts = (double *) realloc (brpts, sizeof (double));
110  *brpts = 0.5 * (x0 + exp (y0));
111  way = TWO_RAND;
112 
113  return CoinMin (brDist [0] = x0 - exp (y0),
114  brDist [1] = log (x0) - y0);
115  }
116 
117  // 2) at least one of them is finite
118 
119  brpts = (double *) realloc (brpts, sizeof (double));
120 
121  if (l <= SQ_COUENNE_EPS) { // u is finite
122 
123  *brpts = obj -> midInterval (exp (y0), l, u, info);
124  way = TWO_RIGHT;
125 
126  brDist [0] = x0 - *brpts;
127  return (brDist [1] = projectSeg (x0, y0, *brpts, log (*brpts), u, log (u), +1)); // exact distance
128  // return CoinMin (x0 - exp (y0), log (x0) - y0);
129  }
130 
131  if (u > COUENNE_INFINITY) { // l is far from zero
132 
133  *brpts = CoinMax (10*x0, obj -> midInterval (x0, l, u, info));
134  way = TWO_LEFT;
135 
136  brDist [1] = log (*brpts) - y0;
137 
138  return (brDist [0] = projectSeg (x0, y0, l, log (l), *brpts, log (*brpts), +1)); // exact distance
139  //return log (x0) - y0;
140  }
141 
142  // both are finite
143 
145 
146  *brpts = obj -> getBrPoint (&ft, x0, l, u, info);
147 
148  // *brpts = midInterval (powNewton (x0, y0, log, inv, oppInvSqr), l, u, info);
149  // WRONG! Local minima may be at bounds
150 
151  // compute distance from current point to new convexification(s) and
152  // to curve. If closer to curve, branch on current point
153 
154  way = TWO_RAND;
155 
156  // exact distance
157  return CoinMin (brDist [0] = projectSeg (x0, y0, l, log (l), *brpts, log (*brpts), +1),
158  brDist [1] = projectSeg (x0, y0, *brpts, log (*brpts), u, log (u), +1));
159 }
void fint fint fint real fint real real real real real real real real real fint real fint fint fint real fint fint fint fint * info
OsiObject for auxiliary variables $w=f(x)$.
CouExpr & log(CouExpr &e)
CouNumber powNewton(CouNumber xc, CouNumber yc, unary_function f, unary_function fp, unary_function fpp)
find proper tangent point to add deepest tangent cut
Definition: powNewton.cpp:29
CouNumber oppInvSqr(register CouNumber x)
derivative of inv (x)
virtual CouNumber selectBranch(const CouenneObject *obj, const OsiBranchingInformation *info, expression *&var, double *&brpts, double *&brDist, int &way)
set up branching object by evaluating many branching points for each expression&#39;s arguments ...
expression * argument_
single argument taken by this expression
CouNumber inv(register CouNumber arg)
the operator itself
ULong * x0
Definition: OSdtoa.cpp:1776
CouExpr & exp(CouExpr &e)
virtual int Index() const
Return index of variable (only valid for exprVar and exprAux)
#define COUENNE_EPS
CouNumber projectSeg(CouNumber x0, CouNumber y0, CouNumber x1, CouNumber y1, CouNumber x2, CouNumber y2, int sign, CouNumber *xp=NULL, CouNumber *yp=NULL)
Compute projection of point (x0, y0) on the segment defined by two points (x1,y1), (x2, y2) – sign provided by parameter sign.
Definition: projections.cpp:78
double CouNumber
main number type in Couenne
#define COUENNE_INFINITY
Expression base class.
#define SQ_COUENNE_EPS