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