1 // (C) Copyright International Business Machines Corporation 2006, 2007
2 // All Rights Reserved.
3 // This code is published under the Eclipse Public License.
4 //
5 // Authors :
6 // Laszlo Ladanyi, International Business Machines Corporation
7 // Pierre Bonami, Carnegie Mellon University
8
9 #include "CoinHelperFunctions.hpp"
10 #include "BCP_lp_node.hpp"
11 #include "BCP_lp.hpp"
12 #include "BCP_lp_functions.hpp"
13 #include "BM.hpp"
14 #include "BonChooseVariable.hpp"
15 #include "BonCurvBranchingSolver.hpp"
16 #include "BonQpBranchingSolver.hpp"
17 #include "BonLpBranchingSolver.hpp"
18 #include "BonOsiTMINLPInterface.hpp"
19 #include "BonIpoptWarmStart.hpp"
20
21 #ifndef BM_DEBUG_PRINT
22 #define BM_DEBUG_PRINT 0
23 #endif
24
25 static bool ifprint = true;
26 static bool ifprint2 = false;
27
28 //#############################################################################
29
30 BCP_branching_decision
select_branching_candidates(const BCP_lp_result & lpres,const BCP_vec<BCP_var * > & vars,const BCP_vec<BCP_cut * > & cuts,const BCP_lp_var_pool & local_var_pool,const BCP_lp_cut_pool & local_cut_pool,BCP_vec<BCP_lp_branching_object * > & cands,bool force_branch)31 BM_lp::select_branching_candidates(const BCP_lp_result& lpres,
32 const BCP_vec<BCP_var*>& vars,
33 const BCP_vec<BCP_cut*>& cuts,
34 const BCP_lp_var_pool& local_var_pool,
35 const BCP_lp_cut_pool& local_cut_pool,
36 BCP_vec<BCP_lp_branching_object*>& cands,
37 bool force_branch)
38 {
39 /* FIXME: this is good only for BB and using NLP to solve a node */
40 bm_stats.incNumberNodeSolves();
41
42 Bonmin::OsiTMINLPInterface* nlp =
43 dynamic_cast<Bonmin::OsiTMINLPInterface*>(getLpProblemPointer()->lp_solver);
44 // If we are doing pure B&B then we have an nlp, and then we check for
45 // consecutive failures.
46 if (nlp) {
47 if (lpres.termcode() & BCP_Abandoned) {
48 if (nlp->isIterationLimitReached()) {
49 print(ifprint, "\
50 BM_lp: At node %i : WARNING: nlp reached iter limit. Will force branching\n",
51 current_index());
52 } else {
53 print(ifprint, "\
54 BM_lp: At node %i : WARNING: nlp is abandoned. Will force branching\n",
55 current_index());
56 }
57 // nlp failed
58 nlp->forceBranchable();
59 numNlpFailed_++;
60 if (numNlpFailed_ >= par.entry(BM_par::NumNlpFailureMax)) {
61 print(ifprint, "WARNING! Too many (%i) NLP failures in a row. Abandoning node.",
62 numNlpFailed_);
63 return BCP_DoNotBranch_Fathomed;
64 }
65 } else {
66 numNlpFailed_ = 0;
67 }
68 }
69
70 OsiBranchingInformation brInfo(nlp, false, true);
71 brInfo.cutoff_ = upper_bound() + get_param(BCP_lp_par::Granularity);
72 brInfo.objectiveValue_ = lpres.objval();
73 brInfo.integerTolerance_ = integerTolerance_;
74 brInfo.timeRemaining_ = get_param(BCP_lp_par::MaxRunTime) - CoinCpuTime();
75 brInfo.numberSolutions_ = 0; /*FIXME*/
76 brInfo.numberBranchingSolutions_ = 0; /*FIXME numBranchingSolutions_;*/
77 brInfo.depth_ = current_level();
78
79 BCP_branching_decision brDecision;
80 if (bonmin_.getAlgorithm() == 0) {
81 /* if pure B&B */
82 brDecision = bbBranch(brInfo, cands);
83 } else {
84 brDecision = hybridBranch();
85 }
86
87 return brDecision;
88 }
89
90 //-----------------------------------------------------------------------------
91
92 BCP_branching_decision
hybridBranch()93 BM_lp::hybridBranch()
94 {
95 // FIXME: most of the pureBB stuff should work here.
96 throw BCP_fatal_error("BM_lp: FIXME: make hybrid work...");
97 }
98
99 /*****************************************************************************/
100
101 void
unpack_pseudo_costs(BCP_buffer & buf)102 BM_lp::unpack_pseudo_costs(BCP_buffer& buf)
103 {
104 Bonmin::BonChooseVariable* choose =
105 dynamic_cast<Bonmin::BonChooseVariable*>(bonmin_.branchingMethod());
106 OsiPseudoCosts& pseudoCosts = choose->pseudoCosts();
107 int numObj = pseudoCosts.numberObjects();
108 double* upTotalChange = pseudoCosts.upTotalChange();
109 int* upNumber = pseudoCosts.upNumber();
110 double* downTotalChange = pseudoCosts.downTotalChange();
111 int* downNumber = pseudoCosts.downNumber();
112
113 buf.unpack(upTotalChange, numObj, false);
114 buf.unpack(upNumber, numObj, false);
115 buf.unpack(downTotalChange, numObj, false);
116 buf.unpack(downNumber, numObj, false);
117 }
118
119
120 //-----------------------------------------------------------------------------
121
122 void
send_pseudo_cost_update(OsiBranchingInformation & branchInfo)123 BM_lp::send_pseudo_cost_update(OsiBranchingInformation& branchInfo)
124 {
125 bm_buf.clear();
126 int itmp;
127 double objchange;
128 itmp = BM_PseudoCostUpdate;
129 bm_buf.pack(itmp);
130 for (int i = 0; i < objNum_; ++i) {
131 const BM_SB_result& sbres = sbResult_[i];
132 if ((sbres.branchEval & 1) != 0 && sbres.status[0] != BCP_Abandoned) {
133 bm_buf.pack(sbres.objInd);
134 itmp = 0;
135 bm_buf.pack(itmp);
136 if (sbres.status[0] == BCP_ProvenOptimal) {
137 objchange = sbres.objval[0] - branchInfo.objectiveValue_;
138 } else { // Must be BCP_ProvenPrimalInf
139 if (branchInfo.cutoff_ < 1e50) {
140 objchange = 2*(branchInfo.cutoff_-branchInfo.objectiveValue_);
141 } else {
142 objchange = 2*fabs(branchInfo.objectiveValue_);
143 }
144 }
145 bm_buf.pack(objchange/sbres.varChange[0]);
146 }
147 if ((sbres.branchEval & 2) != 0 && sbres.status[1] != BCP_Abandoned) {
148 bm_buf.pack(sbres.objInd);
149 itmp = 1;
150 bm_buf.pack(itmp);
151 if (sbres.status[1] == BCP_ProvenOptimal) {
152 objchange = sbres.objval[1] - branchInfo.objectiveValue_;
153 } else { // Must be BCP_ProvenPrimalInf
154 if (branchInfo.cutoff_ < 1e50) {
155 objchange = 2*(branchInfo.cutoff_-branchInfo.objectiveValue_);
156 } else {
157 objchange = 2*fabs(branchInfo.objectiveValue_);
158 }
159 }
160 bm_buf.pack(objchange/sbres.varChange[1]);
161 }
162 }
163 itmp = -1;
164 bm_buf.pack(itmp);
165 send_message(parent(), bm_buf);
166 }
167
168 //-----------------------------------------------------------------------------
169
170 int
sort_objects(OsiBranchingInformation & branchInfo,Bonmin::BonChooseVariable * choose,int & branchNum)171 BM_lp::sort_objects(OsiBranchingInformation& branchInfo,
172 Bonmin::BonChooseVariable* choose, int& branchNum)
173 {
174 const OsiObject* const * objects = branchInfo.solver_->objects();
175 double upMult, downMult;
176 choose->computeMultipliers(upMult, downMult);
177 const double MAXMIN = choose->maxminCrit(&branchInfo);
178
179 /* Order all objects that can be branched on */
180 int lastPriority = objects[objInd_[0]]->priority();
181 int infBlockStart = 0;
182 int feasBlockStart = 0;
183 branchNum = 0;
184 infNum_ = 0;
185 feasNum_ = 0;
186
187 const bool isRoot = (current_index() == 0);
188 int way;
189
190 const bool disregardPriorities = par.entry(BM_par::DisregardPriorities);
191 const bool usePseudoCosts = par.entry(BM_par::UsePseudoCosts);
192
193 for (int i = 0; i < objNum_; ++i) {
194 const int ind = objInd_[i];
195 const OsiObject* object = objects[ind];
196 double value = object->infeasibility(&branchInfo, way);
197 if (value > 0.0) {
198 if (value >= 1e50) { // infeasible
199 return -1;
200 }
201 if (! disregardPriorities) {
202 int priorityLevel = object->priority();
203 if (lastPriority < priorityLevel) {
204 // sort the entries based on their usefulness
205 if (infBlockStart < infNum_) {
206 if (par.entry(BM_par::DecreasingSortInSetupList)) {
207 CoinSort_2(infUseful_ + infBlockStart, infUseful_ + infNum_,
208 infInd_ + infBlockStart,
209 CoinFirstGreater_2<double,int>());
210 } else {
211 CoinSort_2(infUseful_ + infBlockStart, infUseful_ + infNum_,
212 infInd_ + infBlockStart);
213 }
214 infBlockStart = infNum_;
215 }
216 lastPriority = priorityLevel;
217 }
218 }
219 double dummy;
220 infInd_[infNum_] = ind;
221 if (usePseudoCosts) {
222 infUseful_[infNum_] = isRoot ?
223 value : choose->computeUsefulness(MAXMIN, upMult, downMult, value,
224 object, ind, dummy);
225 } else {
226 infUseful_[infNum_] = value;
227 }
228
229 ++infNum_;
230 branchNum += 2;
231
232 } else { /* value == 0.0 */
233 const OsiSOS* sos = dynamic_cast<const OsiSOS*>(object);
234 if (sos) {
235 // if an sos is feasible thne we don't do strong branching on it
236 continue;
237 }
238 const int iCol = object->columnNumber();
239 const double lb = branchInfo.lower_[iCol];
240 const double ub = branchInfo.upper_[iCol];
241 if (fabs(ub - lb) < 0.1) {
242 continue;
243 }
244 value = branchInfo.solution_[iCol];
245 ++branchNum;
246 if (floor(value+0.5) > lb && ceil(value-0.5) < ub) {
247 // The variable is integer, but neither at its lower nor at its upper
248 // bound (the test accounts for tolerances)
249 ++branchNum;
250 }
251 if (! disregardPriorities) {
252 int priorityLevel = object->priority();
253 if (lastPriority < priorityLevel) {
254 // sort the entries based on their usefulness
255 if (feasBlockStart < feasNum_) {
256 if (par.entry(BM_par::DecreasingSortInSetupList)) {
257 CoinSort_2(feasUseful_ + feasBlockStart, feasUseful_ + feasNum_,
258 feasInd_ + feasBlockStart,
259 CoinFirstGreater_2<double,int>());
260 } else {
261 CoinSort_2(feasUseful_ + feasBlockStart, feasUseful_ + feasNum_,
262 feasInd_ + feasBlockStart);
263 }
264 }
265 lastPriority = priorityLevel;
266 }
267 }
268 double dummy;
269 feasInd_[feasNum_] = ind;
270 feasUseful_[feasNum_] = choose->computeUsefulness(MAXMIN,
271 upMult, downMult, value,
272 object, ind, dummy);
273 ++feasNum_;
274 }
275 }
276 if (infBlockStart < infNum_) {
277 if (par.entry(BM_par::DecreasingSortInSetupList)) {
278 CoinSort_2(infUseful_ + infBlockStart, infUseful_ + infNum_,
279 infInd_ + infBlockStart,
280 CoinFirstGreater_2<double,int>());
281 } else {
282 CoinSort_2(infUseful_ + infBlockStart, infUseful_ + infNum_,
283 infInd_ + infBlockStart);
284 }
285 }
286 if (feasBlockStart < feasNum_) {
287 if (par.entry(BM_par::DecreasingSortInSetupList)) {
288 CoinSort_2(feasUseful_ + feasBlockStart, feasUseful_ + feasNum_,
289 feasInd_ + feasBlockStart,
290 CoinFirstGreater_2<double,int>());
291 } else {
292 CoinSort_2(feasUseful_ + feasBlockStart, feasUseful_ + feasNum_,
293 feasInd_ + feasBlockStart);
294 }
295 }
296
297 #if (BM_DEBUG_PRINT != 0)
298 const double t = CoinWallclockTime();
299 printf("LP %.3f: node: %i depth: %i obj: %f infNum: %i feasNum: %i soltime: %.3f\n",
300 t-start_time(), current_index(), current_level(),
301 branchInfo.objectiveValue_,
302 infNum_, feasNum_, t-node_start_time);
303 node_start_time = t;
304 #endif
305
306 return infNum_;
307 }
308
309 //-----------------------------------------------------------------------------
310
311 void
clear_SB_results()312 BM_lp::clear_SB_results()
313 {
314 for (int i = 0; i < objNum_; ++i) {
315 sbResult_[i].branchEval = 0;
316 }
317 bestSbResult_ = NULL;
318 }
319
320 //-----------------------------------------------------------------------------
321
322 void
collect_branch_data(OsiBranchingInformation & branchInfo,OsiSolverInterface * solver,const int branchNum,BM_BranchData * branchData)323 BM_lp::collect_branch_data(OsiBranchingInformation& branchInfo,
324 OsiSolverInterface* solver,
325 const int branchNum,
326 BM_BranchData* branchData)
327 {
328 int i;
329 int b = 0;
330 for (i = 0; i < infNum_; ++i) {
331 /* FIXME: think about SOS */
332 const int objInd = infInd_[i];
333 const int colInd = solver->object(objInd)->columnNumber();
334 const double val = branchInfo.solution_[colInd];
335 branchData[b].changeType = BM_Var_DownBranch;
336 branchData[b].objInd = objInd;
337 branchData[b].colInd = colInd;
338 branchData[b].solval = val;
339 branchData[b].bd = floor(val);
340 BM_SB_result& sbres = sbResult_[objInd];
341 sbres.objInd = objInd;
342 sbres.varChange[0] = val - floor(val);
343 ++b;
344 if (b == branchNum) {
345 return;
346 }
347 branchData[b].changeType = BM_Var_UpBranch;
348 branchData[b].objInd = objInd;
349 branchData[b].colInd = colInd;
350 branchData[b].solval = val;
351 branchData[b].bd = ceil(val);
352 sbres.varChange[1] = ceil(val) - val;
353 ++b;
354 if (b == branchNum) {
355 return;
356 }
357 }
358 for (i = 0; i < feasNum_; ++i) {
359 const int objInd = feasInd_[i];
360 const int colInd = solver->object(objInd)->columnNumber();
361 const double val = branchInfo.solution_[colInd];
362 const double lb = branchInfo.lower_[colInd];
363 const double ub = branchInfo.upper_[colInd];
364 if (floor(val+0.5) > lb) { // not at its lb
365 branchData[b].changeType = BM_Var_DownBranch;
366 branchData[b].objInd = objInd;
367 branchData[b].colInd = colInd;
368 branchData[b].solval = val;
369 branchData[b].bd = floor(val - 0.5);
370 ++b;
371 if (b == branchNum) {
372 return;
373 }
374 }
375 if (ceil(val-0.5) < ub) { // not at its ub
376 branchData[b].changeType = BM_Var_UpBranch;
377 branchData[b].objInd = objInd;
378 branchData[b].colInd = colInd;
379 branchData[b].solval = val;
380 branchData[b].bd = ceil(val + 0.5);
381 ++b;
382 if (b == branchNum) {
383 return;
384 }
385 }
386 }
387 }
388
389 //-----------------------------------------------------------------------------
390
391 void
BM_solve_branches(OsiSolverInterface * solver,const CoinWarmStart * cws,const int numBranch,BM_BranchData * bD)392 BM_solve_branches(OsiSolverInterface* solver, const CoinWarmStart* cws,
393 const int numBranch, BM_BranchData* bD)
394 {
395 for (int i = 0; i < numBranch; ++i) {
396 double t = CoinWallclockTime();
397 const int ind = bD[i].colInd;
398 const int field = bD[i].changeType == BM_Var_UpBranch ? 1 : 0;
399 const double old_lb = solver->getColLower()[ind];
400 const double old_ub = solver->getColUpper()[ind];
401 if (field == 0) {
402 solver->setColUpper(ind, bD[i].bd);
403 } else {
404 solver->setColLower(ind, bD[i].bd);
405 }
406 if (cws) {
407 solver->setWarmStart(cws);
408 solver->resolve();
409 } else {
410 solver->initialSolve();
411 }
412 bD[i].status =
413 (solver->isAbandoned() ? BCP_Abandoned : 0) |
414 (solver->isProvenOptimal() ? BCP_ProvenOptimal : 0) |
415 (solver->isProvenPrimalInfeasible() ? BCP_ProvenPrimalInf : 0);
416 bD[i].objval =
417 (bD[i].status & BCP_ProvenOptimal) != 0 ? solver->getObjValue() : 0.0;
418 bD[i].iter = solver->getIterationCount();
419 solver->setColBounds(ind, old_lb, old_ub);
420 bD[i].time = CoinWallclockTime() - t;
421 }
422 }
423
424 //-----------------------------------------------------------------------------
425
426 void
BM_register_branch_results(const int numBranch,const BM_BranchData * bD,BM_SB_result * sbResults)427 BM_register_branch_results(const int numBranch, const BM_BranchData* bD,
428 BM_SB_result* sbResults)
429 {
430 for (int i = 0; i < numBranch; ++i) {
431 const int field = bD[i].changeType == BM_Var_UpBranch ? 1 : 0;
432 BM_SB_result& sbres = sbResults[bD[i].objInd];
433 sbres.objInd = bD[i].objInd;
434 sbres.branchEval |= field == 0 ? 1 : 2;
435 sbres.status[field] = bD[i].status;
436 sbres.objval[field] = bD[i].objval;
437 sbres.iter[field] = bD[i].iter;
438 sbres.varChange[field] = fabs(bD[i].solval - bD[i].bd);
439 }
440 }
441
442 //-----------------------------------------------------------------------------
443
444 void
do_distributed_SB(OsiBranchingInformation & branchInfo,OsiSolverInterface * solver,const CoinWarmStart * cws,const int branchNum,const int * pids,const int pidNum)445 BM_lp::do_distributed_SB(OsiBranchingInformation& branchInfo,
446 OsiSolverInterface* solver,
447 const CoinWarmStart* cws,
448 const int branchNum,
449 const int* pids, const int pidNum)
450 {
451 const double * clb = solver->getColLower();
452 const double * cub = solver->getColUpper();
453 const int numCols = solver->getNumCols();
454 bm_buf.clear();
455 int tag = BM_StrongBranchRequest;
456 bm_buf.pack(tag);
457 bm_buf.pack(clb, numCols);
458 bm_buf.pack(cub, numCols);
459 bm_buf.pack(branchInfo.objectiveValue_);
460 bm_buf.pack(branchInfo.cutoff_);
461 bool has_ws = cws != NULL;
462 bm_buf.pack(has_ws);
463 if (has_ws) {
464 BCP_lp_prob* p = getLpProblemPointer();
465 CoinWarmStart* cws_tmp = cws->clone(); // the next call destroys it
466 BCP_warmstart* ws = cws ? BCP_lp_convert_CoinWarmStart(*p, cws_tmp) : NULL;
467 BCP_pack_warmstart(ws, bm_buf);
468 delete ws;
469 }
470
471 const int fixed_size = bm_buf.size();
472
473 // collect what we'll need to send off data
474 BM_BranchData* branchData = new BM_BranchData[branchNum];
475 collect_branch_data(branchInfo, solver, branchNum, branchData);
476
477 // We have branchNum branches to process on pidNum+1 (the last is the local
478 // process) processes.
479 int branchLeft = branchNum;
480 BM_BranchData* bD = branchData;
481 for (int pidLeft = pidNum; pidLeft > 0; --pidLeft) {
482 int numSend = branchLeft / pidLeft;
483 if (numSend * pidLeft < branchLeft) {
484 ++numSend;
485 }
486 bm_buf.set_size(fixed_size);
487 // Now pack where we are in branchData and pack numSend branches
488 int location = bD - branchData;
489 bm_buf.pack(location);
490 bm_buf.pack(numSend);
491 for (int s = 0; s < numSend; ++s) {
492 bm_buf.pack(bD[s].changeType);
493 bm_buf.pack(bD[s].objInd);
494 bm_buf.pack(bD[s].colInd);
495 bm_buf.pack(bD[s].solval);
496 bm_buf.pack(bD[s].bd);
497 }
498 send_message(pids[pidLeft-1], bm_buf);
499 bD += numSend;
500 branchLeft -= numSend;
501 }
502 assert(branchNum/(pidNum+1) == branchLeft);
503
504 // Process the leftover branches locally
505 /* FIXME: this assumes that the solver is the NLP solver. Maybe we should use
506 the nlp solver in BM_lp */
507 BM_solve_branches(solver, cws, branchLeft, bD);
508 bm_stats.incNumberSbSolves(branchLeft);
509 BM_register_branch_results(branchLeft, bD, sbResult_);
510
511 // Receive the results from the other processes
512 int numResults = branchLeft;
513 while (numResults < branchNum) {
514 bm_buf.clear();
515 receive_message(BCP_AnyProcess, bm_buf, BCP_Msg_User);
516 bm_buf.unpack(tag);
517 assert(tag == BM_StrongBranchResult);
518 int location;
519 int numRes;
520 bm_buf.unpack(location);
521 bm_buf.unpack(numRes);
522 BM_BranchData* bD = branchData + location;
523 for (int i = 0; i < numRes; ++i) {
524 bm_buf.unpack(bD[i].status);
525 bm_buf.unpack(bD[i].objval);
526 bm_buf.unpack(bD[i].iter);
527 bm_buf.unpack(bD[i].time);
528 }
529 BM_register_branch_results(numRes, bD, sbResult_);
530 numResults += numRes;
531 }
532 }
533
534 //-----------------------------------------------------------------------------
535
536 bool
isBranchFathomable(int status,double obj)537 BM_lp::isBranchFathomable(int status, double obj)
538 {
539 return ( (status & BCP_ProvenPrimalInf) ||
540 ((status & BCP_ProvenOptimal) && over_ub(obj)) );
541 }
542
543 //-----------------------------------------------------------------------------
544
545 int
process_SB_results(OsiBranchingInformation & branchInfo,OsiSolverInterface * solver,Bonmin::BonChooseVariable * choose,OsiBranchingObject * & branchObject)546 BM_lp::process_SB_results(OsiBranchingInformation& branchInfo,
547 OsiSolverInterface* solver,
548 Bonmin::BonChooseVariable* choose,
549 OsiBranchingObject*& branchObject)
550 {
551 int evaluated = 0;
552 int i;
553 #if defined(DEBUG_PRINT)
554 const double t = CoinWallclockTime()-start_time();
555 #endif
556
557 #if 0 // defined(DEBUG_PRINT)
558 const double t = CoinWallclockTime()-start_time();
559 for (i = 0; i < infNum_; ++i) {
560 const BM_SB_result& sbres = sbResult_[infInd_[i]];
561 if (sbres.branchEval == 0) {
562 continue;
563 }
564 printf("LP %.3f: SB: node: %i inf col: %i stati: %i %i, obj: %f %f time: %.3f %.3f\n",
565 t, current_index(), sbres.colInd,
566 sbres.status[0], sbres.status[1],
567 sbres.objval[0], sbres.objval[1], sbres.time[0], sbres.time[1]);
568 }
569 for (i = 0; i < feasNum_; ++i) {
570 const BM_SB_result& sbres = sbResult_[feasInd_[i]];
571 if (sbres.branchEval == 0) {
572 continue;
573 }
574 printf("LP %.3f: SB: node: %i feas col: %i stati: %i %i, obj: %f %f time: %.3f %.3f\n",
575 t, current_index(), sbres.colInd,
576 sbres.status[0], sbres.status[1],
577 sbres.objval[0], sbres.objval[1], sbres.time[0], sbres.time[1]);
578 }
579 #endif
580 // First check if we can fathom the node
581 int listLen=0; // we want this for the bm_stats
582 for (i = 0; i < infNum_; ++i) {
583 const BM_SB_result& sbres = sbResult_[infInd_[i]];
584 if (sbres.branchEval == 0) {
585 continue;
586 }
587 assert(sbres.branchEval == 3);
588 if (isBranchFathomable(sbres.status[0], sbres.objval[0]) &&
589 isBranchFathomable(sbres.status[1], sbres.objval[1])) {
590 #if (BM_DEBUG_PRINT != 0)
591 const double wallclock = CoinWallclockTime();
592 #if 0
593 printf("LP %.3f: SBres: node: %i FATHOM inf/eval/cand: time: %.3f\n",
594 wallclock-start_time(), current_index(),
595 infNum_,
596 wallclock-node_start_time);
597 #else
598 printf("LP %.3f: SBres: node: %i FATHOM time: %.3f\n",
599 wallclock-start_time(), current_index(),
600 wallclock-node_start_time);
601 #endif
602 #endif
603 return -2;
604 }
605 ++listLen;
606 }
607
608 // Nope. Let's see if we can fix anything.
609 // 0: nothing fixed bit 0: fixed but stayd feas bit 1: fixed and lost feas
610 int fixedStat = 0;
611 for (i = 0; i < infNum_; ++i) {
612 const BM_SB_result& sbres = sbResult_[infInd_[i]];
613 if (sbres.branchEval == 0) {
614 continue;
615 }
616 ++evaluated;
617 if (isBranchFathomable(sbres.status[0], sbres.objval[0])) {
618 const int colInd = sbres.colInd;
619 solver->setColLower(colInd, ceil(branchInfo.solution_[colInd]));
620 fixedStat |= 2;
621 bm_stats.incNumberFixed();
622 }
623 if (isBranchFathomable(sbres.status[1], sbres.objval[1])) {
624 const int colInd = sbres.colInd;
625 solver->setColUpper(colInd, floor(branchInfo.solution_[colInd]));
626 fixedStat |= 2;
627 bm_stats.incNumberFixed();
628 }
629 }
630 for (i = 0; i < feasNum_; ++i) {
631 const BM_SB_result& sbres = sbResult_[feasInd_[i]];
632 if (sbres.branchEval == 0) {
633 continue;
634 }
635 ++evaluated;
636 if ( (sbres.branchEval & 1) &&
637 isBranchFathomable(sbres.status[0], sbres.objval[0]) ) {
638 // Down branch evaluated and fathomable
639 const int colInd = sbres.colInd;
640 solver->setColLower(colInd, ceil(branchInfo.solution_[colInd] - 0.5));
641 fixedStat |= 1;
642 bm_stats.incNumberFixed();
643 }
644 if ( (sbres.branchEval & 2) &&
645 isBranchFathomable(sbres.status[1], sbres.objval[1]) ) {
646 // Up branch evaluated and fathomable
647 const int colInd = sbres.colInd;
648 solver->setColUpper(colInd, floor(branchInfo.solution_[colInd] + 0.5));
649 fixedStat |= 1;
650 bm_stats.incNumberFixed();
651 }
652 }
653 if ((fixedStat & 2) != 0) {
654 #if (BM_DEBUG_PRINT != 0)
655 printf("LP %.3f: SBres: node: %i RESOLVE time: %.3f\n",
656 t - start_time(), current_index(), t-node_start_time);
657 node_start_time = t;
658 #endif
659 return -1;
660 }
661
662 // OK, we may have fixed some, but we got to branch, nevertheless.
663 // Try to choose something where both sides have finished properly.
664 // Note: at this point all stati of the inf pieces is either abandoned,
665 // optimal (if primal infeas then we fixed and lost feasibility)
666 const bool preferHigh = par.entry(BM_par::PreferHighCombinationInBranching);
667 int best = -1;
668 int bestWhichWay = 1;
669 double bestTrusted = preferHigh ? -COIN_DBL_MAX : COIN_DBL_MAX;
670 const double MAXMIN = choose->maxminCrit(&branchInfo);
671
672 for (i = 0; i < infNum_; ++i) {
673 const int objInd = infInd_[i];
674 const BM_SB_result& sbres = sbResult_[objInd];
675 if (sbres.branchEval == 0) {
676 // FIXME: Check the pseudocost
677 continue;
678 }
679 if ((sbres.status[0]&BCP_Abandoned) || (sbres.status[1]&BCP_Abandoned)){
680 continue;
681 }
682 double upEst = sbres.objval[1] - branchInfo.objectiveValue_;
683 double downEst = sbres.objval[0] - branchInfo.objectiveValue_;
684 double value =
685 MAXMIN*CoinMin(upEst,downEst) + (1.0-MAXMIN)*CoinMax(upEst,downEst);
686 const bool better = ( ( preferHigh && (value > bestTrusted)) ||
687 ( !preferHigh && (value < bestTrusted)) );
688 if (better) {
689 bestTrusted = value;
690 best = i;
691 bestWhichWay = upEst > downEst ? 0 : 1;
692 // override if there is a preferred way
693 const OsiObject* object = solver->object(objInd);
694 if (object->preferredWay() >= 0 && object->infeasibility()) {
695 bestWhichWay = object->preferredWay();
696 }
697 }
698 }
699 if (best == -1) {
700 // OMG... In *ALL* evaluated candidates at least one side was abandoned...
701 // Oh well, pick the first one not evaluated, or if no such thing then the
702 // first where only one side failed, or if no such thing then anything...
703 // FIXME: this loop will not be needed when pseudocosts are used
704 for (i = 0; i < infNum_; ++i) {
705 const BM_SB_result& sbres = sbResult_[infInd_[i]];
706 if (sbres.branchEval == 0) {
707 best = i;
708 break;
709 }
710 }
711 if (best == -1) {
712 for (i = 0; i < infNum_; ++i) {
713 const BM_SB_result& sbres = sbResult_[infInd_[i]];
714 if ((sbres.status[0] & BCP_Abandoned) == 0 ||
715 (sbres.status[1] & BCP_Abandoned) == 0) {
716 // prefer something where at least one side was not abandoned...
717 best = i;
718 break;
719 }
720 best = i;
721 }
722 }
723 }
724
725 bm_stats.updateStrongBrachingInfo(best, listLen);
726
727 // At this point best is not -1, create the branching object
728 const OsiObject * object = solver->object(infInd_[best]);
729 branchObject = object->createBranch(solver, &branchInfo, bestWhichWay);
730 bestSbResult_ = sbResult_ + infInd_[best];
731 #if (BM_DEBUG_PRINT != 0)
732 const int ind = object->columnNumber();
733 printf("LP %.3f: SBres: node: %i depth: %i BRANCH time: %.3f evaluated: %i bvar: %i val: %f obj0: %f obj1: %f way: %i\n",
734 t-start_time(), current_index(), current_level(), t-node_start_time,
735 evaluated, ind, branchInfo.solution_[ind],
736 bestSbResult_->objval[0], bestSbResult_->objval[1], bestWhichWay);
737 node_start_time = t;
738 #endif
739
740 return (fixedStat & 1) != 0 ? -1 : 0;
741 }
742
743 //-----------------------------------------------------------------------------
744
745 int
try_to_branch(OsiBranchingInformation & branchInfo,OsiSolverInterface * solver,Bonmin::BonChooseVariable * choose,OsiBranchingObject * & branchObject,bool allowVarFix)746 BM_lp::try_to_branch(OsiBranchingInformation& branchInfo,
747 OsiSolverInterface* solver,
748 Bonmin::BonChooseVariable* choose,
749 OsiBranchingObject*& branchObject,
750 bool allowVarFix)
751 {
752 int returnStatus = 0;
753
754 int branchNum;
755 sort_objects(branchInfo, choose, branchNum);
756
757 if (infNum_ == 0) {
758 return 0;
759 }
760
761 const bool do_distributed_branching = true;
762
763 if (do_distributed_branching) {
764
765 bm_buf.clear();
766 bm_buf.pack(branchNum-1);
767 send_message(parent(), bm_buf, BCP_Msg_RequestProcessList);
768 bm_buf.clear();
769 receive_message(parent(), bm_buf, BCP_Msg_ProcessList);
770 int* pids = NULL;
771 int pidNum;
772 bm_buf.unpack(pids, pidNum);
773
774 clear_SB_results();
775
776 // Get the warmstart information
777 CoinWarmStart* cws = solver->getWarmStart();
778 Bonmin::IpoptWarmStart* iws = dynamic_cast<Bonmin::IpoptWarmStart*>(cws);
779 if (iws && iws->empty()) {
780 delete cws;
781 cws = NULL;
782 }
783
784 const int nodeIndex = current_index();
785 if (nodeIndex == 0) {
786 // we are in the root. We want to process all possible branches upto the
787 // parameter value.
788 branchNum = CoinMin(branchNum, par.entry(BM_par::SBNumBranchesInRoot));
789 } else {
790 if (nodeIndex < par.entry(BM_par::SBMaxLevel)) {
791 branchNum = CoinMin(branchNum,
792 CoinMax(pidNum + 1,
793 par.entry(BM_par::SBNumBranchesInTree)));
794 } else {
795 branchNum = pidNum > 0 ? pidNum + 1 : 0;
796 }
797 }
798
799 if (branchNum > 0) {
800 do_distributed_SB(branchInfo, solver, cws, branchNum, pids, pidNum);
801 returnStatus = process_SB_results(branchInfo, solver, choose,
802 branchObject);
803 send_pseudo_cost_update(branchInfo);
804 } else {
805 // We are too deep in the tree, just do something locally (like
806 // pseudocosts...)
807 returnStatus = BCP_lp_user::try_to_branch(branchInfo, solver, choose,
808 branchObject, allowVarFix);
809 }
810
811 delete[] pids;
812 delete cws;
813
814 } else { /* ! do_distributed_branching ===> Do something locally */
815
816 }
817
818 return returnStatus;
819 }
820
821 //-----------------------------------------------------------------------------
822
823 BCP_branching_decision
bbBranch(OsiBranchingInformation & brInfo,BCP_vec<BCP_lp_branching_object * > & cands)824 BM_lp::bbBranch(OsiBranchingInformation& brInfo,
825 BCP_vec<BCP_lp_branching_object*>& cands)
826 {
827 BCP_lp_prob* p = getLpProblemPointer();
828 OsiSolverInterface* osi = p->lp_solver;
829 Bonmin::OsiTMINLPInterface* nlp =
830 dynamic_cast<Bonmin::OsiTMINLPInterface*>(osi);
831 assert(nlp);
832
833 nlp->getDblParam(OsiPrimalTolerance, brInfo.integerTolerance_);
834
835 BCP_branching_decision retCode;
836 OsiBranchingObject* brObj = NULL;
837
838 const int numCols = nlp->getNumCols();
839 double* clb_old = new double[numCols];
840 double* cub_old = new double[numCols];
841 CoinDisjointCopyN(nlp->getColLower(), numCols, clb_old);
842 CoinDisjointCopyN(nlp->getColUpper(), numCols, cub_old);
843
844 Ipopt::SmartPtr<Ipopt::OptionsList> options = bonmin_.options();
845 int numSB = 0;
846 //const bool sbIsSet =
847 options->GetIntegerValue("number_strong_branch",numSB,"bonmin.");
848 int numSBroot = 0;
849 const bool sbRootIsSet =
850 options->GetIntegerValue("number_strong_branch_root",numSBroot,"bonmin.");
851
852 if (sbRootIsSet && brInfo.depth_ == 0) {
853 bonmin_.branchingMethod()->setNumberStrong(numSBroot);
854 } else {
855 bonmin_.branchingMethod()->setNumberStrong(numSB);
856 }
857
858 Bonmin::BonChooseVariable* choose =
859 dynamic_cast<Bonmin::BonChooseVariable*>(bonmin_.branchingMethod());
860 const int brResult = BM_lp::try_to_branch(brInfo, nlp, choose, brObj, true);
861
862 #if 0
863 if (choose->goodSolution()) {
864 /* yipee! a feasible solution! Check that it's really */
865 const double* sol = choose->goodSolution();
866 BM_solution* bmsol = new BM_solution;
867 //Just copy the solution stored in solver to sol
868 double ptol = integerTolerance_;
869 for (int i = 0 ; i < numCols ; i++) {
870 if (fabs(sol[i]) > ptol)
871 bmsol->add_entry(i, sol[i]);
872 }
873 bmsol->setObjective(choose->goodObjectiveValue());
874 choose->clearGoodSolution();
875 send_feasible_solution(bmsol);
876 delete bmsol;
877 }
878 #endif
879
880 switch (brResult) {
881 case -2:
882 // when doing strong branching a candidate has proved that the
883 // problem is infeasible
884 retCode = BCP_DoNotBranch_Fathomed;
885 break;
886 case -1:
887 // OsiChooseVariable::chooseVariable() returned 2, 3, or 4
888 if (!brObj) {
889 // just go back and resolve
890 retCode = BCP_DoNotBranch;
891 } else {
892 // otherwise might as well branch. The forced variable is
893 // unlikely to jump up one more (though who knows...)
894 retCode = BCP_DoBranch;
895 }
896 break;
897 case 0:
898 if (!brObj) {
899 // nothing got fixed, yet couldn't find something to branch on
900 throw BCP_fatal_error("BM: Couldn't branch!\n");
901 }
902 // we've got a branching object
903 retCode = BCP_DoBranch;
904 break;
905 default:
906 throw BCP_fatal_error("\
907 BM: BCP_lp_user::try_to_branch returned with unknown return code.\n");
908 }
909
910 if (brResult < 0) {
911 // If there were some fixings (brResult < 0) then propagate them
912 // where needed
913
914 // FIXME: This is not nice. Meddling w/ BCP internal data. The BCP
915 // user interface should provide a way to change bounds regardless
916 // whether branching is asked for or not.
917 const double* clb = nlp->getColLower();
918 const double* cub = nlp->getColUpper();
919
920 BCP_vec<BCP_var*>& vars = getLpProblemPointer()->node->vars;
921 for (int i = 0; i < numCols; ++i) {
922 if (clb_old[i] != clb[i] || cub_old[i] != cub[i]) {
923 vars[i]->change_bounds(clb[i], cub[i]);
924 nlp->setColBounds(i, clb[i], cub[i]);
925 }
926 }
927 }
928
929 if (retCode == BCP_DoBranch) {
930 // all possibilities are 2-way branches
931 int order[2] = {0, 1};
932 // Now interpret the result (at this point we must have a brObj
933 OsiIntegerBranchingObject* intBrObj =
934 dynamic_cast<OsiIntegerBranchingObject*>(brObj);
935 if (intBrObj) {
936 if (intBrObj->firstBranch() == 1) {
937 order[0] = 1;
938 order[1] = 0;
939 }
940 BCP_lp_integer_branching_object o(intBrObj);
941 cands.push_back(new BCP_lp_branching_object(o, order));
942 if (bestSbResult_) {
943 BCP_vec<double> lb(2, 0.0);
944 lb[0] = bestSbResult_->objval[order[0]];
945 lb[1] = bestSbResult_->objval[order[1]];
946 BCP_vec<int> tc(2, 0);
947 tc[0] = bestSbResult_->status[order[0]];
948 tc[1] = bestSbResult_->status[order[1]];
949 cands.back()->set_presolve_result(lb, tc);
950 }
951 if (par.entry(BM_par::PrintBranchingInfo)) {
952 print(ifprint2, "BM_lp: branching on variable %i value: %f\n",
953 intBrObj->originalObject()->columnNumber(),
954 intBrObj->value());
955 }
956 }
957 OsiSOSBranchingObject* sosBrObj =
958 dynamic_cast<OsiSOSBranchingObject*>(brObj);
959 if (sosBrObj) {
960 if (sosBrObj->firstBranch() == 1) {
961 order[0] = 1;
962 order[1] = 0;
963 }
964 BCP_lp_sos_branching_object o(sosBrObj);
965 cands.push_back(new BCP_lp_branching_object(nlp, o, order));
966 if (bestSbResult_) {
967 BCP_vec<double> lb(2, 0.0);
968 lb[0] = bestSbResult_->objval[order[0]];
969 lb[1] = bestSbResult_->objval[order[1]];
970 BCP_vec<int> tc(2, 0);
971 tc[0] = bestSbResult_->status[order[0]];
972 tc[1] = bestSbResult_->status[order[1]];
973 cands.back()->set_presolve_result(lb, tc);
974 }
975 if (par.entry(BM_par::PrintBranchingInfo)) {
976 print(ifprint2, "BM_lp: branching on SOS %i value: %f\n",
977 sosBrObj->originalObject()->columnNumber(),
978 sosBrObj->value());
979 }
980 }
981 }
982
983 delete brObj;
984 delete[] clb_old;
985 delete[] cub_old;
986
987 return retCode;
988 }
989
990 /****************************************************************************/
991
992 void
set_user_data_for_children(BCP_presolved_lp_brobj * best,const int selected)993 BM_lp::set_user_data_for_children(BCP_presolved_lp_brobj* best,
994 const int selected)
995 {
996 BM_node* data = NULL;
997 data = new BM_node;
998 data->numNlpFailed_ = numNlpFailed_;
999 best->user_data()[0] = data;
1000 data = new BM_node;
1001 data->numNlpFailed_ = numNlpFailed_;
1002 best->user_data()[1] = data;
1003 }
1004
1005 //#############################################################################
1006
1007 void
process_message(BCP_buffer & buf)1008 BM_lp::process_message(BCP_buffer& buf)
1009 {
1010 int msgtag;
1011 buf.unpack(msgtag);
1012 assert(msgtag == BM_StrongBranchRequest);
1013
1014 Bonmin::OsiTMINLPInterface* nlp = bonmin_.nonlinearSolver();
1015
1016 int numCols;
1017 double* clb;
1018 double* cub;
1019 double objvalOrig;
1020 double cutoff;
1021 buf.unpack(clb, numCols);
1022 assert(numCols == nlp->getNumCols());
1023 buf.unpack(cub, numCols);
1024 assert(numCols == nlp->getNumCols());
1025 buf.unpack(objvalOrig);
1026 buf.unpack(cutoff);
1027 bool has_ws;
1028 buf.unpack(has_ws);
1029 BCP_warmstart* ws = has_ws ? BCP_unpack_warmstart(buf) : NULL;
1030 CoinWarmStart* cws = ws ? ws->convert_to_CoinWarmStart() : NULL;
1031 int location; // just a marker that we'll send back
1032 int numBranch;
1033 buf.unpack(location);
1034 buf.unpack(numBranch);
1035
1036 int i;
1037 BM_BranchData* branchData = new BM_BranchData[numBranch];
1038 for (i = 0; i < numBranch; ++i) {
1039 buf.unpack(branchData[i].changeType);
1040 buf.unpack(branchData[i].objInd);
1041 buf.unpack(branchData[i].colInd);
1042 buf.unpack(branchData[i].solval);
1043 buf.unpack(branchData[i].bd);
1044 }
1045
1046 nlp->setColLower(clb);
1047 nlp->setColUpper(cub);
1048 BM_solve_branches(nlp, cws, numBranch, branchData);
1049
1050 bm_buf.clear();
1051 msgtag = BM_StrongBranchResult;
1052 bm_buf.pack(msgtag);
1053 bm_buf.pack(location);
1054 bm_buf.pack(numBranch);
1055 for (i = 0; i < numBranch; ++i) {
1056 const BM_BranchData& bD = branchData[i];
1057 bm_buf.pack(bD.status);
1058 bm_buf.pack(bD.objval);
1059 bm_buf.pack(bD.iter);
1060 bm_buf.pack(bD.time);
1061 }
1062 send_message(buf.sender(), bm_buf, BCP_Msg_User);
1063
1064 bm_buf.clear();
1065 send_message(parent(), bm_buf, BCP_Msg_SBnodeFinished);
1066
1067 delete[] branchData;
1068 delete cws;
1069 delete ws;
1070 delete[] clb;
1071 delete[] cub;
1072
1073 bm_stats.incNumberSbSolves(numBranch);
1074 }
1075