[1271] | 1 | /* $Id: CbcHeuristicLocal.cpp 1820 2012-11-29 13:28:09Z stefan $ */ |
---|
[175] | 2 | // Copyright (C) 2002, 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" |
---|
| 17 | #include "CbcMessage.hpp" |
---|
| 18 | #include "CbcHeuristicLocal.hpp" |
---|
[1802] | 19 | #include "CbcHeuristicFPump.hpp" |
---|
[175] | 20 | #include "CbcBranchActual.hpp" |
---|
| 21 | #include "CbcStrategy.hpp" |
---|
| 22 | #include "CglPreProcess.hpp" |
---|
| 23 | |
---|
| 24 | // Default Constructor |
---|
[1286] | 25 | CbcHeuristicLocal::CbcHeuristicLocal() |
---|
| 26 | : CbcHeuristic() |
---|
[175] | 27 | { |
---|
[1286] | 28 | numberSolutions_ = 0; |
---|
| 29 | swap_ = 0; |
---|
| 30 | used_ = NULL; |
---|
[1613] | 31 | lastRunDeep_ = -1000000; |
---|
[1802] | 32 | switches_ |= 16; // needs a new solution |
---|
[175] | 33 | } |
---|
| 34 | |
---|
| 35 | // Constructor with model - assumed before cuts |
---|
| 36 | |
---|
| 37 | CbcHeuristicLocal::CbcHeuristicLocal(CbcModel & model) |
---|
[1286] | 38 | : CbcHeuristic(model) |
---|
[175] | 39 | { |
---|
[1286] | 40 | numberSolutions_ = 0; |
---|
| 41 | swap_ = 0; |
---|
[1613] | 42 | lastRunDeep_ = -1000000; |
---|
[1802] | 43 | switches_ |= 16; // needs a new solution |
---|
[1286] | 44 | // Get a copy of original matrix |
---|
| 45 | assert(model.solver()); |
---|
| 46 | if (model.solver()->getNumRows()) { |
---|
| 47 | matrix_ = *model.solver()->getMatrixByCol(); |
---|
| 48 | } |
---|
| 49 | int numberColumns = model.solver()->getNumCols(); |
---|
| 50 | used_ = new int[numberColumns]; |
---|
| 51 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
[175] | 52 | } |
---|
| 53 | |
---|
[1286] | 54 | // Destructor |
---|
[175] | 55 | CbcHeuristicLocal::~CbcHeuristicLocal () |
---|
| 56 | { |
---|
[1286] | 57 | delete [] used_; |
---|
[175] | 58 | } |
---|
| 59 | |
---|
| 60 | // Clone |
---|
| 61 | CbcHeuristic * |
---|
| 62 | CbcHeuristicLocal::clone() const |
---|
| 63 | { |
---|
[1286] | 64 | return new CbcHeuristicLocal(*this); |
---|
[175] | 65 | } |
---|
[356] | 66 | // Create C++ lines to get to current state |
---|
[1286] | 67 | void |
---|
| 68 | CbcHeuristicLocal::generateCpp( FILE * fp) |
---|
[356] | 69 | { |
---|
[1286] | 70 | CbcHeuristicLocal other; |
---|
| 71 | fprintf(fp, "0#include \"CbcHeuristicLocal.hpp\"\n"); |
---|
| 72 | fprintf(fp, "3 CbcHeuristicLocal heuristicLocal(*cbcModel);\n"); |
---|
| 73 | CbcHeuristic::generateCpp(fp, "heuristicLocal"); |
---|
| 74 | if (swap_ != other.swap_) |
---|
| 75 | fprintf(fp, "3 heuristicLocal.setSearchType(%d);\n", swap_); |
---|
| 76 | else |
---|
| 77 | fprintf(fp, "4 heuristicLocal.setSearchType(%d);\n", swap_); |
---|
| 78 | fprintf(fp, "3 cbcModel->addHeuristic(&heuristicLocal);\n"); |
---|
[356] | 79 | } |
---|
[175] | 80 | |
---|
[1286] | 81 | // Copy constructor |
---|
[175] | 82 | CbcHeuristicLocal::CbcHeuristicLocal(const CbcHeuristicLocal & rhs) |
---|
[1286] | 83 | : |
---|
| 84 | CbcHeuristic(rhs), |
---|
| 85 | matrix_(rhs.matrix_), |
---|
| 86 | numberSolutions_(rhs.numberSolutions_), |
---|
| 87 | swap_(rhs.swap_) |
---|
[175] | 88 | { |
---|
[1286] | 89 | if (model_ && rhs.used_) { |
---|
| 90 | int numberColumns = model_->solver()->getNumCols(); |
---|
| 91 | used_ = CoinCopyOfArray(rhs.used_, numberColumns); |
---|
| 92 | } else { |
---|
| 93 | used_ = NULL; |
---|
| 94 | } |
---|
[175] | 95 | } |
---|
[640] | 96 | |
---|
[1286] | 97 | // Assignment operator |
---|
| 98 | CbcHeuristicLocal & |
---|
| 99 | CbcHeuristicLocal::operator=( const CbcHeuristicLocal & rhs) |
---|
[640] | 100 | { |
---|
[1286] | 101 | if (this != &rhs) { |
---|
| 102 | CbcHeuristic::operator=(rhs); |
---|
| 103 | matrix_ = rhs.matrix_; |
---|
| 104 | numberSolutions_ = rhs.numberSolutions_; |
---|
| 105 | swap_ = rhs.swap_; |
---|
| 106 | delete [] used_; |
---|
| 107 | if (model_ && rhs.used_) { |
---|
| 108 | int numberColumns = model_->solver()->getNumCols(); |
---|
| 109 | used_ = CoinCopyOfArray(rhs.used_, numberColumns); |
---|
| 110 | } else { |
---|
| 111 | used_ = NULL; |
---|
| 112 | } |
---|
[640] | 113 | } |
---|
[1286] | 114 | return *this; |
---|
[640] | 115 | } |
---|
| 116 | |
---|
[175] | 117 | // Resets stuff if model changes |
---|
[1286] | 118 | void |
---|
[1271] | 119 | CbcHeuristicLocal::resetModel(CbcModel * /*model*/) |
---|
[175] | 120 | { |
---|
[1286] | 121 | //CbcHeuristic::resetModel(model); |
---|
| 122 | delete [] used_; |
---|
| 123 | if (model_ && used_) { |
---|
| 124 | int numberColumns = model_->solver()->getNumCols(); |
---|
| 125 | used_ = new int[numberColumns]; |
---|
| 126 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
| 127 | } else { |
---|
| 128 | used_ = NULL; |
---|
| 129 | } |
---|
[175] | 130 | } |
---|
[1364] | 131 | /* |
---|
| 132 | Run a mini-BaB search after fixing all variables not marked as used by |
---|
| 133 | solution(). (See comments there for semantics.) |
---|
| 134 | |
---|
| 135 | Return values are: |
---|
| 136 | 1: smallBranchAndBound found a solution |
---|
| 137 | 0: everything else |
---|
| 138 | |
---|
| 139 | The degree of overload as return codes from smallBranchAndBound are folded |
---|
| 140 | into 0 is such that it's impossible to distinguish return codes that really |
---|
| 141 | require attention from a simple `nothing of interest'. |
---|
| 142 | */ |
---|
[175] | 143 | // This version fixes stuff and does IP |
---|
[1286] | 144 | int |
---|
[175] | 145 | CbcHeuristicLocal::solutionFix(double & objectiveValue, |
---|
[1286] | 146 | double * newSolution, |
---|
| 147 | const int * /*keep*/) |
---|
[175] | 148 | { |
---|
[1364] | 149 | /* |
---|
| 150 | If when is set to off (0), or set to root (1) and we're not at the root, |
---|
| 151 | return. If this heuristic discovered the current solution, don't continue. |
---|
| 152 | */ |
---|
| 153 | |
---|
[1286] | 154 | numCouldRun_++; |
---|
| 155 | // See if to do |
---|
| 156 | if (!when() || (when() == 1 && model_->phase() != 1)) |
---|
| 157 | return 0; // switched off |
---|
| 158 | // Don't do if it was this heuristic which found solution! |
---|
| 159 | if (this == model_->lastHeuristic()) |
---|
| 160 | return 0; |
---|
[1364] | 161 | /* |
---|
| 162 | Load up a new solver with the solution. |
---|
| 163 | |
---|
| 164 | Why continuousSolver(), as opposed to solver()? |
---|
| 165 | */ |
---|
[1286] | 166 | OsiSolverInterface * newSolver = model_->continuousSolver()->clone(); |
---|
| 167 | const double * colLower = newSolver->getColLower(); |
---|
| 168 | //const double * colUpper = newSolver->getColUpper(); |
---|
[175] | 169 | |
---|
[1286] | 170 | int numberIntegers = model_->numberIntegers(); |
---|
| 171 | const int * integerVariable = model_->integerVariable(); |
---|
[1364] | 172 | /* |
---|
| 173 | The net effect here is that anything that hasn't moved from its lower bound |
---|
| 174 | will be fixed at lower bound. |
---|
[1286] | 175 | |
---|
[1364] | 176 | See comments in solution() w.r.t. asymmetric treatment of upper and lower |
---|
| 177 | bounds. |
---|
| 178 | */ |
---|
| 179 | |
---|
[1286] | 180 | int i; |
---|
| 181 | int nFix = 0; |
---|
| 182 | for (i = 0; i < numberIntegers; i++) { |
---|
| 183 | int iColumn = integerVariable[i]; |
---|
| 184 | const OsiObject * object = model_->object(i); |
---|
| 185 | // get original bounds |
---|
| 186 | double originalLower; |
---|
| 187 | double originalUpper; |
---|
| 188 | getIntegerInformation( object, originalLower, originalUpper); |
---|
| 189 | newSolver->setColLower(iColumn, CoinMax(colLower[iColumn], originalLower)); |
---|
| 190 | if (!used_[iColumn]) { |
---|
| 191 | newSolver->setColUpper(iColumn, colLower[iColumn]); |
---|
| 192 | nFix++; |
---|
| 193 | } |
---|
[175] | 194 | } |
---|
[1364] | 195 | /* |
---|
| 196 | Try a `small' branch-and-bound search. The notion here is that we've fixed a |
---|
| 197 | lot of variables and reduced the amount of `free' problem to a point where a |
---|
| 198 | small BaB search will suffice to fully explore the remaining problem. This |
---|
| 199 | routine will execute integer presolve, then call branchAndBound to do the |
---|
| 200 | actual search. |
---|
| 201 | */ |
---|
[1286] | 202 | int returnCode = 0; |
---|
[1315] | 203 | #ifdef CLP_INVESTIGATE2 |
---|
| 204 | printf("Fixing %d out of %d (%d continuous)\n", |
---|
| 205 | nFix, numberIntegers, newSolver->getNumCols() - numberIntegers); |
---|
| 206 | #endif |
---|
| 207 | if (nFix*10 <= numberIntegers) { |
---|
| 208 | // see if we can fix more |
---|
| 209 | int * which = new int [2*(numberIntegers-nFix)]; |
---|
| 210 | int * sort = which + (numberIntegers - nFix); |
---|
| 211 | int n = 0; |
---|
| 212 | for (i = 0; i < numberIntegers; i++) { |
---|
| 213 | int iColumn = integerVariable[i]; |
---|
| 214 | if (used_[iColumn]) { |
---|
| 215 | which[n] = iColumn; |
---|
| 216 | sort[n++] = used_[iColumn]; |
---|
| 217 | } |
---|
| 218 | } |
---|
| 219 | CoinSort_2(sort, sort + n, which); |
---|
| 220 | // only half fixed in total |
---|
| 221 | n = CoinMin(n, numberIntegers / 2 - nFix); |
---|
| 222 | int allow = CoinMax(numberSolutions_ - 2, sort[0]); |
---|
| 223 | int nFix2 = 0; |
---|
| 224 | for (i = 0; i < n; i++) { |
---|
| 225 | int iColumn = integerVariable[i]; |
---|
| 226 | if (used_[iColumn] <= allow) { |
---|
| 227 | newSolver->setColUpper(iColumn, colLower[iColumn]); |
---|
| 228 | nFix2++; |
---|
| 229 | } else { |
---|
| 230 | break; |
---|
| 231 | } |
---|
| 232 | } |
---|
| 233 | delete [] which; |
---|
| 234 | nFix += nFix2; |
---|
[1359] | 235 | #ifdef CLP_INVESTIGATE2 |
---|
[1315] | 236 | printf("Number fixed increased from %d to %d\n", |
---|
| 237 | nFix - nFix2, nFix); |
---|
[1359] | 238 | #endif |
---|
[1315] | 239 | } |
---|
[1286] | 240 | if (nFix*10 > numberIntegers) { |
---|
| 241 | returnCode = smallBranchAndBound(newSolver, numberNodes_, newSolution, objectiveValue, |
---|
| 242 | objectiveValue, "CbcHeuristicLocal"); |
---|
[1364] | 243 | /* |
---|
| 244 | -2 is return due to user event, and -1 is overloaded with what look to be |
---|
| 245 | two contradictory meanings. |
---|
| 246 | */ |
---|
| 247 | if (returnCode < 0) { |
---|
[1286] | 248 | returnCode = 0; // returned on size |
---|
| 249 | int numberColumns = newSolver->getNumCols(); |
---|
| 250 | int numberContinuous = numberColumns - numberIntegers; |
---|
| 251 | if (numberContinuous > 2*numberIntegers && |
---|
| 252 | nFix*10 < numberColumns) { |
---|
[1271] | 253 | #define LOCAL_FIX_CONTINUOUS |
---|
| 254 | #ifdef LOCAL_FIX_CONTINUOUS |
---|
[1286] | 255 | //const double * colUpper = newSolver->getColUpper(); |
---|
| 256 | const double * colLower = newSolver->getColLower(); |
---|
| 257 | int nAtLb = 0; |
---|
| 258 | //double sumDj=0.0; |
---|
| 259 | const double * dj = newSolver->getReducedCost(); |
---|
| 260 | double direction = newSolver->getObjSense(); |
---|
| 261 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 262 | if (!newSolver->isInteger(iColumn)) { |
---|
| 263 | if (!used_[iColumn]) { |
---|
| 264 | //double djValue = dj[iColumn]*direction; |
---|
| 265 | nAtLb++; |
---|
| 266 | //sumDj += djValue; |
---|
| 267 | } |
---|
| 268 | } |
---|
| 269 | } |
---|
| 270 | if (nAtLb) { |
---|
| 271 | // fix some continuous |
---|
| 272 | double * sort = new double[nAtLb]; |
---|
| 273 | int * which = new int [nAtLb]; |
---|
| 274 | //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6); |
---|
| 275 | int nFix2 = 0; |
---|
| 276 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 277 | if (!newSolver->isInteger(iColumn)) { |
---|
| 278 | if (!used_[iColumn]) { |
---|
| 279 | double djValue = dj[iColumn] * direction; |
---|
| 280 | if (djValue > 1.0e-6) { |
---|
| 281 | sort[nFix2] = -djValue; |
---|
| 282 | which[nFix2++] = iColumn; |
---|
| 283 | } |
---|
| 284 | } |
---|
| 285 | } |
---|
| 286 | } |
---|
| 287 | CoinSort_2(sort, sort + nFix2, which); |
---|
| 288 | int divisor = 2; |
---|
| 289 | nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor); |
---|
| 290 | for (int i = 0; i < nFix2; i++) { |
---|
| 291 | int iColumn = which[i]; |
---|
| 292 | newSolver->setColUpper(iColumn, colLower[iColumn]); |
---|
| 293 | } |
---|
| 294 | delete [] sort; |
---|
| 295 | delete [] which; |
---|
[1271] | 296 | #ifdef CLP_INVESTIGATE2 |
---|
[1286] | 297 | printf("%d integers have zero value, and %d continuous fixed at lb\n", |
---|
| 298 | nFix, nFix2); |
---|
[1271] | 299 | #endif |
---|
[1286] | 300 | returnCode = smallBranchAndBound(newSolver, |
---|
| 301 | numberNodes_, newSolution, |
---|
| 302 | objectiveValue, |
---|
| 303 | objectiveValue, "CbcHeuristicLocal"); |
---|
| 304 | if (returnCode < 0) |
---|
| 305 | returnCode = 0; // returned on size |
---|
| 306 | } |
---|
[1271] | 307 | #endif |
---|
[1286] | 308 | } |
---|
| 309 | } |
---|
[1271] | 310 | } |
---|
[1364] | 311 | /* |
---|
| 312 | If the result is complete exploration with a solution (3) or proven |
---|
| 313 | infeasibility (2), we could generate a cut (the AI folks would call it a |
---|
| 314 | nogood) to prevent us from going down this route in the future. |
---|
| 315 | */ |
---|
[1286] | 316 | if ((returnCode&2) != 0) { |
---|
| 317 | // could add cut |
---|
| 318 | returnCode &= ~2; |
---|
| 319 | } |
---|
[197] | 320 | |
---|
[1286] | 321 | delete newSolver; |
---|
| 322 | return returnCode; |
---|
[175] | 323 | } |
---|
| 324 | /* |
---|
| 325 | First tries setting a variable to better value. If feasible then |
---|
| 326 | tries setting others. If not feasible then tries swaps |
---|
[1364] | 327 | Returns 1 if solution, 0 if not |
---|
| 328 | The main body of this routine implements an O((q^2)/2) brute force search |
---|
| 329 | around the current solution, for q = number of integer variables. Call this |
---|
| 330 | the inc/dec heuristic. For each integer variable x<i>, first decrement the |
---|
| 331 | value. Then, for integer variables x<i+1>, ..., x<q-1>, try increment and |
---|
| 332 | decrement. If one of these permutations produces a better solution, |
---|
| 333 | remember it. Then repeat, with x<i> incremented. If we find a better |
---|
| 334 | solution, update our notion of current solution and continue. |
---|
| 335 | |
---|
| 336 | The net effect is a greedy walk: As each improving pair is found, the |
---|
| 337 | current solution is updated and the search continues from this updated |
---|
| 338 | solution. |
---|
| 339 | |
---|
| 340 | Way down at the end, we call solutionFix, which will create a drastically |
---|
| 341 | restricted problem based on variables marked as used, then do mini-BaC on |
---|
| 342 | the restricted problem. This can occur even if we don't try the inc/dec |
---|
| 343 | heuristic. This would be more obvious if the inc/dec heuristic were broken |
---|
| 344 | out as a separate routine and solutionFix had a name that reflected where |
---|
| 345 | it was headed. |
---|
| 346 | |
---|
| 347 | The return code of 0 is grossly overloaded, because it maps to a return |
---|
| 348 | code of 0 from solutionFix, which is itself grossly overloaded. See |
---|
| 349 | comments in solutionFix and in CbcHeuristic::smallBranchAndBound. |
---|
| 350 | */ |
---|
[175] | 351 | int |
---|
| 352 | CbcHeuristicLocal::solution(double & solutionValue, |
---|
[1286] | 353 | double * betterSolution) |
---|
[175] | 354 | { |
---|
[1364] | 355 | /* |
---|
| 356 | Execute only if a new solution has been discovered since the last time we |
---|
| 357 | were called. |
---|
| 358 | */ |
---|
[175] | 359 | |
---|
[1286] | 360 | numCouldRun_++; |
---|
[1613] | 361 | // See if frequency kills off idea |
---|
| 362 | int swap = swap_%100; |
---|
| 363 | int skip = swap_/100; |
---|
| 364 | int nodeCount = model_->getNodeCount(); |
---|
| 365 | if (nodeCount<lastRunDeep_+skip && nodeCount != lastRunDeep_+1) |
---|
| 366 | return 0; |
---|
| 367 | if (numberSolutions_ == model_->getSolutionCount() && |
---|
| 368 | (numberSolutions_ == howOftenShallow_ || |
---|
| 369 | nodeCount < lastRunDeep_+2*skip)) |
---|
[1286] | 370 | return 0; |
---|
[1613] | 371 | howOftenShallow_ = numberSolutions_; |
---|
[1286] | 372 | numberSolutions_ = model_->getSolutionCount(); |
---|
[1613] | 373 | if (nodeCount<lastRunDeep_+skip ) |
---|
| 374 | return 0; |
---|
| 375 | lastRunDeep_ = nodeCount; |
---|
| 376 | howOftenShallow_ = numberSolutions_; |
---|
| 377 | |
---|
| 378 | if ((swap%10) == 2) { |
---|
| 379 | // try merge |
---|
| 380 | return solutionFix( solutionValue, betterSolution, NULL); |
---|
| 381 | } |
---|
[1364] | 382 | /* |
---|
| 383 | Exclude long (column), thin (row) systems. |
---|
| 384 | |
---|
| 385 | Given the n^2 nature of the search, more than 100,000 columns could get |
---|
| 386 | expensive. But I don't yet see the rationale for the second part of the |
---|
| 387 | condition (cols > 10*rows). And cost is proportional to number of integer |
---|
| 388 | variables --- shouldn't we use that? |
---|
| 389 | |
---|
| 390 | Why wait until we have more than one solution? |
---|
| 391 | */ |
---|
[1286] | 392 | if ((model_->getNumCols() > 100000 && model_->getNumCols() > |
---|
| 393 | 10*model_->getNumRows()) || numberSolutions_ <= 1) |
---|
| 394 | return 0; // probably not worth it |
---|
| 395 | // worth trying |
---|
[175] | 396 | |
---|
[1286] | 397 | OsiSolverInterface * solver = model_->solver(); |
---|
| 398 | const double * rowLower = solver->getRowLower(); |
---|
| 399 | const double * rowUpper = solver->getRowUpper(); |
---|
| 400 | const double * solution = model_->bestSolution(); |
---|
[1364] | 401 | /* |
---|
| 402 | Shouldn't this test be redundant if we've already checked that |
---|
| 403 | numberSolutions_ > 1? Stronger: shouldn't this be an assertion? |
---|
| 404 | */ |
---|
[1286] | 405 | if (!solution) |
---|
| 406 | return 0; // No solution found yet |
---|
| 407 | const double * objective = solver->getObjCoefficients(); |
---|
| 408 | double primalTolerance; |
---|
| 409 | solver->getDblParam(OsiPrimalTolerance, primalTolerance); |
---|
[175] | 410 | |
---|
[1286] | 411 | int numberRows = matrix_.getNumRows(); |
---|
[175] | 412 | |
---|
[1286] | 413 | int numberIntegers = model_->numberIntegers(); |
---|
| 414 | const int * integerVariable = model_->integerVariable(); |
---|
[175] | 415 | |
---|
[1286] | 416 | int i; |
---|
| 417 | double direction = solver->getObjSense(); |
---|
| 418 | double newSolutionValue = model_->getObjValue() * direction; |
---|
| 419 | int returnCode = 0; |
---|
| 420 | numRuns_++; |
---|
| 421 | // Column copy |
---|
| 422 | const double * element = matrix_.getElements(); |
---|
| 423 | const int * row = matrix_.getIndices(); |
---|
| 424 | const CoinBigIndex * columnStart = matrix_.getVectorStarts(); |
---|
| 425 | const int * columnLength = matrix_.getVectorLengths(); |
---|
| 426 | |
---|
| 427 | // Get solution array for heuristic solution |
---|
| 428 | int numberColumns = solver->getNumCols(); |
---|
| 429 | double * newSolution = new double [numberColumns]; |
---|
| 430 | memcpy(newSolution, solution, numberColumns*sizeof(double)); |
---|
[1271] | 431 | #ifdef LOCAL_FIX_CONTINUOUS |
---|
[1286] | 432 | // mark continuous used |
---|
| 433 | const double * columnLower = solver->getColLower(); |
---|
| 434 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 435 | if (!solver->isInteger(iColumn)) { |
---|
| 436 | if (solution[iColumn] > columnLower[iColumn] + 1.0e-8) |
---|
| 437 | used_[iColumn] = numberSolutions_; |
---|
| 438 | } |
---|
[1271] | 439 | } |
---|
| 440 | #endif |
---|
[175] | 441 | |
---|
[1286] | 442 | // way is 1 if down possible, 2 if up possible, 3 if both possible |
---|
| 443 | char * way = new char[numberIntegers]; |
---|
| 444 | // corrected costs |
---|
| 445 | double * cost = new double[numberIntegers]; |
---|
| 446 | // for array to mark infeasible rows after iColumn branch |
---|
| 447 | char * mark = new char[numberRows]; |
---|
| 448 | memset(mark, 0, numberRows); |
---|
| 449 | // space to save values so we don't introduce rounding errors |
---|
| 450 | double * save = new double[numberRows]; |
---|
[1364] | 451 | /* |
---|
| 452 | Force variables within their original bounds, then to the nearest integer. |
---|
| 453 | Overall, we seem to be prepared to cope with noninteger bounds. Is this |
---|
| 454 | necessary? Seems like we'd be better off to force the bounds to integrality |
---|
| 455 | as part of preprocessing. More generally, why do we need to do this? This |
---|
| 456 | solution should have been cleaned and checked when it was accepted as a |
---|
| 457 | solution! |
---|
[175] | 458 | |
---|
[1364] | 459 | Once the value is set, decide whether we can move up or down. |
---|
| 460 | |
---|
| 461 | The only place that used_ is used is in solutionFix; if a variable is not |
---|
| 462 | flagged as used, it will be fixed (at lower bound). Why the asymmetric |
---|
| 463 | treatment? This makes some sense for binary variables (for which there are |
---|
| 464 | only two options). But for general integer variables, why not make a similar |
---|
| 465 | test against the original upper bound? |
---|
| 466 | */ |
---|
| 467 | |
---|
[1286] | 468 | // clean solution |
---|
| 469 | for (i = 0; i < numberIntegers; i++) { |
---|
| 470 | int iColumn = integerVariable[i]; |
---|
| 471 | const OsiObject * object = model_->object(i); |
---|
| 472 | // get original bounds |
---|
| 473 | double originalLower; |
---|
| 474 | double originalUpper; |
---|
| 475 | getIntegerInformation( object, originalLower, originalUpper); |
---|
| 476 | double value = newSolution[iColumn]; |
---|
| 477 | if (value < originalLower) { |
---|
| 478 | value = originalLower; |
---|
| 479 | newSolution[iColumn] = value; |
---|
| 480 | } else if (value > originalUpper) { |
---|
| 481 | value = originalUpper; |
---|
| 482 | newSolution[iColumn] = value; |
---|
| 483 | } |
---|
| 484 | double nearest = floor(value + 0.5); |
---|
| 485 | //assert(fabs(value-nearest)<10.0*primalTolerance); |
---|
| 486 | value = nearest; |
---|
| 487 | newSolution[iColumn] = nearest; |
---|
| 488 | // if away from lower bound mark that fact |
---|
| 489 | if (nearest > originalLower) { |
---|
| 490 | used_[iColumn] = numberSolutions_; |
---|
| 491 | } |
---|
| 492 | cost[i] = direction * objective[iColumn]; |
---|
[1364] | 493 | /* |
---|
| 494 | Given previous computation we're checking that value is at least 1 away |
---|
| 495 | from the original bounds. |
---|
| 496 | */ |
---|
[1286] | 497 | int iway = 0; |
---|
| 498 | |
---|
| 499 | if (value > originalLower + 0.5) |
---|
| 500 | iway = 1; |
---|
| 501 | if (value < originalUpper - 0.5) |
---|
| 502 | iway |= 2; |
---|
| 503 | way[i] = static_cast<char>(iway); |
---|
[175] | 504 | } |
---|
[1364] | 505 | /* |
---|
| 506 | Calculate lhs of each constraint for groomed solution. |
---|
| 507 | */ |
---|
[1286] | 508 | // get row activities |
---|
| 509 | double * rowActivity = new double[numberRows]; |
---|
| 510 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
[175] | 511 | |
---|
[1286] | 512 | for (i = 0; i < numberColumns; i++) { |
---|
| 513 | int j; |
---|
| 514 | double value = newSolution[i]; |
---|
| 515 | if (value) { |
---|
| 516 | for (j = columnStart[i]; |
---|
| 517 | j < columnStart[i] + columnLength[i]; j++) { |
---|
| 518 | int iRow = row[j]; |
---|
| 519 | rowActivity[iRow] += value * element[j]; |
---|
| 520 | } |
---|
| 521 | } |
---|
[175] | 522 | } |
---|
[1364] | 523 | /* |
---|
| 524 | Check that constraints are satisfied. For small infeasibility, force the |
---|
| 525 | activity within bound. Again, why is this necessary if the current solution |
---|
| 526 | was accepted as a valid solution? |
---|
| 527 | |
---|
| 528 | Why are we scanning past the first unacceptable constraint? |
---|
| 529 | */ |
---|
[1286] | 530 | // check was feasible - if not adjust (cleaning may move) |
---|
| 531 | // if very infeasible then give up |
---|
| 532 | bool tryHeuristic = true; |
---|
| 533 | for (i = 0; i < numberRows; i++) { |
---|
| 534 | if (rowActivity[i] < rowLower[i]) { |
---|
| 535 | if (rowActivity[i] < rowLower[i] - 10.0*primalTolerance) |
---|
| 536 | tryHeuristic = false; |
---|
| 537 | rowActivity[i] = rowLower[i]; |
---|
| 538 | } else if (rowActivity[i] > rowUpper[i]) { |
---|
| 539 | if (rowActivity[i] < rowUpper[i] + 10.0*primalTolerance) |
---|
| 540 | tryHeuristic = false; |
---|
| 541 | rowActivity[i] = rowUpper[i]; |
---|
| 542 | } |
---|
[175] | 543 | } |
---|
[1364] | 544 | /* |
---|
| 545 | This bit of code is not quite totally redundant: it'll bail at 10,000 |
---|
| 546 | instead of 100,000. Potentially we can do a lot of work to get here, only |
---|
| 547 | to abandon it. |
---|
| 548 | */ |
---|
[1286] | 549 | // Switch off if may take too long |
---|
| 550 | if (model_->getNumCols() > 10000 && model_->getNumCols() > |
---|
| 551 | 10*model_->getNumRows()) |
---|
| 552 | tryHeuristic = false; |
---|
[1364] | 553 | /* |
---|
| 554 | Try the inc/dec heuristic? |
---|
| 555 | */ |
---|
[1286] | 556 | if (tryHeuristic) { |
---|
| 557 | |
---|
[1613] | 558 | // total change in objective |
---|
| 559 | double totalChange = 0.0; |
---|
| 560 | // local best change in objective |
---|
[1286] | 561 | double bestChange = 0.0; |
---|
[1613] | 562 | // maybe just do 1000 |
---|
| 563 | int maxIntegers = numberIntegers; |
---|
| 564 | if (((swap/10) &1) != 0) { |
---|
| 565 | maxIntegers = CoinMin(1000,numberIntegers); |
---|
| 566 | } |
---|
[1364] | 567 | /* |
---|
| 568 | Outer loop to walk integer variables. Call the current variable x<i>. At the |
---|
| 569 | end of this loop, bestChange will contain the best (negative) change in the |
---|
| 570 | objective for any single pair. |
---|
[1286] | 571 | |
---|
[1364] | 572 | The trouble is, we're limited to monotonically increasing improvement. |
---|
| 573 | Suppose we discover an improvement of 10 for some pair. If, later in the |
---|
| 574 | search, we discover an improvement of 9 for some other pair, we will not use |
---|
| 575 | it. That seems wasteful. |
---|
| 576 | */ |
---|
| 577 | |
---|
[1286] | 578 | for (i = 0; i < numberIntegers; i++) { |
---|
| 579 | int iColumn = integerVariable[i]; |
---|
[1613] | 580 | bestChange = 0.0; |
---|
| 581 | int endInner = CoinMin(numberIntegers,i+maxIntegers); |
---|
[1286] | 582 | |
---|
| 583 | double objectiveCoefficient = cost[i]; |
---|
| 584 | int k; |
---|
| 585 | int j; |
---|
| 586 | int goodK = -1; |
---|
| 587 | int wayK = -1, wayI = -1; |
---|
[1364] | 588 | /* |
---|
| 589 | Try decrementing x<i>. |
---|
| 590 | */ |
---|
[1286] | 591 | if ((way[i]&1) != 0) { |
---|
| 592 | int numberInfeasible = 0; |
---|
[1364] | 593 | /* |
---|
| 594 | Adjust row activities where x<i> has a nonzero coefficient. Save the old |
---|
| 595 | values for restoration. Mark any rows that become infeasible as a result |
---|
| 596 | of the decrement. |
---|
| 597 | */ |
---|
[1286] | 598 | // save row activities and adjust |
---|
| 599 | for (j = columnStart[iColumn]; |
---|
| 600 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 601 | int iRow = row[j]; |
---|
| 602 | save[iRow] = rowActivity[iRow]; |
---|
| 603 | rowActivity[iRow] -= element[j]; |
---|
| 604 | if (rowActivity[iRow] < rowLower[iRow] - primalTolerance || |
---|
| 605 | rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { |
---|
| 606 | // mark row |
---|
| 607 | mark[iRow] = 1; |
---|
| 608 | numberInfeasible++; |
---|
| 609 | } |
---|
| 610 | } |
---|
[1364] | 611 | /* |
---|
| 612 | Run through the remaining integer variables. Try increment and decrement on |
---|
| 613 | each one. If the potential objective change is better than anything we've |
---|
| 614 | seen so far, do a full evaluation of x<k> in that direction. If we can |
---|
| 615 | repair all infeasibilities introduced by pushing x<i> down, we have a |
---|
| 616 | winner. Remember the best variable, and the direction for x<i> and x<k>. |
---|
| 617 | */ |
---|
| 618 | // try down |
---|
[1613] | 619 | for (k = i + 1; k < endInner; k++) { |
---|
[1286] | 620 | if ((way[k]&1) != 0) { |
---|
| 621 | // try down |
---|
| 622 | if (-objectiveCoefficient - cost[k] < bestChange) { |
---|
| 623 | // see if feasible down |
---|
| 624 | bool good = true; |
---|
| 625 | int numberMarked = 0; |
---|
| 626 | int kColumn = integerVariable[k]; |
---|
| 627 | for (j = columnStart[kColumn]; |
---|
| 628 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
---|
| 629 | int iRow = row[j]; |
---|
| 630 | double newValue = rowActivity[iRow] - element[j]; |
---|
| 631 | if (newValue < rowLower[iRow] - primalTolerance || |
---|
| 632 | newValue > rowUpper[iRow] + primalTolerance) { |
---|
| 633 | good = false; |
---|
| 634 | break; |
---|
| 635 | } else if (mark[iRow]) { |
---|
| 636 | // made feasible |
---|
| 637 | numberMarked++; |
---|
| 638 | } |
---|
| 639 | } |
---|
| 640 | if (good && numberMarked == numberInfeasible) { |
---|
| 641 | // better solution |
---|
| 642 | goodK = k; |
---|
| 643 | wayK = -1; |
---|
| 644 | wayI = -1; |
---|
| 645 | bestChange = -objectiveCoefficient - cost[k]; |
---|
| 646 | } |
---|
| 647 | } |
---|
| 648 | } |
---|
| 649 | if ((way[k]&2) != 0) { |
---|
| 650 | // try up |
---|
| 651 | if (-objectiveCoefficient + cost[k] < bestChange) { |
---|
| 652 | // see if feasible up |
---|
| 653 | bool good = true; |
---|
| 654 | int numberMarked = 0; |
---|
| 655 | int kColumn = integerVariable[k]; |
---|
| 656 | for (j = columnStart[kColumn]; |
---|
| 657 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
---|
| 658 | int iRow = row[j]; |
---|
| 659 | double newValue = rowActivity[iRow] + element[j]; |
---|
| 660 | if (newValue < rowLower[iRow] - primalTolerance || |
---|
| 661 | newValue > rowUpper[iRow] + primalTolerance) { |
---|
| 662 | good = false; |
---|
| 663 | break; |
---|
| 664 | } else if (mark[iRow]) { |
---|
| 665 | // made feasible |
---|
| 666 | numberMarked++; |
---|
| 667 | } |
---|
| 668 | } |
---|
| 669 | if (good && numberMarked == numberInfeasible) { |
---|
| 670 | // better solution |
---|
| 671 | goodK = k; |
---|
| 672 | wayK = 1; |
---|
| 673 | wayI = -1; |
---|
| 674 | bestChange = -objectiveCoefficient + cost[k]; |
---|
| 675 | } |
---|
| 676 | } |
---|
| 677 | } |
---|
| 678 | } |
---|
[1364] | 679 | /* |
---|
| 680 | Remove effect of decrementing x<i> by restoring original lhs values. |
---|
| 681 | */ |
---|
[1286] | 682 | // restore row activities |
---|
| 683 | for (j = columnStart[iColumn]; |
---|
| 684 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 685 | int iRow = row[j]; |
---|
| 686 | rowActivity[iRow] = save[iRow]; |
---|
| 687 | mark[iRow] = 0; |
---|
| 688 | } |
---|
| 689 | } |
---|
[1364] | 690 | /* |
---|
| 691 | Try to increment x<i>. Actions as for decrement. |
---|
| 692 | */ |
---|
[1286] | 693 | if ((way[i]&2) != 0) { |
---|
| 694 | int numberInfeasible = 0; |
---|
| 695 | // save row activities and adjust |
---|
| 696 | for (j = columnStart[iColumn]; |
---|
| 697 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 698 | int iRow = row[j]; |
---|
| 699 | save[iRow] = rowActivity[iRow]; |
---|
| 700 | rowActivity[iRow] += element[j]; |
---|
| 701 | if (rowActivity[iRow] < rowLower[iRow] - primalTolerance || |
---|
| 702 | rowActivity[iRow] > rowUpper[iRow] + primalTolerance) { |
---|
| 703 | // mark row |
---|
| 704 | mark[iRow] = 1; |
---|
| 705 | numberInfeasible++; |
---|
| 706 | } |
---|
| 707 | } |
---|
| 708 | // try up |
---|
[1613] | 709 | for (k = i + 1; k < endInner; k++) { |
---|
[1286] | 710 | if ((way[k]&1) != 0) { |
---|
| 711 | // try down |
---|
| 712 | if (objectiveCoefficient - cost[k] < bestChange) { |
---|
| 713 | // see if feasible down |
---|
| 714 | bool good = true; |
---|
| 715 | int numberMarked = 0; |
---|
| 716 | int kColumn = integerVariable[k]; |
---|
| 717 | for (j = columnStart[kColumn]; |
---|
| 718 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
---|
| 719 | int iRow = row[j]; |
---|
| 720 | double newValue = rowActivity[iRow] - element[j]; |
---|
| 721 | if (newValue < rowLower[iRow] - primalTolerance || |
---|
| 722 | newValue > rowUpper[iRow] + primalTolerance) { |
---|
| 723 | good = false; |
---|
| 724 | break; |
---|
| 725 | } else if (mark[iRow]) { |
---|
| 726 | // made feasible |
---|
| 727 | numberMarked++; |
---|
| 728 | } |
---|
| 729 | } |
---|
| 730 | if (good && numberMarked == numberInfeasible) { |
---|
| 731 | // better solution |
---|
| 732 | goodK = k; |
---|
| 733 | wayK = -1; |
---|
| 734 | wayI = 1; |
---|
| 735 | bestChange = objectiveCoefficient - cost[k]; |
---|
| 736 | } |
---|
| 737 | } |
---|
| 738 | } |
---|
| 739 | if ((way[k]&2) != 0) { |
---|
| 740 | // try up |
---|
| 741 | if (objectiveCoefficient + cost[k] < bestChange) { |
---|
| 742 | // see if feasible up |
---|
| 743 | bool good = true; |
---|
| 744 | int numberMarked = 0; |
---|
| 745 | int kColumn = integerVariable[k]; |
---|
| 746 | for (j = columnStart[kColumn]; |
---|
| 747 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
---|
| 748 | int iRow = row[j]; |
---|
| 749 | double newValue = rowActivity[iRow] + element[j]; |
---|
| 750 | if (newValue < rowLower[iRow] - primalTolerance || |
---|
| 751 | newValue > rowUpper[iRow] + primalTolerance) { |
---|
| 752 | good = false; |
---|
| 753 | break; |
---|
| 754 | } else if (mark[iRow]) { |
---|
| 755 | // made feasible |
---|
| 756 | numberMarked++; |
---|
| 757 | } |
---|
| 758 | } |
---|
| 759 | if (good && numberMarked == numberInfeasible) { |
---|
| 760 | // better solution |
---|
| 761 | goodK = k; |
---|
| 762 | wayK = 1; |
---|
| 763 | wayI = 1; |
---|
| 764 | bestChange = objectiveCoefficient + cost[k]; |
---|
| 765 | } |
---|
| 766 | } |
---|
| 767 | } |
---|
| 768 | } |
---|
| 769 | // restore row activities |
---|
| 770 | for (j = columnStart[iColumn]; |
---|
| 771 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 772 | int iRow = row[j]; |
---|
| 773 | rowActivity[iRow] = save[iRow]; |
---|
| 774 | mark[iRow] = 0; |
---|
| 775 | } |
---|
| 776 | } |
---|
[1364] | 777 | /* |
---|
| 778 | We've found a pair x<i> and x<k> which produce a better solution. Update our |
---|
| 779 | notion of current solution to match. |
---|
| 780 | |
---|
| 781 | Why does this not update newSolutionValue? |
---|
| 782 | */ |
---|
[1286] | 783 | if (goodK >= 0) { |
---|
| 784 | // we found something - update solution |
---|
| 785 | for (j = columnStart[iColumn]; |
---|
| 786 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
---|
| 787 | int iRow = row[j]; |
---|
| 788 | rowActivity[iRow] += wayI * element[j]; |
---|
| 789 | } |
---|
| 790 | newSolution[iColumn] += wayI; |
---|
| 791 | int kColumn = integerVariable[goodK]; |
---|
| 792 | for (j = columnStart[kColumn]; |
---|
| 793 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
---|
| 794 | int iRow = row[j]; |
---|
| 795 | rowActivity[iRow] += wayK * element[j]; |
---|
| 796 | } |
---|
| 797 | newSolution[kColumn] += wayK; |
---|
[1364] | 798 | /* |
---|
| 799 | Adjust motion range for x<k>. We may have banged up against a bound with that |
---|
| 800 | last move. |
---|
| 801 | */ |
---|
| 802 | // See if k can go further ? |
---|
[1286] | 803 | const OsiObject * object = model_->object(goodK); |
---|
| 804 | // get original bounds |
---|
| 805 | double originalLower; |
---|
| 806 | double originalUpper; |
---|
| 807 | getIntegerInformation( object, originalLower, originalUpper); |
---|
| 808 | |
---|
| 809 | double value = newSolution[kColumn]; |
---|
| 810 | int iway = 0; |
---|
| 811 | |
---|
| 812 | if (value > originalLower + 0.5) |
---|
| 813 | iway = 1; |
---|
| 814 | if (value < originalUpper - 0.5) |
---|
| 815 | iway |= 2; |
---|
| 816 | way[goodK] = static_cast<char>(iway); |
---|
[1613] | 817 | totalChange += bestChange; |
---|
[1286] | 818 | } |
---|
[175] | 819 | } |
---|
[1364] | 820 | /* |
---|
| 821 | End of loop to try increment/decrement of integer variables. |
---|
| 822 | |
---|
| 823 | newSolutionValue does not necessarily match the current newSolution, and |
---|
| 824 | bestChange simply reflects the best single change. Still, that's sufficient |
---|
| 825 | to indicate that there's been at least one change. Check that we really do |
---|
| 826 | have a valid solution. |
---|
| 827 | */ |
---|
[1613] | 828 | if (totalChange + newSolutionValue < solutionValue) { |
---|
[1286] | 829 | // paranoid check |
---|
| 830 | memset(rowActivity, 0, numberRows*sizeof(double)); |
---|
| 831 | |
---|
| 832 | for (i = 0; i < numberColumns; i++) { |
---|
| 833 | int j; |
---|
| 834 | double value = newSolution[i]; |
---|
| 835 | if (value) { |
---|
| 836 | for (j = columnStart[i]; |
---|
| 837 | j < columnStart[i] + columnLength[i]; j++) { |
---|
| 838 | int iRow = row[j]; |
---|
| 839 | rowActivity[iRow] += value * element[j]; |
---|
| 840 | } |
---|
| 841 | } |
---|
| 842 | } |
---|
| 843 | int numberBad = 0; |
---|
| 844 | double sumBad = 0.0; |
---|
| 845 | // check was approximately feasible |
---|
| 846 | for (i = 0; i < numberRows; i++) { |
---|
| 847 | if (rowActivity[i] < rowLower[i]) { |
---|
| 848 | sumBad += rowLower[i] - rowActivity[i]; |
---|
| 849 | if (rowActivity[i] < rowLower[i] - 10.0*primalTolerance) |
---|
| 850 | numberBad++; |
---|
| 851 | } else if (rowActivity[i] > rowUpper[i]) { |
---|
| 852 | sumBad += rowUpper[i] - rowActivity[i]; |
---|
| 853 | if (rowActivity[i] > rowUpper[i] + 10.0*primalTolerance) |
---|
| 854 | numberBad++; |
---|
| 855 | } |
---|
| 856 | } |
---|
| 857 | if (!numberBad) { |
---|
| 858 | for (i = 0; i < numberIntegers; i++) { |
---|
| 859 | int iColumn = integerVariable[i]; |
---|
| 860 | const OsiObject * object = model_->object(i); |
---|
| 861 | // get original bounds |
---|
| 862 | double originalLower; |
---|
| 863 | double originalUpper; |
---|
| 864 | getIntegerInformation( object, originalLower, originalUpper); |
---|
| 865 | |
---|
| 866 | double value = newSolution[iColumn]; |
---|
| 867 | // if away from lower bound mark that fact |
---|
| 868 | if (value > originalLower) { |
---|
| 869 | used_[iColumn] = numberSolutions_; |
---|
| 870 | } |
---|
| 871 | } |
---|
[1364] | 872 | /* |
---|
| 873 | Copy the solution to the array returned to the client. Grab a basis from |
---|
| 874 | the solver (which, if it exists, is almost certainly infeasible, but it |
---|
| 875 | should be ok for a dual start). The value returned as solutionValue is |
---|
| 876 | conservative because of handling of newSolutionValue and bestChange, as |
---|
| 877 | described above. |
---|
| 878 | */ |
---|
[1286] | 879 | // new solution |
---|
| 880 | memcpy(betterSolution, newSolution, numberColumns*sizeof(double)); |
---|
| 881 | CoinWarmStartBasis * basis = |
---|
| 882 | dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; |
---|
| 883 | if (basis) { |
---|
| 884 | model_->setBestSolutionBasis(* basis); |
---|
| 885 | delete basis; |
---|
| 886 | } |
---|
| 887 | returnCode = 1; |
---|
| 888 | solutionValue = newSolutionValue + bestChange; |
---|
| 889 | } else { |
---|
| 890 | // bad solution - should not happen so debug if see message |
---|
[1641] | 891 | COIN_DETAIL_PRINT(printf("Local search got bad solution with %d infeasibilities summing to %g\n", |
---|
| 892 | numberBad, sumBad)); |
---|
[1286] | 893 | } |
---|
| 894 | } |
---|
[175] | 895 | } |
---|
[1364] | 896 | /* |
---|
| 897 | We're done. Clean up. |
---|
| 898 | */ |
---|
[1286] | 899 | delete [] newSolution; |
---|
| 900 | delete [] rowActivity; |
---|
| 901 | delete [] way; |
---|
| 902 | delete [] cost; |
---|
| 903 | delete [] save; |
---|
| 904 | delete [] mark; |
---|
[1364] | 905 | /* |
---|
| 906 | Do we want to try swapping values between solutions? |
---|
| 907 | swap_ is set elsewhere; it's not adjusted during heuristic execution. |
---|
| 908 | |
---|
| 909 | Again, redundant test. We shouldn't be here if numberSolutions_ = 1. |
---|
| 910 | */ |
---|
[1613] | 911 | if (numberSolutions_ > 1 && (swap%10) == 1) { |
---|
[1286] | 912 | // try merge |
---|
| 913 | int returnCode2 = solutionFix( solutionValue, betterSolution, NULL); |
---|
| 914 | if (returnCode2) |
---|
| 915 | returnCode = 1; |
---|
| 916 | } |
---|
| 917 | return returnCode; |
---|
[175] | 918 | } |
---|
| 919 | // update model |
---|
| 920 | void CbcHeuristicLocal::setModel(CbcModel * model) |
---|
| 921 | { |
---|
[1286] | 922 | model_ = model; |
---|
| 923 | // Get a copy of original matrix |
---|
| 924 | assert(model_->solver()); |
---|
| 925 | if (model_->solver()->getNumRows()) { |
---|
| 926 | matrix_ = *model_->solver()->getMatrixByCol(); |
---|
| 927 | } |
---|
| 928 | delete [] used_; |
---|
| 929 | int numberColumns = model->solver()->getNumCols(); |
---|
| 930 | used_ = new int[numberColumns]; |
---|
| 931 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
[175] | 932 | } |
---|
[1802] | 933 | |
---|
[1100] | 934 | // Default Constructor |
---|
[1802] | 935 | CbcHeuristicProximity::CbcHeuristicProximity() |
---|
| 936 | : CbcHeuristic() |
---|
| 937 | { |
---|
| 938 | feasibilityPump_ = NULL; |
---|
| 939 | numberSolutions_ = 0; |
---|
| 940 | used_ = NULL; |
---|
| 941 | lastRunDeep_ = -1000000; |
---|
| 942 | switches_ |= 16; // needs a new solution |
---|
| 943 | } |
---|
| 944 | |
---|
| 945 | // Constructor with model - assumed before cuts |
---|
| 946 | |
---|
| 947 | CbcHeuristicProximity::CbcHeuristicProximity(CbcModel & model) |
---|
| 948 | : CbcHeuristic(model) |
---|
| 949 | { |
---|
| 950 | feasibilityPump_ = NULL; |
---|
| 951 | numberSolutions_ = 0; |
---|
| 952 | lastRunDeep_ = -1000000; |
---|
| 953 | switches_ |= 16; // needs a new solution |
---|
| 954 | int numberColumns = model.solver()->getNumCols(); |
---|
| 955 | used_ = new int[numberColumns]; |
---|
| 956 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
| 957 | } |
---|
| 958 | |
---|
| 959 | // Destructor |
---|
| 960 | CbcHeuristicProximity::~CbcHeuristicProximity () |
---|
| 961 | { |
---|
| 962 | delete feasibilityPump_; |
---|
| 963 | delete [] used_; |
---|
| 964 | } |
---|
| 965 | |
---|
| 966 | // Clone |
---|
| 967 | CbcHeuristic * |
---|
| 968 | CbcHeuristicProximity::clone() const |
---|
| 969 | { |
---|
| 970 | return new CbcHeuristicProximity(*this); |
---|
| 971 | } |
---|
| 972 | // Create C++ lines to get to current state |
---|
| 973 | void |
---|
| 974 | CbcHeuristicProximity::generateCpp( FILE * fp) |
---|
| 975 | { |
---|
| 976 | CbcHeuristicProximity other; |
---|
| 977 | fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n"); |
---|
| 978 | fprintf(fp, "3 CbcHeuristicProximity heuristicProximity(*cbcModel);\n"); |
---|
| 979 | CbcHeuristic::generateCpp(fp, "heuristicProximity"); |
---|
| 980 | fprintf(fp, "3 cbcModel->addHeuristic(&heuristicProximity);\n"); |
---|
| 981 | } |
---|
| 982 | |
---|
| 983 | // Copy constructor |
---|
| 984 | CbcHeuristicProximity::CbcHeuristicProximity(const CbcHeuristicProximity & rhs) |
---|
| 985 | : |
---|
| 986 | CbcHeuristic(rhs), |
---|
| 987 | numberSolutions_(rhs.numberSolutions_) |
---|
| 988 | { |
---|
| 989 | feasibilityPump_ = NULL; |
---|
| 990 | if (model_ && rhs.used_) { |
---|
| 991 | int numberColumns = model_->solver()->getNumCols(); |
---|
| 992 | used_ = CoinCopyOfArray(rhs.used_, numberColumns); |
---|
| 993 | if (rhs.feasibilityPump_) |
---|
| 994 | feasibilityPump_ = new CbcHeuristicFPump(*rhs.feasibilityPump_); |
---|
| 995 | } else { |
---|
| 996 | used_ = NULL; |
---|
| 997 | } |
---|
| 998 | } |
---|
| 999 | |
---|
| 1000 | // Assignment operator |
---|
| 1001 | CbcHeuristicProximity & |
---|
| 1002 | CbcHeuristicProximity::operator=( const CbcHeuristicProximity & rhs) |
---|
| 1003 | { |
---|
| 1004 | if (this != &rhs) { |
---|
| 1005 | CbcHeuristic::operator=(rhs); |
---|
| 1006 | numberSolutions_ = rhs.numberSolutions_; |
---|
| 1007 | delete [] used_; |
---|
| 1008 | delete feasibilityPump_; |
---|
| 1009 | feasibilityPump_ = NULL; |
---|
| 1010 | if (model_ && rhs.used_) { |
---|
| 1011 | int numberColumns = model_->solver()->getNumCols(); |
---|
| 1012 | used_ = CoinCopyOfArray(rhs.used_, numberColumns); |
---|
| 1013 | if (rhs.feasibilityPump_) |
---|
| 1014 | feasibilityPump_ = new CbcHeuristicFPump(*rhs.feasibilityPump_); |
---|
| 1015 | } else { |
---|
| 1016 | used_ = NULL; |
---|
| 1017 | } |
---|
| 1018 | } |
---|
| 1019 | return *this; |
---|
| 1020 | } |
---|
| 1021 | |
---|
| 1022 | // Resets stuff if model changes |
---|
| 1023 | void |
---|
| 1024 | CbcHeuristicProximity::resetModel(CbcModel * /*model*/) |
---|
| 1025 | { |
---|
| 1026 | //CbcHeuristic::resetModel(model); |
---|
| 1027 | delete [] used_; |
---|
| 1028 | if (model_ && used_) { |
---|
| 1029 | int numberColumns = model_->solver()->getNumCols(); |
---|
| 1030 | used_ = new int[numberColumns]; |
---|
| 1031 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
| 1032 | } else { |
---|
| 1033 | used_ = NULL; |
---|
| 1034 | } |
---|
| 1035 | } |
---|
| 1036 | /* |
---|
| 1037 | Run a mini-BaB search after changing objective |
---|
| 1038 | |
---|
| 1039 | Return values are: |
---|
| 1040 | 1: smallBranchAndBound found a solution |
---|
| 1041 | 0: everything else |
---|
| 1042 | |
---|
| 1043 | The degree of overload as return codes from smallBranchAndBound are folded |
---|
| 1044 | into 0 is such that it's impossible to distinguish return codes that really |
---|
| 1045 | require attention from a simple `nothing of interest'. |
---|
| 1046 | */ |
---|
| 1047 | int |
---|
| 1048 | CbcHeuristicProximity::solution(double & solutionValue, |
---|
| 1049 | double * betterSolution) |
---|
| 1050 | { |
---|
| 1051 | if (feasibilityPumpOptions_ == -3 && numCouldRun_==0 && |
---|
| 1052 | !feasibilityPump_ ) { |
---|
| 1053 | // clone feasibility pump |
---|
| 1054 | for (int i = 0; i < model_->numberHeuristics(); i++) { |
---|
| 1055 | const CbcHeuristicFPump* pump = |
---|
| 1056 | dynamic_cast<const CbcHeuristicFPump*>(model_->heuristic(i)); |
---|
| 1057 | if (pump) { |
---|
| 1058 | feasibilityPump_ = new CbcHeuristicFPump(*pump); |
---|
| 1059 | break; |
---|
| 1060 | } |
---|
| 1061 | } |
---|
| 1062 | } |
---|
| 1063 | /* |
---|
| 1064 | Execute only if a new solution has been discovered since the last time we |
---|
| 1065 | were called. |
---|
| 1066 | */ |
---|
| 1067 | |
---|
| 1068 | numCouldRun_++; |
---|
[1820] | 1069 | if (!model_->bestSolution()) |
---|
| 1070 | return 0; // odd - because in parallel mode |
---|
[1802] | 1071 | int nodeCount = model_->getNodeCount(); |
---|
| 1072 | if (numberSolutions_ == model_->getSolutionCount()) |
---|
| 1073 | return 0; |
---|
| 1074 | numberSolutions_ = model_->getSolutionCount(); |
---|
| 1075 | lastRunDeep_ = nodeCount; |
---|
| 1076 | numRuns_++; |
---|
| 1077 | //howOftenShallow_ = numberSolutions_; |
---|
| 1078 | |
---|
| 1079 | /* |
---|
| 1080 | Load up a new solver with the solution. |
---|
| 1081 | |
---|
| 1082 | Why continuousSolver(), as opposed to solver()? |
---|
| 1083 | */ |
---|
| 1084 | OsiSolverInterface * newSolver = model_->continuousSolver()->clone(); |
---|
| 1085 | int numberColumns=newSolver->getNumCols(); |
---|
| 1086 | double * obj = CoinCopyOfArray(newSolver->getObjCoefficients(),numberColumns); |
---|
| 1087 | int * indices = new int [numberColumns]; |
---|
| 1088 | int n=0; |
---|
| 1089 | for (int i=0;i<numberColumns;i++) { |
---|
| 1090 | if (obj[i]) { |
---|
| 1091 | indices[n]=i; |
---|
| 1092 | obj[n++]=obj[i]; |
---|
| 1093 | } |
---|
| 1094 | } |
---|
| 1095 | double cutoff=model_->getCutoff(); |
---|
| 1096 | assert (cutoff<1.0e20); |
---|
| 1097 | if (model_->getCutoffIncrement()<1.0e-4) |
---|
| 1098 | cutoff -= 0.01; |
---|
| 1099 | double offset; |
---|
| 1100 | newSolver->getDblParam(OsiObjOffset, offset); |
---|
| 1101 | newSolver->setDblParam(OsiObjOffset, 0.0); |
---|
| 1102 | newSolver->addRow(n,indices,obj,-COIN_DBL_MAX,cutoff+offset); |
---|
| 1103 | delete [] indices; |
---|
| 1104 | memset(obj,0,numberColumns*sizeof(double)); |
---|
| 1105 | newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e20); |
---|
| 1106 | int numberIntegers = model_->numberIntegers(); |
---|
| 1107 | const int * integerVariable = model_->integerVariable(); |
---|
| 1108 | const double * solutionIn = model_->bestSolution(); |
---|
| 1109 | for (int i = 0; i < numberIntegers; i++) { |
---|
| 1110 | int iColumn = integerVariable[i]; |
---|
| 1111 | if (fabs(solutionIn[iColumn])<1.0e-5) |
---|
| 1112 | obj[iColumn]=1.0; |
---|
| 1113 | else if (fabs(solutionIn[iColumn]-1.0)<1.0e-5) |
---|
| 1114 | obj[iColumn]=-1.0; |
---|
| 1115 | } |
---|
| 1116 | newSolver->setObjective(obj); |
---|
| 1117 | delete [] obj; |
---|
| 1118 | //newSolver->writeMps("xxxx"); |
---|
| 1119 | char proxPrint[200]; |
---|
| 1120 | sprintf(proxPrint,"Running proximity search for %d nodes",numberNodes_); |
---|
| 1121 | model_->messageHandler()->message(CBC_FPUMP1, model_->messages()) |
---|
| 1122 | << proxPrint |
---|
| 1123 | << CoinMessageEol; |
---|
| 1124 | int maxSolutions = model_->getMaximumSolutions(); |
---|
| 1125 | model_->setMaximumSolutions(1); |
---|
| 1126 | bool pumpAdded = false; |
---|
| 1127 | if (feasibilityPumpOptions_ == -3 && feasibilityPump_) { |
---|
| 1128 | // add back feasibility pump |
---|
| 1129 | pumpAdded = true; |
---|
| 1130 | for (int i = 0; i < model_->numberHeuristics(); i++) { |
---|
| 1131 | const CbcHeuristicFPump* pump = |
---|
| 1132 | dynamic_cast<const CbcHeuristicFPump*>(model_->heuristic(i)); |
---|
| 1133 | if (pump) { |
---|
| 1134 | pumpAdded = false; |
---|
| 1135 | break; |
---|
| 1136 | } |
---|
| 1137 | } |
---|
| 1138 | if (pumpAdded) |
---|
| 1139 | model_->addHeuristic(feasibilityPump_); |
---|
| 1140 | } |
---|
| 1141 | int returnCode = |
---|
| 1142 | smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue, |
---|
| 1143 | 1.0e20, "CbcHeuristicProximity"); |
---|
| 1144 | if (pumpAdded) { |
---|
| 1145 | // take off feasibility pump |
---|
| 1146 | int lastHeuristic = model_->numberHeuristics()-1; |
---|
| 1147 | model_->setNumberHeuristics(lastHeuristic); |
---|
| 1148 | delete model_->heuristic(lastHeuristic); |
---|
| 1149 | } |
---|
| 1150 | model_->setMaximumSolutions(maxSolutions); |
---|
| 1151 | /* |
---|
| 1152 | -2 is return due to user event, and -1 is overloaded with what look to be |
---|
| 1153 | two contradictory meanings. |
---|
| 1154 | */ |
---|
| 1155 | if (returnCode < 0) { |
---|
| 1156 | returnCode = 0; |
---|
| 1157 | } |
---|
| 1158 | /* |
---|
| 1159 | If the result is complete exploration with a solution (3) or proven |
---|
| 1160 | infeasibility (2), we could generate a cut (the AI folks would call it a |
---|
| 1161 | nogood) to prevent us from going down this route in the future. |
---|
| 1162 | */ |
---|
| 1163 | if ((returnCode&2) != 0) { |
---|
| 1164 | // could add cut |
---|
| 1165 | returnCode &= ~2; |
---|
| 1166 | } |
---|
| 1167 | if ((returnCode&1) != 0) { |
---|
| 1168 | // redo objective |
---|
| 1169 | const double * obj = model_->continuousSolver()->getObjCoefficients(); |
---|
| 1170 | solutionValue = - offset; |
---|
| 1171 | for (int i=0;i<numberColumns;i++) { |
---|
| 1172 | solutionValue += obj[i]*betterSolution[i]; |
---|
| 1173 | } |
---|
| 1174 | } |
---|
| 1175 | |
---|
| 1176 | delete newSolver; |
---|
| 1177 | return returnCode; |
---|
| 1178 | } |
---|
| 1179 | // update model |
---|
| 1180 | void CbcHeuristicProximity::setModel(CbcModel * model) |
---|
| 1181 | { |
---|
| 1182 | model_ = model; |
---|
| 1183 | // Get a copy of original matrix |
---|
| 1184 | assert(model_->solver()); |
---|
| 1185 | delete [] used_; |
---|
| 1186 | int numberColumns = model->solver()->getNumCols(); |
---|
| 1187 | used_ = new int[numberColumns]; |
---|
| 1188 | memset(used_, 0, numberColumns*sizeof(int)); |
---|
| 1189 | } |
---|
| 1190 | |
---|
| 1191 | // Default Constructor |
---|
[1286] | 1192 | CbcHeuristicNaive::CbcHeuristicNaive() |
---|
| 1193 | : CbcHeuristic() |
---|
[1100] | 1194 | { |
---|
[1286] | 1195 | large_ = 1.0e6; |
---|
[1100] | 1196 | } |
---|
[175] | 1197 | |
---|
[1100] | 1198 | // Constructor with model - assumed before cuts |
---|
| 1199 | |
---|
| 1200 | CbcHeuristicNaive::CbcHeuristicNaive(CbcModel & model) |
---|
[1286] | 1201 | : CbcHeuristic(model) |
---|
[1100] | 1202 | { |
---|
[1286] | 1203 | large_ = 1.0e6; |
---|
[1100] | 1204 | } |
---|
| 1205 | |
---|
[1286] | 1206 | // Destructor |
---|
[1100] | 1207 | CbcHeuristicNaive::~CbcHeuristicNaive () |
---|
| 1208 | { |
---|
| 1209 | } |
---|
| 1210 | |
---|
| 1211 | // Clone |
---|
| 1212 | CbcHeuristic * |
---|
| 1213 | CbcHeuristicNaive::clone() const |
---|
| 1214 | { |
---|
[1286] | 1215 | return new CbcHeuristicNaive(*this); |
---|
[1100] | 1216 | } |
---|
| 1217 | // Create C++ lines to get to current state |
---|
[1286] | 1218 | void |
---|
| 1219 | CbcHeuristicNaive::generateCpp( FILE * fp) |
---|
[1100] | 1220 | { |
---|
[1286] | 1221 | CbcHeuristicNaive other; |
---|
[1802] | 1222 | fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n"); |
---|
[1286] | 1223 | fprintf(fp, "3 CbcHeuristicNaive naive(*cbcModel);\n"); |
---|
| 1224 | CbcHeuristic::generateCpp(fp, "naive"); |
---|
| 1225 | if (large_ != other.large_) |
---|
| 1226 | fprintf(fp, "3 naive.setLarge(%g);\n", large_); |
---|
| 1227 | else |
---|
| 1228 | fprintf(fp, "4 naive.setLarge(%g);\n", large_); |
---|
| 1229 | fprintf(fp, "3 cbcModel->addHeuristic(&naive);\n"); |
---|
[1100] | 1230 | } |
---|
| 1231 | |
---|
[1286] | 1232 | // Copy constructor |
---|
[1100] | 1233 | CbcHeuristicNaive::CbcHeuristicNaive(const CbcHeuristicNaive & rhs) |
---|
[1286] | 1234 | : |
---|
| 1235 | CbcHeuristic(rhs), |
---|
| 1236 | large_(rhs.large_) |
---|
[1100] | 1237 | { |
---|
| 1238 | } |
---|
| 1239 | |
---|
[1286] | 1240 | // Assignment operator |
---|
| 1241 | CbcHeuristicNaive & |
---|
| 1242 | CbcHeuristicNaive::operator=( const CbcHeuristicNaive & rhs) |
---|
[1100] | 1243 | { |
---|
[1286] | 1244 | if (this != &rhs) { |
---|
| 1245 | CbcHeuristic::operator=(rhs); |
---|
| 1246 | large_ = rhs.large_; |
---|
| 1247 | } |
---|
| 1248 | return *this; |
---|
[1100] | 1249 | } |
---|
| 1250 | |
---|
| 1251 | // Resets stuff if model changes |
---|
[1286] | 1252 | void |
---|
[1100] | 1253 | CbcHeuristicNaive::resetModel(CbcModel * model) |
---|
| 1254 | { |
---|
[1286] | 1255 | CbcHeuristic::resetModel(model); |
---|
[1100] | 1256 | } |
---|
| 1257 | int |
---|
| 1258 | CbcHeuristicNaive::solution(double & solutionValue, |
---|
[1286] | 1259 | double * betterSolution) |
---|
[1100] | 1260 | { |
---|
[1286] | 1261 | numCouldRun_++; |
---|
| 1262 | // See if to do |
---|
| 1263 | bool atRoot = model_->getNodeCount() == 0; |
---|
| 1264 | int passNumber = model_->getCurrentPassNumber(); |
---|
| 1265 | if (!when() || (when() == 1 && model_->phase() != 1) || !atRoot || passNumber != 1) |
---|
| 1266 | return 0; // switched off |
---|
| 1267 | // Don't do if it was this heuristic which found solution! |
---|
| 1268 | if (this == model_->lastHeuristic()) |
---|
| 1269 | return 0; |
---|
| 1270 | numRuns_++; |
---|
| 1271 | double cutoff; |
---|
| 1272 | model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff); |
---|
| 1273 | double direction = model_->solver()->getObjSense(); |
---|
| 1274 | cutoff *= direction; |
---|
| 1275 | cutoff = CoinMin(cutoff, solutionValue); |
---|
| 1276 | OsiSolverInterface * solver = model_->continuousSolver(); |
---|
| 1277 | if (!solver) |
---|
| 1278 | solver = model_->solver(); |
---|
| 1279 | const double * colLower = solver->getColLower(); |
---|
| 1280 | const double * colUpper = solver->getColUpper(); |
---|
| 1281 | const double * objective = solver->getObjCoefficients(); |
---|
[1100] | 1282 | |
---|
[1286] | 1283 | int numberColumns = model_->getNumCols(); |
---|
| 1284 | int numberIntegers = model_->numberIntegers(); |
---|
| 1285 | const int * integerVariable = model_->integerVariable(); |
---|
| 1286 | |
---|
| 1287 | int i; |
---|
| 1288 | bool solutionFound = false; |
---|
| 1289 | CoinWarmStartBasis saveBasis; |
---|
| 1290 | CoinWarmStartBasis * basis = |
---|
| 1291 | dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ; |
---|
| 1292 | if (basis) { |
---|
| 1293 | saveBasis = * basis; |
---|
| 1294 | delete basis; |
---|
[1100] | 1295 | } |
---|
[1286] | 1296 | // First just fix all integers as close to zero as possible |
---|
| 1297 | OsiSolverInterface * newSolver = cloneBut(7); // wassolver->clone(); |
---|
| 1298 | for (i = 0; i < numberIntegers; i++) { |
---|
| 1299 | int iColumn = integerVariable[i]; |
---|
| 1300 | double lower = colLower[iColumn]; |
---|
| 1301 | double upper = colUpper[iColumn]; |
---|
| 1302 | double value; |
---|
| 1303 | if (lower > 0.0) |
---|
| 1304 | value = lower; |
---|
| 1305 | else if (upper < 0.0) |
---|
| 1306 | value = upper; |
---|
| 1307 | else |
---|
| 1308 | value = 0.0; |
---|
| 1309 | newSolver->setColLower(iColumn, value); |
---|
| 1310 | newSolver->setColUpper(iColumn, value); |
---|
[1100] | 1311 | } |
---|
| 1312 | newSolver->initialSolve(); |
---|
| 1313 | if (newSolver->isProvenOptimal()) { |
---|
[1286] | 1314 | double solValue = newSolver->getObjValue() * direction ; |
---|
| 1315 | if (solValue < cutoff) { |
---|
| 1316 | // we have a solution |
---|
| 1317 | solutionFound = true; |
---|
| 1318 | solutionValue = solValue; |
---|
| 1319 | memcpy(betterSolution, newSolver->getColSolution(), |
---|
| 1320 | numberColumns*sizeof(double)); |
---|
[1641] | 1321 | COIN_DETAIL_PRINT(printf("Naive fixing close to zero gave solution of %g\n", solutionValue)); |
---|
[1286] | 1322 | cutoff = solValue - model_->getCutoffIncrement(); |
---|
| 1323 | } |
---|
[1100] | 1324 | } |
---|
[1286] | 1325 | // Now fix all integers as close to zero if zero or large cost |
---|
| 1326 | int nFix = 0; |
---|
| 1327 | for (i = 0; i < numberIntegers; i++) { |
---|
| 1328 | int iColumn = integerVariable[i]; |
---|
| 1329 | double lower = colLower[iColumn]; |
---|
| 1330 | double upper = colUpper[iColumn]; |
---|
| 1331 | double value; |
---|
| 1332 | if (fabs(objective[i]) > 0.0 && fabs(objective[i]) < large_) { |
---|
| 1333 | nFix++; |
---|
| 1334 | if (lower > 0.0) |
---|
| 1335 | value = lower; |
---|
| 1336 | else if (upper < 0.0) |
---|
| 1337 | value = upper; |
---|
| 1338 | else |
---|
| 1339 | value = 0.0; |
---|
| 1340 | newSolver->setColLower(iColumn, value); |
---|
| 1341 | newSolver->setColUpper(iColumn, value); |
---|
| 1342 | } else { |
---|
| 1343 | // set back to original |
---|
| 1344 | newSolver->setColLower(iColumn, lower); |
---|
| 1345 | newSolver->setColUpper(iColumn, upper); |
---|
| 1346 | } |
---|
| 1347 | } |
---|
| 1348 | const double * solution = solver->getColSolution(); |
---|
| 1349 | if (nFix) { |
---|
| 1350 | newSolver->setWarmStart(&saveBasis); |
---|
| 1351 | newSolver->setColSolution(solution); |
---|
| 1352 | newSolver->initialSolve(); |
---|
| 1353 | if (newSolver->isProvenOptimal()) { |
---|
| 1354 | double solValue = newSolver->getObjValue() * direction ; |
---|
| 1355 | if (solValue < cutoff) { |
---|
| 1356 | // try branch and bound |
---|
| 1357 | double * newSolution = new double [numberColumns]; |
---|
[1641] | 1358 | COIN_DETAIL_PRINT(printf("%d fixed after fixing costs\n", nFix)); |
---|
[1286] | 1359 | int returnCode = smallBranchAndBound(newSolver, |
---|
| 1360 | numberNodes_, newSolution, |
---|
| 1361 | solutionValue, |
---|
| 1362 | solutionValue, "CbcHeuristicNaive1"); |
---|
| 1363 | if (returnCode < 0) |
---|
| 1364 | returnCode = 0; // returned on size |
---|
| 1365 | if ((returnCode&2) != 0) { |
---|
| 1366 | // could add cut |
---|
| 1367 | returnCode &= ~2; |
---|
| 1368 | } |
---|
| 1369 | if (returnCode == 1) { |
---|
| 1370 | // solution |
---|
| 1371 | solutionFound = true; |
---|
| 1372 | memcpy(betterSolution, newSolution, |
---|
| 1373 | numberColumns*sizeof(double)); |
---|
[1641] | 1374 | COIN_DETAIL_PRINT(printf("Naive fixing zeros gave solution of %g\n", solutionValue)); |
---|
[1286] | 1375 | cutoff = solutionValue - model_->getCutoffIncrement(); |
---|
| 1376 | } |
---|
| 1377 | delete [] newSolution; |
---|
| 1378 | } |
---|
| 1379 | } |
---|
| 1380 | } |
---|
[1100] | 1381 | #if 1 |
---|
[1286] | 1382 | newSolver->setObjSense(-direction); // maximize |
---|
| 1383 | newSolver->setWarmStart(&saveBasis); |
---|
| 1384 | newSolver->setColSolution(solution); |
---|
| 1385 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 1386 | double value = solution[iColumn]; |
---|
| 1387 | double lower = colLower[iColumn]; |
---|
| 1388 | double upper = colUpper[iColumn]; |
---|
| 1389 | double newLower; |
---|
| 1390 | double newUpper; |
---|
| 1391 | if (newSolver->isInteger(iColumn)) { |
---|
| 1392 | newLower = CoinMax(lower, floor(value) - 2.0); |
---|
| 1393 | newUpper = CoinMin(upper, ceil(value) + 2.0); |
---|
| 1394 | } else { |
---|
| 1395 | newLower = CoinMax(lower, value - 1.0e5); |
---|
| 1396 | newUpper = CoinMin(upper, value + 1.0e-5); |
---|
| 1397 | } |
---|
| 1398 | newSolver->setColLower(iColumn, newLower); |
---|
| 1399 | newSolver->setColUpper(iColumn, newUpper); |
---|
[1100] | 1400 | } |
---|
[1286] | 1401 | newSolver->initialSolve(); |
---|
| 1402 | if (newSolver->isProvenOptimal()) { |
---|
| 1403 | double solValue = newSolver->getObjValue() * direction ; |
---|
| 1404 | if (solValue < cutoff) { |
---|
| 1405 | nFix = 0; |
---|
| 1406 | newSolver->setObjSense(direction); // correct direction |
---|
| 1407 | //const double * thisSolution = newSolver->getColSolution(); |
---|
| 1408 | for (int iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
| 1409 | double value = solution[iColumn]; |
---|
| 1410 | double lower = colLower[iColumn]; |
---|
| 1411 | double upper = colUpper[iColumn]; |
---|
| 1412 | double newLower = lower; |
---|
| 1413 | double newUpper = upper; |
---|
| 1414 | if (newSolver->isInteger(iColumn)) { |
---|
| 1415 | if (value < lower + 1.0e-6) { |
---|
| 1416 | nFix++; |
---|
| 1417 | newUpper = lower; |
---|
| 1418 | } else if (value > upper - 1.0e-6) { |
---|
| 1419 | nFix++; |
---|
| 1420 | newLower = upper; |
---|
| 1421 | } else { |
---|
| 1422 | newLower = CoinMax(lower, floor(value) - 2.0); |
---|
| 1423 | newUpper = CoinMin(upper, ceil(value) + 2.0); |
---|
| 1424 | } |
---|
| 1425 | } |
---|
| 1426 | newSolver->setColLower(iColumn, newLower); |
---|
| 1427 | newSolver->setColUpper(iColumn, newUpper); |
---|
| 1428 | } |
---|
| 1429 | // try branch and bound |
---|
| 1430 | double * newSolution = new double [numberColumns]; |
---|
[1641] | 1431 | COIN_DETAIL_PRINT(printf("%d fixed after maximizing\n", nFix)); |
---|
[1286] | 1432 | int returnCode = smallBranchAndBound(newSolver, |
---|
| 1433 | numberNodes_, newSolution, |
---|
| 1434 | solutionValue, |
---|
| 1435 | solutionValue, "CbcHeuristicNaive1"); |
---|
| 1436 | if (returnCode < 0) |
---|
| 1437 | returnCode = 0; // returned on size |
---|
| 1438 | if ((returnCode&2) != 0) { |
---|
| 1439 | // could add cut |
---|
| 1440 | returnCode &= ~2; |
---|
| 1441 | } |
---|
| 1442 | if (returnCode == 1) { |
---|
| 1443 | // solution |
---|
| 1444 | solutionFound = true; |
---|
| 1445 | memcpy(betterSolution, newSolution, |
---|
| 1446 | numberColumns*sizeof(double)); |
---|
[1641] | 1447 | COIN_DETAIL_PRINT(printf("Naive maximizing gave solution of %g\n", solutionValue)); |
---|
[1286] | 1448 | cutoff = solutionValue - model_->getCutoffIncrement(); |
---|
| 1449 | } |
---|
| 1450 | delete [] newSolution; |
---|
| 1451 | } |
---|
[1100] | 1452 | } |
---|
| 1453 | #endif |
---|
[1286] | 1454 | delete newSolver; |
---|
| 1455 | return solutionFound ? 1 : 0; |
---|
[1100] | 1456 | } |
---|
| 1457 | // update model |
---|
| 1458 | void CbcHeuristicNaive::setModel(CbcModel * model) |
---|
| 1459 | { |
---|
[1286] | 1460 | model_ = model; |
---|
[1100] | 1461 | } |
---|
[1271] | 1462 | // Default Constructor |
---|
[1286] | 1463 | CbcHeuristicCrossover::CbcHeuristicCrossover() |
---|
| 1464 | : CbcHeuristic(), |
---|
| 1465 | numberSolutions_(0), |
---|
| 1466 | useNumber_(3) |
---|
[1271] | 1467 | { |
---|
[1286] | 1468 | setWhen(1); |
---|
[1271] | 1469 | } |
---|
[1100] | 1470 | |
---|
[1271] | 1471 | // Constructor with model - assumed before cuts |
---|
| 1472 | |
---|
| 1473 | CbcHeuristicCrossover::CbcHeuristicCrossover(CbcModel & model) |
---|
[1286] | 1474 | : CbcHeuristic(model), |
---|
| 1475 | numberSolutions_(0), |
---|
| 1476 | useNumber_(3) |
---|
[1271] | 1477 | { |
---|
[1286] | 1478 | setWhen(1); |
---|
| 1479 | for (int i = 0; i < 10; i++) |
---|
| 1480 | random_[i] = model.randomNumberGenerator()->randomDouble(); |
---|
[1271] | 1481 | } |
---|
| 1482 | |
---|
[1286] | 1483 | // Destructor |
---|
[1271] | 1484 | CbcHeuristicCrossover::~CbcHeuristicCrossover () |
---|
| 1485 | { |
---|
| 1486 | } |
---|
| 1487 | |
---|
| 1488 | // Clone |
---|
| 1489 | CbcHeuristic * |
---|
| 1490 | CbcHeuristicCrossover::clone() const |
---|
| 1491 | { |
---|
[1286] | 1492 | return new CbcHeuristicCrossover(*this); |
---|
[1271] | 1493 | } |
---|
| 1494 | // Create C++ lines to get to current state |
---|
[1286] | 1495 | void |
---|
| 1496 | CbcHeuristicCrossover::generateCpp( FILE * fp) |
---|
[1271] | 1497 | { |
---|
[1286] | 1498 | CbcHeuristicCrossover other; |
---|
[1802] | 1499 | fprintf(fp, "0#include \"CbcHeuristicProximity.hpp\"\n"); |
---|
[1286] | 1500 | fprintf(fp, "3 CbcHeuristicCrossover crossover(*cbcModel);\n"); |
---|
| 1501 | CbcHeuristic::generateCpp(fp, "crossover"); |
---|
| 1502 | if (useNumber_ != other.useNumber_) |
---|
| 1503 | fprintf(fp, "3 crossover.setNumberSolutions(%d);\n", useNumber_); |
---|
| 1504 | else |
---|
| 1505 | fprintf(fp, "4 crossover.setNumberSolutions(%d);\n", useNumber_); |
---|
| 1506 | fprintf(fp, "3 cbcModel->addHeuristic(&crossover);\n"); |
---|
[1271] | 1507 | } |
---|
| 1508 | |
---|
[1286] | 1509 | // Copy constructor |
---|
[1271] | 1510 | CbcHeuristicCrossover::CbcHeuristicCrossover(const CbcHeuristicCrossover & rhs) |
---|
[1286] | 1511 | : |
---|
| 1512 | CbcHeuristic(rhs), |
---|
| 1513 | attempts_(rhs.attempts_), |
---|
| 1514 | numberSolutions_(rhs.numberSolutions_), |
---|
| 1515 | useNumber_(rhs.useNumber_) |
---|
[1271] | 1516 | { |
---|
[1286] | 1517 | memcpy(random_, rhs.random_, 10*sizeof(double)); |
---|
[1271] | 1518 | } |
---|
| 1519 | |
---|
[1286] | 1520 | // Assignment operator |
---|
| 1521 | CbcHeuristicCrossover & |
---|
| 1522 | CbcHeuristicCrossover::operator=( const CbcHeuristicCrossover & rhs) |
---|
[1271] | 1523 | { |
---|
[1286] | 1524 | if (this != &rhs) { |
---|
| 1525 | CbcHeuristic::operator=(rhs); |
---|
| 1526 | useNumber_ = rhs.useNumber_; |
---|
| 1527 | attempts_ = rhs.attempts_; |
---|
| 1528 | numberSolutions_ = rhs.numberSolutions_; |
---|
| 1529 | memcpy(random_, rhs.random_, 10*sizeof(double)); |
---|
| 1530 | } |
---|
| 1531 | return *this; |
---|
[1271] | 1532 | } |
---|
| 1533 | |
---|
| 1534 | // Resets stuff if model changes |
---|
[1286] | 1535 | void |
---|
[1271] | 1536 | CbcHeuristicCrossover::resetModel(CbcModel * model) |
---|
| 1537 | { |
---|
[1286] | 1538 | CbcHeuristic::resetModel(model); |
---|
[1271] | 1539 | } |
---|
| 1540 | int |
---|
| 1541 | CbcHeuristicCrossover::solution(double & solutionValue, |
---|
[1286] | 1542 | double * betterSolution) |
---|
[1271] | 1543 | { |
---|
[1286] | 1544 | if (when_ == 0) |
---|
| 1545 | return 0; |
---|
| 1546 | numCouldRun_++; |
---|
| 1547 | bool useBest = (numberSolutions_ != model_->getSolutionCount()); |
---|
| 1548 | if (!useBest && (when_ % 10) == 1) |
---|
| 1549 | return 0; |
---|
| 1550 | numberSolutions_ = model_->getSolutionCount(); |
---|
| 1551 | OsiSolverInterface * continuousSolver = model_->continuousSolver(); |
---|
| 1552 | int useNumber = CoinMin(model_->numberSavedSolutions(), useNumber_); |
---|
| 1553 | if (useNumber < 2 || !continuousSolver) |
---|
| 1554 | return 0; |
---|
| 1555 | // Fix later |
---|
| 1556 | if (!useBest) |
---|
| 1557 | abort(); |
---|
| 1558 | numRuns_++; |
---|
| 1559 | double cutoff; |
---|
| 1560 | model_->solver()->getDblParam(OsiDualObjectiveLimit, cutoff); |
---|
| 1561 | double direction = model_->solver()->getObjSense(); |
---|
| 1562 | cutoff *= direction; |
---|
| 1563 | cutoff = CoinMin(cutoff, solutionValue); |
---|
| 1564 | OsiSolverInterface * solver = cloneBut(2); |
---|
| 1565 | // But reset bounds |
---|
| 1566 | solver->setColLower(continuousSolver->getColLower()); |
---|
| 1567 | solver->setColUpper(continuousSolver->getColUpper()); |
---|
| 1568 | int numberColumns = solver->getNumCols(); |
---|
| 1569 | // Fixed |
---|
| 1570 | double * fixed = new double [numberColumns]; |
---|
| 1571 | for (int i = 0; i < numberColumns; i++) |
---|
| 1572 | fixed[i] = -COIN_DBL_MAX; |
---|
| 1573 | int whichSolution[10]; |
---|
| 1574 | for (int i = 0; i < useNumber; i++) |
---|
| 1575 | whichSolution[i] = i; |
---|
| 1576 | for (int i = 0; i < useNumber; i++) { |
---|
| 1577 | int k = whichSolution[i]; |
---|
| 1578 | const double * solution = model_->savedSolution(k); |
---|
| 1579 | for (int j = 0; j < numberColumns; j++) { |
---|
| 1580 | if (solver->isInteger(j)) { |
---|
| 1581 | if (fixed[j] == -COIN_DBL_MAX) |
---|
| 1582 | fixed[j] = floor(solution[j] + 0.5); |
---|
| 1583 | else if (fabs(fixed[j] - solution[j]) > 1.0e-7) |
---|
| 1584 | fixed[j] = COIN_DBL_MAX; |
---|
| 1585 | } |
---|
| 1586 | } |
---|
[1271] | 1587 | } |
---|
[1286] | 1588 | const double * colLower = solver->getColLower(); |
---|
| 1589 | for (int i = 0; i < numberColumns; i++) { |
---|
| 1590 | if (solver->isInteger(i)) { |
---|
| 1591 | double value = fixed[i]; |
---|
| 1592 | if (value != COIN_DBL_MAX) { |
---|
| 1593 | if (when_ < 10) { |
---|
| 1594 | solver->setColLower(i, value); |
---|
| 1595 | solver->setColUpper(i, value); |
---|
| 1596 | } else if (value == colLower[i]) { |
---|
| 1597 | solver->setColUpper(i, value); |
---|
| 1598 | } |
---|
| 1599 | } |
---|
| 1600 | } |
---|
[1271] | 1601 | } |
---|
[1286] | 1602 | int returnCode = smallBranchAndBound(solver, numberNodes_, betterSolution, |
---|
| 1603 | solutionValue, |
---|
| 1604 | solutionValue, "CbcHeuristicCrossover"); |
---|
| 1605 | if (returnCode < 0) |
---|
| 1606 | returnCode = 0; // returned on size |
---|
| 1607 | if ((returnCode&2) != 0) { |
---|
| 1608 | // could add cut |
---|
| 1609 | returnCode &= ~2; |
---|
| 1610 | } |
---|
[1271] | 1611 | |
---|
[1286] | 1612 | delete solver; |
---|
| 1613 | return returnCode; |
---|
[1271] | 1614 | } |
---|
| 1615 | // update model |
---|
| 1616 | void CbcHeuristicCrossover::setModel(CbcModel * model) |
---|
| 1617 | { |
---|
[1286] | 1618 | model_ = model; |
---|
| 1619 | if (model) { |
---|
| 1620 | for (int i = 0; i < 10; i++) |
---|
| 1621 | random_[i] = model->randomNumberGenerator()->randomDouble(); |
---|
| 1622 | } |
---|
[1271] | 1623 | } |
---|
| 1624 | |
---|
[1286] | 1625 | |
---|