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