source: stable/2.4/Cbc/src/CbcHeuristicFPump.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

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