source: trunk/Cbc/src/CbcHeuristicFPump.cpp @ 1129

Last change on this file since 1129 was 1129, checked in by forrest, 11 years ago

slight modification for rounding heuristic

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