1 /* $Id$
2 *
3 * Name: CouenneBranchingObject.cpp
4 * Authors: Pierre Bonami, IBM Corp.
5 * Pietro Belotti, Carnegie Mellon University
6 * Purpose: Branching object for auxiliary variables
7 *
8 * (C) Carnegie-Mellon University, 2006-11.
9 * This file is licensed under the Eclipse Public License (EPL)
10 */
11
12 #include "CoinHelperFunctions.hpp"
13
14 #include "OsiRowCut.hpp"
15
16 #include "CouenneCutGenerator.hpp"
17
18 #include "CouenneProblem.hpp"
19 #include "CouenneProblemElem.hpp"
20 #include "CouenneObject.hpp"
21 #include "CouenneBranchingObject.hpp"
22
23 using namespace Ipopt;
24 using namespace Couenne;
25
26 // translate changed bound sparse array into a dense one
27 void sparse2dense (int ncols, t_chg_bounds *chg_bds, int *&changed, int &nchanged);
28
29 /** \brief Constructor.
30 *
31 * Get a variable as an argument and set value_ through a call to
32 * operator () of that exprAux.
33 */
34
CouenneBranchingObject(OsiSolverInterface * solver,const OsiObject * originalObject,JnlstPtr jnlst,CouenneCutGenerator * cutGen,CouenneProblem * problem,expression * var,int way,CouNumber brpoint,bool doFBBT,bool doConvCuts)35 CouenneBranchingObject::CouenneBranchingObject (OsiSolverInterface *solver,
36 const OsiObject * originalObject,
37 JnlstPtr jnlst,
38 CouenneCutGenerator *cutGen,
39 CouenneProblem *problem,
40 expression *var,
41 int way,
42 CouNumber brpoint,
43 bool doFBBT, bool doConvCuts):
44
45 OsiTwoWayBranchingObject (solver, originalObject, way, brpoint),
46
47 cutGen_ (cutGen),
48 problem_ (problem),
49 variable_ (var),
50 jnlst_ (jnlst),
51 doFBBT_ (doFBBT),
52 doConvCuts_ (doConvCuts),
53 downEstimate_ (0.),
54 upEstimate_ (0.),
55 simulate_ (false) {
56
57 // This two-way branching rule is only applied when both lower and
58 // upper bound are finite. Otherwise, a CouenneThreeWayBranchObj is
59 // used (see CouenneThreeWayBranchObj.hpp).
60 //
61 // The rule is as follows:
62 //
63 // - if x is well inside the interval (both bounds are infinite or
64 // there is a difference of at least COU_NEAR_BOUND), set
65 // value_ to x;
66 //
67 // - otherwise, try to get far from bounds by setting value_ to a
68 // convex combination of current and midpoint
69 //
70 // TODO: consider branching value that maximizes distance from
71 // current point (how?)
72
73 CouNumber lb, ub;
74
75 variable_ -> getBounds (lb, ub);
76
77 // bounds may have tightened and may exclude value_ now, update it
78
79 value_ = (fabs (brpoint) < 1e27) ? brpoint : (*variable_) ();
80
81 if (lb < -large_bound)
82 if (ub > large_bound) value_ = 0.; // ]-inf,+inf[
83 else value_ = ((value_ < -COUENNE_EPS) ? (AGGR_MUL * (-1+value_)) : // ]-inf,u]
84 (value_ > COUENNE_EPS) ? 0. : -AGGR_MUL);
85 else
86 if (ub > large_bound) value_ = ((value_ > COUENNE_EPS) ? (AGGR_MUL * (1+value_)) : // [l,+inf[
87 (value_ < -COUENNE_EPS) ? 0. : AGGR_MUL);
88 else { // [l,u]
89 double margin = fabs (ub-lb) * closeToBounds;
90 if (value_ < lb + margin) value_ = lb + margin;
91 else if (value_ > ub - margin) value_ = ub - margin;
92 }
93
94 jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING,
95 "Branch: x%-3d will branch on %g (cur. %g) [%g,%g]; firstBranch_ = %d\n",
96 variable_ -> Index (), value_, (*variable_) (), lb, ub, firstBranch_);
97 }
98
99
100 /** \brief Execute the actions required to branch, as specified by the
101 * current state of the branching object, and advance the
102 * object's state.
103 *
104 * Returns change in guessed objective on next branch
105 */
106
branch(OsiSolverInterface * solver)107 double CouenneBranchingObject::branch (OsiSolverInterface * solver) {
108
109 jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING,
110 "::branch: at %.20e\n", value_);
111
112 // way = 0 if "<=" node,
113 // 1 if ">=" node
114
115 int
116 way = (!branchIndex_) ? firstBranch_ : !firstBranch_,
117 index = variable_ -> Index ();
118
119 bool
120 integer = variable_ -> isInteger (),
121 infeasible = false;
122
123 CouNumber
124 l = solver -> getColLower () [index],
125 u = solver -> getColUpper () [index],
126 brpt = value_;
127
128 if (brpt < l) brpt = l;
129 else if (brpt > u) brpt = u;
130
131 // anticipate test on bounds w.r.t. test on integrality.
132
133 if (l < -large_bound) {
134 if (u <= large_bound) // ]-inf,u]
135 brpt = ((u < -COUENNE_EPS) ? CoinMax (CoinMax (brpt, .5 * (l+u)), AGGR_MUL * (-1. + brpt)) :
136 (u > COUENNE_EPS) ? 0. : -AGGR_MUL);
137 else brpt = 0.;
138 } else
139 if (u > large_bound) // [l,+inf[
140 brpt = ((l > COUENNE_EPS) ? CoinMin (CoinMin (brpt, .5 * (l+u)), AGGR_MUL * ( 1. + brpt)) :
141 (l < -COUENNE_EPS) ? 0. : AGGR_MUL);
142 else { // [l,u] (finite)
143
144 CouNumber point = default_alpha * brpt + (1. - default_alpha) * (l + u) / 2.;
145
146 if ((point-l) / (u-l) < closeToBounds) brpt = l + (u-l) * closeToBounds;
147 else if ((u-point) / (u-l) < closeToBounds) brpt = u + (l-u) * closeToBounds;
148 else brpt = point;
149 }
150
151 // If brpt is integer and the variable is constrained to be integer,
152 // there will be a valid but weak branching. Modify brpt depending
153 // on way and on the bounds on the variable, so that the usual
154 // integer branching will be performed.
155
156 if (integer &&
157 ::isInteger (brpt)) {
158
159 // Look at all possible cases (l,u are bounds, b is the branching
160 // point. l,u,b all integer):
161 //
162 // 1) l < b < u: first branch on b +/- 1 depending on branch
163 // direction, right branch on b;
164 //
165 // 2) l <= b < u: LEFT branch on b, RIGHT branch on b+1
166 //
167 // 3) l < b <= u: LEFT branch on b-1, RIGHT branch on b
168
169 // assert ((brpt - l > .5) ||
170 // (u - brpt > .5));
171
172 if ((brpt - l > .5) &&
173 (u - brpt > .5)) {// brpt is integer interior point of [l,u]
174
175 if (!branchIndex_) { // if this is the first branch operation
176 if (!way) brpt -= (1. - COUENNE_EPS);
177 else brpt += (1. - COUENNE_EPS);
178 }
179 else { // adjust brpt so that interval covers integer value
180 // necessary for nvs13.nl
181 if (!way) brpt += COUENNE_EPS;
182 else brpt -= COUENNE_EPS;
183 }
184 }
185 else {
186 if (u - brpt > .5) {
187 if (way) {
188 brpt += (1. - COUENNE_EPS);
189 }
190 else { // adjust brpt so that interval covers integer value
191 brpt += COUENNE_EPS;
192 }
193 }
194 else {
195 if (brpt - l > .5) {
196 if (!way) {
197 brpt -= (1. - COUENNE_EPS);
198 }
199 else { // adjust brpt so that interval covers integer value
200 brpt -= COUENNE_EPS;
201 }
202 }
203 else { // u == l == brpt; must still branch to fix variables in object
204 // but one branch is infeasible
205 if(way) {
206 solver->setColLower(index, u+1); // infeasible
207 }
208 }
209 }
210 }
211 }
212
213 jnlst_ -> Printf (J_ITERSUMMARY, J_BRANCHING, "Branching: x%-3d %c= %.20e\n",
214 index, way ? '>' : '<', integer ? (way ? ceil (brpt) : floor (brpt)) : brpt);
215
216 if (way) {
217 if (brpt < l) jnlst_->Printf (J_STRONGWARNING, J_BRANCHING, "Nonsense up-br: [ %.8f ,(%.8f)] -> %.8f\n",l,u,value_);
218 else if (brpt < l+COUENNE_EPS) jnlst_->Printf (J_STRONGWARNING, J_BRANCHING, "## WEAK up-br: [ %.8f ,(%.8f)] -> %.8f\n",l,u,value_);
219 } else {
220 if (brpt > u) jnlst_->Printf (J_STRONGWARNING, J_BRANCHING, "Nonsense dn-br: [(%.8f), %.8f ] -> %.8f\n",l,u,value_);
221 else if (brpt > u-COUENNE_EPS) jnlst_->Printf (J_STRONGWARNING, J_BRANCHING, "## WEAK dn-br: [(%.8f), %.8f ] -> %.8f\n",l,u,value_);
222 }
223
224 /*
225 double time = CoinCpuTime ();
226 jnlst_ -> Printf (J_VECTOR, J_BRANCHING,"[vbctool] %02d:%02d:%02d.%02d_I x%d%c=%g_[%g,%g]\n",
227 (int) (floor(time) / 3600),
228 (int) (floor(time) / 60) % 60,
229 (int) floor(time) % 60,
230 (int) ((time - floor (time)) * 100),
231 index, way ? '>' : '<', integer ? ((way ? ceil (brpt): floor (brpt))) : brpt,
232 solver -> getColLower () [index],
233 solver -> getColUpper () [index]);
234 */
235
236 t_chg_bounds *chg_bds = NULL;
237
238 // Apply core branching decision (usually new variable bound)
239
240 branchCore (solver, index, way, integer, brpt, chg_bds);
241
242 //CouenneSolverInterface *couenneSolver = dynamic_cast <CouenneSolverInterface *> (solver);
243 //CouenneProblem *p = cutGen_ -> Problem ();
244 //couenneSolver -> CutGen () -> Problem ();
245
246 int
247 nvars = problem_ -> nVars (),
248 objind = problem_ -> Obj (0) -> Body () -> Index ();
249
250 //CouNumber &estimate = way ? upEstimate_ : downEstimate_;
251 CouNumber estimate = 0.;//way ? upEstimate_ : downEstimate_;
252
253 if ((doFBBT_ && problem_ -> doFBBT ()) ||
254 (doConvCuts_ && simulate_ && cutGen_)) {
255
256 problem_ -> domain () -> push (solver); // have to alloc+copy
257
258 if ( doFBBT_ && // this branching object should do FBBT, and
259 problem_ -> doFBBT ()) { // problem allowed to do FBBT
260
261 problem_ -> installCutOff ();
262
263 if (!problem_ -> btCore (chg_bds)) // done FBBT and this branch is infeasible
264 infeasible = true; // ==> report it
265
266 else {
267
268 const double
269 *lb = solver -> getColLower (),
270 *ub = solver -> getColUpper ();
271
272 //CouNumber newEst = problem_ -> Lb (objind) - lb [objind];
273 estimate = CoinMax (0., objind >= 0 ? (problem_ -> Lb (objind) - lb [objind]) : 0.);
274
275 for (int i=0; i<nvars; i++) {
276 if (problem_ -> Lb (i) > lb [i]) solver -> setColLower (i, problem_ -> Lb (i));
277 if (problem_ -> Ub (i) < ub [i]) solver -> setColUpper (i, problem_ -> Ub (i));
278 }
279 }
280 }
281
282 if (!infeasible && doConvCuts_ && simulate_ && cutGen_) {
283
284 // generate convexification cuts before solving new node's LP
285 int nchanged, *changed = NULL;
286 OsiCuts cs;
287
288 // sparsify structure with info on changed bounds and get convexification cuts
289 sparse2dense (nvars, chg_bds, changed, nchanged);
290 cutGen_ -> genRowCuts (*solver, cs, nchanged, changed, chg_bds);
291 free (changed);
292
293 solver -> applyCuts (cs);
294 }
295
296 //delete [] chg_bds;
297
298 problem_ -> domain () -> pop ();
299 }
300
301 if (chg_bds) delete [] chg_bds;
302
303 // next time do other branching
304 branchIndex_++;
305
306 return (infeasible ? COIN_DBL_MAX : estimate); // estimated change in objective function
307 }
308