Changeset 1286 for branches/sandbox/Cbc/src/CbcHeuristicLocal.cpp
 Timestamp:
 Nov 9, 2009 6:33:07 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/sandbox/Cbc/src/CbcHeuristicLocal.cpp
r1271 r1286 20 20 21 21 // Default Constructor 22 CbcHeuristicLocal::CbcHeuristicLocal() 23 :CbcHeuristic()24 { 25 numberSolutions_=0;26 swap_=0;27 used_=NULL;22 CbcHeuristicLocal::CbcHeuristicLocal() 23 : CbcHeuristic() 24 { 25 numberSolutions_ = 0; 26 swap_ = 0; 27 used_ = NULL; 28 28 } 29 29 … … 31 31 32 32 CbcHeuristicLocal::CbcHeuristicLocal(CbcModel & model) 33 :CbcHeuristic(model)34 { 35 numberSolutions_=0;36 swap_=0;37 // Get a copy of original matrix38 assert(model.solver());39 if (model.solver()>getNumRows()) {40 matrix_ = *model.solver()>getMatrixByCol();41 }42 int numberColumns = model.solver()>getNumCols();43 used_ = new int[numberColumns];44 memset(used_,0,numberColumns*sizeof(int));45 } 46 47 // Destructor 33 : CbcHeuristic(model) 34 { 35 numberSolutions_ = 0; 36 swap_ = 0; 37 // Get a copy of original matrix 38 assert(model.solver()); 39 if (model.solver()>getNumRows()) { 40 matrix_ = *model.solver()>getMatrixByCol(); 41 } 42 int numberColumns = model.solver()>getNumCols(); 43 used_ = new int[numberColumns]; 44 memset(used_, 0, numberColumns*sizeof(int)); 45 } 46 47 // Destructor 48 48 CbcHeuristicLocal::~CbcHeuristicLocal () 49 49 { 50 delete [] used_;50 delete [] used_; 51 51 } 52 52 … … 55 55 CbcHeuristicLocal::clone() const 56 56 { 57 return new CbcHeuristicLocal(*this);57 return new CbcHeuristicLocal(*this); 58 58 } 59 59 // Create C++ lines to get to current state 60 void 61 CbcHeuristicLocal::generateCpp( FILE * fp) 62 { 63 CbcHeuristicLocal other;64 fprintf(fp,"0#include \"CbcHeuristicLocal.hpp\"\n");65 fprintf(fp,"3 CbcHeuristicLocal heuristicLocal(*cbcModel);\n");66 CbcHeuristic::generateCpp(fp,"heuristicLocal");67 if (swap_!=other.swap_)68 fprintf(fp,"3 heuristicLocal.setSearchType(%d);\n",swap_);69 else70 fprintf(fp,"4 heuristicLocal.setSearchType(%d);\n",swap_);71 fprintf(fp,"3 cbcModel>addHeuristic(&heuristicLocal);\n");72 } 73 74 // Copy constructor 60 void 61 CbcHeuristicLocal::generateCpp( FILE * fp) 62 { 63 CbcHeuristicLocal other; 64 fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n"); 65 fprintf(fp, "3 CbcHeuristicLocal heuristicLocal(*cbcModel);\n"); 66 CbcHeuristic::generateCpp(fp, "heuristicLocal"); 67 if (swap_ != other.swap_) 68 fprintf(fp, "3 heuristicLocal.setSearchType(%d);\n", swap_); 69 else 70 fprintf(fp, "4 heuristicLocal.setSearchType(%d);\n", swap_); 71 fprintf(fp, "3 cbcModel>addHeuristic(&heuristicLocal);\n"); 72 } 73 74 // Copy constructor 75 75 CbcHeuristicLocal::CbcHeuristicLocal(const CbcHeuristicLocal & rhs) 76 : 77 CbcHeuristic(rhs), 78 matrix_(rhs.matrix_), 79 numberSolutions_(rhs.numberSolutions_), 80 swap_(rhs.swap_) 81 { 82 if (model_&&rhs.used_) { 83 int numberColumns = model_>solver()>getNumCols(); 84 used_ = CoinCopyOfArray(rhs.used_,numberColumns); 85 } else { 86 used_=NULL; 87 } 88 } 89 90 // Assignment operator 91 CbcHeuristicLocal & 92 CbcHeuristicLocal::operator=( const CbcHeuristicLocal& rhs) 93 { 94 if (this!=&rhs) { 95 CbcHeuristic::operator=(rhs); 96 matrix_ = rhs.matrix_; 97 numberSolutions_ = rhs.numberSolutions_; 98 swap_ = rhs.swap_; 76 : 77 CbcHeuristic(rhs), 78 matrix_(rhs.matrix_), 79 numberSolutions_(rhs.numberSolutions_), 80 swap_(rhs.swap_) 81 { 82 if (model_ && rhs.used_) { 83 int numberColumns = model_>solver()>getNumCols(); 84 used_ = CoinCopyOfArray(rhs.used_, numberColumns); 85 } else { 86 used_ = NULL; 87 } 88 } 89 90 // Assignment operator 91 CbcHeuristicLocal & 92 CbcHeuristicLocal::operator=( const CbcHeuristicLocal & rhs) 93 { 94 if (this != &rhs) { 95 CbcHeuristic::operator=(rhs); 96 matrix_ = rhs.matrix_; 97 numberSolutions_ = rhs.numberSolutions_; 98 swap_ = rhs.swap_; 99 delete [] used_; 100 if (model_ && rhs.used_) { 101 int numberColumns = model_>solver()>getNumCols(); 102 used_ = CoinCopyOfArray(rhs.used_, numberColumns); 103 } else { 104 used_ = NULL; 105 } 106 } 107 return *this; 108 } 109 110 // Resets stuff if model changes 111 void 112 CbcHeuristicLocal::resetModel(CbcModel * /*model*/) 113 { 114 //CbcHeuristic::resetModel(model); 99 115 delete [] used_; 100 if (model_&&rhs.used_) { 101 int numberColumns = model_>solver()>getNumCols(); 102 used_ = CoinCopyOfArray(rhs.used_,numberColumns); 116 if (model_ && used_) { 117 int numberColumns = model_>solver()>getNumCols(); 118 used_ = new int[numberColumns]; 119 memset(used_, 0, numberColumns*sizeof(int)); 103 120 } else { 104 used_=NULL; 105 } 106 } 107 return *this; 108 } 109 110 // Resets stuff if model changes 111 void 112 CbcHeuristicLocal::resetModel(CbcModel * /*model*/) 113 { 114 //CbcHeuristic::resetModel(model); 115 delete [] used_; 116 if (model_&&used_) { 117 int numberColumns = model_>solver()>getNumCols(); 118 used_ = new int[numberColumns]; 119 memset(used_,0,numberColumns*sizeof(int)); 120 } else { 121 used_=NULL; 122 } 121 used_ = NULL; 122 } 123 123 } 124 124 // This version fixes stuff and does IP 125 int 125 int 126 126 CbcHeuristicLocal::solutionFix(double & objectiveValue, 127 128 129 { 130 numCouldRun_++;131 // See if to do132 if (!when()(when()==1&&model_>phase()!=1))133 return 0; // switched off134 // Don't do if it was this heuristic which found solution!135 if (this==model_>lastHeuristic())136 return 0;137 OsiSolverInterface * newSolver = model_>continuousSolver()>clone();138 const double * colLower = newSolver>getColLower();139 //const double * colUpper = newSolver>getColUpper();140 141 int numberIntegers = model_>numberIntegers();142 const int * integerVariable = model_>integerVariable();143 144 int i;145 int nFix=0;146 for (i=0;i<numberIntegers;i++) {147 int iColumn=integerVariable[i];148 const OsiObject * object = model_>object(i);149 // get original bounds150 double originalLower;151 double originalUpper;152 getIntegerInformation( object,originalLower, originalUpper);153 newSolver>setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));154 if (!used_[iColumn]) {155 newSolver>setColUpper(iColumn,colLower[iColumn]);156 nFix++;157 }158 }159 int returnCode = 0;160 if (nFix*10>numberIntegers) {161 returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,objectiveValue,162 objectiveValue,"CbcHeuristicLocal");163 if (returnCode<0) {164 returnCode=0; // returned on size165 int numberColumns=newSolver>getNumCols();166 int numberContinuous = numberColumnsnumberIntegers;167 if (numberContinuous>2*numberIntegers&&168 nFix*10<numberColumns) {127 double * newSolution, 128 const int * /*keep*/) 129 { 130 numCouldRun_++; 131 // See if to do 132 if (!when()  (when() == 1 && model_>phase() != 1)) 133 return 0; // switched off 134 // Don't do if it was this heuristic which found solution! 135 if (this == model_>lastHeuristic()) 136 return 0; 137 OsiSolverInterface * newSolver = model_>continuousSolver()>clone(); 138 const double * colLower = newSolver>getColLower(); 139 //const double * colUpper = newSolver>getColUpper(); 140 141 int numberIntegers = model_>numberIntegers(); 142 const int * integerVariable = model_>integerVariable(); 143 144 int i; 145 int nFix = 0; 146 for (i = 0; i < numberIntegers; i++) { 147 int iColumn = integerVariable[i]; 148 const OsiObject * object = model_>object(i); 149 // get original bounds 150 double originalLower; 151 double originalUpper; 152 getIntegerInformation( object, originalLower, originalUpper); 153 newSolver>setColLower(iColumn, CoinMax(colLower[iColumn], originalLower)); 154 if (!used_[iColumn]) { 155 newSolver>setColUpper(iColumn, colLower[iColumn]); 156 nFix++; 157 } 158 } 159 int returnCode = 0; 160 if (nFix*10 > numberIntegers) { 161 returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, objectiveValue, 162 objectiveValue, "CbcHeuristicLocal"); 163 if (returnCode < 0) { 164 returnCode = 0; // returned on size 165 int numberColumns = newSolver>getNumCols(); 166 int numberContinuous = numberColumns  numberIntegers; 167 if (numberContinuous > 2*numberIntegers && 168 nFix*10 < numberColumns) { 169 169 #define LOCAL_FIX_CONTINUOUS 170 170 #ifdef LOCAL_FIX_CONTINUOUS 171 172 173 int nAtLb=0;174 175 176 double direction=newSolver>getObjSense();177 for (int iColumn=0;iColumn<numberColumns;iColumn++) {178 179 180 181 182 183 184 185 186 187 188 189 190 191 int nFix2=0;192 for (int iColumn=0;iColumn<numberColumns;iColumn++) {193 194 195 double djValue = dj[iColumn]*direction;196 if (djValue>1.0e6) {197 sort[nFix2]=djValue;198 which[nFix2++]=iColumn;199 200 201 202 203 CoinSort_2(sort,sort+nFix2,which);204 int divisor=2;205 nFix2=CoinMin(nFix2,(numberColumnsnFix)/divisor);206 for (int i=0;i<nFix2;i++) {207 208 newSolver>setColUpper(iColumn,colLower[iColumn]);209 210 211 171 //const double * colUpper = newSolver>getColUpper(); 172 const double * colLower = newSolver>getColLower(); 173 int nAtLb = 0; 174 //double sumDj=0.0; 175 const double * dj = newSolver>getReducedCost(); 176 double direction = newSolver>getObjSense(); 177 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 178 if (!newSolver>isInteger(iColumn)) { 179 if (!used_[iColumn]) { 180 //double djValue = dj[iColumn]*direction; 181 nAtLb++; 182 //sumDj += djValue; 183 } 184 } 185 } 186 if (nAtLb) { 187 // fix some continuous 188 double * sort = new double[nAtLb]; 189 int * which = new int [nAtLb]; 190 //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e6); 191 int nFix2 = 0; 192 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 193 if (!newSolver>isInteger(iColumn)) { 194 if (!used_[iColumn]) { 195 double djValue = dj[iColumn] * direction; 196 if (djValue > 1.0e6) { 197 sort[nFix2] = djValue; 198 which[nFix2++] = iColumn; 199 } 200 } 201 } 202 } 203 CoinSort_2(sort, sort + nFix2, which); 204 int divisor = 2; 205 nFix2 = CoinMin(nFix2, (numberColumns  nFix) / divisor); 206 for (int i = 0; i < nFix2; i++) { 207 int iColumn = which[i]; 208 newSolver>setColUpper(iColumn, colLower[iColumn]); 209 } 210 delete [] sort; 211 delete [] which; 212 212 #ifdef CLP_INVESTIGATE2 213 214 nFix,nFix2);213 printf("%d integers have zero value, and %d continuous fixed at lb\n", 214 nFix, nFix2); 215 215 #endif 216 217 numberNodes_,newSolution,218 219 objectiveValue,"CbcHeuristicLocal");220 if (returnCode<0) 221 returnCode=0; // returned on size222 216 returnCode = smallBranchAndBound(newSolver, 217 numberNodes_, newSolution, 218 objectiveValue, 219 objectiveValue, "CbcHeuristicLocal"); 220 if (returnCode < 0) 221 returnCode = 0; // returned on size 222 } 223 223 #endif 224 }225 }226 }227 if ((returnCode&2)!=0) {228 // could add cut229 returnCode &= ~2;230 }231 232 delete newSolver;233 return returnCode;224 } 225 } 226 } 227 if ((returnCode&2) != 0) { 228 // could add cut 229 returnCode &= ~2; 230 } 231 232 delete newSolver; 233 return returnCode; 234 234 } 235 235 /* … … 239 239 int 240 240 CbcHeuristicLocal::solution(double & solutionValue, 241 242 { 243 244 numCouldRun_++;245 if (numberSolutions_==model_>getSolutionCount())246 return 0;247 numberSolutions_=model_>getSolutionCount();248 if ((model_>getNumCols()>100000&&model_>getNumCols()>249 10*model_>getNumRows())numberSolutions_<=1)250 return 0; // probably not worth it251 // worth trying252 253 OsiSolverInterface * solver = model_>solver();254 const double * rowLower = solver>getRowLower();255 const double * rowUpper = solver>getRowUpper();256 const double * solution = model_>bestSolution();257 if (!solution)258 return 0; // No solution found yet259 const double * objective = solver>getObjCoefficients();260 double primalTolerance;261 solver>getDblParam(OsiPrimalTolerance,primalTolerance);262 263 int numberRows = matrix_.getNumRows();264 265 int numberIntegers = model_>numberIntegers();266 const int * integerVariable = model_>integerVariable();267 268 int i;269 double direction = solver>getObjSense();270 double newSolutionValue = model_>getObjValue()*direction;271 int returnCode = 0;272 numRuns_++;273 // Column copy274 const double * element = matrix_.getElements();275 const int * row = matrix_.getIndices();276 const CoinBigIndex * columnStart = matrix_.getVectorStarts();277 const int * columnLength = matrix_.getVectorLengths();278 279 // Get solution array for heuristic solution280 int numberColumns = solver>getNumCols();281 double * newSolution = new double [numberColumns];282 memcpy(newSolution,solution,numberColumns*sizeof(double));241 double * betterSolution) 242 { 243 244 numCouldRun_++; 245 if (numberSolutions_ == model_>getSolutionCount()) 246 return 0; 247 numberSolutions_ = model_>getSolutionCount(); 248 if ((model_>getNumCols() > 100000 && model_>getNumCols() > 249 10*model_>getNumRows())  numberSolutions_ <= 1) 250 return 0; // probably not worth it 251 // worth trying 252 253 OsiSolverInterface * solver = model_>solver(); 254 const double * rowLower = solver>getRowLower(); 255 const double * rowUpper = solver>getRowUpper(); 256 const double * solution = model_>bestSolution(); 257 if (!solution) 258 return 0; // No solution found yet 259 const double * objective = solver>getObjCoefficients(); 260 double primalTolerance; 261 solver>getDblParam(OsiPrimalTolerance, primalTolerance); 262 263 int numberRows = matrix_.getNumRows(); 264 265 int numberIntegers = model_>numberIntegers(); 266 const int * integerVariable = model_>integerVariable(); 267 268 int i; 269 double direction = solver>getObjSense(); 270 double newSolutionValue = model_>getObjValue() * direction; 271 int returnCode = 0; 272 numRuns_++; 273 // Column copy 274 const double * element = matrix_.getElements(); 275 const int * row = matrix_.getIndices(); 276 const CoinBigIndex * columnStart = matrix_.getVectorStarts(); 277 const int * columnLength = matrix_.getVectorLengths(); 278 279 // Get solution array for heuristic solution 280 int numberColumns = solver>getNumCols(); 281 double * newSolution = new double [numberColumns]; 282 memcpy(newSolution, solution, numberColumns*sizeof(double)); 283 283 #ifdef LOCAL_FIX_CONTINUOUS 284 // mark continuous used285 const double * columnLower = solver>getColLower();286 for (int iColumn=0;iColumn<numberColumns;iColumn++) {287 if (!solver>isInteger(iColumn)) {288 if (solution[iColumn]>columnLower[iColumn]+1.0e8)289 used_[iColumn]=numberSolutions_;290 }291 }284 // mark continuous used 285 const double * columnLower = solver>getColLower(); 286 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 287 if (!solver>isInteger(iColumn)) { 288 if (solution[iColumn] > columnLower[iColumn] + 1.0e8) 289 used_[iColumn] = numberSolutions_; 290 } 291 } 292 292 #endif 293 293 294 // way is 1 if down possible, 2 if up possible, 3 if both possible295 char * way = new char[numberIntegers];296 // corrected costs297 double * cost = new double[numberIntegers];298 // for array to mark infeasible rows after iColumn branch299 char * mark = new char[numberRows];300 memset(mark,0,numberRows);301 // space to save values so we don't introduce rounding errors302 double * save = new double[numberRows];303 304 // clean solution305 for (i=0;i<numberIntegers;i++) {306 int iColumn = integerVariable[i];307 const OsiObject * object = model_>object(i);308 // get original bounds309 double originalLower;310 double originalUpper;311 getIntegerInformation( object,originalLower, originalUpper);312 double value=newSolution[iColumn];313 if (value<originalLower) {314 value=originalLower;315 newSolution[iColumn]=value;316 } else if (value>originalUpper) {317 value=originalUpper;318 newSolution[iColumn]=value;319 }320 double nearest=floor(value+0.5);321 //assert(fabs(valuenearest)<10.0*primalTolerance);322 value=nearest;323 newSolution[iColumn]=nearest;324 // if away from lower bound mark that fact325 if (nearest>originalLower) {326 used_[iColumn]=numberSolutions_;327 }328 cost[i] = direction*objective[iColumn];329 int iway=0;330 331 if (value>originalLower+0.5)332 iway = 1;333 if (value<originalUpper0.5)334 iway = 2;335 way[i]=static_cast<char>(iway);336 }337 // get row activities338 double * rowActivity = new double[numberRows];339 memset(rowActivity,0,numberRows*sizeof(double));340 341 for (i=0;i<numberColumns;i++) {342 int j;343 double value = newSolution[i];344 if (value) {345 for (j=columnStart[i];346 j<columnStart[i]+columnLength[i];j++) {347 int iRow=row[j];348 rowActivity[iRow] += value*element[j];349 }350 }351 }352 // check was feasible  if not adjust (cleaning may move)353 // if very infeasible then give up354 bool tryHeuristic=true;355 for (i=0;i<numberRows;i++) {356 if(rowActivity[i]<rowLower[i]) {357 if (rowActivity[i]<rowLower[i]10.0*primalTolerance)358 tryHeuristic=false;359 rowActivity[i]=rowLower[i];360 } else if(rowActivity[i]>rowUpper[i]) {361 if (rowActivity[i]<rowUpper[i]+10.0*primalTolerance)362 tryHeuristic=false;363 rowActivity[i]=rowUpper[i];364 }365 }366 // Switch off if may take too long367 if (model_>getNumCols()>10000&&model_>getNumCols()>368 10*model_>getNumRows())369 tryHeuristic=false;370 if (tryHeuristic) {371 372 // best change in objective373 double bestChange=0.0;374 375 for (i=0;i<numberIntegers;i++) {376 int iColumn = integerVariable[i];377 378 double objectiveCoefficient = cost[i];379 int k;380 int j;381 int goodK=1;382 int wayK=1,wayI=1;383 if ((way[i]&1)!=0) {384 int numberInfeasible=0;385 386 for (j=columnStart[iColumn];387 j<columnStart[iColumn]+columnLength[iColumn];j++) {388 389 save[iRow]=rowActivity[iRow];390 391 if(rowActivity[iRow]<rowLower[iRow]primalTolerance392 rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {393 394 mark[iRow]=1;395 396 397 398 399 for (k=i+1;k<numberIntegers;k++) {400 if ((way[k]&1)!=0) {401 402 if (objectiveCoefficientcost[k]<bestChange) {403 404 bool good=true;405 int numberMarked=0;406 407 for (j=columnStart[kColumn];408 j<columnStart[kColumn]+columnLength[kColumn];j++) {409 410 411 if(newValue<rowLower[iRow]primalTolerance412 newValue>rowUpper[iRow]+primalTolerance) {413 good=false;414 415 416 417 418 419 420 if (good&&numberMarked==numberInfeasible) {421 422 goodK=k;423 wayK=1;424 wayI=1;425 bestChange = objectiveCoefficientcost[k];426 427 428 429 if ((way[k]&2)!=0) {430 431 if (objectiveCoefficient+cost[k]<bestChange) {432 433 bool good=true;434 int numberMarked=0;435 436 for (j=columnStart[kColumn];437 j<columnStart[kColumn]+columnLength[kColumn];j++) {438 439 440 if(newValue<rowLower[iRow]primalTolerance441 newValue>rowUpper[iRow]+primalTolerance) {442 good=false;443 444 445 446 447 448 449 if (good&&numberMarked==numberInfeasible) {450 451 goodK=k;452 wayK=1;453 wayI=1;454 bestChange = objectiveCoefficient+cost[k];455 456 457 458 459 460 for (j=columnStart[iColumn];461 j<columnStart[iColumn]+columnLength[iColumn];j++) {462 463 464 mark[iRow]=0;465 466 }467 if ((way[i]&2)!=0) {468 int numberInfeasible=0;469 470 for (j=columnStart[iColumn];471 j<columnStart[iColumn]+columnLength[iColumn];j++) {472 473 save[iRow]=rowActivity[iRow];474 475 if(rowActivity[iRow]<rowLower[iRow]primalTolerance476 rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {477 478 mark[iRow]=1;479 480 481 482 483 for (k=i+1;k<numberIntegers;k++) {484 if ((way[k]&1)!=0) {485 486 if (objectiveCoefficientcost[k]<bestChange) {487 488 bool good=true;489 int numberMarked=0;490 491 for (j=columnStart[kColumn];492 j<columnStart[kColumn]+columnLength[kColumn];j++) {493 494 495 if(newValue<rowLower[iRow]primalTolerance496 newValue>rowUpper[iRow]+primalTolerance) {497 good=false;498 499 500 501 502 503 504 if (good&&numberMarked==numberInfeasible) {505 506 goodK=k;507 wayK=1;508 wayI=1;509 bestChange = objectiveCoefficientcost[k];510 511 512 513 if ((way[k]&2)!=0) {514 515 if (objectiveCoefficient+cost[k]<bestChange) {516 517 bool good=true;518 int numberMarked=0;519 520 for (j=columnStart[kColumn];521 j<columnStart[kColumn]+columnLength[kColumn];j++) {522 523 524 if(newValue<rowLower[iRow]primalTolerance525 newValue>rowUpper[iRow]+primalTolerance) {526 good=false;527 528 529 530 531 532 533 if (good&&numberMarked==numberInfeasible) {534 535 goodK=k;536 wayK=1;537 wayI=1;538 bestChange = objectiveCoefficient+cost[k];539 540 541 542 543 544 for (j=columnStart[iColumn];545 j<columnStart[iColumn]+columnLength[iColumn];j++) {546 547 548 mark[iRow]=0;549 550 }551 if (goodK>=0) {552 553 for (j=columnStart[iColumn];554 j<columnStart[iColumn]+columnLength[iColumn];j++) {555 556 557 558 559 560 for (j=columnStart[kColumn];561 j<columnStart[kColumn]+columnLength[kColumn];j++) {562 563 564 565 566 567 568 569 570 571 getIntegerInformation( object,originalLower, originalUpper); 572 573 double value=newSolution[kColumn];574 int iway=0;575 576 if (value>originalLower+0.5) 577 578 if (value<originalUpper0.5) 579 580 way[goodK]=static_cast<char>(iway);581 }582 }583 if (bestChange+newSolutionValue<solutionValue) {584 // paranoid check585 memset(rowActivity,0,numberRows*sizeof(double));586 587 for (i=0;i<numberColumns;i++) {588 589 590 591 for (j=columnStart[i];592 j<columnStart[i]+columnLength[i];j++) {593 int iRow=row[j];594 rowActivity[iRow] += value*element[j];595 596 597 }598 int numberBad=0;599 double sumBad=0.0;600 // check was approximately feasible601 for (i=0;i<numberRows;i++) {602 if(rowActivity[i]<rowLower[i]) {603 sumBad += rowLower[i]rowActivity[i];604 if (rowActivity[i]<rowLower[i]10.0*primalTolerance)605 numberBad++;606 } else if(rowActivity[i]>rowUpper[i]) {607 sumBad += rowUpper[i]rowActivity[i];608 if (rowActivity[i]>rowUpper[i]+10.0*primalTolerance)609 numberBad++;610 611 }612 if (!numberBad) {613 for (i=0;i<numberIntegers;i++) {614 int iColumn = integerVariable[i];615 const OsiObject * object = model_>object(i);616 // get original bounds617 618 619 getIntegerInformation( object,originalLower, originalUpper); 620 621 double value=newSolution[iColumn];622 // if away from lower bound mark that fact623 if (value>originalLower) {624 used_[iColumn]=numberSolutions_;625 }626 }627 // new solution628 memcpy(betterSolution,newSolution,numberColumns*sizeof(double));629 630 631 632 633 634 635 returnCode=1;636 solutionValue = newSolutionValue + bestChange;637 } else {638 // bad solution  should not happen so debug if see message639 printf("Local search got bad solution with %d infeasibilities summing to %g\n",640 numberBad,sumBad);641 }642 }643 }644 delete [] newSolution;645 delete [] rowActivity;646 delete [] way;647 delete [] cost;648 delete [] save;649 delete [] mark;650 if (numberSolutions_>1&&swap_==1) {651 // try merge652 int returnCode2=solutionFix( solutionValue, betterSolution,NULL);653 if (returnCode2)654 returnCode=1;655 }656 return returnCode;294 // way is 1 if down possible, 2 if up possible, 3 if both possible 295 char * way = new char[numberIntegers]; 296 // corrected costs 297 double * cost = new double[numberIntegers]; 298 // for array to mark infeasible rows after iColumn branch 299 char * mark = new char[numberRows]; 300 memset(mark, 0, numberRows); 301 // space to save values so we don't introduce rounding errors 302 double * save = new double[numberRows]; 303 304 // clean solution 305 for (i = 0; i < numberIntegers; i++) { 306 int iColumn = integerVariable[i]; 307 const OsiObject * object = model_>object(i); 308 // get original bounds 309 double originalLower; 310 double originalUpper; 311 getIntegerInformation( object, originalLower, originalUpper); 312 double value = newSolution[iColumn]; 313 if (value < originalLower) { 314 value = originalLower; 315 newSolution[iColumn] = value; 316 } else if (value > originalUpper) { 317 value = originalUpper; 318 newSolution[iColumn] = value; 319 } 320 double nearest = floor(value + 0.5); 321 //assert(fabs(valuenearest)<10.0*primalTolerance); 322 value = nearest; 323 newSolution[iColumn] = nearest; 324 // if away from lower bound mark that fact 325 if (nearest > originalLower) { 326 used_[iColumn] = numberSolutions_; 327 } 328 cost[i] = direction * objective[iColumn]; 329 int iway = 0; 330 331 if (value > originalLower + 0.5) 332 iway = 1; 333 if (value < originalUpper  0.5) 334 iway = 2; 335 way[i] = static_cast<char>(iway); 336 } 337 // get row activities 338 double * rowActivity = new double[numberRows]; 339 memset(rowActivity, 0, numberRows*sizeof(double)); 340 341 for (i = 0; i < numberColumns; i++) { 342 int j; 343 double value = newSolution[i]; 344 if (value) { 345 for (j = columnStart[i]; 346 j < columnStart[i] + columnLength[i]; j++) { 347 int iRow = row[j]; 348 rowActivity[iRow] += value * element[j]; 349 } 350 } 351 } 352 // check was feasible  if not adjust (cleaning may move) 353 // if very infeasible then give up 354 bool tryHeuristic = true; 355 for (i = 0; i < numberRows; i++) { 356 if (rowActivity[i] < rowLower[i]) { 357 if (rowActivity[i] < rowLower[i]  10.0*primalTolerance) 358 tryHeuristic = false; 359 rowActivity[i] = rowLower[i]; 360 } else if (rowActivity[i] > rowUpper[i]) { 361 if (rowActivity[i] < rowUpper[i] + 10.0*primalTolerance) 362 tryHeuristic = false; 363 rowActivity[i] = rowUpper[i]; 364 } 365 } 366 // Switch off if may take too long 367 if (model_>getNumCols() > 10000 && model_>getNumCols() > 368 10*model_>getNumRows()) 369 tryHeuristic = false; 370 if (tryHeuristic) { 371 372 // best change in objective 373 double bestChange = 0.0; 374 375 for (i = 0; i < numberIntegers; i++) { 376 int iColumn = integerVariable[i]; 377 378 double objectiveCoefficient = cost[i]; 379 int k; 380 int j; 381 int goodK = 1; 382 int wayK = 1, wayI = 1; 383 if ((way[i]&1) != 0) { 384 int numberInfeasible = 0; 385 // save row activities and adjust 386 for (j = columnStart[iColumn]; 387 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 388 int iRow = row[j]; 389 save[iRow] = rowActivity[iRow]; 390 rowActivity[iRow] = element[j]; 391 if (rowActivity[iRow] < rowLower[iRow]  primalTolerance  392 rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { 393 // mark row 394 mark[iRow] = 1; 395 numberInfeasible++; 396 } 397 } 398 // try down 399 for (k = i + 1; k < numberIntegers; k++) { 400 if ((way[k]&1) != 0) { 401 // try down 402 if (objectiveCoefficient  cost[k] < bestChange) { 403 // see if feasible down 404 bool good = true; 405 int numberMarked = 0; 406 int kColumn = integerVariable[k]; 407 for (j = columnStart[kColumn]; 408 j < columnStart[kColumn] + columnLength[kColumn]; j++) { 409 int iRow = row[j]; 410 double newValue = rowActivity[iRow]  element[j]; 411 if (newValue < rowLower[iRow]  primalTolerance  412 newValue > rowUpper[iRow] + primalTolerance) { 413 good = false; 414 break; 415 } else if (mark[iRow]) { 416 // made feasible 417 numberMarked++; 418 } 419 } 420 if (good && numberMarked == numberInfeasible) { 421 // better solution 422 goodK = k; 423 wayK = 1; 424 wayI = 1; 425 bestChange = objectiveCoefficient  cost[k]; 426 } 427 } 428 } 429 if ((way[k]&2) != 0) { 430 // try up 431 if (objectiveCoefficient + cost[k] < bestChange) { 432 // see if feasible up 433 bool good = true; 434 int numberMarked = 0; 435 int kColumn = integerVariable[k]; 436 for (j = columnStart[kColumn]; 437 j < columnStart[kColumn] + columnLength[kColumn]; j++) { 438 int iRow = row[j]; 439 double newValue = rowActivity[iRow] + element[j]; 440 if (newValue < rowLower[iRow]  primalTolerance  441 newValue > rowUpper[iRow] + primalTolerance) { 442 good = false; 443 break; 444 } else if (mark[iRow]) { 445 // made feasible 446 numberMarked++; 447 } 448 } 449 if (good && numberMarked == numberInfeasible) { 450 // better solution 451 goodK = k; 452 wayK = 1; 453 wayI = 1; 454 bestChange = objectiveCoefficient + cost[k]; 455 } 456 } 457 } 458 } 459 // restore row activities 460 for (j = columnStart[iColumn]; 461 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 462 int iRow = row[j]; 463 rowActivity[iRow] = save[iRow]; 464 mark[iRow] = 0; 465 } 466 } 467 if ((way[i]&2) != 0) { 468 int numberInfeasible = 0; 469 // save row activities and adjust 470 for (j = columnStart[iColumn]; 471 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 472 int iRow = row[j]; 473 save[iRow] = rowActivity[iRow]; 474 rowActivity[iRow] += element[j]; 475 if (rowActivity[iRow] < rowLower[iRow]  primalTolerance  476 rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { 477 // mark row 478 mark[iRow] = 1; 479 numberInfeasible++; 480 } 481 } 482 // try up 483 for (k = i + 1; k < numberIntegers; k++) { 484 if ((way[k]&1) != 0) { 485 // try down 486 if (objectiveCoefficient  cost[k] < bestChange) { 487 // see if feasible down 488 bool good = true; 489 int numberMarked = 0; 490 int kColumn = integerVariable[k]; 491 for (j = columnStart[kColumn]; 492 j < columnStart[kColumn] + columnLength[kColumn]; j++) { 493 int iRow = row[j]; 494 double newValue = rowActivity[iRow]  element[j]; 495 if (newValue < rowLower[iRow]  primalTolerance  496 newValue > rowUpper[iRow] + primalTolerance) { 497 good = false; 498 break; 499 } else if (mark[iRow]) { 500 // made feasible 501 numberMarked++; 502 } 503 } 504 if (good && numberMarked == numberInfeasible) { 505 // better solution 506 goodK = k; 507 wayK = 1; 508 wayI = 1; 509 bestChange = objectiveCoefficient  cost[k]; 510 } 511 } 512 } 513 if ((way[k]&2) != 0) { 514 // try up 515 if (objectiveCoefficient + cost[k] < bestChange) { 516 // see if feasible up 517 bool good = true; 518 int numberMarked = 0; 519 int kColumn = integerVariable[k]; 520 for (j = columnStart[kColumn]; 521 j < columnStart[kColumn] + columnLength[kColumn]; j++) { 522 int iRow = row[j]; 523 double newValue = rowActivity[iRow] + element[j]; 524 if (newValue < rowLower[iRow]  primalTolerance  525 newValue > rowUpper[iRow] + primalTolerance) { 526 good = false; 527 break; 528 } else if (mark[iRow]) { 529 // made feasible 530 numberMarked++; 531 } 532 } 533 if (good && numberMarked == numberInfeasible) { 534 // better solution 535 goodK = k; 536 wayK = 1; 537 wayI = 1; 538 bestChange = objectiveCoefficient + cost[k]; 539 } 540 } 541 } 542 } 543 // restore row activities 544 for (j = columnStart[iColumn]; 545 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 546 int iRow = row[j]; 547 rowActivity[iRow] = save[iRow]; 548 mark[iRow] = 0; 549 } 550 } 551 if (goodK >= 0) { 552 // we found something  update solution 553 for (j = columnStart[iColumn]; 554 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 555 int iRow = row[j]; 556 rowActivity[iRow] += wayI * element[j]; 557 } 558 newSolution[iColumn] += wayI; 559 int kColumn = integerVariable[goodK]; 560 for (j = columnStart[kColumn]; 561 j < columnStart[kColumn] + columnLength[kColumn]; j++) { 562 int iRow = row[j]; 563 rowActivity[iRow] += wayK * element[j]; 564 } 565 newSolution[kColumn] += wayK; 566 // See if k can go further ? 567 const OsiObject * object = model_>object(goodK); 568 // get original bounds 569 double originalLower; 570 double originalUpper; 571 getIntegerInformation( object, originalLower, originalUpper); 572 573 double value = newSolution[kColumn]; 574 int iway = 0; 575 576 if (value > originalLower + 0.5) 577 iway = 1; 578 if (value < originalUpper  0.5) 579 iway = 2; 580 way[goodK] = static_cast<char>(iway); 581 } 582 } 583 if (bestChange + newSolutionValue < solutionValue) { 584 // paranoid check 585 memset(rowActivity, 0, numberRows*sizeof(double)); 586 587 for (i = 0; i < numberColumns; i++) { 588 int j; 589 double value = newSolution[i]; 590 if (value) { 591 for (j = columnStart[i]; 592 j < columnStart[i] + columnLength[i]; j++) { 593 int iRow = row[j]; 594 rowActivity[iRow] += value * element[j]; 595 } 596 } 597 } 598 int numberBad = 0; 599 double sumBad = 0.0; 600 // check was approximately feasible 601 for (i = 0; i < numberRows; i++) { 602 if (rowActivity[i] < rowLower[i]) { 603 sumBad += rowLower[i]  rowActivity[i]; 604 if (rowActivity[i] < rowLower[i]  10.0*primalTolerance) 605 numberBad++; 606 } else if (rowActivity[i] > rowUpper[i]) { 607 sumBad += rowUpper[i]  rowActivity[i]; 608 if (rowActivity[i] > rowUpper[i] + 10.0*primalTolerance) 609 numberBad++; 610 } 611 } 612 if (!numberBad) { 613 for (i = 0; i < numberIntegers; i++) { 614 int iColumn = integerVariable[i]; 615 const OsiObject * object = model_>object(i); 616 // get original bounds 617 double originalLower; 618 double originalUpper; 619 getIntegerInformation( object, originalLower, originalUpper); 620 621 double value = newSolution[iColumn]; 622 // if away from lower bound mark that fact 623 if (value > originalLower) { 624 used_[iColumn] = numberSolutions_; 625 } 626 } 627 // new solution 628 memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); 629 CoinWarmStartBasis * basis = 630 dynamic_cast<CoinWarmStartBasis *>(solver>getWarmStart()) ; 631 if (basis) { 632 model_>setBestSolutionBasis(* basis); 633 delete basis; 634 } 635 returnCode = 1; 636 solutionValue = newSolutionValue + bestChange; 637 } else { 638 // bad solution  should not happen so debug if see message 639 printf("Local search got bad solution with %d infeasibilities summing to %g\n", 640 numberBad, sumBad); 641 } 642 } 643 } 644 delete [] newSolution; 645 delete [] rowActivity; 646 delete [] way; 647 delete [] cost; 648 delete [] save; 649 delete [] mark; 650 if (numberSolutions_ > 1 && swap_ == 1) { 651 // try merge 652 int returnCode2 = solutionFix( solutionValue, betterSolution, NULL); 653 if (returnCode2) 654 returnCode = 1; 655 } 656 return returnCode; 657 657 } 658 658 // update model 659 659 void CbcHeuristicLocal::setModel(CbcModel * model) 660 660 { 661 model_ = model;662 // Get a copy of original matrix663 assert(model_>solver());664 if (model_>solver()>getNumRows()) {665 matrix_ = *model_>solver()>getMatrixByCol();666 }667 delete [] used_;668 int numberColumns = model>solver()>getNumCols();669 used_ = new int[numberColumns];670 memset(used_,0,numberColumns*sizeof(int));661 model_ = model; 662 // Get a copy of original matrix 663 assert(model_>solver()); 664 if (model_>solver()>getNumRows()) { 665 matrix_ = *model_>solver()>getMatrixByCol(); 666 } 667 delete [] used_; 668 int numberColumns = model>solver()>getNumCols(); 669 used_ = new int[numberColumns]; 670 memset(used_, 0, numberColumns*sizeof(int)); 671 671 } 672 672 // Default Constructor 673 CbcHeuristicNaive::CbcHeuristicNaive() 674 :CbcHeuristic()675 { 676 large_=1.0e6;673 CbcHeuristicNaive::CbcHeuristicNaive() 674 : CbcHeuristic() 675 { 676 large_ = 1.0e6; 677 677 } 678 678 … … 680 680 681 681 CbcHeuristicNaive::CbcHeuristicNaive(CbcModel & model) 682 :CbcHeuristic(model)683 { 684 large_=1.0e6;685 } 686 687 // Destructor 682 : CbcHeuristic(model) 683 { 684 large_ = 1.0e6; 685 } 686 687 // Destructor 688 688 CbcHeuristicNaive::~CbcHeuristicNaive () 689 689 { … … 694 694 CbcHeuristicNaive::clone() const 695 695 { 696 return new CbcHeuristicNaive(*this);696 return new CbcHeuristicNaive(*this); 697 697 } 698 698 // Create C++ lines to get to current state 699 void 700 CbcHeuristicNaive::generateCpp( FILE * fp) 701 { 702 CbcHeuristicNaive other;703 fprintf(fp,"0#include \"CbcHeuristicLocal.hpp\"\n");704 fprintf(fp,"3 CbcHeuristicNaive naive(*cbcModel);\n");705 CbcHeuristic::generateCpp(fp,"naive");706 if (large_!=other.large_)707 fprintf(fp,"3 naive.setLarge(%g);\n",large_);708 else709 fprintf(fp,"4 naive.setLarge(%g);\n",large_);710 fprintf(fp,"3 cbcModel>addHeuristic(&naive);\n");711 } 712 713 // Copy constructor 699 void 700 CbcHeuristicNaive::generateCpp( FILE * fp) 701 { 702 CbcHeuristicNaive other; 703 fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n"); 704 fprintf(fp, "3 CbcHeuristicNaive naive(*cbcModel);\n"); 705 CbcHeuristic::generateCpp(fp, "naive"); 706 if (large_ != other.large_) 707 fprintf(fp, "3 naive.setLarge(%g);\n", large_); 708 else 709 fprintf(fp, "4 naive.setLarge(%g);\n", large_); 710 fprintf(fp, "3 cbcModel>addHeuristic(&naive);\n"); 711 } 712 713 // Copy constructor 714 714 CbcHeuristicNaive::CbcHeuristicNaive(const CbcHeuristicNaive & rhs) 715 :716 CbcHeuristic(rhs),717 large_(rhs.large_)718 { 719 } 720 721 // Assignment operator 722 CbcHeuristicNaive & 723 CbcHeuristicNaive::operator=( const CbcHeuristicNaive & rhs)724 { 725 if (this!=&rhs) {726 CbcHeuristic::operator=(rhs);727 large_ = rhs.large_;728 }729 return *this;715 : 716 CbcHeuristic(rhs), 717 large_(rhs.large_) 718 { 719 } 720 721 // Assignment operator 722 CbcHeuristicNaive & 723 CbcHeuristicNaive::operator=( const CbcHeuristicNaive & rhs) 724 { 725 if (this != &rhs) { 726 CbcHeuristic::operator=(rhs); 727 large_ = rhs.large_; 728 } 729 return *this; 730 730 } 731 731 732 732 // Resets stuff if model changes 733 void 733 void 734 734 CbcHeuristicNaive::resetModel(CbcModel * model) 735 735 { 736 CbcHeuristic::resetModel(model);736 CbcHeuristic::resetModel(model); 737 737 } 738 738 int 739 739 CbcHeuristicNaive::solution(double & solutionValue, 740 double * betterSolution) 741 { 742 numCouldRun_++; 743 // See if to do 744 bool atRoot = model_>getNodeCount()==0; 745 int passNumber = model_>getCurrentPassNumber(); 746 if (!when()(when()==1&&model_>phase()!=1)!atRootpassNumber!=1) 747 return 0; // switched off 748 // Don't do if it was this heuristic which found solution! 749 if (this==model_>lastHeuristic()) 750 return 0; 751 numRuns_++; 752 double cutoff; 753 model_>solver()>getDblParam(OsiDualObjectiveLimit,cutoff); 754 double direction = model_>solver()>getObjSense(); 755 cutoff *= direction; 756 cutoff = CoinMin(cutoff,solutionValue); 757 OsiSolverInterface * solver = model_>continuousSolver(); 758 if (!solver) 759 solver = model_>solver(); 760 const double * colLower = solver>getColLower(); 761 const double * colUpper = solver>getColUpper(); 762 const double * objective = solver>getObjCoefficients(); 763 764 int numberColumns = model_>getNumCols(); 765 int numberIntegers = model_>numberIntegers(); 766 const int * integerVariable = model_>integerVariable(); 767 768 int i; 769 bool solutionFound=false; 770 CoinWarmStartBasis saveBasis; 771 CoinWarmStartBasis * basis = 772 dynamic_cast<CoinWarmStartBasis *>(solver>getWarmStart()) ; 773 if (basis) { 774 saveBasis = * basis; 775 delete basis; 776 } 777 // First just fix all integers as close to zero as possible 778 OsiSolverInterface * newSolver = cloneBut(7); // wassolver>clone(); 779 for (i=0;i<numberIntegers;i++) { 780 int iColumn=integerVariable[i]; 781 double lower = colLower[iColumn]; 782 double upper = colUpper[iColumn]; 783 double value; 784 if (lower>0.0) 785 value=lower; 786 else if (upper<0.0) 787 value=upper; 788 else 789 value=0.0; 790 newSolver>setColLower(iColumn,value); 791 newSolver>setColUpper(iColumn,value); 792 } 793 newSolver>initialSolve(); 794 if (newSolver>isProvenOptimal()) { 795 double solValue = newSolver>getObjValue()*direction ; 796 if (solValue<cutoff) { 797 // we have a solution 798 solutionFound=true; 799 solutionValue=solValue; 800 memcpy(betterSolution,newSolver>getColSolution(), 801 numberColumns*sizeof(double)); 802 printf("Naive fixing close to zero gave solution of %g\n",solutionValue); 803 cutoff = solValue  model_>getCutoffIncrement(); 804 } 805 } 806 // Now fix all integers as close to zero if zero or large cost 807 int nFix=0; 808 for (i=0;i<numberIntegers;i++) { 809 int iColumn=integerVariable[i]; 810 double lower = colLower[iColumn]; 811 double upper = colUpper[iColumn]; 812 double value; 813 if (fabs(objective[i])>0.0&&fabs(objective[i])<large_) { 814 nFix++; 815 if (lower>0.0) 816 value=lower; 817 else if (upper<0.0) 818 value=upper; 819 else 820 value=0.0; 821 newSolver>setColLower(iColumn,value); 822 newSolver>setColUpper(iColumn,value); 823 } else { 824 // set back to original 825 newSolver>setColLower(iColumn,lower); 826 newSolver>setColUpper(iColumn,upper); 827 } 828 } 829 const double * solution = solver>getColSolution(); 830 if (nFix) { 740 double * betterSolution) 741 { 742 numCouldRun_++; 743 // See if to do 744 bool atRoot = model_>getNodeCount() == 0; 745 int passNumber = model_>getCurrentPassNumber(); 746 if (!when()  (when() == 1 && model_>phase() != 1)  !atRoot  passNumber != 1) 747 return 0; // switched off 748 // Don't do if it was this heuristic which found solution! 749 if (this == model_>lastHeuristic()) 750 return 0; 751 numRuns_++; 752 double cutoff; 753 model_>solver()>getDblParam(OsiDualObjectiveLimit, cutoff); 754 double direction = model_>solver()>getObjSense(); 755 cutoff *= direction; 756 cutoff = CoinMin(cutoff, solutionValue); 757 OsiSolverInterface * solver = model_>continuousSolver(); 758 if (!solver) 759 solver = model_>solver(); 760 const double * colLower = solver>getColLower(); 761 const double * colUpper = solver>getColUpper(); 762 const double * objective = solver>getObjCoefficients(); 763 764 int numberColumns = model_>getNumCols(); 765 int numberIntegers = model_>numberIntegers(); 766 const int * integerVariable = model_>integerVariable(); 767 768 int i; 769 bool solutionFound = false; 770 CoinWarmStartBasis saveBasis; 771 CoinWarmStartBasis * basis = 772 dynamic_cast<CoinWarmStartBasis *>(solver>getWarmStart()) ; 773 if (basis) { 774 saveBasis = * basis; 775 delete basis; 776 } 777 // First just fix all integers as close to zero as possible 778 OsiSolverInterface * newSolver = cloneBut(7); // wassolver>clone(); 779 for (i = 0; i < numberIntegers; i++) { 780 int iColumn = integerVariable[i]; 781 double lower = colLower[iColumn]; 782 double upper = colUpper[iColumn]; 783 double value; 784 if (lower > 0.0) 785 value = lower; 786 else if (upper < 0.0) 787 value = upper; 788 else 789 value = 0.0; 790 newSolver>setColLower(iColumn, value); 791 newSolver>setColUpper(iColumn, value); 792 } 793 newSolver>initialSolve(); 794 if (newSolver>isProvenOptimal()) { 795 double solValue = newSolver>getObjValue() * direction ; 796 if (solValue < cutoff) { 797 // we have a solution 798 solutionFound = true; 799 solutionValue = solValue; 800 memcpy(betterSolution, newSolver>getColSolution(), 801 numberColumns*sizeof(double)); 802 printf("Naive fixing close to zero gave solution of %g\n", solutionValue); 803 cutoff = solValue  model_>getCutoffIncrement(); 804 } 805 } 806 // Now fix all integers as close to zero if zero or large cost 807 int nFix = 0; 808 for (i = 0; i < numberIntegers; i++) { 809 int iColumn = integerVariable[i]; 810 double lower = colLower[iColumn]; 811 double upper = colUpper[iColumn]; 812 double value; 813 if (fabs(objective[i]) > 0.0 && fabs(objective[i]) < large_) { 814 nFix++; 815 if (lower > 0.0) 816 value = lower; 817 else if (upper < 0.0) 818 value = upper; 819 else 820 value = 0.0; 821 newSolver>setColLower(iColumn, value); 822 newSolver>setColUpper(iColumn, value); 823 } else { 824 // set back to original 825 newSolver>setColLower(iColumn, lower); 826 newSolver>setColUpper(iColumn, upper); 827 } 828 } 829 const double * solution = solver>getColSolution(); 830 if (nFix) { 831 newSolver>setWarmStart(&saveBasis); 832 newSolver>setColSolution(solution); 833 newSolver>initialSolve(); 834 if (newSolver>isProvenOptimal()) { 835 double solValue = newSolver>getObjValue() * direction ; 836 if (solValue < cutoff) { 837 // try branch and bound 838 double * newSolution = new double [numberColumns]; 839 printf("%d fixed after fixing costs\n", nFix); 840 int returnCode = smallBranchAndBound(newSolver, 841 numberNodes_, newSolution, 842 solutionValue, 843 solutionValue, "CbcHeuristicNaive1"); 844 if (returnCode < 0) 845 returnCode = 0; // returned on size 846 if ((returnCode&2) != 0) { 847 // could add cut 848 returnCode &= ~2; 849 } 850 if (returnCode == 1) { 851 // solution 852 solutionFound = true; 853 memcpy(betterSolution, newSolution, 854 numberColumns*sizeof(double)); 855 printf("Naive fixing zeros gave solution of %g\n", solutionValue); 856 cutoff = solutionValue  model_>getCutoffIncrement(); 857 } 858 delete [] newSolution; 859 } 860 } 861 } 862 #if 1 863 newSolver>setObjSense(direction); // maximize 831 864 newSolver>setWarmStart(&saveBasis); 832 865 newSolver>setColSolution(solution); 866 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 867 double value = solution[iColumn]; 868 double lower = colLower[iColumn]; 869 double upper = colUpper[iColumn]; 870 double newLower; 871 double newUpper; 872 if (newSolver>isInteger(iColumn)) { 873 newLower = CoinMax(lower, floor(value)  2.0); 874 newUpper = CoinMin(upper, ceil(value) + 2.0); 875 } else { 876 newLower = CoinMax(lower, value  1.0e5); 877 newUpper = CoinMin(upper, value + 1.0e5); 878 } 879 newSolver>setColLower(iColumn, newLower); 880 newSolver>setColUpper(iColumn, newUpper); 881 } 833 882 newSolver>initialSolve(); 834 883 if (newSolver>isProvenOptimal()) { 835 double solValue = newSolver>getObjValue()*direction ; 836 if (solValue<cutoff) { 837 // try branch and bound 838 double * newSolution = new double [numberColumns]; 839 printf("%d fixed after fixing costs\n",nFix); 840 int returnCode = smallBranchAndBound(newSolver, 841 numberNodes_,newSolution, 842 solutionValue, 843 solutionValue,"CbcHeuristicNaive1"); 844 if (returnCode<0) 845 returnCode=0; // returned on size 846 if ((returnCode&2)!=0) { 847 // could add cut 848 returnCode &= ~2; 849 } 850 if (returnCode==1) { 851 // solution 852 solutionFound=true; 853 memcpy(betterSolution,newSolution, 854 numberColumns*sizeof(double)); 855 printf("Naive fixing zeros gave solution of %g\n",solutionValue); 856 cutoff = solutionValue  model_>getCutoffIncrement(); 857 } 858 delete [] newSolution; 859 } 860 } 861 } 862 #if 1 863 newSolver>setObjSense(direction); // maximize 864 newSolver>setWarmStart(&saveBasis); 865 newSolver>setColSolution(solution); 866 for (int iColumn=0;iColumn<numberColumns;iColumn++) { 867 double value = solution[iColumn]; 868 double lower = colLower[iColumn]; 869 double upper = colUpper[iColumn]; 870 double newLower; 871 double newUpper; 872 if (newSolver>isInteger(iColumn)) { 873 newLower = CoinMax(lower,floor(value)2.0); 874 newUpper = CoinMin(upper,ceil(value)+2.0); 875 } else { 876 newLower = CoinMax(lower,value1.0e5); 877 newUpper = CoinMin(upper,value+1.0e5); 878 } 879 newSolver>setColLower(iColumn,newLower); 880 newSolver>setColUpper(iColumn,newUpper); 881 } 882 newSolver>initialSolve(); 883 if (newSolver>isProvenOptimal()) { 884 double solValue = newSolver>getObjValue()*direction ; 885 if (solValue<cutoff) { 886 nFix=0; 887 newSolver>setObjSense(direction); // correct direction 888 //const double * thisSolution = newSolver>getColSolution(); 889 for (int iColumn=0;iColumn<numberColumns;iColumn++) { 890 double value = solution[iColumn]; 891 double lower = colLower[iColumn]; 892 double upper = colUpper[iColumn]; 893 double newLower=lower; 894 double newUpper=upper; 895 if (newSolver>isInteger(iColumn)) { 896 if (value<lower+1.0e6) { 897 nFix++; 898 newUpper=lower; 899 } else if (value>upper1.0e6) { 900 nFix++; 901 newLower=upper; 902 } else { 903 newLower = CoinMax(lower,floor(value)2.0); 904 newUpper = CoinMin(upper,ceil(value)+2.0); 905 } 906 } 907 newSolver>setColLower(iColumn,newLower); 908 newSolver>setColUpper(iColumn,newUpper); 909 } 910 // try branch and bound 911 double * newSolution = new double [numberColumns]; 912 printf("%d fixed after maximizing\n",nFix); 913 int returnCode = smallBranchAndBound(newSolver, 914 numberNodes_,newSolution, 915 solutionValue, 916 solutionValue,"CbcHeuristicNaive1"); 917 if (returnCode<0) 918 returnCode=0; // returned on size 919 if ((returnCode&2)!=0) { 920 // could add cut 921 returnCode &= ~2; 922 } 923 if (returnCode==1) { 924 // solution 925 solutionFound=true; 926 memcpy(betterSolution,newSolution, 927 numberColumns*sizeof(double)); 928 printf("Naive maximizing gave solution of %g\n",solutionValue); 929 cutoff = solutionValue  model_>getCutoffIncrement(); 930 } 931 delete [] newSolution; 932 } 933 } 884 double solValue = newSolver>getObjValue() * direction ; 885 if (solValue < cutoff) { 886 nFix = 0; 887 newSolver>setObjSense(direction); // correct direction 888 //const double * thisSolution = newSolver>getColSolution(); 889 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 890 double value = solution[iColumn]; 891 double lower = colLower[iColumn]; 892 double upper = colUpper[iColumn]; 893 double newLower = lower; 894 double newUpper = upper; 895 if (newSolver>isInteger(iColumn)) { 896 if (value < lower + 1.0e6) { 897 nFix++; 898 newUpper = lower; 899 } else if (value > upper  1.0e6) { 900 nFix++; 901 newLower = upper; 902 } else { 903 newLower = CoinMax(lower, floor(value)  2.0); 904 newUpper = CoinMin(upper, ceil(value) + 2.0); 905 } 906 } 907 newSolver>setColLower(iColumn, newLower); 908 newSolver>setColUpper(iColumn, newUpper); 909 } 910 // try branch and bound 911 double * newSolution = new double [numberColumns]; 912 printf("%d fixed after maximizing\n", nFix); 913 int returnCode = smallBranchAndBound(newSolver, 914 numberNodes_, newSolution, 915 solutionValue, 916 solutionValue, "CbcHeuristicNaive1"); 917 if (returnCode < 0) 918 returnCode = 0; // returned on size 919 if ((returnCode&2) != 0) { 920 // could add cut 921 returnCode &= ~2; 922 } 923 if (returnCode == 1) { 924 // solution 925 solutionFound = true; 926 memcpy(betterSolution, newSolution, 927 numberColumns*sizeof(double)); 928 printf("Naive maximizing gave solution of %g\n", solutionValue); 929 cutoff = solutionValue  model_>getCutoffIncrement(); 930 } 931 delete [] newSolution; 932 } 933 } 934 934 #endif 935 delete newSolver;936 return solutionFound ? 1 : 0;935 delete newSolver; 936 return solutionFound ? 1 : 0; 937 937 } 938 938 // update model 939 939 void CbcHeuristicNaive::setModel(CbcModel * model) 940 940 { 941 model_ = model;941 model_ = model; 942 942 } 943 943 // Default Constructor 944 CbcHeuristicCrossover::CbcHeuristicCrossover() 945 :CbcHeuristic(),946 numberSolutions_(0),947 useNumber_(3)948 { 949 setWhen(1);944 CbcHeuristicCrossover::CbcHeuristicCrossover() 945 : CbcHeuristic(), 946 numberSolutions_(0), 947 useNumber_(3) 948 { 949 setWhen(1); 950 950 } 951 951 … … 953 953 954 954 CbcHeuristicCrossover::CbcHeuristicCrossover(CbcModel & model) 955 :CbcHeuristic(model),956 numberSolutions_(0),957 useNumber_(3)958 { 959 setWhen(1);960 for (int i=0;i<10;i++)961 random_[i]=model.randomNumberGenerator()>randomDouble();962 } 963 964 // Destructor 955 : CbcHeuristic(model), 956 numberSolutions_(0), 957 useNumber_(3) 958 { 959 setWhen(1); 960 for (int i = 0; i < 10; i++) 961 random_[i] = model.randomNumberGenerator()>randomDouble(); 962 } 963 964 // Destructor 965 965 CbcHeuristicCrossover::~CbcHeuristicCrossover () 966 966 { … … 971 971 CbcHeuristicCrossover::clone() const 972 972 { 973 return new CbcHeuristicCrossover(*this);973 return new CbcHeuristicCrossover(*this); 974 974 } 975 975 // Create C++ lines to get to current state 976 void 977 CbcHeuristicCrossover::generateCpp( FILE * fp) 978 { 979 CbcHeuristicCrossover other;980 fprintf(fp,"0#include \"CbcHeuristicLocal.hpp\"\n");981 fprintf(fp,"3 CbcHeuristicCrossover crossover(*cbcModel);\n");982 CbcHeuristic::generateCpp(fp,"crossover");983 if (useNumber_!=other.useNumber_)984 fprintf(fp,"3 crossover.setNumberSolutions(%d);\n",useNumber_);985 else986 fprintf(fp,"4 crossover.setNumberSolutions(%d);\n",useNumber_);987 fprintf(fp,"3 cbcModel>addHeuristic(&crossover);\n");988 } 989 990 // Copy constructor 976 void 977 CbcHeuristicCrossover::generateCpp( FILE * fp) 978 { 979 CbcHeuristicCrossover other; 980 fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n"); 981 fprintf(fp, "3 CbcHeuristicCrossover crossover(*cbcModel);\n"); 982 CbcHeuristic::generateCpp(fp, "crossover"); 983 if (useNumber_ != other.useNumber_) 984 fprintf(fp, "3 crossover.setNumberSolutions(%d);\n", useNumber_); 985 else 986 fprintf(fp, "4 crossover.setNumberSolutions(%d);\n", useNumber_); 987 fprintf(fp, "3 cbcModel>addHeuristic(&crossover);\n"); 988 } 989 990 // Copy constructor 991 991 CbcHeuristicCrossover::CbcHeuristicCrossover(const CbcHeuristicCrossover & rhs) 992 :993 CbcHeuristic(rhs),994 attempts_(rhs.attempts_),995 numberSolutions_(rhs.numberSolutions_),996 useNumber_(rhs.useNumber_)997 { 998 memcpy(random_,rhs.random_,10*sizeof(double));999 } 1000 1001 // Assignment operator 1002 CbcHeuristicCrossover & 1003 CbcHeuristicCrossover::operator=( const CbcHeuristicCrossover & rhs)1004 { 1005 if (this!=&rhs) {1006 CbcHeuristic::operator=(rhs);1007 useNumber_ = rhs.useNumber_;1008 attempts_ = rhs.attempts_;1009 numberSolutions_ = rhs.numberSolutions_;1010 memcpy(random_,rhs.random_,10*sizeof(double));1011 }1012 return *this;992 : 993 CbcHeuristic(rhs), 994 attempts_(rhs.attempts_), 995 numberSolutions_(rhs.numberSolutions_), 996 useNumber_(rhs.useNumber_) 997 { 998 memcpy(random_, rhs.random_, 10*sizeof(double)); 999 } 1000 1001 // Assignment operator 1002 CbcHeuristicCrossover & 1003 CbcHeuristicCrossover::operator=( const CbcHeuristicCrossover & rhs) 1004 { 1005 if (this != &rhs) { 1006 CbcHeuristic::operator=(rhs); 1007 useNumber_ = rhs.useNumber_; 1008 attempts_ = rhs.attempts_; 1009 numberSolutions_ = rhs.numberSolutions_; 1010 memcpy(random_, rhs.random_, 10*sizeof(double)); 1011 } 1012 return *this; 1013 1013 } 1014 1014 1015 1015 // Resets stuff if model changes 1016 void 1016 void 1017 1017 CbcHeuristicCrossover::resetModel(CbcModel * model) 1018 1018 { 1019 CbcHeuristic::resetModel(model);1019 CbcHeuristic::resetModel(model); 1020 1020 } 1021 1021 int 1022 1022 CbcHeuristicCrossover::solution(double & solutionValue, 1023 1024 { 1025 if (when_==0)1026 return 0;1027 numCouldRun_++;1028 bool useBest=(numberSolutions_!=model_>getSolutionCount());1029 if (!useBest&&(when_%10)==1)1030 return 0;1031 numberSolutions_=model_>getSolutionCount();1032 OsiSolverInterface * continuousSolver = model_>continuousSolver();1033 int useNumber =CoinMin(model_>numberSavedSolutions(),useNumber_);1034 if (useNumber<2!continuousSolver)1035 return 0;1036 // Fix later1037 if (!useBest)1038 abort();1039 numRuns_++;1040 double cutoff;1041 model_>solver()>getDblParam(OsiDualObjectiveLimit,cutoff);1042 double direction = model_>solver()>getObjSense();1043 cutoff *= direction;1044 cutoff = CoinMin(cutoff,solutionValue);1045 OsiSolverInterface * solver = cloneBut(2);1046 // But reset bounds1047 solver>setColLower(continuousSolver>getColLower());1048 solver>setColUpper(continuousSolver>getColUpper());1049 int numberColumns = solver>getNumCols();1050 // Fixed1051 double * fixed =new double [numberColumns];1052 for (int i=0;i<numberColumns;i++)1053 fixed[i]=COIN_DBL_MAX;1054 int whichSolution[10];1055 for (int i=0;i<useNumber;i++)1056 whichSolution[i]=i;1057 for (int i=0;i<useNumber;i++) {1058 int k = whichSolution[i];1059 const double * solution = model_>savedSolution(k);1060 for (int j=0;j<numberColumns;j++) {1061 if (solver>isInteger(j)) {1062 if (fixed[j]==COIN_DBL_MAX) 1063 fixed[j]=floor(solution[j]+0.5);1064 else if (fabs(fixed[j]solution[j])>1.0e7)1065 fixed[j]=COIN_DBL_MAX;1066 }1067 }1068 }1069 const double * colLower = solver>getColLower();1070 for (int i=0;i<numberColumns;i++) {1071 if (solver>isInteger(i)) {1072 double value=fixed[i];1073 if (value!=COIN_DBL_MAX) {1074 if (when_<10) {1075 solver>setColLower(i,value);1076 solver>setColUpper(i,value);1077 } else if (value==colLower[i]) {1078 solver>setColUpper(i,value);1079 1080 }1081 }1082 }1083 int returnCode = smallBranchAndBound(solver,numberNodes_,betterSolution,1084 1085 solutionValue,"CbcHeuristicCrossover");1086 if (returnCode<0)1087 returnCode=0; // returned on size1088 if ((returnCode&2)!=0) {1089 // could add cut1090 returnCode &= ~2;1091 }1092 1093 delete solver;1094 return returnCode;1023 double * betterSolution) 1024 { 1025 if (when_ == 0) 1026 return 0; 1027 numCouldRun_++; 1028 bool useBest = (numberSolutions_ != model_>getSolutionCount()); 1029 if (!useBest && (when_ % 10) == 1) 1030 return 0; 1031 numberSolutions_ = model_>getSolutionCount(); 1032 OsiSolverInterface * continuousSolver = model_>continuousSolver(); 1033 int useNumber = CoinMin(model_>numberSavedSolutions(), useNumber_); 1034 if (useNumber < 2  !continuousSolver) 1035 return 0; 1036 // Fix later 1037 if (!useBest) 1038 abort(); 1039 numRuns_++; 1040 double cutoff; 1041 model_>solver()>getDblParam(OsiDualObjectiveLimit, cutoff); 1042 double direction = model_>solver()>getObjSense(); 1043 cutoff *= direction; 1044 cutoff = CoinMin(cutoff, solutionValue); 1045 OsiSolverInterface * solver = cloneBut(2); 1046 // But reset bounds 1047 solver>setColLower(continuousSolver>getColLower()); 1048 solver>setColUpper(continuousSolver>getColUpper()); 1049 int numberColumns = solver>getNumCols(); 1050 // Fixed 1051 double * fixed = new double [numberColumns]; 1052 for (int i = 0; i < numberColumns; i++) 1053 fixed[i] = COIN_DBL_MAX; 1054 int whichSolution[10]; 1055 for (int i = 0; i < useNumber; i++) 1056 whichSolution[i] = i; 1057 for (int i = 0; i < useNumber; i++) { 1058 int k = whichSolution[i]; 1059 const double * solution = model_>savedSolution(k); 1060 for (int j = 0; j < numberColumns; j++) { 1061 if (solver>isInteger(j)) { 1062 if (fixed[j] == COIN_DBL_MAX) 1063 fixed[j] = floor(solution[j] + 0.5); 1064 else if (fabs(fixed[j]  solution[j]) > 1.0e7) 1065 fixed[j] = COIN_DBL_MAX; 1066 } 1067 } 1068 } 1069 const double * colLower = solver>getColLower(); 1070 for (int i = 0; i < numberColumns; i++) { 1071 if (solver>isInteger(i)) { 1072 double value = fixed[i]; 1073 if (value != COIN_DBL_MAX) { 1074 if (when_ < 10) { 1075 solver>setColLower(i, value); 1076 solver>setColUpper(i, value); 1077 } else if (value == colLower[i]) { 1078 solver>setColUpper(i, value); 1079 } 1080 } 1081 } 1082 } 1083 int returnCode = smallBranchAndBound(solver, numberNodes_, betterSolution, 1084 solutionValue, 1085 solutionValue, "CbcHeuristicCrossover"); 1086 if (returnCode < 0) 1087 returnCode = 0; // returned on size 1088 if ((returnCode&2) != 0) { 1089 // could add cut 1090 returnCode &= ~2; 1091 } 1092 1093 delete solver; 1094 return returnCode; 1095 1095 } 1096 1096 // update model 1097 1097 void CbcHeuristicCrossover::setModel(CbcModel * model) 1098 1098 { 1099 model_ = model;1100 if (model) {1101 for (int i=0;i<10;i++)1102 random_[i]=model>randomNumberGenerator()>randomDouble();1103 }1104 } 1105 1106 1099 model_ = model; 1100 if (model) { 1101 for (int i = 0; i < 10; i++) 1102 random_[i] = model>randomNumberGenerator()>randomDouble(); 1103 } 1104 } 1105 1106
Note: See TracChangeset
for help on using the changeset viewer.