Changeset 912 for trunk/Cbc/src/CbcHeuristicDiveCoefficient.cpp
 Timestamp:
 Apr 10, 2008 1:55:10 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Cbc/src/CbcHeuristicDiveCoefficient.cpp
r873 r912 6 6 #endif 7 7 8 //#define PRINT_DEBUG 9 8 10 #include "CbcHeuristicDiveCoefficient.hpp" 9 11 #include "CbcStrategy.hpp" 10 #include "CoinTime.hpp"11 12 12 13 // Default Constructor 13 14 CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient() 14 :CbcHeuristic ()15 :CbcHeuristicDive() 15 16 { 16 // matrix and row copy will automatically be empty17 downLocks_ =NULL;18 upLocks_ = NULL;19 percentageToFix_ = 0.2;20 maxIterations_ = 100;21 maxTime_ = 60;22 17 } 23 18 24 19 // Constructor from model 25 20 CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient(CbcModel & model) 26 :CbcHeuristic (model)21 :CbcHeuristicDive(model) 27 22 { 28 downLocks_ =NULL;29 upLocks_ = NULL;30 // Get a copy of original matrix31 assert(model.solver());32 // model may have empty matrix  wait until setModel33 const CoinPackedMatrix * matrix = model.solver()>getMatrixByCol();34 if (matrix) {35 matrix_ = *matrix;36 }37 percentageToFix_ = 0.2;38 maxIterations_ = 100;39 maxTime_ = 60;40 23 } 41 24 … … 43 26 CbcHeuristicDiveCoefficient::~CbcHeuristicDiveCoefficient () 44 27 { 45 delete [] downLocks_;46 delete [] upLocks_;47 28 } 48 29 … … 62 43 fprintf(fp,"3 CbcHeuristicDiveCoefficient heuristicDiveCoefficient(*cbcModel);\n"); 63 44 CbcHeuristic::generateCpp(fp,"heuristicDiveCoefficient"); 64 if (percentageToFix_!=other.percentageToFix_)65 fprintf(fp,"3 heuristicDiveCoefficient.setPercentageToFix(%.2f);\n",percentageToFix_);66 else67 fprintf(fp,"4 heuristicDiveCoefficient.setPercentageToFix(%.2f);\n",percentageToFix_);68 if (maxIterations_!=other.maxIterations_)69 fprintf(fp,"3 heuristicDiveCoefficient.setMaxIterations(%d);\n",maxIterations_);70 else71 fprintf(fp,"4 heuristicDiveCoefficient.setMaxIterations(%d);\n",maxIterations_);72 if (maxTime_!=other.maxTime_)73 fprintf(fp,"3 heuristicDiveCoefficient.setMaxTime(%.2f);\n",maxTime_);74 else75 fprintf(fp,"4 heuristicDiveCoefficient.setMaxTime(%.2f);\n",maxTime_);76 45 fprintf(fp,"3 cbcModel>addHeuristic(&heuristicDiveCoefficient);\n"); 77 46 } … … 80 49 CbcHeuristicDiveCoefficient::CbcHeuristicDiveCoefficient(const CbcHeuristicDiveCoefficient & rhs) 81 50 : 82 CbcHeuristic(rhs), 83 matrix_(rhs.matrix_), 84 percentageToFix_(rhs.percentageToFix_), 85 maxIterations_(rhs.maxIterations_), 86 maxTime_(rhs.maxTime_) 51 CbcHeuristicDive(rhs) 87 52 { 88 int numberIntegers = model_>numberIntegers();89 if (rhs.downLocks_) {90 int numberIntegers = model_>numberIntegers();91 downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);92 upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);93 } else {94 downLocks_ = NULL;95 upLocks_ = NULL;96 }97 53 } 98 54 … … 102 58 { 103 59 if (this!=&rhs) { 104 CbcHeuristic::operator=(rhs); 105 matrix_ = rhs.matrix_; 106 percentageToFix_ = rhs.percentageToFix_; 107 maxIterations_ = rhs.maxIterations_; 108 maxTime_ = rhs.maxTime_; 109 delete [] downLocks_; 110 delete [] upLocks_; 111 if (rhs.downLocks_) { 112 int numberIntegers = model_>numberIntegers(); 113 downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers); 114 upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers); 115 } else { 116 downLocks_ = NULL; 117 upLocks_ = NULL; 118 } 60 CbcHeuristicDive::operator=(rhs); 119 61 } 120 62 return *this; 121 63 } 122 64 123 // Resets stuff if model changes 124 void 125 CbcHeuristicDiveCoefficient::resetModel(CbcModel * model) 65 void 66 CbcHeuristicDiveCoefficient::selectVariableToBranch(OsiSolverInterface* solver, 67 const double* newSolution, 68 int& bestColumn, 69 int& bestRound) 126 70 { 127 model_=model;128 // Get a copy of original matrix129 assert(model_>solver());130 // model may have empty matrix  wait until setModel131 const CoinPackedMatrix * matrix = model_>solver()>getMatrixByCol();132 if (matrix) {133 matrix_ = *matrix;134 validate();135 }136 }137 138 // See if dive fractional will give better solution139 // Sets value of solution140 // Returns 1 if solution, 0 if not141 int142 CbcHeuristicDiveCoefficient::solution(double & solutionValue,143 double * betterSolution)144 {145 146 // See if to do147 if (!when()(when()%10==1&&model_>phase()!=1)148 (when()%10==2&&(model_>phase()!=2&&model_>phase()!=3)))149 return 0; // switched off150 151 double time1 = CoinCpuTime();152 153 OsiSolverInterface * solver = model_>solver()>clone();154 const double * lower = solver>getColLower();155 const double * upper = solver>getColUpper();156 const double * rowLower = solver>getRowLower();157 const double * rowUpper = solver>getRowUpper();158 const double * solution = solver>getColSolution();159 const double * objective = solver>getObjCoefficients();160 double integerTolerance = model_>getDblParam(CbcModel::CbcIntegerTolerance);161 double primalTolerance;162 solver>getDblParam(OsiPrimalTolerance,primalTolerance);163 164 int numberRows = matrix_.getNumRows();165 assert (numberRows<=solver>getNumRows());166 71 int numberIntegers = model_>numberIntegers(); 167 72 const int * integerVariable = model_>integerVariable(); 168 double direction = solver>getObjSense(); // 1 for min, 1 for max 169 double newSolutionValue = direction*solver>getObjValue(); 170 int returnCode = 0; 171 // Column copy 172 const double * element = matrix_.getElements(); 173 const int * row = matrix_.getIndices(); 174 const CoinBigIndex * columnStart = matrix_.getVectorStarts(); 175 const int * columnLength = matrix_.getVectorLengths(); 73 double integerTolerance = model_>getDblParam(CbcModel::CbcIntegerTolerance); 176 74 177 // Get solution array for heuristic solution 178 int numberColumns = solver>getNumCols(); 179 double * newSolution = new double [numberColumns]; 180 memcpy(newSolution,solution,numberColumns*sizeof(double)); 181 182 // vectors to store the latest variables fixed at their bounds 183 int* columnFixed = new int [numberIntegers]; 184 double* originalBound = new double [numberIntegers]; 185 bool * fixedAtLowerBound = new bool [numberIntegers]; 186 187 const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers); 188 189 // count how many fractional variables 190 int numberFractionalVariables = 0; 75 bestColumn = 1; 76 bestRound = 1; // 1 rounds down, +1 rounds up 77 double bestFraction = DBL_MAX; 78 int bestLocks = COIN_INT_MAX; 191 79 for (int i=0; i<numberIntegers; i++) { 192 80 int iColumn = integerVariable[i]; 193 81 double value=newSolution[iColumn]; 82 double fraction=valuefloor(value); 83 int round = 0; 194 84 if (fabs(floor(value+0.5)value)>integerTolerance) { 195 numberFractionalVariables++; 196 } 197 } 198 199 int iteration = 0; 200 while(numberFractionalVariables) { 201 iteration++; 202 203 // select a fractional variable to bound 204 double bestFraction = DBL_MAX; 205 int bestColumn = 1; 206 int bestRound = 1; // 1 rounds down, +1 rounds up 207 int bestLocks = COIN_INT_MAX; 208 int numberAtBoundFixed = 0; 209 bool canRoundSolution = true; 210 for (int i=0; i<numberIntegers; i++) { 211 int iColumn = integerVariable[i]; 212 double value=newSolution[iColumn]; 213 double fraction=valuefloor(value); 214 int round = 0; 215 if (fabs(floor(value+0.5)value)>integerTolerance) { 216 int nDownLocks = downLocks_[i]; 217 int nUpLocks = upLocks_[i]; 218 if(nDownLocks>0 && nUpLocks>0) { 219 // the variable cannot be rounded 220 canRoundSolution = false; 221 int nLocks = nDownLocks; 222 if(nDownLocks<nUpLocks) 223 round = 1; 224 else if(nDownLocks>nUpLocks) { 225 round = 1; 226 fraction = 1.0  fraction; 227 nLocks = nUpLocks; 228 } 229 else if(fraction < 0.5) 230 round = 1; 231 else { 232 round = 1; 233 fraction = 1.0  fraction; 234 nLocks = nUpLocks; 235 } 236 237 // if variable is not binary, penalize it 238 if(!solver>isBinary(iColumn)) 239 fraction *= 1000.0; 240 241 if(nLocks < bestLocks  (nLocks == bestLocks && 242 fraction < bestFraction)) { 243 bestColumn = iColumn; 244 bestLocks = nLocks; 245 bestFraction = fraction; 246 bestRound = round; 247 } 85 int nDownLocks = downLocks_[i]; 86 int nUpLocks = upLocks_[i]; 87 if(nDownLocks>0 && nUpLocks>0) { 88 // the variable cannot be rounded 89 int nLocks = nDownLocks; 90 if(nDownLocks<nUpLocks) 91 round = 1; 92 else if(nDownLocks>nUpLocks) { 93 round = 1; 94 fraction = 1.0  fraction; 95 nLocks = nUpLocks; 248 96 } 249 } 250 else if(numberAtBoundFixed < maxNumberAtBoundToFix) { 251 // fix the variable at one of its bounds 252 if (fabs(lower[iColumn]value)<=integerTolerance && 253 lower[iColumn] != upper[iColumn]) { 254 columnFixed[numberAtBoundFixed] = iColumn; 255 originalBound[numberAtBoundFixed] = upper[iColumn]; 256 fixedAtLowerBound[numberAtBoundFixed] = true; 257 solver>setColUpper(iColumn, lower[iColumn]); 258 numberAtBoundFixed++; 97 else if(fraction < 0.5) 98 round = 1; 99 else { 100 round = 1; 101 fraction = 1.0  fraction; 102 nLocks = nUpLocks; 259 103 } 260 else if(fabs(upper[iColumn]value)<=integerTolerance && 261 lower[iColumn] != upper[iColumn]) { 262 columnFixed[numberAtBoundFixed] = iColumn; 263 originalBound[numberAtBoundFixed] = lower[iColumn]; 264 fixedAtLowerBound[numberAtBoundFixed] = false; 265 solver>setColLower(iColumn, upper[iColumn]); 266 numberAtBoundFixed++; 104 105 // if variable is not binary, penalize it 106 if(!solver>isBinary(iColumn)) 107 fraction *= 1000.0; 108 109 if(nLocks < bestLocks  (nLocks == bestLocks && 110 fraction < bestFraction)) { 111 bestColumn = iColumn; 112 bestLocks = nLocks; 113 bestFraction = fraction; 114 bestRound = round; 267 115 } 268 116 } 269 117 } 270 271 if(canRoundSolution) {272 // Round all the fractional variables273 for (int i=0; i<numberIntegers; i++) {274 int iColumn = integerVariable[i];275 double value=newSolution[iColumn];276 if (fabs(floor(value+0.5)value)>integerTolerance) {277 if(downLocks_[i]==0  upLocks_[i]==0) {278 if(downLocks_[i]==0 && upLocks_[i]==0) {279 if(direction * objective[iColumn] >= 0.0)280 newSolution[iColumn] = floor(value);281 else282 newSolution[iColumn] = ceil(value);283 }284 else if(downLocks_[i]==0)285 newSolution[iColumn] = floor(value);286 else287 newSolution[iColumn] = ceil(value);288 }289 else290 break;291 }292 }293 break;294 }295 296 297 double originalBoundBestColumn;298 if(bestColumn >= 0) {299 if(bestRound < 0) {300 originalBoundBestColumn = upper[bestColumn];301 solver>setColUpper(bestColumn, floor(newSolution[bestColumn]));302 }303 else {304 originalBoundBestColumn = lower[bestColumn];305 solver>setColLower(bestColumn, ceil(newSolution[bestColumn]));306 }307 } else {308 break;309 }310 int originalBestRound = bestRound;311 while (1) {312 313 solver>resolve();314 315 if(!solver>isProvenOptimal()) {316 if(numberAtBoundFixed > 0) {317 // Remove the bound fix for variables that were at bounds318 for(int i=0; i<numberAtBoundFixed; i++) {319 int iColFixed = columnFixed[i];320 if(fixedAtLowerBound[i])321 solver>setColUpper(iColFixed, originalBound[i]);322 else323 solver>setColLower(iColFixed, originalBound[i]);324 }325 numberAtBoundFixed = 0;326 }327 else if(bestRound == originalBestRound) {328 bestRound *= (1);329 if(bestRound < 0) {330 solver>setColLower(bestColumn, originalBoundBestColumn);331 solver>setColUpper(bestColumn, floor(newSolution[bestColumn]));332 }333 else {334 solver>setColLower(bestColumn, ceil(newSolution[bestColumn]));335 solver>setColUpper(bestColumn, originalBoundBestColumn);336 }337 }338 else339 break;340 }341 else342 break;343 }344 345 if(!solver>isProvenOptimal())346 break;347 348 if(iteration > maxIterations_) {349 break;350 }351 352 if(CoinCpuTime()time1 > maxTime_) {353 break;354 }355 356 memcpy(newSolution,solution,numberColumns*sizeof(double));357 numberFractionalVariables = 0;358 for (int i=0; i<numberIntegers; i++) {359 int iColumn = integerVariable[i];360 double value=newSolution[iColumn];361 if (fabs(floor(value+0.5)value)>integerTolerance) {362 numberFractionalVariables++;363 }364 }365 366 }367 368 369 double * rowActivity = new double[numberRows];370 memset(rowActivity,0,numberRows*sizeof(double));371 372 // recompute new solution value373 double objOffset=0.0;374 solver>getDblParam(OsiObjOffset,objOffset);375 newSolutionValue = objOffset;376 for (int i=0 ; i<numberColumns ; i++ )377 newSolutionValue += objective[i]*newSolution[i];378 newSolutionValue *= direction;379 //printf("new solution value %g %g\n",newSolutionValue,solutionValue);380 if (newSolutionValue<solutionValue) {381 // paranoid check382 memset(rowActivity,0,numberRows*sizeof(double));383 for (int i=0;i<numberColumns;i++) {384 int j;385 double value = newSolution[i];386 if (value) {387 for (j=columnStart[i];388 j<columnStart[i]+columnLength[i];j++) {389 int iRow=row[j];390 rowActivity[iRow] += value*element[j];391 }392 }393 }394 // check was approximately feasible395 bool feasible=true;396 for (int i=0;i<numberRows;i++) {397 if(rowActivity[i]<rowLower[i]) {398 if (rowActivity[i]<rowLower[i]1000.0*primalTolerance)399 feasible = false;400 } else if(rowActivity[i]>rowUpper[i]) {401 if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)402 feasible = false;403 }404 }405 for (int i=0; i<numberIntegers; i++) {406 int iColumn = integerVariable[i];407 double value=newSolution[iColumn];408 if (fabs(floor(value+0.5)value)>integerTolerance) {409 feasible = false;410 break;411 }412 }413 if (feasible) {414 // new solution415 memcpy(betterSolution,newSolution,numberColumns*sizeof(double));416 solutionValue = newSolutionValue;417 //printf("** Solution of %g found by CbcHeuristicDiveCoefficient\n",newSolutionValue);418 returnCode=1;419 } else {420 // Can easily happen421 //printf("Debug CbcHeuristicDiveCoefficient giving bad solution\n");422 }423 }424 425 delete [] newSolution;426 delete [] columnFixed;427 delete [] originalBound;428 delete [] fixedAtLowerBound;429 delete [] rowActivity;430 delete solver;431 return returnCode;432 }433 434 // update model435 void CbcHeuristicDiveCoefficient::setModel(CbcModel * model)436 {437 model_ = model;438 // Get a copy of original matrix (and by row for rounding);439 assert(model_>solver());440 const CoinPackedMatrix * matrix = model_>solver()>getMatrixByCol();441 if (matrix) {442 matrix_ = *matrix;443 // make sure model okay for heuristic444 validate();445 118 } 446 119 } 447 448 // Validate model i.e. sets when_ to 0 if necessary (may be NULL)449 void450 CbcHeuristicDiveCoefficient::validate()451 {452 if (model_&&when()<10) {453 if (model_>numberIntegers()!=454 model_>numberObjects())455 setWhen(0);456 }457 458 int numberIntegers = model_>numberIntegers();459 const int * integerVariable = model_>integerVariable();460 delete [] downLocks_;461 delete [] upLocks_;462 downLocks_ = new unsigned short [numberIntegers];463 upLocks_ = new unsigned short [numberIntegers];464 // Column copy465 const double * element = matrix_.getElements();466 const int * row = matrix_.getIndices();467 const CoinBigIndex * columnStart = matrix_.getVectorStarts();468 const int * columnLength = matrix_.getVectorLengths();469 const double * rowLower = model_>solver()>getRowLower();470 const double * rowUpper = model_>solver()>getRowUpper();471 for (int i=0;i<numberIntegers;i++) {472 int iColumn = integerVariable[i];473 int down=0;474 int up=0;475 if (columnLength[iColumn]>65535) {476 setWhen(0);477 break; // unlikely to work478 }479 for (CoinBigIndex j=columnStart[iColumn];480 j<columnStart[iColumn]+columnLength[iColumn];j++) {481 int iRow=row[j];482 if (rowLower[iRow]>1.0e20&&rowUpper[iRow]<1.0e20) {483 up++;484 down++;485 } else if (element[j]>0.0) {486 if (rowUpper[iRow]<1.0e20)487 up++;488 else489 down++;490 } else {491 if (rowLower[iRow]>1.0e20)492 up++;493 else494 down++;495 }496 }497 downLocks_[i] = (unsigned short) down;498 upLocks_[i] = (unsigned short) up;499 }500 }
Note: See TracChangeset
for help on using the changeset viewer.