Changeset 463 for branches/devel/Cbc/src/CbcHeuristicFPump.cpp
- Timestamp:
- Oct 20, 2006 3:59:21 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/devel/Cbc/src/CbcHeuristicFPump.cpp
r444 r463 27 27 absoluteIncrement_(0.0), 28 28 relativeIncrement_(0.0), 29 initialWeight_(0.0), 30 weightFactor_(0.1), 29 31 maximumPasses_(100), 30 32 maximumRetries_(1), … … 45 47 absoluteIncrement_(0.0), 46 48 relativeIncrement_(0.0), 49 initialWeight_(0.0), 50 weightFactor_(0.1), 47 51 maximumPasses_(100), 48 52 maximumRetries_(1), … … 104 108 fprintf(fp,"4 heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_); 105 109 fprintf(fp,"3 cbcModel->addHeuristic(&heuristicFPump);\n"); 110 if (initialWeight_!=other.initialWeight_) 111 fprintf(fp,"3 heuristicFPump.setInitialWeight(%g);\n",initialWeight_); 112 else 113 fprintf(fp,"4 heuristicFPump.setInitialWeight(%g);\n",initialWeight_); 114 if (weightFactor_!=other.weightFactor_) 115 fprintf(fp,"3 heuristicFPump.setWeightFactor(%g);\n",weightFactor_); 116 else 117 fprintf(fp,"4 heuristicFPump.setWeightFactor(%g);\n",weightFactor_); 106 118 } 107 119 … … 116 128 absoluteIncrement_(rhs.absoluteIncrement_), 117 129 relativeIncrement_(rhs.relativeIncrement_), 130 initialWeight_(rhs.initialWeight_), 131 weightFactor_(rhs.weightFactor_), 118 132 maximumPasses_(rhs.maximumPasses_), 119 133 maximumRetries_(rhs.maximumRetries_), … … 154 168 cutoff *= direction; 155 169 cutoff = CoinMin(cutoff,solutionValue); 156 // space for rounded solution170 // check plausible and space for rounded solution 157 171 int numberColumns = model_->getNumCols(); 172 int numberIntegers = model_->numberIntegers(); 173 const int * integerVariableOrig = model_->integerVariable(); 174 175 // 1. initially check 0-1 176 int i,j; 177 int general=0; 178 int * integerVariable = new int[numberIntegers]; 179 const double * lower = model_->solver()->getColLower(); 180 const double * upper = model_->solver()->getColUpper(); 181 bool doGeneral = (accumulate_&4)!=0; 182 j=0; 183 for (i=0;i<numberIntegers;i++) { 184 int iColumn = integerVariableOrig[i]; 185 #ifndef NDEBUG 186 const OsiObject * object = model_->object(i); 187 const CbcSimpleInteger * integerObject = 188 dynamic_cast<const CbcSimpleInteger *> (object); 189 const OsiSimpleInteger * integerObject2 = 190 dynamic_cast<const OsiSimpleInteger *> (object); 191 assert(integerObject||integerObject2); 192 #endif 193 if (upper[iColumn]-lower[iColumn]>1.000001) { 194 general++; 195 if (doGeneral) 196 integerVariable[j++]=iColumn; 197 } else { 198 integerVariable[j++]=iColumn; 199 } 200 } 201 if (general*3>numberIntegers&&!doGeneral) { 202 delete [] integerVariable; 203 return 0; 204 } 205 int numberIntegersOrig = numberIntegers; 206 numberIntegers = j; 207 if (!doGeneral) 208 general=0; 158 209 double * newSolution = new double [numberColumns]; 159 210 double newSolutionValue=COIN_DBL_MAX; … … 219 270 solver->getDblParam(OsiPrimalTolerance,primalTolerance); 220 271 221 int numberIntegers = model_->numberIntegers();222 const int * integerVariable = model_->integerVariable();223 224 // 1. initially check 0-1225 int i,j;226 int general=0;227 for (i=0;i<numberIntegers;i++) {228 int iColumn = integerVariable[i];229 #ifndef NDEBUG230 const OsiObject * object = model_->object(i);231 const CbcSimpleInteger * integerObject =232 dynamic_cast<const CbcSimpleInteger *> (object);233 const OsiSimpleInteger * integerObject2 =234 dynamic_cast<const OsiSimpleInteger *> (object);235 assert(integerObject||integerObject2);236 #endif237 if (upper[iColumn]-lower[iColumn]>1.000001) {238 general++;239 break;240 }241 }242 if (general*3>numberIntegers) {243 delete solver;244 return 0;245 }246 272 247 273 // 2 space for last rounded solutions … … 271 297 double saveOffset; 272 298 solver->getDblParam(OsiObjOffset,saveOffset); 273 299 // Get amount for original objective 300 double scaleFactor = 0.0; 301 for (i=0;i<numberColumns;i++) 302 scaleFactor += saveObjective[i]*saveObjective[i]; 303 if (scaleFactor) 304 scaleFactor = (initialWeight_*sqrt((double) numberIntegers))/sqrt(scaleFactor); 274 305 // 5. MAIN WHILE LOOP 275 306 bool newLineNeeded=false; … … 291 322 memcpy(newSolution,solution,numberColumns*sizeof(double)); 292 323 int flip; 293 returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip); 324 returnCode = rounds(newSolution,saveObjective,numberIntegers,integerVariable, 325 roundExpensive_,downValue_,&flip); 294 326 if (returnCode) { 295 327 // SOLUTION IS INTEGER … … 312 344 if (general) { 313 345 int numberLeft=0; 314 for (i=0;i<numberIntegers ;i++) {315 int iColumn = integerVariable [i];346 for (i=0;i<numberIntegersOrig;i++) { 347 int iColumn = integerVariableOrig[i]; 316 348 double value = floor(newSolution[iColumn]+0.5); 317 if(solver->isBinary(iColumn) ) {349 if(solver->isBinary(iColumn)||general) { 318 350 solver->setColLower(iColumn,value); 319 351 solver->setColUpper(iColumn,value); … … 358 390 for (i = 0; i <numberIntegers; i++) { 359 391 int iColumn = integerVariable[i]; 360 if(!solver->isBinary(iColumn))361 continue;362 392 if (newSolution[iColumn]!=b[iColumn]) { 363 393 matched=false; … … 374 404 for (i=0;i<numberIntegers;i++) { 375 405 int iColumn = integerVariable[i]; 376 if(!solver->isBinary(iColumn))377 continue;378 406 double value = max(0.0,CoinDrand48()-0.3); 379 407 double difference = fabs(solution[iColumn]-newSolution[iColumn]); 380 408 if (difference+value>0.5) { 381 if (newSolution[iColumn]<lower[iColumn]+primalTolerance) newSolution[iColumn] += 1.0; 382 else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) newSolution[iColumn] -= 1.0; 383 else abort(); 409 if (newSolution[iColumn]<lower[iColumn]+primalTolerance) { 410 newSolution[iColumn] += 1.0; 411 } else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) { 412 newSolution[iColumn] -= 1.0; 413 } else { 414 // general integer 415 if (difference+value>0.75) 416 newSolution[iColumn] += 1.0; 417 else 418 newSolution[iColumn] -= 1.0; 419 } 384 420 } 385 421 } … … 393 429 // 2. update the objective function based on the new rounded solution 394 430 double offset=0.0; 395 for (i=0;i<numberIntegers;i++) { 396 int iColumn = integerVariable[i]; 397 if(!solver->isBinary(iColumn)) 431 double costValue = (1.0-scaleFactor)*solver->getObjSense(); 432 433 for (i=0;i<numberColumns;i++) { 434 // below so we can keep original code and allow for objective 435 int iColumn = i; 436 if(!solver->isBinary(iColumn)&&!general) 398 437 continue; 399 double costValue = 1.0;400 438 // deal with fixed variables (i.e., upper=lower) 401 439 if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) { 402 if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);403 else solver->setObjCoeff(iColumn,costValue);440 //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue); 441 //else solver->setObjCoeff(iColumn,costValue); 404 442 continue; 405 443 } 406 444 if (newSolution[iColumn]<lower[iColumn]+primalTolerance) { 407 solver->setObjCoeff(iColumn,costValue );445 solver->setObjCoeff(iColumn,costValue+scaleFactor*saveObjective[iColumn]); 408 446 } else { 409 447 if (newSolution[iColumn]>upper[iColumn]-primalTolerance) { 410 solver->setObjCoeff(iColumn,-costValue );448 solver->setObjCoeff(iColumn,-costValue+scaleFactor*saveObjective[iColumn]); 411 449 } else { 412 450 abort(); … … 416 454 } 417 455 solver->setDblParam(OsiObjOffset,-offset); 418 if (!general& false) {456 if (!general&&false) { 419 457 // Solve in two goes - first keep satisfied ones fixed 420 458 double * saveLower = new double [numberIntegers]; … … 444 482 memcpy(newSolution,solution,numberColumns*sizeof(double)); 445 483 int flip; 446 returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip); 484 returnCode = rounds(newSolution,saveObjective,numberIntegers,integerVariable, 485 roundExpensive_,downValue_,&flip); 447 486 if (returnCode) { 448 487 // solution - but may not be better … … 465 504 } 466 505 } 467 solver->resolve(); 468 assert (solver->isProvenOptimal()); 506 if (!general) { 507 solver->resolve(); 508 assert (solver->isProvenOptimal()); 509 } else { 510 int * addStart = new int[2*general+1]; 511 int * addIndex = new int[4*general]; 512 double * addElement = new double[4*general]; 513 double * addLower = new double[2*general]; 514 double * addUpper = new double[2*general]; 515 double * obj = new double[general]; 516 int nAdd=0; 517 for (i=0;i<numberIntegers;i++) { 518 int iColumn = integerVariable[i]; 519 if (newSolution[iColumn]>lower[iColumn]+primalTolerance&& 520 newSolution[iColumn]<upper[iColumn]-primalTolerance) { 521 obj[nAdd]=1.0; 522 addLower[nAdd]=0.0; 523 addUpper[nAdd]=COIN_DBL_MAX; 524 nAdd++; 525 } 526 } 527 OsiSolverInterface * solver2 = solver; 528 if (nAdd) { 529 CoinZeroN(addStart,nAdd+1); 530 solver2 = solver->clone(); 531 solver2->addCols(nAdd,addStart,NULL,NULL,addLower,addUpper,obj); 532 // feasible solution 533 double * sol = new double[nAdd+numberColumns]; 534 memcpy(sol,solution,numberColumns*sizeof(double)); 535 // now rows 536 int nAdd=0; 537 int nEl=0; 538 int nAddRow=0; 539 for (i=0;i<numberIntegers;i++) { 540 int iColumn = integerVariable[i]; 541 if (newSolution[iColumn]>lower[iColumn]+primalTolerance&& 542 newSolution[iColumn]<upper[iColumn]-primalTolerance) { 543 addLower[nAddRow]=-newSolution[iColumn];; 544 addUpper[nAddRow]=COIN_DBL_MAX; 545 addIndex[nEl] = iColumn; 546 addElement[nEl++]=-1.0; 547 addIndex[nEl] = numberColumns+nAdd; 548 addElement[nEl++]=1.0; 549 nAddRow++; 550 addStart[nAddRow]=nEl; 551 addLower[nAddRow]=newSolution[iColumn];; 552 addUpper[nAddRow]=COIN_DBL_MAX; 553 addIndex[nEl] = iColumn; 554 addElement[nEl++]=1.0; 555 addIndex[nEl] = numberColumns+nAdd; 556 addElement[nEl++]=1.0; 557 nAddRow++; 558 addStart[nAddRow]=nEl; 559 sol[nAdd+numberColumns] = fabs(sol[iColumn]-newSolution[iColumn]); 560 nAdd++; 561 } 562 } 563 solver2->setColSolution(sol); 564 delete [] sol; 565 solver2->addRows(nAddRow,addStart,addIndex,addElement,addLower,addUpper); 566 } 567 delete [] addStart; 568 delete [] addIndex; 569 delete [] addElement; 570 delete [] addLower; 571 delete [] addUpper; 572 delete [] obj; 573 solver2->resolve(); 574 assert (solver2->isProvenOptimal()); 575 if (nAdd) { 576 solver->setColSolution(solver2->getColSolution()); 577 delete solver2; 578 } 579 } 469 580 if (model_->logLevel()) 470 581 printf("\npass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue()); … … 472 583 473 584 } 585 // reduce scale factor 586 scaleFactor *= weightFactor_; 474 587 } // END WHILE 475 588 if (newLineNeeded&&model_->logLevel()) … … 491 604 for (i=0;i<numberIntegers;i++) { 492 605 int iColumn=integerVariable[i]; 493 const OsiObject * object = model_->object(i); 494 double originalLower; 495 double originalUpper; 496 getIntegerInformation( object,originalLower, originalUpper); 497 assert(colLower[iColumn]==originalLower); 498 newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower)); 499 assert(colUpper[iColumn]==originalUpper); 500 newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper)); 606 //const OsiObject * object = model_->object(i); 607 //double originalLower; 608 //double originalUpper; 609 //getIntegerInformation( object,originalLower, originalUpper); 610 //assert(colLower[iColumn]==originalLower); 611 //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower)); 612 newSolver->setColLower(iColumn,colLower[iColumn]); 613 //assert(colUpper[iColumn]==originalUpper); 614 //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper)); 615 newSolver->setColUpper(iColumn,colUpper[iColumn]); 501 616 if (!usedColumn[iColumn]) { 502 617 double value=lastSolution[iColumn]; … … 598 713 cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement()); 599 714 printf("round again with cutoff of %g\n",cutoff); 600 if ( accumulate_<2&&usedColumn)715 if ((accumulate_&3)<2&&usedColumn) 601 716 memset(usedColumn,0,numberColumns); 602 717 totalNumberPasses += numberPasses; … … 608 723 delete [] lastSolution; 609 724 delete [] newSolution; 725 delete [] integerVariable; 610 726 return finalReturnCode; 611 727 } … … 625 741 CbcHeuristicFPump::rounds(double * solution, 626 742 const double * objective, 743 int numberIntegers, const int * integerVariable, 627 744 bool roundExpensive, double downValue, int *flip) 628 745 { … … 631 748 double primalTolerance ; 632 749 solver->getDblParam(OsiPrimalTolerance,primalTolerance) ; 633 int numberIntegers = model_->numberIntegers();634 const int * integerVariable = model_->integerVariable();635 750 636 751 int i; … … 653 768 for (i=0;i<numberIntegers;i++) { 654 769 int iColumn = integerVariable[i]; 655 if(!solver->isBinary(iColumn))656 continue;657 #ifndef NDEBUG658 const OsiObject * object = model_->object(i);659 const CbcSimpleInteger * integerObject =660 dynamic_cast<const CbcSimpleInteger *> (object);661 const OsiSimpleInteger * integerObject2 =662 dynamic_cast<const OsiSimpleInteger *> (object);663 assert(integerObject||integerObject2);664 #endif665 770 double value=solution[iColumn]; 666 771 double round = floor(value+primalTolerance); … … 691 796 delete [] tmp; 692 797 798 const double * columnLower = solver->getColLower(); 799 const double * columnUpper = solver->getColUpper(); 693 800 if (*flip == 0 && iter != 0) { 694 801 if(model_->logLevel()) 695 802 printf(" -- rand = %4d (%4d) ", nnv, nn); 696 for (i = 0; i < nnv; i++) solution[list[i]] = 1. - solution[list[i]]; 803 for (i = 0; i < nnv; i++) { 804 // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6 805 int index = list[i]; 806 double value = solution[index]; 807 if (value<=1.0) { 808 solution[index] = 1.0-value; 809 } else if (value<columnLower[index]+integerTolerance) { 810 solution[index] = value+1.0; 811 } else if (value>columnUpper[index]-integerTolerance) { 812 solution[index] = value-1.0; 813 } else { 814 solution[index] = value-1.0; 815 } 816 } 697 817 *flip = nnv; 698 818 } else if (model_->logLevel()) {
Note: See TracChangeset
for help on using the changeset viewer.