Changeset 965 for branches/dynamicbranching/dynamicbranching.cpp
 Timestamp:
 Jun 5, 2008 11:13:18 AM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/dynamicbranching/dynamicbranching.cpp
r963 r965 1 #include "CoinTime.hpp" 2 #include "OsiClpSolverInterface.hpp" 3 4 5 // below needed for pathetic branch and bound code 6 #include <vector> 7 #include <map> 8 9 // Trivial class for Branch and Bound 10 11 class DBNodeSimple { 12 13 public: 14 15 // Default Constructor 16 DBNodeSimple (); 17 18 // Constructor from current state (and list of integers) 19 // Also chooses branching variable (if none set to 1) 20 DBNodeSimple (OsiSolverInterface &model, 21 int numberIntegers, int * integer, 22 CoinWarmStart * basis); 23 void gutsOfConstructor (OsiSolverInterface &model, 24 int numberIntegers, int * integer, 25 CoinWarmStart * basis); 26 // Copy constructor 27 DBNodeSimple ( const DBNodeSimple &); 28 29 // Assignment operator 30 DBNodeSimple & operator=( const DBNodeSimple& rhs); 31 32 // Destructor 33 ~DBNodeSimple (); 34 // Work of destructor 35 void gutsOfDestructor(); 36 // Extension  true if other extension of this 37 bool extension(const DBNodeSimple & other, 38 const double * originalLower, 39 const double * originalUpper) const; 40 inline void incrementDescendants() 41 { descendants_++;} 42 // Public data 43 // Basis (should use tree, but not as wasteful as bounds!) 44 CoinWarmStart * basis_; 45 // Objective value (COIN_DBL_MAX) if spare node 46 double objectiveValue_; 47 // Branching variable (0 is first integer) 48 int variable_; 49 // Way to branch  1 down (first), 1 up, 2 down (second), 2 up (second) 50 int way_; 51 // Number of integers (for length of arrays) 52 int numberIntegers_; 53 // Current value 54 double value_; 55 // Number of descendant nodes (so 2 is in interior) 56 int descendants_; 57 // Parent 58 int parent_; 59 // Previous in chain 60 int previous_; 61 // Next in chain 62 int next_; 63 // Now I must use tree 64 // Bounds stored in full (for integers) 65 int * lower_; 66 int * upper_; 67 }; 68 69 70 DBNodeSimple::DBNodeSimple() : 71 basis_(NULL), 72 objectiveValue_(COIN_DBL_MAX), 73 variable_(100), 74 way_(1), 75 numberIntegers_(0), 76 value_(0.5), 77 descendants_(1), 78 parent_(1), 79 previous_(1), 80 next_(1), 81 lower_(NULL), 82 upper_(NULL) 83 { 84 } 85 DBNodeSimple::DBNodeSimple(OsiSolverInterface & model, 86 int numberIntegers, int * integer,CoinWarmStart * basis) 87 { 88 gutsOfConstructor(model,numberIntegers,integer,basis); 89 } 90 void 91 DBNodeSimple::gutsOfConstructor(OsiSolverInterface & model, 92 int numberIntegers, int * integer,CoinWarmStart * basis) 93 { 94 basis_ = basis; 95 variable_=1; 96 way_=1; 97 numberIntegers_=numberIntegers; 98 value_=0.0; 99 descendants_ = 0; 100 parent_ = 1; 101 previous_ = 1; 102 next_ = 1; 103 if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { 104 objectiveValue_ = model.getObjSense()*model.getObjValue(); 105 } else { 106 objectiveValue_ = 1.0e100; 107 lower_ = NULL; 108 upper_ = NULL; 109 return; // node cutoff 110 } 111 lower_ = new int [numberIntegers_]; 112 upper_ = new int [numberIntegers_]; 113 assert (upper_!=NULL); 114 const double * lower = model.getColLower(); 115 const double * upper = model.getColUpper(); 116 const double * solution = model.getColSolution(); 117 int i; 118 // Hard coded integer tolerance 119 #define INTEGER_TOLERANCE 1.0e6 120 ///////// Start of Strong branching code  can be ignored 121 // Number of strong branching candidates 122 #define STRONG_BRANCHING 5 123 #ifdef STRONG_BRANCHING 124 double upMovement[STRONG_BRANCHING]; 125 double downMovement[STRONG_BRANCHING]; 126 double solutionValue[STRONG_BRANCHING]; 127 int chosen[STRONG_BRANCHING]; 128 int iSmallest=0; 129 // initialize distance from integer 130 for (i=0;i<STRONG_BRANCHING;i++) { 131 upMovement[i]=0.0; 132 chosen[i]=1; 133 } 134 variable_=1; 135 // This has hard coded integer tolerance 136 double mostAway=INTEGER_TOLERANCE; 137 int numberAway=0; 138 for (i=0;i<numberIntegers;i++) { 139 int iColumn = integer[i]; 140 lower_[i]=(int)lower[iColumn]; 141 upper_[i]=(int)upper[iColumn]; 142 double value = solution[iColumn]; 143 value = max(value,(double) lower_[i]); 144 value = min(value,(double) upper_[i]); 145 double nearest = floor(value+0.5); 146 if (fabs(valuenearest)>INTEGER_TOLERANCE) 147 numberAway++; 148 if (fabs(valuenearest)>mostAway) { 149 double away = fabs(valuenearest); 150 if (away>upMovement[iSmallest]) { 151 //add to list 152 upMovement[iSmallest]=away; 153 solutionValue[iSmallest]=value; 154 chosen[iSmallest]=i; 155 int j; 156 iSmallest=1; 157 double smallest = 1.0; 158 for (j=0;j<STRONG_BRANCHING;j++) { 159 if (upMovement[j]<smallest) { 160 smallest=upMovement[j]; 161 iSmallest=j; 162 } 163 } 164 } 165 } 166 } 167 int numberStrong=0; 168 for (i=0;i<STRONG_BRANCHING;i++) { 169 if (chosen[i]>=0) { 170 numberStrong ++; 171 variable_ = chosen[i]; 172 } 173 } 174 // out strong branching if bit set 175 OsiClpSolverInterface* clp = 176 dynamic_cast<OsiClpSolverInterface*>(&model); 177 if (clp&&(clp>specialOptions()&16)!=0&&numberStrong>1) { 178 int j; 179 int iBest=1; 180 double best = 0.0; 181 for (j=0;j<STRONG_BRANCHING;j++) { 182 if (upMovement[j]>best) { 183 best=upMovement[j]; 184 iBest=j; 185 } 186 } 187 numberStrong=1; 188 variable_=chosen[iBest]; 189 } 190 if (numberStrong==1) { 191 // just one  makes it easy 192 int iColumn = integer[variable_]; 193 double value = solution[iColumn]; 194 value = max(value,(double) lower_[variable_]); 195 value = min(value,(double) upper_[variable_]); 196 double nearest = floor(value+0.5); 197 value_=value; 198 if (value<=nearest) 199 way_=1; // up 200 else 201 way_=1; // down 202 } else if (numberStrong) { 203 // more than one  choose 204 bool chooseOne=true; 205 model.markHotStart(); 206 for (i=0;i<STRONG_BRANCHING;i++) { 207 int iInt = chosen[i]; 208 if (iInt>=0) { 209 int iColumn = integer[iInt]; 210 double value = solutionValue[i]; // value of variable in original 211 double objectiveChange; 212 value = max(value,(double) lower_[iInt]); 213 value = min(value,(double) upper_[iInt]); 214 215 // try down 216 217 model.setColUpper(iColumn,floor(value)); 218 model.solveFromHotStart(); 219 model.setColUpper(iColumn,upper_[iInt]); 220 if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { 221 objectiveChange = model.getObjSense()*model.getObjValue() 222  objectiveValue_; 223 } else { 224 objectiveChange = 1.0e100; 225 } 226 assert (objectiveChange>1.0e5); 227 objectiveChange = CoinMax(objectiveChange,0.0); 228 downMovement[i]=objectiveChange; 229 230 // try up 231 232 model.setColLower(iColumn,ceil(value)); 233 model.solveFromHotStart(); 234 model.setColLower(iColumn,lower_[iInt]); 235 if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { 236 objectiveChange = model.getObjSense()*model.getObjValue() 237  objectiveValue_; 238 } else { 239 objectiveChange = 1.0e100; 240 } 241 assert (objectiveChange>1.0e5); 242 objectiveChange = CoinMax(objectiveChange,0.0); 243 upMovement[i]=objectiveChange; 244 245 /* Possibilities are: 246 Both sides feasible  store 247 Neither side feasible  set objective high and exit 248 One side feasible  change bounds and resolve 249 */ 250 bool solveAgain=false; 251 if (upMovement[i]<1.0e100) { 252 if(downMovement[i]<1.0e100) { 253 // feasible  no action 254 } else { 255 // up feasible, down infeasible 256 solveAgain = true; 257 model.setColLower(iColumn,ceil(value)); 258 } 259 } else { 260 if(downMovement[i]<1.0e100) { 261 // down feasible, up infeasible 262 solveAgain = true; 263 model.setColUpper(iColumn,floor(value)); 264 } else { 265 // neither side feasible 266 objectiveValue_=1.0e100; 267 chooseOne=false; 268 break; 269 } 270 } 271 if (solveAgain) { 272 // need to solve problem again  signal this 273 variable_ = numberIntegers; 274 chooseOne=false; 275 break; 276 } 277 } 278 } 279 if (chooseOne) { 280 // choose the one that makes most difference both ways 281 double best = 1.0; 282 double best2 = 1.0; 283 for (i=0;i<STRONG_BRANCHING;i++) { 284 int iInt = chosen[i]; 285 if (iInt>=0) { 286 //std::cout<<"Strong branching on " 287 // <<i<<""<<iInt<<" down "<<downMovement[i] 288 // <<" up "<<upMovement[i] 289 // <<" value "<<solutionValue[i] 290 // <<std::endl; 291 bool better = false; 292 if (min(upMovement[i],downMovement[i])>best) { 293 // smaller is better 294 better=true; 295 } else if (min(upMovement[i],downMovement[i])>best1.0e5) { 296 if (max(upMovement[i],downMovement[i])>best2+1.0e5) { 297 // smaller is about same, but larger is better 298 better=true; 299 } 300 } 301 if (better) { 302 best = min(upMovement[i],downMovement[i]); 303 best2 = max(upMovement[i],downMovement[i]); 304 variable_ = iInt; 305 double value = solutionValue[i]; 306 value = max(value,(double) lower_[variable_]); 307 value = min(value,(double) upper_[variable_]); 308 value_=value; 309 if (upMovement[i]<=downMovement[i]) 310 way_=1; // up 311 else 312 way_=1; // down 313 } 314 } 315 } 316 } 317 // Delete the snapshot 318 model.unmarkHotStart(); 319 } 320 ////// End of Strong branching 321 #else 322 variable_=1; 323 // This has hard coded integer tolerance 324 double mostAway=INTEGER_TOLERANCE; 325 int numberAway=0; 326 for (i=0;i<numberIntegers;i++) { 327 int iColumn = integer[i]; 328 lower_[i]=(int)lower[iColumn]; 329 upper_[i]=(int)upper[iColumn]; 330 double value = solution[iColumn]; 331 value = max(value,(double) lower_[i]); 332 value = min(value,(double) upper_[i]); 333 double nearest = floor(value+0.5); 334 if (fabs(valuenearest)>INTEGER_TOLERANCE) 335 numberAway++; 336 if (fabs(valuenearest)>mostAway) { 337 mostAway=fabs(valuenearest); 338 variable_=i; 339 value_=value; 340 if (value<=nearest) 341 way_=1; // up 342 else 343 way_=1; // down 344 } 345 } 346 #endif 347 } 348 349 DBNodeSimple::DBNodeSimple(const DBNodeSimple & rhs) 350 { 351 if (rhs.basis_) 352 basis_=rhs.basis_>clone(); 353 else 354 basis_ = NULL; 355 objectiveValue_=rhs.objectiveValue_; 356 variable_=rhs.variable_; 357 way_=rhs.way_; 358 numberIntegers_=rhs.numberIntegers_; 359 value_=rhs.value_; 360 descendants_ = rhs.descendants_; 361 parent_ = rhs.parent_; 362 previous_ = rhs.previous_; 363 next_ = rhs.next_; 364 lower_=NULL; 365 upper_=NULL; 366 if (rhs.lower_!=NULL) { 367 lower_ = new int [numberIntegers_]; 368 upper_ = new int [numberIntegers_]; 369 assert (upper_!=NULL); 370 memcpy(lower_,rhs.lower_,numberIntegers_*sizeof(int)); 371 memcpy(upper_,rhs.upper_,numberIntegers_*sizeof(int)); 372 } 373 } 374 375 DBNodeSimple & 376 DBNodeSimple::operator=(const DBNodeSimple & rhs) 377 { 378 if (this != &rhs) { 379 gutsOfDestructor(); 380 if (rhs.basis_) 381 basis_=rhs.basis_>clone(); 382 objectiveValue_=rhs.objectiveValue_; 383 variable_=rhs.variable_; 384 way_=rhs.way_; 385 numberIntegers_=rhs.numberIntegers_; 386 value_=rhs.value_; 387 descendants_ = rhs.descendants_; 388 parent_ = rhs.parent_; 389 previous_ = rhs.previous_; 390 next_ = rhs.next_; 391 if (rhs.lower_!=NULL) { 392 lower_ = new int [numberIntegers_]; 393 upper_ = new int [numberIntegers_]; 394 assert (upper_!=NULL); 395 memcpy(lower_,rhs.lower_,numberIntegers_*sizeof(int)); 396 memcpy(upper_,rhs.upper_,numberIntegers_*sizeof(int)); 397 } 398 } 399 return *this; 400 } 401 402 403 DBNodeSimple::~DBNodeSimple () 404 { 405 gutsOfDestructor(); 406 } 407 // Work of destructor 408 void 409 DBNodeSimple::gutsOfDestructor() 410 { 411 delete [] lower_; 412 delete [] upper_; 413 delete basis_; 414 lower_ = NULL; 415 upper_ = NULL; 416 basis_ = NULL; 417 objectiveValue_ = COIN_DBL_MAX; 418 } 419 // Extension  true if other extension of this 420 bool 421 DBNodeSimple::extension(const DBNodeSimple & other, 422 const double * originalLower, 423 const double * originalUpper) const 424 { 425 bool ok=true; 426 for (int i=0;i<numberIntegers_;i++) { 427 if (upper_[i]<originalUpper[i] 428 lower_[i]>originalLower[i]) { 429 if (other.upper_[i]>upper_[i] 430 other.lower_[i]<lower_[i]) { 431 ok=false; 432 break; 433 } 434 } 435 } 436 return ok; 437 } 438 439 #include <vector> 440 #define FUNNY_BRANCHING 1 441 #define FUNNY_TREE 442 #ifndef FUNNY_TREE 443 // Vector of DBNodeSimples 444 typedef std::vector<DBNodeSimple> DBVectorNode; 445 #else 446 // Must code up by hand 447 class DBVectorNode { 448 449 public: 450 451 // Default Constructor 452 DBVectorNode (); 453 454 // Copy constructor 455 DBVectorNode ( const DBVectorNode &); 456 457 // Assignment operator 458 DBVectorNode & operator=( const DBVectorNode& rhs); 459 460 // Destructor 461 ~DBVectorNode (); 462 // Size 463 inline int size() const 464 { return size_sizeDeferred_;} 465 // Push 466 void push_back(const DBNodeSimple & node); 467 // Last one in (or other criterion) 468 DBNodeSimple back() const; 469 // Get rid of last one 470 void pop_back(); 471 // Works out best one 472 int best() const; 473 474 // Public data 475 // Maximum size 476 int maximumSize_; 477 // Current size 478 int size_; 479 // Number still hanging around 480 int sizeDeferred_; 481 // First spare 482 int firstSpare_; 483 // First 484 int first_; 485 // Last 486 int last_; 487 // Chosen one 488 mutable int chosen_; 489 // Nodes 490 DBNodeSimple * nodes_; 491 }; 492 493 494 DBVectorNode::DBVectorNode() : 495 maximumSize_(10), 496 size_(0), 497 sizeDeferred_(0), 498 firstSpare_(0), 499 first_(1), 500 last_(1) 501 { 502 nodes_ = new DBNodeSimple[maximumSize_]; 503 for (int i=0;i<maximumSize_;i++) { 504 nodes_[i].previous_=i1; 505 nodes_[i].next_=i+1; 506 } 507 } 508 509 DBVectorNode::DBVectorNode(const DBVectorNode & rhs) 510 { 511 maximumSize_ = rhs.maximumSize_; 512 size_ = rhs.size_; 513 sizeDeferred_ = rhs.sizeDeferred_; 514 firstSpare_ = rhs.firstSpare_; 515 first_ = rhs.first_; 516 last_ = rhs.last_; 517 nodes_ = new DBNodeSimple[maximumSize_]; 518 for (int i=0;i<maximumSize_;i++) { 519 nodes_[i] = rhs.nodes_[i]; 520 } 521 } 522 523 DBVectorNode & 524 DBVectorNode::operator=(const DBVectorNode & rhs) 525 { 526 if (this != &rhs) { 527 delete [] nodes_; 528 maximumSize_ = rhs.maximumSize_; 529 size_ = rhs.size_; 530 sizeDeferred_ = rhs.sizeDeferred_; 531 firstSpare_ = rhs.firstSpare_; 532 first_ = rhs.first_; 533 last_ = rhs.last_; 534 nodes_ = new DBNodeSimple[maximumSize_]; 535 for (int i=0;i<maximumSize_;i++) { 536 nodes_[i] = rhs.nodes_[i]; 537 } 538 } 539 return *this; 540 } 541 542 543 DBVectorNode::~DBVectorNode () 544 { 545 delete [] nodes_; 546 } 547 // Push 548 void 549 DBVectorNode::push_back(const DBNodeSimple & node) 550 { 551 if (size_==maximumSize_) { 552 assert (firstSpare_==size_); 553 maximumSize_ = (maximumSize_*3)+10; 554 DBNodeSimple * temp = new DBNodeSimple[maximumSize_]; 555 int i; 556 for (i=0;i<size_;i++) { 557 temp[i]=nodes_[i]; 558 } 559 delete [] nodes_; 560 nodes_ = temp; 561 //firstSpare_=size_; 562 int last = 1; 563 for ( i=size_;i<maximumSize_;i++) { 564 nodes_[i].previous_=last; 565 nodes_[i].next_=i+1; 566 last = i; 567 } 568 } 569 assert (firstSpare_<maximumSize_); 570 assert (nodes_[firstSpare_].previous_<0); 571 int next = nodes_[firstSpare_].next_; 572 nodes_[firstSpare_]=node; 573 if (last_>=0) { 574 assert (nodes_[last_].next_==1); 575 nodes_[last_].next_=firstSpare_; 576 } 577 nodes_[firstSpare_].previous_=last_; 578 nodes_[firstSpare_].next_=1; 579 if (last_==1) { 580 assert (first_==1); 581 first_ = firstSpare_; 582 } 583 last_=firstSpare_; 584 if (next>=0&&next<maximumSize_) { 585 firstSpare_ = next; 586 nodes_[firstSpare_].previous_=1; 587 } else { 588 firstSpare_=maximumSize_; 589 } 590 chosen_ = 1; 591 //best(); 592 size_++; 593 assert (node.descendants_<=2); 594 if (node.descendants_==2) 595 sizeDeferred_++; 596 } 597 // Works out best one 598 int 599 DBVectorNode::best() const 600 { 601 // can modify 602 chosen_=1; 603 if (chosen_<0) { 604 chosen_=last_; 605 #if FUNNY_BRANCHING 606 while (nodes_[chosen_].descendants_==2) { 607 chosen_ = nodes_[chosen_].previous_; 608 assert (chosen_>=0); 609 } 610 #endif 611 } 612 return chosen_; 613 } 614 // Last one in (or other criterion) 615 DBNodeSimple 616 DBVectorNode::back() const 617 { 618 assert (last_>=0); 619 return nodes_[best()]; 620 } 621 // Get rid of last one 622 void 623 DBVectorNode::pop_back() 624 { 625 // Temporary until more sophisticated 626 //assert (last_==chosen_); 627 if (nodes_[chosen_].descendants_==2) 628 sizeDeferred_; 629 int previous = nodes_[chosen_].previous_; 630 int next = nodes_[chosen_].next_; 631 nodes_[chosen_].gutsOfDestructor(); 632 if (previous>=0) { 633 nodes_[previous].next_=next; 634 } else { 635 first_ = next; 636 } 637 if (next>=0) { 638 nodes_[next].previous_ = previous; 639 } else { 640 last_ = previous; 641 } 642 nodes_[chosen_].previous_=1; 643 if (firstSpare_>=0) { 644 nodes_[chosen_].next_ = firstSpare_; 645 } else { 646 nodes_[chosen_].next_ = 1; 647 } 648 firstSpare_ = chosen_; 649 chosen_ = 1; 650 assert (size_>0); 651 size_; 652 } 653 #endif 654 // Invoke solver's builtin enumeration algorithm 655 void 656 branchAndBound(OsiSolverInterface & model) { 657 double time1 = CoinCpuTime(); 658 // solve LP 659 model.initialSolve(); 660 int funnyBranching=FUNNY_BRANCHING; 661 662 if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { 663 // Continuous is feasible  find integers 664 int numberIntegers=0; 665 int numberColumns = model.getNumCols(); 666 int iColumn; 667 int i; 668 for (iColumn=0;iColumn<numberColumns;iColumn++) { 669 if( model.isInteger(iColumn)) 670 numberIntegers++; 671 } 672 if (!numberIntegers) { 673 std::cout<<"No integer variables" 674 <<std::endl; 675 return; 676 } 677 int * which = new int[numberIntegers]; // which variables are integer 678 // original bounds 679 int * originalLower = new int[numberIntegers]; 680 int * originalUpper = new int[numberIntegers]; 681 int * relaxedLower = new int[numberIntegers]; 682 int * relaxedUpper = new int[numberIntegers]; 683 { 684 const double * lower = model.getColLower(); 685 const double * upper = model.getColUpper(); 686 numberIntegers=0; 687 for (iColumn=0;iColumn<numberColumns;iColumn++) { 688 if( model.isInteger(iColumn)) { 689 originalLower[numberIntegers]=(int) lower[iColumn]; 690 originalUpper[numberIntegers]=(int) upper[iColumn]; 691 which[numberIntegers++]=iColumn; 692 } 693 } 694 } 695 double direction = model.getObjSense(); 696 // empty tree 697 DBVectorNode branchingTree; 698 699 // Add continuous to it; 700 DBNodeSimple rootNode(model,numberIntegers,which,model.getWarmStart()); 701 // something extra may have been fixed by strong branching 702 // if so go round again 703 while (rootNode.variable_==numberIntegers) { 704 model.resolve(); 705 rootNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart()); 706 } 707 if (rootNode.objectiveValue_<1.0e100) { 708 // push on stack 709 branchingTree.push_back(rootNode); 710 } 711 712 // For printing totals 713 int numberIterations=0; 714 int numberNodes =0; 715 int nRedundantUp=0; 716 int nRedundantDown=0; 717 int nRedundantUp2=0; 718 int nRedundantDown2=0; 719 DBNodeSimple bestNode; 720 ////// Start main while of branch and bound 721 // while until nothing on stack 722 while (branchingTree.size()) { 723 // last node 724 DBNodeSimple node = branchingTree.back(); 725 int kNode = branchingTree.chosen_; 726 branchingTree.pop_back(); 727 assert (node.descendants_<2); 728 numberNodes++; 729 if (node.variable_>=0) { 730 // branch  do bounds 731 for (i=0;i<numberIntegers;i++) { 732 iColumn=which[i]; 733 model.setColBounds( iColumn,node.lower_[i],node.upper_[i]); 734 } 735 // move basis 736 model.setWarmStart(node.basis_); 737 // do branching variable 738 node.incrementDescendants(); 739 if (node.way_<0) { 740 model.setColUpper(which[node.variable_],floor(node.value_)); 741 // now push back node if more to come 742 if (node.way_==1) { 743 node.way_=+2; // Swap direction 744 branchingTree.push_back(node); 745 } else if (funnyBranching) { 746 // put back on tree anyway 747 branchingTree.push_back(node); 748 } 749 } else { 750 model.setColLower(which[node.variable_],ceil(node.value_)); 751 // now push back node if more to come 752 if (node.way_==1) { 753 node.way_=2; // Swap direction 754 branchingTree.push_back(node); 755 } else if (funnyBranching) { 756 // put back on tree anyway 757 branchingTree.push_back(node); 758 } 759 } 760 761 // solve 762 model.resolve(); 763 CoinWarmStart * ws = model.getWarmStart(); 764 const CoinWarmStartBasis* wsb = 765 dynamic_cast<const CoinWarmStartBasis*>(ws); 766 assert (wsb!=NULL); // make sure not volume 767 numberIterations += model.getIterationCount(); 768 // fix on reduced costs 769 int nFixed0=0,nFixed1=0; 770 double cutoff; 771 model.getDblParam(OsiDualObjectiveLimit,cutoff); 772 double gap=(cutoffmodel.getObjValue())*direction+1.0e4; 773 // double gap=(cutoffmodelPtr_>objectiveValue())*direction+1.0e4; 774 if (gap<1.0e10&&model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { 775 const double * dj = model.getReducedCost(); 776 const double * lower = model.getColLower(); 777 const double * upper = model.getColUpper(); 778 for (i=0;i<numberIntegers;i++) { 779 iColumn=which[i]; 780 if (upper[iColumn]>lower[iColumn]) { 781 double djValue = dj[iColumn]*direction; 782 if (wsb>getStructStatus(iColumn)==CoinWarmStartBasis::atLowerBound&& 783 djValue>gap) { 784 nFixed0++; 785 model.setColUpper(iColumn,lower[iColumn]); 786 } else if (wsb>getStructStatus(iColumn)==CoinWarmStartBasis::atUpperBound&& 787 djValue>gap) { 788 nFixed1++; 789 model.setColLower(iColumn,upper[iColumn]); 790 } 791 } 792 } 793 //if (nFixed0+nFixed1) 794 //printf("%d fixed to lower, %d fixed to upper\n",nFixed0,nFixed1); 795 } 796 if (!model.isIterationLimitReached()) { 797 if (model.isProvenOptimal()&&!model.isDualObjectiveLimitReached()) { 798 #if FUNNY_BRANCHING 799 // See if branched variable off bounds 800 const double * dj = model.getReducedCost(); 801 const double * lower = model.getColLower(); 802 const double * upper = model.getColUpper(); 803 const double * solution = model.getColSolution(); 804 // Better to use "natural" value  need flag to say fixed 805 for (i=0;i<numberIntegers;i++) { 806 iColumn=which[i]; 807 relaxedLower[i]=originalLower[i]; 808 relaxedUpper[i]=originalUpper[i]; 809 double djValue = dj[iColumn]*direction; 810 if (djValue>1.0e6) { 811 // wants to go down 812 if (lower[iColumn]>originalLower[i]) { 813 // Lower bound active 814 relaxedLower[i]=(int) lower[iColumn]; 815 } 816 if (upper[iColumn]<originalUpper[i]) { 817 // Upper bound NOT active 818 } 819 } else if (djValue<1.0e6) { 820 // wants to go up 821 if (lower[iColumn]>originalLower[i]) { 822 // Lower bound NOT active 823 } 824 if (upper[iColumn]<originalUpper[i]) { 825 // Upper bound active 826 relaxedUpper[i]=(int) upper[iColumn]; 827 } 828 } 829 } 830 // See if can do anything 831 { 832 /* 833 If kNode is on second branch then 834 a) If other feasible could free up as well 835 b) If other infeasible could do something clever. 836 For now  we have to give up 837 */ 838 int jNode=branchingTree.nodes_[kNode].parent_; 839 bool canDelete = (branchingTree.nodes_[kNode].descendants_<2); 840 while (jNode>=0) { 841 DBNodeSimple & node = branchingTree.nodes_[jNode]; 842 int next = node.parent_; 843 if (node.descendants_<2) { 844 int variable = node.variable_; 845 iColumn=which[variable]; 846 double value = node.value_; 847 double djValue = dj[iColumn]*direction; 848 assert (node.way_==2node.way_==2); 849 // we don't know which branch it was  look at current bounds 850 if (upper[iColumn]<value&&node.lower_[variable]<upper[iColumn]) { 851 // must have been down branch 852 if (djValue>1.0e3solution[iColumn]<upper[iColumn]1.0e5) { 853 if (canDelete) { 854 nRedundantDown++; 855 #if 1 856 printf("%d redundant branch down with value %g current upper %g solution %g dj %g\n", 857 variable,node.value_,upper[iColumn],solution[iColumn],djValue); 858 #endif 859 node.descendants_=2; // ignore 860 branchingTree.sizeDeferred_++; 861 int newUpper = originalUpper[variable]; 862 if (next>=0) { 863 DBNodeSimple & node2 = branchingTree.nodes_[next]; 864 newUpper = node2.upper_[variable]; 865 } 866 if (branchingTree.nodes_[jNode].parent_!=next) 867 assert (newUpper>upper[iColumn]); 868 model.setColUpper(iColumn,newUpper); 869 int kNode2=next; 870 int jNode2=branchingTree.nodes_[kNode].parent_; 871 assert (newUpper>branchingTree.nodes_[kNode].upper_[variable]); 872 branchingTree.nodes_[kNode].upper_[variable]= newUpper; 873 while (jNode2!=kNode2) { 874 DBNodeSimple & node2 = branchingTree.nodes_[jNode2]; 875 int next = node2.parent_; 876 if (next!=kNode2) 877 assert (newUpper>node2.upper_[variable]); 878 node2.upper_[variable]= newUpper; 879 jNode2=next; 880 } 881 } else { 882 // can't delete but can add other way to jNode 883 nRedundantDown2++; 884 DBNodeSimple & node2 = branchingTree.nodes_[kNode]; 885 assert (node2.way_==2node2.way_==2); 886 double value2 = node2.value_; 887 int variable2 = node2.variable_; 888 int iColumn2 = which[variable2]; 889 if (variable != variable2) { 890 if (node2.way_==2&&upper[iColumn2]<value2) { 891 // must have been down branch which was done  carry over 892 int newUpper = (int) floor(value2); 893 assert (newUpper<node.upper_[variable2]); 894 node.upper_[variable2]=newUpper; 895 } else if (node2.way_==2&&lower[iColumn2]>value2) { 896 // must have been up branch which was done  carry over 897 int newLower = (int) ceil(value2); 898 assert (newLower>node.lower_[variable2]); 899 node.lower_[variable2]=newLower; 900 } 901 if (node.lower_[variable2]>node.upper_[variable2]) { 902 // infeasible 903 node.descendants_=2; // ignore 904 branchingTree.sizeDeferred_++; 905 } 906 } 907 } 908 break; 909 } 910 // we don't know which branch it was  look at current bounds 911 } else if (lower[iColumn]>value&&node.upper_[variable]>lower[iColumn]) { 912 // must have been up branch 913 if (djValue<1.0e3solution[iColumn]>lower[iColumn]+1.0e5) { 914 if (canDelete) { 915 nRedundantUp++; 916 #if 1 917 printf("%d redundant branch up with value %g current lower %g solution %g dj %g\n", 918 variable,node.value_,lower[iColumn],solution[iColumn],djValue); 919 #endif 920 node.descendants_=2; // ignore 921 branchingTree.sizeDeferred_++; 922 int newLower = originalLower[variable]; 923 if (next>=0) { 924 DBNodeSimple & node2 = branchingTree.nodes_[next]; 925 newLower = node2.lower_[variable]; 926 } 927 if (branchingTree.nodes_[jNode].parent_!=next) 928 assert (newLower<lower[iColumn]); 929 model.setColLower(iColumn,newLower); 930 int kNode2=next; 931 int jNode2=branchingTree.nodes_[kNode].parent_; 932 assert (newLower<branchingTree.nodes_[kNode].lower_[variable]); 933 branchingTree.nodes_[kNode].lower_[variable]= newLower; 934 while (jNode2!=kNode2) { 935 DBNodeSimple & node2 = branchingTree.nodes_[jNode2]; 936 int next = node2.parent_; 937 if (next!=kNode2) 938 assert (newLower<node2.lower_[variable]); 939 node2.lower_[variable]=newLower; 940 jNode2=next; 941 } 942 } else { 943 // can't delete but can add other way to jNode 944 nRedundantUp2++; 945 DBNodeSimple & node2 = branchingTree.nodes_[kNode]; 946 assert (node2.way_==2node2.way_==2); 947 double value2 = node2.value_; 948 int variable2 = node2.variable_; 949 int iColumn2 = which[variable2]; 950 if (variable != variable2) { 951 if (node2.way_==2&&upper[iColumn2]<value2) { 952 // must have been down branch which was done  carry over 953 int newUpper = (int) floor(value2); 954 assert (newUpper<node.upper_[variable2]); 955 node.upper_[variable2]=newUpper; 956 } else if (node2.way_==2&&lower[iColumn2]>value2) { 957 // must have been up branch which was done  carry over 958 int newLower = (int) ceil(value2); 959 assert (newLower>node.lower_[variable2]); 960 node.lower_[variable2]=newLower; 961 } 962 if (node.lower_[variable2]>node.upper_[variable2]) { 963 // infeasible 964 node.descendants_=2; // ignore 965 branchingTree.sizeDeferred_++; 966 } 967 } 968 } 969 break; 970 } 971 } 972 } else { 973 break; 974 } 975 jNode=next; 976 } 977 } 978 // solve 979 //resolve(); 980 //assert(!getIterationCount()); 981 if ((numberNodes%1000)==0) 982 printf("%d nodes, redundant down %d (%d) up %d (%d) tree size %d\n", 983 numberNodes,nRedundantDown,nRedundantDown2,nRedundantUp,nRedundantUp2,branchingTree.size()); 984 #else 985 if ((numberNodes%1000)==0) 986 printf("%d nodes, tree size %d\n", 987 numberNodes,branchingTree.size()); 988 #endif 989 if (CoinCpuTime()time1>3600.0) { 990 printf("stopping after 3600 seconds\n"); 991 exit(77); 992 } 993 DBNodeSimple newNode(model,numberIntegers,which,ws); 994 // something extra may have been fixed by strong branching 995 // if so go round again 996 while (newNode.variable_==numberIntegers) { 997 model.resolve(); 998 newNode = DBNodeSimple(model,numberIntegers,which,model.getWarmStart()); 999 } 1000 if (newNode.objectiveValue_<1.0e100) { 1001 if (newNode.variable_>=0) 1002 assert (fabs(newNode.value_floor(newNode.value_+0.5))>1.0e6); 1003 newNode.parent_ = kNode; 1004 // push on stack 1005 branchingTree.push_back(newNode); 1006 #if 0 1007 } else { 1008 // integer solution  save 1009 bestNode = node; 1010 // set cutoff (hard coded tolerance) 1011 model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_1.0e5)*direction); 1012 std::cout<<"Integer solution of " 1013 <<bestNode.objectiveValue_ 1014 <<" found after "<<numberIterations 1015 <<" iterations and "<<numberNodes<<" nodes" 1016 <<std::endl; 1017 } 1018 #endif 1019 } 1020 } 1021 } else { 1022 // maximum iterations  exit 1023 std::cout<<"Exiting on maximum iterations" 1024 <<std::endl; 1025 break; 1026 } 1027 } else { 1028 // integer solution  save 1029 bestNode = node; 1030 // set cutoff (hard coded tolerance) 1031 model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_1.0e5)*direction); 1032 std::cout<<"Integer solution of " 1033 <<bestNode.objectiveValue_ 1034 <<" found after "<<numberIterations 1035 <<" iterations and "<<numberNodes<<" nodes" 1036 <<std::endl; 1037 } 1038 } 1039 ////// End main while of branch and bound 1040 std::cout<<"Search took " 1041 <<numberIterations 1042 <<" iterations and "<<numberNodes<<" nodes" 1043 <<std::endl; 1044 if (bestNode.numberIntegers_) { 1045 // we have a solution restore 1046 // do bounds 1047 for (i=0;i<numberIntegers;i++) { 1048 iColumn=which[i]; 1049 model.setColBounds( iColumn,bestNode.lower_[i],bestNode.upper_[i]); 1050 } 1051 // move basis 1052 model.setWarmStart(bestNode.basis_); 1053 // set cutoff so will be good (hard coded tolerance) 1054 model.setDblParam(OsiDualObjectiveLimit,(bestNode.objectiveValue_+1.0e5)*direction); 1055 model.resolve(); 1056 } else { 1057 // modelPtr_>setProblemStatus(1); 1058 } 1059 delete [] which; 1060 delete [] originalLower; 1061 delete [] originalUpper; 1062 delete [] relaxedLower; 1063 delete [] relaxedUpper; 1064 } else { 1065 std::cout<<"The LP relaxation is infeasible" 1066 <<std::endl; 1067 // modelPtr_>setProblemStatus(1); 1068 //throw CoinError("The LP relaxation is infeasible or too expensive", 1069 //"branchAndBound", "OsiClpSolverInterface"); 1070 } 1071 } 1072 1073 1074 1 1075 int main(int argc, char* argv[]) 2 1076 {
Note: See TracChangeset
for help on using the changeset viewer.