source: trunk/Cbc/src/CbcHeuristicFPump.cpp

Last change on this file was 2467, checked in by unxusr, 6 months ago

spaces after angles

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 126.1 KB
Line 
1/* $Id: CbcHeuristicFPump.cpp 2467 2019-01-03 21:26:29Z unxusr $ */
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
33CbcHeuristicFPump::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
56CbcHeuristicFPump::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
80CbcHeuristicFPump::~CbcHeuristicFPump()
81{
82}
83
84// Clone
85CbcHeuristic *
86CbcHeuristicFPump::clone() const
87{
88  return new CbcHeuristicFPump(*this);
89}
90// Create C++ lines to get to current state
91void 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
157CbcHeuristicFPump::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
179CbcHeuristicFPump &
180CbcHeuristicFPump::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
205void 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
214int 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*/
2450int 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
2589void 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*/
2597int 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
3212void 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//-------------------------------------------------------------------
3227CbcDisasterHandler::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//-------------------------------------------------------------------
3242CbcDisasterHandler::CbcDisasterHandler(const CbcDisasterHandler &rhs)
3243  : OsiClpDisasterHandler(rhs)
3244  , cbcModel_(rhs.cbcModel_)
3245{
3246}
3247
3248//-------------------------------------------------------------------
3249// Destructor
3250//-------------------------------------------------------------------
3251CbcDisasterHandler::~CbcDisasterHandler()
3252{
3253}
3254
3255//----------------------------------------------------------------
3256// Assignment operator
3257//-------------------------------------------------------------------
3258CbcDisasterHandler &
3259CbcDisasterHandler::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//-------------------------------------------------------------------
3270ClpDisasterHandler *CbcDisasterHandler::clone() const
3271{
3272  return new CbcDisasterHandler(*this);
3273}
3274// Type of disaster 0 can fix, 1 abort
3275int 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. */
3286void 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*/
Note: See TracBrowser for help on using the repository browser.