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

 1 edited
Legend:
 Unmodified
 Added
 Removed

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