Changeset 1064
 Timestamp:
 Sep 11, 2008 4:56:14 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Cbc/src/CbcHeuristicRandRound.cpp
r1057 r1064 9 9 #include <cmath> 10 10 #include <cfloat> 11 #include <vector.h> 11 12 12 13 #include "OsiSolverInterface.hpp" … … 14 15 #include "CbcMessage.hpp" 15 16 #include "CbcHeuristicRandRound.hpp" 17 #include "OsiClpSolverInterface.hpp" 18 #include "CoinTime.hpp" 16 19 17 20 // Default Constructor … … 74 77 } 75 78 /* 76 Comments needed79 Randomized Rounding Heuristic 77 80 Returns 1 if solution, 0 if not */ 78 81 int … … 80 83 double * betterSolution) 81 84 { 82 83 numCouldRun_++; 84 int returnCode = 0; 85 return returnCode; 85 std::cout << "Lucky you! You're in the Randomized Rounding Heuristic" << std::endl; 86 // The struct should be moved to member data 87 88 typedef struct { 89 int numberSolutions; 90 int maximumSolutions; 91 int numberColumns; 92 double ** solution; 93 int * numberUnsatisfied; 94 } clpSolution; 95 96 double start = CoinCpuTime(); 97 numCouldRun_++; // Todo: Ask JJHF what this for. 98 99 OsiClpSolverInterface * clpSolver 100 = dynamic_cast<OsiClpSolverInterface *> (model_>solver()); 101 assert (clpSolver); 102 ClpSimplex * simplex = clpSolver>getModelPtr(); 103 104 // Initialize the structure holding the solutions 105 clpSolution solutions; 106 // Set typeStruct field of ClpTrustedData struct to one. 107 // This tells Clp it's "Mahdi!" 108 ClpTrustedData trustedSolutions; 109 trustedSolutions.typeStruct = 1; 110 trustedSolutions.data = &solutions; 111 solutions.numberSolutions=0; 112 solutions.maximumSolutions=0; 113 solutions.numberColumns=simplex>numberColumns(); 114 solutions.solution=NULL; 115 solutions.numberUnsatisfied=NULL; 116 simplex>setTrustedUserPointer(&trustedSolutions); // rlh: old was 117 // "userPointer" 118 // Solve from all slack to get some points 119 simplex>allSlackBasis(); 120 simplex>primal(); 121 122 // rlh: Mahdi's code 123 //  124 // Get the problem information 125 126 //  Get the number of cols and rows 127 int numCols = clpSolver>getNumCols(); 128 int numRows = clpSolver>getNumRows(); 129 130 //  Get the right hand side of the rows 131 const double * rhs = clpSolver>getRightHandSide(); 132 133 // Find the integer variables 134 // rlh: consider using columnType 135 // One if not continuous (i.e., bin or gen int) 136 bool * varClassInt = new bool[numCols]; 137 for(int i=0; i<numCols; i++) 138 { 139 if(clpSolver>isContinuous(i)) 140 varClassInt[i] = 0; 141 else 142 varClassInt[i] = 1; 143 } 144 145 // Get the rows sense 146 const char * rowSense; 147 rowSense = clpSolver>getRowSense(); 148 149 // Get the objective coefficients 150 const double *objCoefficients = clpSolver>getObjCoefficients(); 151 double *originalObjCoeff = new double [numCols]; 152 for(int i=0; i<numCols;i++) 153 originalObjCoeff[i] = objCoefficients[i]; 154 155 // Get the matrix of the problem 156 double ** matrix = new double * [numRows]; 157 for(int i = 0; i < numRows; i++) 158 { 159 matrix[i] = new double[numCols]; 160 for(int j = 0; j < numCols; j++) 161 matrix[i][j] = 0; 162 } 163 164 const CoinPackedMatrix* matrixByRow = clpSolver>getMatrixByRow(); 165 const double * matrixElements = matrixByRow>getElements(); 166 const int * matrixIndices = matrixByRow>getIndices(); 167 const int * matrixStarts = matrixByRow>getVectorStarts(); 168 for(int j = 0; j < numRows; j++) 169 { 170 for(int i = matrixStarts[j]; i < matrixStarts[j+1]; i++) 171 { 172 matrix[j][matrixIndices[i]] = matrixElements[i]; 173 } 174 } 175 176 double * newObj = new double [numCols]; 177 srand ( time(NULL) +1); 178 int randNum; 179 180 // Shuffle the rows: 181 // Put the rows in a random order 182 // so that the optimal solution is a different corner point than the 183 // starting point. 184 int * index = new int [numRows]; 185 for(int i=0; i<numRows; i++) 186 index[i]=i; 187 for(int i=0; i<numRows; i++) 188 { 189 int temp = index[i]; 190 int randNumTemp = i+(rand()%(numRowsi)); 191 index[i]=index[randNumTemp]; 192 index[randNumTemp]=temp; 193 } 194 195 // Start finding corner points by iteratively doing the following: 196 //  find randomly tilted objective 197 //  solve 198 for(int i=0; i<numRows; i++) 199 { 200 // TODO: that 10,000 could be a param in the member data 201 if(solutions.numberSolutions > 10000) 202 break; 203 randNum = rand()%10 + 1; 204 randNum = fmod(randNum, 2); 205 for(int j=0; j<numCols; j++) 206 { 207 // for row i and column j vary the coefficient "a bit" 208 if(randNum == 1) 209 // if the element is zero, then round the coefficient down to 0.1 210 if(fabs(matrix[index[i]][j]) < 1e6) 211 newObj[j] = 0.1; 212 else 213 // if the element is nonzero, then increase it "a bit" 214 newObj[j] = matrix[index[i]][j] * 1.1; 215 else 216 // if randnum is 2, then 217 // if the element is zero, then round the coefficient down 218 // to NEGATIVE 0.1 219 if(fabs(matrix[index[i]][j]) < 1e6) 220 newObj[j] = 0.1; 221 else 222 // if the element is nonzero, then DEcrease it "a bit" 223 newObj[j] = matrix[index[i]][j] * 0.9; 224 } 225 // Use the new "tilted" objective 226 clpSolver>setObjective(newObj); 227 228 // Based on the row sense, we decide whether to max or min 229 if(rowSense[i] == 'L') 230 clpSolver>setObjSense(1); 231 else 232 clpSolver>setObjSense(1); 233 234 // Solve with primal simplex 235 clpSolver>getModelPtr()>primal(1); 236 printf(" %d\n", i); 237 } 238 // Iteratively do this process until... 239 // either you reach the max number of corner points (aka 10K) 240 // or all the rows have been used as an objective. 241 242 // Look at solutions 243 int numberSolutions = solutions.numberSolutions; 244 const char * integerInfo = simplex>integerInformation(); 245 const double * columnLower = simplex>columnLower(); 246 const double * columnUpper = simplex>columnUpper(); 247 printf("there are %d solutions\n",numberSolutions); 248 249 // Up to here we have all the corner points 250 // Now we need to do the random walks and roundings 251 252 double ** cornerPoints = new double * [numberSolutions]; 253 for(int j=0;j<numberSolutions;j++) 254 cornerPoints[j] = solutions.solution[j]; 255 256 bool feasibility = 1; 257 double bestObj = 1e30; 258 vector< vector <double> > feasibles; 259 int numFeasibles = 0; 260 261 // Check the feasibility of the corner points 262 int numCornerPoints = numberSolutions; 263 264 rhs = clpSolver>getRightHandSide(); 265 rowSense = clpSolver>getRowSense(); 266 267 for(int i=0; i<numCornerPoints; i++) 268 { 269 //get the objective value for this this point 270 double objValue = 0; 271 for(int k=0; k<numCols; k++) 272 objValue += cornerPoints[i][k] * originalObjCoeff[k]; 273 274 if(objValue < bestObj) 275 { 276 feasibility = 1; 277 for(int j=0; j<numCols; j++) 278 { 279 if(varClassInt[j]) 280 { 281 double closest=floor(cornerPoints[i][j]+0.5); 282 if(fabs(cornerPoints[i][j]  closest)>1e6) 283 { 284 feasibility = 0; 285 break; 286 } 287 } 288 } 289 if(feasibility) 290 { 291 for(int irow = 0; irow < numRows; irow++) 292 { 293 double lhs = 0; 294 for(int j = 0; j <numCols; j++) 295 { 296 lhs += matrix[irow][j] * cornerPoints[i][j]; 297 } 298 if(rowSense[irow] == 'L' && lhs > rhs[irow] + 1e6) 299 { 300 feasibility = 0; 301 break; 302 } 303 if(rowSense[irow] == 'G' && lhs < rhs[irow]  1e6) 304 { 305 feasibility = 0; 306 break; 307 } 308 if(rowSense[irow] == 'E' && (lhs  rhs[irow] > 1e6  lhs  rhs[irow] < 1e6)) 309 { 310 feasibility = 0; 311 break; 312 } 313 } 314 } 315 316 if(feasibility) 317 { 318 numFeasibles++; 319 feasibles.push_back(vector <double> (numCols)); 320 for(int k=0; k<numCols; k++) 321 feasibles[numFeasibles1][k] = cornerPoints[i][k]; 322 printf("obj: %f\n", objValue); 323 if(objValue < bestObj) 324 bestObj = objValue; 325 } 326 } 327 } 328 int numFeasibleCorners; 329 numFeasibleCorners = numFeasibles; 330 //find the center of gravity of the corner points as the first random point 331 double * rp = new double[numCols]; 332 for(int i=0; i<numCols; i++) 333 { 334 rp[i] = 0; 335 for(int j=0; j < numCornerPoints; j++) 336 { 337 rp[i] += cornerPoints[j][i]; 338 } 339 rp[i] = rp[i]/numCornerPoints; 340 } 341 342 // 343 //main loop: 344 // generate the next random point 345 // round the random point 346 // check the feasibility of the random point 347 // 348 349 srand ( time(NULL) +1); 350 int numRandomPoints = 0; 351 while(numRandomPoints < 50000) 352 { 353 numRandomPoints++; 354 //generate the next random point 355 int randomIndex = rand()%numCornerPoints; 356 double random = rand(); 357 random = random/RAND_MAX; 358 for(int i=0; i<numCols; i++) 359 { 360 rp[i] = (random * (cornerPoints[randomIndex][i]  rp[i])) + rp[i]; 361 } 362 363 //CRISP ROUNDING 364 //round the random point just generated 365 double * roundRp = new double[numCols]; 366 for(int i=0; i<numCols; i++) 367 { 368 roundRp[i] = rp[i]; 369 if(varClassInt[i]) 370 { 371 if(rp[i]>=0) 372 { 373 if(fmod(rp[i],1) > 0.5) 374 roundRp[i] = floor(rp[i]) + 1; 375 else 376 roundRp[i] = floor(rp[i]); 377 } 378 else 379 { 380 if(abs(fmod(rp[i],1)) > 0.5) 381 roundRp[i] = floor(rp[i]); 382 else 383 roundRp[i] = floor(rp[i])+1; 384 385 } 386 } 387 } 388 389 390 //SOFT ROUNDING 391 // Look at original files for the "how to" on soft rounding 392 393 394 //check the feasibility of the rounded random point 395 // Check the feasibility 396 // Get the rows sense 397 rowSense = clpSolver>getRowSense(); 398 rhs = clpSolver>getRightHandSide(); 399 400 //get the objective value for this feasible point 401 double objValue = 0; 402 for(int i=0; i<numCols; i++) 403 objValue += roundRp[i] * originalObjCoeff[i]; 404 405 if(objValue<bestObj) 406 { 407 feasibility = 1; 408 for(int i = 0; i < numRows; i++) 409 { 410 double lhs = 0; 411 for(int j = 0; j <numCols; j++) 412 { 413 lhs += matrix[i][j] * roundRp[j]; 414 } 415 if(rowSense[i] == 'L' && lhs > rhs[i] + 1e6) 416 { 417 feasibility = 0; 418 break; 419 } 420 if(rowSense[i] == 'G' && lhs < rhs[i]  1e6) 421 { 422 feasibility = 0; 423 break; 424 } 425 if(rowSense[i] == 'E' && (lhs  rhs[i] > 1e6  lhs  rhs[i] < 1e6)) 426 { 427 feasibility = 0; 428 break; 429 } 430 } 431 if(feasibility) 432 { 433 printf("Feasible Found!!\n"); 434 printf("%.2f\n", CoinCpuTime()start); 435 numFeasibles++; 436 feasibles.push_back(vector <double> (numCols)); 437 for(int i=0; i<numCols; i++) 438 feasibles[numFeasibles1][i] = roundRp[i]; 439 printf("obj: %f\n", objValue); 440 if(objValue < bestObj) 441 bestObj = objValue; 442 } 443 } 444 } 445 printf("Number of Feasible Corners: %d\n", numFeasibleCorners); 446 printf("Number of Feasibles Found: %d\n", numFeasibles); 447 if(numFeasibles > 0) 448 printf("Best Objective: %f\n", bestObj); 449 printf("time: %.2f\n",CoinCpuTime()start); 450 451 if (numFeasibles == 0) return 0; 452 453 // We found something better 454 solutionValue = bestObj; 455 for (int k = 0; k<numCols; k++){ 456 betterSolution[k] = feasibles[numFeasibles1][k]; 457 } 458 std::cout << "See you soon! You're leaving the Randomized Rounding Heuristic" << std::endl; 459 return 1; 460 86 461 } 87 462 // update model
Note: See TracChangeset
for help on using the changeset viewer.