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