Changeset 1582 for trunk/Cbc/src/CbcHeuristicRENS.cpp
- Timestamp:
- Jan 7, 2011 12:16:00 PM (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Cbc/src/CbcHeuristicRENS.cpp
r1573 r1582 90 90 return 0; 91 91 numberTries_++; 92 double saveFractionSmall=fractionSmall_; 92 93 OsiSolverInterface * solver = model_->solver(); 93 94 … … 95 96 const int * integerVariable = model_->integerVariable(); 96 97 97 const double * currentSolution = solver->getColSolution();98 98 OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone(); 99 newSolver->resolve(); 99 const double * currentSolution = newSolver->getColSolution(); 100 int type = rensType_&15; 101 if (type<11) 102 newSolver->resolve(); 100 103 double direction = newSolver->getObjSense(); 101 104 double cutoff=model_->getCutoff(); 102 //newSolver->getDblParam(OsiDualObjectiveLimit, cutoff);105 newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100); 103 106 //cutoff *= direction; 104 107 double gap = cutoff - newSolver->getObjValue() * direction ; 105 108 double tolerance; 106 109 newSolver->getDblParam(OsiDualTolerance, tolerance) ; 107 if ( gap > 0.0 || !newSolver->isProvenOptimal()) {110 if ((gap > 0.0 || !newSolver->isProvenOptimal())&&type<11) { 108 111 gap += 100.0 * tolerance; 109 112 int nFix = newSolver->reducedCostFix(gap); … … 115 118 << CoinMessageEol; 116 119 } 117 } else {120 } else if (type<11) { 118 121 return 0; // finished? 119 122 } 120 123 int numberColumns = solver->getNumCols(); 121 124 double * dj = CoinCopyOfArray(solver->getReducedCost(),numberColumns); 122 int type = rensType_&15;123 125 double djTolerance = (type!=1) ? -1.0e30 : 1.0e-4; 124 126 const double * colLower = newSolver->getColLower(); 125 127 const double * colUpper = newSolver->getColUpper(); 128 double * contribution = NULL; 126 129 int numberFixed = 0; 127 130 if (type==3) { … … 143 146 delete basis; 144 147 } 145 } else if (type>=5&&type<=1 0) {148 } else if (type>=5&&type<=11) { 146 149 /* 5 fix sets at one 147 150 6 fix on dj but leave unfixed SOS slacks … … 149 152 8 fix all at zero but leave unfixed SOS slacks 150 153 9 as 8 but only fix all at zero if just one in set nonzero 151 10 as 7 but pi other way 154 10 fix all "stable" ones 155 11 layered approach 152 156 */ 153 157 // SOS type fixing … … 203 207 } 204 208 } 209 // Just leave one slack in each set 210 { 211 const double * objective = newSolver->getObjCoefficients(); 212 int * best = new int [numberRows]; 213 double * cheapest = new double[numberRows]; 214 for (int i=0;i<numberRows;i++) { 215 best[i]=-1; 216 cheapest[i]=COIN_DBL_MAX; 217 } 218 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 219 if (colUpper[iColumn]>colLower[iColumn]) { 220 if (columnLength[iColumn]==1) { 221 CoinBigIndex j = columnStart[iColumn]; 222 int iRow = row[j]; 223 if (bestDj[iRow]<1.0e30) { 224 double obj = direction*objective[iColumn]; 225 if (obj<cheapest[iRow]) { 226 cheapest[iRow]=obj; 227 best[iRow]=iColumn; 228 } 229 } 230 } 231 } 232 } 233 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 234 if (colUpper[iColumn]>colLower[iColumn]) { 235 if (columnLength[iColumn]==1) { 236 CoinBigIndex j = columnStart[iColumn]; 237 int iRow = row[j]; 238 if (bestDj[iRow]<1.0e30) { 239 if (best[iRow]!=-1&&iColumn!=best[iRow]) { 240 newSolver->setColUpper(iColumn,0.0); 241 } 242 } 243 } 244 } 245 } 246 delete [] best; 247 delete [] cheapest; 248 } 205 249 int nSOS=0; 206 250 double * sort = new double [numberRows]; 207 251 const double * pi = newSolver->getRowPrice(); 252 if (type==11) { 253 contribution = new double [numberRows]; 254 for (int i=0;i<numberRows;i++) { 255 if (bestDj[i]<1.0e30) 256 contribution[i]=0.0; 257 else 258 contribution[i]=-1.0; 259 } 260 } 208 261 for (int i=0;i<numberRows;i++) { 209 262 if (bestDj[i]<1.0e30) { … … 217 270 } 218 271 if (10*nSOS>8*numberRows) { 219 std::sort(sort,sort+nSOS); 220 int last = static_cast<int>(nSOS*0.9*fractionSmall_); 221 double tolerance = sort[last]; 222 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 223 if (colUpper[iColumn]>colLower[iColumn]) { 224 CoinBigIndex j; 225 if (currentSolution[iColumn]<=1.0e-6|| 226 currentSolution[iColumn]>=0.999999) { 227 if (fixSets) { 228 for (j = columnStart[iColumn]; 229 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 230 int iRow = row[j]; 231 double useDj; 232 if (type==5) 233 useDj = bestDj[iRow]; 234 else if (type==7) 235 useDj= -fabs(pi[iRow]); 236 else 237 useDj= fabs(pi[iRow]); 238 if (bestDj[iRow]<1.0e30&&useDj>=tolerance) { 239 numberFixed++; 240 if (currentSolution[iColumn]<=1.0e-6) 241 newSolver->setColUpper(iColumn,0.0); 242 else if (currentSolution[iColumn]>=0.999999) 243 newSolver->setColLower(iColumn,1.0); 272 if (type<10) { 273 std::sort(sort,sort+nSOS); 274 int last = static_cast<int>(nSOS*0.9*fractionSmall_); 275 double tolerance = sort[last]; 276 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 277 if (colUpper[iColumn]>colLower[iColumn]) { 278 CoinBigIndex j; 279 if (currentSolution[iColumn]<=1.0e-6|| 280 currentSolution[iColumn]>=0.999999) { 281 if (fixSets) { 282 for (j = columnStart[iColumn]; 283 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 284 int iRow = row[j]; 285 double useDj; 286 if (type==5) 287 useDj = bestDj[iRow]; 288 else if (type==7) 289 useDj= -fabs(pi[iRow]); 290 else 291 useDj= fabs(pi[iRow]); 292 if (bestDj[iRow]<1.0e30&&useDj>=tolerance) { 293 numberFixed++; 294 if (currentSolution[iColumn]<=1.0e-6) 295 newSolver->setColUpper(iColumn,0.0); 296 else if (currentSolution[iColumn]>=0.999999) 297 newSolver->setColLower(iColumn,1.0); 298 } 299 } 300 } else if (columnLength[iColumn]==1) { 301 // leave more slacks 302 int iRow = row[columnStart[iColumn]]; 303 if (bestDj[iRow]<1.0e30) { 304 // fake dj 305 dj[iColumn] *= 0.000001; 306 } 307 } else if (type==8||type==9) { 308 if (currentSolution[iColumn]<=1.0e-6) { 309 if (type==8) { 310 dj[iColumn] *= 1.0e6; 311 } else { 312 bool fix=false; 313 for (j = columnStart[iColumn]; 314 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 315 int iRow = row[j]; 316 if (bestDj[iRow]<1.0e25) { 317 fix=true; 318 break; 319 } 320 } 321 if (fix) { 322 dj[iColumn] *= 1.0e6; 323 } 324 } 325 } else { 326 dj[iColumn] *= 0.000001; 244 327 } 245 328 } 246 } else if (columnLength[iColumn]==1) {247 // leave more slacks248 int iRow = row[columnStart[iColumn]];249 if (bestDj[iRow]<1.0e30) {250 // fake dj251 dj[iColumn] *= 0.000001;252 }253 } else if (type==8||type==9) {254 if (currentSolution[iColumn]<=1.0e-6) {255 if (type==8) {256 dj[iColumn] *= 1.0e6;257 } else {258 bool fix=false;259 for (j = columnStart[iColumn];260 j < columnStart[iColumn] + columnLength[iColumn]; j++) {261 int iRow = row[j];262 if (bestDj[iRow]<1.0e25) {263 fix=true;264 break;265 }266 }267 if (fix) {268 dj[iColumn] *= 1.0e6;269 }270 }271 } else {272 dj[iColumn] *= 0.000001;273 }274 329 } 275 330 } 276 331 } 277 } 278 if (fixSets) 332 if (fixSets) 333 djTolerance = 1.0e30; 334 } else if (type==10) { 335 double * saveUpper = CoinCopyOfArray(newSolver->getRowUpper(),numberRows); 336 char * mark = new char [numberColumns]; 337 char * nonzero = new char [numberColumns]; 338 double factor=CoinMax(1.000001,fractionSmall_); 339 fractionSmall_ = 0.5; 340 // loosen up 341 for (int i=0;i<numberRows;i++) { 342 if (bestDj[i]>=1.0e30) { 343 newSolver->setRowUpper(i,factor*saveUpper[i]); 344 } 345 } 346 newSolver->resolve(); 347 const double * solution = newSolver->getColSolution(); 348 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 349 mark[iColumn]=0; 350 nonzero[iColumn]=0; 351 if (colUpper[iColumn]>colLower[iColumn]&& 352 solution[iColumn]>0.9999) 353 mark[iColumn]=1; 354 else if (solution[iColumn]>0.00001) 355 nonzero[iColumn]=1; 356 } 357 // slightly small 358 for (int i=0;i<numberRows;i++) { 359 if (bestDj[i]>=1.0e30) { 360 newSolver->setRowUpper(i,saveUpper[i]*0.9999); 361 } 362 } 363 newSolver->resolve(); 364 int nCheck=2; 365 if (newSolver->isProvenOptimal()) { 366 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 367 if (colUpper[iColumn]>colLower[iColumn]&& 368 solution[iColumn]>0.9999) 369 mark[iColumn]++; 370 else if (solution[iColumn]>0.00001) 371 nonzero[iColumn]=1; 372 } 373 } else { 374 nCheck=1; 375 } 376 // correct values 377 for (int i=0;i<numberRows;i++) { 378 if (bestDj[i]>=1.0e30) { 379 newSolver->setRowUpper(i,saveUpper[i]); 380 } 381 } 382 newSolver->resolve(); 383 int nFixed=0; 384 int nFixedToZero=0; 385 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 386 if (colUpper[iColumn]>colLower[iColumn]) { 387 if (solution[iColumn]>0.9999&&mark[iColumn]==nCheck) { 388 newSolver->setColLower(iColumn,1.0); 389 nFixed++; 390 } else if (!mark[iColumn]&&!nonzero[iColumn]&& 391 columnLength[iColumn]>1&&solution[iColumn]<0.00001) { 392 newSolver->setColUpper(iColumn,0.0); 393 nFixedToZero++; 394 } 395 } 396 } 397 char line[100]; 398 sprintf(line,"Heuristic %s fixed %d to one (and %d to zero)", 399 heuristicName(), 400 nFixed,nFixedToZero); 401 model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) 402 << line 403 << CoinMessageEol; 404 delete [] mark; 405 delete []nonzero; 406 delete [] saveUpper; 407 numberFixed=numberColumns; 279 408 djTolerance = 1.0e30; 409 } 280 410 } 281 411 delete basis; 282 412 delete [] sort; 283 413 delete [] bestDj; 414 if (10*nSOS<=8*numberRows) { 415 // give up 416 delete [] contribution; 417 delete newSolver; 418 return 0; 419 } 284 420 } 285 421 } 286 422 // Do dj to get right number 287 if (type==4||type==6|| type>7) {423 if (type==4||type==6||(type>7&&type<10)) { 288 424 double * sort = new double [numberColumns]; 289 425 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { … … 297 433 djTolerance = CoinMax(sort[last],1.0e-5); 298 434 delete [] sort; 435 } else if (type==11) { 436 // Do layered in a different way 437 int numberRows = solver->getNumRows(); 438 // Column copy 439 const CoinPackedMatrix * matrix = newSolver->getMatrixByCol(); 440 const double * element = matrix->getElements(); 441 const int * row = matrix->getIndices(); 442 const CoinBigIndex * columnStart = matrix->getVectorStarts(); 443 const int * columnLength = matrix->getVectorLengths(); 444 int * whichRow = new int[numberRows]; 445 int * whichSet = new int [numberColumns]; 446 int nSOS=0; 447 for (int i=0;i<numberRows;i++) { 448 whichRow[i]=0; 449 if (!contribution[i]) 450 nSOS++; 451 } 452 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 453 whichSet[iColumn]=-2; 454 if (colUpper[iColumn]>colLower[iColumn]) { 455 CoinBigIndex j; 456 double sum=0.0; 457 int iSOS=-1; 458 int n=0; 459 for (j = columnStart[iColumn]; 460 j < columnStart[iColumn] + columnLength[iColumn]; j++) { 461 int iRow = row[j]; 462 if (contribution[iRow]>=0.0) { 463 iSOS=iRow; 464 n++; 465 } else { 466 sum += fabs(element[j]); 467 } 468 } 469 if (n>1) 470 printf("Too many SOS entries (%d) for column %d\n", 471 n,iColumn); 472 if (sum) { 473 assert (iSOS>=0); 474 contribution[iSOS] += sum; 475 whichRow[iSOS]++; 476 whichSet[iColumn]=iSOS; 477 } else { 478 whichSet[iColumn]=iSOS+numberRows; 479 } 480 } 481 } 482 int * chunk = new int [numberRows]; 483 for (int i=0;i<numberRows;i++) { 484 chunk[i]=-1; 485 if (whichRow[i]) { 486 contribution[i]= - contribution[i]/static_cast<double>(whichRow[i]); 487 } else { 488 contribution[i] = COIN_DBL_MAX; 489 } 490 whichRow[i]=i; 491 } 492 newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100); 493 double * saveLower = CoinCopyOfArray(colLower,numberColumns); 494 double * saveUpper = CoinCopyOfArray(colUpper,numberColumns); 495 CoinSort_2(contribution,contribution+numberRows,whichRow); 496 // Set do nothing solution 497 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 498 if(whichSet[iColumn]>=numberRows) 499 newSolver->setColLower(iColumn,1.0); 500 } 501 newSolver->resolve(); 502 int nChunk = (nSOS+9)/10; 503 int nPass=0; 504 int inChunk=0; 505 for (int i=0;i<nSOS;i++) { 506 chunk[whichRow[i]]=nPass; 507 inChunk++; 508 if (inChunk==nChunk) { 509 inChunk=0; 510 // last two together 511 if (i+nChunk<nSOS) 512 nPass++; 513 } 514 } 515 // adjust 516 nPass++; 517 for (int iPass=0;iPass<nPass;iPass++) { 518 // fix last chunk and unfix this chunk 519 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 520 int iSOS = whichSet[iColumn]; 521 if (iSOS>=0) { 522 if (iSOS>=numberRows) 523 iSOS-=numberRows; 524 if (chunk[iSOS]==iPass-1&&betterSolution[iColumn]>0.9999) { 525 newSolver->setColLower(iColumn,1.0); 526 } else if (chunk[iSOS]==iPass) { 527 newSolver->setColLower(iColumn,saveLower[iColumn]); 528 newSolver->setColUpper(iColumn,saveUpper[iColumn]); 529 } 530 } 531 } 532 // solve 533 returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue, 534 model_->getCutoff(), "CbcHeuristicRENS"); 535 if (returnCode < 0) { 536 returnCode = 0; // returned on size 537 break; 538 } else if ((returnCode&1)==0) { 539 // no good 540 break; 541 } 542 } 543 if ((returnCode&2) != 0) { 544 // could add cut 545 returnCode &= ~2; 546 } 547 delete [] chunk; 548 delete [] saveLower; 549 delete [] saveUpper; 550 delete [] whichRow; 551 delete [] whichSet; 552 delete [] contribution; 553 delete newSolver; 554 return returnCode; 299 555 } 300 556 … … 490 746 } 491 747 } 492 748 //delete [] whichRow; 749 //delete [] contribution; 493 750 delete newSolver; 751 fractionSmall_ = saveFractionSmall; 494 752 return returnCode; 495 753 }
Note: See TracChangeset
for help on using the changeset viewer.