- Timestamp:
- Jul 15, 2007 5:28:33 PM (14 years ago)
- Location:
- trunk/Cbc/src
- Files:
-
- 4 added
- 28 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Cbc/src/CbcBranchBase.cpp
r640 r687 179 179 { 180 180 return NULL; 181 } 182 /* Pass in information on branch just done and create CbcObjectUpdateData instance. 183 If object does not need data then backward pointer will be NULL. 184 Assumes can get information from solver */ 185 CbcObjectUpdateData 186 CbcObject::createUpdateInformation(const OsiSolverInterface * solver, 187 const CbcNode * node, 188 const CbcBranchingObject * branchingObject) 189 { 190 return CbcObjectUpdateData(); 181 191 } 182 192 … … 318 328 return *this; 319 329 } 320 321 330 // Default constructor 331 CbcObjectUpdateData::CbcObjectUpdateData() 332 : object_(NULL), 333 way_(0), 334 objectNumber_(-1), 335 change_(0.0), 336 status_(0), 337 intDecrease_(0), 338 branchingValue_(0.0) 339 { 340 } 341 342 // Useful constructor 343 CbcObjectUpdateData::CbcObjectUpdateData (CbcObject * object, 344 int way, 345 double change, 346 int status, 347 int intDecrease, 348 double branchingValue) 349 : object_(object), 350 way_(way), 351 objectNumber_(-1), 352 change_(change), 353 status_(status), 354 intDecrease_(intDecrease), 355 branchingValue_(branchingValue) 356 { 357 } 358 359 // Destructor 360 CbcObjectUpdateData::~CbcObjectUpdateData () 361 { 362 } 363 364 // Copy constructor 365 CbcObjectUpdateData::CbcObjectUpdateData ( const CbcObjectUpdateData & rhs) 366 : object_(rhs.object_), 367 way_(rhs.way_), 368 objectNumber_(rhs.objectNumber_), 369 change_(rhs.change_), 370 status_(rhs.status_), 371 intDecrease_(rhs.intDecrease_), 372 branchingValue_(rhs.branchingValue_) 373 { 374 } 375 376 // Assignment operator 377 CbcObjectUpdateData & 378 CbcObjectUpdateData::operator=( const CbcObjectUpdateData& rhs) 379 { 380 if (this!=&rhs) { 381 object_ = rhs.object_; 382 way_ = rhs.way_; 383 objectNumber_ = rhs.objectNumber_; 384 change_ = rhs.change_; 385 status_ = rhs.status_; 386 intDecrease_ = rhs.intDecrease_; 387 branchingValue_ = rhs.branchingValue_; 388 } 389 return *this; 390 } 391 392 -
trunk/Cbc/src/CbcBranchBase.hpp
r640 r687 15 15 class CbcBranchingObject; 16 16 class OsiChooseVariable; 17 class CbcObjectUpdateData; 17 18 18 19 //############################################################################# … … 190 191 double tolerance) const; 191 192 193 /** Pass in information on branch just done and create CbcObjectUpdateData instance. 194 If object does not need data then backward pointer will be NULL. 195 Assumes can get information from solver */ 196 virtual CbcObjectUpdateData createUpdateInformation(const OsiSolverInterface * solver, 197 const CbcNode * node, 198 const CbcBranchingObject * branchingObject); 199 200 /// Update object by CbcObjectUpdateData 201 virtual void updateInformation(const CbcObjectUpdateData & data) {}; 202 192 203 /// Identifier (normally column number in matrix) 193 204 inline int id() const … … 323 334 {way_=way;}; 324 335 336 /// update model 337 inline void setModel(CbcModel * model) 338 { model_ = model;}; 325 339 /// Return model 326 340 inline CbcModel * model() const … … 489 503 protected: 490 504 }; 505 /* This stores data so an object can be updated 506 */ 507 class CbcObjectUpdateData { 508 509 public: 510 511 /// Default Constructor 512 CbcObjectUpdateData (); 513 514 /// Useful constructor 515 CbcObjectUpdateData (CbcObject * object, 516 int way, 517 double change, 518 int status, 519 int intDecrease_, 520 double branchingValue); 521 522 /// Copy constructor 523 CbcObjectUpdateData ( const CbcObjectUpdateData &); 524 525 /// Assignment operator 526 CbcObjectUpdateData & operator=( const CbcObjectUpdateData& rhs); 527 528 /// Destructor 529 virtual ~CbcObjectUpdateData (); 530 531 532 public: 533 /// data 534 535 /// Object 536 CbcObject * object_; 537 /// Branch as defined by instance of CbcObject 538 int way_; 539 /// Object number 540 int objectNumber_; 541 /// Change in objective 542 double change_; 543 /// Status 0 Optimal, 1 infeasible, 2 unknown 544 int status_; 545 /// Decrease in number unsatisfied 546 int intDecrease_; 547 /// Branching value 548 double branchingValue_; 549 550 }; 491 551 492 552 #endif -
trunk/Cbc/src/CbcBranchDynamic.cpp
r640 r687 267 267 { 268 268 } 269 // Copy some information i.e. just variable stuff 270 void 271 CbcSimpleIntegerDynamicPseudoCost::copySome(CbcSimpleIntegerDynamicPseudoCost * otherObject) 272 { 273 downDynamicPseudoCost_=otherObject->downDynamicPseudoCost_; 274 upDynamicPseudoCost_=otherObject->upDynamicPseudoCost_; 275 sumDownCost_ = otherObject->sumDownCost_; 276 sumUpCost_ = otherObject->sumUpCost_; 277 sumDownChange_ = otherObject->sumDownChange_; 278 sumUpChange_ = otherObject->sumUpChange_; 279 sumDownCostSquared_ = otherObject->sumDownCostSquared_; 280 sumUpCostSquared_ = otherObject->sumUpCostSquared_; 281 sumDownDecrease_ = otherObject->sumDownDecrease_; 282 sumUpDecrease_ = otherObject->sumUpDecrease_; 283 lastDownCost_ = otherObject->lastDownCost_; 284 lastUpCost_ = otherObject->lastUpCost_; 285 lastDownDecrease_ = otherObject->lastDownDecrease_; 286 lastUpDecrease_ = otherObject->lastUpDecrease_; 287 numberTimesDown_ = otherObject->numberTimesDown_; 288 numberTimesUp_ = otherObject->numberTimesUp_; 289 numberTimesDownInfeasible_ = otherObject->numberTimesDownInfeasible_; 290 numberTimesUpInfeasible_ = otherObject->numberTimesUpInfeasible_; 291 numberTimesDownLocalFixed_ = otherObject->numberTimesDownLocalFixed_; 292 numberTimesUpLocalFixed_ = otherObject->numberTimesUpLocalFixed_; 293 numberTimesDownTotalFixed_ = otherObject->numberTimesDownTotalFixed_; 294 numberTimesUpTotalFixed_ = otherObject->numberTimesUpTotalFixed_; 295 numberTimesProbingTotal_ = otherObject->numberTimesProbingTotal_; 296 } 269 297 // Creates a branching object 270 298 CbcBranchingObject * … … 646 674 double downCost = CoinMax((value-below)*downDynamicPseudoCost_,0.0); 647 675 return downCost; 676 } 677 /* Pass in information on branch just done and create CbcObjectUpdateData instance. 678 If object does not need data then backward pointer will be NULL. 679 Assumes can get information from solver */ 680 CbcObjectUpdateData 681 CbcSimpleIntegerDynamicPseudoCost::createUpdateInformation(const OsiSolverInterface * solver, 682 const CbcNode * node, 683 const CbcBranchingObject * branchingObject) 684 { 685 double originalValue=node->objectiveValue(); 686 int originalUnsatisfied = node->numberUnsatisfied(); 687 double objectiveValue = solver->getObjValue()*solver->getObjSense(); 688 int unsatisfied=0; 689 int i; 690 //might be base model - doesn't matter 691 int numberIntegers = model_->numberIntegers();; 692 const double * solution = solver->getColSolution(); 693 //const double * lower = solver->getColLower(); 694 //const double * upper = solver->getColUpper(); 695 double change = CoinMax(0.0,objectiveValue-originalValue); 696 int iStatus; 697 if (solver->isProvenOptimal()) 698 iStatus=0; // optimal 699 else if (solver->isIterationLimitReached() 700 &&!solver->isDualObjectiveLimitReached()) 701 iStatus=2; // unknown 702 else 703 iStatus=1; // infeasible 704 705 bool feasible = iStatus!=1; 706 if (feasible) { 707 double integerTolerance = 708 model_->getDblParam(CbcModel::CbcIntegerTolerance); 709 const int * integerVariable = model_->integerVariable(); 710 for (i=0;i<numberIntegers;i++) { 711 int j=integerVariable[i]; 712 double value = solution[j]; 713 double nearest = floor(value+0.5); 714 if (fabs(value-nearest)>integerTolerance) 715 unsatisfied++; 716 } 717 } 718 int way = branchingObject->way(); 719 way = - way; // because after branch so moved on 720 double value = branchingObject->value(); 721 return CbcObjectUpdateData (this, way, 722 change, iStatus, 723 originalUnsatisfied-unsatisfied,value); 724 } 725 // Update object by CbcObjectUpdateData 726 void 727 CbcSimpleIntegerDynamicPseudoCost::updateInformation(const CbcObjectUpdateData & data) 728 { 729 bool feasible = data.status_!=1; 730 int way = data.way_; 731 double value = data.branchingValue_; 732 double change = data.change_; 733 #define TYPE2 1 734 #define TYPERATIO 0.9 735 if (way<0) { 736 // down 737 if (feasible) { 738 //printf("(down change %g value down %g ",change,value-floor(value)); 739 incrementNumberTimesDown(); 740 addToSumDownChange(1.0e-30+value-floor(value)); 741 addToSumDownDecrease(data.intDecrease_); 742 #if TYPE2==0 743 addToSumDownCost(change/(1.0e-30+(value-floor(value)))); 744 setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown()); 745 #elif TYPE2==1 746 addToSumDownCost(change); 747 setDownDynamicPseudoCost(sumDownCost()/sumDownChange()); 748 #elif TYPE2==2 749 addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value)))); 750 setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown())); 751 #endif 752 } else { 753 //printf("(down infeasible value down %g ",change,value-floor(value)); 754 incrementNumberTimesDownInfeasible(); 755 #if INFEAS==2 756 double distanceToCutoff=0.0; 757 double objectiveValue = model->getCurrentMinimizationObjValue(); 758 distanceToCutoff = model->getCutoff() - originalValue; 759 if (distanceToCutoff<1.0e20) 760 change = distanceToCutoff*2.0; 761 else 762 change = downDynamicPseudoCost()*(value-floor(value))*10.0; 763 change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change); 764 incrementNumberTimesDown(); 765 addToSumDownChange(1.0e-30+value-floor(value)); 766 addToSumDownDecrease(data.intDecrease_); 767 #if TYPE2==0 768 addToSumDownCost(change/(1.0e-30+(value-floor(value)))); 769 setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown()); 770 #elif TYPE2==1 771 addToSumDownCost(change); 772 setDownDynamicPseudoCost(sumDownCost()/sumDownChange()); 773 #elif TYPE2==2 774 addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value)))); 775 setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown())); 776 #endif 777 #endif 778 } 779 } else { 780 // up 781 if (feasible) { 782 //printf("(up change %g value down %g ",change,ceil(value)-value); 783 incrementNumberTimesUp(); 784 addToSumUpChange(1.0e-30+ceil(value)-value); 785 addToSumUpDecrease(data.intDecrease_); 786 #if TYPE2==0 787 addToSumUpCost(change/(1.0e-30+(ceil(value)-value))); 788 setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp()); 789 #elif TYPE2==1 790 addToSumUpCost(change); 791 setUpDynamicPseudoCost(sumUpCost()/sumUpChange()); 792 #elif TYPE2==2 793 addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value))); 794 setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp())); 795 #endif 796 } else { 797 //printf("(up infeasible value down %g ",change,ceil(value)-value); 798 incrementNumberTimesUpInfeasible(); 799 #if INFEAS==2 800 double distanceToCutoff=0.0; 801 double objectiveValue = model->getCurrentMinimizationObjValue(); 802 distanceToCutoff = model->getCutoff() - originalValue; 803 if (distanceToCutoff<1.0e20) 804 change = distanceToCutoff*2.0; 805 else 806 change = upDynamicPseudoCost()*(ceil(value)-value)*10.0; 807 change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change); 808 incrementNumberTimesUp(); 809 addToSumUpChange(1.0e-30+ceil(value)-value); 810 addToSumUpDecrease(data.intDecrease_); 811 #if TYPE2==0 812 addToSumUpCost(change/(1.0e-30+(ceil(value)-value))); 813 setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp()); 814 #elif TYPE2==1 815 addToSumUpCost(change); 816 setUpDynamicPseudoCost(sumUpCost()/sumUpChange()); 817 #elif TYPE2==2 818 addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value))); 819 setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp())); 820 #endif 821 #endif 822 } 823 } 824 //print(1,0.5); 648 825 } 649 826 // Pass in probing information … … 952 1129 int way = object_->way(); 953 1130 double value = object_->value(); 954 #define TYPE2 1955 #define TYPERATIO 0.91131 //#define TYPE2 1 1132 //#define TYPERATIO 0.9 956 1133 if (way<0) { 957 1134 // down … … 1043 1220 } 1044 1221 } 1045 //object->print( );1222 //object->print(1,0.5); 1046 1223 delete object_; 1047 1224 object_=NULL; … … 1061 1238 double changeDown, int numInfDown) 1062 1239 { 1063 int stateOfSearch = thisOne->model()->stateOfSearch(); 1240 CbcModel * model = thisOne->model(); 1241 int stateOfSearch = model->stateOfSearch(); 1064 1242 int betterWay=0; 1065 1243 double value=0.0; 1244 if (!bestObject_) { 1245 bestCriterion_=-1.0; 1246 bestNumberUp_=INT_MAX; 1247 bestNumberDown_=INT_MAX; 1248 } 1066 1249 if (stateOfSearch<=1&&thisOne->model()->currentNode()->depth()>10) { 1067 #if 0 1068 if (!bestObject_) { 1069 bestNumberUp_=INT_MAX; 1070 bestNumberDown_=INT_MAX; 1071 } 1250 #define TRY_STUFF 1 1251 #ifdef TRY_STUFF 1072 1252 // before solution - choose smallest number 1073 1253 // could add in depth as well … … 1108 1288 } 1109 1289 #else 1110 if (!bestObject_) {1111 bestCriterion_=-1.0;1112 }1113 1290 // got a solution 1114 1291 double minValue = CoinMin(changeDown,changeUp); … … 1124 1301 #endif 1125 1302 } else { 1126 if (!bestObject_) { 1127 bestCriterion_=-1.0; 1128 } 1303 #if TRY_STUFF > 1 1304 // Get current number of infeasibilities, cutoff and current objective 1305 CbcNode * node = model->currentNode(); 1306 int numberUnsatisfied = node->numberUnsatisfied(); 1307 double cutoff = model->getCutoff(); 1308 double objectiveValue = node->objectiveValue(); 1309 #endif 1129 1310 // got a solution 1130 1311 double minValue = CoinMin(changeDown,changeUp); 1131 1312 double maxValue = CoinMax(changeDown,changeUp); 1132 1313 // Reduce 1314 #ifdef TRY_STUFF 1315 //maxValue = CoinMin(maxValue,minValue*4.0); 1316 #else 1133 1317 maxValue = CoinMin(maxValue,minValue*2.0); 1318 #endif 1134 1319 value = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue; 1135 if (value>bestCriterion_+1.0e-8) { 1320 double useValue = value; 1321 double useBest = bestCriterion_; 1322 #if TRY_STUFF>1 1323 int thisNumber = CoinMin(numInfUp,numInfDown); 1324 int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_); 1325 double distance = cutoff-objectiveValue; 1326 assert (distance>=0.0); 1327 if (useValue+0.1*distance>useBest&&useValue*1.1>useBest&& 1328 useBest+0.1*distance>useValue&&useBest*1.1>useValue) { 1329 // not much in it - look at unsatisfied 1330 if (thisNumber<numberUnsatisfied||bestNumber<numberUnsatisfied) { 1331 double perInteger = distance/ ((double) numberUnsatisfied); 1332 useValue += thisNumber*perInteger; 1333 useBest += bestNumber*perInteger; 1334 } 1335 } 1336 #endif 1337 if (useValue>useBest+1.0e-8) { 1136 1338 if (changeUp<=changeDown) { 1137 1339 betterWay=1; -
trunk/Cbc/src/CbcBranchDynamic.hpp
r640 r687 64 64 virtual CbcBranchingObject * createBranch(OsiSolverInterface * solver, 65 65 const OsiBranchingInformation * info, int way) ; 66 67 /** Pass in information on branch just done and create CbcObjectUpdateData instance. 68 If object does not need data then backward pointer will be NULL. 69 Assumes can get information from solver */ 70 virtual CbcObjectUpdateData createUpdateInformation(const OsiSolverInterface * solver, 71 const CbcNode * node, 72 const CbcBranchingObject * branchingObject); 73 /// Update object by CbcObjectUpdateData 74 virtual void updateInformation(const CbcObjectUpdateData & data) ; 75 /// Copy some information i.e. just variable stuff 76 void copySome(CbcSimpleIntegerDynamicPseudoCost * otherObject); 66 77 67 78 /** Create an OsiSolverBranch object -
trunk/Cbc/src/CbcCompareBase.hpp
r485 r687 18 18 */ 19 19 #include "CbcNode.hpp" 20 #include "CbcConfig.h" 20 21 21 22 class CbcModel; … … 76 77 { 77 78 assert (x); 79 assert (y); 80 #ifndef CBC_THREAD 78 81 CbcNodeInfo * infoX = x->nodeInfo(); 79 82 assert (infoX); 80 83 int nodeNumberX = infoX->nodeNumber(); 81 assert (y);82 84 CbcNodeInfo * infoY = y->nodeInfo(); 83 85 assert (infoY); … … 85 87 assert (nodeNumberX!=nodeNumberY); 86 88 return (nodeNumberX>nodeNumberY); 89 #else 90 // doesn't work if threaded 91 assert (x!=y); 92 return (x>y); 93 #endif 87 94 }; 88 95 protected: -
trunk/Cbc/src/CbcCountRowCut.cpp
r310 r687 41 41 int whichGenerator) 42 42 : OsiRowCut(rhs), 43 owner_(info),44 ownerCut_(whichOne),43 owner_(info), 44 ownerCut_(whichOne), 45 45 numberPointingToThis_(0), 46 whichCutGenerator_(whichGenerator)46 whichCutGenerator_(whichGenerator) 47 47 { 48 48 #ifdef CHECK_CUT_COUNTS … … 58 58 // Look at owner and delete 59 59 owner_->deleteCut(ownerCut_); 60 ownerCut_=-1234567; 60 61 } 61 62 // Increment number of references … … 63 64 CbcCountRowCut::increment(int change) 64 65 { 66 assert(ownerCut_!=-1234567); 65 67 numberPointingToThis_+=change; 66 68 } … … 70 72 CbcCountRowCut::decrement(int change) 71 73 { 72 assert(numberPointingToThis_>=change); 74 assert(ownerCut_!=-1234567); 75 //assert(numberPointingToThis_>=change); 76 assert(numberPointingToThis_>=0); 77 if(numberPointingToThis_<change) { 78 assert(numberPointingToThis_>0); 79 printf("negative cut count %d - %d\n",numberPointingToThis_, change); 80 change = numberPointingToThis_; 81 } 73 82 numberPointingToThis_-=change; 74 83 return numberPointingToThis_; 75 84 } 85 76 86 // Set information 77 87 void -
trunk/Cbc/src/CbcCutGenerator.cpp
r640 r687 83 83 model_ = rhs.model_; 84 84 generator_=rhs.generator_->clone(); 85 generator_->refreshSolver(model_->solver());85 //generator_->refreshSolver(model_->solver()); 86 86 whenCutGenerator_=rhs.whenCutGenerator_; 87 87 whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_; … … 160 160 */ 161 161 bool 162 CbcCutGenerator::generateCuts( OsiCuts & cs , bool fullScan, CbcNode * node)162 CbcCutGenerator::generateCuts( OsiCuts & cs , bool fullScan, OsiSolverInterface * solver, CbcNode * node) 163 163 { 164 164 int howOften = whenCutGenerator_; … … 172 172 howOften=1; 173 173 bool returnCode=false; 174 OsiSolverInterface * solver = model_->solver();174 //OsiSolverInterface * solver = model_->solver(); 175 175 int depth; 176 176 if (node) … … 210 210 if (!generator) { 211 211 // Pass across model information in case it could be useful 212 //OsiSolverInterface * solver = model_->solver();213 212 //void * saveData = solver->getApplicationData(); 214 213 //solver->setApplicationData(model_); … … 235 234 int numberColumns = solver->getNumCols(); 236 235 double primalTolerance = 1.0e-8; 237 #if 0238 int numberChanged=0,ifCut=0;239 CoinPackedVector lbs;240 CoinPackedVector ubs;241 for (j=0;j<numberColumns;j++) {242 if (solver->isInteger(j)) {243 if (tightUpper[j]<upper[j]) {244 numberChanged++;245 assert (tightUpper[j]==floor(tightUpper[j]+0.5));246 ubs.insert(j,tightUpper[j]);247 if (tightUpper[j]<solution[j]-primalTolerance)248 ifCut=1;249 }250 if (tightLower[j]>lower[j]) {251 numberChanged++;252 assert (tightLower[j]==floor(tightLower[j]+0.5));253 lbs.insert(j,tightLower[j]);254 if (tightLower[j]>solution[j]+primalTolerance)255 ifCut=1;256 }257 } else {258 if (tightUpper[j]==tightLower[j]&&259 upper[j]>lower[j]) {260 // fix261 //solver->setColLower(j,tightLower[j]);262 //solver->setColUpper(j,tightUpper[j]);263 double value = tightUpper[j];264 numberChanged++;265 if (value<upper[j])266 ubs.insert(j,value);267 if (value>lower[j])268 lbs.insert(j,value);269 }270 }271 }272 if (numberChanged) {273 OsiColCut cc;274 cc.setUbs(ubs);275 cc.setLbs(lbs);276 if (ifCut) {277 cc.setEffectiveness(100.0);278 } else {279 cc.setEffectiveness(1.0e-5);280 }281 cs.insert(cc);282 }283 // need to resolve if some bounds changed284 returnCode = !solver->basisIsAvailable();285 assert (!returnCode);286 #else287 236 const char * tightenBounds = generator->tightenBounds(); 288 for (j=0;j<numberColumns;j++) { 289 if (solver->isInteger(j)) { 290 if (tightUpper[j]<upper[j]) { 291 double nearest = floor(tightUpper[j]+0.5); 292 //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible 293 solver->setColUpper(j,nearest); 294 if (nearest<solution[j]-primalTolerance) 295 returnCode=true; 296 } 297 if (tightLower[j]>lower[j]) { 298 double nearest = floor(tightLower[j]+0.5); 299 //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible 300 solver->setColLower(j,nearest); 301 if (nearest>solution[j]+primalTolerance) 302 returnCode=true; 303 } 304 } else { 305 if (upper[j]>lower[j]) { 306 if (tightUpper[j]==tightLower[j]) { 307 // fix 308 solver->setColLower(j,tightLower[j]); 309 solver->setColUpper(j,tightUpper[j]); 310 if (tightLower[j]>solution[j]+primalTolerance|| 311 tightUpper[j]<solution[j]-primalTolerance) 312 returnCode=true; 313 } else if (tightenBounds&&tightenBounds[j]) { 314 solver->setColLower(j,CoinMax(tightLower[j],lower[j])); 315 solver->setColUpper(j,CoinMin(tightUpper[j],upper[j])); 316 if (tightLower[j]>solution[j]+primalTolerance|| 317 tightUpper[j]<solution[j]-primalTolerance) 237 if ((model_->getThreadMode()&2)==0) { 238 for (j=0;j<numberColumns;j++) { 239 if (solver->isInteger(j)) { 240 if (tightUpper[j]<upper[j]) { 241 double nearest = floor(tightUpper[j]+0.5); 242 //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible 243 solver->setColUpper(j,nearest); 244 if (nearest<solution[j]-primalTolerance) 318 245 returnCode=true; 319 246 } 247 if (tightLower[j]>lower[j]) { 248 double nearest = floor(tightLower[j]+0.5); 249 //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible 250 solver->setColLower(j,nearest); 251 if (nearest>solution[j]+primalTolerance) 252 returnCode=true; 253 } 254 } else { 255 if (upper[j]>lower[j]) { 256 if (tightUpper[j]==tightLower[j]) { 257 // fix 258 solver->setColLower(j,tightLower[j]); 259 solver->setColUpper(j,tightUpper[j]); 260 if (tightLower[j]>solution[j]+primalTolerance|| 261 tightUpper[j]<solution[j]-primalTolerance) 262 returnCode=true; 263 } else if (tightenBounds&&tightenBounds[j]) { 264 solver->setColLower(j,CoinMax(tightLower[j],lower[j])); 265 solver->setColUpper(j,CoinMin(tightUpper[j],upper[j])); 266 if (tightLower[j]>solution[j]+primalTolerance|| 267 tightUpper[j]<solution[j]-primalTolerance) 268 returnCode=true; 269 } 270 } 320 271 } 321 } 272 } 273 } else { 274 CoinPackedVector lbs; 275 CoinPackedVector ubs; 276 int numberChanged=0; 277 bool ifCut=false; 278 for (j=0;j<numberColumns;j++) { 279 if (solver->isInteger(j)) { 280 if (tightUpper[j]<upper[j]) { 281 double nearest = floor(tightUpper[j]+0.5); 282 //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible 283 ubs.insert(j,nearest); 284 numberChanged++; 285 if (nearest<solution[j]-primalTolerance) 286 ifCut=true; 287 } 288 if (tightLower[j]>lower[j]) { 289 double nearest = floor(tightLower[j]+0.5); 290 //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible 291 lbs.insert(j,nearest); 292 numberChanged++; 293 if (nearest>solution[j]+primalTolerance) 294 ifCut=true; 295 } 296 } else { 297 if (upper[j]>lower[j]) { 298 if (tightUpper[j]==tightLower[j]) { 299 // fix 300 lbs.insert(j,tightLower[j]); 301 ubs.insert(j,tightUpper[j]); 302 if (tightLower[j]>solution[j]+primalTolerance|| 303 tightUpper[j]<solution[j]-primalTolerance) 304 ifCut=true; 305 } else if (tightenBounds&&tightenBounds[j]) { 306 lbs.insert(j,CoinMax(tightLower[j],lower[j])); 307 ubs.insert(j,CoinMin(tightUpper[j],upper[j])); 308 if (tightLower[j]>solution[j]+primalTolerance|| 309 tightUpper[j]<solution[j]-primalTolerance) 310 ifCut=true; 311 } 312 } 313 } 314 } 315 if (numberChanged) { 316 OsiColCut cc; 317 cc.setUbs(ubs); 318 cc.setLbs(lbs); 319 if (ifCut) { 320 cc.setEffectiveness(100.0); 321 } else { 322 cc.setEffectiveness(1.0e-5); 323 } 324 cs.insert(cc); 325 } 322 326 } 323 327 //if (!solver->basisIsAvailable()) 324 328 //returnCode=true; 325 #endif326 329 #if 0 327 330 // Pass across info to pseudocosts … … 352 355 int k ; 353 356 int nOdd=0; 354 const OsiSolverInterface * solver = model_->solver();357 //const OsiSolverInterface * solver = model_->solver(); 355 358 for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) { 356 359 OsiRowCut & thisCut = cs.rowCut(k) ; … … 382 385 for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 383 386 OsiRowCut thisCut = cs.rowCut(k) ; 387 if (thisCut.lb()>thisCut.ub()|| 388 thisCut.lb()>1.0e8|| 389 thisCut.ub()<-1.0e8) 390 printf("cut from %s has bounds %g and %g!\n", 391 generatorName_,thisCut.lb(),thisCut.ub()); 384 392 if (thisCut.lb()<=thisCut.ub()) { 385 393 /* check size of elements. -
trunk/Cbc/src/CbcCutGenerator.hpp
r640 r687 64 64 If node then can find out depth 65 65 */ 66 bool generateCuts( OsiCuts &cs, bool fullScan, CbcNode * node);66 bool generateCuts( OsiCuts &cs, bool fullScan,OsiSolverInterface * solver,CbcNode * node); 67 67 //@} 68 68 … … 177 177 inline double timeInCutGenerator() const 178 178 { return timeInCutGenerator_;}; 179 inline void incrementTimeInCutGenerator(double value) 180 { timeInCutGenerator_ += value;}; 179 181 /// Get the \c CglCutGenerator corresponding to this \c CbcCutGenerator. 180 182 inline CglCutGenerator * generator() const … … 237 239 inline void setNumberActiveCutsAtRoot(int value) 238 240 { numberActiveCutsAtRoot_ = value;}; 241 /// Set model 242 inline void setModel(CbcModel * model) 243 { model_ = model;}; 239 244 //@} 240 245 -
trunk/Cbc/src/CbcFathom.cpp
r310 r687 9 9 #include <cfloat> 10 10 11 #include "OsiSolverInterface.hpp" 11 #ifdef COIN_HAS_CLP 12 #include "OsiClpSolverInterface.hpp" 13 #endif 12 14 #include "CbcModel.hpp" 13 15 #include "CbcMessage.hpp" … … 45 47 model_ = model; 46 48 } 49 #ifdef COIN_HAS_CLP 50 51 //############################################################################# 52 // Constructors, destructors clone and assignment 53 //############################################################################# 54 55 //------------------------------------------------------------------- 56 // Default Constructor 57 //------------------------------------------------------------------- 58 CbcOsiSolver::CbcOsiSolver () 59 : OsiClpSolverInterface() 60 { 61 cbcModel_ = NULL; 62 } 63 //------------------------------------------------------------------- 64 // Clone 65 //------------------------------------------------------------------- 66 OsiSolverInterface * 67 CbcOsiSolver::clone(bool copyData) const 68 { 69 assert (copyData); 70 return new CbcOsiSolver(*this); 71 } 72 73 74 //------------------------------------------------------------------- 75 // Copy constructor 76 //------------------------------------------------------------------- 77 CbcOsiSolver::CbcOsiSolver ( 78 const CbcOsiSolver & rhs) 79 : OsiClpSolverInterface(rhs) 80 { 81 cbcModel_ = rhs.cbcModel_; 82 } 83 84 //------------------------------------------------------------------- 85 // Destructor 86 //------------------------------------------------------------------- 87 CbcOsiSolver::~CbcOsiSolver () 88 { 89 } 90 91 //------------------------------------------------------------------- 92 // Assignment operator 93 //------------------------------------------------------------------- 94 CbcOsiSolver & 95 CbcOsiSolver::operator=(const CbcOsiSolver& rhs) 96 { 97 if (this != &rhs) { 98 OsiClpSolverInterface::operator=(rhs); 99 cbcModel_ = rhs.cbcModel_; 100 } 101 return *this; 102 } 103 #endif -
trunk/Cbc/src/CbcFathom.hpp
r36 r687 3 3 #ifndef CbcFathom_H 4 4 #define CbcFathom_H 5 #include "CbcConfig.h" 5 6 6 7 class CbcModel; … … 63 64 64 65 }; 66 #ifdef COIN_HAS_CLP 67 #include "OsiClpSolverInterface.hpp" 65 68 69 //############################################################################# 70 71 /** 72 73 This is for codes where solver needs to know about CbcModel 74 */ 75 76 class CbcOsiSolver : public OsiClpSolverInterface { 77 78 public: 79 80 /**@name Constructors and destructors */ 81 //@{ 82 /// Default Constructor 83 CbcOsiSolver (); 84 85 /// Clone 86 virtual OsiSolverInterface * clone(bool copyData=true) const; 87 88 /// Copy constructor 89 CbcOsiSolver (const CbcOsiSolver &); 90 91 /// Assignment operator 92 CbcOsiSolver & operator=(const CbcOsiSolver& rhs); 93 94 /// Destructor 95 virtual ~CbcOsiSolver (); 96 97 //@} 98 99 100 /**@name Sets and Gets */ 101 //@{ 102 /// Set Cbc Model 103 inline void setCbcModel(CbcModel * model) 104 { cbcModel_=model;}; 105 /// Return Cbc Model 106 inline CbcModel * cbcModel() const 107 { return cbcModel_;}; 108 //@} 109 110 //--------------------------------------------------------------------------- 111 112 protected: 113 114 115 /**@name Private member data */ 116 //@{ 117 /// Pointer back to CbcModel 118 CbcModel * cbcModel_; 119 //@} 120 }; 66 121 #endif 122 #endif -
trunk/Cbc/src/CbcHeuristic.cpp
r640 r687 178 178 // probably faster to use a basis to get integer solutions 179 179 model.setSpecialOptions(2); 180 #ifdef CBC_THREAD 181 if (model_->getNumberThreads()>0&&(model_->getThreadMode()&1)!=0) { 182 // See if at root node 183 bool atRoot = model_->getNodeCount()==0; 184 int passNumber = model_->getCurrentPassNumber(); 185 if (atRoot&&passNumber==1) 186 model.setNumberThreads(model_->getNumberThreads()); 187 } 188 #endif 180 189 model.branchAndBound(); 181 190 if (logLevel>1) … … 302 311 (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3))) 303 312 return 0; // switched off 304 305 313 OsiSolverInterface * solver = model_->solver(); 306 314 const double * lower = solver->getColLower(); … … 315 323 316 324 int numberRows = matrix_.getNumRows(); 317 325 assert (numberRows<=solver->getNumRows()); 318 326 int numberIntegers = model_->numberIntegers(); 319 327 const int * integerVariable = model_->integerVariable(); … … 322 330 double newSolutionValue = direction*solver->getObjValue(); 323 331 int returnCode = 0; 324 325 332 // Column copy 326 333 const double * element = matrix_.getElements(); … … 893 900 } 894 901 } 902 // Just in case of some stupidity 903 double objOffset=0.0; 904 solver->getDblParam(OsiObjOffset,objOffset); 905 newSolutionValue = -objOffset; 906 for ( i=0 ; i<numberColumns ; i++ ) 907 newSolutionValue += objective[i]*newSolution[i]; 908 newSolutionValue *= direction; 909 //printf("new solution value %g %g\n",newSolutionValue,solutionValue); 895 910 if (newSolutionValue<solutionValue) { 896 911 // paranoid check -
trunk/Cbc/src/CbcHeuristic.hpp
r640 r687 79 79 inline int numberNodes() const 80 80 { return numberNodes_;}; 81 /// Just set model - do not do anything else 82 inline void setModelOnly(CbcModel * model) 83 { model_ = model;}; 84 81 85 82 86 /// Sets fraction of new(rows+columns)/old(rows+columns) before doing small branch and bound (default 1.0) -
trunk/Cbc/src/CbcHeuristicFPump.cpp
r640 r687 175 175 double * betterSolution) 176 176 { 177 #define LEN_PRINT 132 178 char pumpPrint[LEN_PRINT]; 179 pumpPrint[0]='\0'; 177 180 if (!when()||(when()==1&&model_->phase()!=1)) 178 181 return 0; // switched off … … 364 367 } 365 368 if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break; 366 numberPasses++;367 369 memcpy(newSolution,solution,numberColumns*sizeof(double)); 368 370 int flip; 369 returnCode = rounds(solver,newSolution,saveObjective,numberIntegers,integerVariable, 370 roundExpensive_,defaultRounding_,&flip); 371 returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable, 372 pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip); 373 numberPasses++; 371 374 if (returnCode) { 372 375 // SOLUTION IS INTEGER … … 382 385 newSolutionValue += saveObjective[i]*newSolution[i]; 383 386 newSolutionValue *= direction; 384 if (model_->logLevel()) 385 printf(" - solution found of %g",newSolutionValue); 387 sprintf(pumpPrint+strlen(pumpPrint)," - solution found of %g",newSolutionValue); 386 388 newLineNeeded=false; 387 389 if (newSolutionValue<solutionValue) { … … 416 418 solutionFound=true; 417 419 if (general&&saveValue!=newSolutionValue) 418 printf(" - cleaned solution of %g\n",solutionValue); 419 else 420 printf("\n"); 420 sprintf(pumpPrint+strlen(pumpPrint)," - cleaned solution of %g",solutionValue); 421 if (pumpPrint[0]!='\0') 422 model_->messageHandler()->message(CBC_FPUMP1,model_->messages()) 423 << pumpPrint 424 <<CoinMessageEol; 425 pumpPrint[0]='\0'; 421 426 } else { 422 if (model_->logLevel()) 423 printf(" - not good enough after small branch and bound\n"); 427 sprintf(pumpPrint+strlen(pumpPrint)," - not good enough after mini branch and bound"); 428 model_->messageHandler()->message(CBC_FPUMP1,model_->messages()) 429 << pumpPrint 430 <<CoinMessageEol; 431 pumpPrint[0]='\0'; 424 432 } 425 433 } else { 426 if (model_->logLevel()) 427 printf(" - not good enough\n"); 434 sprintf(pumpPrint+strlen(pumpPrint)," - not good enough"); 435 model_->messageHandler()->message(CBC_FPUMP1,model_->messages()) 436 << pumpPrint 437 <<CoinMessageEol; 438 pumpPrint[0]='\0'; 428 439 newLineNeeded=false; 429 440 returnCode=0; … … 448 459 if (matched || numberPasses%100 == 0) { 449 460 // perturbation 450 if (model_->logLevel()) 451 printf("Perturbation applied"); 461 sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied"); 452 462 newLineNeeded=true; 453 463 for (i=0;i<numberIntegers;i++) { … … 531 541 memcpy(newSolution,solution,numberColumns*sizeof(double)); 532 542 int flip; 533 returnCode = rounds(solver,newSolution,saveObjective,numberIntegers,integerVariable, 534 roundExpensive_,defaultRounding_,&flip); 543 returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable, 544 pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip); 545 numberPasses++; 535 546 if (returnCode) { 536 547 // solution - but may not be better … … 540 551 newSolutionValue += saveObjective[i]*newSolution[i]; 541 552 newSolutionValue *= direction; 542 if (model_->logLevel()) 543 printf(" - intermediate solution found of %g",newSolutionValue); 553 sprintf(pumpPrint+strlen(pumpPrint)," - intermediate solution found of %g",newSolutionValue); 544 554 if (newSolutionValue<solutionValue) { 545 555 memcpy(betterSolution,newSolution,numberColumns*sizeof(double)); … … 636 646 } 637 647 } 638 if (model_->logLevel()) 639 printf("\npass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue()); 648 if (pumpPrint[0]!='\0') 649 model_->messageHandler()->message(CBC_FPUMP1,model_->messages()) 650 << pumpPrint 651 <<CoinMessageEol; 652 pumpPrint[0]='\0'; 653 sprintf(pumpPrint+strlen(pumpPrint),"Pass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue()); 640 654 if (closestSolution&&solver->getObjValue()<closestObjectiveValue) { 641 655 int i; … … 656 670 scaleFactor *= weightFactor_; 657 671 } // END WHILE 658 if (newLineNeeded&&model_->logLevel()) 659 printf(" - no solution found\n"); 672 if (!solutionFound) { 673 sprintf(pumpPrint+strlen(pumpPrint),"No solution found this major pass"); 674 model_->messageHandler()->message(CBC_FPUMP1,model_->messages()) 675 << pumpPrint 676 <<CoinMessageEol; 677 pumpPrint[0]='\0'; 678 } 660 679 delete solver; 661 680 for ( j=0;j<NUMBER_OLD;j++) … … 724 743 newSolver->initialSolve(); 725 744 if (!newSolver->isProvenOptimal()) { 726 newSolver->writeMps("bad.mps");745 //newSolver->writeMps("bad.mps"); 727 746 assert (newSolver->isProvenOptimal()); 728 747 break; 729 748 } 730 printf("%d integers at bound fixed and %d continuous",749 sprintf(pumpPrint+strlen(pumpPrint),"Before mini branch and bound, %d integers at bound fixed and %d continuous", 731 750 nFix,nFixC); 732 if (nFixC2+nFixI==0) 733 printf("\n"); 734 else 735 printf(" of which %d were internal integer and %d internal continuous\n", 736 nFixI,nFixC2); 751 if (nFixC2+nFixI!=0) 752 sprintf(pumpPrint+strlen(pumpPrint)," of which %d were internal integer and %d internal continuous", 753 nFixI,nFixC2); 754 model_->messageHandler()->message(CBC_FPUMP1,model_->messages()) 755 << pumpPrint 756 <<CoinMessageEol; 757 pumpPrint[0]='\0'; 737 758 double saveValue = newSolutionValue; 738 759 returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue, 739 newSolutionValue,"CbcHeuristicLocalAfterFPump");760 cutoff,"CbcHeuristicLocalAfterFPump"); 740 761 if ((returnCode&2)!=0) { 741 762 // could add cut … … 743 764 } 744 765 if (returnCode) { 745 printf("old sol of %g new of %g\n",saveValue,newSolutionValue); 766 sprintf(pumpPrint+strlen(pumpPrint),"Mini branch and bound improved solution from %g to %g",saveValue,newSolutionValue); 767 model_->messageHandler()->message(CBC_FPUMP1,model_->messages()) 768 << pumpPrint 769 <<CoinMessageEol; 770 pumpPrint[0]='\0'; 746 771 memcpy(betterSolution,newSolution,numberColumns*sizeof(double)); 747 772 if (fixContinuous) { … … 763 788 double value = newSolver->getObjValue()*newSolver->getObjSense(); 764 789 if (value<newSolutionValue) { 765 printf("freeing continuous gives a solution of %g\n", value); 790 sprintf(pumpPrint+strlen(pumpPrint),"Freeing continuous variables gives a solution of %g", value); 791 model_->messageHandler()->message(CBC_FPUMP1,model_->messages()) 792 << pumpPrint 793 <<CoinMessageEol; 794 pumpPrint[0]='\0'; 766 795 newSolutionValue=value; 767 796 memcpy(betterSolution,newSolver->getColSolution(),numberColumns*sizeof(double)); 768 797 } 769 798 } else { 770 newSolver->writeMps("bad3.mps");799 //newSolver->writeMps("bad3.mps"); 771 800 } 772 801 } … … 786 815 double gap = relativeIncrement_*fabs(solutionValue); 787 816 cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement()); 788 printf("round again with cutoff of %g\n",cutoff); 817 sprintf(pumpPrint+strlen(pumpPrint),"Round again with cutoff of %g",cutoff); 818 model_->messageHandler()->message(CBC_FPUMP1,model_->messages()) 819 << pumpPrint 820 <<CoinMessageEol; 821 pumpPrint[0]='\0'; 789 822 if ((accumulate_&3)<2&&usedColumn) 790 823 memset(usedColumn,0,numberColumns); 791 totalNumberPasses += numberPasses ;824 totalNumberPasses += numberPasses-1; 792 825 } else { 793 826 break; … … 817 850 newSolver->addRow(numberIntegersOrig,closestSolution, 818 851 lastSolution,-COIN_DBL_MAX,rhs+10.0); 819 double saveValue = newSolutionValue;852 //double saveValue = newSolutionValue; 820 853 //newSolver->writeMps("sub"); 821 854 int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue, … … 857 890 const double * objective, 858 891 int numberIntegers, const int * integerVariable, 892 char * pumpPrint, int & iter, 859 893 bool roundExpensive, double downValue, int *flip) 860 894 { … … 865 899 int i; 866 900 867 static int iter = 0;868 901 int numberColumns = model_->getNumCols(); 869 902 // tmp contains the current obj coefficients … … 905 938 906 939 if (nnv > nn) nnv = nn; 907 if (iter != 0 &&model_->logLevel())908 printf("up = %5d , down = %5d", flip_up, flip_down); fflush(stdout);940 if (iter != 0) 941 sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down); 909 942 *flip = flip_up + flip_down; 910 943 delete [] tmp; … … 913 946 const double * columnUpper = solver->getColUpper(); 914 947 if (*flip == 0 && iter != 0) { 915 if(model_->logLevel()) 916 printf(" -- rand = %4d (%4d) ", nnv, nn); 948 sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn); 917 949 for (i = 0; i < nnv; i++) { 918 950 // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6 … … 930 962 } 931 963 *flip = nnv; 932 } else if (model_->logLevel()){933 printf(" ");964 } else { 965 //sprintf(pumpPrint+strlen(pumpPrint)," "); 934 966 } 935 967 delete [] list; delete [] val; 936 iter++;968 //iter++; 937 969 938 970 const double * rowLower = solver->getRowLower(); -
trunk/Cbc/src/CbcHeuristicFPump.hpp
r640 r687 165 165 int rounds(OsiSolverInterface * solver,double * solution, const double * objective, 166 166 int numberIntegers, const int * integerVariable, 167 char * pumpPrint,int & passNumber, 167 168 bool roundExpensive=false, 168 169 double downValue=0.5, int *flip=0); -
trunk/Cbc/src/CbcHeuristicGreedy.cpp
r640 r687 764 764 memcpy(betterSolution,newSolution,numberColumns*sizeof(double)); 765 765 solutionValue = newSolutionValue; 766 model_->messageHandler()->message(CBC_HEURISTIC_SOLUTION,model_->messages())767 << newSolutionValue768 << "CbcHeuristicGreedy"<<model_->getCurrentSeconds()769 <<CoinMessageEol;770 766 returnCode=1; 771 767 } -
trunk/Cbc/src/CbcHeuristicLocal.cpp
r640 r687 549 549 returnCode=1; 550 550 solutionValue = newSolutionValue + bestChange; 551 if (bestChange<-1.0e-1)552 model_->messageHandler()->message(CBC_HEURISTIC_SOLUTION,model_->messages())553 << solutionValue554 << "CbcHeuristicLocal"<<model_->getCurrentSeconds()555 <<CoinMessageEol;556 551 } else { 557 552 // bad solution - should not happen so debug if see message -
trunk/Cbc/src/CbcLinked.cpp
r640 r687 3 3 #include "CbcConfig.h" 4 4 5 #ifndef COIN_HAS_LINK6 #ifdef COIN_HAS_ASL7 #define COIN_HAS_LINK8 #endif9 #endif10 5 #include "CoinTime.hpp" 11 6 … … 13 8 #include "CoinModel.hpp" 14 9 #include "ClpSimplex.hpp" 10 #ifdef COIN_HAS_ASL 11 #include "ClpAmplObjective.hpp" 12 #endif 15 13 // returns jColumn (-2 if linear term, -1 if unknown) and coefficient 16 14 static … … 92 90 return jColumn; 93 91 } 94 #i fdef COIN_HAS_LINK92 #include "ClpQuadraticObjective.hpp" 95 93 #include <cassert> 96 94 #if defined(_MSC_VER) … … 105 103 #include "ClpPackedMatrix.hpp" 106 104 #include "CoinTime.hpp" 107 #include "ClpQuadraticObjective.hpp"108 105 #include "CbcModel.hpp" 109 106 #include "CbcCutGenerator.hpp" … … 180 177 //sprintf(temp,"cc%d",iPass); 181 178 //writeMps(temp); 179 //writeMps("tight"); 180 //exit(33); 182 181 //printf("wrote cc%d\n",iPass); 183 182 OsiClpSolverInterface::initialSolve(); … … 185 184 if (modelPtr_->status()==0&&(secondaryStatus==2||secondaryStatus==4)) 186 185 modelPtr_->cleanup(1); 187 if (!isProvenOptimal())188 186 //if (!isProvenOptimal()) 187 //writeMps("yy"); 189 188 if (isProvenOptimal()&&quadraticModel_&&modelPtr_->numberColumns()==quadraticModel_->numberColumns()) { 190 189 // see if qp can get better solution … … 233 232 int numberGenerators = cbcModel_->numberCutGenerators(); 234 233 int iGenerator; 234 cbcModel_->lockThread(); 235 235 for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) { 236 236 CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator); … … 263 263 } 264 264 } 265 cbcModel_->unlockThread(); 265 266 } 266 267 } … … 461 462 bestSolution_ = CoinCopyOfArray(modelPtr_->getColSolution(),modelPtr_->getNumCols()); 462 463 bestObjectiveValue_ = value; 464 if (maxIts<=10000&&cbcModel_) { 465 OsiSolverLink * solver2 = dynamic_cast<OsiSolverLink *> (cbcModel_->solver()); 466 assert (solver2); 467 if (solver2!=this) { 468 // in strong branching - need to store in original solver 469 if (value<solver2->bestObjectiveValue_-1.0e-3) { 470 delete [] solver2->bestSolution_; 471 solver2->bestSolution_ = CoinCopyOfArray(bestSolution_,modelPtr_->getNumCols()); 472 solver2->bestObjectiveValue_ = value; 473 } 474 } 475 } 463 476 // If model has stored then add cut (if convex) 464 477 if (cbcModel_&&(specialOptions2_&4)!=0&&quadraticModel_) { … … 470 483 CglStored * gen2 = dynamic_cast<CglStored *> (gen); 471 484 if (gen2) { 485 cbcModel_->lockThread(); 472 486 // add OA cut 473 double offset ;487 double offset=0.0; 474 488 int numberColumns = quadraticModel_->numberColumns(); 475 489 double * gradient = new double [numberColumns+1]; … … 492 506 int yColumn = obj->yColumn(); 493 507 if (xColumn!=yColumn) { 494 double coefficient = 2.0*obj->coefficient();508 double coefficient = /* 2.0* */obj->coefficient(); 495 509 gradient[xColumn] += coefficient*solution[yColumn]; 496 510 gradient[yColumn] += coefficient*solution[xColumn]; … … 520 534 delete [] gradient; 521 535 delete [] column; 536 cbcModel_->unlockThread(); 522 537 break; 523 538 } … … 634 649 delete [] lower2; 635 650 delete [] upper2; 636 if (isProvenOptimal())637 651 //if (isProvenOptimal()) 652 //writeMps("zz"); 638 653 } 639 654 } … … 652 667 CglTemporary * gen2 = dynamic_cast<CglTemporary *> (gen); 653 668 if (gen2) { 669 cbcModel_->lockThread(); 654 670 const double * solution = getColSolution(); 655 671 // add OA cut 656 double offset ;672 double offset=0.0; 657 673 int numberColumns = quadraticModel_->numberColumns(); 658 674 double * gradient = new double [numberColumns+1]; … … 677 693 int yColumn = obj->yColumn(); 678 694 if (xColumn!=yColumn) { 679 double coefficient = 2.0*obj->coefficient();695 double coefficient = /* 2.0* */obj->coefficient(); 680 696 gradient[xColumn] += coefficient*solution[yColumn]; 681 697 gradient[yColumn] += coefficient*solution[xColumn]; … … 710 726 delete [] gradient; 711 727 delete [] column; 728 cbcModel_->unlockThread(); 712 729 break; 713 730 } … … 735 752 //const double * columnLower = modelPtr_->columnLower(); 736 753 //const double * columnUpper = modelPtr_->columnUpper(); 754 cbcModel_->lockThread(); 737 755 for (int iNon=0;iNon<numberNonLinearRows_;iNon++) { 738 756 int iRow = rowNonLinear_[iNon]; … … 754 772 int yColumn = obj->yColumn(); 755 773 if (xColumn!=yColumn) { 756 double coefficient = 2.0*obj->coefficient();774 double coefficient = /* 2.0* */obj->coefficient(); 757 775 gradient[xColumn] += coefficient*solution[yColumn]; 758 776 gradient[yColumn] += coefficient*solution[xColumn]; … … 791 809 gen2->addCut(offset-1.0e-7,COIN_DBL_MAX,n,column,gradient); 792 810 } 811 cbcModel_->unlockThread(); 793 812 delete [] gradient; 794 813 delete [] column; … … 812 831 //------------------------------------------------------------------- 813 832 OsiSolverLink::OsiSolverLink () 814 : OsiClpSolverInterface()833 : CbcOsiSolver() 815 834 { 816 835 gutsOfDestructor(true); … … 893 912 */ 894 913 OsiSolverLink::OsiSolverLink ( CoinModel & coinModel) 895 : OsiClpSolverInterface()914 : CbcOsiSolver() 896 915 { 897 916 gutsOfDestructor(true); … … 899 918 } 900 919 // need bounds 901 static void fakeBounds(OsiSolverInterface * solver,int column,double maximumValue) 920 static void fakeBounds(OsiSolverInterface * solver,int column,double maximumValue, 921 CoinModel * model1, CoinModel * model2) 902 922 { 903 923 double lo = solver->getColLower()[column]; 904 if (lo<-maximumValue) 924 if (lo<-maximumValue) { 905 925 solver->setColLower(column,-maximumValue); 926 if (model1) 927 model1->setColLower(column,-maximumValue); 928 if (model2) 929 model2->setColLower(column,-maximumValue); 930 } 906 931 double up = solver->getColUpper()[column]; 907 if (up>maximumValue) 932 if (up>maximumValue) { 908 933 solver->setColUpper(column,maximumValue); 909 } 910 void OsiSolverLink::load ( CoinModel & coinModel, bool tightenBounds,int logLevel) 934 if (model1) 935 model1->setColUpper(column,maximumValue); 936 if (model2) 937 model2->setColUpper(column,maximumValue); 938 } 939 } 940 void OsiSolverLink::load ( CoinModel & coinModelOriginal, bool tightenBounds,int logLevel) 911 941 { 912 942 // first check and set up arrays 913 int numberColumns = coinModel .numberColumns();914 int numberRows = coinModel .numberRows();943 int numberColumns = coinModelOriginal.numberColumns(); 944 int numberRows = coinModelOriginal.numberRows(); 915 945 // List of nonlinear entries 916 946 int * which = new int[numberColumns]; … … 920 950 int numberErrors=0; 921 951 for (iColumn=0;iColumn<numberColumns;iColumn++) { 922 CoinModelLink triple=coinModel .firstInColumn(iColumn);952 CoinModelLink triple=coinModelOriginal.firstInColumn(iColumn); 923 953 bool linear=true; 924 954 int n=0; 925 955 // See if quadratic objective 926 const char * expr = coinModel .getColumnObjectiveAsString(iColumn);956 const char * expr = coinModelOriginal.getColumnObjectiveAsString(iColumn); 927 957 if (strcmp(expr,"Numeric")) { 928 958 linear=false; … … 930 960 while (triple.row()>=0) { 931 961 int iRow = triple.row(); 932 const char * expr = coinModel .getElementAsString(iRow,iColumn);962 const char * expr = coinModelOriginal.getElementAsString(iRow,iColumn); 933 963 if (strcmp(expr,"Numeric")) { 934 964 linear=false; 935 965 } 936 triple=coinModel .next(triple);966 triple=coinModelOriginal.next(triple); 937 967 n++; 938 968 } … … 944 974 if (!numberVariables_) { 945 975 delete [] which; 946 coinModel_ = coinModel ;976 coinModel_ = coinModelOriginal; 947 977 int nInt=0; 948 978 for (iColumn=0;iColumn<numberColumns;iColumn++) { … … 951 981 } 952 982 printf("There are %d integers\n",nInt); 953 loadFromCoinModel(coinModel ,true);983 loadFromCoinModel(coinModelOriginal,true); 954 984 OsiObject ** objects = new OsiObject * [nInt]; 955 985 nInt=0; … … 968 998 return; 969 999 } else { 970 coinModel_ = coinModel ;1000 coinModel_ = coinModelOriginal; 971 1001 // arrays for tightening bounds 972 1002 int * freeRow = new int [numberRows]; … … 1004 1034 ifFirst=false; 1005 1035 } 1006 coinModel .setObjective(iColumn,linearTerm);1036 coinModelOriginal.setObjective(iColumn,linearTerm); 1007 1037 } 1008 1038 } … … 1039 1069 ifFirst=false; 1040 1070 } 1041 coinModel .setElement(iRow,iColumn,linearTerm);1071 coinModelOriginal.setElement(iRow,iColumn,linearTerm); 1042 1072 } 1043 1073 triple=coinModel_.next(triple); … … 1053 1083 } 1054 1084 printf("There are %d bilinear and %d integers\n",nBi,nInt); 1055 loadFromCoinModel(coinModel,true); 1085 loadFromCoinModel(coinModelOriginal,true); 1086 CoinModel coinModel = coinModelOriginal; 1056 1087 if (tightenBounds&&numberColumns<100) { 1057 1088 // first fake bounds 1058 1089 for (iColumn=0;iColumn<numberColumns;iColumn++) { 1059 1090 if (tryColumn[iColumn]) { 1060 fakeBounds(this,iColumn,defaultBound_ );1091 fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_); 1061 1092 } 1062 1093 } … … 1091 1122 columnLower[iColumn]=value; 1092 1123 coinModel_.setColumnLower(iColumn,value); 1124 coinModel.setColumnLower(iColumn,value); 1093 1125 } 1094 1126 objective[iColumn]=-1.0; … … 1102 1134 columnUpper[iColumn]=value; 1103 1135 coinModel_.setColumnUpper(iColumn,value); 1136 coinModel.setColumnUpper(iColumn,value); 1104 1137 } 1105 1138 objective[iColumn]=0.0; … … 1127 1160 // add in objective as constraint 1128 1161 objectiveVariable_= numberColumns; 1129 objectiveRow_ = modelPtr_->numberRows();1130 addCol(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX,1.0);1162 objectiveRow_ = coinModel.numberRows(); 1163 coinModel.addColumn(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX,1.0); 1131 1164 int * column = new int[numberColumns+1]; 1132 1165 double * element = new double[numberColumns+1]; 1133 double * objective = modelPtr_->objective();1166 double * objective = coinModel.objectiveArray(); 1134 1167 int n=0; 1135 int starts[2]={0,0};1136 1168 for (int i=0;i<numberColumns;i++) { 1137 1169 if (objective[i]) { … … 1143 1175 column[n]=objectiveVariable_; 1144 1176 element[n++]=-1.0; 1145 starts[1]=n; 1146 double offset = - modelPtr_->objectiveOffset(); 1177 double offset = - coinModel.objectiveOffset(); 1147 1178 assert (!offset); // get sign right if happens 1148 modelPtr_->setObjectiveOffset(0.0);1179 coinModel.setObjectiveOffset(0.0); 1149 1180 double lowerBound = -COIN_DBL_MAX; 1150 addRows(1,starts,column,element,&lowerBound,&offset);1181 coinModel.addRow(n,column,element,lowerBound,offset); 1151 1182 delete [] column; 1152 1183 delete [] element; … … 1159 1190 double * sort = new double [nBi]; 1160 1191 nBi=nInt; 1192 const OsiObject ** justBi = const_cast<const OsiObject **> (objects+nInt); 1161 1193 for (iColumn=0;iColumn<numberColumns;iColumn++) { 1162 1194 if (quadraticObjective) … … 1166 1198 if (strcmp(expr,"Numeric")) { 1167 1199 // need bounds 1168 fakeBounds(this,iColumn,defaultBound_ );1200 fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_); 1169 1201 // value*x*y 1170 1202 char temp[20000]; … … 1179 1211 if (quadraticObjective) { 1180 1212 columnQuadratic[numberQuadratic]=jColumn; 1181 elementQuadratic[numberQuadratic++]=2.0*value; // convention 1213 if (jColumn==iColumn) 1214 elementQuadratic[numberQuadratic++]=2.0*value; // convention 1215 else 1216 elementQuadratic[numberQuadratic++]=1.0*value; // convention 1182 1217 } 1183 1218 // need bounds 1184 fakeBounds(this,jColumn,defaultBound_ );1219 fakeBounds(this,jColumn,defaultBound_,&coinModel,&coinModel_); 1185 1220 double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0; 1186 1221 if (meshI) … … 1209 1244 meshJ=0.0; 1210 1245 } 1211 OsiBiLinear * newObj = new OsiBiLinear( this,iColumn,jColumn,objectiveRow_,value,meshI,meshJ,1212 nBi-nInt, (const OsiObject **) (objects+nInt));1246 OsiBiLinear * newObj = new OsiBiLinear(&coinModel,iColumn,jColumn,objectiveRow_,value,meshI,meshJ, 1247 nBi-nInt,justBi); 1213 1248 newObj->setPriority(biLinearPriority_); 1214 1249 objects[nBi++] = newObj; … … 1240 1275 if (strcmp("Numeric",el)) { 1241 1276 // need bounds 1242 fakeBounds(this,iColumn,defaultBound_ );1277 fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_); 1243 1278 // value*x*y 1244 1279 char temp[20000]; … … 1252 1287 if (jColumn>=0) { 1253 1288 // need bounds 1254 fakeBounds(this,jColumn,defaultBound_ );1289 fakeBounds(this,jColumn,defaultBound_,&coinModel,&coinModel_); 1255 1290 double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0; 1256 1291 if (meshI) … … 1279 1314 meshJ=0.0; 1280 1315 } 1281 OsiBiLinear * newObj = new OsiBiLinear( this,iColumn,jColumn,iRow,value,meshI,meshJ,1282 nBi-nInt, (const OsiObject **) (objects+nInt));1316 OsiBiLinear * newObj = new OsiBiLinear(&coinModel,iColumn,jColumn,iRow,value,meshI,meshJ, 1317 nBi-nInt,justBi); 1283 1318 newObj->setPriority(biLinearPriority_); 1284 1319 objects[nBi++] = newObj; … … 1308 1343 stats[0],stats[1],stats[2],nBi-nInt-nDiff); 1309 1344 } 1345 // reload with all bilinear stuff 1346 loadFromCoinModel(coinModel,true); 1347 //exit(77); 1310 1348 nInt=0; 1311 1349 for (iColumn=0;iColumn<numberColumns;iColumn++) { … … 1458 1496 } 1459 1497 } 1498 delete [] which; 1499 if ((specialOptions2_&16)!=0) 1500 addTighterConstraints(); 1501 } 1502 // Add reformulated bilinear constraints 1503 void 1504 OsiSolverLink::addTighterConstraints() 1505 { 1506 // This is first attempt - for now get working on trimloss 1507 int numberW=0; 1508 int * xW = new int[numberObjects_]; 1509 int * yW = new int[numberObjects_]; 1510 // Points to firstlambda 1511 int * wW = new int[numberObjects_]; 1512 // Coefficient 1513 double * alphaW = new double[numberObjects_]; 1514 // Objects 1515 OsiBiLinear ** objW = new OsiBiLinear * [numberObjects_]; 1516 int numberColumns = getNumCols(); 1517 int firstLambda=numberColumns; 1518 // set up list (better to rethink and do properly as column ordered) 1519 int * list = new int[numberColumns]; 1520 memset(list,0,numberColumns*sizeof(int)); 1521 int i; 1522 for ( i =0;i<numberObjects_;i++) { 1523 OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]); 1524 if (obj) { 1525 //obj->setBranchingStrategy(4); // ***** temp 1526 objW[numberW]=obj; 1527 xW[numberW]=obj->xColumn(); 1528 yW[numberW]=obj->yColumn(); 1529 list[xW[numberW]]=1; 1530 list[yW[numberW]]=1; 1531 wW[numberW]=obj->firstLambda(); 1532 firstLambda = CoinMin(firstLambda,obj->firstLambda()); 1533 alphaW[numberW]=obj->coefficient(); 1534 //assert (alphaW[numberW]==1.0); // fix when occurs 1535 numberW++; 1536 } 1537 } 1538 int nList = 0; 1539 for (i=0;i<numberColumns;i++) { 1540 if (list[i]) 1541 list[nList++]=i; 1542 } 1543 // set up mark array 1544 char * mark = new char [firstLambda*firstLambda]; 1545 memset(mark,0,firstLambda*firstLambda); 1546 for (i=0;i<numberW;i++) { 1547 int x = xW[i]; 1548 int y = yW[i]; 1549 mark[x*firstLambda+y]=1; 1550 mark[y*firstLambda+x]=1; 1551 } 1552 int numberRows2 = originalRowCopy_->getNumRows(); 1553 int * addColumn = new int [numberColumns]; 1554 double * addElement = new double [numberColumns]; 1555 int * addW = new int [numberColumns]; 1556 assert (objectiveRow_<0); // fix when occurs 1557 for (int iRow=0;iRow<numberRows2;iRow++) { 1558 for (int iList=0;iList<nList;iList++) { 1559 int kColumn = list[iList]; 1560 const double * columnLower = getColLower(); 1561 //const double * columnUpper = getColUpper(); 1562 const double * rowLower = getRowLower(); 1563 const double * rowUpper = getRowUpper(); 1564 const CoinPackedMatrix * rowCopy = getMatrixByRow(); 1565 const double * element = rowCopy->getElements(); 1566 const int * column = rowCopy->getIndices(); 1567 const CoinBigIndex * rowStart = rowCopy->getVectorStarts(); 1568 const int * rowLength = rowCopy->getVectorLengths(); 1569 CoinBigIndex j; 1570 int numberElements = rowLength[iRow]; 1571 int n=0; 1572 for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) { 1573 int iColumn = column[j]; 1574 if (iColumn>=firstLambda) { 1575 // no good 1576 n=-1; 1577 break; 1578 } 1579 if (mark[iColumn*firstLambda+kColumn]) 1580 n++; 1581 } 1582 if (n==numberElements) { 1583 printf("can add row %d\n",iRow); 1584 assert (columnLower[kColumn]>=0); // might be able to fix 1585 n=0; 1586 for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) { 1587 int xColumn=kColumn; 1588 int yColumn = column[j]; 1589 int k; 1590 for (k=0;k<numberW;k++) { 1591 if ((xW[k]==yColumn&&yW[k]==xColumn)|| 1592 (yW[k]==yColumn&&xW[k]==xColumn)) 1593 break; 1594 } 1595 assert (k<numberW); 1596 if (xW[k]!=xColumn) { 1597 int temp=xColumn; 1598 xColumn=yColumn; 1599 yColumn=temp; 1600 } 1601 addW[n/4]=k; 1602 int start = wW[k]; 1603 double value = element[j]; 1604 for (int kk=0;kk<4;kk++) { 1605 // Dummy value 1606 addElement[n]= value; 1607 addColumn[n++]=start+kk; 1608 } 1609 } 1610 addColumn[n++] = kColumn; 1611 double lo = rowLower[iRow]; 1612 double up = rowUpper[iRow]; 1613 if (lo>-1.0e20) { 1614 // and tell object 1615 for (j=0;j<n-1;j+=4) { 1616 int iObject = addW[j/4]; 1617 objW[iObject]->addExtraRow(matrix_->getNumRows(),addElement[j]); 1618 } 1619 addElement[n-1]=-lo; 1620 if (lo==up) 1621 addRow(n,addColumn,addElement,0.0,0.0); 1622 else 1623 addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX); 1624 matrix_->appendRow(n,addColumn,addElement); 1625 } 1626 if (up<1.0e20&&up>lo) { 1627 // and tell object 1628 for (j=0;j<n-1;j+=4) { 1629 int iObject = addW[j/4]; 1630 objW[iObject]->addExtraRow(matrix_->getNumRows(),addElement[j]); 1631 } 1632 addElement[n-1]=-up; 1633 addRow(n,addColumn,addElement,-COIN_DBL_MAX,0.0); 1634 matrix_->appendRow(n,addColumn,addElement); 1635 } 1636 } 1637 } 1638 } 1460 1639 #if 0 1461 //int fake[10]={3,4,2,5,5,3,1,11,2,0}; 1462 int fake[10]={11,5,5,4,3,3,2,2,1,0}; 1463 for (int kk=0;kk<10;kk++) { 1464 setColUpper(kk,fake[kk]); 1465 setColLower(kk,fake[kk]); 1640 // possibly do bounds 1641 for (int iColumn=0;iColumn<firstLambda;iColumn++) { 1642 for (int iList=0;iList<nList;iList++) { 1643 int kColumn = list[iList]; 1644 const double * columnLower = getColLower(); 1645 const double * columnUpper = getColUpper(); 1646 if (mark[iColumn*firstLambda+kColumn]) { 1647 printf("can add column %d\n",iColumn); 1648 assert (columnLower[kColumn]>=0); // might be able to fix 1649 int xColumn=kColumn; 1650 int yColumn = iColumn; 1651 int k; 1652 for (k=0;k<numberW;k++) { 1653 if ((xW[k]==yColumn&&yW[k]==xColumn)|| 1654 (yW[k]==yColumn&&xW[k]==xColumn)) 1655 break; 1656 } 1657 assert (k<numberW); 1658 if (xW[k]!=xColumn) { 1659 int temp=xColumn; 1660 xColumn=yColumn; 1661 yColumn=temp; 1662 } 1663 int start = wW[k]; 1664 int n=0; 1665 for (int kk=0;kk<4;kk++) { 1666 // Dummy value 1667 addElement[n]= 1.0e-19; 1668 addColumn[n++]=start+kk; 1669 } 1670 // Tell object about this 1671 objW[k]->addExtraRow(matrix_->getNumRows(),1.0); 1672 addColumn[n++] = kColumn; 1673 double lo = columnLower[iColumn]; 1674 double up = columnUpper[iColumn]; 1675 if (lo>-1.0e20) { 1676 addElement[n-1]=-lo; 1677 if (lo==up) 1678 addRow(n,addColumn,addElement,0.0,0.0); 1679 else 1680 addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX); 1681 matrix_->appendRow(n,addColumn,addElement); 1682 } 1683 if (up<1.0e20&&up>lo) { 1684 addElement[n-1]=-up; 1685 addRow(n,addColumn,addElement,-COIN_DBL_MAX,0.0); 1686 matrix_->appendRow(n,addColumn,addElement); 1687 } 1688 } 1689 } 1466 1690 } 1467 1691 #endif 1468 delete [] which; 1692 delete [] xW; 1693 delete [] yW; 1694 delete [] wW; 1695 delete [] alphaW; 1696 delete [] addColumn; 1697 delete [] addElement; 1698 delete [] addW; 1699 delete [] mark; 1700 delete [] list; 1701 delete [] objW; 1469 1702 } 1470 1703 // Set all biLinear priorities on x-x variables … … 1491 1724 objNew->setYOtherSatisfied(oldSatisfied); 1492 1725 objNew->setYMeshSize(meshSize); 1493 objNew->setXYSatisfied(0. 5*meshSize);1726 objNew->setXYSatisfied(0.25*meshSize); 1494 1727 objNew->setPriority(value); 1495 1728 objNew->setBranchingStrategy(8); … … 1502 1735 delete [] newObject; 1503 1736 } 1737 /* Set options and priority on all or some biLinear variables 1738 1 - on I-I 1739 2 - on I-x 1740 4 - on x-x 1741 or combinations. 1742 -1 means leave (for priority value and strategy value) 1743 */ 1744 void 1745 OsiSolverLink::setBranchingStrategyOnVariables(int strategyValue, int priorityValue, 1746 int mode) 1747 { 1748 int i; 1749 for ( i =0;i<numberObjects_;i++) { 1750 OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]); 1751 if (obj) { 1752 bool change=false; 1753 if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0&&(mode&4)!=0) 1754 change=true; 1755 else if (((obj->xMeshSize()==1.0&&obj->yMeshSize()<1.0)|| 1756 (obj->xMeshSize()<1.0&&obj->yMeshSize()==1.0))&&(mode&2)!=0) 1757 change=true; 1758 else if (obj->xMeshSize()==1.0&&obj->yMeshSize()==1.0&&(mode&1)!=0) 1759 change=true; 1760 else if (obj->xMeshSize()>1.0||obj->yMeshSize()>1.0) 1761 abort(); 1762 if (change) { 1763 if (strategyValue>=0) 1764 obj->setBranchingStrategy(strategyValue); 1765 if (priorityValue>=0) 1766 obj->setPriority(priorityValue); 1767 } 1768 } 1769 } 1770 } 1771 1504 1772 // Say convex (should work it out) 1505 1773 void … … 1712 1980 } 1713 1981 model.initialSolve(); 1714 model.writeMps("bad.mps");1982 //model.writeMps("bad.mps"); 1715 1983 // redo number of columns 1716 1984 numberColumns = model.numberColumns(); … … 2079 2347 char temp[20]; 2080 2348 sprintf(temp,"pass%d.mps",iPass); 2081 model.writeMps(temp);2349 //model.writeMps(temp); 2082 2350 #ifdef CLP_DEBUG 2083 2351 if (model.status()) { … … 2138 2406 } 2139 2407 model.primal(1); 2140 model.writeMps("xx.mps");2408 //model.writeMps("xx.mps"); 2141 2409 iPass--; 2142 2410 goodMove=-1; … … 2270 2538 cbcModel->cutGenerator(6)->setTiming(true); 2271 2539 // For now - switch off most heuristics (because CglPreProcess is bad with QP) 2272 #if 02540 #if 1 2273 2541 CbcHeuristicFPump heuristicFPump(*cbcModel); 2274 2542 heuristicFPump.setWhen(13); … … 2291 2559 2292 2560 CbcRounding rounding(*cbcModel); 2561 rounding.setHeuristicName("rounding"); 2293 2562 cbcModel->addHeuristic(&rounding); 2294 2563 … … 2329 2598 // if convex 2330 2599 if ((specialOptions2()&4)!=0) { 2600 if (cbcModel_) 2601 cbcModel_->lockThread(); 2331 2602 // add OA cut 2332 2603 double offset; … … 2350 2621 delete [] gradient; 2351 2622 delete [] column; 2623 if (cbcModel_) 2624 cbcModel_->unlockThread(); 2352 2625 } 2353 2626 delete qp; … … 2369 2642 double * solution = CoinCopyOfArray(temp->primalColumnSolution(),numberColumns); 2370 2643 delete temp; 2644 if (mode==0) { 2645 return solution; 2646 } else if (mode==2) { 2647 const double * lower = getColLower(); 2648 const double * upper = getColUpper(); 2649 for (int iObject =0;iObject<numberObjects_;iObject++) { 2650 OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]); 2651 if (obj&&(obj->priority()<biLinearPriority_||biLinearPriority_<=0)) { 2652 int iColumn = obj->columnNumber(); 2653 double value = solution[iColumn]; 2654 value = floor(value+0.5); 2655 if (fabs(value-solution[iColumn])>0.01) { 2656 setColLower(iColumn,CoinMax(lower[iColumn],value-CoinMax(defaultBound_,0.0))); 2657 setColUpper(iColumn,CoinMin(upper[iColumn],value+CoinMax(defaultBound_,1.0))); 2658 } else { 2659 // could fix to integer 2660 setColLower(iColumn,CoinMax(lower[iColumn],value-CoinMax(defaultBound_,0.0))); 2661 setColUpper(iColumn,CoinMin(upper[iColumn],value+CoinMax(defaultBound_,0.0))); 2662 } 2663 } 2664 } 2665 return solution; 2666 } 2371 2667 OsiClpSolverInterface newSolver; 2372 2668 if (mode==1) { … … 2376 2672 tempModel = coinModel_; 2377 2673 // solve modified problem 2674 char * mark = new char[numberColumns]; 2675 memset(mark,0,numberColumns); 2378 2676 for (int iObject =0;iObject<numberObjects_;iObject++) { 2379 2677 OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]); … … 2383 2681 value = ceil(value-1.0e-7); 2384 2682 tempModel.associateElement(coinModel_.columnName(iColumn),value); 2385 } 2386 } 2683 mark[iColumn]=1; 2684 } 2685 OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (object_[iObject]); 2686 if (objB) { 2687 // if one or both continuous then fix one 2688 if (objB->xMeshSize()<1.0) { 2689 int xColumn = objB->xColumn(); 2690 double value = solution[xColumn]; 2691 tempModel.associateElement(coinModel_.columnName(xColumn),value); 2692 mark[xColumn]=1; 2693 } else if (objB->yMeshSize()<1.0) { 2694 int yColumn = objB->yColumn(); 2695 double value = solution[yColumn]; 2696 tempModel.associateElement(coinModel_.columnName(yColumn),value); 2697 mark[yColumn]=1; 2698 } 2699 } 2700 } 2701 CoinModel * reOrdered = tempModel.reorder(mark); 2702 assert (reOrdered); 2703 tempModel=*reOrdered; 2704 delete reOrdered; 2705 delete [] mark; 2387 2706 newSolver.loadFromCoinModel(tempModel,true); 2388 2707 for (int iObject =0;iObject<numberObjects_;iObject++) { … … 2394 2713 newSolver.setColLower(iColumn,value); 2395 2714 newSolver.setColUpper(iColumn,value); 2715 } 2716 OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (object_[iObject]); 2717 if (objB) { 2718 // if one or both continuous then fix one 2719 if (objB->xMeshSize()<1.0) { 2720 int xColumn = objB->xColumn(); 2721 double value = solution[xColumn]; 2722 newSolver.setColLower(xColumn,value); 2723 newSolver.setColUpper(xColumn,value); 2724 } else if (objB->yMeshSize()<1.0) { 2725 int yColumn = objB->yColumn(); 2726 double value = solution[yColumn]; 2727 newSolver.setColLower(yColumn,value); 2728 newSolver.setColUpper(yColumn,value); 2729 } 2396 2730 } 2397 2731 } … … 2491 2825 clpModel->setLogLevel(saveLogLevel); 2492 2826 returnCode=-1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl; 2493 clpModel->writeMps("infeas2.mps");2827 //clpModel->writeMps("infeas2.mps"); 2494 2828 } else { 2495 2829 clpModel->setLogLevel(saveLogLevel); … … 2710 3044 { 2711 3045 assert (copyData); 2712 return new OsiSolverLink(*this); 3046 OsiSolverLink * newModel = new OsiSolverLink(*this); 3047 return newModel; 2713 3048 } 2714 3049 … … 2719 3054 OsiSolverLink::OsiSolverLink ( 2720 3055 const OsiSolverLink & rhs) 2721 : OsiClpSolverInterface(rhs)3056 : CbcOsiSolver(rhs) 2722 3057 { 2723 3058 gutsOfDestructor(true); … … 2743 3078 if (this != &rhs) { 2744 3079 gutsOfDestructor(); 2745 OsiClpSolverInterface::operator=(rhs);3080 CbcOsiSolver::operator=(rhs); 2746 3081 gutsOfCopy(rhs); 2747 3082 } … … 2771 3106 convex_ = NULL; 2772 3107 whichNonLinear_ = NULL; 2773 cbcModel_ = NULL;2774 3108 info_ = NULL; 2775 3109 fixVariables_=NULL; … … 2801 3135 biLinearPriority_ = rhs.biLinearPriority_; 2802 3136 numberFix_ = rhs.numberFix_; 2803 cbcModel_ = rhs.cbcModel_;2804 3137 if (numberVariables_) { 2805 3138 if (rhs.matrix_) … … 3094 3427 } 3095 3428 } 3096 newSolver.writeMps("xx");3429 //newSolver.writeMps("xx"); 3097 3430 CbcModel model(newSolver); 3098 3431 // Now do requested saves and modifications … … 3318 3651 setRowBounds(i,-COIN_DBL_MAX,COIN_DBL_MAX); 3319 3652 initialSolve(); 3320 if (!isProvenOptimal())3321 3653 //if (!isProvenOptimal()) 3654 //getModelPtr()->writeMps("bad.mps"); 3322 3655 if (isProvenOptimal()) { 3323 3656 delete [] bestSolution_; … … 4302 4635 xyRow_(-1), 4303 4636 convexity_(-1), 4637 numberExtraRows_(0), 4638 multiplier_(NULL), 4639 extraRow_(NULL), 4304 4640 chosen_(-1) 4305 4641 { … … 4330 4666 xyRow_(xyRow), 4331 4667 convexity_(-1), 4668 numberExtraRows_(0), 4669 multiplier_(NULL), 4670 extraRow_(NULL), 4332 4671 chosen_(-1) 4333 4672 { … … 4517 4856 } 4518 4857 } 4858 // Useful constructor 4859 OsiBiLinear::OsiBiLinear (CoinModel * coinModel, int xColumn, 4860 int yColumn, int xyRow, double coefficient, 4861 double xMesh, double yMesh, 4862 int numberExistingObjects,const OsiObject ** objects ) 4863 : OsiObject2(), 4864 coefficient_(coefficient), 4865 xMeshSize_(xMesh), 4866 yMeshSize_(yMesh), 4867 xSatisfied_(1.0e-6), 4868 ySatisfied_(1.0e-6), 4869 xOtherSatisfied_(0.0), 4870 yOtherSatisfied_(0.0), 4871 xySatisfied_(1.0e-6), 4872 xyBranchValue_(0.0), 4873 xColumn_(xColumn), 4874 yColumn_(yColumn), 4875 firstLambda_(-1), 4876 branchingStrategy_(0), 4877 boundType_(0), 4878 xRow_(-1), 4879 yRow_(-1), 4880 xyRow_(xyRow), 4881 convexity_(-1), 4882 numberExtraRows_(0), 4883 multiplier_(NULL), 4884 extraRow_(NULL), 4885 chosen_(-1) 4886 { 4887 double columnLower[4]; 4888 double columnUpper[4]; 4889 double objective[4]; 4890 double rowLower[3]; 4891 double rowUpper[3]; 4892 CoinBigIndex starts[5]; 4893 int index[16]; 4894 double element[16]; 4895 int i; 4896 starts[0]=0; 4897 // rows 4898 int numberRows = coinModel->numberRows(); 4899 // convexity 4900 rowLower[0]=1.0; 4901 rowUpper[0]=1.0; 4902 convexity_ = numberRows; 4903 starts[1]=0; 4904 // x 4905 rowLower[1]=0.0; 4906 rowUpper[1]=0.0; 4907 index[0]=xColumn_; 4908 element[0]=-1.0; 4909 xRow_ = numberRows+1; 4910 starts[2]=1; 4911 int nAdd=2; 4912 if (xColumn_!=yColumn_) { 4913 rowLower[2]=0.0; 4914 rowUpper[2]=0.0; 4915 index[1]=yColumn; 4916 element[1]=-1.0; 4917 nAdd=3; 4918 yRow_ = numberRows+2; 4919 starts[3]=2; 4920 } else { 4921 yRow_=-1; 4922 branchingStrategy_=1; 4923 } 4924 // may be objective 4925 assert (xyRow_>=-1); 4926 for (i=0;i<nAdd;i++) { 4927 CoinBigIndex iStart = starts[i]; 4928 coinModel->addRow(starts[i+1]-iStart,index+iStart,element+iStart,rowLower[i],rowUpper[i]); 4929 } 4930 int n=0; 4931 // order is LxLy, LxUy, UxLy and UxUy 4932 firstLambda_ = coinModel->numberColumns(); 4933 // bit sloppy as theoretically could be infeasible but otherwise need to do more work 4934 double xB[2]; 4935 double yB[2]; 4936 const double * lower = coinModel->columnLowerArray(); 4937 const double * upper = coinModel->columnUpperArray(); 4938 xB[0]=lower[xColumn_]; 4939 xB[1]=upper[xColumn_]; 4940 yB[0]=lower[yColumn_]; 4941 yB[1]=upper[yColumn_]; 4942 if (xMeshSize_!=floor(xMeshSize_)) { 4943 // not integral 4944 xSatisfied_ = CoinMax(xSatisfied_,0.51*xMeshSize_); 4945 if (!yMeshSize_) { 4946 xySatisfied_ = CoinMax(xySatisfied_,xSatisfied_*CoinMax(fabs(yB[0]),fabs(yB[1]))); 4947 } 4948 } 4949 if (yMeshSize_!=floor(yMeshSize_)) { 4950 // not integral 4951 ySatisfied_ = CoinMax(ySatisfied_,0.51*yMeshSize_); 4952 if (!xMeshSize_) { 4953 xySatisfied_ = CoinMax(xySatisfied_,ySatisfied_*CoinMax(fabs(xB[0]),fabs(xB[1]))); 4954 } 4955 } 4956 // adjust 4957 double distance; 4958 double steps; 4959 if (xMeshSize_) { 4960 distance = xB[1]-xB[0]; 4961 steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_); 4962 distance = xB[0]+xMeshSize_*steps; 4963 if (fabs(xB[1]-distance)>xSatisfied_) { 4964 printf("bad x mesh %g %g %g -> %g\n",xB[0],xMeshSize_,xB[1],distance); 4965 //double newValue = CoinMax(fabs(xB[1]-distance),xMeshSize_); 4966 //printf("xSatisfied increased to %g\n",newValue); 4967 //xSatisfied_ = newValue; 4968 //xB[1]=distance; 4969 //coinModel->setColUpper(xColumn_,distance); 4970 } 4971 } 4972 if (yMeshSize_) { 4973 distance = yB[1]-yB[0]; 4974 steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_); 4975 distance = yB[0]+yMeshSize_*steps; 4976 if (fabs(yB[1]-distance)>ySatisfied_) { 4977 printf("bad y mesh %g %g %g -> %g\n",yB[0],yMeshSize_,yB[1],distance); 4978 //double newValue = CoinMax(fabs(yB[1]-distance),yMeshSize_); 4979 //printf("ySatisfied increased to %g\n",newValue); 4980 //ySatisfied_ = newValue; 4981 //yB[1]=distance; 4982 //coinModel->setColUpper(yColumn_,distance); 4983 } 4984 } 4985 for (i=0;i<4;i++) { 4986 double x = (i<2) ? xB[0] : xB[1]; 4987 double y = ((i&1)==0) ? yB[0] : yB[1]; 4988 columnLower[i]=0.0; 4989 columnUpper[i]=2.0; 4990 objective[i]=0.0; 4991 double value; 4992 // xy 4993 value=coefficient_*x*y; 4994 if (xyRow_>=0) { 4995 if (fabs(value)<1.0e-19) 4996 value = 1.0e-19; 4997 element[n]=value; 4998 index[n++]=xyRow_; 4999 } else { 5000 objective[i]=value; 5001 } 5002 // convexity 5003 value=1.0; 5004 element[n]=value; 5005 index[n++]=0+numberRows; 5006 // x 5007 value=x; 5008 if (fabs(value)<1.0e-19) 5009 value = 1.0e-19; 5010 element[n]=value; 5011 index[n++]=1+numberRows; 5012 if (xColumn_!=yColumn_) { 5013 // y 5014 value=y; 5015 if (fabs(value)<1.0e-19) 5016 value = 1.0e-19; 5017 element[n]=value; 5018 index[n++]=2+numberRows; 5019 } 5020 starts[i+1]=n; 5021 } 5022 for (i=0;i<4;i++) { 5023 CoinBigIndex iStart = starts[i]; 5024 coinModel->addColumn(starts[i+1]-iStart,index+iStart,element+iStart,columnLower[i], 5025 columnUpper[i],objective[i]); 5026 } 5027 // At least one has to have a mesh 5028 if (!xMeshSize_&&(!yMeshSize_||yRow_<0)) { 5029 printf("one of x and y must have a mesh size\n"); 5030 abort(); 5031 } else if (yRow_>=0) { 5032 if (!xMeshSize_) 5033 branchingStrategy_ = 2; 5034 else if (!yMeshSize_) 5035 branchingStrategy_ = 1; 5036 } 5037 // Now add constraints to link in x and or y to existing ones. 5038 bool xDone=false; 5039 bool yDone=false; 5040 // order is LxLy, LxUy, UxLy and UxUy 5041 for (i=numberExistingObjects-1;i>=0;i--) { 5042 const OsiObject * obj = objects[i]; 5043 const OsiBiLinear * obj2 = 5044 dynamic_cast <const OsiBiLinear *>(obj) ; 5045 if (obj2) { 5046 if (xColumn_==obj2->xColumn_&&!xDone) { 5047 // make sure y equal 5048 double rhs=0.0; 5049 CoinBigIndex starts[2]; 5050 int index[4]; 5051 double element[4]= {1.0,1.0,-1.0,-1.0}; 5052 starts[0]=0; 5053 starts[1]=4; 5054 index[0]=firstLambda_+0; 5055 index[1]=firstLambda_+1; 5056 index[2]=obj2->firstLambda_+0; 5057 index[3]=obj2->firstLambda_+1; 5058 coinModel->addRow(4,index,element,rhs,rhs); 5059 xDone=true; 5060 } 5061 if (yColumn_==obj2->yColumn_&&yRow_>=0&&!yDone) { 5062 // make sure x equal 5063 double rhs=0.0; 5064 CoinBigIndex starts[2]; 5065 int index[4]; 5066 double element[4]= {1.0,1.0,-1.0,-1.0}; 5067 starts[0]=0; 5068 starts[1]=4; 5069 index[0]=firstLambda_+0; 5070 index[1]=firstLambda_+2; 5071 index[2]=obj2->firstLambda_+0; 5072 index[3]=obj2->firstLambda_+2; 5073 coinModel->addRow(4,index,element,rhs,rhs); 5074 yDone=true; 5075 } 5076 } 5077 } 5078 } 4519 5079 4520 5080 // Copy constructor … … 4539 5099 xyRow_(rhs.xyRow_), 4540 5100 convexity_(rhs.convexity_), 5101 numberExtraRows_(rhs.numberExtraRows_), 5102 multiplier_(NULL), 5103 extraRow_(NULL), 4541 5104 chosen_(rhs.chosen_) 4542 5105 { 5106 if (numberExtraRows_) { 5107 multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_); 5108 extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_); 5109 } 4543 5110 } 4544 5111 … … 4574 5141 xyRow_ = rhs.xyRow_; 4575 5142 convexity_ = rhs.convexity_; 5143 numberExtraRows_ = rhs.numberExtraRows_; 5144 delete [] multiplier_; 5145 delete [] extraRow_; 5146 if (numberExtraRows_) { 5147 multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_); 5148 extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_); 5149 } else { 5150 multiplier_ = NULL; 5151 extraRow_ = NULL; 5152 } 4576 5153 chosen_ = rhs.chosen_; 4577 5154 } … … 4582 5159 OsiBiLinear::~OsiBiLinear () 4583 5160 { 4584 } 4585 5161 delete [] multiplier_; 5162 delete [] extraRow_; 5163 } 5164 // Adds in data for extra row with variable coefficients 5165 void 5166 OsiBiLinear::addExtraRow(int row, double multiplier) 5167 { 5168 int * tempI = new int [numberExtraRows_+1]; 5169 double * tempD = new double [numberExtraRows_+1]; 5170 memcpy(tempI,extraRow_,numberExtraRows_*sizeof(int)); 5171 memcpy(tempD,multiplier_,numberExtraRows_*sizeof(double)); 5172 tempI[numberExtraRows_]=row; 5173 tempD[numberExtraRows_]=multiplier; 5174 if (numberExtraRows_) 5175 assert (row>tempI[numberExtraRows_-1]); 5176 numberExtraRows_++; 5177 delete [] extraRow_; 5178 extraRow_ = tempI; 5179 delete [] multiplier_; 5180 multiplier_ = tempD; 5181 } 5182 static bool testCoarse=true; 4586 5183 // Infeasibility - large is 0.5 4587 5184 double … … 4657 5254 double steps; 4658 5255 bool xSatisfied; 4659 double xNew ;5256 double xNew=xB[0]; 4660 5257 if (xMeshSize_) { 4661 5258 if (x<0.5*(xB[0]+xB[1])) { … … 4673 5270 } 4674 5271 // but if first coarse grid then only if gap small 4675 if ( false&&(branchingStrategy_&8)!=0&&xSatisfied&&5272 if (testCoarse&&(branchingStrategy_&8)!=0&&xSatisfied&& 4676 5273 xB[1]-xB[0]>=xMeshSize_) { 4677 5274 // but allow if fine grid would allow … … 4687 5284 } 4688 5285 bool ySatisfied; 4689 double yNew ;5286 double yNew=yB[0]; 4690 5287 if (yMeshSize_) { 4691 5288 if (y<0.5*(yB[0]+yB[1])) { … … 4703 5300 } 4704 5301 // but if first coarse grid then only if gap small 4705 if ( false&&(branchingStrategy_&8)!=0&&ySatisfied&&5302 if (testCoarse&&(branchingStrategy_&8)!=0&&ySatisfied&& 4706 5303 yB[1]-yB[0]>=yMeshSize_) { 4707 5304 // but allow if fine grid would allow … … 4765 5362 xyLambda /= coefficient_; 4766 5363 } 5364 if (0) { 5365 // only true with positive values 5366 // see if all convexification constraints OK with true 5367 assert (xyTrue+1.0e-5>xB[0]*y+yB[0]*x - xB[0]*yB[0]); 5368 assert (xyTrue+1.0e-5>xB[1]*y+yB[1]*x - xB[1]*yB[1]); 5369 assert (xyTrue-1.0e-5<xB[1]*y+yB[0]*x - xB[1]*yB[0]); 5370 assert (xyTrue-1.0e-5<xB[0]*y+yB[1]*x - xB[0]*yB[1]); 5371 // see if all convexification constraints OK with lambda version 5372 #if 1 5373 assert (xyLambda+1.0e-5>xB[0]*y+yB[0]*x - xB[0]*yB[0]); 5374 assert (xyLambda+1.0e-5>xB[1]*y+yB[1]*x - xB[1]*yB[1]); 5375 assert (xyLambda-1.0e-5<xB[1]*y+yB[0]*x - xB[1]*yB[0]); 5376 assert (xyLambda-1.0e-5<xB[0]*y+yB[1]*x - xB[0]*yB[1]); 5377 #endif 5378 // see if other bound stuff true 5379 assert (xyLambda+1.0e-5>xB[0]*y); 5380 assert (xyLambda+1.0e-5>yB[0]*x); 5381 assert (xyLambda-1.0e-5<xB[1]*y); 5382 assert (xyLambda-1.0e-5<yB[1]*x); 5383 #define SIZE 2 5384 if (yColumn_==xColumn_+SIZE) { 5385 #if SIZE==6 5386 double bMax = 2200.0; 5387 double bMin = bMax - 100.0; 5388 double b[] = {330.0,360.0,380.0,430.0,490.0,530.0}; 5389 #elif SIZE==2 5390 double bMax = 1900.0; 5391 double bMin = bMax - 200.0; 5392 double b[] = {460.0,570.0}; 5393 #else 5394 abort(); 5395 #endif 5396 double sum =0.0; 5397 double sum2 =0.0; 5398 int m=xColumn_; 5399 double x = info->solution_[m]; 5400 double xB[2]; 5401 double yB[2]; 5402 xB[0]=info->lower_[m]; 5403 xB[1]=info->upper_[m]; 5404 for (int i=0;i<SIZE*SIZE;i+=SIZE) { 5405 int n = i+SIZE+m; 5406 double y = info->solution_[n]; 5407 yB[0]=info->lower_[n]; 5408 yB[1]=info->upper_[n]; 5409 int firstLambda=SIZE*SIZE+2*SIZE+4*i+4*m; 5410 double xyLambda=0.0; 5411 for (int j=0;j<4;j++) { 5412 int iX = j>>1; 5413 int iY = j&1; 5414 xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda+j]; 5415 } 5416 sum += xyLambda*b[i/SIZE]; 5417 double xyTrue = x*y; 5418 sum2 += xyTrue*b[i/SIZE]; 5419 } 5420 if (sum>bMax*x+1.0e-5||sum<bMin*x-1.0e-5) { 5421 //if (sum<bMax*x+1.0e-5&&sum>bMin*x-1.0e-5) { 5422 printf("bmin*x %g b*w %g bmax*x %g (true) %g\n",bMin*x,sum,bMax*x,sum2); 5423 printf("m %d lb %g value %g up %g\n", 5424 m,xB[0],x,xB[1]); 5425 sum=0.0; 5426 for (int i=0;i<SIZE*SIZE;i+=SIZE) { 5427 int n = i+SIZE+m; 5428 double y = info->solution_[n]; 5429 yB[0]=info->lower_[n]; 5430 yB[1]=info->upper_[n]; 5431 printf("n %d lb %g value %g up %g\n", 5432 n,yB[0],y,yB[1]); 5433 int firstLambda=SIZE*SIZE+2*SIZE+4*i+m*4; 5434 double xyLambda=0.0; 5435 for (int j=0;j<4;j++) { 5436 int iX = j>>1; 5437 int iY = j&1; 5438 xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda+j]; 5439 printf("j %d l %d new xylambda %g ",j,firstLambda+j,xyLambda); 5440 } 5441 sum += xyLambda*b[i/SIZE]; 5442 printf(" - sum now %g\n",sum); 5443 } 5444 } 5445 if (sum2>bMax*x+1.0e-5||sum2<bMin*x-1.0e-5) { 5446 printf("bmin*x %g b*x*y %g bmax*x %g (estimate) %g\n",bMin*x,sum2,bMax*x,sum); 5447 printf("m %d lb %g value %g up %g\n", 5448 m,xB[0],x,xB[1]); 5449 sum2=0.0; 5450 for (int i=0;i<SIZE*SIZE;i+=SIZE) { 5451 int n = i+SIZE+m; 5452 double y = info->solution_[n]; 5453 yB[0]=info->lower_[n]; 5454 yB[1]=info->upper_[n]; 5455 printf("n %d lb %g value %g up %g\n", 5456 n,yB[0],y,yB[1]); 5457 double xyTrue = x*y; 5458 sum2 += xyTrue*b[i/SIZE]; 5459 printf("xyTrue %g - sum now %g\n",xyTrue,sum2); 5460 } 5461 } 5462 } 5463 } 4767 5464 // If pseudo shadow prices then see what would happen 4768 5465 //double pseudoEstimate = 0.0; … … 4781 5478 // == row so move x and y not xy 4782 5479 } 5480 } 5481 if ((branchingStrategy_&16)!=0) { 5482 // always treat as satisfied!! 5483 xSatisfied=true; 5484 ySatisfied=true; 5485 xyTrue=xyLambda; 4783 5486 } 4784 5487 if ( !xSatisfied) { … … 4852 5555 } 4853 5556 } 5557 if (testCoarse&&(branchingStrategy_&8)!=0&&xB[1]-xB[0]<1.0001*xSatisfied_&& 5558 yB[1]-yB[0]<1.0001*ySatisfied_) 5559 feasible=true; 4854 5560 if (feasible) { 4855 5561 if (xB[1]-xB[0]>=xSatisfied_&&xMeshSize_) { … … 4938 5644 } 4939 5645 whichWay=whichWay_; 5646 //if (infeasibility_&&priority_==10) 5647 //printf("x %d %g %g %g, y %d %g %g %g\n",xColumn_,xB[0],x,xB[1],yColumn_,yB[0],y,yB[1]); 4940 5648 return infeasibility_; 4941 5649 } 4942 5650 // Sets infeasibility and other when pseudo shadow prices 5651 void 5652 OsiBiLinear::getPseudoShadow(const OsiBranchingInformation * info) 5653 { 5654 // order is LxLy, LxUy, UxLy and UxUy 5655 double xB[2]; 5656 double yB[2]; 5657 xB[0]=info->lower_[xColumn_]; 5658 xB[1]=info->upper_[xColumn_]; 5659 yB[0]=info->lower_[yColumn_]; 5660 yB[1]=info->upper_[yColumn_]; 5661 double x = info->solution_[xColumn_]; 5662 x = CoinMax(x,xB[0]); 5663 x = CoinMin(x,xB[1]); 5664 double y = info->solution_[yColumn_]; 5665 y = CoinMax(y,yB[0]); 5666 y = CoinMin(y,yB[1]); 5667 int j; 5668 double xyTrue = x*y; 5669 double xyLambda = 0.0; 5670 if ((branchingStrategy_&4)==0) { 5671 for (j=0;j<4;j++) { 5672 int iX = j>>1; 5673 int iY = j&1; 5674 xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda_+j]; 5675 } 5676 } else { 5677 if (xyRow_>=0) { 5678 const double * element = info->elementByColumn_; 5679 const int * row = info->row_; 5680 const CoinBigIndex * columnStart = info->columnStart_; 5681 const int * columnLength = info->columnLength_; 5682 for (j=0;j<4;j++) { 5683 int iColumn = firstLambda_+j; 5684 int iStart = columnStart[iColumn]; 5685 int iEnd = iStart + columnLength[iColumn]; 5686 int k=iStart; 5687 double sol = info->solution_[iColumn]; 5688 for (;k<iEnd;k++) { 5689 if (xyRow_==row[k]) 5690 xyLambda += element[k]*sol; 5691 } 5692 } 5693 } else { 5694 // objective 5695 const double * objective = info->objective_; 5696 for (j=0;j<4;j++) { 5697 int iColumn = firstLambda_+j; 5698 double sol = info->solution_[iColumn]; 5699 xyLambda += objective[iColumn]*sol; 5700 } 5701 } 5702 xyLambda /= coefficient_; 5703 } 5704 assert (info->defaultDual_>=0.0); 5705 // If we move to xy then we move by coefficient * (xyTrue-xyLambda) on row xyRow_ 5706 double movement = xyTrue-xyLambda; 5707 infeasibility_=0.0; 5708 const double * pi = info->pi_; 5709 const double * activity = info->rowActivity_; 5710 const double * lower = info->rowLower_; 5711 const double * upper = info->rowUpper_; 5712 double tolerance = info->primalTolerance_; 5713 double direction = info->direction_; 5714 if (xyRow_>=0) { 5715 assert (!boundType_); 5716 if (lower[xyRow_]<-1.0e20) 5717 assert (pi[xyRow_]<=1.0e-3); 5718 if (upper[xyRow_]>1.0e20) 5719 assert (pi[xyRow_]>=-1.0e-3); 5720 double valueP = pi[xyRow_]*direction; 5721 // if move makes infeasible then make at least default 5722 double newValue = activity[xyRow_] + movement*coefficient_; 5723 if (newValue>upper[xyRow_]+tolerance||newValue<lower[xyRow_]-tolerance) 5724 infeasibility_ += fabs(movement*coefficient_)*CoinMax(fabs(valueP),info->defaultDual_); 5725 } else { 5726 // objective 5727 assert (movement>-1.0e-7); 5728 infeasibility_ += movement; 5729 } 5730 for (int i=0;i<numberExtraRows_;i++) { 5731 int iRow = extraRow_[i]; 5732 if (lower[iRow]<-1.0e20) 5733 assert (pi[iRow]<=1.0e-3); 5734 if (upper[iRow]>1.0e20) 5735 assert (pi[iRow]>=-1.0e-3); 5736 double valueP = pi[iRow]*direction; 5737 // if move makes infeasible then make at least default 5738 double newValue = activity[iRow] + movement*multiplier_[i]; 5739 if (newValue>upper[iRow]+tolerance||newValue<lower[iRow]-tolerance) 5740 infeasibility_ += fabs(movement*multiplier_[i])*CoinMax(fabs(valueP),info->defaultDual_); 5741 } 5742 if (infeasibility_<1.0e-7) 5743 infeasibility_=0.0; 5744 otherInfeasibility_ = CoinMax(1.0e-12,infeasibility_*10.0); 5745 } 4943 5746 // This looks at solution and sets bounds to contain solution 4944 5747 double … … 5476 6279 const int * row = matrix->getIndices(); 5477 6280 const CoinBigIndex * columnStart = matrix->getVectorStarts(); 5478 //const int * columnLength = matrix->getVectorLengths();6281 const int * columnLength = matrix->getVectorLengths(); 5479 6282 // order is LxLy, LxUy, UxLy and UxUy 5480 6283 double xB[2]; … … 5497 6300 double y = yB[iY]; 5498 6301 CoinBigIndex k = columnStart[j+firstLambda_]; 6302 CoinBigIndex last = k+columnLength[j+firstLambda_]; 5499 6303 double value; 5500 6304 // xy … … 5531 6335 element[k++]=value; 5532 6336 numberUpdated++; 6337 } 6338 // Do extra rows 6339 for (int i=0;i<numberExtraRows_;i++) { 6340 int iRow = extraRow_[i]; 6341 for (;k<last;k++) { 6342 if (row[k]==iRow) 6343 break; 6344 } 6345 assert (k<last); 6346 element[k++] = x*y*multiplier_[i]; 5533 6347 } 5534 6348 } … … 6099 6913 if (newUp>upper[iRow]+tolerance||newUp<lower[iRow]-tolerance) 6100 6914 u = CoinMax(u,info->defaultDual_); 6101 upEstimate += u*upMovement ;6915 upEstimate += u*upMovement*fabs(el2); 6102 6916 // if down makes infeasible then make at least default 6103 6917 double newDown = activity[iRow] - downMovement*el2; 6104 6918 if (newDown>upper[iRow]+tolerance||newDown<lower[iRow]-tolerance) 6105 6919 d = CoinMax(d,info->defaultDual_); 6106 downEstimate += d*downMovement ;6920 downEstimate += d*downMovement*fabs(el2); 6107 6921 } 6108 6922 if (downEstimate>=upEstimate) { … … 6692 7506 return nelCreate; 6693 7507 } 6694 #endif6695 7508 #include "ClpConstraint.hpp" 6696 7509 #include "ClpConstraintLinear.hpp" … … 6705 7518 int mode) 6706 7519 { 7520 #ifdef COIN_HAS_ASL 7521 // matrix etc will be changed 7522 CoinModel coinModel2 = coinModel; 7523 if (coinModel2.moreInfo()) { 7524 // for now just ampl objective 7525 ClpSimplex * model = new ClpSimplex(); 7526 model->loadProblem(coinModel2); 7527 int numberConstraints; 7528 ClpConstraint ** constraints=NULL; 7529 int type = model->loadNonLinear(coinModel2.moreInfo(), 7530 numberConstraints,constraints); 7531 if (type==1||type==3) { 7532 int returnCode = model->nonlinearSLP(numberPasses,deltaTolerance); 7533 assert (!returnCode); 7534 } else if (type==2||type==4) { 7535 int returnCode = model->nonlinearSLP(numberConstraints,constraints, 7536 numberPasses,deltaTolerance); 7537 assert (!returnCode); 7538 } else { 7539 printf("error or linear - fix %d\n",type); 7540 } 7541 //exit(66); 7542 return model; 7543 } 6707 7544 // first check and set up arrays 6708 7545 int numberColumns = coinModel.numberColumns(); 6709 7546 int numberRows = coinModel.numberRows(); 6710 7547 // List of nonlinear rows 6711 int * which = new int[numberRows +1];7548 int * which = new int[numberRows]; 6712 7549 bool testLinear=false; 6713 7550 int numberConstraints=0; 6714 7551 int iColumn; 6715 // matrix etc will be changed6716 CoinModel coinModel2 = coinModel;6717 7552 bool linearObjective=true; 6718 7553 int maximumQuadraticElements=0; … … 6749 7584 for (iColumn=0;iColumn<numberColumns;iColumn++) 6750 7585 coinModel2.setObjective(iColumn,0.0); 6751 which[numberConstraints++]=-1;6752 7586 } 6753 7587 int iRow; … … 6798 7632 ClpSimplex * model = new ClpSimplex(); 6799 7633 // return if nothing 6800 if (!numberConstraints ) {7634 if (!numberConstraints&&linearObjective) { 6801 7635 delete [] which; 6802 7636 model->loadProblem(coinModel); … … 6814 7648 int saveNumber=numberConstraints; 6815 7649 numberConstraints=0; 7650 ClpQuadraticObjective * quadObj = NULL; 6816 7651 if (!linearObjective) { 6817 7652 int numberQuadratic=0; 6818 int largestColumn=-1;6819 7653 CoinZeroN(linearTerm,numberColumns); 6820 7654 for (iColumn=0;iColumn<numberColumns;iColumn++) { … … 6823 7657 const char * expr = coinModel.getColumnObjectiveAsString(iColumn); 6824 7658 if (strcmp(expr,"Numeric")) { 6825 largestColumn = CoinMax(largestColumn,iColumn);6826 7659 // value*x*y 6827 7660 char temp[20000]; … … 6835 7668 if (jColumn>=0) { 6836 7669 columnQuadratic[numberQuadratic]=jColumn; 6837 elementQuadratic[numberQuadratic++]=2.0*value; // convention 6838 largestColumn = CoinMax(largestColumn,jColumn); 7670 if (jColumn!=iColumn) 7671 elementQuadratic[numberQuadratic++]=1.0*value; // convention 7672 else if (jColumn==iColumn) 7673 elementQuadratic[numberQuadratic++]=2.0*value; // convention 6839 7674 } else if (jColumn==-2) { 6840 7675 linearTerm[iColumn] = value; 6841 // and put in as row -16842 columnQuadratic[numberQuadratic]=-1;6843 elementQuadratic[numberQuadratic++]=value;6844 largestColumn = CoinMax(largestColumn,iColumn);6845 7676 } else { 6846 7677 printf("bad nonlinear term %s\n",temp); … … 6852 7683 // linear part 6853 7684 linearTerm[iColumn] = coinModel.getColumnObjective(iColumn); 6854 // and put in as row -16855 columnQuadratic[numberQuadratic]=-1;6856 elementQuadratic[numberQuadratic++]=linearTerm[iColumn];6857 if (linearTerm[iColumn])6858 largestColumn = CoinMax(largestColumn,iColumn);6859 7685 } 6860 7686 } 6861 7687 startQuadratic[numberColumns] = numberQuadratic; 6862 // here we create ClpConstraint 6863 constraints[numberConstraints++] = new ClpConstraintQuadratic(-1, largestColumn+1, numberColumns, 6864 startQuadratic,columnQuadratic,elementQuadratic); 7688 quadObj = new ClpQuadraticObjective(linearTerm, numberColumns, 7689 startQuadratic,columnQuadratic,elementQuadratic); 6865 7690 } 6866 7691 int iConstraint; … … 6949 7774 delete [] which; 6950 7775 model->loadProblem(coinModel2); 6951 int returnCode = model->nonlinearSLP(numberConstraints, constraints, 6952 numberPasses,deltaTolerance); 6953 // See if any integers 6954 for (iConstraint=0;iConstraint<saveNumber;iConstraint++) 6955 delete constraints[iConstraint]; 6956 7776 if (quadObj) 7777 model->setObjective(quadObj); 7778 delete quadObj; 7779 int returnCode; 7780 if (numberConstraints) { 7781 returnCode = model->nonlinearSLP(numberConstraints, constraints, 7782 numberPasses,deltaTolerance); 7783 for (iConstraint=0;iConstraint<saveNumber;iConstraint++) 7784 delete constraints[iConstraint]; 7785 } else { 7786 returnCode = model->nonlinearSLP(numberPasses,deltaTolerance); 7787 } 6957 7788 delete [] constraints; 6958 7789 assert (!returnCode); 6959 7790 return model; 6960 } 7791 #else 7792 printf("loadNonLinear needs ampl\n"); 7793 abort(); 7794 return NULL; 7795 #endif 7796 } 7797 OsiChooseStrongSubset::OsiChooseStrongSubset() : 7798 OsiChooseStrong(), 7799 numberObjectsToUse_(0) 7800 { 7801 } 7802 7803 OsiChooseStrongSubset::OsiChooseStrongSubset(const OsiSolverInterface * solver) : 7804 OsiChooseStrong(solver), 7805 numberObjectsToUse_(-1) 7806 { 7807 } 7808 7809 OsiChooseStrongSubset::OsiChooseStrongSubset(const OsiChooseStrongSubset & rhs) 7810 : OsiChooseStrong(rhs) 7811 { 7812 numberObjectsToUse_ = -1; 7813 } 7814 7815 OsiChooseStrongSubset & 7816 OsiChooseStrongSubset::operator=(const OsiChooseStrongSubset & rhs) 7817 { 7818 if (this != &rhs) { 7819 OsiChooseStrong::operator=(rhs); 7820 numberObjectsToUse_ = -1; 7821 } 7822 return *this; 7823 } 7824 7825 7826 OsiChooseStrongSubset::~OsiChooseStrongSubset () 7827 { 7828 } 7829 7830 // Clone 7831 OsiChooseVariable * 7832 OsiChooseStrongSubset::clone() const 7833 { 7834 return new OsiChooseStrongSubset(*this); 7835 } 7836 // Initialize 7837 int 7838 OsiChooseStrongSubset::setupList ( OsiBranchingInformation *info, bool initialize) 7839 { 7840 assert (solver_==info->solver_); 7841 // Only has to work with Clp 7842 OsiSolverInterface * solverA = const_cast<OsiSolverInterface *> (solver_); 7843 OsiSolverLink * solver = dynamic_cast<OsiSolverLink *> (solverA); 7844 assert (solver); 7845 int numberObjects = solver->numberObjects(); 7846 if (numberObjects>numberObjects_) { 7847 // redo useful arrays 7848 delete [] upTotalChange_; 7849 delete [] downTotalChange_; 7850 delete [] upNumber_; 7851 delete [] downNumber_; 7852 numberObjects_ = solver->numberObjects(); 7853 upTotalChange_ = new double [numberObjects_]; 7854 downTotalChange_ = new double [numberObjects_]; 7855 upNumber_ = new int [numberObjects_]; 7856 downNumber_ = new int [numberObjects_]; 7857 CoinZeroN(upTotalChange_,numberObjects_); 7858 CoinZeroN(downTotalChange_,numberObjects_); 7859 CoinZeroN(upNumber_,numberObjects_); 7860 CoinZeroN(downNumber_,numberObjects_); 7861 } 7862 if (numberObjectsToUse_<0) { 7863 // Sort objects so bilinear at end 7864 OsiObject ** sorted = new OsiObject * [numberObjects]; 7865 OsiObject ** objects = solver->objects(); 7866 numberObjects_=0; 7867 int numberBiLinear=0; 7868 int i; 7869 for (i=0;i<numberObjects;i++) { 7870 OsiObject * obj = objects[i]; 7871 OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj); 7872 if (!objB) 7873 objects[numberObjects_++]=obj; 7874 else 7875 sorted[numberBiLinear++]=obj; 7876 } 7877 numberObjectsToUse_ = numberObjects_; 7878 for (i=0;i<numberBiLinear;i++) 7879 objects[numberObjects_++]=sorted[i]; 7880 delete [] sorted; 7881 // See if any master objects 7882 for (i=0;i<numberObjectsToUse_;i++) { 7883 OsiUsesBiLinear * obj = dynamic_cast<OsiUsesBiLinear *> (objects[i]); 7884 if (obj) 7885 obj->addBiLinearObjects(solver); 7886 } 7887 } 7888 solver->setNumberObjects(numberObjectsToUse_); 7889 numberObjects_=numberObjectsToUse_; 7890 // Use shadow prices 7891 info->defaultDual_=0.0; 7892 int numberUnsatisfied=OsiChooseStrong::setupList ( info, initialize); 7893 solver->setNumberObjects(numberObjects); 7894 numberObjects_=numberObjects; 7895 return numberUnsatisfied; 7896 } 7897 /* Choose a variable 7898 Returns - 7899 -1 Node is infeasible 7900 0 Normal termination - we have a candidate 7901 1 All looks satisfied - no candidate 7902 2 We can change the bound on a variable - but we also have a strong branching candidate 7903 3 We can change the bound on a variable - but we have a non-strong branching candidate 7904 4 We can change the bound on a variable - no other candidates 7905 We can pick up branch from whichObject() and whichWay() 7906 We can pick up a forced branch (can change bound) from whichForcedObject() and whichForcedWay() 7907 If we have a solution then we can pick up from goodObjectiveValue() and goodSolution() 7908 */ 7909 int 7910 OsiChooseStrongSubset::chooseVariable( OsiSolverInterface * solver, OsiBranchingInformation *info, bool fixVariables) 7911 { 7912 int numberObjects = solver->numberObjects(); 7913 solver->setNumberObjects(numberObjectsToUse_); 7914 numberObjects_=numberObjectsToUse_; 7915 // Use shadow prices 7916 info->defaultDual_=0.0; 7917 int returnCode=OsiChooseStrong::chooseVariable(solver,info,fixVariables); 7918 solver->setNumberObjects(numberObjects); 7919 numberObjects_=numberObjects; 7920 return returnCode; 7921 } 7922 /** Default Constructor 7923 7924 Equivalent to an unspecified binary variable. 7925 */ 7926 OsiUsesBiLinear::OsiUsesBiLinear () 7927 : OsiSimpleInteger(), 7928 numberBiLinear_(0), 7929 type_(0), 7930 objects_(NULL) 7931 { 7932 } 7933 7934 /** Useful constructor 7935 7936 Loads actual upper & lower bounds for the specified variable. 7937 */ 7938 OsiUsesBiLinear::OsiUsesBiLinear (const OsiSolverInterface * solver, int iColumn, int type) 7939 : OsiSimpleInteger(solver,iColumn), 7940 numberBiLinear_(0), 7941 type_(type), 7942 objects_(NULL) 7943 { 7944 if (type_) { 7945 assert(originalLower_==floor(originalLower_+0.5)); 7946 assert(originalUpper_==floor(originalUpper_+0.5)); 7947 } 7948 } 7949 7950 7951 // Useful constructor - passed solver index and original bounds 7952 OsiUsesBiLinear::OsiUsesBiLinear ( int iColumn, double lower, double upper, int type) 7953 : OsiSimpleInteger(iColumn,lower,upper), 7954 numberBiLinear_(0), 7955 type_(type), 7956 objects_(NULL) 7957 { 7958 if (type_) { 7959 assert(originalLower_==floor(originalLower_+0.5)); 7960 assert(originalUpper_==floor(originalUpper_+0.5)); 7961 } 7962 } 7963 7964 // Useful constructor - passed simple integer 7965 OsiUsesBiLinear::OsiUsesBiLinear ( const OsiSimpleInteger &rhs, int type) 7966 : OsiSimpleInteger(rhs), 7967 numberBiLinear_(0), 7968 type_(type), 7969 objects_(NULL) 7970 { 7971 if (type_) { 7972 assert(originalLower_==floor(originalLower_+0.5)); 7973 assert(originalUpper_==floor(originalUpper_+0.5)); 7974 } 7975 } 7976 7977 // Copy constructor 7978 OsiUsesBiLinear::OsiUsesBiLinear ( const OsiUsesBiLinear & rhs) 7979 :OsiSimpleInteger(rhs), 7980 numberBiLinear_(0), 7981 type_(rhs.type_), 7982 objects_(NULL) 7983 7984 { 7985 } 7986 7987 // Clone 7988 OsiObject * 7989 OsiUsesBiLinear::clone() const 7990 { 7991 return new OsiUsesBiLinear(*this); 7992 } 7993 7994 // Assignment operator 7995 OsiUsesBiLinear & 7996 OsiUsesBiLinear::operator=( const OsiUsesBiLinear& rhs) 7997 { 7998 if (this!=&rhs) { 7999 OsiSimpleInteger::operator=(rhs); 8000 delete [] objects_; 8001 numberBiLinear_ = 0; 8002 type_ = rhs.type_; 8003 objects_ = NULL; 8004 } 8005 return *this; 8006 } 8007 8008 // Destructor 8009 OsiUsesBiLinear::~OsiUsesBiLinear () 8010 { 8011 delete [] objects_; 8012 } 8013 // Infeasibility - large is 0.5 8014 double 8015 OsiUsesBiLinear::infeasibility(const OsiBranchingInformation * info, int & whichWay) const 8016 { 8017 assert (type_==0); // just continuous for now 8018 double value = info->solution_[columnNumber_]; 8019 value = CoinMax(value, info->lower_[columnNumber_]); 8020 value = CoinMin(value, info->upper_[columnNumber_]); 8021 infeasibility_ = 0.0; 8022 for (int i=0;i<numberBiLinear_;i++) { 8023 OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (objects_[i]); 8024 assert (obj); 8025 obj->getPseudoShadow(info); 8026 infeasibility_ += objects_[i]->infeasibility(); 8027 } 8028 bool satisfied=false; 8029 whichWay=-1; 8030 if (infeasibility_<=info->integerTolerance_) { 8031 otherInfeasibility_ = 1.0; 8032 satisfied=true; 8033 infeasibility_ = 0.0; 8034 } else { 8035 otherInfeasibility_ = 10.0*infeasibility_; 8036 if (value-info->lower_[columnNumber_]> 8037 info->upper_[columnNumber_]-value) 8038 whichWay=1; 8039 else 8040 whichWay=-1; 8041 } 8042 if (preferredWay_>=0&&!satisfied) 8043 whichWay = preferredWay_; 8044 whichWay_=whichWay; 8045 return infeasibility_; 8046 } 8047 // Creates a branching object 8048 OsiBranchingObject * 8049 OsiUsesBiLinear::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const 8050 { 8051 double value = info->solution_[columnNumber_]; 8052 value = CoinMax(value, info->lower_[columnNumber_]); 8053 value = CoinMin(value, info->upper_[columnNumber_]); 8054 assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]); 8055 double nearest = floor(value+0.5); 8056 double integerTolerance = info->integerTolerance_; 8057 if (fabs(value-nearest)<integerTolerance) { 8058 // adjust value 8059 if (nearest!=info->upper_[columnNumber_]) 8060 value = nearest+2.0*integerTolerance; 8061 else 8062 value = nearest-2.0*integerTolerance; 8063 } 8064 OsiBranchingObject * branch = new OsiIntegerBranchingObject(solver,this,way, 8065 value,value,value); 8066 return branch; 8067 } 8068 // This looks at solution and sets bounds to contain solution 8069 /** More precisely: it first forces the variable within the existing 8070 bounds, and then tightens the bounds to fix the variable at the 8071 nearest integer value. 8072 */ 8073 double 8074 OsiUsesBiLinear::feasibleRegion(OsiSolverInterface * solver, 8075 const OsiBranchingInformation * info) const 8076 { 8077 double value = info->solution_[columnNumber_]; 8078 double newValue = CoinMax(value, info->lower_[columnNumber_]); 8079 newValue = CoinMin(newValue, info->upper_[columnNumber_]); 8080 solver->setColLower(columnNumber_,newValue); 8081 solver->setColUpper(columnNumber_,newValue); 8082 return fabs(value-newValue); 8083 } 8084 // Add all bi-linear objects 8085 void 8086 OsiUsesBiLinear::addBiLinearObjects(OsiSolverLink * solver) 8087 { 8088 delete [] objects_; 8089 numberBiLinear_=0; 8090 OsiObject ** objects = solver->objects(); 8091 int i; 8092 int numberObjects = solver->numberObjects(); 8093 for (i=0;i<numberObjects;i++) { 8094 OsiObject * obj = objects[i]; 8095 OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj); 8096 if (objB) { 8097 if (objB->xColumn()==columnNumber_||objB->yColumn()==columnNumber_) 8098 numberBiLinear_++; 8099 } 8100 } 8101 if (numberBiLinear_) { 8102 objects_ = new OsiObject * [numberBiLinear_]; 8103 numberBiLinear_=0; 8104 for (i=0;i<numberObjects;i++) { 8105 OsiObject * obj = objects[i]; 8106 OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj); 8107 if (objB) { 8108 if (objB->xColumn()==columnNumber_||objB->yColumn()==columnNumber_) 8109 objects_[numberBiLinear_++] = obj;; 8110 } 8111 } 8112 } else { 8113 objects_=NULL; 8114 } 8115 } -
trunk/Cbc/src/CbcLinked.hpp
r640 r687 8 8 CglTemporary 9 9 */ 10 #ifndef COIN_HAS_LINK11 #ifdef COIN_HAS_ASL12 #define COIN_HAS_LINK13 #endif14 #endif15 10 #include "CoinModel.hpp" 16 #ifdef COIN_HAS_LINK17 11 #include "OsiClpSolverInterface.hpp" 12 #include "OsiChooseVariable.hpp" 13 #include "CbcFathom.hpp" 18 14 class CbcModel; 19 15 class CoinPackedMatrix; … … 21 17 class OsiObject; 22 18 class CglStored; 23 //#############################################################################24 25 19 /** 26 20 … … 29 23 */ 30 24 31 class OsiSolverLink : public OsiClpSolverInterface{25 class OsiSolverLink : public CbcOsiSolver { 32 26 33 27 public: … … 63 57 heuristic solution 64 58 Returns solution array 59 mode - 60 0 just get continuous 61 1 round and try normal bab 62 2 use defaultBound_ to bound integer variables near current solution 65 63 */ 66 64 double * heuristicSolution(int numberPasses,double deltaTolerance,int mode); … … 112 110 /// Analyze constraints to see which are convex (quadratic) 113 111 void analyzeObjects(); 112 /// Add reformulated bilinear constraints 113 void addTighterConstraints(); 114 114 /// Objective value of best solution found internally 115 115 inline double bestObjectiveValue() const … … 163 163 inline int integerPriority() const 164 164 { return integerPriority_;}; 165 /// Objective transfer variable if one 166 inline int objectiveVariable() const 167 { return objectiveVariable_;} 165 168 /// Set biLinear priority 166 169 inline void setBiLinearPriority(int value) … … 169 172 inline int biLinearPriority() const 170 173 { return biLinearPriority_;}; 171 /// Set Cbc Model172 inline void setCbcModel(CbcModel * model)173 { cbcModel_=model;};174 174 /// Return CoinModel 175 175 inline const CoinModel * coinModel() const … … 177 177 /// Set all biLinear priorities on x-x variables 178 178 void setBiLinearPriorities(int value, double meshSize=1.0); 179 /** Set options and priority on all or some biLinear variables 180 1 - on I-I 181 2 - on I-x 182 4 - on x-x 183 or combinations. 184 -1 means leave (for priority value and strategy value) 185 */ 186 void setBranchingStrategyOnVariables(int strategyValue, int priorityValue=-1, 187 int mode=7); 179 188 /// Set all mesh sizes on x-x variables 180 189 void setMeshSizes(double value); … … 212 221 /// Copy of quadratic model if one 213 222 ClpSimplex * quadraticModel_; 214 /// Pointer back to CbcModel215 CbcModel * cbcModel_;216 223 /// Number of rows with nonLinearities 217 224 int numberNonLinearRows_; … … 236 243 1 bit (2) - quadratic only in objective (add OA cuts) 237 244 2 bit (4) - convex 238 4 bit (8) - try adding OA cuts 245 3 bit (8) - try adding OA cuts 246 4 bit (16) - add linearized constraints 239 247 */ 240 248 int specialOptions2_; … … 681 689 */ 682 690 OsiBiLinear (OsiSolverInterface * solver, int xColumn, 691 int yColumn, int xyRow, double coefficient, 692 double xMesh, double yMesh, 693 int numberExistingObjects=0,const OsiObject ** objects=NULL ); 694 695 /** Useful constructor - 696 This Adds in rows and variables to construct valid Linked Ordered Set 697 Adds extra constraints to match other x/y 698 So note not const model 699 */ 700 OsiBiLinear (CoinModel * coinModel, int xColumn, 683 701 int yColumn, int xyRow, double coefficient, 684 702 double xMesh, double yMesh, … … 793 811 next bit 794 812 8 set to say don't use in feasible region 813 next bit 814 16 set to say - Always satisfied !! 795 815 */ 796 816 inline int branchingStrategy() const … … 820 840 /// Compute lambdas (third entry in each .B is current value) (nonzero if bad) 821 841 double computeLambdas(const double xB[3], const double yB[3],const double xybar[4],double lambda[4]) const; 842 /// Adds in data for extra row with variable coefficients 843 void addExtraRow(int row, double multiplier); 844 /// Sets infeasibility and other when pseudo shadow prices 845 void getPseudoShadow(const OsiBranchingInformation * info); 822 846 823 847 protected: … … 857 881 next bit 858 882 8 set to say don't use in feasible region 883 next bit 884 16 set to say - Always satisfied !! 859 885 */ 860 886 int branchingStrategy_; … … 875 901 /// Convexity row 876 902 int convexity_; 903 /// Number of extra rows (coefficients to be modified) 904 int numberExtraRows_; 905 /// Multiplier for coefficient on row 906 double * multiplier_; 907 /// Row number 908 int * extraRow_; 877 909 /// Which chosen -1 none, 0 x, 1 y 878 910 mutable short chosen_; … … 972 1004 int numberPoints_; 973 1005 }; 974 /// Define a single integer class - but one where you ke p branching until fixed even if satsified1006 /// Define a single integer class - but one where you keep branching until fixed even if satisfied 975 1007 976 1008 … … 1013 1045 /// data 1014 1046 1047 }; 1048 /** Define a single variable class which is involved with OsiBiLinear objects. 1049 This is used so can make better decision on where to branch as it can look at 1050 all objects. 1051 1052 This version sees if it can re-use code from OsiSimpleInteger 1053 even if not an integer variable. If not then need to duplicate code. 1054 */ 1055 1056 1057 class OsiUsesBiLinear : public OsiSimpleInteger { 1058 1059 public: 1060 1061 /// Default Constructor 1062 OsiUsesBiLinear (); 1063 1064 /// Useful constructor - passed solver index 1065 OsiUsesBiLinear (const OsiSolverInterface * solver, int iColumn, int type); 1066 1067 /// Useful constructor - passed solver index and original bounds 1068 OsiUsesBiLinear (int iColumn, double lower, double upper, int type); 1069 1070 /// Useful constructor - passed simple integer 1071 OsiUsesBiLinear (const OsiSimpleInteger & rhs, int type); 1072 1073 /// Copy constructor 1074 OsiUsesBiLinear ( const OsiUsesBiLinear & rhs); 1075 1076 /// Clone 1077 virtual OsiObject * clone() const; 1078 1079 /// Assignment operator 1080 OsiUsesBiLinear & operator=( const OsiUsesBiLinear& rhs); 1081 1082 /// Destructor 1083 virtual ~OsiUsesBiLinear (); 1084 1085 /// Infeasibility - large is 0.5 1086 virtual double infeasibility(const OsiBranchingInformation * info, int & whichWay) const; 1087 /** Creates a branching object 1088 1089 The preferred direction is set by \p way, 0 for down, 1 for up. 1090 */ 1091 virtual OsiBranchingObject * createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const; 1092 /** Set bounds to fix the variable at the current value. 1093 1094 Given an current value, set the lower and upper bounds to fix the 1095 variable. Returns amount it had to move variable. 1096 */ 1097 virtual double feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const; 1098 /// Add all bi-linear objects 1099 void addBiLinearObjects(OsiSolverLink * solver); 1100 protected: 1101 /// data 1102 /// Number of bilinear objects (maybe could be more general) 1103 int numberBiLinear_; 1104 /// Type of variable - 0 continuous, 1 integer 1105 int type_; 1106 /// Objects 1107 OsiObject ** objects_; 1108 }; 1109 /** This class chooses a variable to branch on 1110 1111 This is just as OsiChooseStrong but it fakes it so only 1112 first so many are looked at in this phase 1113 1114 */ 1115 1116 class OsiChooseStrongSubset : public OsiChooseStrong { 1117 1118 public: 1119 1120 /// Default Constructor 1121 OsiChooseStrongSubset (); 1122 1123 /// Constructor from solver (so we can set up arrays etc) 1124 OsiChooseStrongSubset (const OsiSolverInterface * solver); 1125 1126 /// Copy constructor 1127 OsiChooseStrongSubset (const OsiChooseStrongSubset &); 1128 1129 /// Assignment operator 1130 OsiChooseStrongSubset & operator= (const OsiChooseStrongSubset& rhs); 1131 1132 /// Clone 1133 virtual OsiChooseVariable * clone() const; 1134 1135 /// Destructor 1136 virtual ~OsiChooseStrongSubset (); 1137 1138 /** Sets up strong list and clears all if initialize is true. 1139 Returns number of infeasibilities. 1140 If returns -1 then has worked out node is infeasible! 1141 */ 1142 virtual int setupList ( OsiBranchingInformation *info, bool initialize); 1143 /** Choose a variable 1144 Returns - 1145 -1 Node is infeasible 1146 0 Normal termination - we have a candidate 1147 1 All looks satisfied - no candidate 1148 2 We can change the bound on a variable - but we also have a strong branching candidate 1149 3 We can change the bound on a variable - but we have a non-strong branching candidate 1150 4 We can change the bound on a variable - no other candidates 1151 We can pick up branch from bestObjectIndex() and bestWhichWay() 1152 We can pick up a forced branch (can change bound) from firstForcedObjectIndex() and firstForcedWhichWay() 1153 If we have a solution then we can pick up from goodObjectiveValue() and goodSolution() 1154 If fixVariables is true then 2,3,4 are all really same as problem changed 1155 */ 1156 virtual int chooseVariable( OsiSolverInterface * solver, OsiBranchingInformation *info, bool fixVariables); 1157 1158 /// Number of objects to use 1159 inline int numberObjectsToUse() const 1160 { return numberObjectsToUse_;}; 1161 /// Set number of objects to use 1162 inline void setNumberObjectsToUse(int value) 1163 { numberObjectsToUse_ = value;}; 1164 1165 protected: 1166 // Data 1167 /// Number of objects to be used (and set in solver) 1168 int numberObjectsToUse_; 1015 1169 }; 1016 1170 … … 1147 1301 //@} 1148 1302 }; 1149 #endif1150 1303 class ClpSimplex; 1151 1304 /** Return an approximate solution to a CoinModel. -
trunk/Cbc/src/CbcMessage.cpp
r640 r687 31 31 {CBC_STATUS,10,1,"After %d nodes, %d on tree, %g best solution, best possible %g (%.2f seconds)"}, 32 32 {CBC_GAP,11,1,"Exiting as integer gap of %g less than %g or %g%%"}, 33 {CBC_ROUNDING,12,1,"Integer solution of %g found by heuristicafter %d iterations and %d nodes (%.2f seconds)"},33 {CBC_ROUNDING,12,1,"Integer solution of %g found by %s after %d iterations and %d nodes (%.2f seconds)"}, 34 34 {CBC_TREE_SOL,24,1,"Integer solution of %g found by subtree after %d iterations and %d nodes (%.2f seconds)"}, 35 35 {CBC_ROOT,13,1,"At root node, %d cuts changed objective from %g to %g in %d passes"}, … … 49 49 {CBC_START_SUB,28,1,"Starting sub-tree for %s - maximum nodes %d"}, 50 50 {CBC_END_SUB,29,1,"Ending sub-tree for %s"}, 51 {CBC_ HEURISTIC_SOLUTION,30,1,"solution of %g found by %s after %2.f seconds"},51 {CBC_THREAD_STATS,30,1,"%s%? %d used %d times, waiting to start %g, %?%g cpu time,%? %g waiting for threads, %? %d locks, %g locked, %g waiting for locks"}, 52 52 {CBC_CUTS_STATS,31,1,"%d added rows had average density of %g"}, 53 53 {CBC_STRONG_STATS,32,1,"Strong branching done %d times (%d iterations), fathomed %d nodes and fixed %d variables"}, … … 56 56 {CBC_HEURISTICS_OFF,36,1,"Heuristics switched off as %d branching objects are of wrong type"}, 57 57 {CBC_STATUS2,37,1,"%d nodes, %d on tree, best %g - possible %g depth %d unsat %d its %d (%.2f seconds)"}, 58 {CBC_FPUMP1,38,1,"%s"}, 59 {CBC_FPUMP2,39,2,"%s"}, 58 60 {CBC_DUMMY_END,999999,0,""} 59 61 }; -
trunk/Cbc/src/CbcMessage.hpp
r640 r687 54 54 CBC_START_SUB, 55 55 CBC_END_SUB, 56 CBC_ HEURISTIC_SOLUTION,56 CBC_THREAD_STATS, 57 57 CBC_CUTS_STATS, 58 58 CBC_STRONG_STATS, … … 61 61 CBC_HEURISTICS_OFF, 62 62 CBC_STATUS2, 63 CBC_FPUMP1, 64 CBC_FPUMP2, 63 65 CBC_DUMMY_END 64 66 }; -
trunk/Cbc/src/CbcModel.cpp
r650 r687 14 14 //#define NODE_LOG 15 15 //#define GLOBAL_CUTS_JUST_POINTERS 16 #ifndef CLP_FAST_CODE 17 #ifdef NDEBUG 18 #undef NDEBUG 19 #endif 20 #endif 16 21 #include <cassert> 17 22 #include <cmath> … … 36 41 #include "CbcBranchDynamic.hpp" 37 42 #include "CbcHeuristic.hpp" 43 #include "CbcHeuristicFPump.hpp" 38 44 #include "CbcModel.hpp" 39 45 #include "CbcTreeLocal.hpp" … … 48 54 #include "CbcCutGenerator.hpp" 49 55 #include "CbcFeasibilityBase.hpp" 56 #include "CbcFathom.hpp" 50 57 // include Probing 51 58 #include "CglProbing.hpp" … … 60 67 #include "CbcCompareActual.hpp" 61 68 #include "CbcTree.hpp" 69 //#define CBC_THREAD 70 #ifdef CBC_THREAD 71 #include <pthread.h> 72 #ifndef CLP_FAST_CODE 73 //#define CBC_THREAD_DEBUG 1 74 #endif 75 #ifdef CBC_THREAD_DEBUG 76 #ifdef NDEBUG 77 #undef NDEBUG 78 #undef assert 79 # define assert(expression) { \ 80 if (!(expression)) { \ 81 throw CoinError(__STRING(expression), __PRETTY_FUNCTION__, \ 82 "", __FILE__, __LINE__); \ 83 } \ 84 } 85 #endif 86 #endif 87 // To Pass across to doOneNode 88 typedef struct { 89 CbcModel * baseModel; 90 CbcModel * thisModel; 91 CbcNode * node; // filled in every time 92 CbcNode * createdNode; // filled in every time on return 93 pthread_t threadIdOfBase; 94 pthread_mutex_t * mutex; // for locking data 95 pthread_mutex_t * mutex2; // for waking up threads 96 pthread_cond_t * condition2; // for waking up thread 97 int returnCode; // -1 available, 0 busy, 1 finished , 2?? 98 double timeLocked; 99 double timeWaitingToLock; 100 double timeWaitingToStart; 101 double timeInThread; 102 int numberTimesLocked; 103 int numberTimesUnlocked; 104 int numberTimesWaitingToStart; 105 int saveStuff[2]; 106 struct timespec absTime; 107 bool locked; 108 #if CBC_THREAD_DEBUG 109 int threadNumber; 110 #endif 111 } threadStruct; 112 static void * doNodesThread(void * voidInfo); 113 static void * doCutsThread(void * voidInfo); 114 #endif 62 115 /* Various functions local to CbcModel.cpp */ 63 116 … … 449 502 system held by the solver. 450 503 */ 451 452 504 void CbcModel::branchAndBound(int doStatistics) 453 505 … … 461 513 strongInfo_[2]=0; 462 514 numberStrongIterations_ = 0; 515 #ifndef NDEBUG 516 { 517 int i; 518 int n = solver_->getNumCols(); 519 const double *lower = solver_->getColLower() ; 520 const double *upper = solver_->getColUpper() ; 521 for (i=0;i<n;i++) { 522 assert (lower[i]<1.0e10); 523 assert (upper[i]>-1.0e10); 524 } 525 n = solver_->getNumRows(); 526 lower = solver_->getRowLower() ; 527 upper = solver_->getRowUpper() ; 528 for (i=0;i<n;i++) { 529 assert (lower[i]<1.0e10); 530 assert (upper[i]>-1.0e10); 531 } 532 } 533 #endif 463 534 // original solver (only set if pre-processing) 464 535 OsiSolverInterface * originalSolver=NULL; … … 466 537 OsiObject ** originalObject = NULL; 467 538 // Set up strategies 539 #if 0 540 std::string problemName ; 541 solver_->getStrParam(OsiProbName,problemName) ; 542 if (!strcmp(problemName.c_str(),"PP08A")) solver_->activateRowCutDebugger(problemName.c_str()) ; 543 #endif 468 544 if (strategy_) { 469 545 // May do preprocessing … … 533 609 numberObjects_= numberNewIntegers+numberOldIntegers+numberOldOther+nNonInt; 534 610 object_ = new OsiObject * [numberObjects_]; 611 delete [] integerVariable_; 535 612 integerVariable_ = new int [numberNewIntegers+numberOldIntegers]; 536 613 /* … … 943 1020 whichGenerator_ = new int[maximumWhich_] ; 944 1021 memset(whichGenerator_,0,maximumWhich_*sizeof(int)); 945 int currentNumberCuts = 0 ;946 1022 maximumNumberCuts_ = 0 ; 947 1023 currentNumberCuts_ = 0 ; … … 1005 1081 maximumDepthActual_=0; 1006 1082 numberDJFixed_=0.0; 1083 // Do heuristics 1084 { 1085 double * newSolution = new double [numberColumns] ; 1086 double heuristicValue = getCutoff() ; 1087 int found = -1; // no solution found 1088 1089 currentPassNumber_ = 1; // so root heuristics will run 1090 int i; 1091 for (i = 0;i<numberHeuristics_;i++) { 1092 // see if heuristic will do anything 1093 double saveValue = heuristicValue ; 1094 int ifSol = heuristic_[i]->solution(heuristicValue, 1095 newSolution); 1096 if (ifSol>0) { 1097 // better solution found 1098 found = i ; 1099 incrementUsed(newSolution); 1100 // increment number of solutions so other heuristics can test 1101 numberSolutions_++; 1102 numberHeuristicSolutions_++; 1103 } else { 1104 heuristicValue = saveValue ; 1105 } 1106 } 1107 currentPassNumber_ = 0; 1108 /* 1109 Did any of the heuristics turn up a new solution? Record it before we free 1110 the vector. 1111 */ 1112 if (found >= 0) { 1113 // For compiler error on terra cluster! 1114 if (found<numberHeuristics_) 1115 lastHeuristic_ = heuristic_[found]; 1116 else 1117 lastHeuristic_ = heuristic_[0]; 1118 setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ; 1119 CbcTreeLocal * tree 1120 = dynamic_cast<CbcTreeLocal *> (tree_); 1121 if (tree) 1122 tree->passInSolution(bestSolution_,heuristicValue); 1123 } 1124 for (i = 0;i<numberHeuristics_;i++) { 1125 // delete FPump 1126 CbcHeuristicFPump * pump 1127 = dynamic_cast<CbcHeuristicFPump *> (heuristic_[i]); 1128 if (pump) { 1129 delete pump; 1130 numberHeuristics_ --; 1131 for (int j=i;j<numberHeuristics_;j++) 1132 heuristic_[j] = heuristic_[j+1]; 1133 } 1134 } 1135 delete [] newSolution ; 1136 } 1007 1137 statistics_ = NULL; 1008 1138 // Do on switch … … 1037 1167 feasible = solveWithCuts(cuts, 1, 1038 1168 NULL); 1039 #if 01040 currentNumberCuts_ = cuts.sizeRowCuts();1041 if (currentNumberCuts_ >= maximumNumberCuts_) {1042 maximumNumberCuts_ = currentNumberCuts;1043 delete [] addedCuts_;1044 addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];1045 }1046 #endif1047 1169 } 1048 1170 } … … 1064 1186 *dblParam_[CbcAllowableFractionGap]); 1065 1187 if (bestObjective_-bestPossibleObjective_ < testGap && getCutoffIncrement()>=0.0) { 1066 if (bestPossibleObjective_<getCutoff()) {1188 if (bestPossibleObjective_<getCutoff()) 1067 1189 stoppedOnGap_ = true ; 1068 messageHandler()->message(CBC_GAP,messages())1069 << bestObjective_-bestPossibleObjective_1070 << dblParam_[CbcAllowableGap]1071 << dblParam_[CbcAllowableFractionGap]*100.01072 << CoinMessageEol ;1073 secondaryStatus_ = 2;1074 }1075 1190 feasible = false; 1076 1191 } … … 1127 1242 } 1128 1243 } 1129 #if 0 1130 while (anyAction == -1) 1131 { 1132 // Set objective value (not so obvious if NLP etc) 1133 setObjectiveValue(newNode,NULL); 1134 if (numberBeforeTrust_==0 ) { 1135 anyAction = newNode->chooseBranch(this,NULL,numberPassesLeft) ; 1136 } else { 1137 OsiSolverBranch * branches=NULL; 1138 anyAction = newNode->chooseDynamicBranch(this,NULL,branches,numberPassesLeft) ; 1139 if (anyAction==-3) 1140 anyAction = newNode->chooseBranch(this,NULL,numberPassesLeft) ; // dynamic did nothing 1141 } 1142 numberPassesLeft--; 1143 if (anyAction == -1) 1144 { 1145 //if (solverCharacteristics_->solverType()!=4) 1146 feasible = resolve(NULL,10) != 0 ; 1147 //else 1148 //feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,NULL); 1149 if (problemFeasibility_->feasible(this,0)<0) { 1150 feasible=false; // pretend infeasible 1151 } 1152 reducedCostFix() ; 1153 resolved = true ; 1154 # ifdef CBC_DEBUG 1155 printf("Resolve (root) as something fixed, Obj value %g %d rows\n", 1156 solver_->getObjValue(), 1157 solver_->getNumRows()) ; 1158 # endif 1159 if (!feasible) anyAction = -2 ; } 1160 if (anyAction == -2||newNode->objectiveValue() >= cutoff) 1161 { delete newNode ; 1162 newNode = NULL ; 1163 feasible = false ; } } 1164 #else 1165 OsiSolverBranch * branches = NULL; 1166 anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts,resolved, 1167 NULL,NULL,NULL,branches); 1168 if (anyAction == -2||newNode->objectiveValue() >= cutoff) { 1169 if (anyAction != -2) { 1170 // zap parent nodeInfo 1244 OsiSolverBranch * branches = NULL; 1245 anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts,resolved, 1246 NULL,NULL,NULL,branches); 1247 if (anyAction == -2||newNode->objectiveValue() >= cutoff) { 1248 if (anyAction != -2) { 1249 // zap parent nodeInfo 1171 1250 #ifdef COIN_DEVELOP 1172 printf("zapping CbcNodeInfo %x\n",newNode->nodeInfo()->parent()); 1173 #endif 1174 if (newNode->nodeInfo()) 1175 newNode->nodeInfo()->nullParent(); 1176 } 1177 delete newNode ; 1178 newNode = NULL ; 1179 feasible = false ; 1180 } 1181 #endif 1251 printf("zapping CbcNodeInfo %x\n",newNode->nodeInfo()->parent()); 1252 #endif 1253 if (newNode->nodeInfo()) 1254 newNode->nodeInfo()->nullParent(); 1255 } 1256 delete newNode ; 1257 newNode = NULL ; 1258 feasible = false ; 1259 } 1182 1260 } 1183 1261 /* … … 1188 1266 assert (!newNode || newNode->objectiveValue() <= cutoff) ; 1189 1267 // Save address of root node as we don't want to delete it 1190 CbcNode * rootNode = newNode;1191 1268 // initialize for print out 1192 1269 int lastDepth=0; … … 1204 1281 a valid solution for use by setBestSolution(). 1205 1282 */ 1206 CoinWarmStartBasis *lastws = 0;1283 CoinWarmStartBasis *lastws = NULL ; 1207 1284 if (feasible && newNode->branchingObject()) 1208 1285 { if (resolved) … … 1291 1368 */ 1292 1369 numberLongStrong_=0; 1370 double totalTime = 0.0; 1371 #ifdef CBC_THREAD 1372 CbcNode * createdNode=NULL; 1373 CbcModel ** threadModel = NULL; 1374 pthread_t * threadId = NULL; 1375 int * threadCount = NULL; 1376 pthread_mutex_t mutex; 1377 pthread_cond_t condition_main; 1378 pthread_mutex_t condition_mutex; 1379 pthread_mutex_t * mutex2 = NULL; 1380 pthread_cond_t * condition2 = NULL; 1381 threadStruct * threadInfo = NULL; 1382 bool locked=false; 1383 int threadStats[6]; 1384 memset(threadStats,0,sizeof(threadStats)); 1385 double timeWaiting=0.0; 1386 // For now just one model 1387 if (numberThreads_) { 1388 threadId = new pthread_t [numberThreads_]; 1389 threadCount = new int [numberThreads_]; 1390 CoinZeroN(threadCount,numberThreads_); 1391 pthread_mutex_init(&mutex,NULL); 1392 pthread_cond_init(&condition_main,NULL); 1393 pthread_mutex_init(&condition_mutex,NULL); 1394 threadModel = new CbcModel * [numberThreads_+1]; 1395 threadInfo = new threadStruct [numberThreads_+1]; 1396 mutex2 = new pthread_mutex_t [numberThreads_]; 1397 condition2 = new pthread_cond_t [numberThreads_]; 1398 // we don't want a strategy object 1399 CbcStrategy * saveStrategy = strategy_; 1400 strategy_ = NULL; 1401 for (int i=0;i<numberThreads_;i++) { 1402 pthread_mutex_init(mutex2+i,NULL); 1403 pthread_cond_init(condition2+i,NULL); 1404 threadId[i]=0; 1405 threadInfo[i].baseModel=this; 1406 threadModel[i]=new CbcModel(*this); 1407 #ifdef COIN_HAS_CLP 1408 // Solver may need to know about model 1409 CbcModel * thisModel = threadModel[i]; 1410 CbcOsiSolver * solver = 1411 dynamic_cast<CbcOsiSolver *>(thisModel->solver()) ; 1412 if (solver) 1413 solver->setCbcModel(thisModel); 1414 #endif 1415 mutex_ = (void *) (threadInfo+i); 1416 threadModel[i]->moveToModel(this,-1); 1417 threadInfo[i].thisModel=threadModel[i]; 1418 threadInfo[i].node=NULL; 1419 threadInfo[i].createdNode=NULL; 1420 threadInfo[i].threadIdOfBase=pthread_self(); 1421 threadInfo[i].mutex=&mutex; 1422 threadInfo[i].mutex2=mutex2+i; 1423 threadInfo[i].condition2=condition2+i; 1424 threadInfo[i].returnCode=-1; 1425 threadInfo[i].timeLocked=0.0; 1426 threadInfo[i].timeWaitingToLock=0.0; 1427 threadInfo[i].timeWaitingToStart=0.0; 1428 threadInfo[i].timeInThread=0.0; 1429 threadInfo[i].numberTimesLocked=0; 1430 threadInfo[i].numberTimesUnlocked=0; 1431 threadInfo[i].numberTimesWaitingToStart=0; 1432 threadInfo[i].locked=false; 1433 #if CBC_THREAD_DEBUG 1434 threadInfo[i].threadNumber=i+2; 1435 #endif 1436 pthread_create(threadId+i,NULL,doNodesThread,threadInfo+i); 1437 } 1438 strategy_ = saveStrategy; 1439 // Do a partial one for base model 1440 threadInfo[numberThreads_].baseModel=this; 1441 threadModel[numberThreads_]=this; 1442 mutex_ = (void *) (threadInfo+numberThreads_); 1443 threadInfo[numberThreads_].node=NULL; 1444 threadInfo[numberThreads_].mutex=&mutex; 1445 threadInfo[numberThreads_].condition2=&condition_main; 1446 threadInfo[numberThreads_].mutex2=&condition_mutex; 1447 threadInfo[numberThreads_].timeLocked=0.0; 1448 threadInfo[numberThreads_].timeWaitingToLock=0.0; 1449 threadInfo[numberThreads_].numberTimesLocked=0; 1450 threadInfo[numberThreads_].numberTimesUnlocked=0; 1451 threadInfo[numberThreads_].locked=false; 1452 #if CBC_THREAD_DEBUG 1453 threadInfo[numberThreads_].threadNumber=1; 1454 #endif 1455 } 1456 #endif 1293 1457 /* 1294 1458 At last, the actual branch-and-cut search loop, which will iterate until … … 1301 1465 than the current objective cutoff. 1302 1466 */ 1303 while (!tree_->empty()) { 1467 while (true) { 1468 #ifdef CBC_THREAD 1469 if (!locked) { 1470 lockThread(); 1471 locked=true; 1472 } 1473 #endif 1474 if (tree_->empty()) { 1475 #ifdef CBC_THREAD 1476 if (numberThreads_) { 1477 #ifdef COIN_DEVELOP 1478 printf("empty\n"); 1479 #endif 1480 // may still be outstanding nodes 1481 int iThread; 1482 for (iThread=0;iThread<numberThreads_;iThread++) { 1483 if (threadId[iThread]) { 1484 if (threadInfo[iThread].returnCode==0) 1485 break; 1486 } 1487 } 1488 if (iThread<numberThreads_) { 1489 #ifdef COIN_DEVELOP 1490 printf("waiting for thread %d code 0\n",iThread); 1491 #endif 1492 unlockThread(); 1493 locked = false; 1494 pthread_cond_signal(threadInfo[iThread].condition2); // unlock in case 1495 while (true) { 1496 pthread_mutex_lock(&condition_mutex); 1497 struct timespec absTime; 1498 clock_gettime(CLOCK_REALTIME,&absTime); 1499 double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec; 1500 absTime.tv_nsec += 1000000; // millisecond 1501 if (absTime.tv_nsec>=1000000000) { 1502 absTime.tv_nsec -= 1000000000; 1503 absTime.tv_sec++; 1504 } 1505 pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime); 1506 clock_gettime(CLOCK_REALTIME,&absTime); 1507 double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec; 1508 timeWaiting += time2-time; 1509 pthread_mutex_unlock(&condition_mutex); 1510 if (threadInfo[iThread].returnCode!=0) 1511 break; 1512 pthread_cond_signal(threadInfo[iThread].condition2); // unlock 1513 } 1514 threadModel[iThread]->moveToModel(this,1); 1515 assert (threadInfo[iThread].returnCode==1); 1516 // say available 1517 threadInfo[iThread].returnCode=-1; 1518 threadStats[4]++; 1519 #ifdef COIN_DEVELOP 1520 printf("thread %d code now -1\n",iThread); 1521 #endif 1522 continue; 1523 } else { 1524 #ifdef COIN_DEVELOP 1525 printf("no threads at code 0 \n"); 1526 #endif 1527 // now check if any have just finished 1528 for (iThread=0;iThread<numberThreads_;iThread++) { 1529 if (threadId[iThread]) { 1530 if (threadInfo[iThread].returnCode==1) 1531 break; 1532 } 1533 } 1534 if (iThread<numberThreads_) { 1535 unlockThread(); 1536 locked = false; 1537 threadModel[iThread]->moveToModel(this,1); 1538 assert (threadInfo[iThread].returnCode==1); 1539 // say available 1540 threadInfo[iThread].returnCode=-1; 1541 threadStats[4]++; 1542 #ifdef COIN_DEVELOP 1543 printf("thread %d code now -1\n",iThread); 1544 #endif 1545 continue; 1546 } 1547 } 1548 if (!tree_->empty()) { 1549 #ifdef COIN_DEVELOP 1550 printf("tree not empty!!!!!!\n"); 1551 #endif 1552 continue; 1553 } 1554 for (iThread=0;iThread<numberThreads_;iThread++) { 1555 if (threadId[iThread]) { 1556 if (threadInfo[iThread].returnCode!=-1) { 1557 printf("bad end of tree\n"); 1558 abort(); 1559 } 1560 } 1561 } 1562 #ifdef COIN_DEVELOP 1563 printf("finished ************\n"); 1564 #endif 1565 } 1566 unlockThread(); 1567 locked=false; // not needed as break 1568 #endif 1569 break; 1570 } 1571 #ifdef CBC_THREAD 1572 unlockThread(); 1573 locked = false; 1574 #endif 1575 /* 1576 Check for abort on limits: node count, solution count, time, integrality gap. 1577 */ 1578 totalTime = getCurrentSeconds() ; 1579 if (!(numberNodes_ < intParam_[CbcMaxNumNode] && 1580 numberSolutions_ < intParam_[CbcMaxNumSol] && 1581 totalTime < dblParam_[CbcMaximumSeconds] && 1582 !stoppedOnGap_&&!eventHappened_)) { 1583 // out of loop 1584 break; 1585 } 1304 1586 #ifdef BONMIN 1305 1587 assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible()); … … 1308 1590 double newCutoff = getCutoff(); 1309 1591 if (analyzeResults_) { 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1592 // see if we could fix any (more) 1593 int n=0; 1594 double * newLower = analyzeResults_; 1595 double * objLower = newLower+numberIntegers_; 1596 double * newUpper = objLower+numberIntegers_; 1597 double * objUpper = newUpper+numberIntegers_; 1598 for (int i=0;i<numberIntegers_;i++) { 1599 if (objLower[i]>newCutoff) { 1600 n++; 1601 if (objUpper[i]>newCutoff) { 1602 newCutoff = -COIN_DBL_MAX; 1603 break; 1604 } 1605 } else if (objUpper[i]>newCutoff) { 1606 n++; 1607 } 1608 } 1609 if (newCutoff==-COIN_DBL_MAX) { 1610 printf("Root analysis says finished\n"); 1611 } else if (n>numberFixedNow_) { 1612 printf("%d more fixed by analysis - now %d\n",n-numberFixedNow_,n); 1613 numberFixedNow_=n; 1614 } 1333 1615 } 1334 1616 if (eventHandler) { 1335 if (!eventHandler->event(CbcEventHandler::solution)) { 1336 eventHappened_=true; // exit 1337 } 1338 } 1617 if (!eventHandler->event(CbcEventHandler::solution)) { 1618 eventHappened_=true; // exit 1619 } 1620 } 1621 lockThread(); 1339 1622 // Do from deepest 1340 1623 tree_->cleanTree(this, newCutoff,bestPossibleObjective_) ; 1341 1624 nodeCompare_->newSolution(this) ; 1342 1625 nodeCompare_->newSolution(this,continuousObjective_, 1343 1626 continuousInfeasibilities_) ; 1344 1627 tree_->setComparison(*nodeCompare_) ; 1345 if (tree_->empty()) 1346 break; // finished 1628 if (tree_->empty()) { 1629 unlockThread(); 1630 // For threads we need to check further 1631 //break; // finished 1632 continue; 1633 } 1634 unlockThread(); 1347 1635 } 1348 1636 cutoff = getCutoff() ; 1349 1637 /* 1350 Periodic activities: Opportunities to1638 Periodic activities: Opportunities to 1351 1639 + tweak the nodeCompare criteria, 1352 1640 + check if we've closed the integrality gap enough to quit, … … 1354 1642 */ 1355 1643 if ((numberNodes_%1000) == 0) { 1644 lockThread(); 1356 1645 bool redoTree=nodeCompare_->every1000Nodes(this, numberNodes_) ; 1357 1646 #ifdef CHECK_CUT_SIZE … … 1361 1650 if (redoTree) 1362 1651 tree_->setComparison(*nodeCompare_) ; 1652 unlockThread(); 1363 1653 } 1364 1654 if (saveCompare&&!hotstartSolution_) { … … 1368 1658 saveCompare=NULL; 1369 1659 // redo tree 1660 lockThread(); 1370 1661 tree_->setComparison(*nodeCompare_) ; 1662 unlockThread(); 1371 1663 } 1372 1664 if ((numberNodes_%printFrequency_) == 0) { 1665 lockThread(); 1373 1666 int j ; 1374 1667 int nNodes = tree_->size() ; … … 1379 1672 bestPossibleObjective_ = node->objectiveValue() ; 1380 1673 } 1674 unlockThread(); 1381 1675 if (!intParam_[CbcPrinting]) { 1382 1676 messageHandler()->message(CBC_STATUS,messages()) … … 1405 1699 stoppedOnGap_ = true ; 1406 1700 } 1407 1701 1408 1702 # ifdef CHECK_NODE_FULL 1409 1703 verifyTreeNodes(tree_,*this) ; … … 1412 1706 verifyCutCounts(tree_,*this) ; 1413 1707 # endif 1414 1415 1708 /* 1416 1709 Now we come to the meat of the loop. To create the active subproblem, we'll … … 1423 1716 if (!node) 1424 1717 continue; 1718 #ifndef CBC_THREAD 1719 int currentNumberCuts = 0 ; 1425 1720 currentNode_=node; // so can be accessed elsewhere 1426 1721 #ifdef CBC_DEBUG … … 1429 1724 node->guessedObjectiveValue()); 1430 1725 #endif 1726 #if NEW_UPDATE_OBJECT==0 1431 1727 // Save clone in branching decision 1432 1728 if(branchingMethod_) 1433 1729 branchingMethod_->saveBranchingObject(node->modifiableBranchingObject()); 1434 bool nodeOnTree=false; // Node has been popped 1730 #endif 1435 1731 // Say not on optimal path 1436 1732 bool onOptimalPath=false; … … 1460 1756 CbcNodeInfo * nodeInfo = node->nodeInfo() ; 1461 1757 newNode = NULL ; 1758 int branchesLeft=0; 1462 1759 if (!addCuts(node,lastws,numberFixedNow_>numberFixedAtRoot_)) 1463 1760 { int i ; … … 1471 1768 solverCharacteristics_->setBeforeUpper(upperBefore); 1472 1769 } 1473 bool deleteNode ;1474 1770 if (messageHandler()->logLevel()>2) 1475 1771 node->modifiableBranchingObject()->print(); 1476 int retCode;1477 1772 if (!useOsiBranching) 1478 retCode= node->branch(NULL); // old way1773 branchesLeft = node->branch(NULL); // old way 1479 1774 else 1480 retCode = node->branch(solver_); // new way 1481 if (retCode) 1482 { 1775 branchesLeft = node->branch(solver_); // new way 1776 if (branchesLeft) { 1483 1777 // set nodenumber correctly 1484 1778 node->nodeInfo()->setNodeNumber(numberNodes2_); … … 1497 1791 } 1498 1792 numberNodes2_++; 1499 nodeOnTree=true; // back on tree1500 deleteNode = false ;1793 //nodeOnTree=true; // back on tree 1794 //deleteNode = false ; 1501 1795 # ifdef CHECK_NODE 1502 1796 printf("Node %x pushed back on tree - %d left, %d count\n",node, … … 1504 1798 nodeInfo->numberPointingToThis()) ; 1505 1799 # endif 1506 } 1507 else 1508 { deleteNode = true ; 1509 if (!nodeInfo->numberBranchesLeft()) 1510 nodeInfo->allBranchesGone(); // can clean up 1511 } 1512 1800 } else { 1801 //deleteNode = true ; 1802 if (!nodeInfo->numberBranchesLeft()) 1803 nodeInfo->allBranchesGone(); // can clean up 1804 } 1513 1805 if ((specialOptions_&1)!=0) { 1514 1806 /* … … 1524 1816 } 1525 1817 } 1818 1526 1819 /* 1527 1820 Reoptimize, possibly generating cuts and/or using heuristics to find … … 1597 1890 } 1598 1891 /* 1599 Check for abort on limits: node count, solution count, time, integrality gap.1600 */1601 double totalTime = CoinCpuTime()-dblParam_[CbcStartSeconds] ;1602 if (numberNodes_ < intParam_[CbcMaxNumNode] &&1603 numberSolutions_ < intParam_[CbcMaxNumSol] &&1604 totalTime < dblParam_[CbcMaximumSeconds] &&1605 !stoppedOnGap_&&!eventHappened_)1606 {1607 /*1608 1892 Are we still feasible? If so, create a node and do the work to attach a 1609 1893 branching object, reoptimising as needed if chooseBranch() identifies … … 1645 1929 int numberPassesLeft=5; 1646 1930 checkingNode=true; 1647 #if 01648 while (anyAction == -1)1649 {1650 // Set objective value (not so obvious if NLP etc)1651 setObjectiveValue(newNode,node);1652 if (numberBeforeTrust_==0 ) {1653 anyAction = newNode->chooseBranch(this,node,numberPassesLeft) ;1654 } else {1655 OsiSolverBranch * branches=NULL;1656 anyAction = newNode->chooseDynamicBranch(this,node,branches,numberPassesLeft) ;1657 if (anyAction==-3)1658 anyAction = newNode->chooseBranch(this,node,numberPassesLeft) ; // dynamic did nothing1659 }1660 if (solverCharacteristics_ &&1661 solverCharacteristics_->solutionAddsCuts() && // we are in some OA based bab1662 feasible && (newNode->numberUnsatisfied()==0) //solution has become integer feasible during strong branching1663 )1664 {1665 //in the present case we need to check here integer infeasibility if the node is not fathomed we will have to do the loop1666 // again1667 std::cout<<solver_<<std::endl;1668 resolve(solver_);1669 double objval = solver_->getObjValue();1670 setBestSolution(CBC_SOLUTION, objval,1671 solver_->getColSolution()) ;1672 lastHeuristic_ = NULL;1673 int easy=2;1674 if (!solverCharacteristics_->mipFeasible())//did we prove that the node could be pruned?1675 feasible = false;1676 // Reset the bound now1677 solverCharacteristics_->setMipBound(-COIN_DBL_MAX);1678 1679 1680 //numberPassesLeft++;1681 solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;1682 feasible &= resolve(node ? node->nodeInfo() : NULL,11) != 0 ;1683 solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;1684 resolved = true ;1685 if (problemFeasibility_->feasible(this,0)<0) {1686 feasible=false; // pretend infeasible1687 }1688 if(feasible)1689 anyAction = -1;1690 else1691 anyAction = -2;1692 }1693 /*1694 Yep, false positives for sure. And no easy way to distinguish honest1695 infeasibility from `found a solution and tightened objective target.'1696 1697 if (onOptimalPath)1698 assert (anyAction!=-2); // can be useful but gives false positives on strong1699 */1700 numberPassesLeft--;1701 if (numberPassesLeft<=-1) {1702 if (!numberLongStrong_)1703 messageHandler()->message(CBC_WARNING_STRONG,1704 messages()) << CoinMessageEol ;1705 numberLongStrong_++;1706 }1707 if (anyAction == -1)1708 {1709 // can do quick optimality check1710 int easy=2;1711 solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;1712 //if (solverCharacteristics_->solverType()!=4)1713 feasible = resolve(node ? node->nodeInfo() : NULL,11) != 0 ;1714 //else1715 //feasible = solveWithCuts(cuts,maximumCutPasses_,node);1716 solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;1717 resolved = true ;1718 if (problemFeasibility_->feasible(this,0)<0) {1719 feasible=false; // pretend infeasible1720 }1721 if (feasible)1722 {1723 // Set objective value (not so obvious if NLP etc)1724 setObjectiveValue(newNode,node);1725 reducedCostFix() ;1726 if (newNode->objectiveValue() >= getCutoff())1727 anyAction=-2;1728 }1729 else1730 { anyAction = -2 ; } } }1731 if (anyAction >= 0)1732 { if (resolved)1733 { bool needValidSolution = (newNode->branchingObject() == NULL) ;1734 takeOffCuts(cuts,needValidSolution,NULL) ;1735 # ifdef CHECK_CUT_COUNTS1736 { printf("Number of rows after chooseBranch fix (node)"1737 "(active only) %d\n",1738 numberRowsAtContinuous_+numberNewCuts_+1739 numberOldActiveCuts_) ;1740 const CoinWarmStartBasis* debugws =1741 dynamic_cast<const CoinWarmStartBasis*>1742 (solver_->getWarmStart()) ;1743 debugws->print() ;1744 delete debugws ; }1745 # endif1746 }1747 newNode->createInfo(this,node,lastws,lowerBefore,upperBefore,1748 numberOldActiveCuts_,numberNewCuts_) ;1749 if (newNode->numberUnsatisfied()) {1750 maximumDepthActual_ = CoinMax(maximumDepthActual_,newNode->depth());1751 newNode->initializeInfo() ;1752 newNode->nodeInfo()->addCuts(cuts,newNode->numberBranches(),1753 whichGenerator_) ; } } }1754 else {1755 anyAction = -2 ;1756 // Reset bound anyway (no harm if not odd)1757 solverCharacteristics_->setMipBound(-COIN_DBL_MAX);1758 node->nodeInfo()->decrement();1759 }1760 // May have slipped through i.e. anyAction == 0 and objective above cutoff1761 // I think this will screw up cut reference counts if executed.1762 // We executed addCuts just above. (lh)1763 if ( anyAction >=0 ) {1764 assert (newNode);1765 if (newNode->objectiveValue() >= getCutoff()) {1766 anyAction = -2; // say bad after all1767 #ifdef COIN_DEVELOP1768 printf("zapping2 CbcNodeInfo %x\n",newNode->nodeInfo()->parent());1769 #endif1770 // zap parent nodeInfo1771 if (newNode->nodeInfo())1772 newNode->nodeInfo()->nullParent();1773 }1774 }1775 /*1776 If we end up infeasible, we can delete the new node immediately. Since this1777 node won't be needing the cuts we collected, decrement the reference counts.1778 If we are feasible, then we'll be placing this node into the live set, so1779 increment the reference count in the current (parent) nodeInfo.1780 */1781 if (anyAction == -2)1782 { delete newNode ;1783 newNode = NULL ;1784 // say strong doing well1785 if (checkingNode)1786 setSpecialOptions(specialOptions_|8);1787 for (i = 0 ; i < currentNumberCuts_ ; i++)1788 { if (addedCuts_[i])1789 { if (!addedCuts_[i]->decrement(1))1790 delete addedCuts_[i] ; } } }1791 else1792 { nodeInfo->increment() ;1793 if ((numberNodes_%20)==0) {1794 // say strong not doing as well1795 setSpecialOptions(specialOptions_&~8);1796 }1797 #else1798 1931 OsiSolverBranch * branches=NULL; 1799 1932 // point to useful information … … 1823 1956 } 1824 1957 } 1825 #endif1826 1958 } 1827 1959 /* … … 1866 1998 generate=false; // only special cuts 1867 1999 if (generate) { 1868 generator_[i]->generateCuts(theseCuts,true, NULL) ;2000 generator_[i]->generateCuts(theseCuts,true,solver_,NULL) ; 1869 2001 int numberRowCutsAfter = theseCuts.sizeRowCuts() ; 1870 2002 if (numberRowCutsAfter) { … … 1913 2045 heurValue = saveValue ; } } 1914 2046 if (found >= 0) { 2047 lastHeuristic_ = heuristic_[found]; 1915 2048 setBestSolution(CBC_ROUNDING,heurValue,newSolution) ; 1916 lastHeuristic_ = heuristic_[found];1917 2049 } 1918 2050 delete [] newSolution ; … … 1964 2096 set. 1965 2097 */ 1966 if (deleteNode) 1967 delete node ; } 1968 /* 1969 End of the non-abort actions. The next block of code is executed if we've 1970 aborted because we hit one of the limits. Clean up by deleting the live set 1971 and break out of the node processing loop. Note that on an abort, node may 1972 have been pushed back onto the tree for further processing, in which case 1973 it'll be deleted in cleanTree. We need to check. 1974 */ 2098 if (branchesLeft) 2099 { 2100 } 1975 2101 else 1976 2102 { 1977 if (tree_->size()) 1978 tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ; 1979 delete nextRowCut_; 1980 // We need to get rid of node if is has already been popped from tree 1981 if (!nodeOnTree&&!stoppedOnGap_&&node!=rootNode) 1982 delete node; 1983 if (stoppedOnGap_) 1984 { messageHandler()->message(CBC_GAP,messages()) 1985 << bestObjective_-bestPossibleObjective_ 1986 << dblParam_[CbcAllowableGap] 1987 << dblParam_[CbcAllowableFractionGap]*100.0 1988 << CoinMessageEol ; 1989 secondaryStatus_ = 2; 1990 status_ = 0 ; } 1991 else 1992 if (isNodeLimitReached()) 1993 { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ; 1994 secondaryStatus_ = 3; 1995 status_ = 1 ; } 1996 else 1997 if (totalTime >= dblParam_[CbcMaximumSeconds]) 1998 { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ; 1999 secondaryStatus_ = 4; 2000 status_ = 1 ; } 2001 else 2002 if (eventHappened_) 2003 { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ; 2004 secondaryStatus_ = 5; 2005 status_ = 5 ; } 2006 else 2007 { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ; 2008 secondaryStatus_ = 6; 2009 status_ = 1 ; } 2010 break ; } 2103 if (!nodeInfo->numberBranchesLeft()) 2104 nodeInfo->allBranchesGone(); // can clean up 2105 delete node ; } 2106 } else { 2107 // add cuts found to be infeasible (on bound)! 2108 abort(); 2109 delete node; 2110 } 2011 2111 /* 2012 2112 Delete cuts to get back to the original system. … … 2022 2122 { delRows[i] = i+numberRowsAtContinuous_ ; } 2023 2123 solver_->deleteRows(numberToDelete,delRows) ; 2024 delete [] delRows ; } } 2025 /* 2026 This node fathomed when addCuts atttempted to revive it. Toss it. 2027 */ 2028 else 2029 { delete node ; } } 2124 delete [] delRows ; } 2125 #else // end of not CBC_THREAD 2126 if (!numberThreads_) { 2127 doOneNode(this,node,createdNode); 2128 } else { 2129 threadStats[0]++; 2130 //need to think 2131 int iThread; 2132 // Start one off if any available 2133 for (iThread=0;iThread<numberThreads_;iThread++) { 2134 if (threadInfo[iThread].returnCode==-1) 2135 break; 2136 } 2137 if (iThread<numberThreads_) { 2138 threadInfo[iThread].node=node; 2139 assert (threadInfo[iThread].returnCode==-1); 2140 // say in use 2141 threadInfo[iThread].returnCode=0; 2142 threadModel[iThread]->moveToModel(this,0); 2143 pthread_cond_signal(threadInfo[iThread].condition2); // unlock 2144 threadCount[iThread]++; 2145 } 2146 lockThread(); 2147 locked=true; 2148 // see if any finished 2149 for (iThread=0;iThread<numberThreads_;iThread++) { 2150 if (threadInfo[iThread].returnCode>0) 2151 break; 2152 } 2153 unlockThread(); 2154 locked=false; 2155 if (iThread<numberThreads_) { 2156 threadModel[iThread]->moveToModel(this,1); 2157 assert (threadInfo[iThread].returnCode==1); 2158 // say available 2159 threadInfo[iThread].returnCode=-1; 2160 // carry on 2161 threadStats[3]++; 2162 } else { 2163 // Start one off if any available 2164 for (iThread=0;iThread<numberThreads_;iThread++) { 2165 if (threadInfo[iThread].returnCode==-1) 2166 break; 2167 } 2168 if (iThread<numberThreads_) { 2169 lockThread(); 2170 locked=true; 2171 // If any on tree get 2172 if (!tree_->empty()) { 2173 //node = tree_->bestNode(cutoff) ; 2174 //assert (node); 2175 threadStats[1]++; 2176 continue; // ** get another node 2177 } 2178 unlockThread(); 2179 locked=false; 2180 } 2181 // wait (for debug could sleep and use test) 2182 bool finished=false; 2183 while (!finished) { 2184 pthread_mutex_lock(&condition_mutex); 2185 struct timespec absTime; 2186 clock_gettime(CLOCK_REALTIME,&absTime); 2187 double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec; 2188 absTime.tv_nsec += 1000000; // millisecond 2189 if (absTime.tv_nsec>=1000000000) { 2190 absTime.tv_nsec -= 1000000000; 2191 absTime.tv_sec++; 2192 } 2193 pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime); 2194 clock_gettime(CLOCK_REALTIME,&absTime); 2195 double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec; 2196 timeWaiting += time2-time; 2197 pthread_mutex_unlock(&condition_mutex); 2198 for (iThread=0;iThread<numberThreads_;iThread++) { 2199 if (threadInfo[iThread].returnCode>0) { 2200 finished=true; 2201 break; 2202 } else if (threadInfo[iThread].returnCode==0) { 2203 pthread_cond_signal(threadInfo[iThread].condition2); // unlock 2204 } 2205 } 2206 } 2207 assert (iThread<numberThreads_); 2208 threadModel[iThread]->moveToModel(this,1); 2209 node = threadInfo[iThread].node; 2210 threadInfo[iThread].node=NULL; 2211 assert (threadInfo[iThread].returnCode==1); 2212 // say available 2213 threadInfo[iThread].returnCode=-1; 2214 // carry on 2215 threadStats[2]++; 2216 } 2217 } 2218 //lastDepth=node->depth(); 2219 //lastUnsatisfied=node->numberUnsatisfied(); 2220 #endif // end of CBC_THREAD 2221 } 2222 #ifdef CBC_THREAD 2223 if (numberThreads_) { 2224 //printf("stats "); 2225 //for (unsigned int j=0;j<sizeof(threadStats)/sizeof(int);j++) 2226 //printf("%d ",threadStats[j]); 2227 //printf("\n"); 2228 int i; 2229 // Seems to be bug in CoinCpu on Linux - does threads as well despite documentation 2230 double time=0.0; 2231 for (i=0;i<numberThreads_;i++) 2232 time += threadInfo[i].timeInThread; 2233 bool goodTimer = time<(getCurrentSeconds()); 2234 //bool stopped = (!(numberNodes_ < intParam_[CbcMaxNumNode] && 2235 // numberSolutions_ < intParam_[CbcMaxNumSol] && 2236 // totalTime < dblParam_[CbcMaximumSeconds] && 2237 // !stoppedOnGap_&&!eventHappened_)); 2238 for (i=0;i<numberThreads_;i++) { 2239 while (threadInfo[i].returnCode==0) { 2240 pthread_cond_signal(threadInfo[i].condition2); // unlock 2241 pthread_mutex_lock(&condition_mutex); 2242 struct timespec absTime; 2243 clock_gettime(CLOCK_REALTIME,&absTime); 2244 absTime.tv_nsec += 1000000; // millisecond 2245 if (absTime.tv_nsec>=1000000000) { 2246 absTime.tv_nsec -= 1000000000; 2247 absTime.tv_sec++; 2248 } 2249 pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime); 2250 clock_gettime(CLOCK_REALTIME,&absTime); 2251 pthread_mutex_unlock(&condition_mutex); 2252 } 2253 pthread_cond_signal(threadInfo[i].condition2); // unlock 2254 pthread_mutex_lock(&condition_mutex); // not sure necessary but have had one hang on interrupt 2255 threadModel[i]->numberThreads_=0; // say exit 2256 threadInfo[i].returnCode=0; 2257 pthread_mutex_unlock(&condition_mutex); 2258 pthread_cond_signal(threadInfo[i].condition2); // unlock 2259 //if (!stopped) 2260 pthread_join(threadId[i],NULL); 2261 //else 2262 //pthread_kill(threadId[i]); // kill rather than try and synchronize 2263 threadModel[i]->moveToModel(this,2); 2264 pthread_mutex_destroy (threadInfo[i].mutex2); 2265 pthread_cond_destroy (threadInfo[i].condition2); 2266 assert (threadInfo[i].numberTimesLocked==threadInfo[i].numberTimesUnlocked); 2267 handler_->message(CBC_THREAD_STATS,messages_) 2268 <<"Thread"; 2269 handler_->printing(true) 2270 <<i<<threadCount[i]<<threadInfo[i].timeWaitingToStart; 2271 handler_->printing(goodTimer)<<threadInfo[i].timeInThread; 2272 handler_->printing(false)<<0.0; 2273 handler_->printing(true)<<threadInfo[i].numberTimesLocked 2274 <<threadInfo[i].timeLocked<<threadInfo[i].timeWaitingToLock 2275 <<CoinMessageEol; 2276 } 2277 assert (threadInfo[numberThreads_].numberTimesLocked==threadInfo[numberThreads_].numberTimesUnlocked); 2278 handler_->message(CBC_THREAD_STATS,messages_) 2279 <<"Main thread"; 2280 handler_->printing(false)<<0<<0<<0.0; 2281 handler_->printing(false)<<0.0; 2282 handler_->printing(true)<<timeWaiting; 2283 handler_->printing(true)<<threadInfo[numberThreads_].numberTimesLocked 2284 <<threadInfo[numberThreads_].timeLocked<<threadInfo[numberThreads_].timeWaitingToLock 2285 <<CoinMessageEol; 2286 pthread_mutex_destroy (&mutex); 2287 pthread_cond_destroy (&condition_main); 2288 pthread_mutex_destroy (&condition_mutex); 2289 // delete models (here in case some point to others) 2290 for (i=0;i<numberThreads_;i++) { 2291 delete threadModel[i]; 2292 } 2293 delete [] mutex2; 2294 delete [] condition2; 2295 delete [] threadId; 2296 delete [] threadInfo; 2297 delete [] threadModel; 2298 delete [] threadCount; 2299 mutex_=NULL; 2300 // adjust time to allow for children on some systems 2301 dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren(); 2302 } 2303 #endif 2304 /* 2305 End of the non-abort actions. The next block of code is executed if we've 2306 aborted because we hit one of the limits. Clean up by deleting the live set 2307 and break out of the node processing loop. Note that on an abort, node may 2308 have been pushed back onto the tree for further processing, in which case 2309 it'll be deleted in cleanTree. We need to check. 2310 */ 2311 if (!(numberNodes_ < intParam_[CbcMaxNumNode] && 2312 numberSolutions_ < intParam_[CbcMaxNumSol] && 2313 totalTime < dblParam_[CbcMaximumSeconds] && 2314 !stoppedOnGap_&&!eventHappened_)) { 2315 if (tree_->size()) 2316 tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ; 2317 delete nextRowCut_; 2318 if (stoppedOnGap_) 2319 { messageHandler()->message(CBC_GAP,messages()) 2320 << bestObjective_-bestPossibleObjective_ 2321 << dblParam_[CbcAllowableGap] 2322 << dblParam_[CbcAllowableFractionGap]*100.0 2323 << CoinMessageEol ; 2324 secondaryStatus_ = 2; 2325 status_ = 0 ; } 2326 else 2327 if (isNodeLimitReached()) 2328 { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ; 2329 secondaryStatus_ = 3; 2330 status_ = 1 ; } 2331 else 2332 if (totalTime >= dblParam_[CbcMaximumSeconds]) 2333 { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ; 2334 secondaryStatus_ = 4; 2335 status_ = 1 ; } 2336 else 2337 if (eventHappened_) 2338 { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ; 2339 secondaryStatus_ = 5; 2340 status_ = 5 ; } 2341 else 2342 { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ; 2343 secondaryStatus_ = 6; 2344 status_ = 1 ; } 2345 } 2030 2346 /* 2031 2347 That's it, we've exhausted the search tree, or broken out of the loop because … … 2464 2780 subTreeModel_(NULL), 2465 2781 numberStoppedSubTrees_(0), 2782 mutex_(NULL), 2466 2783 presolve_(0), 2467 2784 numberStrong_(5), … … 2511 2828 searchStrategy_(-1), 2512 2829 numberStrongIterations_(0), 2513 resolveAfterTakeOffCuts_(true) 2830 resolveAfterTakeOffCuts_(true), 2831 #if NEW_UPDATE_OBJECT>1 2832 numberUpdateItems_(0), 2833 maximumNumberUpdateItems_(0), 2834 updateItems_(NULL), 2835 #endif 2836 numberThreads_(0), 2837 threadMode_(0) 2514 2838 { 2515 2839 memset(intParam_,0,sizeof(intParam_)); … … 2594 2918 subTreeModel_(NULL), 2595 2919 numberStoppedSubTrees_(0), 2920 mutex_(NULL), 2596 2921 presolve_(0), 2597 2922 numberStrong_(5), … … 2641 2966 searchStrategy_(-1), 2642 2967 numberStrongIterations_(0), 2643 resolveAfterTakeOffCuts_(true) 2968 resolveAfterTakeOffCuts_(true), 2969 #if NEW_UPDATE_OBJECT>1 2970 numberUpdateItems_(0), 2971 maximumNumberUpdateItems_(0), 2972 updateItems_(NULL), 2973 #endif 2974 numberThreads_(0), 2975 threadMode_(0) 2644 2976 { 2645 2977 memset(intParam_,0,sizeof(intParam_)); … … 2780 3112 numberIntegers_++; 2781 3113 } 3114 delete [] integerVariable_; 2782 3115 if (numberIntegers_) { 2783 delete [] integerVariable_;2784 3116 integerVariable_ = new int [numberIntegers_]; 2785 3117 numberIntegers_=0; … … 2819 3151 subTreeModel_(rhs.subTreeModel_), 2820 3152 numberStoppedSubTrees_(rhs.numberStoppedSubTrees_), 3153 mutex_(NULL), 2821 3154 presolve_(rhs.presolve_), 2822 3155 numberStrong_(rhs.numberStrong_), … … 2855 3188 searchStrategy_(rhs.searchStrategy_), 2856 3189 numberStrongIterations_(rhs.numberStrongIterations_), 2857 resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_) 3190 resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_), 3191 #if NEW_UPDATE_OBJECT>1 3192 numberUpdateItems_(rhs.numberUpdateItems_), 3193 maximumNumberUpdateItems_(rhs.maximumNumberUpdateItems_), 3194 updateItems_(NULL), 3195 #endif 3196 numberThreads_(rhs.numberThreads_), 3197 threadMode_(rhs.threadMode_) 2858 3198 { 2859 3199 memcpy(intParam_,rhs.intParam_,sizeof(intParam_)); … … 2907 3247 object_ = new OsiObject * [numberObjects_]; 2908 3248 int i; 2909 for (i=0;i<numberObjects_;i++) 3249 for (i=0;i<numberObjects_;i++) { 2910 3250 object_[i]=(rhs.object_[i])->clone(); 3251 CbcObject * obj = dynamic_cast <CbcObject *>(object_[i]) ; 3252 // Could be OsiObjects 3253 if (obj) 3254 obj->setModel(this); 3255 } 2911 3256 } else { 2912 3257 object_=NULL; … … 2932 3277 originalColumns_=NULL; 2933 3278 } 3279 #if NEW_UPDATE_OBJECT>1 3280 if (maximumNumberUpdateItems_) { 3281 updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_]; 3282 for (int i=0;i<maximumNumberUpdateItems_;i++) 3283 updateItems_[i] = rhs.updateItems_[i]; 3284 } 3285 #endif 2934 3286 if (maximumWhich_&&rhs.whichGenerator_) 2935 3287 whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_); … … 3105 3457 subTreeModel_ = rhs.subTreeModel_; 3106 3458 numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_; 3459 mutex_ = NULL; 3107 3460 presolve_ = rhs.presolve_; 3108 3461 numberStrong_ = rhs.numberStrong_; … … 3156 3509 numberNewCuts_ = rhs.numberNewCuts_; 3157 3510 resolveAfterTakeOffCuts_=rhs.resolveAfterTakeOffCuts_; 3511 #if NEW_UPDATE_OBJECT>1 3512 numberUpdateItems_ = rhs.numberUpdateItems_; 3513 maximumNumberUpdateItems_ = rhs.maximumNumberUpdateItems_; 3514 delete [] updateItems_; 3515 if (maximumNumberUpdateItems_) { 3516 updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_]; 3517 for (i=0;i<maximumNumberUpdateItems_;i++) 3518 updateItems_[i] = rhs.updateItems_[i]; 3519 } else { 3520 updateItems_ = NULL; 3521 } 3522 #endif 3523 numberThreads_ = rhs.numberThreads_; 3524 threadMode_ = rhs.threadMode_; 3158 3525 sizeMiniTree_ = rhs.sizeMiniTree_; 3159 3526 searchStrategy_ = rhs.searchStrategy_; … … 3329 3696 originalColumns_=NULL; 3330 3697 delete strategy_; 3698 #if NEW_UPDATE_OBJECT>1 3699 delete [] updateItems_; 3700 updateItems_=NULL; 3701 numberUpdateItems_=0; 3702 maximumNumberUpdateItems_=0; 3703 #endif 3331 3704 gutsOfDestructor2(); 3332 3705 } … … 3686 4059 */ 3687 4060 currentNumberCuts_=currentNumberCuts; 3688 if (currentNumberCuts > =maximumNumberCuts_) {4061 if (currentNumberCuts > maximumNumberCuts_) { 3689 4062 maximumNumberCuts_ = currentNumberCuts; 3690 4063 delete [] addedCuts_; … … 3784 4157 reconstructed basis in the solver. 3785 4158 */ 3786 if (node->objectiveValue() < cutoff )4159 if (node->objectiveValue() < cutoff||numberThreads_) 3787 4160 { 3788 4161 # ifdef CBC_CHECK_BASIS … … 3864 4237 set. Adjust the cut reference counts to reflect that we no longer need to 3865 4238 explore the remaining branch arms, hence they will no longer reference any 3866 cuts. Cuts whose reference count falls to zero are deleted. 4239 cuts. Cuts whose reference count falls to zero are deleted. 3867 4240 */ 3868 4241 else 3869 4242 { int i; 3870 int numberLeft = nodeInfo->numberBranchesLeft(); 3871 for (i = 0 ; i < currentNumberCuts ; i++) 3872 { if (addedCuts_[i]) 3873 { if (!addedCuts_[i]->decrement(numberLeft)) 3874 { delete addedCuts_[i]; 3875 addedCuts_[i] = NULL; } } } 4243 if (currentNumberCuts) { 4244 lockThread(); 4245 int numberLeft = nodeInfo->numberBranchesLeft(); 4246 for (i = 0 ; i < currentNumberCuts ; i++) 4247 { if (addedCuts_[i]) 4248 { if (!addedCuts_[i]->decrement(numberLeft)) 4249 { delete addedCuts_[i]; 4250 addedCuts_[i] = NULL; } } } 4251 unlockThread(); 4252 } 3876 4253 return 1 ; } 3877 4254 } … … 3988 4365 saveClpOptions = clpSolver->specialOptions(); 3989 4366 # endif 4367 #ifdef CBC_THREAD 4368 CbcModel ** threadModel = NULL; 4369 pthread_t * threadId = NULL; 4370 pthread_cond_t condition_main; 4371 pthread_mutex_t condition_mutex; 4372 pthread_mutex_t * mutex2 = NULL; 4373 pthread_cond_t * condition2 = NULL; 4374 threadStruct * threadInfo = NULL; 4375 void * saveMutex = NULL; 4376 if (numberThreads_&&(threadMode_&2)!=0&&!numberNodes_) { 4377 threadId = new pthread_t [numberThreads_]; 4378 pthread_cond_init(&condition_main,NULL); 4379 pthread_mutex_init(&condition_mutex,NULL); 4380 threadModel = new CbcModel * [numberThreads_]; 4381 threadInfo = new threadStruct [numberThreads_+1]; 4382 mutex2 = new pthread_mutex_t [numberThreads_]; 4383 condition2 = new pthread_cond_t [numberThreads_]; 4384 saveMutex = mutex_; 4385 for (int i=0;i<numberThreads_;i++) { 4386 pthread_mutex_init(mutex2+i,NULL); 4387 pthread_cond_init(condition2+i,NULL); 4388 threadId[i]=0; 4389 threadModel[i]=new CbcModel; 4390 threadModel[i]->generator_ = new CbcCutGenerator * [1]; 4391 delete threadModel[i]->solver_; 4392 threadModel[i]->solver_=NULL; 4393 threadModel[i]->numberThreads_=numberThreads_; 4394 mutex_ = (void *) (threadInfo+i); 4395 threadInfo[i].thisModel=(CbcModel *) threadModel[i]; 4396 threadInfo[i].baseModel=this; 4397 threadInfo[i].threadIdOfBase=pthread_self(); 4398 threadInfo[i].mutex2=mutex2+i; 4399 threadInfo[i].condition2=condition2+i; 4400 threadInfo[i].returnCode=-1; 4401 pthread_create(threadId+i,NULL,doCutsThread,threadInfo+i); 4402 } 4403 // Do a partial one for base model 4404 threadInfo[numberThreads_].baseModel=this; 4405 mutex_ = (void *) (threadInfo+numberThreads_); 4406 threadInfo[numberThreads_].condition2=&condition_main; 4407 threadInfo[numberThreads_].mutex2=&condition_mutex; 4408 } 4409 #endif 3990 4410 bool feasible = true ; 3991 4411 int lastNumberCuts = 0 ; … … 4026 4446 //printf("zero iterations on first solve of branch\n"); 4027 4447 #endif 4028 if (node&& !node->nodeInfo()->numberBranchesLeft())4448 if (node&&node->nodeInfo()&&!node->nodeInfo()->numberBranchesLeft()) 4029 4449 node->nodeInfo()->allBranchesGone(); // can clean up 4030 4450 feasible = returnCode != 0 ; … … 4035 4455 } 4036 4456 4457 #if NEW_UPDATE_OBJECT==0 4037 4458 // Update branching information if wanted 4038 4459 if(node &&branchingMethod_) 4039 4460 branchingMethod_->updateInformation(solver_,node); 4461 #elif NEW_UPDATE_OBJECT<2 4462 // Update branching information if wanted 4463 if(node &&branchingMethod_) { 4464 OsiBranchingObject * bobj = node->modifiableBranchingObject(); 4465 CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (bobj); 4466 if (cbcobj) { 4467 CbcObject * object = cbcobj->object(); 4468 CbcObjectUpdateData update = object->createUpdateInformation(solver_,node,cbcobj); 4469 object->updateInformation(update); 4470 } else { 4471 branchingMethod_->updateInformation(solver_,node); 4472 } 4473 } 4474 #else 4475 // Update branching information if wanted 4476 if(node &&branchingMethod_) { 4477 OsiBranchingObject * bobj = node->modifiableBranchingObject(); 4478 CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (bobj); 4479 if (cbcobj) { 4480 CbcObject * object = cbcobj->object(); 4481 CbcObjectUpdateData update = object->createUpdateInformation(solver_,node,cbcobj); 4482 // have to compute object number as not saved 4483 CbcSimpleInteger * simpleObject = 4484 dynamic_cast <CbcSimpleInteger *>(object) ; 4485 int iObject; 4486 int iColumn = simpleObject->columnNumber(); 4487 for (iObject = 0 ; iObject < numberObjects_ ; iObject++) { 4488 simpleObject = 4489 dynamic_cast <CbcSimpleInteger *>(object_[iObject]) ; 4490 if (simpleObject->columnNumber()==iColumn) 4491 break; 4492 } 4493 assert (iObject<numberObjects_); 4494 update.objectNumber_ = iObject; 4495 addUpdateInformation(update); 4496 } else { 4497 branchingMethod_->updateInformation(solver_,node); 4498 } 4499 } 4500 #endif 4040 4501 4041 4502 #ifdef CBC_DEBUG … … 4078 4539 sumChangeObjective1_ += solver_->getObjValue()*solver_->getObjSense() 4079 4540 - objectiveValue ; 4080 if ( CoinCpuTime()-dblParam_[CbcStartSeconds]> dblParam_[CbcMaximumSeconds] )4541 if ( getCurrentSeconds() > dblParam_[CbcMaximumSeconds] ) 4081 4542 numberTries=0; // exit 4082 4543 //if ((numberNodes_%100)==0) … … 4244 4705 whichGenerator_[numberViolated++]=-2; 4245 4706 } 4246 double * newSolution = new double [numberColumns] ;4247 double heuristicValue = getCutoff() ;4248 int found = -1; // no solution found4249 4707 4250 4708 // reset probing info 4251 4709 if (probingInfo_) 4252 4710 probingInfo_->initializeFixing(); 4253 for (int i = 0;i<numberCutGenerators_+numberHeuristics_;i++) { 4254 int numberRowCutsBefore = theseCuts.sizeRowCuts() ; 4255 int numberColumnCutsBefore = theseCuts.sizeColCuts() ; 4256 if (i<numberCutGenerators_) { 4257 bool generate = generator_[i]->normal(); 4258 // skip if not optimal and should be (maybe a cut generator has fixed variables) 4259 if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable()) 4260 generate=false; 4261 if (generator_[i]->switchedOff()) 4262 generate=false;; 4711 int i; 4712 if ((threadMode_&2)==0||numberNodes_) { 4713 for (i = 0;i<numberCutGenerators_;i++) { 4714 int numberRowCutsBefore = theseCuts.sizeRowCuts() ; 4715 int numberColumnCutsBefore = theseCuts.sizeColCuts() ; 4716 bool generate = generator_[i]->normal(); 4717 // skip if not optimal and should be (maybe a cut generator has fixed variables) 4718 if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable()) 4719 generate=false; 4720 if (generator_[i]->switchedOff()) 4721 generate=false;; 4263 4722 if (generate) { 4264 4723 bool mustResolve = 4265 generator_[i]->generateCuts(theseCuts,fullScan, node) ;4724 generator_[i]->generateCuts(theseCuts,fullScan,solver_,node) ; 4266 4725 int numberRowCutsAfter = theseCuts.sizeRowCuts() ; 4267 4268 4269 4726 if(numberRowCutsBefore < numberRowCutsAfter && 4727 generator_[i]->mustCallAgain()) 4728 keepGoing=true; // say must go round 4270 4729 // Check last cut to see if infeasible 4271 4730 if(numberRowCutsBefore < numberRowCutsAfter) { 4272 4731 const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ; 4273 4732 if (thisCut->lb()>thisCut->ub()) { … … 4277 4736 } 4278 4737 #ifdef CBC_DEBUG 4279 { 4280 int numberRowCutsAfter = theseCuts.sizeRowCuts() ; 4281 int k ; 4282 for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 4283 OsiRowCut thisCut = theseCuts.rowCut(k) ; 4284 /* check size of elements. 4285 We can allow smaller but this helps debug generators as it 4286 is unsafe to have small elements */ 4287 int n=thisCut.row().getNumElements(); 4288 const int * column = thisCut.row().getIndices(); 4289 const double * element = thisCut.row().getElements(); 4290 //assert (n); 4291 for (int i=0;i<n;i++) { 4292 double value = element[i]; 4293 assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20); 4294 } 4295 } 4296 } 4738 { 4739 int k ; 4740 for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 4741 OsiRowCut thisCut = theseCuts.rowCut(k) ; 4742 /* check size of elements. 4743 We can allow smaller but this helps debug generators as it 4744 is unsafe to have small elements */ 4745 int n=thisCut.row().getNumElements(); 4746 const int * column = thisCut.row().getIndices(); 4747 const double * element = thisCut.row().getElements(); 4748 //assert (n); 4749 for (int i=0;i<n;i++) { 4750 double value = element[i]; 4751 assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20); 4752 } 4753 } 4754 } 4297 4755 #endif 4298 4756 if (mustResolve) { 4299 4300 4301 4302 4303 4304 4305 4306 4307 4308 4309 4310 4311 4757 int returncode = resolve(node ? node->nodeInfo() : NULL,2); 4758 feasible = returnCode != 0 ; 4759 if (returncode<0) 4760 numberTries=0; 4761 if ((specialOptions_&1)!=0) { 4762 debugger = solver_->getRowCutDebugger() ; 4763 if (debugger) 4764 onOptimalPath = (debugger->onOptimalPath(*solver_)) ; 4765 else 4766 onOptimalPath=false; 4767 if (onOptimalPath && !solver_->isDualObjectiveLimitReached()) 4768 assert(feasible) ; 4769 } 4312 4770 if (!feasible) 4313 4771 break ; 4314 4772 } 4315 4773 } 4774 int numberRowCutsAfter = theseCuts.sizeRowCuts() ; 4775 int numberColumnCutsAfter = theseCuts.sizeColCuts() ; 4776 4777 if ((specialOptions_&1)!=0) { 4778 if (onOptimalPath) { 4779 int k ; 4780 for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 4781 OsiRowCut thisCut = theseCuts.rowCut(k) ; 4782 if(debugger->invalidCut(thisCut)) { 4783 solver_->writeMps("badCut"); 4784 #ifdef NDEBUG 4785 printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n", 4786 i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore); 4787 abort(); 4788 #endif 4789 } 4790 assert(!debugger->invalidCut(thisCut)) ; 4791 } 4792 } 4793 } 4794 /* 4795 The cut generator has done its thing, and maybe it generated some 4796 cuts. Do a bit of bookkeeping: load 4797 whichGenerator[i] with the index of the generator responsible for a cut, 4798 and place cuts flagged as global in the global cut pool for the model. 4799 4800 lastNumberCuts is the sum of cuts added in previous iterations; it's the 4801 offset to the proper starting position in whichGenerator. 4802 */ 4803 int numberBefore = 4804 numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ; 4805 int numberAfter = 4806 numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ; 4807 // possibly extend whichGenerator 4808 resizeWhichGenerator(numberBefore, numberAfter); 4809 int j ; 4810 if (fullScan) { 4811 // counts 4812 countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ; 4813 } 4814 countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ; 4815 4816 bool dodgyCuts=false; 4817 for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) { 4818 const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ; 4819 if (thisCut->lb()>1.0e10||thisCut->ub()<-1.0e10) { 4820 dodgyCuts=true; 4821 break; 4822 } 4823 whichGenerator_[numberBefore++] = i ; 4824 if (thisCut->lb()>thisCut->ub()) 4825 violated=-2; // sub-problem is infeasible 4826 if (thisCut->globallyValid()) { 4827 // add to global list 4828 OsiRowCut newCut(*thisCut); 4829 newCut.setGloballyValid(2); 4830 newCut.mutableRow().setTestForDuplicateIndex(false); 4831 globalCuts_.insert(newCut) ; 4832 } 4833 } 4834 if (dodgyCuts) { 4835 for (int k=numberRowCutsAfter-1;k>=j;k--) { 4836 const OsiRowCut * thisCut = theseCuts.rowCutPtr(k) ; 4837 if (thisCut->lb()>thisCut->ub()) 4838 violated=-2; // sub-problem is infeasible 4839 if (thisCut->lb()>1.0e10||thisCut->ub()<-1.0e10) 4840 theseCuts.eraseRowCut(k); 4841 } 4842 numberRowCutsAfter = theseCuts.sizeRowCuts() ; 4843 for (;j<numberRowCutsAfter;j++) { 4844 const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ; 4845 whichGenerator_[numberBefore++] = i ; 4846 if (thisCut->globallyValid()) { 4847 // add to global list 4848 OsiRowCut newCut(*thisCut); 4849 newCut.setGloballyValid(2); 4850 newCut.mutableRow().setTestForDuplicateIndex(false); 4851 globalCuts_.insert(newCut) ; 4852 } 4853 } 4854 } 4855 for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) { 4856 whichGenerator_[numberBefore++] = i ; 4857 const OsiColCut * thisCut = theseCuts.colCutPtr(j) ; 4858 if (thisCut->globallyValid()) { 4859 // add to global list 4860 OsiColCut newCut(*thisCut); 4861 newCut.setGloballyValid(2); 4862 globalCuts_.insert(newCut) ; 4863 } 4864 } 4865 } 4866 // Add in any violated saved cuts 4867 if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) { 4868 int numberOld = theseCuts.sizeRowCuts(); 4869 int numberCuts = slackCuts.sizeRowCuts() ; 4870 int i; 4871 // possibly extend whichGenerator 4872 resizeWhichGenerator(numberOld, numberOld+numberCuts); 4873 for ( i = 0;i<numberCuts;i++) { 4874 const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ; 4875 if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) { 4876 if (messageHandler()->logLevel()>2) 4877 printf("Old cut added - violation %g\n", 4878 thisCut->violated(cbcColSolution_)) ; 4879 whichGenerator_[numberOld++]=-1; 4880 theseCuts.insert(*thisCut) ; 4881 } 4882 } 4883 } 4884 } else { 4885 // do cuts independently 4886 OsiCuts * eachCuts = new OsiCuts [numberCutGenerators_];; 4887 #ifdef CBC_THREAD 4888 if (!threadModel) { 4889 #endif 4890 // generate cuts 4891 for (i = 0;i<numberCutGenerators_;i++) { 4892 bool generate = generator_[i]->normal(); 4893 // skip if not optimal and should be (maybe a cut generator has fixed variables) 4894 if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable()) 4895 generate=false; 4896 if (generator_[i]->switchedOff()) 4897 generate=false;; 4898 if (generate) 4899 generator_[i]->generateCuts(eachCuts[i],fullScan,solver_,node) ; 4900 } 4901 #ifdef CBC_THREAD 4316 4902 } else { 4903 for (i=0;i<numberThreads_;i++) { 4904 // set solver here after cloning 4905 threadModel[i]->solver_=solver_->clone(); 4906 threadModel[i]->numberNodes_ = (fullScan) ? 1 : 0; 4907 } 4908 // generate cuts 4909 for (i = 0;i<numberCutGenerators_;i++) { 4910 bool generate = generator_[i]->normal(); 4911 // skip if not optimal and should be (maybe a cut generator has fixed variables) 4912 if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable()) 4913 generate=false; 4914 if (generator_[i]->switchedOff()) 4915 generate=false;; 4916 if (generate) { 4917 bool finished=false; 4918 int iThread=-1; 4919 // see if any available 4920 for (iThread=0;iThread<numberThreads_;iThread++) { 4921 if (threadInfo[iThread].returnCode) { 4922 finished=true; 4923 break; 4924 } else if (threadInfo[iThread].returnCode==0) { 4925 pthread_cond_signal(threadInfo[iThread].condition2); // unlock 4926 } 4927 } 4928 while (!finished) { 4929 pthread_mutex_lock(&condition_mutex); 4930 struct timespec absTime; 4931 clock_gettime(CLOCK_REALTIME,&absTime); 4932 absTime.tv_nsec += 1000000; // millisecond 4933 if (absTime.tv_nsec>=1000000000) { 4934 absTime.tv_nsec -= 1000000000; 4935 absTime.tv_sec++; 4936 } 4937 pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime); 4938 pthread_mutex_unlock(&condition_mutex); 4939 for (iThread=0;iThread<numberThreads_;iThread++) { 4940 if (threadInfo[iThread].returnCode>0) { 4941 finished=true; 4942 break; 4943 } else if (threadInfo[iThread].returnCode==0) { 4944 pthread_cond_signal(threadInfo[iThread].condition2); // unlock 4945 } 4946 } 4947 } 4948 assert (iThread<numberThreads_); 4949 assert (threadInfo[iThread].returnCode); 4950 threadModel[iThread]->generator_[0]=generator_[i]; 4951 threadModel[iThread]->object_ = (OsiObject **) (eachCuts+i); 4952 // allow to start 4953 threadInfo[iThread].returnCode=0; 4954 pthread_cond_signal(threadInfo[iThread].condition2); // unlock 4955 } 4956 } 4957 // wait 4958 for (int iThread=0;iThread<numberThreads_;iThread++) { 4959 if (threadInfo[iThread].returnCode==0) { 4960 bool finished=false; 4961 pthread_cond_signal(threadInfo[iThread].condition2); // unlock 4962 while (!finished) { 4963 pthread_mutex_lock(&condition_mutex); 4964 struct timespec absTime; 4965 clock_gettime(CLOCK_REALTIME,&absTime); 4966 absTime.tv_nsec += 1000000; // millisecond 4967 if (absTime.tv_nsec>=1000000000) { 4968 absTime.tv_nsec -= 1000000000; 4969 absTime.tv_sec++; 4970 } 4971 pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime); 4972 pthread_mutex_unlock(&condition_mutex); 4973 if (threadInfo[iThread].returnCode>0) { 4974 finished=true; 4975 break; 4976 } else if (threadInfo[iThread].returnCode==0) { 4977 pthread_cond_signal(threadInfo[iThread].condition2); // unlock 4978 } 4979 } 4980 } 4981 assert (threadInfo[iThread].returnCode); 4982 // say available 4983 threadInfo[iThread].returnCode=-1; 4984 delete threadModel[iThread]->solver_; 4985 threadModel[iThread]->solver_=NULL; 4986 } 4987 } 4988 #endif 4989 // Now put together 4990 for (i = 0;i<numberCutGenerators_;i++) { 4991 // add column cuts 4992 int numberColumnCutsBefore = theseCuts.sizeColCuts() ; 4993 int numberColumnCuts = eachCuts[i].sizeColCuts(); 4994 int numberColumnCutsAfter = numberColumnCutsBefore 4995 + numberColumnCuts; 4996 int j; 4997 for (j=0;j<numberColumnCuts;j++) { 4998 theseCuts.insert(eachCuts[i].colCut(j)); 4999 } 5000 int numberRowCutsBefore = theseCuts.sizeRowCuts() ; 5001 int numberRowCuts = eachCuts[i].sizeRowCuts(); 5002 int numberRowCutsAfter = numberRowCutsBefore 5003 + numberRowCuts; 5004 if (numberRowCuts) { 5005 for (j=0;j<numberRowCuts;j++) { 5006 const OsiRowCut * thisCut = eachCuts[i].rowCutPtr(j) ; 5007 if (thisCut->lb()<=1.0e10&&thisCut->ub()>=-1.0e10) 5008 theseCuts.insert(eachCuts[i].rowCut(j)); 5009 } 5010 if (generator_[i]->mustCallAgain()) 5011 keepGoing=true; // say must go round 5012 // Check last cut to see if infeasible 5013 const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ; 5014 if (thisCut->lb()>thisCut->ub()) { 5015 feasible = false; // sub-problem is infeasible 5016 break; 5017 } 5018 } 5019 #ifdef CBC_DEBUG 5020 { 5021 int k ; 5022 for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 5023 OsiRowCut thisCut = theseCuts.rowCut(k) ; 5024 /* check size of elements. 5025 We can allow smaller but this helps debug generators as it 5026 is unsafe to have small elements */ 5027 int n=thisCut.row().getNumElements(); 5028 const int * column = thisCut.row().getIndices(); 5029 const double * element = thisCut.row().getElements(); 5030 //assert (n); 5031 for (int i=0;i<n;i++) { 5032 double value = element[i]; 5033 assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20); 5034 } 5035 } 5036 } 5037 #endif 5038 if ((specialOptions_&1)!=0) { 5039 if (onOptimalPath) { 5040 int k ; 5041 for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 5042 OsiRowCut thisCut = theseCuts.rowCut(k) ; 5043 if(debugger->invalidCut(thisCut)) { 5044 solver_->writeMps("badCut"); 5045 #ifdef NDEBUG 5046 printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n", 5047 i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore); 5048 abort(); 5049 #endif 5050 } 5051 assert(!debugger->invalidCut(thisCut)) ; 5052 } 5053 } 5054 } 5055 /* 5056 The cut generator has done its thing, and maybe it generated some 5057 cuts. Do a bit of bookkeeping: load 5058 whichGenerator[i] with the index of the generator responsible for a cut, 5059 and place cuts flagged as global in the global cut pool for the model. 5060 5061 lastNumberCuts is the sum of cuts added in previous iterations; it's the 5062 offset to the proper starting position in whichGenerator. 5063 */ 5064 int numberBefore = 5065 numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ; 5066 int numberAfter = 5067 numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ; 5068 // possibly extend whichGenerator 5069 resizeWhichGenerator(numberBefore, numberAfter); 5070 if (fullScan) { 5071 // counts 5072 countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ; 5073 } 5074 countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ; 5075 5076 for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) { 5077 whichGenerator_[numberBefore++] = i ; 5078 const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ; 5079 if (thisCut->lb()>thisCut->ub()) 5080 violated=-2; // sub-problem is infeasible 5081 if (thisCut->globallyValid()) { 5082 // add to global list 5083 OsiRowCut newCut(*thisCut); 5084 newCut.setGloballyValid(2); 5085 newCut.mutableRow().setTestForDuplicateIndex(false); 5086 globalCuts_.insert(newCut) ; 5087 } 5088 } 5089 for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) { 5090 whichGenerator_[numberBefore++] = i ; 5091 const OsiColCut * thisCut = theseCuts.colCutPtr(j) ; 5092 if (thisCut->globallyValid()) { 5093 // add to global list 5094 OsiColCut newCut(*thisCut); 5095 newCut.setGloballyValid(2); 5096 globalCuts_.insert(newCut) ; 5097 } 5098 } 5099 } 5100 // Add in any violated saved cuts 5101 if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) { 5102 int numberOld = theseCuts.sizeRowCuts(); 5103 int numberCuts = slackCuts.sizeRowCuts() ; 5104 int i; 5105 // possibly extend whichGenerator 5106 resizeWhichGenerator(numberOld, numberOld+numberCuts); 5107 for ( i = 0;i<numberCuts;i++) { 5108 const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ; 5109 if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) { 5110 if (messageHandler()->logLevel()>2) 5111 printf("Old cut added - violation %g\n", 5112 thisCut->violated(cbcColSolution_)) ; 5113 whichGenerator_[numberOld++]=-1; 5114 theseCuts.insert(*thisCut) ; 5115 } 5116 } 5117 } 5118 delete [] eachCuts; 5119 } 5120 //if (!feasible) 5121 //break; 5122 /* 5123 End of the loop to exercise each generator - try heuristics 5124 - unless at root node and first pass 5125 */ 5126 if (numberNodes_||currentPassNumber_!=1) { 5127 double * newSolution = new double [numberColumns] ; 5128 double heuristicValue = getCutoff() ; 5129 int found = -1; // no solution found 5130 for (i = 0;i<numberHeuristics_;i++) { 4317 5131 // see if heuristic will do anything 4318 5132 double saveValue = heuristicValue ; 4319 5133 int ifSol = 4320 heuristic_[i -numberCutGenerators_]->solution(heuristicValue,4321 4322 5134 heuristic_[i]->solution(heuristicValue, 5135 newSolution, 5136 theseCuts) ; 4323 5137 if (ifSol>0) { 4324 5138 // better solution found 4325 found = i -numberCutGenerators_;5139 found = i ; 4326 5140 incrementUsed(newSolution); 4327 5141 } else if (ifSol<0) { … … 4329 5143 } 4330 5144 } 4331 int numberRowCutsAfter = theseCuts.sizeRowCuts() ; 4332 int numberColumnCutsAfter = theseCuts.sizeColCuts() ; 4333 4334 if ((specialOptions_&1)!=0) { 4335 if (onOptimalPath) { 4336 int k ; 4337 for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) { 4338 OsiRowCut thisCut = theseCuts.rowCut(k) ; 4339 if(debugger->invalidCut(thisCut)) { 4340 solver_->writeMps("badCut"); 4341 } 4342 assert(!debugger->invalidCut(thisCut)) ; 4343 } 4344 } 4345 } 4346 /* 4347 The cut generator/heuristic has done its thing, and maybe it generated some 4348 cuts and/or a new solution. Do a bit of bookkeeping: load 4349 whichGenerator[i] with the index of the generator responsible for a cut, 4350 and place cuts flagged as global in the global cut pool for the model. 4351 4352 lastNumberCuts is the sum of cuts added in previous iterations; it's the 4353 offset to the proper starting position in whichGenerator. 4354 */ 4355 int numberBefore = 4356 numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ; 4357 int numberAfter = 4358 numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ; 4359 // possibly extend whichGenerator 4360 resizeWhichGenerator(numberBefore, numberAfter); 4361 int j ; 4362 if (fullScan) { 4363 // counts 4364 countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ; 4365 } 4366 countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ; 4367 4368 for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) { 4369 whichGenerator_[numberBefore++] = i ; 4370 const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ; 4371 if (thisCut->lb()>thisCut->ub()) 4372 violated=-2; // sub-problem is infeasible 4373 if (thisCut->globallyValid()) { 4374 // add to global list 4375 OsiRowCut newCut(*thisCut); 4376 newCut.setGloballyValid(2); 4377 newCut.mutableRow().setTestForDuplicateIndex(false); 4378 globalCuts_.insert(newCut) ; 4379 } 4380 } 4381 for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) { 4382 whichGenerator_[numberBefore++] = i ; 4383 const OsiColCut * thisCut = theseCuts.colCutPtr(j) ; 4384 if (thisCut->globallyValid()) { 4385 // add to global list 4386 OsiColCut newCut(*thisCut); 4387 newCut.setGloballyValid(2); 4388 globalCuts_.insert(newCut) ; 4389 } 4390 } 4391 } 4392 // If at root node and first pass do heuristics without cuts 4393 if (!numberNodes_&¤tPassNumber_==1) { 4394 // Save number solutions 4395 int saveNumberSolutions = numberSolutions_; 4396 int saveNumberHeuristicSolutions = numberHeuristicSolutions_; 4397 // make sure bounds tight 4398 { 4399 for (int i = 0;i < numberObjects_;i++) 4400 object_[i]->resetBounds(solver_) ; 4401 } 4402 for (int i = 0;i<numberHeuristics_;i++) { 4403 // see if heuristic will do anything 4404 double saveValue = heuristicValue ; 4405 int ifSol = heuristic_[i]->solution(heuristicValue, 4406 newSolution); 4407 if (ifSol>0) { 4408 // better solution found 4409 found = i ; 4410 incrementUsed(newSolution); 4411 // increment number of solutions so other heuristics can test 4412 numberSolutions_++; 4413 numberHeuristicSolutions_++; 4414 } else { 4415 heuristicValue = saveValue ; 4416 } 4417 } 4418 // Restore number solutions 4419 numberSolutions_ = saveNumberSolutions; 4420 numberHeuristicSolutions_ = saveNumberHeuristicSolutions; 4421 } 4422 // Add in any violated saved cuts 4423 if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) { 4424 int numberOld = theseCuts.sizeRowCuts(); 4425 int numberCuts = slackCuts.sizeRowCuts() ; 4426 int i; 4427 // possibly extend whichGenerator 4428 resizeWhichGenerator(numberOld, numberOld+numberCuts); 4429 for ( i = 0;i<numberCuts;i++) { 4430 const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ; 4431 if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) { 4432 if (messageHandler()->logLevel()>2) 4433 printf("Old cut added - violation %g\n", 4434 thisCut->violated(cbcColSolution_)) ; 4435 whichGenerator_[numberOld++]=-1; 4436 theseCuts.insert(*thisCut) ; 4437 } 4438 } 4439 } 4440 /* 4441 End of the loop to exercise each generator/heuristic. 4442 5145 /* 4443 5146 Did any of the heuristics turn up a new solution? Record it before we free 4444 5147 the vector. 4445 5148 */ 4446 if (found >= 0) {4447 4448 4449 setBestSolution(CBC_ROUNDING,heuristicValue,newSolution);4450 lastHeuristic_ = heuristic_[found];4451 5149 if (found >= 0) { 5150 phase_=4; 5151 incrementUsed(newSolution); 5152 lastHeuristic_ = heuristic_[found]; 5153 setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ; 5154 CbcTreeLocal * tree 4452 5155 = dynamic_cast<CbcTreeLocal *> (tree_); 4453 if (tree) 4454 tree->passInSolution(bestSolution_,heuristicValue); 4455 } 4456 delete [] newSolution ; 5156 if (tree) 5157 tree->passInSolution(bestSolution_,heuristicValue); 5158 } 5159 delete [] newSolution ; 5160 } 4457 5161 4458 5162 #if 0 … … 4616 5320 } 4617 5321 feasible = resolve(node ? node->nodeInfo() : NULL,2) ; 4618 if ( CoinCpuTime()-dblParam_[CbcStartSeconds]> dblParam_[CbcMaximumSeconds] )5322 if ( getCurrentSeconds() > dblParam_[CbcMaximumSeconds] ) 4619 5323 numberTries=0; // exit 4620 5324 # ifdef CBC_DEBUG … … 4689 5393 if (!feasible) 4690 5394 { int i ; 4691 for (i = 0;i<currentNumberCuts_;i++) { 4692 // take off node 4693 if (addedCuts_[i]) { 4694 if (!addedCuts_[i]->decrement()) 4695 delete addedCuts_[i] ; 4696 addedCuts_[i] = NULL ; 5395 if (currentNumberCuts_) { 5396 lockThread(); 5397 for (i = 0;i<currentNumberCuts_;i++) { 5398 // take off node 5399 if (addedCuts_[i]) { 5400 if (!addedCuts_[i]->decrement()) 5401 delete addedCuts_[i] ; 5402 addedCuts_[i] = NULL ; 5403 } 4697 5404 } 4698 5405 } … … 4830 5537 phase_=4; 4831 5538 incrementUsed(newSolution); 5539 lastHeuristic_ = heuristic_[found]; 4832 5540 setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ; 4833 lastHeuristic_ = heuristic_[found];4834 5541 } 4835 5542 delete [] newSolution ; … … 5165 5872 clpSolver->setSpecialOptions(saveClpOptions); 5166 5873 # endif 5874 #ifdef CBC_THREAD 5875 if (threadModel) { 5876 // stop threads 5877 int i; 5878 for (i=0;i<numberThreads_;i++) { 5879 while (threadInfo[i].returnCode==0) { 5880 pthread_cond_signal(threadInfo[i].condition2); // unlock 5881 pthread_mutex_lock(&condition_mutex); 5882 struct timespec absTime;