1 /*===========================================================================*
2 * This file is part of the Bcps Linear Solver (BLIS). *
3 * *
4 * BLIS 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 * Ted Ralphs, Lehigh University *
11 * *
12 * Conceptual Design: *
13 * *
14 * Yan Xu, Lehigh University *
15 * Ted Ralphs, Lehigh University *
16 * Laszlo Ladanyi, IBM T.J. Watson Research Center *
17 * Matthew Saltzman, Clemson University *
18 * *
19 * *
20 * Copyright (C) 2001-2017, Lehigh University, Yan Xu, and Ted Ralphs. *
21 * All Rights Reserved. *
22 *===========================================================================*/
23
24 #include "CoinHelperFunctions.hpp"
25 #include "CoinTime.hpp"
26 #include "CoinWarmStartBasis.hpp"
27
28 #include "Alps.h"
29 #include "AlpsKnowledgeBroker.h"
30
31 #include "BlisBranchStrategyStrong.h"
32 #include "BlisSolution.h"
33 #include "BlisObjectInt.h"
34
35 //#############################################################################
36
37 // Copy constructor
BlisBranchStrategyStrong(const BlisBranchStrategyStrong & rhs)38 BlisBranchStrategyStrong::BlisBranchStrategyStrong (
39 const BlisBranchStrategyStrong & rhs
40 )
41 :
42 BcpsBranchStrategy()
43 {
44 bestChangeUp_ = rhs.bestChangeUp_;
45 bestNumberUp_ = rhs.bestNumberUp_;
46 bestChangeDown_ = rhs.bestChangeDown_;
47 bestNumberDown_ = rhs.bestNumberDown_;
48 }
49
50 //#############################################################################
51
52 /** Create a set of candidate branching objects. */
53 int
createCandBranchObjects(int numPassesLeft)54 BlisBranchStrategyStrong::createCandBranchObjects(int numPassesLeft)
55 {
56 int bStatus = 0;
57 int i, j, pass;
58
59 int col, ind;
60 int numInfs = 0;
61 int numIntegerInfs = 0; // For integer objects.
62 int numObjectInfs = 0; // For non-integer objects.
63 double sumDeg = 0.0; // For solution estimate.
64 double lpX;
65
66 BlisObjectInt *intObject = NULL;
67
68 BlisModel *model = dynamic_cast<BlisModel *>(model_);
69 OsiSolverInterface * solver = model->solver();
70
71 int numCols = model->getNumCols();
72 int numObjects = model->numObjects();
73 bool beforeSolution = (model->getNumSolutions() == 0);
74
75 int givenStrongLen = dynamic_cast<BlisParams*>
76 (model->BlisPar())->entry(BlisParams::strongCandSize);
77 int strongLen = givenStrongLen;
78 int maxStrongLen = CoinMax(CoinMin(givenStrongLen, numObjects), 1);
79
80 BlisStrong * candStrongs = new BlisStrong [maxStrongLen];
81
82 //------------------------------------------------------
83 // Backup solver status
84 //------------------------------------------------------
85
86 // Objective value.
87 double objValue = solver->getObjSense() * solver->getObjValue();
88
89 // Column bounds.
90 const double * lower = solver->getColLower();
91 const double * upper = solver->getColUpper();
92 double * saveUpper = new double[numCols];
93 double * saveLower = new double[numCols];
94 for (i = 0; i < numCols; ++i) {
95 saveLower[i] = lower[i];
96 saveUpper[i] = upper[i];
97 }
98
99 // Primal solution.
100 double * saveSolution = new double[numCols];
101 memcpy(saveSolution, solver->getColSolution(),numCols * sizeof(double));
102
103 //------------------------------------------------------
104 // Select a set of objects based on feasibility.
105 // NOTE: we might go round this loop twice if we are feed in
106 // a "feasible" solution.
107 //------------------------------------------------------
108
109 for (pass = 0; pass < 2; ++pass) {
110
111 // Compute how many infeasible objects.
112 // NOTE: it set model->savedLpSolution
113 model->feasibleSolution(numIntegerInfs, numObjectInfs);
114
115 sumDeg = 0.0;
116 numInfs = 0;
117
118 int preferDir;
119 int leastFrac = 0;
120 double infeasibility;
121 double minInf = ALPS_ZERO;
122 BlisObjectInt * object = NULL;
123
124 for (i = 0; i < maxStrongLen; ++i) {
125 candStrongs[i].bObject = NULL;
126 }
127
128 strongLen = 0;
129
130 for (i = 0; i < numObjects; ++i) {
131
132 // TODO: currently all integer object.
133 object = dynamic_cast<BlisObjectInt *>(model->objects(i));
134 infeasibility = object->infeasibility(model, preferDir);
135
136 if (infeasibility) {
137 ++numInfs;
138
139 // Increase estimated degradation to solution
140 sumDeg += object->pseudocost().getScore();
141
142 // Check for suitability based on infeasibility.
143 if (infeasibility > minInf) {
144
145 if (candStrongs[leastFrac].bObject) {
146 // The slot already has one, free it.
147 delete candStrongs[leastFrac].bObject;
148 }
149
150 // Create new branching object.
151 candStrongs[leastFrac].bObject =
152 object->createBranchObject(model, preferDir);
153
154 candStrongs[leastFrac].bObject->setUpScore(infeasibility);
155 candStrongs[leastFrac].bObject->setObjectIndex(i);
156 candStrongs[leastFrac].objectIndex = i;
157
158 strongLen = CoinMax(strongLen, leastFrac + 1);
159
160 // Find an empty or the worst slot.
161 leastFrac = -1;
162 minInf = ALPS_INFINITY;
163
164 for(j = 0; j < maxStrongLen; ++j) {
165 if (!candStrongs[j].bObject) {
166 // j is an empty slots.
167 minInf = ALPS_ZERO;
168 leastFrac = j;
169 break;
170 }
171 else if(candStrongs[j].bObject->getUpScore() < minInf){
172 minInf = candStrongs[j].bObject->getUpScore();
173 leastFrac = j;
174 }
175 }
176 }
177 }
178 }
179
180 if (numInfs) {
181 #ifdef BLIS_DEBUG_MORE
182 std::cout << "STRONG: numInfs = " << numInfs
183 << std::endl;
184 #endif
185 model->setSolEstimate(objValue + sumDeg);
186
187 break;
188 }
189 else if (pass == 0) {
190 // The first pass and ip feasible
191
192 #ifdef BLIS_DEBUG_MORE
193 std::cout << "STRONG: given a feasible sol" << std::endl;
194 #endif
195
196 bool roundAgain = false;
197
198 CoinWarmStartBasis * ws =
199 dynamic_cast<CoinWarmStartBasis*>(solver->getWarmStart());
200 if (!ws) break;
201
202 // Force solution values within bounds
203 for (i = 0; i < numCols; ++i) {
204 double lpX = saveSolution[i];
205 if (lpX < lower[i]) {
206 saveSolution[i] = lower[i];
207 roundAgain = true;
208 ws->setStructStatus(i, CoinWarmStartBasis::atLowerBound);
209 }
210 else if (lpX > upper[i]) {
211 saveSolution[i] = upper[i];
212 roundAgain = true;
213 ws->setStructStatus(i, CoinWarmStartBasis::atUpperBound);
214 }
215 }
216
217 if (roundAgain) {
218 // Need resolve, then do the second round selection.
219 solver->setWarmStart(ws);
220 delete ws;
221
222 // Resolve.
223 solver->resolve();
224
225 // Save new lp solution.
226 memcpy(saveSolution,
227 solver->getColSolution(),
228 numCols * sizeof(double));
229
230 if (!solver->isProvenOptimal()) {
231 // Become infeasible, can do nothing.
232 bStatus = -2;
233 break;
234 }
235 }
236 else {
237 delete ws;
238 break;
239 }
240 }
241 }
242
243 //------------------------------------------------------
244 // Now we have a set of candidate branching object,
245 // evaluate them.
246 //------------------------------------------------------
247
248 int bestBO = 0;
249 double bestScore = 0.0;
250
251 for (i = 0; i < strongLen; ++i) {
252
253 candStrongs[i].numIntInfUp = numInfs;
254 candStrongs[i].numIntInfDown = numInfs;
255
256 if ( !model->objects(candStrongs[i].bObject->getObjectIndex())->
257 boundBranch(model) ) {
258 // Weild branching: not by modifying variable bounds.
259 // Exit selection routine.
260 strongLen = 0;
261 }
262
263 // Find the most fractional object in case of doing simple branch
264 if (candStrongs[i].bObject->getUpScore() > bestScore) {
265 bestScore = candStrongs[i].bObject->getUpScore();
266 bestBO = i;
267 }
268 }
269
270 // If we have hit max time don't do strong branching
271 double timeLimit = model->AlpsPar()->entry(AlpsParams::timeLimit);
272 bool maxTimeReached = (CoinCpuTime() - model->startTime_ > timeLimit);
273
274 #ifdef BLIS_DEBUG_MORE
275 printf("1. strongLen = %d, maxTimeReached %d, numPassesLeft %d\n",
276 strongLen, maxTimeReached, numPassesLeft);
277 #endif
278
279 if (strongLen <= 0 || maxTimeReached) {
280
281 //--------------------------------------------------
282 // Simple max infeasibility branching.
283 //--------------------------------------------------
284
285 #ifdef BLIS_DEBUG
286 std::cout << "NOT STRONG: maxTimeReached=" << maxTimeReached
287 << "; numPassesLeft=" << numPassesLeft
288 << std::endl;
289 #endif
290
291 // Create candidate set.
292 numBranchObjects_ = 1;
293 branchObjects_ = new BcpsBranchObject * [1];
294 branchObjects_[0] = candStrongs[bestBO].bObject;
295 candStrongs[bestBO].bObject = NULL;
296 strongLen = 0;
297 }
298 else {
299
300 //--------------------------------------------------------
301 /* Strong branching */
302 //--------------------------------------------------------
303
304 // set true to say look at all even if some fixed (experiment)
305 bool solveAll = false;
306 int saveLimit;
307
308 CoinWarmStart * ws = solver->getWarmStart();
309 solver->getIntParam(OsiMaxNumIterationHotStart, saveLimit);
310 if (beforeSolution) {
311 solver->setIntParam(OsiMaxNumIterationHotStart, 10000);
312 }
313
314 // Mark hot start
315 solver->markHotStart();
316
317 #ifdef BLIS_DEBUG_MORE
318 printf("BEFORE LOOP: strongLen = %d\n",strongLen);
319 #endif
320
321 for (i = 0; i < strongLen; ++i) {
322 double objChange;
323 double newObjValue = ALPS_DBL_MAX;
324
325 // status is 0 finished, 1 infeasible and other
326 int lpStatus;
327
328 //----------------------------------------------
329 // Branching down.
330 //----------------------------------------------
331
332 candStrongs[i].bObject->setDirection(-1);
333 candStrongs[i].bObject->branch();
334 solver->solveFromHotStart();
335
336 if (solver->isProvenOptimal()) {
337 lpStatus = 0; // optimal
338 }
339 else if (solver->isIterationLimitReached()
340 &&!solver->isDualObjectiveLimitReached()) {
341 lpStatus = 2; // unknown
342 }
343 else {
344 lpStatus = 1; // infeasible
345 }
346
347 newObjValue = solver->getObjSense() * solver->getObjValue();
348 objChange = newObjValue-objValue ;
349
350 #ifdef BLIS_DEBUG_MORE
351 std::cout << "Down: lpStatus = " << lpStatus << std::endl;
352 #endif
353
354 if (lpStatus == 0) {
355
356 // Update pseudocost
357 ind = candStrongs[i].objectIndex;
358 intObject = dynamic_cast<BlisObjectInt *>(model->objects(ind));
359 col = intObject->columnIndex();
360 lpX = saveSolution[col];
361 intObject->pseudocost().update(-1, objChange, lpX);
362
363 candStrongs[i].finishedDown = true ;
364 if (newObjValue >= model->getCutoff()) {
365 objChange = ALPS_DBL_MAX; // say infeasible
366 }
367 else {
368 if(model->feasibleSolution(candStrongs[i].numIntInfDown)){
369 printf("STRONG:easy:down:found a feasible solution\n");
370 }
371
372 // See if integer solution
373 if (model->feasibleSolution(candStrongs[i].numIntInfDown,
374 candStrongs[i].numObjInfDown)) {
375 #ifdef BLIS_DEBUG_MORE
376 printf("STRONG:down:found a feasible solution\n");
377 #endif
378
379 model->setBestSolution(BLIS_SOL_STRONG,
380 newObjValue,
381 solver->getColSolution());
382 (model->getKnowledgeBroker())->
383 getNodeSelection()->setWeight(0.0);
384 BlisSolution* ksol =
385 new BlisSolution(solver->getNumCols(),
386 solver->getColSolution(),
387 newObjValue);
388 model->getKnowledgeBroker()->
389 addKnowledge(AlpsKnowledgeTypeSolution,
390 ksol,
391 newObjValue);
392
393 objChange = ALPS_DBL_MAX ;
394 }
395 }
396 }
397 else if (lpStatus == 1) {
398 objChange = ALPS_DBL_MAX ;
399 }
400 else {
401 // Can't say much as we did not finish
402 candStrongs[i].finishedDown = false ;
403 }
404
405 candStrongs[i].bObject->setDownScore(objChange);
406
407 // restore bounds
408 int numDiff = 0;
409 for (j = 0; j < numCols; ++j) {
410 if (saveLower[j] != lower[j]) {
411 solver->setColLower(j, saveLower[j]);
412 ++numDiff;
413 }
414 if (saveUpper[j] != upper[j]) {
415 solver->setColUpper(j, saveUpper[j]);
416 ++numDiff;
417 }
418 }
419
420 #ifdef BLIS_DEBUG_MORE
421 std::cout << "numDiff = " << numDiff << std::endl;
422 #endif
423
424 //----------------------------------------------
425 // Branching up.
426 //----------------------------------------------
427
428 candStrongs[i].bObject->branch();
429 solver->solveFromHotStart();
430
431 if (solver->isProvenOptimal()) {
432 lpStatus = 0; // optimal
433 }
434 else if (solver->isIterationLimitReached()
435 &&!solver->isDualObjectiveLimitReached()) {
436 lpStatus = 2; // unknown
437 }
438 else {
439 lpStatus = 1; // infeasible
440 }
441
442 newObjValue = solver->getObjSense() * solver->getObjValue();
443 objChange = newObjValue - objValue ;
444
445 #ifdef BLIS_DEBUG_MORE
446 std::cout << "STRONG: Up: lpStatus = " << lpStatus << std::endl;
447 #endif
448
449 if (lpStatus == 0) {
450 // Update pseudocost
451 ind = candStrongs[i].objectIndex;
452 intObject = dynamic_cast<BlisObjectInt *>(model->objects(ind));
453 col = intObject->columnIndex();
454 lpX = saveSolution[col];
455 intObject->pseudocost().update(1, objChange, lpX);
456
457 candStrongs[i].finishedUp = true ;
458 if (newObjValue >= model->getCutoff()) {
459 objChange = ALPS_DBL_MAX; // Cutoff
460 }
461 else {
462 // See if integer solution
463 if (model->feasibleSolution(candStrongs[i].numIntInfUp,
464 candStrongs[i].numObjInfUp)) {
465
466 #ifdef BLIS_DEBUG_MORE
467 printf("STRONG:up:found a feasible solution\n");
468 #endif
469
470 model->setBestSolution(BLIS_SOL_STRONG,
471 newObjValue,
472 solver->getColSolution());
473
474 model->getKnowledgeBroker()->
475 getNodeSelection()->setWeight(0.0);
476
477 BlisSolution* ksol =
478 new BlisSolution(solver->getNumCols(),
479 solver->getColSolution(),
480 newObjValue);
481
482 model->getKnowledgeBroker()->
483 addKnowledge(AlpsKnowledgeTypeSolution,
484 ksol,
485 newObjValue);
486
487 objChange = ALPS_DBL_MAX;
488 }
489 }
490 }
491 else if (lpStatus == 1) {
492 objChange = ALPS_DBL_MAX;
493 }
494 else {
495 // Can't say much as we did not finish
496 candStrongs[i].finishedUp = false ;
497 }
498 candStrongs[i].bObject->setUpScore(objChange);
499
500 // restore bounds
501 for (j = 0; j < numCols; ++j) {
502 if (saveLower[j] != lower[j]) {
503 solver->setColLower(j,saveLower[j]);
504 }
505 if (saveUpper[j] != upper[j]) {
506 solver->setColUpper(j,saveUpper[j]);
507 }
508 }
509
510 //----------------------------------------------
511 // End of evaluation for this branching object.
512 // Possibilities are:
513 // 1) Both sides below cutoff; this variable is a
514 // candidate for branching.
515 // 2) Both sides infeasible or above the obj cutoff:
516 // no further action here. Break from the evaluation loop and
517 // assume the node will be purged by the caller.
518 // 3) One side below cutoff: Install the branch (i.e., fix the
519 // variable). Break from the evaluation loop and assume the
520 // node will be reoptimised by the caller.
521 //----------------------------------------------
522
523 if (candStrongs[i].bObject->getUpScore() < ALPS_INFINITY) {
524 if(candStrongs[i].bObject->getDownScore() < ALPS_INFINITY) {
525 // feasible - no action
526 }
527 else {
528 // up feasible, down infeasible
529 bStatus = -1;
530 if (!solveAll) {
531 candStrongs[i].bObject->setDirection(1);
532 candStrongs[i].bObject->branch();
533 break;
534 }
535 }
536 }
537 else {
538 if(candStrongs[i].bObject->getDownScore() < ALPS_INFINITY) {
539 // down feasible, up infeasible
540 bStatus = -1;
541 if (!solveAll) {
542 candStrongs[i].bObject->setDirection(-1);
543 candStrongs[i].bObject->branch();
544 break;
545 }
546 }
547 else {
548 // neither side feasible
549 bStatus = -2;
550 break;
551 }
552 }
553 }// EOF the loop for checking each candiate
554
555
556 //--------------------------------------------------
557 // Unmark hotstart and reset lp solver.
558 //--------------------------------------------------
559
560 solver->unmarkHotStart();
561 solver->setIntParam(OsiMaxNumIterationHotStart, saveLimit);
562 solver->setWarmStart(ws);
563 delete ws;
564 }
565
566
567 if (bStatus >= 0) {
568 // Store the set of candidate branching objects.
569 numBranchObjects_ = strongLen;
570
571 branchObjects_ = new BcpsBranchObject* [strongLen];
572 for (i = 0; i < strongLen; ++i) {
573 branchObjects_[i] = candStrongs[i].bObject;
574 candStrongs[i].bObject = NULL;
575 }
576 }
577
578 // Restore solution
579 solver->setColSolution(saveSolution);
580
581 // Cleanup.
582 delete [] saveSolution;
583 delete [] saveLower;
584 delete [] saveUpper;
585
586 for (i = 0; i < maxStrongLen; ++i) {
587 if (candStrongs[i].bObject) {
588 delete candStrongs[i].bObject;
589 }
590 }
591 delete [] candStrongs;
592
593
594 return bStatus;
595 }
596
597 //#############################################################################
598
599 /** Strong method to compare object thisOne to the bestObject_.
600 Compare based on number of infeasibility objects (numInfUp, numInfDn)
601 until a solution is found by search, then swith to change
602 (changeUp, changeDn) in objective.
603 The lesser of the MIN(numInfUp, numInfDown), the better.
604 The larger of the MIN(changeUp, changeDown), the better. */
605 int
betterBranchObject(BcpsBranchObject * thisOne,BcpsBranchObject * bestSoFar)606 BlisBranchStrategyStrong::betterBranchObject(BcpsBranchObject * thisOne,
607 BcpsBranchObject * bestSoFar)
608 {
609 int betterDirection = 0;
610 double bestChange;
611
612 //BlisModel *model = dynamic_cast<BlisModel *>(model_);
613
614 if (!bestSoFar) {
615 bestChange = -1.0;
616 }
617 else {
618 bestChange = ALPS_MIN(bestChangeUp_, bestChangeDown_);
619 }
620
621 double upCost = thisOne->getUpScore();
622 double downCost = thisOne->getDownScore();
623
624 if (upCost <= downCost) {
625 if (upCost > bestChange) {
626 betterDirection = 1; // Branching up.
627 }
628 }
629 else {
630 if (downCost > bestChange) {
631 betterDirection = -1; // Branching down
632 }
633 }
634
635 if (betterDirection) {
636 bestChangeUp_ = upCost;
637 bestChangeDown_ = downCost;
638 }
639
640 return betterDirection;
641 }
642
643 //#############################################################################
644