[1271] | 1 | /* $Id: CbcHeuristicFPump.cpp 1573 2011-01-05 01:12:36Z tkr $ */ |
---|
[175] | 2 | // Copyright (C) 2004, International Business Machines |
---|
| 3 | // Corporation and others. All Rights Reserved. |
---|
[1573] | 4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
---|
| 5 | |
---|
[175] | 6 | #if defined(_MSC_VER) |
---|
| 7 | // Turn off compiler warning about long names |
---|
| 8 | # pragma warning(disable:4786) |
---|
| 9 | #endif |
---|
| 10 | #include <cassert> |
---|
[904] | 11 | #include <cstdlib> |
---|
[175] | 12 | #include <cmath> |
---|
| 13 | #include <cfloat> |
---|
| 14 | |
---|
| 15 | #include "OsiSolverInterface.hpp" |
---|
| 16 | #include "CbcModel.hpp" |
---|
[1088] | 17 | #ifdef COIN_HAS_CLP |
---|
| 18 | #include "OsiClpSolverInterface.hpp" |
---|
| 19 | #endif |
---|
[175] | 20 | #include "CbcMessage.hpp" |
---|
| 21 | #include "CbcHeuristicFPump.hpp" |
---|
| 22 | #include "CbcBranchActual.hpp" |
---|
[1271] | 23 | #include "CbcBranchDynamic.hpp" |
---|
[175] | 24 | #include "CoinHelperFunctions.hpp" |
---|
[738] | 25 | #include "CoinWarmStartBasis.hpp" |
---|
[175] | 26 | #include "CoinTime.hpp" |
---|
[747] | 27 | #include "CbcEventHandler.hpp" |
---|
[175] | 28 | |
---|
| 29 | |
---|
| 30 | // Default Constructor |
---|
[1286] | 31 | CbcHeuristicFPump::CbcHeuristicFPump() |
---|
| 32 | : CbcHeuristic(), |
---|
| 33 | startTime_(0.0), |
---|
| 34 | maximumTime_(0.0), |
---|
| 35 | fakeCutoff_(COIN_DBL_MAX), |
---|
| 36 | absoluteIncrement_(0.0), |
---|
| 37 | relativeIncrement_(0.0), |
---|
| 38 | defaultRounding_(0.49999), |
---|
| 39 | initialWeight_(0.0), |
---|
| 40 | weightFactor_(0.1), |
---|
| 41 | artificialCost_(COIN_DBL_MAX), |
---|
| 42 | iterationRatio_(0.0), |
---|
| 43 | reducedCostMultiplier_(1.0), |
---|
| 44 | maximumPasses_(100), |
---|
| 45 | maximumRetries_(1), |
---|
| 46 | accumulate_(0), |
---|
| 47 | fixOnReducedCosts_(1), |
---|
| 48 | roundExpensive_(false) |
---|
[175] | 49 | { |
---|
[1286] | 50 | setWhen(1); |
---|
[175] | 51 | } |
---|
| 52 | |
---|
| 53 | // Constructor from model |
---|
| 54 | CbcHeuristicFPump::CbcHeuristicFPump(CbcModel & model, |
---|
[1286] | 55 | double downValue, bool roundExpensive) |
---|
| 56 | : CbcHeuristic(model), |
---|
| 57 | startTime_(0.0), |
---|
| 58 | maximumTime_(0.0), |
---|
| 59 | fakeCutoff_(COIN_DBL_MAX), |
---|
| 60 | absoluteIncrement_(0.0), |
---|
| 61 | relativeIncrement_(0.0), |
---|
| 62 | defaultRounding_(downValue), |
---|
| 63 | initialWeight_(0.0), |
---|
| 64 | weightFactor_(0.1), |
---|
| 65 | artificialCost_(COIN_DBL_MAX), |
---|
| 66 | iterationRatio_(0.0), |
---|
| 67 | reducedCostMultiplier_(1.0), |
---|
| 68 | maximumPasses_(100), |
---|
| 69 | maximumRetries_(1), |
---|
| 70 | accumulate_(0), |
---|
| 71 | fixOnReducedCosts_(1), |
---|
| 72 | roundExpensive_(roundExpensive) |
---|
[175] | 73 | { |
---|
[1286] | 74 | setWhen(1); |
---|
[175] | 75 | } |
---|
| 76 | |
---|
[1286] | 77 | // Destructor |
---|
[175] | 78 | CbcHeuristicFPump::~CbcHeuristicFPump () |
---|
| 79 | { |
---|
| 80 | } |
---|
| 81 | |
---|
| 82 | // Clone |
---|
| 83 | CbcHeuristic * |
---|
| 84 | CbcHeuristicFPump::clone() const |
---|
| 85 | { |
---|
[1286] | 86 | return new CbcHeuristicFPump(*this); |
---|
[175] | 87 | } |
---|
[356] | 88 | // Create C++ lines to get to current state |
---|
[1286] | 89 | void |
---|
| 90 | CbcHeuristicFPump::generateCpp( FILE * fp) |
---|
[356] | 91 | { |
---|
[1286] | 92 | CbcHeuristicFPump other; |
---|
| 93 | fprintf(fp, "0#include \"CbcHeuristicFPump.hpp\"\n"); |
---|
| 94 | fprintf(fp, "3 CbcHeuristicFPump heuristicFPump(*cbcModel);\n"); |
---|
| 95 | CbcHeuristic::generateCpp(fp, "heuristicFPump"); |
---|
| 96 | if (maximumPasses_ != other.maximumPasses_) |
---|
| 97 | fprintf(fp, "3 heuristicFPump.setMaximumPasses(%d);\n", maximumPasses_); |
---|
| 98 | else |
---|
| 99 | fprintf(fp, "4 heuristicFPump.setMaximumPasses(%d);\n", maximumPasses_); |
---|
| 100 | if (maximumRetries_ != other.maximumRetries_) |
---|
| 101 | fprintf(fp, "3 heuristicFPump.setMaximumRetries(%d);\n", maximumRetries_); |
---|
| 102 | else |
---|
| 103 | fprintf(fp, "4 heuristicFPump.setMaximumRetries(%d);\n", maximumRetries_); |
---|
| 104 | if (accumulate_ != other.accumulate_) |
---|
| 105 | fprintf(fp, "3 heuristicFPump.setAccumulate(%d);\n", accumulate_); |
---|
| 106 | else |
---|
| 107 | fprintf(fp, "4 heuristicFPump.setAccumulate(%d);\n", accumulate_); |
---|
| 108 | if (fixOnReducedCosts_ != other.fixOnReducedCosts_) |
---|
| 109 | fprintf(fp, "3 heuristicFPump.setFixOnReducedCosts(%d);\n", fixOnReducedCosts_); |
---|
| 110 | else |
---|
| 111 | fprintf(fp, "4 heuristicFPump.setFixOnReducedCosts(%d);\n", fixOnReducedCosts_); |
---|
| 112 | if (maximumTime_ != other.maximumTime_) |
---|
| 113 | fprintf(fp, "3 heuristicFPump.setMaximumTime(%g);\n", maximumTime_); |
---|
| 114 | else |
---|
| 115 | fprintf(fp, "4 heuristicFPump.setMaximumTime(%g);\n", maximumTime_); |
---|
| 116 | if (fakeCutoff_ != other.fakeCutoff_) |
---|
| 117 | fprintf(fp, "3 heuristicFPump.setFakeCutoff(%g);\n", fakeCutoff_); |
---|
| 118 | else |
---|
| 119 | fprintf(fp, "4 heuristicFPump.setFakeCutoff(%g);\n", fakeCutoff_); |
---|
| 120 | if (absoluteIncrement_ != other.absoluteIncrement_) |
---|
| 121 | fprintf(fp, "3 heuristicFPump.setAbsoluteIncrement(%g);\n", absoluteIncrement_); |
---|
| 122 | else |
---|
| 123 | fprintf(fp, "4 heuristicFPump.setAbsoluteIncrement(%g);\n", absoluteIncrement_); |
---|
| 124 | if (relativeIncrement_ != other.relativeIncrement_) |
---|
| 125 | fprintf(fp, "3 heuristicFPump.setRelativeIncrement(%g);\n", relativeIncrement_); |
---|
| 126 | else |
---|
| 127 | fprintf(fp, "4 heuristicFPump.setRelativeIncrement(%g);\n", relativeIncrement_); |
---|
| 128 | if (defaultRounding_ != other.defaultRounding_) |
---|
| 129 | fprintf(fp, "3 heuristicFPump.setDefaultRounding(%g);\n", defaultRounding_); |
---|
| 130 | else |
---|
| 131 | fprintf(fp, "4 heuristicFPump.setDefaultRounding(%g);\n", defaultRounding_); |
---|
| 132 | if (initialWeight_ != other.initialWeight_) |
---|
| 133 | fprintf(fp, "3 heuristicFPump.setInitialWeight(%g);\n", initialWeight_); |
---|
| 134 | else |
---|
| 135 | fprintf(fp, "4 heuristicFPump.setInitialWeight(%g);\n", initialWeight_); |
---|
| 136 | if (weightFactor_ != other.weightFactor_) |
---|
| 137 | fprintf(fp, "3 heuristicFPump.setWeightFactor(%g);\n", weightFactor_); |
---|
| 138 | else |
---|
| 139 | fprintf(fp, "4 heuristicFPump.setWeightFactor(%g);\n", weightFactor_); |
---|
| 140 | if (artificialCost_ != other.artificialCost_) |
---|
| 141 | fprintf(fp, "3 heuristicFPump.setArtificialCost(%g);\n", artificialCost_); |
---|
| 142 | else |
---|
| 143 | fprintf(fp, "4 heuristicFPump.setArtificialCost(%g);\n", artificialCost_); |
---|
| 144 | if (iterationRatio_ != other.iterationRatio_) |
---|
| 145 | fprintf(fp, "3 heuristicFPump.setIterationRatio(%g);\n", iterationRatio_); |
---|
| 146 | else |
---|
| 147 | fprintf(fp, "4 heuristicFPump.setIterationRatio(%g);\n", iterationRatio_); |
---|
| 148 | if (reducedCostMultiplier_ != other.reducedCostMultiplier_) |
---|
| 149 | fprintf(fp, "3 heuristicFPump.setReducedCostMultiplier(%g);\n", reducedCostMultiplier_); |
---|
| 150 | else |
---|
| 151 | fprintf(fp, "4 heuristicFPump.setReducedCostMultiplier(%g);\n", reducedCostMultiplier_); |
---|
| 152 | fprintf(fp, "3 cbcModel->addHeuristic(&heuristicFPump);\n"); |
---|
[356] | 153 | } |
---|
[175] | 154 | |
---|
[1286] | 155 | // Copy constructor |
---|
[175] | 156 | CbcHeuristicFPump::CbcHeuristicFPump(const CbcHeuristicFPump & rhs) |
---|
[1286] | 157 | : |
---|
| 158 | CbcHeuristic(rhs), |
---|
| 159 | startTime_(rhs.startTime_), |
---|
| 160 | maximumTime_(rhs.maximumTime_), |
---|
| 161 | fakeCutoff_(rhs.fakeCutoff_), |
---|
| 162 | absoluteIncrement_(rhs.absoluteIncrement_), |
---|
| 163 | relativeIncrement_(rhs.relativeIncrement_), |
---|
| 164 | defaultRounding_(rhs.defaultRounding_), |
---|
| 165 | initialWeight_(rhs.initialWeight_), |
---|
| 166 | weightFactor_(rhs.weightFactor_), |
---|
| 167 | artificialCost_(rhs.artificialCost_), |
---|
| 168 | iterationRatio_(rhs.iterationRatio_), |
---|
| 169 | reducedCostMultiplier_(rhs.reducedCostMultiplier_), |
---|
| 170 | maximumPasses_(rhs.maximumPasses_), |
---|
| 171 | maximumRetries_(rhs.maximumRetries_), |
---|
| 172 | accumulate_(rhs.accumulate_), |
---|
| 173 | fixOnReducedCosts_(rhs.fixOnReducedCosts_), |
---|
| 174 | roundExpensive_(rhs.roundExpensive_) |
---|
[175] | 175 | { |
---|
| 176 | } |
---|
[640] | 177 | |
---|
[1286] | 178 | // Assignment operator |
---|
| 179 | CbcHeuristicFPump & |
---|
| 180 | CbcHeuristicFPump::operator=( const CbcHeuristicFPump & rhs) |
---|
[640] | 181 | { |
---|
[1286] | 182 | if (this != &rhs) { |
---|
| 183 | CbcHeuristic::operator=(rhs); |
---|
| 184 | startTime_ = rhs.startTime_; |
---|
| 185 | maximumTime_ = rhs.maximumTime_; |
---|
| 186 | fakeCutoff_ = rhs.fakeCutoff_; |
---|
| 187 | absoluteIncrement_ = rhs.absoluteIncrement_; |
---|
| 188 | relativeIncrement_ = rhs.relativeIncrement_; |
---|
| 189 | defaultRounding_ = rhs.defaultRounding_; |
---|
| 190 | initialWeight_ = rhs.initialWeight_; |
---|
| 191 | weightFactor_ = rhs.weightFactor_; |
---|
| 192 | artificialCost_ = rhs.artificialCost_; |
---|
| 193 | iterationRatio_ = rhs.iterationRatio_; |
---|
| 194 | reducedCostMultiplier_ = rhs.reducedCostMultiplier_; |
---|
| 195 | maximumPasses_ = rhs.maximumPasses_; |
---|
| 196 | maximumRetries_ = rhs.maximumRetries_; |
---|
| 197 | accumulate_ = rhs.accumulate_; |
---|
| 198 | fixOnReducedCosts_ = rhs.fixOnReducedCosts_; |
---|
| 199 | roundExpensive_ = rhs.roundExpensive_; |
---|
| 200 | } |
---|
| 201 | return *this; |
---|
[640] | 202 | } |
---|
| 203 | |
---|
[175] | 204 | // Resets stuff if model changes |
---|
[1286] | 205 | void |
---|
[1271] | 206 | CbcHeuristicFPump::resetModel(CbcModel * ) |
---|
[175] | 207 | { |
---|
| 208 | } |
---|
| 209 | |
---|
| 210 | /**************************BEGIN MAIN PROCEDURE ***********************************/ |
---|
| 211 | |
---|
| 212 | // See if feasibility pump will give better solution |
---|
| 213 | // Sets value of solution |
---|
| 214 | // Returns 1 if solution, 0 if not |
---|
| 215 | int |
---|
| 216 | CbcHeuristicFPump::solution(double & solutionValue, |
---|
[1286] | 217 | double * betterSolution) |
---|
[175] | 218 | { |
---|
[1286] | 219 | startTime_ = CoinCpuTime(); |
---|
| 220 | numCouldRun_++; |
---|
| 221 | double incomingObjective = solutionValue; |
---|
[698] | 222 | #define LEN_PRINT 200 |
---|
[1286] | 223 | char pumpPrint[LEN_PRINT]; |
---|
| 224 | pumpPrint[0] = '\0'; |
---|
[1364] | 225 | /* |
---|
| 226 | Decide if we want to run. Standard values for when are described in |
---|
| 227 | CbcHeuristic.hpp. If we're off, or running only at root and this isn't the |
---|
| 228 | root, bail out. |
---|
| 229 | |
---|
| 230 | The double test (against phase, then atRoot and passNumber) has a fair bit |
---|
| 231 | of redundancy, but the results will differ depending on whether we're |
---|
| 232 | actually at the root of the main search tree or at the root of a small tree |
---|
| 233 | (recursive call to branchAndBound). |
---|
| 234 | |
---|
| 235 | FPump also supports some exotic values (11 -- 15) for when, described in |
---|
| 236 | CbcHeuristicFPump.hpp. |
---|
| 237 | */ |
---|
[1286] | 238 | if (!when() || (when() == 1 && model_->phase() != 1)) |
---|
| 239 | return 0; // switched off |
---|
| 240 | // See if at root node |
---|
| 241 | bool atRoot = model_->getNodeCount() == 0; |
---|
| 242 | int passNumber = model_->getCurrentPassNumber(); |
---|
| 243 | // just do once |
---|
| 244 | if (!atRoot) |
---|
| 245 | return 0; |
---|
| 246 | int options = feasibilityPumpOptions_; |
---|
| 247 | if ((options % 1000000) > 0) { |
---|
| 248 | int kOption = options / 1000000; |
---|
| 249 | options = options % 1000000; |
---|
| 250 | /* |
---|
| 251 | Add 10 to do even if solution |
---|
| 252 | 1 - do after cuts |
---|
| 253 | 2 - do after cuts (not before) |
---|
| 254 | 3 - not used do after every cut round (and after cuts) |
---|
| 255 | k not used do after every (k-2)th round |
---|
| 256 | */ |
---|
| 257 | if (kOption < 10 && model_->getSolutionCount()) |
---|
| 258 | return 0; |
---|
| 259 | if (model_->getSolutionCount()) |
---|
| 260 | kOption = kOption % 10; |
---|
| 261 | bool good; |
---|
| 262 | if (kOption == 1) { |
---|
| 263 | good = (passNumber == 999999); |
---|
| 264 | } else if (kOption == 2) { |
---|
| 265 | good = (passNumber == 999999); |
---|
| 266 | passNumber = 2; // so won't run before |
---|
| 267 | //} else if (kOption==3) { |
---|
| 268 | //good = true; |
---|
| 269 | } else { |
---|
| 270 | //good = (((passNumber-1)%(kOption-2))==0); |
---|
| 271 | good = false; |
---|
| 272 | } |
---|
| 273 | if (passNumber != 1 && !good) |
---|
| 274 | return 0; |
---|
[1271] | 275 | } else { |
---|
[1286] | 276 | if (passNumber != 1) |
---|
| 277 | return 0; |
---|
[1271] | 278 | } |
---|
[1286] | 279 | // loop round doing repeated pumps |
---|
| 280 | double cutoff; |
---|
| 281 | model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff); |
---|
| 282 | double direction = model_->solver()->getObjSense(); |
---|
| 283 | cutoff *= direction; |
---|
| 284 | int numberBandBsolutions = 0; |
---|
| 285 | double firstCutoff = fabs(cutoff); |
---|
| 286 | cutoff = CoinMin(cutoff, solutionValue); |
---|
| 287 | // check plausible and space for rounded solution |
---|
| 288 | int numberColumns = model_->getNumCols(); |
---|
| 289 | int numberIntegers = model_->numberIntegers(); |
---|
| 290 | const int * integerVariableOrig = model_->integerVariable(); |
---|
| 291 | double iterationLimit = -1.0; |
---|
| 292 | //iterationRatio_=1.0; |
---|
| 293 | if (iterationRatio_ > 0.0) |
---|
| 294 | iterationLimit = (2 * model_->solver()->getNumRows() + 2 * numberColumns) * |
---|
| 295 | iterationRatio_; |
---|
| 296 | int totalNumberIterations = 0; |
---|
| 297 | int averageIterationsPerTry = -1; |
---|
| 298 | int numberIterationsLastPass = 0; |
---|
| 299 | // 1. initially check 0-1 |
---|
[1364] | 300 | /* |
---|
| 301 | I'm skeptical of the above comment, but it's likely accurate as the default. |
---|
| 302 | Bit 4 or bit 8 needs to be set in order to consider working with general |
---|
| 303 | integers. |
---|
| 304 | */ |
---|
[1286] | 305 | int i, j; |
---|
| 306 | int general = 0; |
---|
| 307 | int * integerVariable = new int[numberIntegers]; |
---|
| 308 | const double * lower = model_->solver()->getColLower(); |
---|
| 309 | const double * upper = model_->solver()->getColUpper(); |
---|
| 310 | bool doGeneral = (accumulate_ & 4) != 0; |
---|
| 311 | j = 0; |
---|
[1364] | 312 | /* |
---|
| 313 | Scan the objects, recording the columns and counting general integers. |
---|
| 314 | |
---|
| 315 | Seems like the NDEBUG tests could be made into an applicability test. If |
---|
| 316 | a scan of the objects reveals complex objects, just clean up and return |
---|
| 317 | failure. |
---|
| 318 | */ |
---|
[1286] | 319 | for (i = 0; i < numberIntegers; i++) { |
---|
| 320 | int iColumn = integerVariableOrig[i]; |
---|
[356] | 321 | #ifndef NDEBUG |
---|
[1286] | 322 | const OsiObject * object = model_->object(i); |
---|
| 323 | const CbcSimpleInteger * integerObject = |
---|
| 324 | dynamic_cast<const CbcSimpleInteger *> (object); |
---|
| 325 | const OsiSimpleInteger * integerObject2 = |
---|
| 326 | dynamic_cast<const OsiSimpleInteger *> (object); |
---|
| 327 | assert(integerObject || integerObject2); |
---|
[201] | 328 | #endif |
---|
[1286] | 329 | if (upper[iColumn] - lower[iColumn] > 1.000001) { |
---|
| 330 | general++; |
---|
| 331 | if (doGeneral) |
---|
| 332 | integerVariable[j++] = iColumn; |
---|
| 333 | } else { |
---|
| 334 | integerVariable[j++] = iColumn; |
---|
| 335 | } |
---|
[175] | 336 | } |
---|
[1364] | 337 | /* |
---|
| 338 | If 2/3 of integers are general integers, and we're not going to work with |
---|
| 339 | them, might as well go home. |
---|
| 340 | |
---|
| 341 | The else case is unclear to me. We reach it if general integers are less than |
---|
| 342 | 2/3 of the total, or if either of bit 4 or 8 is set. But only bit 8 is used |
---|
| 343 | in the decision. (Let manyGen = 1 if more than 2/3 of integers are general |
---|
| 344 | integers. Then a k-map on manyGen, bit4, and bit8 shows it clearly.) |
---|
| 345 | |
---|
| 346 | So there's something odd here. In the case where bit4 = 1 and bit8 = 0, |
---|
| 347 | we've included general integers in integerVariable, but we're not going to |
---|
| 348 | process them. |
---|
| 349 | */ |
---|
[1286] | 350 | if (general*3 > 2*numberIntegers && !doGeneral) { |
---|
| 351 | delete [] integerVariable; |
---|
| 352 | return 0; |
---|
| 353 | } else if ((accumulate_&4) == 0) { |
---|
| 354 | doGeneral = false; |
---|
| 355 | j = 0; |
---|
| 356 | for (i = 0; i < numberIntegers; i++) { |
---|
| 357 | int iColumn = integerVariableOrig[i]; |
---|
| 358 | if (upper[iColumn] - lower[iColumn] < 1.000001) |
---|
| 359 | integerVariable[j++] = iColumn; |
---|
| 360 | } |
---|
[1053] | 361 | } |
---|
[1286] | 362 | if (!general) |
---|
| 363 | doGeneral = false; |
---|
[1053] | 364 | #ifdef CLP_INVESTIGATE |
---|
[1286] | 365 | if (doGeneral) |
---|
| 366 | printf("DOing general with %d out of %d\n", general, numberIntegers); |
---|
[1053] | 367 | #endif |
---|
[1364] | 368 | /* |
---|
| 369 | This `closest solution' will satisfy integrality, but violate some other |
---|
| 370 | constraints? |
---|
| 371 | */ |
---|
[1286] | 372 | // For solution closest to feasible if none found |
---|
| 373 | int * closestSolution = general ? NULL : new int[numberIntegers]; |
---|
| 374 | double closestObjectiveValue = COIN_DBL_MAX; |
---|
[1271] | 375 | |
---|
[1286] | 376 | int numberIntegersOrig = numberIntegers; |
---|
| 377 | numberIntegers = j; |
---|
| 378 | double * newSolution = new double [numberColumns]; |
---|
| 379 | double newSolutionValue = COIN_DBL_MAX; |
---|
| 380 | bool solutionFound = false; |
---|
| 381 | int * usedColumn = NULL; |
---|
| 382 | double * lastSolution = NULL; |
---|
| 383 | int fixContinuous = 0; |
---|
| 384 | bool fixInternal = false; |
---|
| 385 | bool allSlack = false; |
---|
| 386 | if (when_ >= 21 && when_ <= 25) { |
---|
| 387 | when_ -= 10; |
---|
| 388 | allSlack = true; |
---|
[753] | 389 | } |
---|
[1286] | 390 | double time1 = CoinCpuTime(); |
---|
[1364] | 391 | /* |
---|
| 392 | Obtain a relaxed lp solution. |
---|
| 393 | */ |
---|
[1286] | 394 | model_->solver()->resolve(); |
---|
| 395 | if (!model_->solver()->isProvenOptimal()) { |
---|
| 396 | // presumably max time or some such |
---|
| 397 | return 0; |
---|
[1115] | 398 | } |
---|
[1286] | 399 | numRuns_++; |
---|
| 400 | if (cutoff < 1.0e50 && false) { |
---|
| 401 | // Fix on djs |
---|
| 402 | double direction = model_->solver()->getObjSense() ; |
---|
| 403 | double gap = cutoff - model_->solver()->getObjValue() * direction ; |
---|
| 404 | double tolerance; |
---|
| 405 | model_->solver()->getDblParam(OsiDualTolerance, tolerance) ; |
---|
| 406 | if (gap > 0.0) { |
---|
| 407 | gap += 100.0 * tolerance; |
---|
| 408 | int nFix = model_->solver()->reducedCostFix(gap); |
---|
| 409 | printf("dj fixing fixed %d variables\n", nFix); |
---|
| 410 | } |
---|
| 411 | } |
---|
[1364] | 412 | /* |
---|
| 413 | I have no idea why we're doing this, except perhaps that saveBasis will be |
---|
| 414 | automagically deleted on exit from the routine. |
---|
| 415 | */ |
---|
[1286] | 416 | CoinWarmStartBasis saveBasis; |
---|
| 417 | CoinWarmStartBasis * basis = |
---|
| 418 | dynamic_cast<CoinWarmStartBasis *>(model_->solver()->getWarmStart()) ; |
---|
| 419 | if (basis) { |
---|
| 420 | saveBasis = * basis; |
---|
| 421 | delete basis; |
---|
| 422 | } |
---|
| 423 | double continuousObjectiveValue = model_->solver()->getObjValue() * model_->solver()->getObjSense(); |
---|
| 424 | double * firstPerturbedObjective = NULL; |
---|
| 425 | double * firstPerturbedSolution = NULL; |
---|
| 426 | double firstPerturbedValue = COIN_DBL_MAX; |
---|
| 427 | if (when_ >= 11 && when_ <= 15) { |
---|
| 428 | fixInternal = when_ > 11 && when_ < 15; |
---|
| 429 | if (when_ < 13) |
---|
| 430 | fixContinuous = 0; |
---|
| 431 | else if (when_ != 14) |
---|
| 432 | fixContinuous = 1; |
---|
| 433 | else |
---|
| 434 | fixContinuous = 2; |
---|
| 435 | when_ = 1; |
---|
| 436 | if ((accumulate_&1) != 0) { |
---|
| 437 | usedColumn = new int [numberColumns]; |
---|
| 438 | for (int i = 0; i < numberColumns; i++) |
---|
| 439 | usedColumn[i] = -1; |
---|
| 440 | } |
---|
| 441 | lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(), numberColumns); |
---|
| 442 | } |
---|
| 443 | int finalReturnCode = 0; |
---|
| 444 | int totalNumberPasses = 0; |
---|
| 445 | int numberTries = 0; |
---|
| 446 | CoinWarmStartBasis bestBasis; |
---|
| 447 | bool exitAll = false; |
---|
| 448 | //double saveBestObjective = model_->getMinimizationObjValue(); |
---|
| 449 | int numberSolutions = 0; |
---|
| 450 | OsiSolverInterface * solver = NULL; |
---|
| 451 | double artificialFactor = 0.00001; |
---|
| 452 | // also try rounding! |
---|
| 453 | double * roundingSolution = new double[numberColumns]; |
---|
| 454 | double roundingObjective = solutionValue; |
---|
| 455 | CbcRounding roundingHeuristic(*model_); |
---|
| 456 | int dualPass = 0; |
---|
| 457 | int secondPassOpt = 0; |
---|
[1271] | 458 | #define RAND_RAND |
---|
| 459 | #ifdef RAND_RAND |
---|
[1286] | 460 | int offRandom = 0; |
---|
[1271] | 461 | #endif |
---|
[1286] | 462 | int maximumAllowed = -1; |
---|
| 463 | bool moreIterations = false; |
---|
| 464 | if (options > 0) { |
---|
| 465 | if (options >= 1000) |
---|
| 466 | maximumAllowed = options / 1000; |
---|
| 467 | int options2 = (options % 1000) / 100; |
---|
[1271] | 468 | #ifdef RAND_RAND |
---|
[1286] | 469 | offRandom = options2 & 1; |
---|
[1271] | 470 | #endif |
---|
[1286] | 471 | moreIterations = (options2 & 2) != 0; |
---|
| 472 | secondPassOpt = (options / 10) % 10; |
---|
| 473 | /* 1 to 7 - re-use solution |
---|
| 474 | 8 use dual and current solution(ish) |
---|
| 475 | 9 use dual and allslack |
---|
| 476 | 1 - primal and mod obj |
---|
| 477 | 2 - dual and mod obj |
---|
| 478 | 3 - primal and no mod obj |
---|
| 479 | add 3 to redo current solution |
---|
| 480 | */ |
---|
| 481 | if (secondPassOpt >= 8) { |
---|
| 482 | dualPass = secondPassOpt - 7; |
---|
| 483 | secondPassOpt = 0; |
---|
| 484 | } |
---|
[1053] | 485 | } |
---|
[1286] | 486 | // Number of passes to do |
---|
| 487 | int maximumPasses = maximumPasses_; |
---|
[1271] | 488 | #ifdef COIN_HAS_CLP |
---|
[1286] | 489 | if (maximumPasses == 30) { |
---|
| 490 | OsiClpSolverInterface * clpSolver |
---|
| 491 | = dynamic_cast<OsiClpSolverInterface *> (model_->solver()); |
---|
| 492 | if (clpSolver && clpSolver->fakeObjective()) |
---|
| 493 | maximumPasses = 100; // feasibility problem? |
---|
| 494 | } |
---|
[1271] | 495 | #endif |
---|
| 496 | #ifdef RAND_RAND |
---|
[1286] | 497 | double * randomFactor = new double [numberColumns]; |
---|
| 498 | for (int i = 0; i < numberColumns; i++) { |
---|
| 499 | double value = floor(1.0e3 * randomNumberGenerator_.randomDouble()); |
---|
| 500 | randomFactor[i] = 1.0 + value * 1.0e-4; |
---|
| 501 | } |
---|
[1271] | 502 | #endif |
---|
[1286] | 503 | // guess exact multiple of objective |
---|
| 504 | double exactMultiple = model_->getCutoffIncrement(); |
---|
| 505 | exactMultiple *= 2520; |
---|
| 506 | if (fabs(exactMultiple / 0.999 - floor(exactMultiple / 0.999 + 0.5)) < 1.0e-9) |
---|
| 507 | exactMultiple /= 2520.0 * 0.999; |
---|
| 508 | else if (fabs(exactMultiple - floor(exactMultiple + 0.5)) < 1.0e-9) |
---|
| 509 | exactMultiple /= 2520.0; |
---|
| 510 | else |
---|
| 511 | exactMultiple = 0.0; |
---|
| 512 | // check for rounding errors (only for integral case) |
---|
| 513 | if (fabs(exactMultiple - floor(exactMultiple + 0.5)) < 1.0e-8) |
---|
| 514 | exactMultiple = floor(exactMultiple + 0.5); |
---|
| 515 | //printf("exact multiple %g\n",exactMultiple); |
---|
| 516 | // Clone solver for rounding |
---|
| 517 | OsiSolverInterface * clonedSolver = cloneBut(2); // wasmodel_->solver()->clone(); |
---|
| 518 | while (!exitAll) { |
---|
| 519 | // Cutoff rhs |
---|
| 520 | double useRhs = COIN_DBL_MAX; |
---|
| 521 | double useOffset = 0.0; |
---|
| 522 | int numberPasses = 0; |
---|
| 523 | artificialFactor *= 10.0; |
---|
| 524 | int lastMove = (!numberTries) ? -10 : 1000000; |
---|
| 525 | double lastSumInfeas = COIN_DBL_MAX; |
---|
| 526 | numberTries++; |
---|
| 527 | // Clone solver - otherwise annoys root node computations |
---|
| 528 | solver = cloneBut(2); // was model_->solver()->clone(); |
---|
[1053] | 529 | #ifdef COIN_HAS_CLP |
---|
[1286] | 530 | { |
---|
| 531 | OsiClpSolverInterface * clpSolver |
---|
| 532 | = dynamic_cast<OsiClpSolverInterface *> (solver); |
---|
| 533 | if (clpSolver) { |
---|
| 534 | // better to clean up using primal? |
---|
| 535 | ClpSimplex * lp = clpSolver->getModelPtr(); |
---|
| 536 | int options = lp->specialOptions(); |
---|
| 537 | lp->setSpecialOptions(options | 8192); |
---|
| 538 | //lp->setSpecialOptions(options|0x01000000); |
---|
[1271] | 539 | #ifdef CLP_INVESTIGATE |
---|
[1286] | 540 | clpSolver->setHintParam(OsiDoReducePrint, false, OsiHintTry); |
---|
| 541 | lp->setLogLevel(CoinMax(1, lp->logLevel())); |
---|
[1271] | 542 | #endif |
---|
[1286] | 543 | } |
---|
| 544 | } |
---|
[1053] | 545 | #endif |
---|
[1286] | 546 | if (CoinMin(fakeCutoff_, cutoff) < 1.0e50) { |
---|
| 547 | // Fix on djs |
---|
| 548 | double direction = solver->getObjSense() ; |
---|
| 549 | double gap = CoinMin(fakeCutoff_, cutoff) - solver->getObjValue() * direction ; |
---|
| 550 | double tolerance; |
---|
| 551 | solver->getDblParam(OsiDualTolerance, tolerance) ; |
---|
| 552 | if (gap > 0.0 && (fixOnReducedCosts_ == 1 || (numberTries == 1 && fixOnReducedCosts_ == 2))) { |
---|
| 553 | gap += 100.0 * tolerance; |
---|
| 554 | gap *= reducedCostMultiplier_; |
---|
| 555 | int nFix = solver->reducedCostFix(gap); |
---|
| 556 | if (nFix) { |
---|
| 557 | sprintf(pumpPrint, "Reduced cost fixing fixed %d variables on major pass %d", nFix, numberTries); |
---|
| 558 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 559 | << pumpPrint |
---|
| 560 | << CoinMessageEol; |
---|
| 561 | //pumpPrint[0]='\0'; |
---|
| 562 | } |
---|
| 563 | } |
---|
| 564 | } |
---|
| 565 | // if cutoff exists then add constraint |
---|
| 566 | bool useCutoff = (fabs(cutoff) < 1.0e20 && (fakeCutoff_ != COIN_DBL_MAX || numberTries > 1)); |
---|
| 567 | // but there may be a close one |
---|
| 568 | if (firstCutoff < 2.0*solutionValue && numberTries == 1 && CoinMin(cutoff, fakeCutoff_) < 1.0e20) |
---|
| 569 | useCutoff = true; |
---|
| 570 | if (useCutoff) { |
---|
| 571 | double rhs = CoinMin(cutoff, fakeCutoff_); |
---|
| 572 | const double * objective = solver->getObjCoefficients(); |
---|
| 573 | int numberColumns = solver->getNumCols(); |
---|
| 574 | int * which = new int[numberColumns]; |
---|
| 575 | double * els = new double[numberColumns]; |
---|
| 576 | int nel = 0; |
---|
| 577 | for (int i = 0; i < numberColumns; i++) { |
---|
| 578 | double value = objective[i]; |
---|
| 579 | if (value) { |
---|
| 580 | which[nel] = i; |
---|
| 581 | els[nel++] = direction * value; |
---|
| 582 | } |
---|
| 583 | } |
---|
| 584 | solver->getDblParam(OsiObjOffset, useOffset); |
---|
[753] | 585 | #ifdef COIN_DEVELOP |
---|
[1286] | 586 | if (useOffset) |
---|
| 587 | printf("CbcHeuristicFPump obj offset %g\n", useOffset); |
---|
[753] | 588 | #endif |
---|
[1286] | 589 | useOffset *= direction; |
---|
| 590 | // Tweak rhs and save |
---|
| 591 | useRhs = rhs; |
---|
[1393] | 592 | #ifdef JJF_ZERO |
---|
[1286] | 593 | double tempValue = 60.0 * useRhs; |
---|
| 594 | if (fabs(tempValue - floor(tempValue + 0.5)) < 1.0e-7 && rhs != fakeCutoff_) { |
---|
| 595 | // add a little |
---|
| 596 | useRhs += 1.0e-5; |
---|
| 597 | } |
---|
[1117] | 598 | #endif |
---|
[1286] | 599 | solver->addRow(nel, which, els, -COIN_DBL_MAX, useRhs + useOffset); |
---|
| 600 | delete [] which; |
---|
| 601 | delete [] els; |
---|
| 602 | bool takeHint; |
---|
| 603 | OsiHintStrength strength; |
---|
| 604 | solver->getHintParam(OsiDoDualInResolve, takeHint, strength); |
---|
| 605 | solver->setHintParam(OsiDoDualInResolve, true, OsiHintDo); |
---|
| 606 | solver->resolve(); |
---|
| 607 | solver->setHintParam(OsiDoDualInResolve, takeHint, strength); |
---|
| 608 | if (!solver->isProvenOptimal()) { |
---|
| 609 | // presumably max time or some such |
---|
| 610 | exitAll = true; |
---|
| 611 | break; |
---|
| 612 | } |
---|
| 613 | } |
---|
| 614 | solver->setDblParam(OsiDualObjectiveLimit, 1.0e50); |
---|
| 615 | solver->resolve(); |
---|
| 616 | // Solver may not be feasible |
---|
| 617 | if (!solver->isProvenOptimal()) { |
---|
| 618 | exitAll = true; |
---|
| 619 | break; |
---|
| 620 | } |
---|
| 621 | const double * lower = solver->getColLower(); |
---|
| 622 | const double * upper = solver->getColUpper(); |
---|
| 623 | const double * solution = solver->getColSolution(); |
---|
| 624 | if (lastSolution) |
---|
| 625 | memcpy(lastSolution, solution, numberColumns*sizeof(double)); |
---|
| 626 | double primalTolerance; |
---|
| 627 | solver->getDblParam(OsiPrimalTolerance, primalTolerance); |
---|
| 628 | |
---|
| 629 | // 2 space for last rounded solutions |
---|
[640] | 630 | #define NUMBER_OLD 4 |
---|
[1286] | 631 | double ** oldSolution = new double * [NUMBER_OLD]; |
---|
| 632 | for (j = 0; j < NUMBER_OLD; j++) { |
---|
| 633 | oldSolution[j] = new double[numberColumns]; |
---|
| 634 | for (i = 0; i < numberColumns; i++) oldSolution[j][i] = -COIN_DBL_MAX; |
---|
| 635 | } |
---|
| 636 | |
---|
| 637 | // 3. Replace objective with an initial 0-valued objective |
---|
| 638 | double * saveObjective = new double [numberColumns]; |
---|
| 639 | memcpy(saveObjective, solver->getObjCoefficients(), numberColumns*sizeof(double)); |
---|
| 640 | for (i = 0; i < numberColumns; i++) { |
---|
| 641 | solver->setObjCoeff(i, 0.0); |
---|
| 642 | } |
---|
| 643 | bool finished = false; |
---|
| 644 | double direction = solver->getObjSense(); |
---|
| 645 | int returnCode = 0; |
---|
| 646 | bool takeHint; |
---|
| 647 | OsiHintStrength strength; |
---|
| 648 | solver->getHintParam(OsiDoDualInResolve, takeHint, strength); |
---|
| 649 | solver->setHintParam(OsiDoDualInResolve, false); |
---|
| 650 | //solver->messageHandler()->setLogLevel(0); |
---|
| 651 | |
---|
| 652 | // 4. Save objective offset so we can see progress |
---|
| 653 | double saveOffset; |
---|
| 654 | solver->getDblParam(OsiObjOffset, saveOffset); |
---|
| 655 | // Get amount for original objective |
---|
| 656 | double scaleFactor = 0.0; |
---|
[833] | 657 | #ifdef COIN_DEVELOP |
---|
[1286] | 658 | double largestCost = 0.0; |
---|
| 659 | int nArtificial = 0; |
---|
[833] | 660 | #endif |
---|
[1286] | 661 | for (i = 0; i < numberColumns; i++) { |
---|
| 662 | double value = saveObjective[i]; |
---|
| 663 | scaleFactor += value * value; |
---|
[833] | 664 | #ifdef COIN_DEVELOP |
---|
[1286] | 665 | largestCost = CoinMax(largestCost, fabs(value)); |
---|
| 666 | if (value*direction >= artificialCost_) |
---|
| 667 | nArtificial++; |
---|
[833] | 668 | #endif |
---|
[1286] | 669 | } |
---|
| 670 | if (scaleFactor) |
---|
| 671 | scaleFactor = (initialWeight_ * sqrt(static_cast<double> (numberIntegers))) / sqrt(scaleFactor); |
---|
[1271] | 672 | #ifdef CLP_INVESTIGATE |
---|
[833] | 673 | #ifdef COIN_DEVELOP |
---|
[1286] | 674 | if (scaleFactor || nArtificial) |
---|
| 675 | printf("Using %g fraction of original objective (decay %g) - largest %g - %d artificials\n", scaleFactor, weightFactor_, |
---|
| 676 | largestCost, nArtificial); |
---|
[1271] | 677 | #else |
---|
[1286] | 678 | if (scaleFactor) |
---|
| 679 | printf("Using %g fraction of original objective (decay %g)\n", |
---|
| 680 | scaleFactor, weightFactor_); |
---|
[833] | 681 | #endif |
---|
[1271] | 682 | #endif |
---|
[1286] | 683 | // This is an array of sums of infeasibilities so can see if "bobbling" |
---|
[1127] | 684 | #define SIZE_BOBBLE 20 |
---|
[1286] | 685 | double saveSumInf[SIZE_BOBBLE]; |
---|
| 686 | CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX); |
---|
| 687 | // 0 before doing anything |
---|
| 688 | int bobbleMode = 0; |
---|
| 689 | // 5. MAIN WHILE LOOP |
---|
| 690 | //bool newLineNeeded=false; |
---|
[1364] | 691 | /* |
---|
| 692 | finished occurs exactly twice in this routine: immediately above, where it's |
---|
| 693 | set to false, and here in the loop condition. |
---|
| 694 | */ |
---|
[1286] | 695 | while (!finished) { |
---|
| 696 | double newTrueSolutionValue = 0.0; |
---|
| 697 | double newSumInfeas = 0.0; |
---|
| 698 | int newNumberInfeas = 0; |
---|
| 699 | returnCode = 0; |
---|
| 700 | if (model_->maximumSecondsReached()) { |
---|
| 701 | exitAll = true; |
---|
| 702 | break; |
---|
| 703 | } |
---|
| 704 | // see what changed |
---|
| 705 | if (usedColumn) { |
---|
| 706 | for (i = 0; i < numberColumns; i++) { |
---|
| 707 | if (fabs(solution[i] - lastSolution[i]) > 1.0e-8) |
---|
| 708 | usedColumn[i] = numberPasses; |
---|
| 709 | lastSolution[i] = solution[i]; |
---|
| 710 | } |
---|
| 711 | } |
---|
| 712 | if (averageIterationsPerTry >= 0) { |
---|
| 713 | int n = totalNumberIterations - numberIterationsLastPass; |
---|
| 714 | double perPass = totalNumberIterations / |
---|
| 715 | (totalNumberPasses + numberPasses + 1.0e-5); |
---|
| 716 | perPass /= (solver->getNumRows() + numberColumns); |
---|
| 717 | double test = moreIterations ? 0.3 : 0.05; |
---|
| 718 | if (n > CoinMax(20000, 3*averageIterationsPerTry) |
---|
| 719 | && (switches_&2) == 0 && maximumPasses<200 && perPass>test) { |
---|
| 720 | exitAll = true; |
---|
| 721 | } |
---|
| 722 | } |
---|
| 723 | // Exit on exact total number if maximumPasses large |
---|
| 724 | if ((maximumPasses >= 200 || (switches_&2) != 0) |
---|
| 725 | && numberPasses + totalNumberPasses >= |
---|
| 726 | maximumPasses) |
---|
| 727 | exitAll = true; |
---|
| 728 | bool exitThis = false; |
---|
| 729 | if (iterationLimit < 0.0) { |
---|
| 730 | if (numberPasses >= maximumPasses) { |
---|
| 731 | // If going well then keep going if maximumPasses small |
---|
| 732 | if (lastMove < numberPasses - 4 || lastMove == 1000000) |
---|
| 733 | exitThis = true; |
---|
| 734 | if (maximumPasses > 20 || numberPasses >= 40) |
---|
| 735 | exitThis = true; |
---|
| 736 | } |
---|
| 737 | } |
---|
| 738 | if (iterationLimit > 0.0 && totalNumberIterations > iterationLimit |
---|
| 739 | && numberPasses > 15) { |
---|
| 740 | // exiting on iteration count |
---|
| 741 | exitAll = true; |
---|
| 742 | } else if (maximumPasses<30 && numberPasses>100) { |
---|
| 743 | // too many passes anyway |
---|
| 744 | exitAll = true; |
---|
| 745 | } |
---|
| 746 | if (maximumTime_ > 0.0 && CoinCpuTime() >= startTime_ + maximumTime_) { |
---|
| 747 | exitAll = true; |
---|
| 748 | // force exit |
---|
| 749 | switches_ |= 2048; |
---|
| 750 | } |
---|
| 751 | if (exitAll || exitThis) |
---|
| 752 | break; |
---|
| 753 | memcpy(newSolution, solution, numberColumns*sizeof(double)); |
---|
| 754 | int flip; |
---|
| 755 | if (numberPasses == 0 && false) { |
---|
| 756 | // always use same seed |
---|
| 757 | randomNumberGenerator_.setSeed(987654321); |
---|
| 758 | } |
---|
| 759 | returnCode = rounds(solver, newSolution,/*saveObjective,*/ |
---|
| 760 | numberIntegers, integerVariable, |
---|
| 761 | /*pumpPrint,*/numberPasses, |
---|
| 762 | /*roundExpensive_,*/defaultRounding_, &flip); |
---|
| 763 | if (numberPasses == 0 && false) { |
---|
| 764 | // Make sure random will be different |
---|
| 765 | for (i = 1; i < numberTries; i++) |
---|
| 766 | randomNumberGenerator_.randomDouble(); |
---|
| 767 | } |
---|
| 768 | numberPasses++; |
---|
| 769 | if (returnCode) { |
---|
| 770 | // SOLUTION IS INTEGER |
---|
| 771 | // Put back correct objective |
---|
| 772 | for (i = 0; i < numberColumns; i++) |
---|
| 773 | solver->setObjCoeff(i, saveObjective[i]); |
---|
[640] | 774 | |
---|
[1286] | 775 | // solution - but may not be better |
---|
| 776 | // Compute using dot product |
---|
| 777 | solver->setDblParam(OsiObjOffset, saveOffset); |
---|
| 778 | newSolutionValue = -saveOffset; |
---|
| 779 | for ( i = 0 ; i < numberColumns ; i++ ) |
---|
| 780 | newSolutionValue += saveObjective[i] * newSolution[i]; |
---|
| 781 | newSolutionValue *= direction; |
---|
| 782 | sprintf(pumpPrint, "Solution found of %g", newSolutionValue); |
---|
| 783 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 784 | << pumpPrint |
---|
| 785 | << CoinMessageEol; |
---|
| 786 | //newLineNeeded=false; |
---|
| 787 | if (newSolutionValue < solutionValue) { |
---|
| 788 | double saveValue = solutionValue; |
---|
| 789 | if (!doGeneral) { |
---|
| 790 | int numberLeft = 0; |
---|
| 791 | for (i = 0; i < numberIntegersOrig; i++) { |
---|
| 792 | int iColumn = integerVariableOrig[i]; |
---|
| 793 | double value = floor(newSolution[iColumn] + 0.5); |
---|
| 794 | if (solver->isBinary(iColumn)) { |
---|
| 795 | solver->setColLower(iColumn, value); |
---|
| 796 | solver->setColUpper(iColumn, value); |
---|
| 797 | } else { |
---|
| 798 | if (fabs(value - newSolution[iColumn]) > 1.0e-7) |
---|
| 799 | numberLeft++; |
---|
| 800 | } |
---|
| 801 | } |
---|
| 802 | if (numberLeft) { |
---|
| 803 | sprintf(pumpPrint, "Branch and bound needed to clear up %d general integers", numberLeft); |
---|
| 804 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 805 | << pumpPrint |
---|
| 806 | << CoinMessageEol; |
---|
| 807 | returnCode = smallBranchAndBound(solver, numberNodes_, newSolution, newSolutionValue, |
---|
| 808 | solutionValue, "CbcHeuristicFpump"); |
---|
| 809 | if (returnCode < 0) { |
---|
| 810 | if (returnCode == -2) |
---|
| 811 | exitAll = true; |
---|
| 812 | returnCode = 0; // returned on size or event |
---|
| 813 | } |
---|
| 814 | if ((returnCode&2) != 0) { |
---|
| 815 | // could add cut |
---|
| 816 | returnCode &= ~2; |
---|
| 817 | } |
---|
| 818 | if (returnCode != 1) |
---|
| 819 | newSolutionValue = saveValue; |
---|
| 820 | if (returnCode && newSolutionValue < saveValue) |
---|
| 821 | numberBandBsolutions++; |
---|
| 822 | } |
---|
| 823 | } |
---|
| 824 | if (returnCode && newSolutionValue < saveValue) { |
---|
| 825 | memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); |
---|
| 826 | solutionFound = true; |
---|
| 827 | if (exitNow(newSolutionValue)) |
---|
| 828 | exitAll = true; |
---|
| 829 | CoinWarmStartBasis * basis = |
---|
| 830 | dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; |
---|
| 831 | if (basis) { |
---|
| 832 | bestBasis = * basis; |
---|
| 833 | delete basis; |
---|
| 834 | int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution); |
---|
| 835 | if (action == 0) { |
---|
| 836 | double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns); |
---|
| 837 | double saveObjectiveValue = model_->getMinimizationObjValue(); |
---|
| 838 | model_->setBestSolution(betterSolution, numberColumns, newSolutionValue); |
---|
| 839 | if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue()) |
---|
| 840 | model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue); |
---|
| 841 | delete [] saveOldSolution; |
---|
| 842 | } |
---|
| 843 | if (action == 0 || model_->maximumSecondsReached()) { |
---|
| 844 | exitAll = true; // exit |
---|
| 845 | break; |
---|
| 846 | } |
---|
| 847 | } |
---|
| 848 | if ((accumulate_&1) != 0) { |
---|
| 849 | model_->incrementUsed(betterSolution); // for local search |
---|
| 850 | numberSolutions++; |
---|
| 851 | } |
---|
| 852 | solutionValue = newSolutionValue; |
---|
| 853 | solutionFound = true; |
---|
| 854 | if (general && saveValue != newSolutionValue) { |
---|
| 855 | sprintf(pumpPrint, "Cleaned solution of %g", solutionValue); |
---|
| 856 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 857 | << pumpPrint |
---|
| 858 | << CoinMessageEol; |
---|
| 859 | } |
---|
| 860 | if (exitNow(newSolutionValue)) |
---|
| 861 | exitAll = true; |
---|
| 862 | } else { |
---|
| 863 | sprintf(pumpPrint, "Mini branch and bound could not fix general integers"); |
---|
| 864 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 865 | << pumpPrint |
---|
| 866 | << CoinMessageEol; |
---|
| 867 | } |
---|
| 868 | } else { |
---|
| 869 | sprintf(pumpPrint, "After further testing solution no better than previous of %g", solutionValue); |
---|
| 870 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 871 | << pumpPrint |
---|
| 872 | << CoinMessageEol; |
---|
| 873 | //newLineNeeded=false; |
---|
| 874 | returnCode = 0; |
---|
| 875 | } |
---|
| 876 | break; |
---|
| 877 | } else { |
---|
| 878 | // SOLUTION IS not INTEGER |
---|
| 879 | // 1. check for loop |
---|
| 880 | bool matched; |
---|
| 881 | for (int k = NUMBER_OLD - 1; k > 0; k--) { |
---|
| 882 | double * b = oldSolution[k]; |
---|
| 883 | matched = true; |
---|
| 884 | for (i = 0; i < numberIntegers; i++) { |
---|
| 885 | int iColumn = integerVariable[i]; |
---|
| 886 | if (newSolution[iColumn] != b[iColumn]) { |
---|
| 887 | matched = false; |
---|
| 888 | break; |
---|
| 889 | } |
---|
| 890 | } |
---|
| 891 | if (matched) break; |
---|
| 892 | } |
---|
| 893 | int numberPerturbed = 0; |
---|
| 894 | if (matched || numberPasses % 100 == 0) { |
---|
| 895 | // perturbation |
---|
| 896 | //sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied"); |
---|
| 897 | //newLineNeeded=true; |
---|
| 898 | double factorX[10] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0}; |
---|
| 899 | double factor = 1.0; |
---|
| 900 | double target = -1.0; |
---|
| 901 | double * randomX = new double [numberIntegers]; |
---|
| 902 | for (i = 0; i < numberIntegers; i++) |
---|
| 903 | randomX[i] = CoinMax(0.0, randomNumberGenerator_.randomDouble() - 0.3); |
---|
| 904 | for (int k = 0; k < 10; k++) { |
---|
[753] | 905 | #ifdef COIN_DEVELOP_x |
---|
[1286] | 906 | printf("kpass %d\n", k); |
---|
[750] | 907 | #endif |
---|
[1286] | 908 | int numberX[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; |
---|
| 909 | for (i = 0; i < numberIntegers; i++) { |
---|
| 910 | int iColumn = integerVariable[i]; |
---|
| 911 | double value = randomX[i]; |
---|
| 912 | double difference = fabs(solution[iColumn] - newSolution[iColumn]); |
---|
| 913 | for (int j = 0; j < 10; j++) { |
---|
| 914 | if (difference + value*factorX[j] > 0.5) |
---|
| 915 | numberX[j]++; |
---|
| 916 | } |
---|
| 917 | } |
---|
| 918 | if (target < 0.0) { |
---|
| 919 | if (numberX[9] <= 200) |
---|
| 920 | break; // not very many changes |
---|
| 921 | target = CoinMax(200.0, CoinMin(0.05 * numberX[9], 1000.0)); |
---|
| 922 | } |
---|
| 923 | int iX = -1; |
---|
| 924 | int iBand = -1; |
---|
| 925 | for (i = 0; i < 10; i++) { |
---|
[753] | 926 | #ifdef COIN_DEVELOP_x |
---|
[1286] | 927 | printf("** %d changed at %g\n", numberX[i], factorX[i]); |
---|
[750] | 928 | #endif |
---|
[1286] | 929 | if (numberX[i] >= target && numberX[i] < 2.0*target && iX < 0) |
---|
| 930 | iX = i; |
---|
| 931 | if (iBand<0 && numberX[i]>target) { |
---|
| 932 | iBand = i; |
---|
| 933 | factor = factorX[i]; |
---|
| 934 | } |
---|
| 935 | } |
---|
| 936 | if (iX >= 0) { |
---|
| 937 | factor = factorX[iX]; |
---|
| 938 | break; |
---|
| 939 | } else { |
---|
| 940 | assert (iBand >= 0); |
---|
| 941 | double hi = factor; |
---|
| 942 | double lo = (iBand > 0) ? factorX[iBand-1] : 0.0; |
---|
| 943 | double diff = (hi - lo) / 9.0; |
---|
| 944 | for (i = 0; i < 10; i++) { |
---|
| 945 | factorX[i] = lo; |
---|
| 946 | lo += diff; |
---|
| 947 | } |
---|
| 948 | } |
---|
| 949 | } |
---|
| 950 | for (i = 0; i < numberIntegers; i++) { |
---|
| 951 | int iColumn = integerVariable[i]; |
---|
| 952 | double value = randomX[i]; |
---|
| 953 | double difference = fabs(solution[iColumn] - newSolution[iColumn]); |
---|
| 954 | if (difference + value*factor > 0.5) { |
---|
| 955 | numberPerturbed++; |
---|
| 956 | if (newSolution[iColumn] < lower[iColumn] + primalTolerance) { |
---|
| 957 | newSolution[iColumn] += 1.0; |
---|
| 958 | } else if (newSolution[iColumn] > upper[iColumn] - primalTolerance) { |
---|
| 959 | newSolution[iColumn] -= 1.0; |
---|
| 960 | } else { |
---|
| 961 | // general integer |
---|
| 962 | if (difference + value > 0.75) |
---|
| 963 | newSolution[iColumn] += 1.0; |
---|
| 964 | else |
---|
| 965 | newSolution[iColumn] -= 1.0; |
---|
| 966 | } |
---|
| 967 | } |
---|
| 968 | } |
---|
| 969 | delete [] randomX; |
---|
| 970 | } else { |
---|
| 971 | for (j = NUMBER_OLD - 1; j > 0; j--) { |
---|
| 972 | for (i = 0; i < numberColumns; i++) oldSolution[j][i] = oldSolution[j-1][i]; |
---|
| 973 | } |
---|
| 974 | for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j]; |
---|
| 975 | } |
---|
| 976 | |
---|
| 977 | // 2. update the objective function based on the new rounded solution |
---|
| 978 | double offset = 0.0; |
---|
| 979 | double costValue = (1.0 - scaleFactor) * solver->getObjSense(); |
---|
| 980 | int numberChanged = 0; |
---|
| 981 | const double * oldObjective = solver->getObjCoefficients(); |
---|
| 982 | for (i = 0; i < numberColumns; i++) { |
---|
| 983 | // below so we can keep original code and allow for objective |
---|
| 984 | int iColumn = i; |
---|
| 985 | // Special code for "artificials" |
---|
| 986 | if (direction*saveObjective[iColumn] >= artificialCost_) { |
---|
| 987 | //solver->setObjCoeff(iColumn,scaleFactor*saveObjective[iColumn]); |
---|
| 988 | solver->setObjCoeff(iColumn, (artificialFactor*saveObjective[iColumn]) / artificialCost_); |
---|
| 989 | } |
---|
| 990 | if (!solver->isBinary(iColumn) && !doGeneral) |
---|
| 991 | continue; |
---|
| 992 | // deal with fixed variables (i.e., upper=lower) |
---|
| 993 | if (fabs(lower[iColumn] - upper[iColumn]) < primalTolerance || !solver->isInteger(iColumn)) { |
---|
| 994 | //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue); |
---|
| 995 | //else solver->setObjCoeff(iColumn,costValue); |
---|
| 996 | continue; |
---|
| 997 | } |
---|
| 998 | double newValue = 0.0; |
---|
| 999 | if (newSolution[iColumn] < lower[iColumn] + primalTolerance) { |
---|
| 1000 | newValue = costValue + scaleFactor * saveObjective[iColumn]; |
---|
| 1001 | } else { |
---|
| 1002 | if (newSolution[iColumn] > upper[iColumn] - primalTolerance) { |
---|
| 1003 | newValue = -costValue + scaleFactor * saveObjective[iColumn]; |
---|
| 1004 | } |
---|
| 1005 | } |
---|
[1271] | 1006 | #ifdef RAND_RAND |
---|
[1286] | 1007 | if (!offRandom) |
---|
| 1008 | newValue *= randomFactor[iColumn]; |
---|
[1271] | 1009 | #endif |
---|
[1286] | 1010 | if (newValue != oldObjective[iColumn]) { |
---|
| 1011 | numberChanged++; |
---|
| 1012 | } |
---|
| 1013 | solver->setObjCoeff(iColumn, newValue); |
---|
| 1014 | offset += costValue * newSolution[iColumn]; |
---|
| 1015 | } |
---|
| 1016 | solver->setDblParam(OsiObjOffset, -offset); |
---|
| 1017 | if (!general && false) { |
---|
| 1018 | // Solve in two goes - first keep satisfied ones fixed |
---|
| 1019 | double * saveLower = new double [numberIntegers]; |
---|
| 1020 | double * saveUpper = new double [numberIntegers]; |
---|
| 1021 | for (i = 0; i < numberIntegers; i++) { |
---|
| 1022 | int iColumn = integerVariable[i]; |
---|
| 1023 | saveLower[i] = COIN_DBL_MAX; |
---|
| 1024 | saveUpper[i] = -COIN_DBL_MAX; |
---|
| 1025 | if (solution[iColumn] < lower[iColumn] + primalTolerance) { |
---|
| 1026 | saveUpper[i] = upper[iColumn]; |
---|
| 1027 | solver->setColUpper(iColumn, lower[iColumn]); |
---|
| 1028 | } else if (solution[iColumn] > upper[iColumn] - primalTolerance) { |
---|
| 1029 | saveLower[i] = lower[iColumn]; |
---|
| 1030 | solver->setColLower(iColumn, upper[iColumn]); |
---|
| 1031 | } |
---|
| 1032 | } |
---|
| 1033 | solver->resolve(); |
---|
| 1034 | if (!solver->isProvenOptimal()) { |
---|
| 1035 | // presumably max time or some such |
---|
| 1036 | exitAll = true; |
---|
| 1037 | break; |
---|
| 1038 | } |
---|
| 1039 | for (i = 0; i < numberIntegers; i++) { |
---|
| 1040 | int iColumn = integerVariable[i]; |
---|
| 1041 | if (saveLower[i] != COIN_DBL_MAX) |
---|
| 1042 | solver->setColLower(iColumn, saveLower[i]); |
---|
| 1043 | if (saveUpper[i] != -COIN_DBL_MAX) |
---|
| 1044 | solver->setColUpper(iColumn, saveUpper[i]); |
---|
| 1045 | saveUpper[i] = -COIN_DBL_MAX; |
---|
| 1046 | } |
---|
| 1047 | memcpy(newSolution, solution, numberColumns*sizeof(double)); |
---|
| 1048 | int flip; |
---|
| 1049 | returnCode = rounds(solver, newSolution,/*saveObjective,*/ |
---|
| 1050 | numberIntegers, integerVariable, |
---|
| 1051 | /*pumpPrint,*/numberPasses, |
---|
| 1052 | /*roundExpensive_,*/defaultRounding_, &flip); |
---|
| 1053 | numberPasses++; |
---|
| 1054 | if (returnCode) { |
---|
| 1055 | // solution - but may not be better |
---|
| 1056 | // Compute using dot product |
---|
| 1057 | double newSolutionValue = -saveOffset; |
---|
| 1058 | for ( i = 0 ; i < numberColumns ; i++ ) |
---|
| 1059 | newSolutionValue += saveObjective[i] * newSolution[i]; |
---|
| 1060 | newSolutionValue *= direction; |
---|
| 1061 | sprintf(pumpPrint, "Intermediate solution found of %g", newSolutionValue); |
---|
| 1062 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1063 | << pumpPrint |
---|
| 1064 | << CoinMessageEol; |
---|
| 1065 | if (newSolutionValue < solutionValue) { |
---|
| 1066 | memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); |
---|
| 1067 | CoinWarmStartBasis * basis = |
---|
| 1068 | dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; |
---|
| 1069 | solutionFound = true; |
---|
| 1070 | if (exitNow(newSolutionValue)) |
---|
| 1071 | exitAll = true; |
---|
| 1072 | if (basis) { |
---|
| 1073 | bestBasis = * basis; |
---|
| 1074 | delete basis; |
---|
| 1075 | int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution); |
---|
| 1076 | if (!action) { |
---|
| 1077 | double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns); |
---|
| 1078 | double saveObjectiveValue = model_->getMinimizationObjValue(); |
---|
| 1079 | model_->setBestSolution(betterSolution, numberColumns, newSolutionValue); |
---|
| 1080 | if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue()) |
---|
| 1081 | model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue); |
---|
| 1082 | delete [] saveOldSolution; |
---|
| 1083 | } |
---|
| 1084 | if (!action || model_->maximumSecondsReached()) { |
---|
| 1085 | exitAll = true; // exit |
---|
| 1086 | break; |
---|
| 1087 | } |
---|
| 1088 | } |
---|
| 1089 | if ((accumulate_&1) != 0) { |
---|
| 1090 | model_->incrementUsed(betterSolution); // for local search |
---|
| 1091 | numberSolutions++; |
---|
| 1092 | } |
---|
| 1093 | solutionValue = newSolutionValue; |
---|
| 1094 | solutionFound = true; |
---|
| 1095 | if (exitNow(newSolutionValue)) |
---|
| 1096 | exitAll = true; |
---|
| 1097 | } else { |
---|
| 1098 | returnCode = 0; |
---|
| 1099 | } |
---|
| 1100 | } |
---|
| 1101 | } |
---|
| 1102 | int numberIterations = 0; |
---|
| 1103 | if (!doGeneral) { |
---|
| 1104 | // faster to do from all slack!!!! |
---|
| 1105 | if (allSlack) { |
---|
| 1106 | CoinWarmStartBasis dummy; |
---|
| 1107 | solver->setWarmStart(&dummy); |
---|
| 1108 | } |
---|
[750] | 1109 | #ifdef COIN_DEVELOP |
---|
[1286] | 1110 | printf("%d perturbed out of %d columns (%d changed)\n", numberPerturbed, numberColumns, numberChanged); |
---|
[750] | 1111 | #endif |
---|
[1286] | 1112 | bool takeHint; |
---|
| 1113 | OsiHintStrength strength; |
---|
| 1114 | solver->getHintParam(OsiDoDualInResolve, takeHint, strength); |
---|
| 1115 | if (dualPass && numberChanged > 2) { |
---|
| 1116 | solver->setHintParam(OsiDoDualInResolve, true); // dual may be better |
---|
| 1117 | if (dualPass == 1 && 2*numberChanged < numberColumns && |
---|
| 1118 | (numberChanged < 5000 || 6*numberChanged < numberColumns)) { |
---|
| 1119 | // but we need to make infeasible |
---|
| 1120 | CoinWarmStartBasis * basis = |
---|
| 1121 | dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; |
---|
| 1122 | if (basis) { |
---|
| 1123 | // modify |
---|
| 1124 | const double * lower = solver->getColLower(); |
---|
| 1125 | const double * upper = solver->getColUpper(); |
---|
| 1126 | double * solution = CoinCopyOfArray(solver->getColSolution(), |
---|
| 1127 | numberColumns); |
---|
| 1128 | const double * objective = solver->getObjCoefficients(); |
---|
| 1129 | int nChanged = 0; |
---|
| 1130 | for (i = 0; i < numberIntegersOrig; i++) { |
---|
| 1131 | int iColumn = integerVariableOrig[i]; |
---|
[1271] | 1132 | #ifdef RAND_RAND |
---|
[1286] | 1133 | if (nChanged > numberChanged) |
---|
| 1134 | break; |
---|
[1271] | 1135 | #endif |
---|
[1286] | 1136 | if (objective[iColumn] > 0.0) { |
---|
| 1137 | if (basis->getStructStatus(iColumn) == |
---|
| 1138 | CoinWarmStartBasis::atUpperBound) { |
---|
| 1139 | solution[iColumn] = lower[iColumn]; |
---|
| 1140 | basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound); |
---|
| 1141 | nChanged++; |
---|
| 1142 | } |
---|
| 1143 | } else if (objective[iColumn] < 0.0) { |
---|
| 1144 | if (basis->getStructStatus(iColumn) == |
---|
| 1145 | CoinWarmStartBasis::atLowerBound) { |
---|
| 1146 | solution[iColumn] = upper[iColumn]; |
---|
| 1147 | basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound); |
---|
| 1148 | nChanged++; |
---|
| 1149 | } |
---|
| 1150 | } |
---|
| 1151 | } |
---|
| 1152 | if (!nChanged) { |
---|
| 1153 | for (i = 0; i < numberIntegersOrig; i++) { |
---|
| 1154 | int iColumn = integerVariableOrig[i]; |
---|
| 1155 | if (objective[iColumn] > 0.0) { |
---|
| 1156 | if (basis->getStructStatus(iColumn) == |
---|
| 1157 | CoinWarmStartBasis::basic) { |
---|
| 1158 | solution[iColumn] = lower[iColumn]; |
---|
| 1159 | basis->setStructStatus(iColumn, CoinWarmStartBasis::atLowerBound); |
---|
| 1160 | break; |
---|
| 1161 | } |
---|
| 1162 | } else if (objective[iColumn] < 0.0) { |
---|
| 1163 | if (basis->getStructStatus(iColumn) == |
---|
| 1164 | CoinWarmStartBasis::basic) { |
---|
| 1165 | solution[iColumn] = upper[iColumn]; |
---|
| 1166 | basis->setStructStatus(iColumn, CoinWarmStartBasis::atUpperBound); |
---|
| 1167 | break; |
---|
| 1168 | } |
---|
| 1169 | } |
---|
| 1170 | } |
---|
| 1171 | } |
---|
| 1172 | solver->setColSolution(solution); |
---|
| 1173 | delete [] solution; |
---|
| 1174 | solver->setWarmStart(basis); |
---|
| 1175 | delete basis; |
---|
| 1176 | } |
---|
| 1177 | } else { |
---|
| 1178 | // faster to do from all slack!!!! ??? |
---|
| 1179 | CoinWarmStartBasis dummy; |
---|
| 1180 | solver->setWarmStart(&dummy); |
---|
| 1181 | } |
---|
| 1182 | } |
---|
| 1183 | if (numberTries > 1 && numberPasses == 1 && firstPerturbedObjective) { |
---|
| 1184 | // Modify to use convex combination |
---|
| 1185 | // use basis from first time |
---|
| 1186 | solver->setWarmStart(&saveBasis); |
---|
| 1187 | // and objective |
---|
| 1188 | if (secondPassOpt < 3 || (secondPassOpt >= 4 && secondPassOpt < 6)) |
---|
| 1189 | solver->setObjective(firstPerturbedObjective); |
---|
| 1190 | // and solution |
---|
| 1191 | solver->setColSolution(firstPerturbedSolution); |
---|
| 1192 | //if (secondPassOpt==2||secondPassOpt==5|| |
---|
| 1193 | if (firstPerturbedValue > cutoff) |
---|
| 1194 | solver->setHintParam(OsiDoDualInResolve, true); // dual may be better |
---|
| 1195 | } |
---|
| 1196 | solver->resolve(); |
---|
| 1197 | if (!solver->isProvenOptimal()) { |
---|
| 1198 | // presumably max time or some such |
---|
| 1199 | exitAll = true; |
---|
| 1200 | break; |
---|
| 1201 | } |
---|
| 1202 | solver->setHintParam(OsiDoDualInResolve, takeHint); |
---|
| 1203 | newTrueSolutionValue = -saveOffset; |
---|
| 1204 | newSumInfeas = 0.0; |
---|
| 1205 | newNumberInfeas = 0; |
---|
| 1206 | { |
---|
| 1207 | const double * newSolution = solver->getColSolution(); |
---|
| 1208 | for ( i = 0 ; i < numberColumns ; i++ ) { |
---|
| 1209 | if (solver->isInteger(i)) { |
---|
| 1210 | double value = newSolution[i]; |
---|
| 1211 | double nearest = floor(value + 0.5); |
---|
| 1212 | newSumInfeas += fabs(value - nearest); |
---|
| 1213 | if (fabs(value - nearest) > 1.0e-6) |
---|
| 1214 | newNumberInfeas++; |
---|
| 1215 | } |
---|
| 1216 | newTrueSolutionValue += saveObjective[i] * newSolution[i]; |
---|
| 1217 | } |
---|
| 1218 | newTrueSolutionValue *= direction; |
---|
| 1219 | if (numberPasses == 1 && secondPassOpt) { |
---|
| 1220 | if (numberTries == 1 || secondPassOpt > 3) { |
---|
| 1221 | // save basis |
---|
| 1222 | CoinWarmStartBasis * basis = |
---|
| 1223 | dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; |
---|
| 1224 | if (basis) { |
---|
| 1225 | saveBasis = * basis; |
---|
| 1226 | delete basis; |
---|
| 1227 | } |
---|
| 1228 | delete [] firstPerturbedObjective; |
---|
| 1229 | delete [] firstPerturbedSolution; |
---|
| 1230 | firstPerturbedObjective = CoinCopyOfArray(solver->getObjCoefficients(), numberColumns); |
---|
| 1231 | firstPerturbedSolution = CoinCopyOfArray(solver->getColSolution(), numberColumns); |
---|
| 1232 | firstPerturbedValue = newTrueSolutionValue; |
---|
| 1233 | } |
---|
| 1234 | } |
---|
| 1235 | if (newNumberInfeas && newNumberInfeas < 15) { |
---|
[1393] | 1236 | #ifdef JJF_ZERO |
---|
[1286] | 1237 | roundingObjective = solutionValue; |
---|
| 1238 | OsiSolverInterface * saveSolver = model_->swapSolver(solver); |
---|
| 1239 | double * currentObjective = |
---|
| 1240 | CoinCopyOfArray(solver->getObjCoefficients(), numberColumns); |
---|
| 1241 | solver->setObjective(saveObjective); |
---|
| 1242 | double saveOffset2; |
---|
| 1243 | solver->getDblParam(OsiObjOffset, saveOffset2); |
---|
| 1244 | solver->setDblParam(OsiObjOffset, saveOffset); |
---|
| 1245 | int ifSol = roundingHeuristic.solution(roundingObjective, roundingSolution); |
---|
| 1246 | solver->setObjective(currentObjective); |
---|
| 1247 | solver->setDblParam(OsiObjOffset, saveOffset2); |
---|
| 1248 | delete [] currentObjective; |
---|
| 1249 | model_->swapSolver(saveSolver); |
---|
| 1250 | if (ifSol > 0) |
---|
| 1251 | abort(); |
---|
[1013] | 1252 | #endif |
---|
[1286] | 1253 | int numberRows = solver->getNumRows(); |
---|
| 1254 | double * rowActivity = new double[numberRows]; |
---|
| 1255 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
| 1256 | int * which = new int[newNumberInfeas]; |
---|
| 1257 | int * stack = new int[newNumberInfeas+1]; |
---|
| 1258 | double * baseValue = new double[newNumberInfeas]; |
---|
| 1259 | int * whichRow = new int[numberRows]; |
---|
| 1260 | double * rowValue = new double[numberRows]; |
---|
| 1261 | memset(rowValue, 0, numberRows*sizeof(double)); |
---|
| 1262 | int nRow = 0; |
---|
| 1263 | // Column copy |
---|
| 1264 | const double * element = solver->getMatrixByCol()->getElements(); |
---|
| 1265 | const int * row = solver->getMatrixByCol()->getIndices(); |
---|
| 1266 | const CoinBigIndex * columnStart = solver->getMatrixByCol()->getVectorStarts(); |
---|
| 1267 | const int * columnLength = solver->getMatrixByCol()->getVectorLengths(); |
---|
| 1268 | int n = 0; |
---|
| 1269 | double contrib = 0.0; |
---|
| 1270 | for ( i = 0 ; i < numberColumns ; i++ ) { |
---|
| 1271 | double value = newSolution[i]; |
---|
| 1272 | if (solver->isInteger(i)) { |
---|
| 1273 | double nearest = floor(value + 0.5); |
---|
| 1274 | if (fabs(value - nearest) > 1.0e-6) { |
---|
| 1275 | //printf("Column %d value %g\n",i,value); |
---|
| 1276 | for (CoinBigIndex j = columnStart[i]; |
---|
| 1277 | j < columnStart[i] + columnLength[i]; j++) { |
---|
| 1278 | int iRow = row[j]; |
---|
| 1279 | //printf("row %d element %g\n",iRow,element[j]); |
---|
| 1280 | if (!rowValue[iRow]) { |
---|
| 1281 | rowValue[iRow] = 1.0; |
---|
| 1282 | whichRow[nRow++] = iRow; |
---|
| 1283 | } |
---|
| 1284 | } |
---|
| 1285 | baseValue[n] = floor(value); |
---|
| 1286 | contrib += saveObjective[i] * value; |
---|
| 1287 | value = 0.0; |
---|
| 1288 | stack[n] = 0; |
---|
| 1289 | which[n++] = i; |
---|
| 1290 | } |
---|
| 1291 | } |
---|
| 1292 | for (CoinBigIndex j = columnStart[i]; |
---|
| 1293 | j < columnStart[i] + columnLength[i]; j++) { |
---|
| 1294 | int iRow = row[j]; |
---|
| 1295 | rowActivity[iRow] += value * element[j]; |
---|
| 1296 | } |
---|
| 1297 | } |
---|
| 1298 | if (newNumberInfeas < 15) { |
---|
| 1299 | stack[n] = newNumberInfeas + 100; |
---|
| 1300 | int iStack = n; |
---|
| 1301 | memset(rowValue, 0, numberRows*sizeof(double)); |
---|
| 1302 | const double * rowLower = solver->getRowLower(); |
---|
| 1303 | const double * rowUpper = solver->getRowUpper(); |
---|
| 1304 | while (iStack >= 0) { |
---|
| 1305 | double contrib2 = 0.0; |
---|
| 1306 | // Could do faster |
---|
| 1307 | for (int k = 0 ; k < n ; k++ ) { |
---|
| 1308 | i = which[k]; |
---|
| 1309 | double value = baseValue[k] + stack[k]; |
---|
| 1310 | contrib2 += saveObjective[i] * value; |
---|
| 1311 | for (CoinBigIndex j = columnStart[i]; |
---|
| 1312 | j < columnStart[i] + columnLength[i]; j++) { |
---|
| 1313 | int iRow = row[j]; |
---|
| 1314 | rowValue[iRow] += value * element[j]; |
---|
| 1315 | } |
---|
| 1316 | } |
---|
| 1317 | // check if feasible |
---|
| 1318 | bool feasible = true; |
---|
| 1319 | for (int k = 0; k < nRow; k++) { |
---|
| 1320 | i = whichRow[k]; |
---|
| 1321 | double value = rowValue[i] + rowActivity[i]; |
---|
| 1322 | rowValue[i] = 0.0; |
---|
| 1323 | if (value < rowLower[i] - 1.0e-7 || |
---|
| 1324 | value > rowUpper[i] + 1.0e-7) |
---|
| 1325 | feasible = false; |
---|
| 1326 | } |
---|
| 1327 | if (feasible) { |
---|
| 1328 | double newObj = newTrueSolutionValue * direction; |
---|
| 1329 | newObj += contrib2 - contrib; |
---|
| 1330 | newObj *= direction; |
---|
[1127] | 1331 | #ifdef COIN_DEVELOP |
---|
[1286] | 1332 | printf("FFFeasible! - obj %g\n", newObj); |
---|
[1127] | 1333 | #endif |
---|
[1286] | 1334 | if (newObj < roundingObjective - 1.0e-6) { |
---|
[1127] | 1335 | #ifdef COIN_DEVELOP |
---|
[1286] | 1336 | printf("FBetter\n"); |
---|
[1127] | 1337 | #endif |
---|
[1286] | 1338 | roundingObjective = newObj; |
---|
| 1339 | memcpy(roundingSolution, newSolution, numberColumns*sizeof(double)); |
---|
| 1340 | for (int k = 0 ; k < n ; k++ ) { |
---|
| 1341 | i = which[k]; |
---|
| 1342 | double value = baseValue[k] + stack[k]; |
---|
| 1343 | roundingSolution[i] = value; |
---|
| 1344 | } |
---|
| 1345 | } |
---|
| 1346 | } |
---|
| 1347 | while (iStack >= 0 && stack[iStack]) { |
---|
| 1348 | stack[iStack]--; |
---|
| 1349 | iStack--; |
---|
| 1350 | } |
---|
| 1351 | if (iStack >= 0) { |
---|
| 1352 | stack[iStack] = 1; |
---|
| 1353 | iStack = n; |
---|
| 1354 | stack[n] = 1; |
---|
| 1355 | } |
---|
| 1356 | } |
---|
| 1357 | } |
---|
| 1358 | delete [] rowActivity; |
---|
| 1359 | delete [] which; |
---|
| 1360 | delete [] stack; |
---|
| 1361 | delete [] baseValue; |
---|
| 1362 | delete [] whichRow; |
---|
| 1363 | delete [] rowValue; |
---|
| 1364 | } |
---|
| 1365 | } |
---|
| 1366 | if (true) { |
---|
| 1367 | OsiSolverInterface * saveSolver = model_->swapSolver(clonedSolver); |
---|
| 1368 | clonedSolver->setColSolution(solver->getColSolution()); |
---|
| 1369 | CbcRounding heuristic1(*model_); |
---|
| 1370 | heuristic1.setHeuristicName("rounding in feaspump!"); |
---|
| 1371 | heuristic1.setWhen(1); |
---|
| 1372 | roundingObjective = CoinMin(roundingObjective, solutionValue); |
---|
| 1373 | double testSolutionValue = newTrueSolutionValue; |
---|
| 1374 | int returnCode = heuristic1.solution(roundingObjective, |
---|
| 1375 | roundingSolution, |
---|
| 1376 | testSolutionValue) ; |
---|
| 1377 | if (returnCode == 1) { |
---|
[1127] | 1378 | #ifdef COIN_DEVELOP |
---|
[1286] | 1379 | printf("rounding obj of %g?\n", roundingObjective); |
---|
[1127] | 1380 | #endif |
---|
[1286] | 1381 | //roundingObjective = newSolutionValue; |
---|
| 1382 | //} else { |
---|
| 1383 | //roundingObjective = COIN_DBL_MAX; |
---|
| 1384 | } |
---|
| 1385 | model_->swapSolver(saveSolver); |
---|
| 1386 | } |
---|
| 1387 | if (!solver->isProvenOptimal()) { |
---|
| 1388 | // presumably max time or some such |
---|
| 1389 | exitAll = true; |
---|
| 1390 | break; |
---|
| 1391 | } |
---|
| 1392 | // in case very dubious solver |
---|
| 1393 | lower = solver->getColLower(); |
---|
| 1394 | upper = solver->getColUpper(); |
---|
| 1395 | solution = solver->getColSolution(); |
---|
| 1396 | numberIterations = solver->getIterationCount(); |
---|
| 1397 | } else { |
---|
| 1398 | int * addStart = new int[2*general+1]; |
---|
| 1399 | int * addIndex = new int[4*general]; |
---|
| 1400 | double * addElement = new double[4*general]; |
---|
| 1401 | double * addLower = new double[2*general]; |
---|
| 1402 | double * addUpper = new double[2*general]; |
---|
| 1403 | double * obj = new double[general]; |
---|
| 1404 | int nAdd = 0; |
---|
| 1405 | for (i = 0; i < numberIntegers; i++) { |
---|
| 1406 | int iColumn = integerVariable[i]; |
---|
| 1407 | if (newSolution[iColumn] > lower[iColumn] + primalTolerance && |
---|
| 1408 | newSolution[iColumn] < upper[iColumn] - primalTolerance) { |
---|
| 1409 | assert (upper[iColumn] - lower[iColumn] > 1.00001); |
---|
| 1410 | obj[nAdd] = 1.0; |
---|
| 1411 | addLower[nAdd] = 0.0; |
---|
| 1412 | addUpper[nAdd] = COIN_DBL_MAX; |
---|
| 1413 | nAdd++; |
---|
| 1414 | } |
---|
| 1415 | } |
---|
| 1416 | OsiSolverInterface * solver2 = solver; |
---|
| 1417 | if (nAdd) { |
---|
| 1418 | CoinZeroN(addStart, nAdd + 1); |
---|
| 1419 | solver2 = solver->clone(); |
---|
| 1420 | solver2->addCols(nAdd, addStart, NULL, NULL, addLower, addUpper, obj); |
---|
| 1421 | // feasible solution |
---|
| 1422 | double * sol = new double[nAdd+numberColumns]; |
---|
| 1423 | memcpy(sol, solution, numberColumns*sizeof(double)); |
---|
| 1424 | // now rows |
---|
| 1425 | int nAdd = 0; |
---|
| 1426 | int nEl = 0; |
---|
| 1427 | int nAddRow = 0; |
---|
| 1428 | for (i = 0; i < numberIntegers; i++) { |
---|
| 1429 | int iColumn = integerVariable[i]; |
---|
| 1430 | if (newSolution[iColumn] > lower[iColumn] + primalTolerance && |
---|
| 1431 | newSolution[iColumn] < upper[iColumn] - primalTolerance) { |
---|
| 1432 | addLower[nAddRow] = -newSolution[iColumn];; |
---|
| 1433 | addUpper[nAddRow] = COIN_DBL_MAX; |
---|
| 1434 | addIndex[nEl] = iColumn; |
---|
| 1435 | addElement[nEl++] = -1.0; |
---|
| 1436 | addIndex[nEl] = numberColumns + nAdd; |
---|
| 1437 | addElement[nEl++] = 1.0; |
---|
| 1438 | nAddRow++; |
---|
| 1439 | addStart[nAddRow] = nEl; |
---|
| 1440 | addLower[nAddRow] = newSolution[iColumn];; |
---|
| 1441 | addUpper[nAddRow] = COIN_DBL_MAX; |
---|
| 1442 | addIndex[nEl] = iColumn; |
---|
| 1443 | addElement[nEl++] = 1.0; |
---|
| 1444 | addIndex[nEl] = numberColumns + nAdd; |
---|
| 1445 | addElement[nEl++] = 1.0; |
---|
| 1446 | nAddRow++; |
---|
| 1447 | addStart[nAddRow] = nEl; |
---|
| 1448 | sol[nAdd+numberColumns] = fabs(sol[iColumn] - newSolution[iColumn]); |
---|
| 1449 | nAdd++; |
---|
| 1450 | } |
---|
| 1451 | } |
---|
| 1452 | solver2->setColSolution(sol); |
---|
| 1453 | delete [] sol; |
---|
| 1454 | solver2->addRows(nAddRow, addStart, addIndex, addElement, addLower, addUpper); |
---|
| 1455 | } |
---|
| 1456 | delete [] addStart; |
---|
| 1457 | delete [] addIndex; |
---|
| 1458 | delete [] addElement; |
---|
| 1459 | delete [] addLower; |
---|
| 1460 | delete [] addUpper; |
---|
| 1461 | delete [] obj; |
---|
| 1462 | solver2->resolve(); |
---|
| 1463 | if (!solver2->isProvenOptimal()) { |
---|
| 1464 | // presumably max time or some such |
---|
| 1465 | exitAll = true; |
---|
| 1466 | break; |
---|
| 1467 | } |
---|
| 1468 | //assert (solver2->isProvenOptimal()); |
---|
| 1469 | if (nAdd) { |
---|
| 1470 | solver->setColSolution(solver2->getColSolution()); |
---|
| 1471 | numberIterations = solver2->getIterationCount(); |
---|
| 1472 | delete solver2; |
---|
| 1473 | } else { |
---|
| 1474 | numberIterations = solver->getIterationCount(); |
---|
| 1475 | } |
---|
| 1476 | lower = solver->getColLower(); |
---|
| 1477 | upper = solver->getColUpper(); |
---|
| 1478 | solution = solver->getColSolution(); |
---|
| 1479 | newTrueSolutionValue = -saveOffset; |
---|
| 1480 | newSumInfeas = 0.0; |
---|
| 1481 | newNumberInfeas = 0; |
---|
| 1482 | { |
---|
| 1483 | const double * newSolution = solver->getColSolution(); |
---|
| 1484 | for ( i = 0 ; i < numberColumns ; i++ ) { |
---|
| 1485 | if (solver->isInteger(i)) { |
---|
| 1486 | double value = newSolution[i]; |
---|
| 1487 | double nearest = floor(value + 0.5); |
---|
| 1488 | newSumInfeas += fabs(value - nearest); |
---|
| 1489 | if (fabs(value - nearest) > 1.0e-6) |
---|
| 1490 | newNumberInfeas++; |
---|
| 1491 | } |
---|
| 1492 | newTrueSolutionValue += saveObjective[i] * newSolution[i]; |
---|
| 1493 | } |
---|
| 1494 | newTrueSolutionValue *= direction; |
---|
| 1495 | } |
---|
| 1496 | } |
---|
| 1497 | if (lastMove != 1000000) { |
---|
| 1498 | if (newSumInfeas < lastSumInfeas) { |
---|
| 1499 | lastMove = numberPasses; |
---|
| 1500 | lastSumInfeas = newSumInfeas; |
---|
| 1501 | } else if (newSumInfeas > lastSumInfeas + 1.0e-5) { |
---|
| 1502 | lastMove = 1000000; // going up |
---|
| 1503 | } |
---|
| 1504 | } |
---|
| 1505 | totalNumberIterations += numberIterations; |
---|
| 1506 | if (solver->getNumRows() < 3000) |
---|
| 1507 | sprintf(pumpPrint, "Pass %3d: suminf. %10.5f (%d) obj. %g iterations %d", |
---|
| 1508 | numberPasses + totalNumberPasses, |
---|
| 1509 | newSumInfeas, newNumberInfeas, |
---|
| 1510 | newTrueSolutionValue, numberIterations); |
---|
| 1511 | else |
---|
| 1512 | sprintf(pumpPrint, "Pass %3d: (%.2f seconds) suminf. %10.5f (%d) obj. %g iterations %d", numberPasses + totalNumberPasses, |
---|
| 1513 | model_->getCurrentSeconds(), newSumInfeas, newNumberInfeas, |
---|
| 1514 | newTrueSolutionValue, numberIterations); |
---|
| 1515 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1516 | << pumpPrint |
---|
| 1517 | << CoinMessageEol; |
---|
| 1518 | if (closestSolution && solver->getObjValue() < closestObjectiveValue) { |
---|
| 1519 | int i; |
---|
| 1520 | const double * objective = solver->getObjCoefficients(); |
---|
| 1521 | for (i = 0; i < numberIntegersOrig; i++) { |
---|
| 1522 | int iColumn = integerVariableOrig[i]; |
---|
| 1523 | if (objective[iColumn] > 0.0) |
---|
| 1524 | closestSolution[i] = 0; |
---|
| 1525 | else |
---|
| 1526 | closestSolution[i] = 1; |
---|
| 1527 | } |
---|
| 1528 | closestObjectiveValue = solver->getObjValue(); |
---|
| 1529 | } |
---|
| 1530 | // See if we need to think about changing rhs |
---|
| 1531 | if ((switches_&12) != 0 && useRhs < 1.0e50) { |
---|
| 1532 | double oldRhs = useRhs; |
---|
| 1533 | bool trying = false; |
---|
| 1534 | if ((switches_&4) != 0 && numberPasses && (numberPasses % 50) == 0) { |
---|
| 1535 | if (solutionValue > 1.0e20) { |
---|
| 1536 | // only if no genuine solution |
---|
| 1537 | double gap = useRhs - continuousObjectiveValue; |
---|
| 1538 | useRhs += 0.1 * gap; |
---|
| 1539 | if (exactMultiple) { |
---|
| 1540 | useRhs = exactMultiple * ceil(useRhs / exactMultiple); |
---|
| 1541 | useRhs = CoinMax(useRhs, oldRhs + exactMultiple); |
---|
| 1542 | } |
---|
| 1543 | trying = true; |
---|
| 1544 | } |
---|
| 1545 | } |
---|
| 1546 | if ((switches_&8) != 0) { |
---|
| 1547 | // Put in new suminf and check |
---|
| 1548 | double largest = newSumInfeas; |
---|
| 1549 | double smallest = newSumInfeas; |
---|
| 1550 | for (int i = 0; i < SIZE_BOBBLE - 1; i++) { |
---|
| 1551 | double value = saveSumInf[i+1]; |
---|
| 1552 | saveSumInf[i] = value; |
---|
| 1553 | largest = CoinMax(largest, value); |
---|
| 1554 | smallest = CoinMin(smallest, value); |
---|
| 1555 | } |
---|
| 1556 | saveSumInf[SIZE_BOBBLE-1] = newSumInfeas; |
---|
| 1557 | if (smallest*1.5 > largest && smallest > 2.0) { |
---|
| 1558 | if (bobbleMode == 0) { |
---|
| 1559 | // go closer |
---|
| 1560 | double gap = oldRhs - continuousObjectiveValue; |
---|
| 1561 | useRhs -= 0.4 * gap; |
---|
| 1562 | if (exactMultiple) { |
---|
| 1563 | double value = floor(useRhs / exactMultiple); |
---|
| 1564 | useRhs = CoinMin(value * exactMultiple, oldRhs - exactMultiple); |
---|
| 1565 | } |
---|
| 1566 | if (useRhs < continuousObjectiveValue) { |
---|
| 1567 | // skip decrease |
---|
| 1568 | bobbleMode = 1; |
---|
| 1569 | useRhs = oldRhs; |
---|
| 1570 | } |
---|
| 1571 | } |
---|
| 1572 | if (bobbleMode) { |
---|
| 1573 | trying = true; |
---|
| 1574 | // weaken |
---|
| 1575 | if (solutionValue < 1.0e20) { |
---|
| 1576 | double gap = solutionValue - oldRhs; |
---|
| 1577 | useRhs += 0.3 * gap; |
---|
| 1578 | } else { |
---|
| 1579 | double gap = oldRhs - continuousObjectiveValue; |
---|
| 1580 | useRhs += 0.05 * gap; |
---|
| 1581 | } |
---|
| 1582 | if (exactMultiple) { |
---|
| 1583 | double value = ceil(useRhs / exactMultiple); |
---|
| 1584 | useRhs = CoinMin(value * exactMultiple, |
---|
| 1585 | solutionValue - exactMultiple); |
---|
| 1586 | } |
---|
| 1587 | } |
---|
| 1588 | bobbleMode++; |
---|
| 1589 | // reset |
---|
| 1590 | CoinFillN(saveSumInf, SIZE_BOBBLE, COIN_DBL_MAX); |
---|
| 1591 | } |
---|
| 1592 | } |
---|
| 1593 | if (useRhs != oldRhs) { |
---|
| 1594 | // tidy up |
---|
| 1595 | if (exactMultiple) { |
---|
| 1596 | double value = floor(useRhs / exactMultiple); |
---|
| 1597 | double bestPossible = ceil(continuousObjectiveValue / exactMultiple); |
---|
| 1598 | useRhs = CoinMax(value, bestPossible) * exactMultiple; |
---|
| 1599 | } else { |
---|
| 1600 | useRhs = CoinMax(useRhs, continuousObjectiveValue); |
---|
| 1601 | } |
---|
| 1602 | int k = solver->getNumRows() - 1; |
---|
| 1603 | solver->setRowUpper(k, useRhs + useOffset); |
---|
| 1604 | bool takeHint; |
---|
| 1605 | OsiHintStrength strength; |
---|
| 1606 | solver->getHintParam(OsiDoDualInResolve, takeHint, strength); |
---|
| 1607 | if (useRhs < oldRhs) { |
---|
| 1608 | solver->setHintParam(OsiDoDualInResolve, true); |
---|
| 1609 | solver->resolve(); |
---|
| 1610 | } else if (useRhs > oldRhs) { |
---|
| 1611 | solver->setHintParam(OsiDoDualInResolve, false); |
---|
| 1612 | solver->resolve(); |
---|
| 1613 | } |
---|
| 1614 | solver->setHintParam(OsiDoDualInResolve, takeHint); |
---|
| 1615 | if (!solver->isProvenOptimal()) { |
---|
| 1616 | // presumably max time or some such |
---|
| 1617 | exitAll = true; |
---|
| 1618 | break; |
---|
| 1619 | } |
---|
| 1620 | } else if (trying) { |
---|
| 1621 | // doesn't look good |
---|
| 1622 | break; |
---|
| 1623 | } |
---|
| 1624 | } |
---|
| 1625 | } |
---|
| 1626 | // reduce scale factor |
---|
| 1627 | scaleFactor *= weightFactor_; |
---|
| 1628 | } // END WHILE |
---|
| 1629 | // see if rounding worked! |
---|
| 1630 | if (roundingObjective < solutionValue) { |
---|
| 1631 | if (roundingObjective < solutionValue - 1.0e-6*fabs(roundingObjective)) { |
---|
| 1632 | sprintf(pumpPrint, "Rounding solution of %g is better than previous of %g\n", |
---|
| 1633 | roundingObjective, solutionValue); |
---|
| 1634 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1635 | << pumpPrint |
---|
| 1636 | << CoinMessageEol; |
---|
| 1637 | } |
---|
| 1638 | solutionValue = roundingObjective; |
---|
| 1639 | newSolutionValue = solutionValue; |
---|
| 1640 | memcpy(betterSolution, roundingSolution, numberColumns*sizeof(double)); |
---|
| 1641 | solutionFound = true; |
---|
| 1642 | if (exitNow(roundingObjective)) |
---|
| 1643 | exitAll = true; |
---|
| 1644 | } |
---|
| 1645 | if (!solutionFound) { |
---|
| 1646 | sprintf(pumpPrint, "No solution found this major pass"); |
---|
| 1647 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1648 | << pumpPrint |
---|
| 1649 | << CoinMessageEol; |
---|
| 1650 | } |
---|
| 1651 | //} |
---|
| 1652 | delete solver; |
---|
| 1653 | solver = NULL; |
---|
| 1654 | for ( j = 0; j < NUMBER_OLD; j++) |
---|
| 1655 | delete [] oldSolution[j]; |
---|
| 1656 | delete [] oldSolution; |
---|
| 1657 | delete [] saveObjective; |
---|
| 1658 | if (usedColumn && !exitAll) { |
---|
| 1659 | OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone(); |
---|
| 1660 | const double * colLower = newSolver->getColLower(); |
---|
| 1661 | const double * colUpper = newSolver->getColUpper(); |
---|
| 1662 | bool stopBAB = false; |
---|
| 1663 | int allowedPass = -1; |
---|
| 1664 | if (maximumAllowed > 0) |
---|
| 1665 | allowedPass = CoinMax(numberPasses - maximumAllowed, -1); |
---|
| 1666 | while (!stopBAB) { |
---|
| 1667 | stopBAB = true; |
---|
| 1668 | int i; |
---|
| 1669 | int nFix = 0; |
---|
| 1670 | int nFixI = 0; |
---|
| 1671 | int nFixC = 0; |
---|
| 1672 | int nFixC2 = 0; |
---|
| 1673 | for (i = 0; i < numberIntegersOrig; i++) { |
---|
| 1674 | int iColumn = integerVariableOrig[i]; |
---|
| 1675 | //const OsiObject * object = model_->object(i); |
---|
| 1676 | //double originalLower; |
---|
| 1677 | //double originalUpper; |
---|
| 1678 | //getIntegerInformation( object,originalLower, originalUpper); |
---|
| 1679 | //assert(colLower[iColumn]==originalLower); |
---|
| 1680 | //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower)); |
---|
| 1681 | newSolver->setColLower(iColumn, colLower[iColumn]); |
---|
| 1682 | //assert(colUpper[iColumn]==originalUpper); |
---|
| 1683 | //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper)); |
---|
| 1684 | newSolver->setColUpper(iColumn, colUpper[iColumn]); |
---|
| 1685 | if (usedColumn[iColumn] <= allowedPass) { |
---|
| 1686 | double value = lastSolution[iColumn]; |
---|
| 1687 | double nearest = floor(value + 0.5); |
---|
| 1688 | if (fabs(value - nearest) < 1.0e-7) { |
---|
| 1689 | if (nearest == colLower[iColumn]) { |
---|
| 1690 | newSolver->setColUpper(iColumn, colLower[iColumn]); |
---|
| 1691 | nFix++; |
---|
| 1692 | } else if (nearest == colUpper[iColumn]) { |
---|
| 1693 | newSolver->setColLower(iColumn, colUpper[iColumn]); |
---|
| 1694 | nFix++; |
---|
| 1695 | } else if (fixInternal) { |
---|
| 1696 | newSolver->setColLower(iColumn, nearest); |
---|
| 1697 | newSolver->setColUpper(iColumn, nearest); |
---|
| 1698 | nFix++; |
---|
| 1699 | nFixI++; |
---|
| 1700 | } |
---|
| 1701 | } |
---|
| 1702 | } |
---|
| 1703 | } |
---|
| 1704 | if (fixContinuous) { |
---|
| 1705 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 1706 | if (!newSolver->isInteger(iColumn) && usedColumn[iColumn] <= allowedPass) { |
---|
| 1707 | double value = lastSolution[iColumn]; |
---|
| 1708 | if (value < colLower[iColumn] + 1.0e-8) { |
---|
| 1709 | newSolver->setColUpper(iColumn, colLower[iColumn]); |
---|
| 1710 | nFixC++; |
---|
| 1711 | } else if (value > colUpper[iColumn] - 1.0e-8) { |
---|
| 1712 | newSolver->setColLower(iColumn, colUpper[iColumn]); |
---|
| 1713 | nFixC++; |
---|
| 1714 | } else if (fixContinuous == 2) { |
---|
| 1715 | newSolver->setColLower(iColumn, value); |
---|
| 1716 | newSolver->setColUpper(iColumn, value); |
---|
| 1717 | nFixC++; |
---|
| 1718 | nFixC2++; |
---|
| 1719 | } |
---|
| 1720 | } |
---|
| 1721 | } |
---|
| 1722 | } |
---|
| 1723 | newSolver->initialSolve(); |
---|
| 1724 | if (!newSolver->isProvenOptimal()) { |
---|
| 1725 | //newSolver->writeMps("bad.mps"); |
---|
| 1726 | //assert (newSolver->isProvenOptimal()); |
---|
| 1727 | exitAll = true; |
---|
| 1728 | break; |
---|
| 1729 | } |
---|
| 1730 | sprintf(pumpPrint, "Before mini branch and bound, %d integers at bound fixed and %d continuous", |
---|
| 1731 | nFix, nFixC); |
---|
| 1732 | if (nFixC2 + nFixI != 0) |
---|
| 1733 | sprintf(pumpPrint + strlen(pumpPrint), " of which %d were internal integer and %d internal continuous", |
---|
| 1734 | nFixI, nFixC2); |
---|
| 1735 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1736 | << pumpPrint |
---|
| 1737 | << CoinMessageEol; |
---|
| 1738 | double saveValue = newSolutionValue; |
---|
| 1739 | if (newSolutionValue - model_->getCutoffIncrement() |
---|
| 1740 | > continuousObjectiveValue - 1.0e-7) { |
---|
| 1741 | double saveFraction = fractionSmall_; |
---|
| 1742 | if (numberTries > 1 && !numberBandBsolutions) |
---|
| 1743 | fractionSmall_ *= 0.5; |
---|
| 1744 | // Give branch and bound a bit more freedom |
---|
| 1745 | double cutoff2 = newSolutionValue + |
---|
| 1746 | CoinMax(model_->getCutoffIncrement(), 1.0e-3); |
---|
| 1747 | int returnCode2 = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue, |
---|
| 1748 | cutoff2, "CbcHeuristicLocalAfterFPump"); |
---|
| 1749 | fractionSmall_ = saveFraction; |
---|
| 1750 | if (returnCode2 < 0) { |
---|
| 1751 | if (returnCode2 == -2) { |
---|
| 1752 | exitAll = true; |
---|
| 1753 | returnCode = 0; |
---|
| 1754 | } else { |
---|
| 1755 | returnCode2 = 0; // returned on size - try changing |
---|
| 1756 | //#define ROUND_AGAIN |
---|
[1053] | 1757 | #ifdef ROUND_AGAIN |
---|
[1286] | 1758 | if (numberTries == 1 && numberPasses > 20 && allowedPass < numberPasses - 1) { |
---|
| 1759 | allowedPass = (numberPasses + allowedPass) >> 1; |
---|
| 1760 | sprintf(pumpPrint, |
---|
| 1761 | "Fixing all variables which were last changed on pass %d and trying again", |
---|
| 1762 | allowedPass); |
---|
| 1763 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1764 | << pumpPrint |
---|
| 1765 | << CoinMessageEol; |
---|
| 1766 | stopBAB = false; |
---|
| 1767 | continue; |
---|
| 1768 | } |
---|
[1053] | 1769 | #endif |
---|
[1286] | 1770 | } |
---|
| 1771 | } |
---|
| 1772 | if ((returnCode2&2) != 0) { |
---|
| 1773 | // could add cut |
---|
| 1774 | returnCode2 &= ~2; |
---|
| 1775 | } |
---|
| 1776 | if (returnCode2) |
---|
| 1777 | numberBandBsolutions++; |
---|
| 1778 | } else { |
---|
| 1779 | // no need |
---|
| 1780 | exitAll = true; |
---|
| 1781 | //returnCode=0; |
---|
| 1782 | } |
---|
| 1783 | // recompute solution value |
---|
| 1784 | if (returnCode && true) { |
---|
| 1785 | delete newSolver; |
---|
| 1786 | newSolver = cloneBut(3); // was model_->continuousSolver()->clone(); |
---|
| 1787 | newSolutionValue = -saveOffset; |
---|
| 1788 | double newSumInfeas = 0.0; |
---|
| 1789 | const double * obj = newSolver->getObjCoefficients(); |
---|
| 1790 | for (int i = 0 ; i < numberColumns ; i++ ) { |
---|
| 1791 | if (newSolver->isInteger(i)) { |
---|
| 1792 | double value = newSolution[i]; |
---|
| 1793 | double nearest = floor(value + 0.5); |
---|
| 1794 | newSumInfeas += fabs(value - nearest); |
---|
| 1795 | } |
---|
| 1796 | newSolutionValue += obj[i] * newSolution[i]; |
---|
| 1797 | } |
---|
| 1798 | newSolutionValue *= direction; |
---|
| 1799 | } |
---|
| 1800 | bool gotSolution = false; |
---|
| 1801 | if (returnCode && newSolutionValue < saveValue) { |
---|
| 1802 | sprintf(pumpPrint, "Mini branch and bound improved solution from %g to %g (%.2f seconds)", |
---|
| 1803 | saveValue, newSolutionValue, model_->getCurrentSeconds()); |
---|
| 1804 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1805 | << pumpPrint |
---|
| 1806 | << CoinMessageEol; |
---|
| 1807 | memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); |
---|
| 1808 | gotSolution = true; |
---|
| 1809 | if (fixContinuous && nFixC + nFixC2 > 0) { |
---|
| 1810 | // may be able to do even better |
---|
| 1811 | int nFixed = 0; |
---|
| 1812 | const double * lower = model_->solver()->getColLower(); |
---|
| 1813 | const double * upper = model_->solver()->getColUpper(); |
---|
| 1814 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 1815 | double value = newSolution[iColumn]; |
---|
| 1816 | if (newSolver->isInteger(iColumn)) { |
---|
| 1817 | value = floor(newSolution[iColumn] + 0.5); |
---|
| 1818 | newSolver->setColLower(iColumn, value); |
---|
| 1819 | newSolver->setColUpper(iColumn, value); |
---|
| 1820 | nFixed++; |
---|
| 1821 | } else { |
---|
| 1822 | newSolver->setColLower(iColumn, lower[iColumn]); |
---|
| 1823 | newSolver->setColUpper(iColumn, upper[iColumn]); |
---|
| 1824 | if (value < lower[iColumn]) |
---|
| 1825 | value = lower[iColumn]; |
---|
| 1826 | else if (value > upper[iColumn]) |
---|
| 1827 | value = upper[iColumn]; |
---|
| 1828 | |
---|
| 1829 | } |
---|
| 1830 | newSolution[iColumn] = value; |
---|
| 1831 | } |
---|
| 1832 | newSolver->setColSolution(newSolution); |
---|
| 1833 | //#define CLP_INVESTIGATE2 |
---|
[1271] | 1834 | #ifdef CLP_INVESTIGATE2 |
---|
[1286] | 1835 | { |
---|
| 1836 | // check |
---|
| 1837 | // get row activities |
---|
| 1838 | int numberRows = newSolver->getNumRows(); |
---|
| 1839 | double * rowActivity = new double[numberRows]; |
---|
| 1840 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
| 1841 | newSolver->getMatrixByCol()->times(newSolution, rowActivity) ; |
---|
| 1842 | double largestInfeasibility = primalTolerance; |
---|
| 1843 | double sumInfeasibility = 0.0; |
---|
| 1844 | int numberBadRows = 0; |
---|
| 1845 | const double * rowLower = newSolver->getRowLower(); |
---|
| 1846 | const double * rowUpper = newSolver->getRowUpper(); |
---|
| 1847 | for (i = 0 ; i < numberRows ; i++) { |
---|
| 1848 | double value; |
---|
| 1849 | value = rowLower[i] - rowActivity[i]; |
---|
| 1850 | if (value > primalTolerance) { |
---|
| 1851 | numberBadRows++; |
---|
| 1852 | largestInfeasibility = CoinMax(largestInfeasibility, value); |
---|
| 1853 | sumInfeasibility += value; |
---|
| 1854 | } |
---|
| 1855 | value = rowActivity[i] - rowUpper[i]; |
---|
| 1856 | if (value > primalTolerance) { |
---|
| 1857 | numberBadRows++; |
---|
| 1858 | largestInfeasibility = CoinMax(largestInfeasibility, value); |
---|
| 1859 | sumInfeasibility += value; |
---|
| 1860 | } |
---|
| 1861 | } |
---|
| 1862 | printf("%d bad rows, largest inf %g sum %g\n", |
---|
| 1863 | numberBadRows, largestInfeasibility, sumInfeasibility); |
---|
| 1864 | delete [] rowActivity; |
---|
| 1865 | } |
---|
[1271] | 1866 | #endif |
---|
| 1867 | #ifdef COIN_HAS_CLP |
---|
[1286] | 1868 | OsiClpSolverInterface * clpSolver |
---|
| 1869 | = dynamic_cast<OsiClpSolverInterface *> (newSolver); |
---|
| 1870 | if (clpSolver) { |
---|
| 1871 | ClpSimplex * simplex = clpSolver->getModelPtr(); |
---|
| 1872 | //simplex->writeBasis("test.bas",true,2); |
---|
| 1873 | //simplex->writeMps("test.mps",2,1); |
---|
| 1874 | if (nFixed*3 > numberColumns*2) |
---|
| 1875 | simplex->allSlackBasis(); // may as well go from all slack |
---|
| 1876 | simplex->primal(1); |
---|
| 1877 | clpSolver->setWarmStart(NULL); |
---|
| 1878 | } |
---|
[1271] | 1879 | #endif |
---|
[1286] | 1880 | newSolver->initialSolve(); |
---|
| 1881 | if (newSolver->isProvenOptimal()) { |
---|
| 1882 | double value = newSolver->getObjValue() * newSolver->getObjSense(); |
---|
| 1883 | if (value < newSolutionValue) { |
---|
| 1884 | //newSolver->writeMps("query","mps"); |
---|
[1393] | 1885 | #ifdef JJF_ZERO |
---|
[1286] | 1886 | { |
---|
| 1887 | double saveOffset; |
---|
| 1888 | newSolver->getDblParam(OsiObjOffset, saveOffset); |
---|
| 1889 | const double * obj = newSolver->getObjCoefficients(); |
---|
| 1890 | double newTrueSolutionValue = -saveOffset; |
---|
| 1891 | double newSumInfeas = 0.0; |
---|
| 1892 | int numberColumns = newSolver->getNumCols(); |
---|
| 1893 | const double * solution = newSolver->getColSolution(); |
---|
| 1894 | for (int i = 0 ; i < numberColumns ; i++ ) { |
---|
| 1895 | if (newSolver->isInteger(i)) { |
---|
| 1896 | double value = solution[i]; |
---|
| 1897 | double nearest = floor(value + 0.5); |
---|
| 1898 | newSumInfeas += fabs(value - nearest); |
---|
| 1899 | } |
---|
| 1900 | if (solution[i]) |
---|
| 1901 | printf("%d obj %g val %g - total %g\n", i, obj[i], solution[i], |
---|
| 1902 | newTrueSolutionValue); |
---|
| 1903 | newTrueSolutionValue += obj[i] * solution[i]; |
---|
| 1904 | } |
---|
| 1905 | printf("obj %g - inf %g\n", newTrueSolutionValue, |
---|
| 1906 | newSumInfeas); |
---|
| 1907 | } |
---|
[1053] | 1908 | #endif |
---|
[1286] | 1909 | sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value); |
---|
| 1910 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1911 | << pumpPrint |
---|
| 1912 | << CoinMessageEol; |
---|
| 1913 | newSolutionValue = value; |
---|
| 1914 | memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double)); |
---|
| 1915 | } |
---|
| 1916 | } else { |
---|
| 1917 | //newSolver->writeMps("bad3.mps"); |
---|
| 1918 | exitAll = true; |
---|
| 1919 | break; |
---|
| 1920 | } |
---|
| 1921 | } |
---|
| 1922 | } else { |
---|
| 1923 | sprintf(pumpPrint, "Mini branch and bound did not improve solution (%.2f seconds)", |
---|
| 1924 | model_->getCurrentSeconds()); |
---|
| 1925 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1926 | << pumpPrint |
---|
| 1927 | << CoinMessageEol; |
---|
| 1928 | if (returnCode && newSolutionValue < saveValue + 1.0e-3 && nFixC + nFixC2) { |
---|
| 1929 | // may be able to do better |
---|
| 1930 | const double * lower = model_->solver()->getColLower(); |
---|
| 1931 | const double * upper = model_->solver()->getColUpper(); |
---|
| 1932 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 1933 | if (newSolver->isInteger(iColumn)) { |
---|
| 1934 | double value = floor(newSolution[iColumn] + 0.5); |
---|
| 1935 | newSolver->setColLower(iColumn, value); |
---|
| 1936 | newSolver->setColUpper(iColumn, value); |
---|
| 1937 | } else { |
---|
| 1938 | newSolver->setColLower(iColumn, lower[iColumn]); |
---|
| 1939 | newSolver->setColUpper(iColumn, upper[iColumn]); |
---|
| 1940 | } |
---|
| 1941 | } |
---|
| 1942 | newSolver->initialSolve(); |
---|
| 1943 | if (newSolver->isProvenOptimal()) { |
---|
| 1944 | double value = newSolver->getObjValue() * newSolver->getObjSense(); |
---|
| 1945 | if (value < saveValue) { |
---|
| 1946 | sprintf(pumpPrint, "Freeing continuous variables gives a solution of %g", value); |
---|
| 1947 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1948 | << pumpPrint |
---|
| 1949 | << CoinMessageEol; |
---|
| 1950 | newSolutionValue = value; |
---|
| 1951 | memcpy(betterSolution, newSolver->getColSolution(), numberColumns*sizeof(double)); |
---|
| 1952 | gotSolution = true; |
---|
| 1953 | } |
---|
| 1954 | } |
---|
| 1955 | } |
---|
| 1956 | } |
---|
| 1957 | if (gotSolution) { |
---|
| 1958 | if ((accumulate_&1) != 0) { |
---|
| 1959 | model_->incrementUsed(betterSolution); // for local search |
---|
| 1960 | numberSolutions++; |
---|
| 1961 | } |
---|
| 1962 | solutionValue = newSolutionValue; |
---|
| 1963 | solutionFound = true; |
---|
| 1964 | if (exitNow(newSolutionValue)) |
---|
| 1965 | exitAll = true; |
---|
| 1966 | CoinWarmStartBasis * basis = |
---|
| 1967 | dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ; |
---|
| 1968 | if (basis) { |
---|
| 1969 | bestBasis = * basis; |
---|
| 1970 | delete basis; |
---|
| 1971 | int action = model_->dealWithEventHandler(CbcEventHandler::heuristicSolution, newSolutionValue, betterSolution); |
---|
| 1972 | if (action == 0) { |
---|
| 1973 | double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(), numberColumns); |
---|
| 1974 | double saveObjectiveValue = model_->getMinimizationObjValue(); |
---|
| 1975 | model_->setBestSolution(betterSolution, numberColumns, newSolutionValue); |
---|
| 1976 | if (saveOldSolution && saveObjectiveValue < model_->getMinimizationObjValue()) |
---|
| 1977 | model_->setBestSolution(saveOldSolution, numberColumns, saveObjectiveValue); |
---|
| 1978 | delete [] saveOldSolution; |
---|
| 1979 | } |
---|
| 1980 | if (!action || model_->getCurrentSeconds() > model_->getMaximumSeconds()) { |
---|
| 1981 | exitAll = true; // exit |
---|
| 1982 | break; |
---|
| 1983 | } |
---|
| 1984 | } |
---|
| 1985 | } |
---|
| 1986 | } // end stopBAB while |
---|
| 1987 | delete newSolver; |
---|
| 1988 | } |
---|
| 1989 | if (solutionFound) finalReturnCode = 1; |
---|
| 1990 | cutoff = CoinMin(cutoff, solutionValue - model_->getCutoffIncrement()); |
---|
| 1991 | if (numberTries >= maximumRetries_ || !solutionFound || exitAll || cutoff < continuousObjectiveValue + 1.0e-7) { |
---|
| 1992 | break; |
---|
| 1993 | } else { |
---|
| 1994 | solutionFound = false; |
---|
| 1995 | if (absoluteIncrement_ > 0.0 || relativeIncrement_ > 0.0) { |
---|
| 1996 | double gap = relativeIncrement_ * fabs(solutionValue); |
---|
| 1997 | double change = CoinMax(gap, absoluteIncrement_); |
---|
| 1998 | cutoff = CoinMin(cutoff, solutionValue - change); |
---|
| 1999 | } else { |
---|
| 2000 | //double weights[10]={0.1,0.1,0.2,0.2,0.2,0.3,0.3,0.3,0.4,0.5}; |
---|
| 2001 | double weights[10] = {0.1, 0.2, 0.3, 0.3, 0.4, 0.4, 0.4, 0.5, 0.5, 0.6}; |
---|
| 2002 | cutoff -= weights[CoinMin(numberTries-1, 9)] * (cutoff - continuousObjectiveValue); |
---|
| 2003 | } |
---|
| 2004 | // But round down |
---|
| 2005 | if (exactMultiple) |
---|
| 2006 | cutoff = exactMultiple * floor(cutoff / exactMultiple); |
---|
| 2007 | if (cutoff < continuousObjectiveValue) |
---|
| 2008 | break; |
---|
| 2009 | sprintf(pumpPrint, "Round again with cutoff of %g", cutoff); |
---|
| 2010 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 2011 | << pumpPrint |
---|
| 2012 | << CoinMessageEol; |
---|
| 2013 | if ((accumulate_&3) < 2 && usedColumn) { |
---|
| 2014 | for (int i = 0; i < numberColumns; i++) |
---|
| 2015 | usedColumn[i] = -1; |
---|
| 2016 | } |
---|
| 2017 | averageIterationsPerTry = totalNumberIterations / numberTries; |
---|
| 2018 | numberIterationsLastPass = totalNumberIterations; |
---|
| 2019 | totalNumberPasses += numberPasses - 1; |
---|
| 2020 | } |
---|
[175] | 2021 | } |
---|
[1364] | 2022 | /* |
---|
| 2023 | End of the `exitAll' loop. |
---|
| 2024 | */ |
---|
[1271] | 2025 | #ifdef RAND_RAND |
---|
[1286] | 2026 | delete [] randomFactor; |
---|
[1271] | 2027 | #endif |
---|
[1286] | 2028 | delete solver; // probably NULL but do anyway |
---|
| 2029 | if (!finalReturnCode && closestSolution && closestObjectiveValue <= 10.0 && |
---|
| 2030 | usedColumn && !model_->maximumSecondsReached()) { |
---|
| 2031 | // try a bit of branch and bound |
---|
| 2032 | OsiSolverInterface * newSolver = cloneBut(1); // was model_->continuousSolver()->clone(); |
---|
| 2033 | const double * colLower = newSolver->getColLower(); |
---|
| 2034 | const double * colUpper = newSolver->getColUpper(); |
---|
| 2035 | int i; |
---|
| 2036 | double rhs = 0.0; |
---|
| 2037 | for (i = 0; i < numberIntegersOrig; i++) { |
---|
| 2038 | int iColumn = integerVariableOrig[i]; |
---|
| 2039 | int direction = closestSolution[i]; |
---|
| 2040 | closestSolution[i] = iColumn; |
---|
| 2041 | if (direction == 0) { |
---|
| 2042 | // keep close to LB |
---|
| 2043 | rhs += colLower[iColumn]; |
---|
| 2044 | lastSolution[i] = 1.0; |
---|
| 2045 | } else { |
---|
| 2046 | // keep close to UB |
---|
| 2047 | rhs -= colUpper[iColumn]; |
---|
| 2048 | lastSolution[i] = -1.0; |
---|
| 2049 | } |
---|
| 2050 | } |
---|
| 2051 | newSolver->addRow(numberIntegersOrig, closestSolution, |
---|
| 2052 | lastSolution, -COIN_DBL_MAX, rhs + 10.0); |
---|
| 2053 | //double saveValue = newSolutionValue; |
---|
| 2054 | //newSolver->writeMps("sub"); |
---|
| 2055 | int returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, newSolutionValue, |
---|
| 2056 | newSolutionValue, "CbcHeuristicLocalAfterFPump"); |
---|
| 2057 | if (returnCode < 0) |
---|
| 2058 | returnCode = 0; // returned on size |
---|
| 2059 | if ((returnCode&2) != 0) { |
---|
| 2060 | // could add cut |
---|
| 2061 | returnCode &= ~2; |
---|
| 2062 | } |
---|
| 2063 | if (returnCode) { |
---|
| 2064 | //printf("old sol of %g new of %g\n",saveValue,newSolutionValue); |
---|
| 2065 | memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); |
---|
| 2066 | //abort(); |
---|
| 2067 | solutionValue = newSolutionValue; |
---|
| 2068 | solutionFound = true; |
---|
| 2069 | if (exitNow(newSolutionValue)) |
---|
| 2070 | exitAll = true; |
---|
| 2071 | } |
---|
| 2072 | delete newSolver; |
---|
[640] | 2073 | } |
---|
[1286] | 2074 | delete clonedSolver; |
---|
| 2075 | delete [] roundingSolution; |
---|
| 2076 | delete [] usedColumn; |
---|
| 2077 | delete [] lastSolution; |
---|
| 2078 | delete [] newSolution; |
---|
| 2079 | delete [] closestSolution; |
---|
| 2080 | delete [] integerVariable; |
---|
| 2081 | delete [] firstPerturbedObjective; |
---|
| 2082 | delete [] firstPerturbedSolution; |
---|
| 2083 | if (solutionValue == incomingObjective) |
---|
| 2084 | sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting - took %.2f seconds", |
---|
| 2085 | model_->getCurrentSeconds(), CoinCpuTime() - time1); |
---|
| 2086 | else |
---|
| 2087 | sprintf(pumpPrint, "After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds", |
---|
| 2088 | model_->getCurrentSeconds(), solutionValue, CoinCpuTime() - time1); |
---|
| 2089 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 2090 | << pumpPrint |
---|
| 2091 | << CoinMessageEol; |
---|
| 2092 | if (bestBasis.getNumStructural()) |
---|
| 2093 | model_->setBestSolutionBasis(bestBasis); |
---|
| 2094 | //model_->setMinimizationObjValue(saveBestObjective); |
---|
| 2095 | if ((accumulate_&1) != 0 && numberSolutions > 1 && !model_->getSolutionCount()) { |
---|
| 2096 | model_->setSolutionCount(1); // for local search |
---|
| 2097 | model_->setNumberHeuristicSolutions(1); |
---|
[640] | 2098 | } |
---|
[1286] | 2099 | #ifdef COIN_DEVELOP |
---|
| 2100 | { |
---|
| 2101 | double ncol = model_->solver()->getNumCols(); |
---|
| 2102 | double nrow = model_->solver()->getNumRows(); |
---|
| 2103 | printf("XXX total iterations %d ratios - %g %g %g\n", |
---|
| 2104 | totalNumberIterations, |
---|
| 2105 | static_cast<double> (totalNumberIterations) / nrow, |
---|
| 2106 | static_cast<double> (totalNumberIterations) / ncol, |
---|
| 2107 | static_cast<double> (totalNumberIterations) / (2*nrow + 2*ncol)); |
---|
[640] | 2108 | } |
---|
[931] | 2109 | #endif |
---|
[1286] | 2110 | return finalReturnCode; |
---|
[175] | 2111 | } |
---|
| 2112 | |
---|
| 2113 | /**************************END MAIN PROCEDURE ***********************************/ |
---|
| 2114 | |
---|
| 2115 | // update model |
---|
| 2116 | void CbcHeuristicFPump::setModel(CbcModel * model) |
---|
| 2117 | { |
---|
[1286] | 2118 | model_ = model; |
---|
[175] | 2119 | } |
---|
| 2120 | |
---|
| 2121 | /* Rounds solution - down if < downValue |
---|
| 2122 | returns 1 if current is a feasible solution |
---|
| 2123 | */ |
---|
[1286] | 2124 | int |
---|
| 2125 | CbcHeuristicFPump::rounds(OsiSolverInterface * solver, double * solution, |
---|
| 2126 | //const double * objective, |
---|
| 2127 | int numberIntegers, const int * integerVariable, |
---|
| 2128 | /*char * pumpPrint,*/ int iter, |
---|
| 2129 | /*bool roundExpensive,*/ double downValue, int *flip) |
---|
[175] | 2130 | { |
---|
[1286] | 2131 | double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance); |
---|
| 2132 | double primalTolerance ; |
---|
| 2133 | solver->getDblParam(OsiPrimalTolerance, primalTolerance) ; |
---|
[175] | 2134 | |
---|
[1286] | 2135 | int i; |
---|
[175] | 2136 | |
---|
[1286] | 2137 | const double * cost = solver->getObjCoefficients(); |
---|
| 2138 | int flip_up = 0; |
---|
| 2139 | int flip_down = 0; |
---|
| 2140 | double v = randomNumberGenerator_.randomDouble() * 20.0; |
---|
| 2141 | int nn = 10 + static_cast<int> (v); |
---|
| 2142 | int nnv = 0; |
---|
| 2143 | int * list = new int [nn]; |
---|
| 2144 | double * val = new double [nn]; |
---|
| 2145 | for (i = 0; i < nn; i++) val[i] = .001; |
---|
[175] | 2146 | |
---|
[1286] | 2147 | const double * rowLower = solver->getRowLower(); |
---|
| 2148 | const double * rowUpper = solver->getRowUpper(); |
---|
| 2149 | int numberRows = solver->getNumRows(); |
---|
| 2150 | if (false && (iter&1) != 0) { |
---|
| 2151 | // Do set covering variables |
---|
| 2152 | const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow(); |
---|
| 2153 | const double * elementByRow = matrixByRow->getElements(); |
---|
| 2154 | const int * column = matrixByRow->getIndices(); |
---|
| 2155 | const CoinBigIndex * rowStart = matrixByRow->getVectorStarts(); |
---|
| 2156 | const int * rowLength = matrixByRow->getVectorLengths(); |
---|
| 2157 | for (i = 0; i < numberRows; i++) { |
---|
| 2158 | if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) { |
---|
| 2159 | bool cover = true; |
---|
| 2160 | double largest = 0.0; |
---|
| 2161 | int jColumn = -1; |
---|
| 2162 | for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) { |
---|
| 2163 | int iColumn = column[k]; |
---|
| 2164 | if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) { |
---|
| 2165 | cover = false; |
---|
| 2166 | break; |
---|
| 2167 | } else { |
---|
| 2168 | if (solution[iColumn]) { |
---|
| 2169 | double value = solution[iColumn] * |
---|
| 2170 | (randomNumberGenerator_.randomDouble() + 5.0); |
---|
| 2171 | if (value > largest) { |
---|
| 2172 | largest = value; |
---|
| 2173 | jColumn = iColumn; |
---|
| 2174 | } |
---|
| 2175 | } |
---|
| 2176 | } |
---|
| 2177 | } |
---|
| 2178 | if (cover) { |
---|
| 2179 | for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) { |
---|
| 2180 | int iColumn = column[k]; |
---|
| 2181 | if (iColumn == jColumn) |
---|
| 2182 | solution[iColumn] = 1.0; |
---|
| 2183 | else |
---|
| 2184 | solution[iColumn] = 0.0; |
---|
| 2185 | } |
---|
| 2186 | } |
---|
| 2187 | } |
---|
| 2188 | } |
---|
| 2189 | } |
---|
| 2190 | int numberColumns = solver->getNumCols(); |
---|
[1393] | 2191 | #ifdef JJF_ZERO |
---|
[1053] | 2192 | // Do set covering variables |
---|
| 2193 | const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow(); |
---|
| 2194 | const double * elementByRow = matrixByRow->getElements(); |
---|
| 2195 | const int * column = matrixByRow->getIndices(); |
---|
| 2196 | const CoinBigIndex * rowStart = matrixByRow->getVectorStarts(); |
---|
| 2197 | const int * rowLength = matrixByRow->getVectorLengths(); |
---|
[1286] | 2198 | double * sortTemp = new double[numberColumns]; |
---|
| 2199 | int * whichTemp = new int [numberColumns]; |
---|
| 2200 | char * rowTemp = new char [numberRows]; |
---|
| 2201 | memset(rowTemp, 0, numberRows); |
---|
| 2202 | for (i = 0; i < numberColumns; i++) |
---|
| 2203 | whichTemp[i] = -1; |
---|
| 2204 | int nSOS = 0; |
---|
| 2205 | for (i = 0; i < numberRows; i++) { |
---|
| 2206 | if (rowLower[i] == 1.0 && rowUpper[i] == 1.0) { |
---|
| 2207 | bool cover = true; |
---|
| 2208 | for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) { |
---|
| 2209 | int iColumn = column[k]; |
---|
| 2210 | if (elementByRow[k] != 1.0 || !solver->isInteger(iColumn)) { |
---|
| 2211 | cover = false; |
---|
| 2212 | break; |
---|
| 2213 | } |
---|
| 2214 | } |
---|
| 2215 | if (cover) { |
---|
| 2216 | rowTemp[i] = 1; |
---|
| 2217 | nSOS++; |
---|
| 2218 | for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) { |
---|
| 2219 | int iColumn = column[k]; |
---|
| 2220 | double value = solution[iColumn]; |
---|
| 2221 | whichTemp[iColumn] = iColumn; |
---|
| 2222 | } |
---|
| 2223 | } |
---|
| 2224 | } |
---|
[1013] | 2225 | } |
---|
[1286] | 2226 | if (nSOS) { |
---|
| 2227 | // Column copy |
---|
| 2228 | const CoinPackedMatrix * matrixByColumn = solver->getMatrixByCol(); |
---|
| 2229 | //const double * element = matrixByColumn->getElements(); |
---|
| 2230 | const int * row = matrixByColumn->getIndices(); |
---|
| 2231 | const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts(); |
---|
| 2232 | const int * columnLength = matrixByColumn->getVectorLengths(); |
---|
| 2233 | int nLook = 0; |
---|
| 2234 | for (i = 0; i < numberColumns; i++) { |
---|
| 2235 | if (whichTemp[i] >= 0) { |
---|
| 2236 | whichTemp[nLook] = i; |
---|
| 2237 | double value = solution[i]; |
---|
| 2238 | if (value < 0.5) |
---|
| 2239 | value *= (0.1 * randomNumberGenerator_.randomDouble() + 0.3); |
---|
| 2240 | sortTemp[nLook++] = -value; |
---|
| 2241 | } |
---|
| 2242 | } |
---|
| 2243 | CoinSort_2(sortTemp, sortTemp + nLook, whichTemp); |
---|
| 2244 | double smallest = 1.0; |
---|
| 2245 | int nFix = 0; |
---|
| 2246 | int nOne = 0; |
---|
| 2247 | for (int j = 0; j < nLook; j++) { |
---|
| 2248 | int jColumn = whichTemp[j]; |
---|
| 2249 | double thisValue = solution[jColumn]; |
---|
| 2250 | if (!thisValue) |
---|
| 2251 | continue; |
---|
| 2252 | if (thisValue == 1.0) |
---|
| 2253 | nOne++; |
---|
| 2254 | smallest = CoinMin(smallest, thisValue); |
---|
| 2255 | solution[jColumn] = 1.0; |
---|
| 2256 | double largest = 0.0; |
---|
| 2257 | for (CoinBigIndex jEl = columnStart[jColumn]; |
---|
| 2258 | jEl < columnStart[jColumn] + columnLength[jColumn]; jEl++) { |
---|
| 2259 | int jRow = row[jEl]; |
---|
| 2260 | if (rowTemp[jRow]) { |
---|
| 2261 | for (CoinBigIndex k = rowStart[jRow]; k < rowStart[jRow] + rowLength[jRow]; k++) { |
---|
| 2262 | int iColumn = column[k]; |
---|
| 2263 | if (solution[iColumn]) { |
---|
| 2264 | if (iColumn != jColumn) { |
---|
| 2265 | double value = solution[iColumn]; |
---|
| 2266 | if (value > largest) |
---|
| 2267 | largest = value; |
---|
| 2268 | solution[iColumn] = 0.0; |
---|
| 2269 | } |
---|
| 2270 | } |
---|
| 2271 | } |
---|
| 2272 | } |
---|
| 2273 | } |
---|
| 2274 | if (largest > thisValue) |
---|
| 2275 | printf("%d was at %g - chosen over a value of %g\n", |
---|
| 2276 | jColumn, thisValue, largest); |
---|
| 2277 | nFix++; |
---|
| 2278 | } |
---|
| 2279 | printf("%d fixed out of %d (%d at one already)\n", |
---|
| 2280 | nFix, nLook, nOne); |
---|
[1013] | 2281 | } |
---|
[1286] | 2282 | delete [] sortTemp; |
---|
| 2283 | delete [] whichTemp; |
---|
| 2284 | delete [] rowTemp; |
---|
[1013] | 2285 | #endif |
---|
[1286] | 2286 | const double * columnLower = solver->getColLower(); |
---|
| 2287 | const double * columnUpper = solver->getColUpper(); |
---|
| 2288 | // Check if valid with current solution (allow for 0.99999999s) |
---|
| 2289 | for (i = 0; i < numberIntegers; i++) { |
---|
| 2290 | int iColumn = integerVariable[i]; |
---|
| 2291 | double value = solution[iColumn]; |
---|
| 2292 | double round = floor(value + 0.5); |
---|
| 2293 | if (fabs(value - round) > primalTolerance) |
---|
| 2294 | break; |
---|
[1053] | 2295 | } |
---|
[1286] | 2296 | if (i == numberIntegers) { |
---|
| 2297 | // may be able to use solution even if 0.99999's |
---|
| 2298 | double * saveLower = CoinCopyOfArray(columnLower, numberColumns); |
---|
| 2299 | double * saveUpper = CoinCopyOfArray(columnUpper, numberColumns); |
---|
| 2300 | double * saveSolution = CoinCopyOfArray(solution, numberColumns); |
---|
| 2301 | double * tempSolution = CoinCopyOfArray(solution, numberColumns); |
---|
| 2302 | CoinWarmStartBasis * saveBasis = |
---|
| 2303 | dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; |
---|
| 2304 | for (i = 0; i < numberIntegers; i++) { |
---|
| 2305 | int iColumn = integerVariable[i]; |
---|
| 2306 | double value = solution[iColumn]; |
---|
| 2307 | double round = floor(value + 0.5); |
---|
| 2308 | solver->setColLower(iColumn, round); |
---|
| 2309 | solver->setColUpper(iColumn, round); |
---|
| 2310 | tempSolution[iColumn] = round; |
---|
| 2311 | } |
---|
| 2312 | solver->setColSolution(tempSolution); |
---|
| 2313 | delete [] tempSolution; |
---|
| 2314 | solver->resolve(); |
---|
| 2315 | solver->setColLower(saveLower); |
---|
| 2316 | solver->setColUpper(saveUpper); |
---|
| 2317 | solver->setWarmStart(saveBasis); |
---|
| 2318 | delete [] saveLower; |
---|
| 2319 | delete [] saveUpper; |
---|
| 2320 | delete saveBasis; |
---|
| 2321 | if (!solver->isProvenOptimal()) { |
---|
| 2322 | solver->setColSolution(saveSolution); |
---|
| 2323 | } |
---|
| 2324 | delete [] saveSolution; |
---|
| 2325 | if (solver->isProvenOptimal()) { |
---|
| 2326 | // feasible |
---|
| 2327 | delete [] list; |
---|
| 2328 | delete [] val; |
---|
| 2329 | return 1; |
---|
| 2330 | } |
---|
[1127] | 2331 | } |
---|
[1286] | 2332 | //double * saveSolution = CoinCopyOfArray(solution,numberColumns); |
---|
| 2333 | // return rounded solution |
---|
| 2334 | for (i = 0; i < numberIntegers; i++) { |
---|
| 2335 | int iColumn = integerVariable[i]; |
---|
| 2336 | double value = solution[iColumn]; |
---|
| 2337 | double round = floor(value + primalTolerance); |
---|
| 2338 | if (value - round > downValue) round += 1.; |
---|
[1393] | 2339 | #ifndef JJF_ONE |
---|
[1286] | 2340 | if (round < integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++; |
---|
| 2341 | if (round > 1. - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++; |
---|
[1013] | 2342 | #else |
---|
[1286] | 2343 | if (round < columnLower[iColumn] + integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++; |
---|
| 2344 | if (round > columnUpper[iColumn] - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++; |
---|
[1013] | 2345 | #endif |
---|
[1286] | 2346 | if (flip_up + flip_down == 0) { |
---|
| 2347 | for (int k = 0; k < nn; k++) { |
---|
| 2348 | if (fabs(value - round) > val[k]) { |
---|
| 2349 | nnv++; |
---|
| 2350 | for (int j = nn - 2; j >= k; j--) { |
---|
| 2351 | val[j+1] = val[j]; |
---|
| 2352 | list[j+1] = list[j]; |
---|
| 2353 | } |
---|
| 2354 | val[k] = fabs(value - round); |
---|
| 2355 | list[k] = iColumn; |
---|
| 2356 | break; |
---|
| 2357 | } |
---|
| 2358 | } |
---|
| 2359 | } |
---|
| 2360 | solution[iColumn] = round; |
---|
[175] | 2361 | } |
---|
| 2362 | |
---|
[1286] | 2363 | if (nnv > nn) nnv = nn; |
---|
| 2364 | //if (iter != 0) |
---|
| 2365 | //sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down); |
---|
| 2366 | *flip = flip_up + flip_down; |
---|
[175] | 2367 | |
---|
[1286] | 2368 | if (*flip == 0 && iter != 0) { |
---|
| 2369 | //sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn); |
---|
| 2370 | for (i = 0; i < nnv; i++) { |
---|
| 2371 | // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6 |
---|
| 2372 | int index = list[i]; |
---|
| 2373 | double value = solution[index]; |
---|
| 2374 | if (value <= 1.0) { |
---|
| 2375 | solution[index] = 1.0 - value; |
---|
| 2376 | } else if (value < columnLower[index] + integerTolerance) { |
---|
| 2377 | solution[index] = value + 1.0; |
---|
| 2378 | } else if (value > columnUpper[index] - integerTolerance) { |
---|
| 2379 | solution[index] = value - 1.0; |
---|
| 2380 | } else { |
---|
| 2381 | solution[index] = value - 1.0; |
---|
| 2382 | } |
---|
| 2383 | } |
---|
| 2384 | *flip = nnv; |
---|
| 2385 | } else { |
---|
| 2386 | //sprintf(pumpPrint+strlen(pumpPrint)," "); |
---|
[1100] | 2387 | } |
---|
[1286] | 2388 | delete [] list; |
---|
| 2389 | delete [] val; |
---|
| 2390 | //iter++; |
---|
| 2391 | |
---|
| 2392 | // get row activities |
---|
| 2393 | double * rowActivity = new double[numberRows]; |
---|
| 2394 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
| 2395 | solver->getMatrixByCol()->times(solution, rowActivity) ; |
---|
| 2396 | double largestInfeasibility = primalTolerance; |
---|
| 2397 | double sumInfeasibility = 0.0; |
---|
| 2398 | int numberBadRows = 0; |
---|
| 2399 | for (i = 0 ; i < numberRows ; i++) { |
---|
| 2400 | double value; |
---|
| 2401 | value = rowLower[i] - rowActivity[i]; |
---|
| 2402 | if (value > primalTolerance) { |
---|
| 2403 | numberBadRows++; |
---|
| 2404 | largestInfeasibility = CoinMax(largestInfeasibility, value); |
---|
| 2405 | sumInfeasibility += value; |
---|
| 2406 | } |
---|
| 2407 | value = rowActivity[i] - rowUpper[i]; |
---|
| 2408 | if (value > primalTolerance) { |
---|
| 2409 | numberBadRows++; |
---|
| 2410 | largestInfeasibility = CoinMax(largestInfeasibility, value); |
---|
| 2411 | sumInfeasibility += value; |
---|
| 2412 | } |
---|
[1100] | 2413 | } |
---|
[1393] | 2414 | #ifdef JJF_ZERO |
---|
[1286] | 2415 | if (largestInfeasibility > primalTolerance && numberBadRows*10 < numberRows) { |
---|
| 2416 | // Can we improve by flipping |
---|
| 2417 | for (int iPass = 0; iPass < 10; iPass++) { |
---|
| 2418 | int numberColumns = solver->getNumCols(); |
---|
| 2419 | const CoinPackedMatrix * matrixByCol = solver->getMatrixByCol(); |
---|
| 2420 | const double * element = matrixByCol->getElements(); |
---|
| 2421 | const int * row = matrixByCol->getIndices(); |
---|
| 2422 | const CoinBigIndex * columnStart = matrixByCol->getVectorStarts(); |
---|
| 2423 | const int * columnLength = matrixByCol->getVectorLengths(); |
---|
| 2424 | double oldSum = sumInfeasibility; |
---|
| 2425 | // First improve by moving continuous ones |
---|
| 2426 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 2427 | if (!solver->isInteger(iColumn)) { |
---|
| 2428 | double solValue = solution[iColumn]; |
---|
| 2429 | double thetaUp = columnUpper[iColumn] - solValue; |
---|
| 2430 | double improvementUp = 0.0; |
---|
| 2431 | if (thetaUp > primalTolerance) { |
---|
| 2432 | // can go up |
---|
| 2433 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
| 2434 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 2435 | int iRow = row[j]; |
---|
| 2436 | double distanceUp = rowUpper[iRow] - rowActivity[iRow]; |
---|
| 2437 | double distanceDown = rowLower[iRow] - rowActivity[iRow]; |
---|
| 2438 | double el = element[j]; |
---|
| 2439 | if (el > 0.0) { |
---|
| 2440 | // positive element |
---|
| 2441 | if (distanceUp > 0.0) { |
---|
| 2442 | if (thetaUp*el > distanceUp) |
---|
| 2443 | thetaUp = distanceUp / el; |
---|
| 2444 | } else { |
---|
| 2445 | improvementUp -= el; |
---|
| 2446 | } |
---|
| 2447 | if (distanceDown > 0.0) { |
---|
| 2448 | if (thetaUp*el > distanceDown) |
---|
| 2449 | thetaUp = distanceDown / el; |
---|
| 2450 | improvementUp += el; |
---|
| 2451 | } |
---|
| 2452 | } else { |
---|
| 2453 | // negative element |
---|
| 2454 | if (distanceDown < 0.0) { |
---|
| 2455 | if (thetaUp*el < distanceDown) |
---|
| 2456 | thetaUp = distanceDown / el; |
---|
| 2457 | } else { |
---|
| 2458 | improvementUp += el; |
---|
| 2459 | } |
---|
| 2460 | if (distanceUp < 0.0) { |
---|
| 2461 | if (thetaUp*el < distanceUp) |
---|
| 2462 | thetaUp = distanceUp / el; |
---|
| 2463 | improvementUp -= el; |
---|
| 2464 | } |
---|
| 2465 | } |
---|
| 2466 | } |
---|
| 2467 | } |
---|
| 2468 | double thetaDown = solValue - columnLower[iColumn]; |
---|
| 2469 | double improvementDown = 0.0; |
---|
| 2470 | if (thetaDown > primalTolerance) { |
---|
| 2471 | // can go down |
---|
| 2472 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
| 2473 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 2474 | int iRow = row[j]; |
---|
| 2475 | double distanceUp = rowUpper[iRow] - rowActivity[iRow]; |
---|
| 2476 | double distanceDown = rowLower[iRow] - rowActivity[iRow]; |
---|
| 2477 | double el = -element[j]; // not change in sign form up |
---|
| 2478 | if (el > 0.0) { |
---|
| 2479 | // positive element |
---|
| 2480 | if (distanceUp > 0.0) { |
---|
| 2481 | if (thetaDown*el > distanceUp) |
---|
| 2482 | thetaDown = distanceUp / el; |
---|
| 2483 | } else { |
---|
| 2484 | improvementDown -= el; |
---|
| 2485 | } |
---|
| 2486 | if (distanceDown > 0.0) { |
---|
| 2487 | if (thetaDown*el > distanceDown) |
---|
| 2488 | thetaDown = distanceDown / el; |
---|
| 2489 | improvementDown += el; |
---|
| 2490 | } |
---|
| 2491 | } else { |
---|
| 2492 | // negative element |
---|
| 2493 | if (distanceDown < 0.0) { |
---|
| 2494 | if (thetaDown*el < distanceDown) |
---|
| 2495 | thetaDown = distanceDown / el; |
---|
| 2496 | } else { |
---|
| 2497 | improvementDown += el; |
---|
| 2498 | } |
---|
| 2499 | if (distanceUp < 0.0) { |
---|
| 2500 | if (thetaDown*el < distanceUp) |
---|
| 2501 | thetaDown = distanceUp / el; |
---|
| 2502 | improvementDown -= el; |
---|
| 2503 | } |
---|
| 2504 | } |
---|
| 2505 | } |
---|
| 2506 | if (thetaUp < 1.0e-8) |
---|
| 2507 | improvementUp = 0.0; |
---|
| 2508 | if (thetaDown < 1.0e-8) |
---|
| 2509 | improvementDown = 0.0; |
---|
| 2510 | double theta; |
---|
| 2511 | if (improvementUp >= improvementDown) { |
---|
| 2512 | theta = thetaUp; |
---|
| 2513 | } else { |
---|
| 2514 | improvementUp = improvementDown; |
---|
| 2515 | theta = -thetaDown; |
---|
| 2516 | } |
---|
| 2517 | if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) { |
---|
| 2518 | // Could move |
---|
| 2519 | double oldSum = 0.0; |
---|
| 2520 | double newSum = 0.0; |
---|
| 2521 | solution[iColumn] += theta; |
---|
| 2522 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
| 2523 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 2524 | int iRow = row[j]; |
---|
| 2525 | double lower = rowLower[iRow]; |
---|
| 2526 | double upper = rowUpper[iRow]; |
---|
| 2527 | double value = rowActivity[iRow]; |
---|
| 2528 | if (value > upper) |
---|
| 2529 | oldSum += value - upper; |
---|
| 2530 | else if (value < lower) |
---|
| 2531 | oldSum += lower - value; |
---|
| 2532 | value += theta * element[j]; |
---|
| 2533 | rowActivity[iRow] = value; |
---|
| 2534 | if (value > upper) |
---|
| 2535 | newSum += value - upper; |
---|
| 2536 | else if (value < lower) |
---|
| 2537 | newSum += lower - value; |
---|
| 2538 | } |
---|
| 2539 | assert (newSum <= oldSum); |
---|
| 2540 | sumInfeasibility += newSum - oldSum; |
---|
| 2541 | } |
---|
| 2542 | } |
---|
| 2543 | } |
---|
| 2544 | } |
---|
| 2545 | // Now flip some integers? |
---|
[1393] | 2546 | #ifdef JJF_ZERO |
---|
[1286] | 2547 | for (i = 0; i < numberIntegers; i++) { |
---|
| 2548 | int iColumn = integerVariable[i]; |
---|
| 2549 | double solValue = solution[iColumn]; |
---|
| 2550 | assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8); |
---|
| 2551 | double improvementUp = 0.0; |
---|
| 2552 | if (columnUpper[iColumn] >= solValue + 1.0) { |
---|
| 2553 | // can go up |
---|
| 2554 | double oldSum = 0.0; |
---|
| 2555 | double newSum = 0.0; |
---|
| 2556 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
| 2557 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 2558 | int iRow = row[j]; |
---|
| 2559 | double lower = rowLower[iRow]; |
---|
| 2560 | double upper = rowUpper[iRow]; |
---|
| 2561 | double value = rowActivity[iRow]; |
---|
| 2562 | if (value > upper) |
---|
| 2563 | oldSum += value - upper; |
---|
| 2564 | else if (value < lower) |
---|
| 2565 | oldSum += lower - value; |
---|
| 2566 | value += element[j]; |
---|
| 2567 | if (value > upper) |
---|
| 2568 | newSum += value - upper; |
---|
| 2569 | else if (value < lower) |
---|
| 2570 | newSum += lower - value; |
---|
| 2571 | } |
---|
| 2572 | improvementUp = oldSum - newSum; |
---|
| 2573 | } |
---|
| 2574 | double improvementDown = 0.0; |
---|
| 2575 | if (columnLower[iColumn] <= solValue - 1.0) { |
---|
| 2576 | // can go down |
---|
| 2577 | double oldSum = 0.0; |
---|
| 2578 | double newSum = 0.0; |
---|
| 2579 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
| 2580 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 2581 | int iRow = row[j]; |
---|
| 2582 | double lower = rowLower[iRow]; |
---|
| 2583 | double upper = rowUpper[iRow]; |
---|
| 2584 | double value = rowActivity[iRow]; |
---|
| 2585 | if (value > upper) |
---|
| 2586 | oldSum += value - upper; |
---|
| 2587 | else if (value < lower) |
---|
| 2588 | oldSum += lower - value; |
---|
| 2589 | value -= element[j]; |
---|
| 2590 | if (value > upper) |
---|
| 2591 | newSum += value - upper; |
---|
| 2592 | else if (value < lower) |
---|
| 2593 | newSum += lower - value; |
---|
| 2594 | } |
---|
| 2595 | improvementDown = oldSum - newSum; |
---|
| 2596 | } |
---|
| 2597 | double theta; |
---|
| 2598 | if (improvementUp >= improvementDown) { |
---|
| 2599 | theta = 1.0; |
---|
| 2600 | } else { |
---|
| 2601 | improvementUp = improvementDown; |
---|
| 2602 | theta = -1.0; |
---|
| 2603 | } |
---|
| 2604 | if (improvementUp > 1.0e-8 && fabs(theta) > 1.0e-8) { |
---|
| 2605 | // Could move |
---|
| 2606 | double oldSum = 0.0; |
---|
| 2607 | double newSum = 0.0; |
---|
| 2608 | solution[iColumn] += theta; |
---|
| 2609 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
| 2610 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 2611 | int iRow = row[j]; |
---|
| 2612 | double lower = rowLower[iRow]; |
---|
| 2613 | double upper = rowUpper[iRow]; |
---|
| 2614 | double value = rowActivity[iRow]; |
---|
| 2615 | if (value > upper) |
---|
| 2616 | oldSum += value - upper; |
---|
| 2617 | else if (value < lower) |
---|
| 2618 | oldSum += lower - value; |
---|
| 2619 | value += theta * element[j]; |
---|
| 2620 | rowActivity[iRow] = value; |
---|
| 2621 | if (value > upper) |
---|
| 2622 | newSum += value - upper; |
---|
| 2623 | else if (value < lower) |
---|
| 2624 | newSum += lower - value; |
---|
| 2625 | } |
---|
| 2626 | assert (newSum <= oldSum); |
---|
| 2627 | sumInfeasibility += newSum - oldSum; |
---|
| 2628 | } |
---|
| 2629 | } |
---|
[1100] | 2630 | #else |
---|
[1286] | 2631 | int bestColumn = -1; |
---|
| 2632 | double bestImprovement = primalTolerance; |
---|
| 2633 | double theta = 0.0; |
---|
| 2634 | for (i = 0; i < numberIntegers; i++) { |
---|
| 2635 | int iColumn = integerVariable[i]; |
---|
| 2636 | double solValue = solution[iColumn]; |
---|
| 2637 | assert (fabs(solValue - floor(solValue + 0.5)) < 1.0e-8); |
---|
| 2638 | double improvementUp = 0.0; |
---|
| 2639 | if (columnUpper[iColumn] >= solValue + 1.0) { |
---|
| 2640 | // can go up |
---|
| 2641 | double oldSum = 0.0; |
---|
| 2642 | double newSum = 0.0; |
---|
| 2643 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
| 2644 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 2645 | int iRow = row[j]; |
---|
| 2646 | double lower = rowLower[iRow]; |
---|
| 2647 | double upper = rowUpper[iRow]; |
---|
| 2648 | double value = rowActivity[iRow]; |
---|
| 2649 | if (value > upper) |
---|
| 2650 | oldSum += value - upper; |
---|
| 2651 | else if (value < lower) |
---|
| 2652 | oldSum += lower - value; |
---|
| 2653 | value += element[j]; |
---|
| 2654 | if (value > upper) |
---|
| 2655 | newSum += value - upper; |
---|
| 2656 | else if (value < lower) |
---|
| 2657 | newSum += lower - value; |
---|
| 2658 | } |
---|
| 2659 | improvementUp = oldSum - newSum; |
---|
| 2660 | } |
---|
| 2661 | double improvementDown = 0.0; |
---|
| 2662 | if (columnLower[iColumn] <= solValue - 1.0) { |
---|
| 2663 | // can go down |
---|
| 2664 | double oldSum = 0.0; |
---|
| 2665 | double newSum = 0.0; |
---|
| 2666 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
| 2667 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 2668 | int iRow = row[j]; |
---|
| 2669 | double lower = rowLower[iRow]; |
---|
| 2670 | double upper = rowUpper[iRow]; |
---|
| 2671 | double value = rowActivity[iRow]; |
---|
| 2672 | if (value > upper) |
---|
| 2673 | oldSum += value - upper; |
---|
| 2674 | else if (value < lower) |
---|
| 2675 | oldSum += lower - value; |
---|
| 2676 | value -= element[j]; |
---|
| 2677 | if (value > upper) |
---|
| 2678 | newSum += value - upper; |
---|
| 2679 | else if (value < lower) |
---|
| 2680 | newSum += lower - value; |
---|
| 2681 | } |
---|
| 2682 | improvementDown = oldSum - newSum; |
---|
| 2683 | } |
---|
| 2684 | double improvement = CoinMax(improvementUp, improvementDown); |
---|
| 2685 | if (improvement > bestImprovement) { |
---|
| 2686 | bestImprovement = improvement; |
---|
| 2687 | bestColumn = iColumn; |
---|
| 2688 | if (improvementUp > improvementDown) |
---|
| 2689 | theta = 1.0; |
---|
| 2690 | else |
---|
| 2691 | theta = -1.0; |
---|
| 2692 | } |
---|
| 2693 | } |
---|
| 2694 | if (bestColumn >= 0) { |
---|
| 2695 | // Could move |
---|
| 2696 | int iColumn = bestColumn; |
---|
| 2697 | double oldSum = 0.0; |
---|
| 2698 | double newSum = 0.0; |
---|
| 2699 | solution[iColumn] += theta; |
---|
| 2700 | for (CoinBigIndex j = columnStart[iColumn]; |
---|
| 2701 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 2702 | int iRow = row[j]; |
---|
| 2703 | double lower = rowLower[iRow]; |
---|
| 2704 | double upper = rowUpper[iRow]; |
---|
| 2705 | double value = rowActivity[iRow]; |
---|
| 2706 | if (value > upper) |
---|
| 2707 | oldSum += value - upper; |
---|
| 2708 | else if (value < lower) |
---|
| 2709 | oldSum += lower - value; |
---|
| 2710 | value += theta * element[j]; |
---|
| 2711 | rowActivity[iRow] = value; |
---|
| 2712 | if (value > upper) |
---|
| 2713 | newSum += value - upper; |
---|
| 2714 | else if (value < lower) |
---|
| 2715 | newSum += lower - value; |
---|
| 2716 | } |
---|
| 2717 | assert (newSum <= oldSum); |
---|
| 2718 | sumInfeasibility += newSum - oldSum; |
---|
| 2719 | } |
---|
[1100] | 2720 | #endif |
---|
[1286] | 2721 | if (oldSum <= sumInfeasibility + primalTolerance) |
---|
| 2722 | break; // no good |
---|
| 2723 | } |
---|
[1100] | 2724 | } |
---|
[1286] | 2725 | //delete [] saveSolution; |
---|
[1100] | 2726 | #endif |
---|
[1286] | 2727 | delete [] rowActivity; |
---|
| 2728 | return (largestInfeasibility > primalTolerance) ? 0 : 1; |
---|
[175] | 2729 | } |
---|
| 2730 | // Set maximum Time (default off) - also sets starttime to current |
---|
[1286] | 2731 | void |
---|
[175] | 2732 | CbcHeuristicFPump::setMaximumTime(double value) |
---|
| 2733 | { |
---|
[1286] | 2734 | startTime_ = CoinCpuTime(); |
---|
| 2735 | maximumTime_ = value; |
---|
[175] | 2736 | } |
---|
| 2737 | |
---|
[1088] | 2738 | # ifdef COIN_HAS_CLP |
---|
[1286] | 2739 | |
---|
[1088] | 2740 | //############################################################################# |
---|
| 2741 | // Constructors / Destructor / Assignment |
---|
| 2742 | //############################################################################# |
---|
| 2743 | |
---|
| 2744 | //------------------------------------------------------------------- |
---|
[1286] | 2745 | // Default Constructor |
---|
[1088] | 2746 | //------------------------------------------------------------------- |
---|
[1286] | 2747 | CbcDisasterHandler::CbcDisasterHandler (CbcModel * model) |
---|
| 2748 | : OsiClpDisasterHandler(), |
---|
| 2749 | cbcModel_(model) |
---|
[1088] | 2750 | { |
---|
[1286] | 2751 | if (model) { |
---|
| 2752 | osiModel_ |
---|
| 2753 | = dynamic_cast<OsiClpSolverInterface *> (model->solver()); |
---|
| 2754 | if (osiModel_) |
---|
| 2755 | setSimplex(osiModel_->getModelPtr()); |
---|
| 2756 | } |
---|
[1088] | 2757 | } |
---|
| 2758 | |
---|
| 2759 | //------------------------------------------------------------------- |
---|
[1286] | 2760 | // Copy constructor |
---|
[1088] | 2761 | //------------------------------------------------------------------- |
---|
[1286] | 2762 | CbcDisasterHandler::CbcDisasterHandler (const CbcDisasterHandler & rhs) |
---|
| 2763 | : OsiClpDisasterHandler(rhs), |
---|
| 2764 | cbcModel_(rhs.cbcModel_) |
---|
| 2765 | { |
---|
[1088] | 2766 | } |
---|
| 2767 | |
---|
| 2768 | |
---|
| 2769 | //------------------------------------------------------------------- |
---|
[1286] | 2770 | // Destructor |
---|
[1088] | 2771 | //------------------------------------------------------------------- |
---|
| 2772 | CbcDisasterHandler::~CbcDisasterHandler () |
---|
| 2773 | { |
---|
| 2774 | } |
---|
| 2775 | |
---|
| 2776 | //---------------------------------------------------------------- |
---|
[1286] | 2777 | // Assignment operator |
---|
[1088] | 2778 | //------------------------------------------------------------------- |
---|
| 2779 | CbcDisasterHandler & |
---|
[1286] | 2780 | CbcDisasterHandler::operator=(const CbcDisasterHandler & rhs) |
---|
[1088] | 2781 | { |
---|
[1286] | 2782 | if (this != &rhs) { |
---|
| 2783 | OsiClpDisasterHandler::operator=(rhs); |
---|
| 2784 | cbcModel_ = rhs.cbcModel_; |
---|
| 2785 | } |
---|
| 2786 | return *this; |
---|
[1088] | 2787 | } |
---|
| 2788 | //------------------------------------------------------------------- |
---|
| 2789 | // Clone |
---|
| 2790 | //------------------------------------------------------------------- |
---|
| 2791 | ClpDisasterHandler * CbcDisasterHandler::clone() const |
---|
| 2792 | { |
---|
[1286] | 2793 | return new CbcDisasterHandler(*this); |
---|
[1088] | 2794 | } |
---|
| 2795 | // Type of disaster 0 can fix, 1 abort |
---|
[1286] | 2796 | int |
---|
[1088] | 2797 | CbcDisasterHandler::typeOfDisaster() |
---|
| 2798 | { |
---|
[1286] | 2799 | if (!cbcModel_->parentModel() && |
---|
| 2800 | (cbcModel_->specialOptions()&2048) == 0) { |
---|
| 2801 | return 0; |
---|
| 2802 | } else { |
---|
| 2803 | if (cbcModel_->parentModel()) |
---|
| 2804 | cbcModel_->setMaximumNodes(0); |
---|
| 2805 | return 1; |
---|
| 2806 | } |
---|
[1088] | 2807 | } |
---|
| 2808 | /* set model. */ |
---|
[1286] | 2809 | void |
---|
[1088] | 2810 | CbcDisasterHandler::setCbcModel(CbcModel * model) |
---|
| 2811 | { |
---|
[1286] | 2812 | cbcModel_ = model; |
---|
| 2813 | if (model) { |
---|
| 2814 | osiModel_ |
---|
| 2815 | = dynamic_cast<OsiClpSolverInterface *> (model->solver()); |
---|
| 2816 | if (osiModel_) |
---|
| 2817 | setSimplex(osiModel_->getModelPtr()); |
---|
| 2818 | else |
---|
| 2819 | setSimplex(NULL); |
---|
| 2820 | } |
---|
[1088] | 2821 | } |
---|
| 2822 | #endif |
---|
[1286] | 2823 | |
---|