1 // Copyright (C) 2006, International Business Machines
2 // Corporation and others. All Rights Reserved.
3 // This code is licensed under the terms of the Eclipse Public License (EPL).
4
5 #include <string>
6 #include <cassert>
7 #include <cfloat>
8 #include <cmath>
9 #include "CoinPragma.hpp"
10 #include "OsiSolverInterface.hpp"
11 #include "OsiAuxInfo.hpp"
12 #include "OsiSolverBranch.hpp"
13 #include "CoinWarmStartBasis.hpp"
14 #include "CoinPackedMatrix.hpp"
15 #include "CoinTime.hpp"
16 #include "CoinSort.hpp"
17 #include "CoinFinite.hpp"
18 #include "OsiChooseVariable.hpp"
19 using namespace std;
20
OsiChooseVariable()21 OsiChooseVariable::OsiChooseVariable() :
22 goodObjectiveValue_(COIN_DBL_MAX),
23 upChange_(0.0),
24 downChange_(0.0),
25 goodSolution_(NULL),
26 list_(NULL),
27 useful_(NULL),
28 solver_(NULL),
29 status_(-1),
30 bestObjectIndex_(-1),
31 bestWhichWay_(-1),
32 firstForcedObjectIndex_(-1),
33 firstForcedWhichWay_(-1),
34 numberUnsatisfied_(0),
35 numberStrong_(0),
36 numberOnList_(0),
37 numberStrongDone_(0),
38 numberStrongIterations_(0),
39 numberStrongFixed_(0),
40 trustStrongForBound_(true),
41 trustStrongForSolution_(true)
42 {
43 }
44
OsiChooseVariable(const OsiSolverInterface * solver)45 OsiChooseVariable::OsiChooseVariable(const OsiSolverInterface * solver) :
46 goodObjectiveValue_(COIN_DBL_MAX),
47 upChange_(0.0),
48 downChange_(0.0),
49 goodSolution_(NULL),
50 solver_(solver),
51 status_(-1),
52 bestObjectIndex_(-1),
53 bestWhichWay_(-1),
54 firstForcedObjectIndex_(-1),
55 firstForcedWhichWay_(-1),
56 numberUnsatisfied_(0),
57 numberStrong_(0),
58 numberOnList_(0),
59 numberStrongDone_(0),
60 numberStrongIterations_(0),
61 numberStrongFixed_(0),
62 trustStrongForBound_(true),
63 trustStrongForSolution_(true)
64 {
65 // create useful arrays
66 int numberObjects = solver_->numberObjects();
67 list_ = new int [numberObjects];
68 useful_ = new double [numberObjects];
69 }
70
OsiChooseVariable(const OsiChooseVariable & rhs)71 OsiChooseVariable::OsiChooseVariable(const OsiChooseVariable & rhs)
72 {
73 goodObjectiveValue_ = rhs.goodObjectiveValue_;
74 upChange_ = rhs.upChange_;
75 downChange_ = rhs.downChange_;
76 status_ = rhs.status_;
77 bestObjectIndex_ = rhs.bestObjectIndex_;
78 bestWhichWay_ = rhs.bestWhichWay_;
79 firstForcedObjectIndex_ = rhs.firstForcedObjectIndex_;
80 firstForcedWhichWay_ = rhs.firstForcedWhichWay_;
81 numberUnsatisfied_ = rhs.numberUnsatisfied_;
82 numberStrong_ = rhs.numberStrong_;
83 numberOnList_ = rhs.numberOnList_;
84 numberStrongDone_ = rhs.numberStrongDone_;
85 numberStrongIterations_ = rhs.numberStrongIterations_;
86 numberStrongFixed_ = rhs.numberStrongFixed_;
87 trustStrongForBound_ = rhs.trustStrongForBound_;
88 trustStrongForSolution_ = rhs.trustStrongForSolution_;
89 solver_ = rhs.solver_;
90 if (solver_) {
91 int numberObjects = solver_->numberObjects();
92 int numberColumns = solver_->getNumCols();
93 if (rhs.goodSolution_) {
94 goodSolution_ = CoinCopyOfArray(rhs.goodSolution_,numberColumns);
95 } else {
96 goodSolution_ = NULL;
97 }
98 list_ = CoinCopyOfArray(rhs.list_,numberObjects);
99 useful_ = CoinCopyOfArray(rhs.useful_,numberObjects);
100 } else {
101 goodSolution_ = NULL;
102 list_ = NULL;
103 useful_ = NULL;
104 }
105 }
106
107 OsiChooseVariable &
operator =(const OsiChooseVariable & rhs)108 OsiChooseVariable::operator=(const OsiChooseVariable & rhs)
109 {
110 if (this != &rhs) {
111 delete [] goodSolution_;
112 delete [] list_;
113 delete [] useful_;
114 goodObjectiveValue_ = rhs.goodObjectiveValue_;
115 upChange_ = rhs.upChange_;
116 downChange_ = rhs.downChange_;
117 status_ = rhs.status_;
118 bestObjectIndex_ = rhs.bestObjectIndex_;
119 bestWhichWay_ = rhs.bestWhichWay_;
120 firstForcedObjectIndex_ = rhs.firstForcedObjectIndex_;
121 firstForcedWhichWay_ = rhs.firstForcedWhichWay_;
122 numberUnsatisfied_ = rhs.numberUnsatisfied_;
123 numberStrong_ = rhs.numberStrong_;
124 numberOnList_ = rhs.numberOnList_;
125 numberStrongDone_ = rhs.numberStrongDone_;
126 numberStrongIterations_ = rhs.numberStrongIterations_;
127 numberStrongFixed_ = rhs.numberStrongFixed_;
128 trustStrongForBound_ = rhs.trustStrongForBound_;
129 trustStrongForSolution_ = rhs.trustStrongForSolution_;
130 solver_ = rhs.solver_;
131 if (solver_) {
132 int numberObjects = solver_->numberObjects();
133 int numberColumns = solver_->getNumCols();
134 if (rhs.goodSolution_) {
135 goodSolution_ = CoinCopyOfArray(rhs.goodSolution_,numberColumns);
136 } else {
137 goodSolution_ = NULL;
138 }
139 list_ = CoinCopyOfArray(rhs.list_,numberObjects);
140 useful_ = CoinCopyOfArray(rhs.useful_,numberObjects);
141 } else {
142 goodSolution_ = NULL;
143 list_ = NULL;
144 useful_ = NULL;
145 }
146 }
147 return *this;
148 }
149
150
~OsiChooseVariable()151 OsiChooseVariable::~OsiChooseVariable ()
152 {
153 delete [] goodSolution_;
154 delete [] list_;
155 delete [] useful_;
156 }
157
158 // Clone
159 OsiChooseVariable *
clone() const160 OsiChooseVariable::clone() const
161 {
162 return new OsiChooseVariable(*this);
163 }
164 // Set solver and redo arrays
165 void
setSolver(const OsiSolverInterface * solver)166 OsiChooseVariable::setSolver (const OsiSolverInterface * solver)
167 {
168 solver_ = solver;
169 delete [] list_;
170 delete [] useful_;
171 // create useful arrays
172 int numberObjects = solver_->numberObjects();
173 list_ = new int [numberObjects];
174 useful_ = new double [numberObjects];
175 }
176
177
178 // Initialize
179 int
setupList(OsiBranchingInformation * info,bool initialize)180 OsiChooseVariable::setupList ( OsiBranchingInformation *info, bool initialize)
181 {
182 if (initialize) {
183 status_=-2;
184 delete [] goodSolution_;
185 bestObjectIndex_=-1;
186 numberStrongDone_=0;
187 numberStrongIterations_ = 0;
188 numberStrongFixed_ = 0;
189 goodSolution_ = NULL;
190 goodObjectiveValue_ = COIN_DBL_MAX;
191 }
192 numberOnList_=0;
193 numberUnsatisfied_=0;
194 int numberObjects = solver_->numberObjects();
195 assert (numberObjects);
196 double check = 0.0;
197 int checkIndex=0;
198 int bestPriority=COIN_INT_MAX;
199 // pretend one strong even if none
200 int maximumStrong= numberStrong_ ? CoinMin(numberStrong_,numberObjects) : 1;
201 int putOther = numberObjects;
202 int i;
203 for (i=0;i<maximumStrong;i++) {
204 list_[i]=-1;
205 useful_[i]=0.0;
206 }
207 OsiObject ** object = info->solver_->objects();
208 // Say feasible
209 bool feasible = true;
210 for ( i=0;i<numberObjects;i++) {
211 int way;
212 double value = object[i]->infeasibility(info,way);
213 if (value>0.0) {
214 numberUnsatisfied_++;
215 if (value==COIN_DBL_MAX) {
216 // infeasible
217 feasible=false;
218 break;
219 }
220 int priorityLevel = object[i]->priority();
221 // Better priority? Flush choices.
222 if (priorityLevel<bestPriority) {
223 for (int j=0;j<maximumStrong;j++) {
224 if (list_[j]>=0) {
225 int iObject = list_[j];
226 list_[j]=-1;
227 useful_[j]=0.0;
228 list_[--putOther]=iObject;
229 }
230 }
231 bestPriority = priorityLevel;
232 check=0.0;
233 }
234 if (priorityLevel==bestPriority) {
235 if (value>check) {
236 //add to list
237 int iObject = list_[checkIndex];
238 if (iObject>=0)
239 list_[--putOther]=iObject; // to end
240 list_[checkIndex]=i;
241 useful_[checkIndex]=value;
242 // find worst
243 check=COIN_DBL_MAX;
244 for (int j=0;j<maximumStrong;j++) {
245 if (list_[j]>=0) {
246 if (useful_[j]<check) {
247 check=useful_[j];
248 checkIndex=j;
249 }
250 } else {
251 check=0.0;
252 checkIndex = j;
253 break;
254 }
255 }
256 } else {
257 // to end
258 list_[--putOther]=i;
259 }
260 } else {
261 // to end
262 list_[--putOther]=i;
263 }
264 }
265 }
266 // Get list
267 numberOnList_=0;
268 if (feasible) {
269 for (i=0;i<maximumStrong;i++) {
270 if (list_[i]>=0) {
271 list_[numberOnList_]=list_[i];
272 useful_[numberOnList_++]=-useful_[i];
273 }
274 }
275 if (numberOnList_) {
276 // Sort
277 CoinSort_2(useful_,useful_+numberOnList_,list_);
278 // move others
279 i = numberOnList_;
280 for (;putOther<numberObjects;putOther++)
281 list_[i++]=list_[putOther];
282 assert (i==numberUnsatisfied_);
283 if (!numberStrong_)
284 numberOnList_=0;
285 }
286 } else {
287 // not feasible
288 numberUnsatisfied_=-1;
289 }
290 return numberUnsatisfied_;
291 }
292 /* Choose a variable
293 Returns -
294 -1 Node is infeasible
295 0 Normal termination - we have a candidate
296 1 All looks satisfied - no candidate
297 2 We can change the bound on a variable - but we also have a strong branching candidate
298 3 We can change the bound on a variable - but we have a non-strong branching candidate
299 4 We can change the bound on a variable - no other candidates
300 We can pick up branch from whichObject() and whichWay()
301 We can pick up a forced branch (can change bound) from whichForcedObject() and whichForcedWay()
302 If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
303 */
304 int
chooseVariable(OsiSolverInterface * solver,OsiBranchingInformation *,bool)305 OsiChooseVariable::chooseVariable( OsiSolverInterface * solver, OsiBranchingInformation *, bool )
306 {
307 if (numberUnsatisfied_) {
308 bestObjectIndex_=list_[0];
309 bestWhichWay_ = solver->object(bestObjectIndex_)->whichWay();
310 firstForcedObjectIndex_ = -1;
311 firstForcedWhichWay_ =-1;
312 return 0;
313 } else {
314 return 1;
315 }
316 }
317 // Returns true if solution looks feasible against given objects
318 bool
feasibleSolution(const OsiBranchingInformation * info,const double * solution,int numberObjects,const OsiObject ** objects)319 OsiChooseVariable::feasibleSolution(const OsiBranchingInformation * info,
320 const double * solution,
321 int numberObjects,
322 const OsiObject ** objects)
323 {
324 bool satisfied=true;
325 const double * saveSolution = info->solution_;
326 info->solution_ = solution;
327 for (int i=0;i<numberObjects;i++) {
328 double value = objects[i]->checkInfeasibility(info);
329 if (value>0.0) {
330 satisfied=false;
331 break;
332 }
333 }
334 info->solution_ = saveSolution;
335 return satisfied;
336 }
337 // Saves a good solution
338 void
saveSolution(const OsiSolverInterface * solver)339 OsiChooseVariable::saveSolution(const OsiSolverInterface * solver)
340 {
341 delete [] goodSolution_;
342 int numberColumns = solver->getNumCols();
343 goodSolution_ = CoinCopyOfArray(solver->getColSolution(),numberColumns);
344 goodObjectiveValue_ = solver->getObjSense()*solver->getObjValue();
345 }
346 // Clears out good solution after use
347 void
clearGoodSolution()348 OsiChooseVariable::clearGoodSolution()
349 {
350 delete [] goodSolution_;
351 goodSolution_ = NULL;
352 goodObjectiveValue_ = COIN_DBL_MAX;
353 }
354
355 /* This is a utility function which does strong branching on
356 a list of objects and stores the results in OsiHotInfo.objects.
357 On entry the object sequence is stored in the OsiHotInfo object
358 and maybe more.
359 It returns -
360 -1 - one branch was infeasible both ways
361 0 - all inspected - nothing can be fixed
362 1 - all inspected - some can be fixed (returnCriterion==0)
363 2 - may be returning early - one can be fixed (last one done) (returnCriterion==1)
364 3 - returning because max time
365 */
366 int
doStrongBranching(OsiSolverInterface * solver,OsiBranchingInformation * info,int numberToDo,int returnCriterion)367 OsiChooseStrong::doStrongBranching( OsiSolverInterface * solver,
368 OsiBranchingInformation *info,
369 int numberToDo, int returnCriterion)
370 {
371
372 // Might be faster to extend branch() to return bounds changed
373 double * saveLower = NULL;
374 double * saveUpper = NULL;
375 int numberColumns = solver->getNumCols();
376 solver->markHotStart();
377 const double * lower = info->lower_;
378 const double * upper = info->upper_;
379 saveLower = CoinCopyOfArray(info->lower_,numberColumns);
380 saveUpper = CoinCopyOfArray(info->upper_,numberColumns);
381 numResults_=0;
382 int returnCode=0;
383 double timeStart = CoinCpuTime();
384 for (int iDo=0;iDo<numberToDo;iDo++) {
385 OsiHotInfo * result = results_ + iDo;
386 // For now just 2 way
387 OsiBranchingObject * branch = result->branchingObject();
388 assert (branch->numberBranches()==2);
389 /*
390 Try the first direction. Each subsequent call to branch() performs the
391 specified branch and advances the branch object state to the next branch
392 alternative.)
393 */
394 OsiSolverInterface * thisSolver = solver;
395 if (branch->boundBranch()) {
396 // ordinary
397 branch->branch(solver);
398 // maybe we should check bounds for stupidities here?
399 solver->solveFromHotStart() ;
400 } else {
401 // adding cuts or something
402 thisSolver = solver->clone();
403 branch->branch(thisSolver);
404 // set hot start iterations
405 int limit;
406 thisSolver->getIntParam(OsiMaxNumIterationHotStart,limit);
407 thisSolver->setIntParam(OsiMaxNumIteration,limit);
408 thisSolver->resolve();
409 }
410 // can check if we got solution
411 // status is 0 finished, 1 infeasible and 2 unfinished and 3 is solution
412 int status0 = result->updateInformation(thisSolver,info,this);
413 numberStrongIterations_ += thisSolver->getIterationCount();
414 if (status0==3) {
415 // new solution already saved
416 if (trustStrongForSolution_) {
417 info->cutoff_ = goodObjectiveValue_;
418 status0=0;
419 }
420 }
421 if (solver!=thisSolver)
422 delete thisSolver;
423 // Restore bounds
424 for (int j=0;j<numberColumns;j++) {
425 if (saveLower[j] != lower[j])
426 solver->setColLower(j,saveLower[j]);
427 if (saveUpper[j] != upper[j])
428 solver->setColUpper(j,saveUpper[j]);
429 }
430 /*
431 Try the next direction
432 */
433 thisSolver = solver;
434 if (branch->boundBranch()) {
435 // ordinary
436 branch->branch(solver);
437 // maybe we should check bounds for stupidities here?
438 solver->solveFromHotStart() ;
439 } else {
440 // adding cuts or something
441 thisSolver = solver->clone();
442 branch->branch(thisSolver);
443 // set hot start iterations
444 int limit;
445 thisSolver->getIntParam(OsiMaxNumIterationHotStart,limit);
446 thisSolver->setIntParam(OsiMaxNumIteration,limit);
447 thisSolver->resolve();
448 }
449 // can check if we got solution
450 // status is 0 finished, 1 infeasible and 2 unfinished and 3 is solution
451 int status1 = result->updateInformation(thisSolver,info,this);
452 numberStrongDone_++;
453 numberStrongIterations_ += thisSolver->getIterationCount();
454 if (status1==3) {
455 // new solution already saved
456 if (trustStrongForSolution_) {
457 info->cutoff_ = goodObjectiveValue_;
458 status1=0;
459 }
460 }
461 if (solver!=thisSolver)
462 delete thisSolver;
463 // Restore bounds
464 for (int j=0;j<numberColumns;j++) {
465 if (saveLower[j] != lower[j])
466 solver->setColLower(j,saveLower[j]);
467 if (saveUpper[j] != upper[j])
468 solver->setColUpper(j,saveUpper[j]);
469 }
470 /*
471 End of evaluation for this candidate variable. Possibilities are:
472 * Both sides below cutoff; this variable is a candidate for branching.
473 * Both sides infeasible or above the objective cutoff: no further action
474 here. Break from the evaluation loop and assume the node will be purged
475 by the caller.
476 * One side below cutoff: Install the branch (i.e., fix the variable). Possibly break
477 from the evaluation loop and assume the node will be reoptimised by the
478 caller.
479 */
480 numResults_++;
481 if (status0==1&&status1==1) {
482 // infeasible
483 returnCode=-1;
484 break; // exit loop
485 } else if (status0==1||status1==1) {
486 numberStrongFixed_++;
487 if (!returnCriterion) {
488 returnCode=1;
489 } else {
490 returnCode=2;
491 break;
492 }
493 }
494 bool hitMaxTime = ( CoinCpuTime()-timeStart > info->timeRemaining_);
495 if (hitMaxTime) {
496 returnCode=3;
497 break;
498 }
499 }
500 delete [] saveLower;
501 delete [] saveUpper;
502 // Delete the snapshot
503 solver->unmarkHotStart();
504 return returnCode;
505 }
506
507 // Given a candidate fill in useful information e.g. estimates
508 void
updateInformation(const OsiBranchingInformation * info,int,OsiHotInfo * hotInfo)509 OsiChooseVariable::updateInformation(const OsiBranchingInformation *info,
510 int , OsiHotInfo * hotInfo)
511 {
512 int index = hotInfo->whichObject();
513 assert (index<solver_->numberObjects());
514 //assert (branch<2);
515 OsiObject ** object = info->solver_->objects();
516 upChange_ = object[index]->upEstimate();
517 downChange_ = object[index]->downEstimate();
518 }
519 #if 1
520 // Given a branch fill in useful information e.g. estimates
521 void
updateInformation(int index,int branch,double,double,int)522 OsiChooseVariable::updateInformation( int index, int branch,
523 double , double ,
524 int )
525 {
526 assert (index<solver_->numberObjects());
527 assert (branch<2);
528 OsiObject ** object = solver_->objects();
529 if (branch)
530 upChange_ = object[index]->upEstimate();
531 else
532 downChange_ = object[index]->downEstimate();
533 }
534 #endif
535
536 //##############################################################################
537
538 void
gutsOfDelete()539 OsiPseudoCosts::gutsOfDelete()
540 {
541 if (numberObjects_ > 0) {
542 numberObjects_ = 0;
543 numberBeforeTrusted_ = 0;
544 delete[] upTotalChange_; upTotalChange_ = NULL;
545 delete[] downTotalChange_; downTotalChange_ = NULL;
546 delete[] upNumber_; upNumber_ = NULL;
547 delete[] downNumber_; downNumber_ = NULL;
548 }
549 }
550
551 void
gutsOfCopy(const OsiPseudoCosts & rhs)552 OsiPseudoCosts::gutsOfCopy(const OsiPseudoCosts& rhs)
553 {
554 numberObjects_ = rhs.numberObjects_;
555 numberBeforeTrusted_ = rhs.numberBeforeTrusted_;
556 if (numberObjects_ > 0) {
557 upTotalChange_ = CoinCopyOfArray(rhs.upTotalChange_,numberObjects_);
558 downTotalChange_ = CoinCopyOfArray(rhs.downTotalChange_,numberObjects_);
559 upNumber_ = CoinCopyOfArray(rhs.upNumber_,numberObjects_);
560 downNumber_ = CoinCopyOfArray(rhs.downNumber_,numberObjects_);
561 }
562 }
563
OsiPseudoCosts()564 OsiPseudoCosts::OsiPseudoCosts() :
565 upTotalChange_(NULL),
566 downTotalChange_(NULL),
567 upNumber_(NULL),
568 downNumber_(NULL),
569 numberObjects_(0),
570 numberBeforeTrusted_(0)
571 {
572 }
573
~OsiPseudoCosts()574 OsiPseudoCosts::~OsiPseudoCosts()
575 {
576 gutsOfDelete();
577 }
578
OsiPseudoCosts(const OsiPseudoCosts & rhs)579 OsiPseudoCosts::OsiPseudoCosts(const OsiPseudoCosts& rhs) :
580 upTotalChange_(NULL),
581 downTotalChange_(NULL),
582 upNumber_(NULL),
583 downNumber_(NULL),
584 numberObjects_(0),
585 numberBeforeTrusted_(0)
586 {
587 gutsOfCopy(rhs);
588 }
589
590 OsiPseudoCosts&
operator =(const OsiPseudoCosts & rhs)591 OsiPseudoCosts::operator=(const OsiPseudoCosts& rhs)
592 {
593 if (this != &rhs) {
594 gutsOfDelete();
595 gutsOfCopy(rhs);
596 }
597 return *this;
598 }
599
600 void
initialize(int n)601 OsiPseudoCosts::initialize(int n)
602 {
603 gutsOfDelete();
604 numberObjects_ = n;
605 if (numberObjects_ > 0) {
606 upTotalChange_ = new double [numberObjects_];
607 downTotalChange_ = new double [numberObjects_];
608 upNumber_ = new int [numberObjects_];
609 downNumber_ = new int [numberObjects_];
610 CoinZeroN(upTotalChange_,numberObjects_);
611 CoinZeroN(downTotalChange_,numberObjects_);
612 CoinZeroN(upNumber_,numberObjects_);
613 CoinZeroN(downNumber_,numberObjects_);
614 }
615 }
616
617
618 //##############################################################################
619
OsiChooseStrong()620 OsiChooseStrong::OsiChooseStrong() :
621 OsiChooseVariable(),
622 shadowPriceMode_(0),
623 pseudoCosts_(),
624 results_(NULL),
625 numResults_(0)
626 {
627 }
628
OsiChooseStrong(const OsiSolverInterface * solver)629 OsiChooseStrong::OsiChooseStrong(const OsiSolverInterface * solver) :
630 OsiChooseVariable(solver),
631 shadowPriceMode_(0),
632 pseudoCosts_(),
633 results_(NULL),
634 numResults_(0)
635 {
636 // create useful arrays
637 pseudoCosts_.initialize(solver_->numberObjects());
638 }
639
OsiChooseStrong(const OsiChooseStrong & rhs)640 OsiChooseStrong::OsiChooseStrong(const OsiChooseStrong & rhs) :
641 OsiChooseVariable(rhs),
642 shadowPriceMode_(rhs.shadowPriceMode_),
643 pseudoCosts_(rhs.pseudoCosts_),
644 results_(NULL),
645 numResults_(0)
646 {
647 }
648
649 OsiChooseStrong &
operator =(const OsiChooseStrong & rhs)650 OsiChooseStrong::operator=(const OsiChooseStrong & rhs)
651 {
652 if (this != &rhs) {
653 OsiChooseVariable::operator=(rhs);
654 shadowPriceMode_ = rhs.shadowPriceMode_;
655 pseudoCosts_ = rhs.pseudoCosts_;
656 delete[] results_;
657 results_ = NULL;
658 numResults_ = 0;
659 }
660 return *this;
661 }
662
663
~OsiChooseStrong()664 OsiChooseStrong::~OsiChooseStrong ()
665 {
666 delete[] results_;
667 }
668
669 // Clone
670 OsiChooseVariable *
clone() const671 OsiChooseStrong::clone() const
672 {
673 return new OsiChooseStrong(*this);
674 }
675 #define MAXMIN_CRITERION 0.85
676 // Initialize
677 int
setupList(OsiBranchingInformation * info,bool initialize)678 OsiChooseStrong::setupList ( OsiBranchingInformation *info, bool initialize)
679 {
680 if (initialize) {
681 status_=-2;
682 delete [] goodSolution_;
683 bestObjectIndex_=-1;
684 numberStrongDone_=0;
685 numberStrongIterations_ = 0;
686 numberStrongFixed_ = 0;
687 goodSolution_ = NULL;
688 goodObjectiveValue_ = COIN_DBL_MAX;
689 }
690 numberOnList_=0;
691 numberUnsatisfied_=0;
692 int numberObjects = solver_->numberObjects();
693 if (numberObjects>pseudoCosts_.numberObjects()) {
694 // redo useful arrays
695 pseudoCosts_.initialize(numberObjects);
696 }
697 double check = -COIN_DBL_MAX;
698 int checkIndex=0;
699 int bestPriority=COIN_INT_MAX;
700 int maximumStrong= CoinMin(numberStrong_,numberObjects) ;
701 int putOther = numberObjects;
702 int i;
703 for (i=0;i<numberObjects;i++) {
704 list_[i]=-1;
705 useful_[i]=0.0;
706 }
707 OsiObject ** object = info->solver_->objects();
708 // Get average pseudo costs and see if pseudo shadow prices possible
709 int shadowPossible=shadowPriceMode_;
710 if (shadowPossible) {
711 for ( i=0;i<numberObjects;i++) {
712 if ( !object[i]->canHandleShadowPrices()) {
713 shadowPossible=0;
714 break;
715 }
716 }
717 if (shadowPossible) {
718 int numberRows = solver_->getNumRows();
719 const double * pi = info->pi_;
720 double sumPi=0.0;
721 for (i=0;i<numberRows;i++)
722 sumPi += fabs(pi[i]);
723 sumPi /= static_cast<double> (numberRows);
724 // and scale back
725 sumPi *= 0.01;
726 info->defaultDual_ = sumPi; // switch on
727 int numberColumns = solver_->getNumCols();
728 int size = CoinMax(numberColumns,2*numberRows);
729 info->usefulRegion_ = new double [size];
730 CoinZeroN(info->usefulRegion_,size);
731 info->indexRegion_ = new int [size];
732 }
733 }
734 double sumUp=0.0;
735 double numberUp=0.0;
736 double sumDown=0.0;
737 double numberDown=0.0;
738 const double* upTotalChange = pseudoCosts_.upTotalChange();
739 const double* downTotalChange = pseudoCosts_.downTotalChange();
740 const int* upNumber = pseudoCosts_.upNumber();
741 const int* downNumber = pseudoCosts_.downNumber();
742 const int numberBeforeTrusted = pseudoCosts_.numberBeforeTrusted();
743 for ( i=0;i<numberObjects;i++) {
744 sumUp += upTotalChange[i];
745 numberUp += upNumber[i];
746 sumDown += downTotalChange[i];
747 numberDown += downNumber[i];
748 }
749 double upMultiplier=(1.0+sumUp)/(1.0+numberUp);
750 double downMultiplier=(1.0+sumDown)/(1.0+numberDown);
751 // Say feasible
752 bool feasible = true;
753 #if 0
754 int pri[]={10,1000,10000};
755 int priCount[]={0,0,0};
756 #endif
757 for ( i=0;i<numberObjects;i++) {
758 int way;
759 double value = object[i]->infeasibility(info,way);
760 if (value>0.0) {
761 numberUnsatisfied_++;
762 if (value==COIN_DBL_MAX) {
763 // infeasible
764 feasible=false;
765 break;
766 }
767 int priorityLevel = object[i]->priority();
768 #if 0
769 for (int k=0;k<3;k++) {
770 if (priorityLevel==pri[k])
771 priCount[k]++;
772 }
773 #endif
774 // Better priority? Flush choices.
775 if (priorityLevel<bestPriority) {
776 for (int j=maximumStrong-1;j>=0;j--) {
777 if (list_[j]>=0) {
778 int iObject = list_[j];
779 list_[j]=-1;
780 useful_[j]=0.0;
781 list_[--putOther]=iObject;
782 }
783 }
784 maximumStrong = CoinMin(maximumStrong,putOther);
785 bestPriority = priorityLevel;
786 check=-COIN_DBL_MAX;
787 checkIndex=0;
788 }
789 if (priorityLevel==bestPriority) {
790 // Modify value
791 sumUp = upTotalChange[i]+1.0e-30;
792 numberUp = upNumber[i];
793 sumDown = downTotalChange[i]+1.0e-30;
794 numberDown = downNumber[i];
795 double upEstimate = object[i]->upEstimate();
796 double downEstimate = object[i]->downEstimate();
797 if (shadowPossible<2) {
798 upEstimate = numberUp ? ((upEstimate*sumUp)/numberUp) : (upEstimate*upMultiplier);
799 if (numberUp<numberBeforeTrusted)
800 upEstimate *= (numberBeforeTrusted+1.0)/(numberUp+1.0);
801 downEstimate = numberDown ? ((downEstimate*sumDown)/numberDown) : (downEstimate*downMultiplier);
802 if (numberDown<numberBeforeTrusted)
803 downEstimate *= (numberBeforeTrusted+1.0)/(numberDown+1.0);
804 } else {
805 // use shadow prices always
806 }
807 value = MAXMIN_CRITERION*CoinMin(upEstimate,downEstimate) + (1.0-MAXMIN_CRITERION)*CoinMax(upEstimate,downEstimate);
808 if (value>check) {
809 //add to list
810 int iObject = list_[checkIndex];
811 if (iObject>=0) {
812 assert (list_[putOther-1]<0);
813 list_[--putOther]=iObject; // to end
814 }
815 list_[checkIndex]=i;
816 assert (checkIndex<putOther);
817 useful_[checkIndex]=value;
818 // find worst
819 check=COIN_DBL_MAX;
820 maximumStrong = CoinMin(maximumStrong,putOther);
821 for (int j=0;j<maximumStrong;j++) {
822 if (list_[j]>=0) {
823 if (useful_[j]<check) {
824 check=useful_[j];
825 checkIndex=j;
826 }
827 } else {
828 check=0.0;
829 checkIndex = j;
830 break;
831 }
832 }
833 } else {
834 // to end
835 assert (list_[putOther-1]<0);
836 list_[--putOther]=i;
837 maximumStrong = CoinMin(maximumStrong,putOther);
838 }
839 } else {
840 // worse priority
841 // to end
842 assert (list_[putOther-1]<0);
843 list_[--putOther]=i;
844 maximumStrong = CoinMin(maximumStrong,putOther);
845 }
846 }
847 }
848 #if 0
849 printf("%d at %d, %d at %d and %d at %d\n",priCount[0],pri[0],
850 priCount[1],pri[1],priCount[2],pri[2]);
851 #endif
852 // Get list
853 numberOnList_=0;
854 if (feasible) {
855 for (i=0;i<CoinMin(maximumStrong,putOther);i++) {
856 if (list_[i]>=0) {
857 list_[numberOnList_]=list_[i];
858 useful_[numberOnList_++]=-useful_[i];
859 }
860 }
861 if (numberOnList_) {
862 // Sort
863 CoinSort_2(useful_,useful_+numberOnList_,list_);
864 // move others
865 i = numberOnList_;
866 for (;putOther<numberObjects;putOther++)
867 list_[i++]=list_[putOther];
868 assert (i==numberUnsatisfied_);
869 if (!numberStrong_)
870 numberOnList_=0;
871 }
872 } else {
873 // not feasible
874 numberUnsatisfied_=-1;
875 }
876 // Get rid of any shadow prices info
877 info->defaultDual_ = -1.0; // switch off
878 delete [] info->usefulRegion_;
879 delete [] info->indexRegion_;
880 return numberUnsatisfied_;
881 }
882
883 void
resetResults(int num)884 OsiChooseStrong::resetResults(int num)
885 {
886 delete[] results_;
887 numResults_ = 0;
888 results_ = new OsiHotInfo[num];
889 }
890
891 /* Choose a variable
892 Returns -
893 -1 Node is infeasible
894 0 Normal termination - we have a candidate
895 1 All looks satisfied - no candidate
896 2 We can change the bound on a variable - but we also have a strong branching candidate
897 3 We can change the bound on a variable - but we have a non-strong branching candidate
898 4 We can change the bound on a variable - no other candidates
899 We can pick up branch from whichObject() and whichWay()
900 We can pick up a forced branch (can change bound) from whichForcedObject() and whichForcedWay()
901 If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
902 */
903 int
chooseVariable(OsiSolverInterface * solver,OsiBranchingInformation * info,bool fixVariables)904 OsiChooseStrong::chooseVariable( OsiSolverInterface * solver, OsiBranchingInformation *info, bool fixVariables)
905 {
906 if (numberUnsatisfied_) {
907 const double* upTotalChange = pseudoCosts_.upTotalChange();
908 const double* downTotalChange = pseudoCosts_.downTotalChange();
909 const int* upNumber = pseudoCosts_.upNumber();
910 const int* downNumber = pseudoCosts_.downNumber();
911 int numberBeforeTrusted = pseudoCosts_.numberBeforeTrusted();
912 // Somehow we can get here with it 0 !
913 if (!numberBeforeTrusted) {
914 numberBeforeTrusted=5;
915 pseudoCosts_.setNumberBeforeTrusted(numberBeforeTrusted);
916 }
917
918 int numberLeft = CoinMin(numberStrong_-numberStrongDone_,numberUnsatisfied_);
919 int numberToDo=0;
920 resetResults(numberLeft);
921 int returnCode=0;
922 bestObjectIndex_ = -1;
923 bestWhichWay_ = -1;
924 firstForcedObjectIndex_ = -1;
925 firstForcedWhichWay_ =-1;
926 double bestTrusted=-COIN_DBL_MAX;
927 for (int i=0;i<numberLeft;i++) {
928 int iObject = list_[i];
929 if (upNumber[iObject]<numberBeforeTrusted||downNumber[iObject]<numberBeforeTrusted) {
930 results_[numberToDo++] = OsiHotInfo(solver, info,
931 solver->objects(), iObject);
932 } else {
933 const OsiObject * obj = solver->object(iObject);
934 double upEstimate = (upTotalChange[iObject]*obj->upEstimate())/upNumber[iObject];
935 double downEstimate = (downTotalChange[iObject]*obj->downEstimate())/downNumber[iObject];
936 double value = MAXMIN_CRITERION*CoinMin(upEstimate,downEstimate) + (1.0-MAXMIN_CRITERION)*CoinMax(upEstimate,downEstimate);
937 if (value > bestTrusted) {
938 bestObjectIndex_=iObject;
939 bestWhichWay_ = upEstimate>downEstimate ? 0 : 1;
940 bestTrusted = value;
941 }
942 }
943 }
944 int numberFixed=0;
945 if (numberToDo) {
946 returnCode = doStrongBranching(solver, info, numberToDo, 1);
947 if (returnCode>=0&&returnCode<=2) {
948 if (returnCode) {
949 returnCode=4;
950 if (bestObjectIndex_>=0)
951 returnCode=3;
952 }
953 for (int i=0;i<numResults_;i++) {
954 int iObject = results_[i].whichObject();
955 double upEstimate;
956 if (results_[i].upStatus()!=1) {
957 assert (results_[i].upStatus()>=0);
958 upEstimate = results_[i].upChange();
959 } else {
960 // infeasible - just say expensive
961 if (info->cutoff_<1.0e50)
962 upEstimate = 2.0*(info->cutoff_-info->objectiveValue_);
963 else
964 upEstimate = 2.0*fabs(info->objectiveValue_);
965 if (firstForcedObjectIndex_ <0) {
966 firstForcedObjectIndex_ = iObject;
967 firstForcedWhichWay_ =0;
968 }
969 numberFixed++;
970 if (fixVariables) {
971 const OsiObject * obj = solver->object(iObject);
972 OsiBranchingObject * branch = obj->createBranch(solver,info,0);
973 branch->branch(solver);
974 delete branch;
975 }
976 }
977 double downEstimate;
978 if (results_[i].downStatus()!=1) {
979 assert (results_[i].downStatus()>=0);
980 downEstimate = results_[i].downChange();
981 } else {
982 // infeasible - just say expensive
983 if (info->cutoff_<1.0e50)
984 downEstimate = 2.0*(info->cutoff_-info->objectiveValue_);
985 else
986 downEstimate = 2.0*fabs(info->objectiveValue_);
987 if (firstForcedObjectIndex_ <0) {
988 firstForcedObjectIndex_ = iObject;
989 firstForcedWhichWay_ =1;
990 }
991 numberFixed++;
992 if (fixVariables) {
993 const OsiObject * obj = solver->object(iObject);
994 OsiBranchingObject * branch = obj->createBranch(solver,info,1);
995 branch->branch(solver);
996 delete branch;
997 }
998 }
999 double value = MAXMIN_CRITERION*CoinMin(upEstimate,downEstimate) + (1.0-MAXMIN_CRITERION)*CoinMax(upEstimate,downEstimate);
1000 if (value>bestTrusted) {
1001 bestTrusted = value;
1002 bestObjectIndex_ = iObject;
1003 bestWhichWay_ = upEstimate>downEstimate ? 0 : 1;
1004 // but override if there is a preferred way
1005 const OsiObject * obj = solver->object(iObject);
1006 if (obj->preferredWay()>=0&&obj->infeasibility())
1007 bestWhichWay_ = obj->preferredWay();
1008 if (returnCode)
1009 returnCode=2;
1010 }
1011 }
1012 } else if (returnCode==3) {
1013 // max time - just choose one
1014 bestObjectIndex_ = list_[0];
1015 bestWhichWay_ = 0;
1016 returnCode=0;
1017 }
1018 } else {
1019 bestObjectIndex_=list_[0];
1020 }
1021 if ( bestObjectIndex_ >=0 ) {
1022 OsiObject * obj = solver->objects()[bestObjectIndex_];
1023 obj->setWhichWay( bestWhichWay_);
1024 }
1025 if (numberFixed==numberUnsatisfied_&&numberFixed)
1026 returnCode=4;
1027 return returnCode;
1028 } else {
1029 return 1;
1030 }
1031 }
1032 // Given a candidate fill in useful information e.g. estimates
1033 void
updateInformation(const OsiBranchingInformation * info,int branch,OsiHotInfo * hotInfo)1034 OsiPseudoCosts::updateInformation(const OsiBranchingInformation *info,
1035 int branch, OsiHotInfo * hotInfo)
1036 {
1037 int index = hotInfo->whichObject();
1038 assert (index<info->solver_->numberObjects());
1039 const OsiObject * object = info->solver_->object(index);
1040 assert (object->upEstimate()>0.0&&object->downEstimate()>0.0);
1041 assert (branch<2);
1042 if (branch) {
1043 if (hotInfo->upStatus()!=1) {
1044 assert (hotInfo->upStatus()>=0);
1045 upTotalChange_[index] += hotInfo->upChange()/object->upEstimate();
1046 upNumber_[index]++;
1047 } else {
1048 #if 0
1049 // infeasible - just say expensive
1050 if (info->cutoff_<1.0e50)
1051 upTotalChange_[index] += 2.0*(info->cutoff_-info->objectiveValue_)/object->upEstimate();
1052 else
1053 upTotalChange_[index] += 2.0*fabs(info->objectiveValue_)/object->upEstimate();
1054 #endif
1055 }
1056 } else {
1057 if (hotInfo->downStatus()!=1) {
1058 assert (hotInfo->downStatus()>=0);
1059 downTotalChange_[index] += hotInfo->downChange()/object->downEstimate();
1060 downNumber_[index]++;
1061 } else {
1062 #if 0
1063 // infeasible - just say expensive
1064 if (info->cutoff_<1.0e50)
1065 downTotalChange_[index] += 2.0*(info->cutoff_-info->objectiveValue_)/object->downEstimate();
1066 else
1067 downTotalChange_[index] += 2.0*fabs(info->objectiveValue_)/object->downEstimate();
1068 #endif
1069 }
1070 }
1071 }
1072 #if 1
1073 // Given a branch fill in useful information e.g. estimates
1074 void
updateInformation(int index,int branch,double changeInObjective,double changeInValue,int status)1075 OsiPseudoCosts::updateInformation(int index, int branch,
1076 double changeInObjective,
1077 double changeInValue,
1078 int status)
1079 {
1080 //assert (index<solver_->numberObjects());
1081 assert (branch<2);
1082 assert (changeInValue>0.0);
1083 assert (branch<2);
1084 if (branch) {
1085 if (status!=1) {
1086 assert (status>=0);
1087 upTotalChange_[index] += changeInObjective/changeInValue;
1088 upNumber_[index]++;
1089 }
1090 } else {
1091 if (status!=1) {
1092 assert (status>=0);
1093 downTotalChange_[index] += changeInObjective/changeInValue;
1094 downNumber_[index]++;
1095 }
1096 }
1097 }
1098 #endif
1099
OsiHotInfo()1100 OsiHotInfo::OsiHotInfo() :
1101 originalObjectiveValue_(COIN_DBL_MAX),
1102 changes_(NULL),
1103 iterationCounts_(NULL),
1104 statuses_(NULL),
1105 branchingObject_(NULL),
1106 whichObject_(-1)
1107 {
1108 }
1109
OsiHotInfo(OsiSolverInterface * solver,const OsiBranchingInformation * info,const OsiObject * const * objects,int whichObject)1110 OsiHotInfo::OsiHotInfo(OsiSolverInterface * solver,
1111 const OsiBranchingInformation * info,
1112 const OsiObject * const * objects,
1113 int whichObject) :
1114 originalObjectiveValue_(COIN_DBL_MAX),
1115 whichObject_(whichObject)
1116 {
1117 originalObjectiveValue_ = info->objectiveValue_;
1118 const OsiObject * object = objects[whichObject_];
1119 // create object - "down" first
1120 branchingObject_ = object->createBranch(solver,info,0);
1121 // create arrays
1122 int numberBranches = branchingObject_->numberBranches();
1123 changes_ = new double [numberBranches];
1124 iterationCounts_ = new int [numberBranches];
1125 statuses_ = new int [numberBranches];
1126 CoinZeroN(changes_,numberBranches);
1127 CoinZeroN(iterationCounts_,numberBranches);
1128 CoinFillN(statuses_,numberBranches,-1);
1129 }
1130
OsiHotInfo(const OsiHotInfo & rhs)1131 OsiHotInfo::OsiHotInfo(const OsiHotInfo & rhs)
1132 {
1133 originalObjectiveValue_ = rhs.originalObjectiveValue_;
1134 whichObject_ = rhs.whichObject_;
1135 if (rhs.branchingObject_) {
1136 branchingObject_ = rhs.branchingObject_->clone();
1137 int numberBranches = branchingObject_->numberBranches();
1138 changes_ = CoinCopyOfArray(rhs.changes_,numberBranches);
1139 iterationCounts_ = CoinCopyOfArray(rhs.iterationCounts_,numberBranches);
1140 statuses_ = CoinCopyOfArray(rhs.statuses_,numberBranches);
1141 } else {
1142 branchingObject_ = NULL;
1143 changes_ = NULL;
1144 iterationCounts_ = NULL;
1145 statuses_ = NULL;
1146 }
1147 }
1148
1149 OsiHotInfo &
operator =(const OsiHotInfo & rhs)1150 OsiHotInfo::operator=(const OsiHotInfo & rhs)
1151 {
1152 if (this != &rhs) {
1153 delete branchingObject_;
1154 delete [] changes_;
1155 delete [] iterationCounts_;
1156 delete [] statuses_;
1157 originalObjectiveValue_ = rhs.originalObjectiveValue_;
1158 whichObject_ = rhs.whichObject_;
1159 if (rhs.branchingObject_) {
1160 branchingObject_ = rhs.branchingObject_->clone();
1161 int numberBranches = branchingObject_->numberBranches();
1162 changes_ = CoinCopyOfArray(rhs.changes_,numberBranches);
1163 iterationCounts_ = CoinCopyOfArray(rhs.iterationCounts_,numberBranches);
1164 statuses_ = CoinCopyOfArray(rhs.statuses_,numberBranches);
1165 } else {
1166 branchingObject_ = NULL;
1167 changes_ = NULL;
1168 iterationCounts_ = NULL;
1169 statuses_ = NULL;
1170 }
1171 }
1172 return *this;
1173 }
1174
1175
~OsiHotInfo()1176 OsiHotInfo::~OsiHotInfo ()
1177 {
1178 delete branchingObject_;
1179 delete [] changes_;
1180 delete [] iterationCounts_;
1181 delete [] statuses_;
1182 }
1183
1184 // Clone
1185 OsiHotInfo *
clone() const1186 OsiHotInfo::clone() const
1187 {
1188 return new OsiHotInfo(*this);
1189 }
1190 /* Fill in useful information after strong branch
1191 */
updateInformation(const OsiSolverInterface * solver,const OsiBranchingInformation * info,OsiChooseVariable * choose)1192 int OsiHotInfo::updateInformation( const OsiSolverInterface * solver, const OsiBranchingInformation * info,
1193 OsiChooseVariable * choose)
1194 {
1195 int iBranch = branchingObject_->branchIndex()-1;
1196 assert (iBranch>=0&&iBranch<branchingObject_->numberBranches());
1197 iterationCounts_[iBranch] += solver->getIterationCount();
1198 int status;
1199 if (solver->isProvenOptimal())
1200 status=0; // optimal
1201 else if (solver->isIterationLimitReached()
1202 &&!solver->isDualObjectiveLimitReached())
1203 status=2; // unknown
1204 else
1205 status=1; // infeasible
1206 // Could do something different if we can't trust
1207 double newObjectiveValue = solver->getObjSense()*solver->getObjValue();
1208 changes_[iBranch] =CoinMax(0.0,newObjectiveValue-originalObjectiveValue_);
1209 // we might have got here by primal
1210 if (choose->trustStrongForBound()) {
1211 if (!status&&newObjectiveValue>=info->cutoff_) {
1212 status=1; // infeasible
1213 changes_[iBranch] = 1.0e100;
1214 }
1215 }
1216 statuses_[iBranch] = status;
1217 if (!status&&choose->trustStrongForSolution()&&newObjectiveValue<choose->goodObjectiveValue()) {
1218 // check if solution
1219 const OsiSolverInterface * saveSolver = info->solver_;
1220 info->solver_=solver;
1221 const double * saveLower = info->lower_;
1222 info->lower_ = solver->getColLower();
1223 const double * saveUpper = info->upper_;
1224 info->upper_ = solver->getColUpper();
1225 // also need to make sure bounds OK as may not be info solver
1226 #if 0
1227 if (saveSolver->getMatrixByCol()) {
1228 const CoinBigIndex * columnStart = info->columnStart_;
1229 assert (saveSolver->getMatrixByCol()->getVectorStarts()==columnStart);
1230 }
1231 #endif
1232 if (choose->feasibleSolution(info,solver->getColSolution(),solver->numberObjects(),
1233 const_cast<const OsiObject **> (solver->objects()))) {
1234 // put solution somewhere
1235 choose->saveSolution(solver);
1236 status=3;
1237 }
1238 info->solver_=saveSolver;
1239 info->lower_ = saveLower;
1240 info->upper_ = saveUpper;
1241 }
1242 // Now update - possible strong branching info
1243 choose->updateInformation( info,iBranch,this);
1244 return status;
1245 }
1246