Changeset 862 for trunk/Cbc/src/CbcModel.cpp
 Timestamp:
 Jan 22, 2008 4:23:57 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Cbc/src/CbcModel.cpp
r857 r862 1194 1194 } 1195 1195 } 1196 1196 1197 solverCharacteristics_>setSolver(solver_); 1197 1198 // Set so we can tell we are in initial phase in resolve … … 2517 2518 double heurValue = getCutoff() ; 2518 2519 int iHeur ; 2519 for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++) 2520 {double saveValue = heurValue ;2520 for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++) { 2521 double saveValue = heurValue ; 2521 2522 int ifSol = heuristic_[iHeur]>solution(heurValue,newSolution) ; 2522 2523 if (ifSol > 0) { … … 2524 2525 found = iHeur ; 2525 2526 incrementUsed(newSolution); 2526 } 2527 else 2528 if (ifSol < 0) // just returning an estimate 2529 { estValue = CoinMin(heurValue,estValue) ; 2530 heurValue = saveValue ; } } 2531 if (found >= 0) { 2532 lastHeuristic_ = heuristic_[found]; 2533 setBestSolution(CBC_ROUNDING,heurValue,newSolution) ; 2534 } 2527 lastHeuristic_ = heuristic_[found]; 2528 setBestSolution(CBC_ROUNDING,heurValue,newSolution) ; 2529 } else if (ifSol < 0) { 2530 // just returning an estimate 2531 estValue = CoinMin(heurValue,estValue) ; 2532 heurValue = saveValue ; 2533 } 2534 } 2535 2535 delete [] newSolution ; 2536 2536 newNode>setGuessedObjectiveValue(estValue) ; … … 5108 5108 numberFixed++ ; } } } 5109 5109 numberDJFixed_ += numberFixed; 5110 5111 5110 return numberFixed; } 5112 5111 … … 5636 5635 printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n", 5637 5636 i,generator_[i]>cutGeneratorName(),knumberRowCutsBefore); 5637 const double *lower = getColLower() ; 5638 const double *upper = getColUpper() ; 5639 int numberColumns = solver_>getNumCols(); 5640 for (int i=0;i<numberColumns;i++) 5641 printf("%d bounds %g,%g\n",i,lower[i],upper[i]); 5638 5642 abort(); 5639 5643 #endif … … 5897 5901 printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n", 5898 5902 i,generator_[i]>cutGeneratorName(),knumberRowCutsBefore); 5903 const double *lower = getColLower() ; 5904 const double *upper = getColUpper() ; 5905 int numberColumns = solver_>getNumCols(); 5906 for (int i=0;i<numberColumns;i++) 5907 printf("%d bounds %g,%g\n",i,lower[i],upper[i]); 5899 5908 abort(); 5900 5909 #endif … … 5990 5999 found = i ; 5991 6000 incrementUsed(newSolution); 6001 lastHeuristic_ = heuristic_[found]; 6002 setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ; 5992 6003 } else if (ifSol<0) { 5993 6004 heuristicValue = saveValue ; … … 6000 6011 if (found >= 0) { 6001 6012 phase_=4; 6002 incrementUsed(newSolution);6003 lastHeuristic_ = heuristic_[found];6004 setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;6005 6013 CbcTreeLocal * tree 6006 6014 = dynamic_cast<CbcTreeLocal *> (tree_); … … 6354 6362 because a generator indicated we're infeasible. 6355 6363 */ 6356 if (feasible && solver_>isProvenOptimal()) 6357 reducedCostFix() ;6364 if (feasible && solver_>isProvenOptimal()) 6365 reducedCostFix() ; 6358 6366 // If at root node do heuristics 6359 6367 if (!numberNodes_) { … … 6393 6401 found = i ; 6394 6402 incrementUsed(newSolution); 6403 lastHeuristic_ = heuristic_[found]; 6404 setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ; 6395 6405 } else { 6396 6406 heuristicValue = saveValue ; … … 6400 6410 if (found >= 0) { 6401 6411 phase_=4; 6402 incrementUsed(newSolution);6403 lastHeuristic_ = heuristic_[found];6404 setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;6405 6412 } 6406 6413 delete [] newSolution ; … … 7022 7029 if (feasible) 7023 7030 { 7024 resolve(solver_) ; 7025 numberIterations_ += solver_>getIterationCount() ; 7026 feasible = (solver_>isProvenOptimal() && 7027 !solver_>isDualObjectiveLimitReached()) ; 7031 bool onOptimalPath=false; 7032 if ((specialOptions_&1)!=0) { 7033 const OsiRowCutDebugger *debugger = solver_>getRowCutDebugger() ; 7034 if (debugger) { 7035 onOptimalPath=true; 7036 printf("On optimal path d\n") ; 7037 } 7038 } 7039 int nTightened=0; 7040 #if 1 7041 if (!currentNode_(currentNode_>depth()&2)!=0) 7042 nTightened=tightenBounds(); 7043 if (nTightened) { 7044 //printf("%d bounds tightened\n",nTightened); 7045 if ((specialOptions_&1)!=0&&onOptimalPath) { 7046 const OsiRowCutDebugger *debugger = solver_>getRowCutDebugger() ; 7047 if (!debugger) { 7048 // tighten did something??? 7049 solver_>writeMps("infeas4"); 7050 printf("Not on optimalpath aaaa\n"); 7051 abort(); 7052 } 7053 } 7054 } 7055 #endif 7056 if (nTightened>=0) { 7057 resolve(solver_) ; 7058 numberIterations_ += solver_>getIterationCount() ; 7059 feasible = (solver_>isProvenOptimal() && 7060 !solver_>isDualObjectiveLimitReached()) ; 7061 if ((specialOptions_&1)!=0&&onOptimalPath) { 7062 if (!solver_>getRowCutDebugger()) { 7063 // tighten did something??? 7064 solver_>writeMps("infeas4"); 7065 assert (solver_>getRowCutDebugger()) ; 7066 printf("Not on optimalpath e\n"); 7067 abort(); 7068 } 7069 } 7070 } else { 7071 feasible=false; 7072 } 7028 7073 } 7029 7074 if (0&&feasible) { … … 9147 9192 setCutoff(saveCutoff); 9148 9193 return true; 9194 } 9195 // Tighten bounds  lightweight 9196 int 9197 CbcModel::tightenBounds() 9198 { 9199 //CoinPackedMatrix matrixByRow(*solver_>getMatrixByRow()); 9200 int numberRows = solver_>getNumRows(); 9201 int numberColumns = solver_>getNumCols(); 9202 9203 int iRow,iColumn; 9204 9205 // Row copy 9206 //const double * elementByRow = matrixByRow.getElements(); 9207 //const int * column = matrixByRow.getIndices(); 9208 //const CoinBigIndex * rowStart = matrixByRow.getVectorStarts(); 9209 //const int * rowLength = matrixByRow.getVectorLengths(); 9210 9211 const double * columnUpper = solver_>getColUpper(); 9212 const double * columnLower = solver_>getColLower(); 9213 const double * rowUpper = solver_>getRowUpper(); 9214 const double * rowLower = solver_>getRowLower(); 9215 9216 // Column copy of matrix 9217 const double * element = solver_>getMatrixByCol()>getElements(); 9218 const int * row = solver_>getMatrixByCol()>getIndices(); 9219 const CoinBigIndex * columnStart = solver_>getMatrixByCol()>getVectorStarts(); 9220 const int * columnLength = solver_>getMatrixByCol()>getVectorLengths(); 9221 const double *objective = solver_>getObjCoefficients() ; 9222 double direction = solver_>getObjSense(); 9223 double * down = new double [numberRows]; 9224 double * up = new double [numberRows]; 9225 double * sum = new double [numberRows]; 9226 int * type = new int [numberRows]; 9227 CoinZeroN(down,numberRows); 9228 CoinZeroN(up,numberRows); 9229 CoinZeroN(sum,numberRows); 9230 CoinZeroN(type,numberRows); 9231 double infinity = solver_>getInfinity(); 9232 for (iColumn=0;iColumn<numberColumns;iColumn++) { 9233 CoinBigIndex start = columnStart[iColumn]; 9234 CoinBigIndex end = start + columnLength[iColumn]; 9235 double lower = columnLower[iColumn]; 9236 double upper = columnUpper[iColumn]; 9237 if (lower==upper) { 9238 for (CoinBigIndex j=start;j<end;j++) { 9239 int iRow = row[j]; 9240 double value = element[j]; 9241 sum[iRow]+=2.0*fabs(value*lower); 9242 if ((type[iRow]&1)==0) 9243 down[iRow] += value*lower; 9244 if ((type[iRow]&2)==0) 9245 up[iRow] += value*lower; 9246 } 9247 } else { 9248 for (CoinBigIndex j=start;j<end;j++) { 9249 int iRow = row[j]; 9250 double value = element[j]; 9251 if (value>0.0) { 9252 if ((type[iRow]&1)==0) { 9253 if (lower!=infinity) { 9254 down[iRow] += value*lower; 9255 sum[iRow]+=fabs(value*lower); 9256 } else { 9257 type[iRow] = 1; 9258 } 9259 } 9260 if ((type[iRow]&2)==0) { 9261 if (upper!=infinity) { 9262 up[iRow] += value*upper; 9263 sum[iRow]+=fabs(value*upper); 9264 } else { 9265 type[iRow] = 2; 9266 } 9267 } 9268 } else { 9269 if ((type[iRow]&1)==0) { 9270 if (upper!=infinity) { 9271 down[iRow] += value*upper; 9272 sum[iRow]+=fabs(value*upper); 9273 } else { 9274 type[iRow] = 1; 9275 } 9276 } 9277 if ((type[iRow]&2)==0) { 9278 if (lower!=infinity) { 9279 up[iRow] += value*lower; 9280 sum[iRow]+=fabs(value*lower); 9281 } else { 9282 type[iRow] = 2; 9283 } 9284 } 9285 } 9286 } 9287 } 9288 } 9289 int nTightened=0; 9290 double tolerance = 1.0e6; 9291 for (iRow=0;iRow<numberRows;iRow++) { 9292 if ((type[iRow]&1)!=0) 9293 down[iRow]=infinity; 9294 if (down[iRow]>rowUpper[iRow]) { 9295 if (down[iRow]>rowUpper[iRow]+tolerance+1.0e8*sum[iRow]) { 9296 // infeasible 9297 #ifdef COIN_DEVELOP 9298 printf("infeasible on row %d\n",iRow); 9299 #endif 9300 nTightened=1; 9301 break; 9302 } else { 9303 down[iRow]=rowUpper[iRow]; 9304 } 9305 } 9306 if ((type[iRow]&2)!=0) 9307 up[iRow]=infinity; 9308 if (up[iRow]<rowLower[iRow]) { 9309 if (up[iRow]<rowLower[iRow]tolerance1.0e8*sum[iRow]) { 9310 // infeasible 9311 #ifdef COIN_DEVELOP 9312 printf("infeasible on row %d\n",iRow); 9313 #endif 9314 nTightened=1; 9315 break; 9316 } else { 9317 up[iRow]=rowLower[iRow]; 9318 } 9319 } 9320 } 9321 if (nTightened) 9322 numberColumns=0; // so will skip 9323 for (iColumn=0;iColumn<numberColumns;iColumn++) { 9324 CoinBigIndex start = columnStart[iColumn]; 9325 CoinBigIndex end = start + columnLength[iColumn]; 9326 double lower = columnLower[iColumn]; 9327 double upper = columnUpper[iColumn]; 9328 double gap = upperlower; 9329 if (!gap) 9330 continue; 9331 int canGo=0; 9332 if (integerInfo_[iColumn]) { 9333 if (lower!=floor(lower+0.5)) { 9334 #ifdef COIN_DEVELOP 9335 printf("increasing lower bound on %d from %g to %g\n",iColumn, 9336 lower,ceil(lower)); 9337 #endif 9338 lower=ceil(lower); 9339 gap=upperlower; 9340 solver_>setColLower(iColumn,lower); 9341 } 9342 if (upper!=floor(upper+0.5)) { 9343 #ifdef COIN_DEVELOP 9344 printf("decreasing upper bound on %d from %g to %g\n",iColumn, 9345 upper,floor(upper)); 9346 #endif 9347 upper=floor(upper); 9348 gap=upperlower; 9349 solver_>setColUpper(iColumn,upper); 9350 } 9351 double newLower=lower; 9352 double newUpper=upper; 9353 for (CoinBigIndex j=start;j<end;j++) { 9354 int iRow = row[j]; 9355 double value = element[j]; 9356 if (value>0.0) { 9357 if ((type[iRow]&1)==0) { 9358 // has to be at most something 9359 if (down[iRow] + value*gap > rowUpper[iRow]+tolerance) { 9360 double newGap = (rowUpper[iRow]down[iRow])/value; 9361 // adjust 9362 newGap += 1.0e10*sum[iRow]; 9363 newGap = floor(newGap); 9364 if (lower+newGap<newUpper) 9365 newUpper=lower+newGap; 9366 } 9367 } 9368 if (down[iRow]<rowLower[iRow]) 9369 canGo =1; // can't go down without affecting result 9370 if ((type[iRow]&2)==0) { 9371 // has to be at least something 9372 if (up[iRow]  value*gap < rowLower[iRow]tolerance) { 9373 double newGap = (up[iRow]rowLower[iRow])/value; 9374 // adjust 9375 newGap += 1.0e10*sum[iRow]; 9376 newGap = floor(newGap); 9377 if (uppernewGap>newLower) 9378 newLower=uppernewGap; 9379 } 9380 } 9381 if (up[iRow]>rowUpper[iRow]) 9382 canGo =2; // can't go up without affecting result 9383 } else { 9384 if ((type[iRow]&1)==0) { 9385 // has to be at least something 9386 if (down[iRow]  value*gap > rowUpper[iRow]+tolerance) { 9387 double newGap = (rowUpper[iRow]down[iRow])/value; 9388 // adjust 9389 newGap += 1.0e10*sum[iRow]; 9390 newGap = floor(newGap); 9391 if (uppernewGap>newLower) 9392 newLower=uppernewGap; 9393 } 9394 } 9395 if (up[iRow]>rowUpper[iRow]) 9396 canGo =1; // can't go down without affecting result 9397 if ((type[iRow]&2)==0) { 9398 // has to be at most something 9399 if (up[iRow] + value*gap < rowLower[iRow]tolerance) { 9400 double newGap = (up[iRow]rowLower[iRow])/value; 9401 // adjust 9402 newGap += 1.0e10*sum[iRow]; 9403 newGap = floor(newGap); 9404 if (lower+newGap<newUpper) 9405 newUpper=lower+newGap; 9406 } 9407 } 9408 if (down[iRow]<rowLower[iRow]) 9409 canGo =2; // can't go up without affecting result 9410 } 9411 } 9412 if (newUpper<uppernewLower>lower) { 9413 nTightened++; 9414 if (newLower>newUpper) { 9415 // infeasible 9416 #if COIN_DEVELOP>1 9417 printf("infeasible on column %d\n",iColumn); 9418 #endif 9419 nTightened=1; 9420 break; 9421 } else { 9422 solver_>setColLower(iColumn,newLower); 9423 solver_>setColUpper(iColumn,newUpper); 9424 } 9425 for (CoinBigIndex j=start;j<end;j++) { 9426 int iRow = row[j]; 9427 double value = element[j]; 9428 if (value>0.0) { 9429 if ((type[iRow]&1)==0) { 9430 down[iRow] += value*(newLowerlower); 9431 } 9432 if ((type[iRow]&2)==0) { 9433 up[iRow] += value*(newUpperupper); 9434 } 9435 } else { 9436 if ((type[iRow]&1)==0) { 9437 down[iRow] += value*(newUpperupper); 9438 } 9439 if ((type[iRow]&2)==0) { 9440 up[iRow] += value*(newLowerlower); 9441 } 9442 } 9443 } 9444 } else { 9445 if (canGo!=3) { 9446 double objValue = direction*objective[iColumn]; 9447 if (objValue>=0.0&&(canGo&1)==0) { 9448 #if COIN_DEVELOP>1 9449 printf("dual fix down on column %d\n",iColumn); 9450 #endif 9451 nTightened++;; 9452 solver_>setColUpper(iColumn,lower); 9453 } else if (objValue<=0.0&&(canGo&2)==0) { 9454 #if COIN_DEVELOP>1 9455 printf("dual fix up on column %d\n",iColumn); 9456 #endif 9457 nTightened++;; 9458 solver_>setColLower(iColumn,upper); 9459 } 9460 } 9461 } 9462 } else { 9463 // just do dual tests 9464 for (CoinBigIndex j=start;j<end;j++) { 9465 int iRow = row[j]; 9466 double value = element[j]; 9467 if (value>0.0) { 9468 if (down[iRow]<rowLower[iRow]) 9469 canGo =1; // can't go down without affecting result 9470 if (up[iRow]>rowUpper[iRow]) 9471 canGo =2; // can't go up without affecting result 9472 } else { 9473 if (up[iRow]>rowUpper[iRow]) 9474 canGo =1; // can't go down without affecting result 9475 if (down[iRow]<rowLower[iRow]) 9476 canGo =2; // can't go up without affecting result 9477 } 9478 } 9479 if (canGo!=3) { 9480 double objValue = direction*objective[iColumn]; 9481 if (objValue>=0.0&&(canGo&1)==0) { 9482 #if COIN_DEVELOP>1 9483 printf("dual fix down on continuous column %d\n",iColumn); 9484 #endif 9485 nTightened++;; 9486 solver_>setColUpper(iColumn,lower); 9487 } else if (objValue<=0.0&&(canGo&2)==0) { 9488 #if COIN_DEVELOP>1 9489 printf("dual fix up on continuoe column %d\n",iColumn); 9490 #endif 9491 nTightened++;; 9492 solver_>setColLower(iColumn,upper); 9493 } 9494 } 9495 } 9496 } 9497 delete [] type; 9498 delete [] down; 9499 delete [] up; 9500 delete [] sum; 9501 return nTightened; 9149 9502 } 9150 9503 /* … … 10868 11221 numberSolutions_++; 10869 11222 numberHeuristicSolutions_++; 11223 lastHeuristic_ = heuristic_[i]; 11224 setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ; 10870 11225 } else { 10871 11226 heuristicValue = saveValue ; … … 10878 11233 */ 10879 11234 if (found >= 0) { 10880 // For compiler error on terra cluster!10881 if (found<numberHeuristics_)10882 lastHeuristic_ = heuristic_[found];10883 else10884 lastHeuristic_ = heuristic_[0];10885 setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;10886 11235 CbcTreeLocal * tree 10887 11236 = dynamic_cast<CbcTreeLocal *> (tree_); … … 12795 13144 CbcIntegerBranchingObject * objectI = dynamic_cast<CbcIntegerBranchingObject *> (obj); 12796 13145 //const OsiObject * object2 = obj>orig 13146 #ifndef NDEBUG 12797 13147 const CbcSimpleInteger * object2 = dynamic_cast<const CbcSimpleInteger *> (objectI>object()); 12798 13148 assert (object2); 12799 13149 assert (iColumn == object2>columnNumber()); 13150 #endif 12800 13151 double bounds[2]; 12801 13152 bounds[0]=lower; … … 12807 13158 nNode; 12808 13159 walkback_[nNode]>applyBounds(iColumn,lower,upper,force); 13160 #if 0 12809 13161 CbcNode * nodeLook = walkback_[nNode]>mutableOwner(); 12810 if (nodeLook &&0) {13162 if (nodeLook) { 12811 13163 const OsiBranchingObject * obj = nodeLook>branchingObject(); 12812 13164 const CbcIntegerBranchingObject * objectI = dynamic_cast<const CbcIntegerBranchingObject *> (obj); … … 12817 13169 assert (iColumn!=iColumn2); 12818 13170 } 13171 #endif 12819 13172 } 12820 13173 }
Note: See TracChangeset
for help on using the changeset viewer.