1 /* $Id$
2 *
3 * Name: branchExprQuad.cpp
4 * Author: Pietro Belotti
5 * Purpose: return branch data for quadratic forms
6 *
7 * (C) Carnegie-Mellon University, 2007-10.
8 * This file is licensed under the Eclipse Public License (EPL)
9 */
10
11 #include "CoinHelperFunctions.hpp"
12
13 #include "CouenneExprQuad.hpp"
14 #include "CouenneObject.hpp"
15 #include "CouenneBranchingObject.hpp"
16
17 using namespace Couenne;
18
19 //#define DEBUG
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 exprQuad::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 int ind = -1;
31
32 // use a combination of eigenvectors and bounds
33
34 CouNumber delta = (*(obj -> Reference ())) () - (*this) ();
35
36 /*printf ("infeasibility: ");
37 obj -> Reference () -> print ();
38 printf (" [%g=%g] := ",
39 (*(obj -> Reference ())) (), info -> solution_ [obj -> Reference () -> Index ()]);
40 print (); printf (" [%g]\n", (*this) ());*/
41
42 brpts = (double *) realloc (brpts, sizeof (double));
43 brDist = (double *) realloc (brDist, 2*sizeof (double));
44
45 way = TWO_RAND;
46
47 // depending on where the current point is w.r.t. the curve,
48 // branching is done on the eigenvector corresponding to the minimum
49 // or maximum eigenvalue
50
51 std::vector <std::pair <CouNumber,
52 std::vector <std::pair <exprVar *, CouNumber> > > >::iterator fi = eigen_.begin ();
53
54 std::vector <std::pair <CouNumber,
55 std::vector <std::pair <exprVar *, CouNumber> > > >::reverse_iterator ri = eigen_.rbegin ();
56
57 CouNumber max_span = -COUENNE_INFINITY;
58
59 bool changed_sign = false;
60
61 /////////////////////////////////////////////////////////
62
63 for (;(((delta < 0.) && (fi != eigen_. end ())) || // && (fi -> first < 0.) ||
64 ((delta > 0.) && (ri != eigen_. rend ()))); // && (ri -> first > 0.));
65 ++fi, ++ri) {
66
67 std::vector <std::pair <exprVar *, CouNumber> > &ev =
68 (delta < 0.) ? fi -> second : ri -> second;
69
70 if (((delta < 0.) && (fi -> first > 0.)) ||
71 ((delta > 0.) && (ri -> first < 0.))) {
72
73 if (max_span > 0.) break; // if found a variable already, return
74 changed_sign = true; // otherwise, keep in mind we are on
75 // the wrong eigenvalues' sign
76 }
77
78 for (std::vector <std::pair <exprVar *, CouNumber> >::iterator j = ev.begin ();
79 j != ev.end (); ++j) {
80
81 int index = j -> first -> Index ();
82
83 CouNumber
84 lb = info -> lower_ [index],
85 ub = info -> upper_ [index];
86
87 // only accept variable if not fixed
88 if (fabs (ub-lb) > COUENNE_EPS) {
89
90 // no variable was found but the eigenvalue is already
91 // positive (negative)
92 // changed_sign &&
93 // (max_span < 0.))
94
95 CouNumber span = -1;
96
97 if ((lb < -COUENNE_INFINITY) ||
98 (ub > COUENNE_INFINITY) ||
99 ((span = (ub-lb) * fabs (j -> second)) > max_span + COUENNE_EPS)) {
100
101 ind = index;
102 var = j -> first;
103
104 *brpts = obj -> midInterval (info -> solution_ [index], lb, ub, info);
105
106 if (changed_sign)
107 break;
108
109 if (span >= 0) max_span = span; // finite, keep searching
110 else break; // span unbounded, stop searching
111 }
112 }
113 }
114 }
115
116 if ((eigen_.size () == 0) // if no eigenvalue has been computed yet
117 || (ind < 0)) { // or no index was found, pick largest interval
118
119 CouNumber max_span = -COUENNE_INFINITY;
120
121 for (std::map <exprVar *, std::pair <CouNumber, CouNumber> >::iterator i = bounds_. begin ();
122 i != bounds_. end (); ++i) {
123
124 CouNumber
125 lb = i -> second.first,
126 ub = i -> second.second,
127 span = ub - lb;
128
129 if ((span > COUENNE_EPS) &&
130 (span > max_span)) {
131
132 var = i -> first;
133 ind = var -> Index ();
134 }
135 }
136
137 if (ind < 0) {
138
139 var = obj -> Reference ();
140 ind = var -> Index ();
141
142 *brpts = obj -> midInterval (info -> solution_ [ind],
143 info -> lower_ [ind],
144 info -> upper_ [ind], info);
145 }
146 else *brpts = obj -> midInterval (info -> solution_ [ind],
147 info -> lower_ [ind],
148 info -> upper_ [ind], info);
149 //return fabs (delta);
150 }
151
152 return (brDist [0] = brDist [1] = fabs (delta));
153 }
154