1 /*===========================================================================*
2 * This file is part of the Abstract Library for Parallel Search (ALPS). *
3 * *
4 * ALPS is distributed under the Eclipse Public License as part of the *
5 * COIN-OR repository (http://www.coin-or.org). *
6 * *
7 * Authors: *
8 * *
9 * Yan Xu, Lehigh University *
10 * Aykut Bulut, Lehigh University *
11 * Ted Ralphs, Lehigh University *
12 * *
13 * Conceptual Design: *
14 * *
15 * Yan Xu, Lehigh University *
16 * Ted Ralphs, Lehigh University *
17 * Laszlo Ladanyi, IBM T.J. Watson Research Center *
18 * Matthew Saltzman, Clemson University *
19 * *
20 * *
21 * Copyright (C) 2001-2019, Lehigh University, Yan Xu, Aykut Bulut, and *
22 * Ted Ralphs. *
23 * All Rights Reserved. *
24 *===========================================================================*/
25
26
27 #include <cassert>
28 #include <iostream>
29 #include <utility>
30 #include <cmath>
31
32 #include "CoinUtility.hpp"
33 #include "OsiSolverInterface.hpp"
34 #include "OsiClpSolverInterface.hpp"
35 #include "OsiRowCut.hpp"
36 #include "OsiColCut.hpp"
37 #include "OsiRowCutDebugger.hpp"
38 #include "OsiCuts.hpp"
39
40 #include "AlpsKnowledge.h"
41 #include "AlpsEnumProcessT.h"
42
43 #include "AbcBranchActual.h"
44 #include "AbcMessage.h"
45 #include "AbcParams.h"
46 #include "AbcTreeNode.h"
47 #include "AbcSolution.h"
48
49 //#############################################################################
50
51 /// Infeasibility - large is 0.5
checkInfeasibility(AbcModel * model,const int columnNumber,int & preferredWay,double & otherWay)52 static double checkInfeasibility(AbcModel* model, const int columnNumber,
53 int & preferredWay, double & otherWay)
54 {
55 OsiSolverInterface * solver = model->solver();
56 const double * solution = model->currentSolution();
57 const double * lower = solver->getColLower();
58 const double * upper = solver->getColUpper();
59 double value = solution[columnNumber];
60 value = std::max(value, lower[columnNumber]);
61 value = std::min(value, upper[columnNumber]);
62 /*printf("%d %g %g %g %g\n",columnNumber_,value,lower[columnNumber_],
63 solution[columnNumber_],upper[columnNumber_]);*/
64 double nearest = floor(value + 0.5);
65 double integerTolerance =
66 model->getDblParam(AbcModel::AbcIntegerTolerance);
67 if (nearest > value)
68 preferredWay = 1;
69 else
70 preferredWay = -1;
71 otherWay = 1.0 - fabs(value - nearest);
72 if (fabs(value - nearest) <= integerTolerance)
73 return 0.0;
74 else
75 return fabs(value - nearest);
76 }
77
78 //#############################################################################
79
80 AlpsTreeNode*
createNewTreeNode(AlpsNodeDesc * & desc) const81 AbcTreeNode::createNewTreeNode(AlpsNodeDesc *&desc) const
82 {
83 // Create a new tree node
84 AbcNodeDesc* d = dynamic_cast<AbcNodeDesc*>(desc);
85 AbcTreeNode* node = new AbcTreeNode(d);
86 desc = 0;
87 return(node);
88 }
89
90 //#############################################################################
91
92 int
process(bool isRoot,bool rampUp)93 AbcTreeNode::process(bool isRoot, bool rampUp)
94 {
95 bool betterSolution = false;
96 double bestValue = broker()->getIncumbentValue();
97 double parentObjValue = getObjValue();
98 double primalTolerance = 1.0e-7;
99 bool cutDuringRampup;
100 AlpsProcessType myType = broker()->getProcType();
101
102 AbcModel *model=dynamic_cast<AbcModel *>(broker()->getModel());
103
104 cutDuringRampup = true;
105
106 if (rampUp) {
107 if (myType == AlpsProcessTypeHub || myType == AlpsProcessTypeMaster){
108 cutDuringRampup=model->AbcPar()->entry(AbcParams::cutDuringRampup);
109 }
110 }
111
112
113 // Check if can quit early. Becare the root.
114 // REMOVE
115 //if (parentObjValue < 1.0e99) bestValue = 3983;
116
117 if (parentObjValue - primalTolerance > bestValue) {
118 setStatus(AlpsNodeStatusFathomed);
119 return false;
120 }
121
122 int i = -1;
123 AbcNodeDesc *desc = dynamic_cast<AbcNodeDesc*>(desc_);
124 AbcModel *m = dynamic_cast<AbcModel*>(broker()->getModel());
125
126 // Update cutoff if recv a better solution from other process
127 if (bestValue < m->getCutoff()) {
128 double cutoff = bestValue -
129 m->getDblParam(AbcModel::AbcCutoffIncrement);
130 m->setCutoff(cutoff);
131 }
132
133 //-------------------------------------------------------------------------
134 // Report search status
135 int nodeInterval = model->AbcPar()->entry(AbcParams::statusInterval);
136
137 assert(nodeInterval > 0);
138
139 if(m->getNodeCount() % nodeInterval == 0){
140 const int nodeLeft = broker()->updateNumNodesLeft();
141 m->messageHandler()->message(ABC_STATUS, m->messages())
142 << broker()->getProcRank() << (m->getNodeCount())
143 << nodeLeft << m->getCutoff() << parentObjValue << CoinMessageEol;
144 }
145
146 //-------------------------------------------------------------------------
147 // Load and solve the lp relaxation.
148 // (1) LP infeasible
149 // a. set status to be fathom
150 // (2) LP feasible
151 // a. MILP feasible. Check whether need update incumbent.
152 // b. LP feasible but not MIP feasible. Check whether can be
153 // fathomed, if not, choose a branch variable
154 //-------------------------------------------------------------------------
155 const int numCols = m->solver()->getNumCols();
156 const double *lbs = desc->lowerBounds();
157 const double *ubs = desc->upperBounds();
158
159 for(i = 0; i < numCols; ++i) {
160 m->solver()->setColBounds(i, lbs[i], ubs[i]);
161 }
162
163 //-------------------------------------------------------------
164 OsiCuts cuts;
165 int maximumCutPassesAtRoot = m->getMaximumCutPassesAtRoot();
166 int numberOldActiveCuts = 0;
167 int numberNewCuts = 0;
168 int maximumWhich = 1000;
169 int *whichGenerator = new int[maximumWhich];
170 int heurFound = -1;
171
172 bool feasible = m->solveWithCuts(cuts, maximumCutPassesAtRoot,
173 NULL, numberOldActiveCuts,
174 numberNewCuts, maximumWhich,
175 whichGenerator,
176 cutDuringRampup,
177 heurFound);
178
179 m->setCurrentNumberCuts(m->currentNumberCuts() + numberNewCuts);
180
181 // if ( (heurFound >= 0) && (m->getObjValue() < bestValue) ) {
182 if (heurFound >= 0) {
183 betterSolution = true;
184 AbcSolution* sol = new AbcSolution(numCols,
185 m->bestSolution(),
186 m->getObjValue());
187 broker()->addKnowledge(AlpsKnowledgeTypeSolution, sol,
188 m->getObjValue());
189 }
190
191
192 //-------------------------------------------------------------
193
194 if (!feasible) {
195 setStatus(AlpsNodeStatusFathomed);
196 setObjValue(-ALPS_OBJ_MAX); // Remove it as soon as possilbe
197 }
198 else {
199 double val = (m->getCurrentObjValue()) * (m->getObjSense());
200 double xS = desc->getBranchedOnValue();
201 int bDir = desc->getBranchedDir();
202 int bInd = desc->getBranchedOn();
203
204 /* Update pseudo: not for the root */
205 if ((m->numberStrong() == 0) && (parent_ != 0) && (bInd >= 0)) {
206 /* FIXME: not sure which bInd < 0 */
207 // std::cout << "bInd=" << bInd << "; bVal=" << xS << std::endl;
208 (m->getPseudoList()[bInd])->update(bDir, quality_, val, xS);
209 }
210
211 int numberInfeasibleIntegers;
212 if(m->feasibleSolution(numberInfeasibleIntegers)) { // MIP feasible
213 if (val < bestValue) {
214 betterSolution = true;
215 m->setBestSolution(ABC_BRANCHSOL,
216 val, m->getColSolution());
217 AbcSolution* ksol = new AbcSolution(numCols,
218 m->getColSolution(),
219 val);
220 broker()->addKnowledge(AlpsKnowledgeTypeSolution,
221 ksol, val);
222 }
223 setStatus(AlpsNodeStatusFathomed);
224 }
225 else {
226 if (val < bestValue) {
227 bool strongFound = false;
228 int action = -1;
229 while (action == -1) {
230 if(broker()->getProcRank() == -1) {
231 std::cout << "*** I AM RANK ONE: before choose:action = " << action
232 << std::endl;
233 }
234 action = chooseBranch(m, strongFound);
235 if ( (strongFound) &&
236 (m->getObjValue() < bestValue) ) {
237 betterSolution = true;
238 AbcSolution *sol = new AbcSolution(numCols,
239 m->bestSolution(),
240 m->getObjValue());
241 broker()->addKnowledge(
242 AlpsKnowledgeTypeSolution, sol, m->getObjValue());
243
244 }
245 if (action == -1) {
246 feasible = m->resolve();
247 //resolved = true ;
248 #if defined(ABC_DEBUG_MORE)
249 printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
250 m->solver()->getObjValue(),
251 m->solver()->getNumRows());
252 #endif
253 if (!feasible) action = -2;
254 }
255 if (action == -2) {
256 feasible = false;
257 }
258 if(broker()->getProcRank() == -1) {
259 std::cout << "*** I AM RANK ONE: action = " << action
260 << std::endl;
261 }
262 }
263 assert(action != -1);
264
265
266 if (action >= 0) {
267 const double * newLbs = m->getColLower();
268 const double * newUbs = m->getColUpper();
269 desc->setLowerBounds(newLbs, numCols);
270 desc->setUpperBounds(newUbs, numCols);
271 #if defined(ABC_DEBUG_MORE)
272 std::cout << "SetPregnant: branchedOn = " << branchedOn_
273 << "; index = " << index_ << std::endl;
274 #endif
275 setStatus(AlpsNodeStatusPregnant);
276 }
277 else if (action == -2) {
278 setStatus(AlpsNodeStatusFathomed);
279 }
280 else {
281 throw CoinError("No branch object found", "process",
282 "AbcTreeNode");
283 }
284
285 }
286 else {
287 setStatus(AlpsNodeStatusFathomed);
288 }
289 }
290 quality_ = val;
291 }
292
293 m->takeOffCuts();
294
295 delete [] whichGenerator;
296
297 return betterSolution;
298 }
299
300 //#############################################################################
301
302 std::vector< CoinTriple<AlpsNodeDesc*, AlpsNodeStatus, double> >
branch()303 AbcTreeNode::branch()
304 {
305 AbcNodeDesc* desc =
306 dynamic_cast<AbcNodeDesc*>(desc_);
307 AbcModel* m = dynamic_cast<AbcModel*>(broker()->getModel());
308
309 double* oldLbs = desc->lowerBounds();
310 double* oldUbs = desc->upperBounds();
311 const int numCols = m->getNumCols();
312 assert(oldLbs && oldUbs && numCols);
313
314 if ( (branchedOn_ < 0) || (branchedOn_ >= numCols) ) {
315 std::cout << "AbcError: branchedOn_ = "<< branchedOn_ << "; numCols = "
316 << numCols << "; index_ = " << index_ << std::endl;
317 throw CoinError("branch index is out of range",
318 "branch", "AbcTreeNode");
319 }
320
321 double* newLbs = new double[numCols];
322 double* newUbs = new double[numCols];
323 std::copy(oldLbs, oldLbs + numCols, newLbs);
324 std::copy(oldUbs, oldUbs + numCols, newUbs);
325 //memcpy(newLbs, oldLbs, sizeof(double) * numCols);
326 //memcpy(newUbs, oldUbs, sizeof(double) * numCols);
327
328 std::vector< CoinTriple<AlpsNodeDesc*, AlpsNodeStatus, double> > newNodes;
329
330 double objVal;
331 // if (getParent() == 0) { //ROOT
332 objVal = getObjValue();
333 //}
334 //else {
335 //objVal = std::min(getParent()->getObjValue(), getObjValue());
336 //}
337
338 // Branch down
339 newLbs[branchedOn_] = oldLbs[branchedOn_];
340 newUbs[branchedOn_] = floor(branchedOnVal_);//floor(branchedOnVal_+1.0e-5);
341 AbcNodeDesc* child;
342 assert(branchedOn_ >= 0);
343 child = new AbcNodeDesc();
344 child->setBroker(broker());
345 child->setLowerBounds(newLbs, numCols);
346 child->setUpperBounds(newUbs, numCols);
347
348 child->setBranchedOn(branchedOn_);
349 child->setBranchedOnValue(branchedOnVal_);
350 child->setBranchedDir(-1);
351 newNodes.push_back(CoinMakeTriple(static_cast<AlpsNodeDesc *>(child),
352 AlpsNodeStatusCandidate,
353 objVal));
354
355 // Branch up
356 newUbs[branchedOn_] = oldUbs[branchedOn_];
357 newLbs[branchedOn_] = ceil(branchedOnVal_);//ceil(branchedOnVal_ - 1.0e-5);
358 child = 0;
359 child = new AbcNodeDesc();
360 child->setBroker(broker());
361 child->setLowerBounds(newLbs, numCols);
362 child->setUpperBounds(newUbs, numCols);
363 child->setBranchedOn(branchedOn_);
364 child->setBranchedOnValue(branchedOnVal_);
365 child->setBranchedDir(1);
366 newNodes.push_back(CoinMakeTriple(static_cast<AlpsNodeDesc *>(child),
367 AlpsNodeStatusCandidate,
368 objVal));
369 if (newLbs != 0) {
370 delete [] newLbs;
371 newLbs = 0;
372 }
373 if (newUbs != 0) {
374 delete [] newUbs;
375 newUbs = 0;
376 }
377 return newNodes;
378 }
379
380 //#############################################################################
381 #if 0
382 // Right now just choose the variable farest to be integral
383 int
384 AbcTreeNode::chooseBranch(AbcModel* model, bool& strongFound)
385 {
386 strongFound = false;
387 int i = -1;
388 const int numInts = model->numberIntegers();
389 const int* intVar = model->integerVariable();
390 const double* sol = model->solver()->getColSolution();
391 double maxDist = model->getDblParam(AbcModel::AbcIntegerTolerance);
392 double value = 0.0;
393 double nearest = 0.0;
394 double distance = 0.0;
395
396 branchedOn_ = -1;
397
398 for (i = 0; i < numInts; ++i) {
399 value = sol[intVar[i]];
400 double nearest = floor(value + 0.5);
401 double distance = fabs(value - nearest);
402
403 if (distance > maxDist) {
404 maxDist = distance;
405 branchedOn_ = intVar[i];
406 branchedOnVal_ = value;
407 }
408 }
409
410 if (branchedOn_ == -1) {
411 return -2; // Should not happend
412 }
413 else {
414 return 2; // Find a branch variable
415 }
416 }
417 #endif
418
419 //#############################################################################
420 //#############################################################################
421
422 //todo(aykut) This can be improved using AlpsTreeNode::encode to pack the
423 //AlpsTreeNode fields.
encode(AlpsEncoded * encoded) const424 AlpsReturnStatus AbcTreeNode::encode(AlpsEncoded * encoded) const {
425 #if defined(ABC_DEBUG_MORE)
426 std::cout << "AbcTreeNode::encode()--start to encode" << std::endl;
427 #endif
428 AbcNodeDesc* desc = dynamic_cast<AbcNodeDesc*>(desc_);
429 AbcModel* model = dynamic_cast<AbcModel*>(broker_->getModel());
430
431 int numCols = model->getNumCols();
432 assert(numCols);
433
434 const double * lb = desc->lowerBounds();
435 const double * ub = desc->upperBounds();
436
437 encoded->writeRep(explicit_);
438 encoded->writeRep(numCols);
439 encoded->writeRep(lb, numCols);
440 encoded->writeRep(ub, numCols);
441 encoded->writeRep(index_);
442 encoded->writeRep(depth_);
443 encoded->writeRep(quality_);
444 encoded->writeRep(parentIndex_);
445 encoded->writeRep(numChildren_);
446 encoded->writeRep(status_);
447 encoded->writeRep(sentMark_);
448 encoded->writeRep(branchedOn_);
449 encoded->writeRep(branchedOnVal_);
450
451 #if defined(ABC_DEBUG_MORE)
452 std::cout << "numCols = " << numCols << "; ";
453 for (int i = 0; i < numCols; ++i) {
454 std::cout << "[" << lb[i] << "," << ub[i] <<"], ";
455 }
456 std::cout << std::endl;
457 std::cout << "index_ = " << index_ << "; ";
458 std::cout << "depth_ = " << depth_ << "; ";
459 std::cout << "quality_ = " << quality_ << "; ";
460 std::cout << "parentIndex_ = " << parentIndex_ << "; ";
461 std::cout << "numChildren_ = " << numChildren_ << std::endl;
462 #endif
463
464 return AlpsReturnStatusOk;
465 }
466
decode(AlpsEncoded & encoded) const467 AlpsKnowledge * AbcTreeNode::decode(AlpsEncoded& encoded) const {
468 int expli;
469 int numCols;
470 double* lb;
471 double* ub;
472 int index;
473 int depth;
474 double objValue;
475 int parentIndex;
476 int numChildren;
477 AlpsNodeStatus nodeStatus;
478 int sentMark;
479 int branchedOn;
480 double branchedOnVal;
481
482 encoded.readRep(expli); // Check whether has full or diff
483 encoded.readRep(numCols);
484 encoded.readRep(lb, numCols);
485 encoded.readRep(ub, numCols);
486 encoded.readRep(index);
487 encoded.readRep(depth);
488 encoded.readRep(objValue);
489 encoded.readRep(parentIndex);
490 encoded.readRep(numChildren);
491 encoded.readRep(nodeStatus);
492 encoded.readRep(sentMark);
493 encoded.readRep(branchedOn);
494 encoded.readRep(branchedOnVal);
495
496
497 if (nodeStatus == AlpsNodeStatusPregnant) {
498 assert(branchedOn >= 0 && branchedOn < numCols);
499 }
500
501 AbcNodeDesc* nodeDesc = new AbcNodeDesc();
502 //nodeDesc->setBroker(desc_->broker());
503 nodeDesc->setLowerBounds(lb, numCols);
504 nodeDesc->setUpperBounds(ub, numCols);
505
506 // nodeDesc->setModel(getKnowledgeBroker()->getDataPool()->getModel());
507 AbcTreeNode* treeNode = new AbcTreeNode(nodeDesc);
508 treeNode->setIndex(index);
509 treeNode->setDepth(depth);
510 treeNode->setObjValue(objValue);
511 treeNode->setParentIndex(parentIndex);
512 treeNode->setParent(0);
513 treeNode->setNumChildren(numChildren);
514 treeNode->setStatus(nodeStatus);
515 treeNode->setSentMark(sentMark);
516 treeNode->setBranchedOn(branchedOn);
517 treeNode->setBranchedOnValue(branchedOnVal);
518
519 if(!lb) {
520 delete [] lb;
521 lb = 0;
522 }
523 if(!ub) {
524 delete [] ub;
525 ub = 0;
526 }
527
528 #if defined(ABC_DEBUG_MORE)
529 std::cout << "numCols = " << numCols << "; ";
530 std::cout << "index = " << index << "; ";
531 std::cout << "depth = " << depth << "; ";
532 std::cout << "objValue = " << objValue<<"; ";
533 std::cout << "parentIndex = " << parentIndex << "; ";
534 std::cout << "status = " << nodeStatus << "; ";
535 std::cout << "branchedOn = " << branchedOn << "; ";
536 std::cout << "branchedOnVal = " << branchedOnVal << "; ";
537 std::cout << "numChildren = " << numChildren << std::endl;
538 #endif
539
540 return treeNode;
541 }
542
543 //#############################################################################
544 // This function tests each integer using strong branching and selects
545 // the one with the least objective degradation. If strong branching is
546 // disabled, the most infeasible integer is choosen
547 // Returns: 2 Find a branching integer
548 // -1 Monotone
549 // -2 Infeasible
chooseBranch(AbcModel * model,bool & strongFound)550 int AbcTreeNode::chooseBranch(AbcModel *model, bool& strongFound)
551 {
552 int numberIntegerInfeasibilities;
553 int anyAction = 0;
554
555 strongFound = false;
556
557 // Make sure we are talking about the same solver
558 bool feasible = model->resolve();
559 if (!feasible) {
560 anyAction = -2;
561 return anyAction;
562 }
563
564 OsiSolverInterface * solver = model->solver();
565 model->feasibleSolution(numberIntegerInfeasibilities);
566 if (numberIntegerInfeasibilities <= 0) {
567 double objectiveValue =
568 solver->getObjSense() * solver->getObjValue();
569 strongFound = model->setBestSolution(ABC_STRONGSOL,
570 objectiveValue,
571 solver->getColSolution());
572 anyAction = -2;
573 return anyAction;
574 }
575
576 double saveObjectiveValue = solver->getObjValue();
577 double objectiveValue = solver->getObjSense() * saveObjectiveValue;
578
579 if(broker()->getProcRank() == -1) {
580 std::cout << "*** I AM RANK FIVE: obj = "<< objectiveValue
581 << ", cutoff = "<< model->getCutoff() << std::endl;
582 }
583
584 if (objectiveValue >= model->getCutoff()) {
585 anyAction = -2;
586 return anyAction;
587 }
588
589 const double * lower = solver->getColLower();
590 const double * upper = solver->getColUpper();
591
592 double integerTolerance =
593 model->getDblParam(AbcModel::AbcIntegerTolerance);
594 int i, j;
595 bool beforeSolution = model->getSolutionCount()==0;
596 int numberStrong = model->numberStrong();
597 int numberObjects = model->numberIntegers();
598 int maximumStrong = std::max(std::min(model->numberStrong(), numberObjects), 1);
599 int numberColumns = model->getNumCols();
600 double * saveUpper = new double[numberColumns];
601 double * saveLower = new double[numberColumns];
602
603 // Save solution in case heuristics need good solution later
604 double * saveSolution = new double[numberColumns];
605 memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
606
607 // Get a branching decision object.
608 AbcBranchDecision * decision = model->branchingMethod();
609 if (!decision) {
610 decision = new AbcBranchDefaultDecision();
611 }
612
613 typedef struct {
614 int possibleBranch; // branching integer variable
615 double upMovement; // cost going up (and initial away from feasible)
616 double downMovement; // cost going down
617 int numIntInfeasUp; // without odd ones
618 int numObjInfeasUp; // just odd ones
619 bool finishedUp; // true if solver finished
620 int numItersUp; // number of iterations in solver
621 int numIntInfeasDown; // without odd ones
622 int numObjInfeasDown; // just odd ones
623 bool finishedDown; // true if solver finished
624 int numItersDown; // number of iterations in solver
625 int objectNumber; // Which object it is
626 } Strong;
627
628 Strong * choice = new Strong[maximumStrong];
629 for (i = 0; i < numberColumns; ++i) {
630 saveLower[i] = lower[i];
631 saveUpper[i] = upper[i];
632 }
633
634 //assert(model->getForcePriority() < 0); // For hot start
635
636 double saveOtherWay = 0.0; // just for non-strong branching
637 numberUnsatisfied_ = 0;
638 int bestPriority = INT_MAX;
639
640 // Scan for branching objects that indicate infeasibility.
641 // Choose the best maximumStrong candidates, using priority as the
642 // first criteria, then integer infeasibility.
643 // The algorithm is to fill the choice array with a set of good
644 // candidates (by infeasibility) with priority bestPriority. Finding
645 // a candidate with priority better (less) than bestPriority flushes
646 // the choice array. (This serves as initialization when the first
647 // candidate is found.) When a candidate is added,
648 // it replaces the candidate with the smallest infeasibility (tracked by
649 // iSmallest)
650 int iSmallest;
651 double mostAway;
652 for (i = 0; i < maximumStrong; ++i) choice[i].possibleBranch = -1;
653 numberStrong = 0;
654 for (i = 0; i < numberObjects; ++i) {
655 const int object = model->integerVariable()[i];
656 int preferredWay;
657 double otherWay;
658 double infeasibility = checkInfeasibility(model, object,
659 preferredWay, otherWay);
660 if (infeasibility > integerTolerance) {
661 if (model->numberStrong() == 0) { /* pseudo cost branching */
662 model->getPseudoIndices()[numberUnsatisfied_] = object;
663 }
664 ++numberUnsatisfied_;
665 int priorityLevel = model->priority(i);
666 if (priorityLevel < bestPriority) {
667 int j;
668 iSmallest = 0;
669 for (j = 0; j < maximumStrong; ++j) {
670 choice[j].upMovement = 0.0;
671 choice[j].possibleBranch = -1;
672 }
673 bestPriority = priorityLevel;
674 mostAway = integerTolerance;
675 numberStrong = 0;
676 }
677 else if (priorityLevel > bestPriority) {
678 continue;
679 }
680
681 // Check for suitability based on infeasibility.
682 if (infeasibility > mostAway) {
683 choice[iSmallest].upMovement = infeasibility;
684 choice[iSmallest].downMovement = otherWay;
685 saveOtherWay = otherWay;
686 choice[iSmallest].possibleBranch = object;
687 numberStrong = std::max(numberStrong, iSmallest + 1);
688 choice[iSmallest].objectNumber = i;
689 int j;
690 iSmallest = -1;
691 mostAway = 1.0e50;
692 for (j = 0; j < maximumStrong; ++j) {
693 if (choice[j].upMovement < mostAway) {
694 mostAway = choice[j].upMovement;
695 iSmallest = j;
696 }
697 }
698 }
699 }
700 }
701
702 #if 0
703 if (numberStrong == 0) { // FIXME: Should not happen
704 bool fea = model->feasibleSolution(numberIntegerInfeasibilities);
705 if (fea) {
706 double objectiveValue =
707 solver->getObjSense() * solver->getObjValue();
708 strongFound = model->setBestSolution(ABC_STRONGSOL,
709 objectiveValue,
710 solver->getColSolution());
711 }
712 //abort();
713 if (!model->branchingMethod()) delete decision;
714 delete [] choice;
715 delete [] saveLower;
716 delete [] saveUpper;
717 // restore solution
718 solver->setColSolution(saveSolution);
719 delete [] saveSolution;
720 anyAction = -2;
721 return anyAction;
722 }
723 #endif
724
725 // Some solvers can do the strong branching calculations faster if
726 // they do them all at once. At present only Clp does for ordinary
727 // integers but I think this coding would be easy to modify
728 bool allNormal = true; // to say if we can do fast strong branching
729 for (i = 0; i < numberStrong; ++i) {
730 choice[i].numIntInfeasUp = numberUnsatisfied_;
731 choice[i].numIntInfeasDown = numberUnsatisfied_;
732 }
733
734 // Setup for strong branching involves saving the current basis (for
735 // restoration afterwards) and setting up for hot starts.
736 if (model->numberStrong() >0) {
737 // set true to say look at all even if some fixed (experiment)
738 bool solveAll = false;
739 // Save basis
740 CoinWarmStart * ws = solver->getWarmStart();
741 // save limit
742 int saveLimit;
743 solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
744 if (beforeSolution) // go to end
745 solver->setIntParam(OsiMaxNumIterationHotStart, 10000);
746
747 // If we are doing all strong branching in one go then we create
748 // new arrays to store information. If clp NULL then doing old way.
749 // Going down -
750 // outputSolution[2*i] is final solution.
751 // outputStuff[2*i] is status (0 - finished, 1 infeas, other unknown)
752 // outputStuff[2*i+numberStrong] is number iterations
753 // On entry newUpper[i] is new upper bound, on exit obj change
754 // Going up -
755 // outputSolution[2*i+1] is final solution.
756 // outputStuff[2*i+1] is status (0 - finished, 1 infeas, other unknown
757 // outputStuff[2*i+1+numberStrong] is number iterations
758 // On entry newLower[i] is new lower bound, on exit obj change
759 OsiClpSolverInterface * osiclp =
760 dynamic_cast< OsiClpSolverInterface*>(solver);
761 ClpSimplex * clp=NULL;
762 double * newLower = NULL;
763 double * newUpper = NULL;
764 double ** outputSolution = NULL;
765 int * outputStuff = NULL;
766 int saveLogLevel;
767
768 //allNormal=false;
769 if (osiclp && allNormal) {
770 clp = osiclp->getModelPtr();
771 saveLogLevel = clp->logLevel();
772 int saveMaxIts = clp->maximumIterations();
773 clp->setMaximumIterations(saveLimit);
774 clp->setLogLevel(0);
775 newLower = new double[numberStrong];
776 newUpper = new double[numberStrong];
777 outputSolution = new double * [2 * numberStrong];
778 outputStuff = new int [4 * numberStrong];
779 int * which = new int [numberStrong];
780 for (i = 0; i < numberStrong; ++i) {
781 int iObject = choice[i].objectNumber;
782 int iSequence = model->integerVariable()[iObject];
783 newLower[i] = ceil(saveSolution[iSequence]);
784 newUpper[i] = floor(saveSolution[iSequence]);
785 which[i] = iSequence;
786 outputSolution[2*i] = new double [numberColumns];
787 outputSolution[2*i+1] = new double [numberColumns];
788 }
789 clp->strongBranching(numberStrong, which,
790 newLower, newUpper,
791 outputSolution, outputStuff,
792 outputStuff + 2 * numberStrong,
793 !solveAll, false);
794 clp->setMaximumIterations(saveMaxIts);
795 delete [] which;
796 }
797 else { // Doing normal way
798 solver->markHotStart();
799 }
800
801 //---------------------------------------------------------------------
802 // Open a loop to do the strong branching LPs. For each candidate
803 // variable, solve an LP with the variable forced down, then up.
804 // If a direction turns out to be infeasible or monotonic (i.e.,
805 // over the dual objective cutoff), force the objective change to
806 // be big (1.0e100). If we determine the problem is infeasible,
807 // or find a monotone variable, escape the loop.
808
809 for (i = 0; i < numberStrong; ++i) {
810 double objectiveChange ;
811 double newObjectiveValue = 1.0e100;
812 // status is 0 finished, 1 infeasible and other
813 int iStatus;
814 int colInd = choice[i].possibleBranch;
815
816 // Try the down direction first.
817 if (!clp) {
818 solver->setColUpper(colInd, floor(saveSolution[colInd]));
819 solver->solveFromHotStart();
820 // restore bounds
821 for (int j = 0; j < numberColumns; ++j) {
822 if (saveLower[j] != lower[j])
823 solver->setColLower(j, saveLower[j]);
824 if (saveUpper[j] != upper[j])
825 solver->setColUpper(j, saveUpper[j]);
826 }
827
828 if (solver->isProvenOptimal())
829 iStatus = 0; // optimal
830 else if (solver->isIterationLimitReached()
831 &&!solver->isDualObjectiveLimitReached())
832 iStatus = 2; // unknown
833 else
834 iStatus = 1; // infeasible
835 newObjectiveValue =
836 solver->getObjSense() * solver->getObjValue();
837 choice[i].numItersDown = solver->getIterationCount();
838 //objectiveChange = newObjectiveValue - objectiveValue ;
839 }
840 else {
841 iStatus = outputStuff[2*i];
842 choice[i].numItersDown = outputStuff[2*numberStrong + 2*i];
843 newObjectiveValue = objectiveValue + newUpper[i];
844 solver->setColSolution(outputSolution[2*i]);
845 }
846 objectiveChange = newObjectiveValue - objectiveValue;
847 if (!iStatus) {
848 choice[i].finishedDown = true;
849 if (newObjectiveValue > model->getCutoff()) {
850 objectiveChange = 1.0e100; // discard it
851 }
852 else {
853 // See if integer solution
854 if (model->feasibleSolution(choice[i].numIntInfeasDown)){
855 strongFound =
856 model->setBestSolution(ABC_STRONGSOL,
857 newObjectiveValue,
858 solver->getColSolution());
859 if (newObjectiveValue > model->getCutoff())
860 objectiveChange = 1.0e100; // discard it
861 }
862 }
863 }
864 else if (iStatus == 1) {
865 objectiveChange = 1.0e100;
866 }
867 else { // Can't say much as we did not finish
868 choice[i].finishedDown = false ;
869 }
870 choice[i].downMovement = objectiveChange ;
871
872 // repeat the whole exercise, forcing the variable up
873 if (!clp) {
874 solver->setColLower(colInd, ceil(saveSolution[colInd]));
875 solver->solveFromHotStart();
876 // restore bounds
877 for (int j = 0; j < numberColumns; j++) {
878 if (saveLower[j] != lower[j])
879 solver->setColLower(j, saveLower[j]);
880 if (saveUpper[j] != upper[j])
881 solver->setColUpper(j, saveUpper[j]);
882 }
883 if (solver->isProvenOptimal())
884 iStatus = 0; // optimal
885 else if (solver->isIterationLimitReached()
886 &&!solver->isDualObjectiveLimitReached())
887 iStatus = 2; // unknown
888 else
889 iStatus = 1; // infeasible
890 newObjectiveValue =
891 solver->getObjSense() * solver->getObjValue();
892 choice[i].numItersUp = solver->getIterationCount();
893 objectiveChange = newObjectiveValue - objectiveValue ;
894 }
895 else {
896 iStatus = outputStuff[2 * i + 1];
897 choice[i].numItersUp = outputStuff[2*numberStrong + 2*i + 1];
898 newObjectiveValue = objectiveValue + newLower[i];
899 solver->setColSolution(outputSolution[2*i + 1]);
900 }
901 objectiveChange = newObjectiveValue - objectiveValue;
902 if (!iStatus) {
903 choice[i].finishedUp = true;
904 if (newObjectiveValue > model->getCutoff()) {
905 objectiveChange = 1.0e100;
906 }
907 else {
908 // See if integer solution
909 if (model->feasibleSolution(choice[i].numIntInfeasUp)){
910 strongFound =
911 model->setBestSolution(ABC_STRONGSOL,
912 newObjectiveValue,
913 solver->getColSolution());
914 if (newObjectiveValue > model->getCutoff())
915 objectiveChange = 1.0e100;
916 }
917 }
918 }
919 else if (iStatus == 1) {
920 objectiveChange = 1.0e100;
921 }
922 else {
923 // Can't say much as we did not finish
924 choice[i].finishedUp = false;
925 }
926 choice[i].upMovement = objectiveChange;
927
928 // End of evaluation for this candidate variable. Possibilities
929 // are:
930 // * Both sides below cutoff; this variable is a candidate
931 // for branching.
932 // * Both sides infeasible or above the objective cutoff:
933 // no further action here. Break from the evaluation loop and
934 // assume the node will be purged by the caller.
935 // * One side below cutoff: Install the branch (i.e., fix the
936 // variable). Break from the evaluation loop and assume
937 // the node will be reoptimised by the caller.
938
939 AbcNodeDesc* desc = dynamic_cast<AbcNodeDesc*>(desc_);
940 int monoIndex;
941
942 if (choice[i].upMovement < 1.0e100) {
943 if(choice[i].downMovement < 1.0e100) {
944 anyAction = 2; // Both
945 }
946 else { // Only up
947 anyAction = -1;
948 monoIndex = choice[i].possibleBranch;
949 solver->setColLower(monoIndex,
950 ceil(saveSolution[monoIndex]));
951 desc->setLowerBound(monoIndex,
952 ceil(saveSolution[monoIndex]));
953 break;
954 }
955 }
956 else {
957 if(choice[i].downMovement < 1.0e100) {
958 anyAction = -1; // Only down
959 monoIndex = choice[i].possibleBranch;
960 solver->setColUpper(monoIndex,
961 floor(saveSolution[monoIndex]));
962 desc->setUpperBound(monoIndex,
963 floor(saveSolution[monoIndex]));
964 break;
965 }
966 else { // neither
967 anyAction = -2;
968 break;
969 }
970 }
971 }
972
973 if (!clp) {
974 solver->unmarkHotStart();
975 } else {
976 clp->setLogLevel(saveLogLevel);
977 delete [] newLower;
978 delete [] newUpper;
979 delete [] outputStuff;
980 for (int i = 0; i < 2*numberStrong; ++i)
981 delete [] outputSolution[i];
982 delete [] outputSolution;
983 }
984
985 solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
986 solver->setWarmStart(ws); // restore basis
987 delete ws;
988
989 if(broker()->getProcRank() == -1) {
990 std::cout << "*** I AM RANK ONE: chooseBranch():anyAction = "
991 << anyAction << "numberStrong = "<< numberStrong
992 << std::endl;
993 }
994
995 // Sift through the candidates for the best one.
996 // QUERY: Setting numberNodes looks to be a distributed noop.
997 // numberNodes is local to this code block. Perhaps should be
998 // numberNodes_ from model?
999 // Unclear what this calculation is doing.
1000 if (anyAction > 0) {
1001 int numberNodes = model->getNodeCount();
1002 // get average cost per iteration and assume stopped ones
1003 // would stop after 50% more iterations at average cost??? !!! ???
1004 double averageCostPerIteration = 0.0;
1005 double totalNumberIterations = 1.0;
1006 int smallestNumberInfeasibilities = INT_MAX;
1007 for (i = 0; i < numberStrong; ++i) {
1008 totalNumberIterations += choice[i].numItersDown +
1009 choice[i].numItersUp;
1010 averageCostPerIteration += choice[i].downMovement +
1011 choice[i].upMovement;
1012 smallestNumberInfeasibilities =
1013 std::min(std::min(choice[i].numIntInfeasDown,
1014 choice[i].numIntInfeasUp),
1015 smallestNumberInfeasibilities);
1016 }
1017 if (smallestNumberInfeasibilities >= numberIntegerInfeasibilities)
1018 numberNodes = 1000000; // switch off search for better solution
1019 numberNodes = 1000000; // switch off anyway
1020 averageCostPerIteration /= totalNumberIterations;
1021 // all feasible - choose best bet
1022
1023 // New method does all at once so it can be more sophisticated
1024 // in deciding how to balance actions.
1025 // But it does need arrays
1026 double * changeUp = new double [numberStrong];
1027 int * numberInfeasibilitiesUp = new int [numberStrong];
1028 double * changeDown = new double [numberStrong];
1029 int * numberInfeasibilitiesDown = new int [numberStrong];
1030 int * objects = new int [numberStrong];
1031
1032 for (i = 0; i < numberStrong; ++i) {
1033
1034 int iColumn = choice[i].possibleBranch;
1035
1036 model->messageHandler()->message(ABC_STRONG,model->messages())
1037 << i << iColumn
1038 << choice[i].downMovement << choice[i].numIntInfeasDown
1039 << choice[i].upMovement << choice[i].numIntInfeasUp
1040 << choice[i].possibleBranch
1041 << CoinMessageEol;
1042
1043 changeUp[i] = choice[i].upMovement;
1044 numberInfeasibilitiesUp[i] = choice[i].numIntInfeasUp;
1045 changeDown[i] = choice[i].downMovement;
1046 numberInfeasibilitiesDown[i] = choice[i].numIntInfeasDown;
1047 objects[i] = choice[i].possibleBranch;
1048 }
1049 int whichObject =
1050 decision->bestBranch(model,
1051 objects, numberStrong,
1052 numberUnsatisfied_, changeUp,
1053 numberInfeasibilitiesUp, changeDown,
1054 numberInfeasibilitiesDown,objectiveValue);
1055
1056 // move branching object and make sure it will not be deleted
1057 if (whichObject >= 0) {
1058 branchedOn_ = objects[whichObject];
1059 branchedOnVal_= saveSolution[branchedOn_];
1060 assert(branchedOn_ >= 0 && branchedOn_ < numberColumns);
1061
1062 #if 0 // FIXME
1063 std::cout << "AbcStrongBranch: col index : ";
1064 for (i = 0; i < numberStrong; ++i) {
1065 std::cout<< objects[i] << " ";
1066 }
1067 std::cout <<"\nAbcStrongBranch: branchedOn_ = " << branchedOn_
1068 << "; branchedOnVal_ = " << branchedOnVal_
1069 << std::endl;
1070 #endif
1071 }
1072 else {
1073 std::cout << "AbcError: whichObject = " << whichObject
1074 << "; numberStrong = " << numberStrong << std::endl;
1075 throw CoinError("bestBranch find nothing",
1076 "chooseBranch",
1077 "AbcTreeNode");
1078 }
1079
1080 delete [] changeUp;
1081 delete [] numberInfeasibilitiesUp;
1082 delete [] changeDown;
1083 delete [] numberInfeasibilitiesDown;
1084 delete [] objects;
1085 }
1086
1087 if (anyAction == 0){
1088 throw CoinError("Can't find candidates",
1089 "chooseBranch",
1090 "AbcTreeNode");
1091 }
1092 }
1093 else { // pseudocost branching
1094 #if 1
1095 int object;
1096 int mostInd = -1;
1097 double mostDeg = -1.0;
1098 double deg, upCost, downCost, fraction;
1099 AbcPseudocost *pseudoC;
1100
1101 for(j = 0; j < numberUnsatisfied_; ++j) {
1102 //for(j = 0; j < numberStrong; ++j) {
1103 object = model->getPseudoIndices()[j];
1104 //object = choice[j].possibleBranch;
1105 fraction = saveSolution[object] - floor(saveSolution[object]);
1106 pseudoC = model->getPseudoList()[object];
1107 upCost = pseudoC->upCost_ * (1.0 - fraction);
1108 downCost = pseudoC->downCost_ * fraction;
1109 deg = 4.0 * std::min(upCost, downCost) + 1.0 * std::max(upCost, downCost);
1110 if (deg > mostDeg) {
1111 mostDeg = deg;
1112 mostInd = object;
1113 }
1114 }
1115
1116 branchedOn_ = mostInd;
1117 #else
1118 branchedOn_ = choice[0].possibleBranch;
1119 #endif
1120 branchedOnVal_= saveSolution[branchedOn_];
1121 assert( (branchedOn_ >= 0) && (branchedOn_ < numberColumns));
1122 }
1123
1124 if (!model->branchingMethod()) delete decision;
1125
1126 delete [] choice;
1127 delete [] saveLower;
1128 delete [] saveUpper;
1129
1130 // restore solution
1131 solver->setColSolution(saveSolution);
1132 delete [] saveSolution;
1133
1134 return anyAction;
1135 }
1136