1 /* $Id$
2 *
3 * Name: branchExprTrilinear.cpp
4 * Author: Pietro Belotti
5 * Purpose: return branch data for trilinear terms
6 *
7 * This file is licensed under the Eclipse Public License (EPL)
8 */
9
10 #include "CouennePrecisions.hpp"
11 #include "CouenneTypes.hpp"
12 #include "CouenneObject.hpp"
13
14 #include "CouenneExprTrilinear.hpp"
15 #include "CouenneFunTriplets.hpp"
16 #include "CouenneProjections.hpp"
17
18 using namespace Couenne;
19
20 /// set up branching object by evaluating many branching points for
21 /// each expression's arguments
selectBranch(const CouenneObject * obj,const OsiBranchingInformation * info,expression * & var,double * & brpts,double * & brDist,int & way)22 CouNumber exprTrilinear::selectBranch (const CouenneObject *obj,
23 const OsiBranchingInformation *info,
24 expression *&var,
25 double * &brpts,
26 double * &brDist, // distance of current LP point
27 // to new convexifications
28 int &way) {
29
30 if (brDist) {free (brDist); brDist = NULL;} // clear it, computeMulBrDist will fill it
31
32 int
33 xi = arglist_ [0] -> Index (),
34 yi = arglist_ [1] -> Index (),
35 zi = arglist_ [2] -> Index ();
36
37 assert ((xi >= 0) && (yi >= 0) && (zi >= 0));
38
39 CouNumber
40 xl = info -> lower_ [xi], yl = info -> lower_ [yi], zl = info -> lower_ [zi],
41 xu = info -> upper_ [xi], yu = info -> upper_ [yi], zu = info -> upper_ [zi];
42
43 brpts = (double *) realloc (brpts, sizeof (double));
44 brDist = (double *) realloc (brDist, 2 * sizeof (double));
45
46 // case 0: term is fixed
47
48 if ((fabs (xu - xl) < COUENNE_EPS) &&
49 (fabs (xu - xl) < COUENNE_EPS) &&
50 (fabs (xu - xl) < COUENNE_EPS)) {
51
52 var = NULL;
53 return 0.;
54 }
55
56 // a very simple branching scheme:
57 //
58 // if any bound interval is (-inf,+inf), break it by branching at zero
59 // otherwise, if any bound interval is [a,+inf) or (-inf,b],
60 // break it by branching at zero (if interval includes zero) or
61 // at 2a-1 (resp. 2b+1)
62 // otherwise, branch on the largest bound, in the middle
63
64 if ((xl < -COUENNE_INFINITY) && (xu > COUENNE_INFINITY)) {*brpts = 0.; brDist [0] = brDist [1] = 1.; var = arglist_ [0]; return 1.;}
65 if ((yl < -COUENNE_INFINITY) && (yu > COUENNE_INFINITY)) {*brpts = 0.; brDist [0] = brDist [1] = 1.; var = arglist_ [1]; return 1.;}
66 if ((zl < -COUENNE_INFINITY) && (zu > COUENNE_INFINITY)) {*brpts = 0.; brDist [0] = brDist [1] = 1.; var = arglist_ [2]; return 1.;}
67
68
69 #define SETBNDS(l,u,ind) { \
70 \
71 if (l < -COUENNE_INFINITY) {\
72 if (u > 1.) {*brpts = 0.; brDist [0] = brDist [1] = 1.; var = arglist_ [ind]; return 1.;}\
73 else {*brpts = 2*-fabs (u) - 1.; brDist [0] = brDist [1] = 1.; var = arglist_ [ind]; return 1.;}\
74 }\
75 \
76 if (u > COUENNE_INFINITY) {\
77 if (l < -1.) {*brpts = 0.; brDist [0] = brDist [1] = 1.; var = arglist_ [ind]; return 1.;}\
78 else {*brpts = 2*fabs (u) + 1.; brDist [0] = brDist [1] = 1.; var = arglist_ [ind]; return 1.;}\
79 }\
80 }
81
82 SETBNDS (xl, xu, 0);
83 SETBNDS (xl, xu, 1);
84 SETBNDS (xl, xu, 2);
85
86 // bounds all finite, choose largest bound interval width
87 if ((xu - xl > yu - yl) && (xu - xl > zu - zl)) {*brpts = .5 * (xl + xu); brDist [0] = brDist [1] = 1.; var = arglist_ [0]; return 1.;}
88 else if ((yu - yl > zu - zl)) {*brpts = .5 * (yl + yu); brDist [0] = brDist [1] = 1.; var = arglist_ [1]; return 1.;}
89 else {*brpts = .5 * (zl + zu); brDist [0] = brDist [1] = 1.; var = arglist_ [2]; return 1.;}
90 }
91