Changeset 1065
 Timestamp:
 Sep 12, 2008 5:08:57 PM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Cbc/src/CbcHeuristicPivotAndFix.cpp
r1058 r1065 9 9 #include <cmath> 10 10 #include <cfloat> 11 #include <vector> 11 12 12 13 #include "OsiSolverInterface.hpp" … … 14 15 #include "CbcMessage.hpp" 15 16 #include "CbcHeuristicPivotAndFix.hpp" 17 #include "OsiClpSolverInterface.hpp" 18 #include "CoinTime.hpp" 19 20 //#define FORNOW 16 21 17 22 // Default Constructor … … 81 86 { 82 87 83 numCouldRun_++; 84 printf("entered pivot and fix\n"); 85 int returnCode = 0; 86 return returnCode; 87 } 88 numCouldRun_++; // Todo: Ask JJHF what this for. 89 std::cout << "Entering PivotandFix Heuristic" << std::endl; 90 91 #ifdef FORNOW 92 std::cout << "Lucky you! You're in the PivotandFix Heuristic" << std::endl; 93 // The struct should be moved to member data 94 95 96 97 typedef struct { 98 int numberSolutions; 99 int maximumSolutions; 100 int numberColumns; 101 double ** solution; 102 int * numberUnsatisfied; 103 } clpSolution; 104 105 double start = CoinCpuTime(); 106 107 OsiClpSolverInterface * clpSolverOriginal 108 = dynamic_cast<OsiClpSolverInterface *> (model_>solver()); 109 assert (clpSolverOriginal); 110 OsiClpSolverInterface *clpSolver(clpSolverOriginal); 111 112 ClpSimplex * simplex = clpSolver>getModelPtr(); 113 114 // Initialize the structure holding the solutions 115 clpSolution solutions; 116 // Set typeStruct field of ClpTrustedData struct to one. 117 // This tells Clp it's "Mahdi!" 118 ClpTrustedData trustedSolutions; 119 trustedSolutions.typeStruct = 1; 120 trustedSolutions.data = &solutions; 121 solutions.numberSolutions=0; 122 solutions.maximumSolutions=0; 123 solutions.numberColumns=simplex>numberColumns(); 124 solutions.solution=NULL; 125 solutions.numberUnsatisfied=NULL; 126 simplex>setTrustedUserPointer(&trustedSolutions); 127 128 // Solve from all slack to get some points 129 simplex>allSlackBasis(); 130 simplex>primal(); 131 132 //  133 // Get the problem information 134 //  get the number of cols and rows 135 int numCols = clpSolver>getNumCols(); 136 int numRows = clpSolver>getNumRows(); 137 138 //  get the right hand side of the rows 139 const double * rhs = clpSolver>getRightHandSide(); 140 141 //  find the integer variables 142 bool * varClassInt = new bool[numCols]; 143 int numInt = 0; 144 for(int i=0; i<numCols; i++) 145 { 146 if(clpSolver>isContinuous(i)) 147 varClassInt[i] = 0; 148 else 149 { 150 varClassInt[i] = 1; 151 numInt++; 152 } 153 } 154 155 // Get the rows sense 156 const char * rowSense; 157 rowSense = clpSolver>getRowSense(); 158 159 // Get the objective coefficients 160 const double *objCoefficients = clpSolver>getObjCoefficients(); 161 double *originalObjCoeff = new double [numCols]; 162 for(int i=0; i<numCols;i++) 163 originalObjCoeff[i] = objCoefficients[i]; 164 165 // Get the matrix of the problem 166 double ** matrix = new double * [numRows]; 167 for(int i = 0; i < numRows; i++) 168 { 169 matrix[i] = new double[numCols]; 170 for(int j = 0; j < numCols; j++) 171 matrix[i][j] = 0; 172 } 173 const CoinPackedMatrix* matrixByRow = clpSolver>getMatrixByRow(); 174 const double * matrixElements = matrixByRow>getElements(); 175 const int * matrixIndices = matrixByRow>getIndices(); 176 const int * matrixStarts = matrixByRow>getVectorStarts(); 177 for(int j = 0; j < numRows; j++) 178 { 179 for(int i = matrixStarts[j]; i < matrixStarts[j+1]; i++) 180 { 181 matrix[j][matrixIndices[i]] = matrixElements[i]; 182 } 183 } 184 185 // The newObj is the randomly perturbed constraint used to find new 186 // corner points 187 double * newObj = new double [numCols]; 188 189 // Set the random seed 190 srand ( time(NULL) +1); 191 int randNum; 192 193 // We're going to add a new row to the LP formulation 194 // after finding each new solution. 195 // Adding a new row requires the new elements and the new indices. 196 // The elements are original objective function coefficients. 197 // The indicies are the (dense) columns indices stored in addRowIndex. 198 // The rhs is the value of the new solution stored in solutionValue. 199 int * addRowIndex = new int[numCols]; 200 for(int i=0;i<numCols;i++) 201 addRowIndex[i]=i; 202 203 // The number of feasible solutions found by the PF heuristic. 204 // This controls the return code of the solution() method. 205 int numFeasibles =0; 206 207 // Shuffle the rows 208 int * index = new int [numRows]; 209 for(int i=0; i<numRows; i++) 210 index[i]=i; 211 for(int i=0; i<numRows; i++) 212 { 213 int temp = index[i]; 214 int randNumTemp = i+(rand()%(numRowsi)); 215 index[i]=index[randNumTemp]; 216 index[randNumTemp]=temp; 217 } 218 219 // In the clpSolution struct, we store a lot of column solutions. 220 // For each perturb objective, we store the solution from each 221 // iteration of the LP solve. 222 // For each perturb objective, we look at the collection of 223 // solutions to do something extremly intelligent :) 224 // We could (and should..and will :) wipe out the block of 225 // solutions when we're done with them. But for now, we just move on 226 // and store the next block of solutions for the next (perturbed) 227 // objective. 228 // The variable startIndex tells us where the new block begins. 229 int startIndex = 0; 230 231 // At most "fixThreshold" number of integer variables can be unsatisfied 232 // for calling smallBranchAndBound(). 233 // The PF Heuristic only fixes fixThreshold number of variables to 234 // their integer values. Not more. Not less. The reason is to give 235 // the smallBB some opportunity to find better solutions. If we fix 236 // everything it might be too many (leading the heuristic to come up 237 // with infeasibility rather than a useful result). 238 // (This is an important paramater. And it is dynamically set.) 239 double fixThreshold; 240 /* 241 if(numInt > 400) 242 fixThreshold = 17*sqrt(numInt); 243 if(numInt<=400 && numInt>100) 244 fixThreshold = 5*sqrt(numInt); 245 if(numInt<=100) 246 fixThreshold = 4*sqrt(numInt); 247 */ 248 // Initialize fixThreshold based on the number of integer 249 // variables 250 if(numInt<=100) 251 fixThreshold = .35*numInt; 252 if(numInt>100 && numInt<1000) 253 fixThreshold = .85*numInt; 254 if(numInt>=1000) 255 fixThreshold = .1*numInt; 256 257 // Whenever the dynamic system for changing fixThreshold 258 // kicks in, it changes the parameter by the 259 // fixThresholdChange amount. 260 // (The 25% should be member data and tuned. Another paper!) 261 double fixThresholdChange = 0.25 * fixThreshold; 262 263 // maxNode is the maximum number of nodes we allow smallBB to 264 // search. It's initialized to 400 and changed dynamically. 265 // The 400 should be member data, if we become virtuous. 266 int maxNode = 400; 267 268 // We control the decision to change maxNode through the boolean 269 // variable changeMaxNode. The boolean variable is initialized to 270 // true and gets set to false under a condition (and is never true 271 // again.) 272 // It's flipped off and stays off (in the current incarnation of PF) 273 bool changeMaxNode = 1; 274 275 // The sumReturnCode is used for the dynamic system that sets 276 // fixThreshold and changeMaxNode. 277 // 278 // We track what's happening in sumReturnCode. There are 8 switches. 279 // The first 5 switches corresponds to a return code for smallBB. 280 // 281 // We want to know how many times we consecutively get the same 282 // return code. 283 // 284 // If "good" return codes are happening often enough, we're happy. 285 // 286 // If a "bad" returncodes happen consecutively, we want to 287 // change something. 288 // 289 // The switch 5 is the number of times PF didn't call smallBB 290 // becuase the number of integer variables that took integer values 291 // was less than fixThreshold. 292 // 293 // The swicth 6 was added for a brilliant idea...to be announced 294 // later (another paper!) 295 // 296 // The switch 7 is the one that changes the max node. Read the 297 // code. (Todo: Verbalize the brilliant idea for the masses.) 298 // 299 int sumReturnCode[8]; 300 /* 301 sumReturnCode[0] ~ 1 > problem too big for smallBB 302 sumReturnCode[1] ~ 0 > smallBB not finshed and no soln 303 sumReturnCode[2] ~ 1 > smallBB not finshed and there is a soln 304 sumReturnCode[3] ~ 2 > smallBB finished and no soln 305 sumReturnCode[4] ~ 3 > smallBB finished and there is a soln 306 sumReturnCode[5] ~ didn't call smallBranchAndBound too few to fix 307 sumReturnCode[6] ~ didn't call smallBranchAndBound too many unsatisfied 308 sumReturnCode[7] ~ the same as sumReturnCode[1] but becomes zero just if the returnCode is not 0 309 */ 310 311 for(int i=0; i<8; i++) 312 sumReturnCode[i]=0; 313 int * colIndex = new int[numCols]; 314 for(int i=0; i<numCols; i++) 315 colIndex[i]=i; 316 double cutoff = COIN_DBL_MAX; 317 bool didMiniBB; 318 319 // Main loop 320 for(int i=0; i<numRows; i++) 321 { 322 // track the number of minibb for the dynamic threshold setting 323 didMiniBB=0; 324 325 for(int k=startIndex; k<solutions.numberSolutions; k++) 326 //if the point has 0 unsatisfied variables; make sure it is 327 //feasible. Check integer feasiblity and constraints. 328 if(solutions.numberUnsatisfied[k] == 0) 329 { 330 double feasibility=1; 331 //check integer feasibility 332 for(int icol=0; icol<numCols; icol++) 333 { 334 double closest = floor(solutions.solution[k][icol]+0.5); 335 if(varClassInt[icol]&&(fabs(solutions.solution[k][icol]closest)>1e6)) 336 { 337 feasibility = 0; 338 break; 339 } 340 } 341 //check if the solution satisfies the constraints 342 for(int irow = 0; irow < numRows; irow++) 343 { 344 double lhs = 0; 345 for(int j = 0; j <numCols; j++) 346 lhs += matrix[irow][j] * solutions.solution[k][j]; 347 if(rowSense[irow] == 'L' && lhs > rhs[irow] + 1e6) 348 { 349 feasibility = 0; 350 break; 351 } 352 if(rowSense[irow] == 'G' && lhs < rhs[irow]  1e6) 353 { 354 feasibility = 0; 355 break; 356 } 357 if(rowSense[irow] == 'E' && (lhs  rhs[irow] > 1e6  lhs  rhs[irow] < 1e6)) 358 { 359 feasibility = 0; 360 break; 361 } 362 } 363 364 //if feasible, find the objective value and set the cutoff 365 // for the smallBB and add a new constraint to the LP 366 // (and update the best solution found so far for the 367 // return arguments) 368 if(feasibility) 369 { 370 double objectiveValue=0; 371 for(int j=0; j<numCols; j++) 372 objectiveValue+=solutions.solution[k][j]*originalObjCoeff[j]; 373 cutoff=objectiveValue; 374 clpSolver>addRow(numCols, addRowIndex, originalObjCoeff,COIN_DBL_MAX, cutoff); 375 376 // Todo: pick up the best solution in the block (not 377 // the last). 378 solutionValue = objectiveValue; 379 for(int m=0; m<numCols; m++) 380 betterSolution[m] = solutions.solution[k][m]; 381 numFeasibles++; 382 } 383 } 384 385 // Go through the block of solution and decide if to call smallBB 386 for(int k=startIndex; k<solutions.numberSolutions; k++) 387 { 388 if(solutions.numberUnsatisfied[k] <= fixThreshold) 389 { 390 // get new copy 391 OsiSolverInterface * newSolver; 392 newSolver = new OsiClpSolverInterface(*clpSolver); 393 newSolver>setObjSense(1); 394 newSolver>setObjective(originalObjCoeff); 395 int numberColumns = newSolver>getNumCols(); 396 int numFixed=0; 397 398 // Fix the first fixThreshold number of integer vars 399 // that are satisfied 400 for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) 401 { 402 if (newSolver>isInteger(iColumn)) 403 { 404 double value = solutions.solution[k][iColumn]; 405 double intValue = floor(value+0.5); 406 if (fabs(valueintValue)<1.0e5) 407 { 408 newSolver>setColLower(iColumn,intValue); 409 newSolver>setColUpper(iColumn,intValue); 410 numFixed++; 411 if(numFixed>numIntfixThreshold) 412 break; 413 } 414 } 415 } 416 printf("numFixed: %d\n", numFixed); 417 printf("fixThreshold: %f\n", fixThreshold); 418 printf("numInt: %d\n", numInt); 419 double *newSolution = new double[numCols]; 420 double newSolutionValue; 421 422 // Call smallBB on the modified problem 423 int returnCode = smallBranchAndBound(newSolver, maxNode, newSolution, 424 newSolutionValue, cutoff, "mini"); 425 426 // If smallBB found a solution, update the better 427 // solution and solutionValue (we gave smallBB our 428 // cutoff, so it only finds improving solutions) 429 if(returnCode == 1  returnCode == 3) 430 { 431 numFeasibles ++; 432 solutionValue = newSolutionValue; 433 for(int m=0; m<numCols; m++) 434 betterSolution[m] = newSolution[m]; 435 printf("cutoff: %f\n", newSolutionValue); 436 printf("time: %.2lf\n", CoinCpuTime()start); 437 } 438 didMiniBB=1; 439 printf("returnCode: %d\n", returnCode); 440 441 //Update sumReturnCode array 442 for(int iRC=0; iRC<6; iRC++) 443 { 444 if(iRC == returnCode+1) 445 sumReturnCode[iRC]++; 446 else 447 sumReturnCode[iRC]=0; 448 } 449 if(returnCode!=0) 450 sumReturnCode[7]=0; 451 else 452 sumReturnCode[7]++; 453 if(returnCode==1  returnCode==3) 454 { 455 cutoff = newSolutionValue; 456 clpSolver>addRow(numCols, addRowIndex, originalObjCoeff,COIN_DBL_MAX, cutoff); 457 printf("******************\n\n*****************\n"); 458 } 459 break; 460 } 461 } 462 463 if(!didMiniBB && solutions.numberSolutions  startIndex > 0) 464 { 465 sumReturnCode[5]++; 466 for(int iRC=0; iRC<5; iRC++) 467 sumReturnCode[iRC]=0; 468 } 469 470 //Change "fixThreshold" if needed 471 // using the data we've recorded in sumReturnCode 472 if(sumReturnCode[1]>=3) 473 fixThreshold = fixThresholdChange; 474 if(sumReturnCode[7]>=3 && changeMaxNode) 475 { 476 maxNode *= 5; 477 changeMaxNode=0; 478 } 479 if(sumReturnCode[3]>=3 && fixThreshold < 0.95 * numInt) 480 fixThreshold += fixThresholdChange; 481 if(sumReturnCode[5]>=4) 482 fixThreshold += fixThresholdChange; 483 if(sumReturnCode[0]>3) 484 fixThreshold = fixThresholdChange; 485 486 startIndex = solutions.numberSolutions; 487 488 //Check if the maximum iterations limit is reached 489 // rlh: Ask John how this is working with the change to trustedUserPtr. 490 if(solutions.numberSolutions > 20000) 491 break; 492 493 // The first time in this loop PF solves orig LP. 494 495 //Generate the random objective function 496 randNum = rand()%10 + 1; 497 randNum = fmod(randNum, 2); 498 for(int j=0; j<numCols; j++) 499 { 500 if(randNum == 1) 501 if(fabs(matrix[index[i]][j]) < 1e6) 502 newObj[j] = 0.1; 503 else 504 newObj[j] = matrix[index[i]][j] * 1.1; 505 else 506 if(fabs(matrix[index[i]][j]) < 1e6) 507 newObj[j] = 0.1; 508 else 509 newObj[j] = matrix[index[i]][j] * 0.9; 510 } 511 clpSolver>setObjective(newObj); 512 if(rowSense[i] == 'L') 513 clpSolver>setObjSense(1); 514 else 515 // Todo #1: We don't need to solve the LPs to optimality. 516 // We just need corner points. 517 // There's a problem in stopping Clp that needs to be looked 518 // into. So for now, we solve optimality. 519 clpSolver>setObjSense(1); 520 // simplex>setMaximumIterations(100); 521 clpSolver>getModelPtr()>primal(1); 522 // simplex>setMaximumIterations(100000); 523 printf("cutoff: %f\n", cutoff); 524 printf("time: %.2f\n", CoinCpuTime()start); 525 for(int iRC=0; iRC<8; iRC++) 526 printf("%d ",sumReturnCode[iRC]); 527 printf("\nfixThreshold: %f\n", fixThreshold); 528 printf("numInt: %d\n", numInt); 529 printf("\n %d\n", i); 530 531 //temp: 532 if(i>3) break; 533 534 } 535 536 printf("Best Feasible Found: %f\n", cutoff); 537 printf("Total time: %.2f\n", CoinCpuTime()start); 538 539 if (numFeasibles == 0) { 540 return 0; 541 } 542 543 544 545 // We found something better 546 std::cout << "See you soon! You're leaving the PivotandFix Heuristic" << std::endl; 547 std::cout << std::endl; 548 549 return 1; 550 #endif 551 552 return 0; 553 554 } 555 556 557 558 88 559 // update model 89 560 void CbcHeuristicPivotAndFix::setModel(CbcModel * model)
Note: See TracChangeset
for help on using the changeset viewer.