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

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

chnages to try and make faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 82.2 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      bool exitThis=false;
579      if (iterationLimit<0.0) {
580        if (numberPasses>=maximumPasses_) {
581          // If going well then keep going if maximumPasses_ small
582          if (lastMove<numberPasses-4||lastMove==1000000)
583            exitThis=true;
584          if (maximumPasses_>20||numberPasses>=40)
585            exitThis=true;
586        }
587      } 
588      if (iterationLimit>0.0&&totalNumberIterations>iterationLimit
589          &&numberPasses>15) {
590          // exiting on iteration count
591        exitAll=true;
592      } else if (maximumPasses_<30&&numberPasses>100) {
593        // too many passes anyway
594        exitAll=true;
595      }
596      if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) 
597        exitAll=true;
598      if (exitAll||exitThis)
599        break;
600      memcpy(newSolution,solution,numberColumns*sizeof(double));
601      int flip;
602      if (numberPasses==0&&false) {
603        // always use same seed
604        randomNumberGenerator_.setSeed(987654321);
605      }
606      returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
607                          pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
608      if (numberPasses==0&&false) {
609        // Make sure random will be different
610        for (i=1;i<numberTries;i++)
611          randomNumberGenerator_.randomDouble();
612      }
613      numberPasses++;
614      if (returnCode) {
615        // SOLUTION IS INTEGER
616        // Put back correct objective
617        for (i=0;i<numberColumns;i++)
618          solver->setObjCoeff(i,saveObjective[i]);
619
620        // solution - but may not be better
621        // Compute using dot product
622        solver->setDblParam(OsiObjOffset,saveOffset);
623        newSolutionValue = -saveOffset;
624        for (  i=0 ; i<numberColumns ; i++ )
625          newSolutionValue += saveObjective[i]*newSolution[i];
626        newSolutionValue *= direction;
627        sprintf(pumpPrint,"Solution found of %g",newSolutionValue);
628        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
629          << pumpPrint
630          <<CoinMessageEol;
631        //newLineNeeded=false;
632        if (newSolutionValue<solutionValue) {
633          double saveValue = solutionValue;
634          if (!doGeneral) {
635            int numberLeft=0;
636            for (i=0;i<numberIntegersOrig;i++) {
637              int iColumn = integerVariableOrig[i];
638              double value = floor(newSolution[iColumn]+0.5);
639              if(solver->isBinary(iColumn)) {
640                solver->setColLower(iColumn,value);
641                solver->setColUpper(iColumn,value);
642              } else {
643                if (fabs(value-newSolution[iColumn])>1.0e-7) 
644                  numberLeft++;
645              }
646            }
647            if (numberLeft) {
648              sprintf(pumpPrint,"Branch and bound needed to clear up %d general integers",numberLeft);
649              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
650                << pumpPrint
651                <<CoinMessageEol;
652              returnCode = smallBranchAndBound(solver,numberNodes_,newSolution,newSolutionValue,
653                                               solutionValue,"CbcHeuristicFpump");
654              if (returnCode<0) {
655                if (returnCode==-2)
656                  exitAll=true;
657                returnCode=0; // returned on size or event
658              }
659              if ((returnCode&2)!=0) {
660                // could add cut
661                returnCode &= ~2;
662              }
663              if (returnCode!=1)
664                newSolutionValue=saveValue;
665              if (returnCode&&newSolutionValue<saveValue) 
666                numberBandBsolutions++;
667            }
668          }
669          if (returnCode&&newSolutionValue<saveValue) {
670            memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
671            solutionFound=true;
672            if (exitNow(newSolutionValue))
673              exitAll=true;
674            CoinWarmStartBasis * basis =
675              dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
676            if (basis) {
677              bestBasis = * basis;
678              delete basis;
679              CbcEventHandler * handler = model_->getEventHandler();
680              if (handler) {
681                double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
682                double saveObjectiveValue = model_->getMinimizationObjValue();
683                model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
684                int action = handler->event(CbcEventHandler::heuristicSolution);
685                if (saveOldSolution&&saveObjectiveValue<model_->getMinimizationObjValue())
686                  model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
687                delete [] saveOldSolution;
688                if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
689                  exitAll=true; // exit
690                  break;
691                }
692              }
693            }
694            if ((accumulate_&1)!=0) {
695              model_->incrementUsed(betterSolution); // for local search
696              numberSolutions++;
697            }
698            solutionValue=newSolutionValue;
699            solutionFound=true;
700            if (general&&saveValue!=newSolutionValue) {
701              sprintf(pumpPrint,"Cleaned solution of %g",solutionValue);
702              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
703                << pumpPrint
704                <<CoinMessageEol;
705            }
706            if (exitNow(newSolutionValue))
707              exitAll=true;
708          } else {
709            sprintf(pumpPrint,"Mini branch and bound could not fix general integers");
710            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
711              << pumpPrint
712              <<CoinMessageEol;
713          }
714        } else {
715          sprintf(pumpPrint,"After further testing solution no better than previous of %g",solutionValue);
716          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
717            << pumpPrint
718            <<CoinMessageEol;
719          //newLineNeeded=false;
720          returnCode=0;
721        }     
722        break;
723      } else {
724        // SOLUTION IS not INTEGER
725        // 1. check for loop
726        bool matched;
727        for (int k = NUMBER_OLD-1; k > 0; k--) {
728          double * b = oldSolution[k];
729          matched = true;
730          for (i = 0; i <numberIntegers; i++) {
731            int iColumn = integerVariable[i];
732            if (newSolution[iColumn]!=b[iColumn]) {
733              matched=false;
734              break;
735            }
736          }
737          if (matched) break;
738        }
739        int numberPerturbed=0;
740        if (matched || numberPasses%100 == 0) {
741          // perturbation
742          //sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
743          //newLineNeeded=true;
744          double factorX[10]={0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0};
745          double factor=1.0;
746          double target=-1.0;
747          double * randomX = new double [numberIntegers];
748          for (i=0;i<numberIntegers;i++) 
749            randomX[i] = CoinMax(0.0,randomNumberGenerator_.randomDouble()-0.3);
750          for (int k=0;k<10;k++) {
751#ifdef COIN_DEVELOP_x
752            printf("kpass %d\n",k);
753#endif
754            int numberX[10]={0,0,0,0,0,0,0,0,0,0};
755            for (i=0;i<numberIntegers;i++) {
756              int iColumn = integerVariable[i];
757              double value = randomX[i];
758              double difference = fabs(solution[iColumn]-newSolution[iColumn]);
759              for (int j=0;j<10;j++) {
760                if (difference+value*factorX[j]>0.5) 
761                  numberX[j]++;
762              }
763            }
764            if (target<0.0) {
765              if (numberX[9]<=200)
766                break; // not very many changes
767              target=CoinMax(200.0,CoinMin(0.05*numberX[9],1000.0));
768            }
769            int iX=-1;
770            int iBand=-1;
771            for (i=0;i<10;i++) {
772#ifdef COIN_DEVELOP_x
773              printf("** %d changed at %g\n",numberX[i],factorX[i]);
774#endif
775              if (numberX[i]>=target&&numberX[i]<2.0*target&&iX<0)
776                iX=i;
777              if (iBand<0&&numberX[i]>target) {
778                iBand=i;
779                factor=factorX[i];
780              }
781            }
782            if (iX>=0) {
783              factor=factorX[iX];
784              break;
785            } else {
786              assert (iBand>=0);
787              double hi = factor;
788              double lo = (iBand>0) ? factorX[iBand-1] : 0.0;
789              double diff = (hi-lo)/9.0;
790              for (i=0;i<10;i++) {
791                factorX[i]=lo;
792                lo += diff;
793              }
794            }
795          }
796          for (i=0;i<numberIntegers;i++) {
797            int iColumn = integerVariable[i];
798            double value = randomX[i];
799            double difference = fabs(solution[iColumn]-newSolution[iColumn]);
800            if (difference+value*factor>0.5) {
801              numberPerturbed++;
802              if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
803                newSolution[iColumn] += 1.0;
804              } else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
805                newSolution[iColumn] -= 1.0;
806              } else {
807                // general integer
808                if (difference+value>0.75)
809                  newSolution[iColumn] += 1.0;
810                else
811                  newSolution[iColumn] -= 1.0;
812              }
813            }
814          }
815          delete [] randomX;
816        } else {
817          for (j=NUMBER_OLD-1;j>0;j--) {
818            for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
819          }
820          for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
821        }
822       
823        // 2. update the objective function based on the new rounded solution
824        double offset=0.0;
825        double costValue = (1.0-scaleFactor)*solver->getObjSense();
826        int numberChanged=0;
827        const double * oldObjective = solver->getObjCoefficients();
828        for (i=0;i<numberColumns;i++) {
829          // below so we can keep original code and allow for objective
830          int iColumn = i;
831          // Special code for "artificials"
832          if (direction*saveObjective[iColumn]>=artificialCost_) {
833            //solver->setObjCoeff(iColumn,scaleFactor*saveObjective[iColumn]);
834            solver->setObjCoeff(iColumn,(artificialFactor*saveObjective[iColumn])/artificialCost_);
835          }
836          if(!solver->isBinary(iColumn)&&!doGeneral)
837            continue;
838          // deal with fixed variables (i.e., upper=lower)
839          if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance||!solver->isInteger(iColumn)) {
840            //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
841            //else                                       solver->setObjCoeff(iColumn,costValue);
842            continue;
843          }
844          double newValue=0.0;
845          if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
846            newValue = costValue+scaleFactor*saveObjective[iColumn];
847          } else {
848            if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
849              newValue = -costValue+scaleFactor*saveObjective[iColumn];
850            }
851          }
852          if (newValue!=oldObjective[iColumn])
853            numberChanged++;
854          solver->setObjCoeff(iColumn,newValue);
855          offset += costValue*newSolution[iColumn];
856        }
857        solver->setDblParam(OsiObjOffset,-offset);
858        if (!general&&false) {
859          // Solve in two goes - first keep satisfied ones fixed
860          double * saveLower = new double [numberIntegers];
861          double * saveUpper = new double [numberIntegers];
862          for (i=0;i<numberIntegers;i++) {
863            int iColumn = integerVariable[i];
864            saveLower[i]=COIN_DBL_MAX;
865            saveUpper[i]=-COIN_DBL_MAX;
866            if (solution[iColumn]<lower[iColumn]+primalTolerance) {
867              saveUpper[i]=upper[iColumn];
868              solver->setColUpper(iColumn,lower[iColumn]);
869            } else if (solution[iColumn]>upper[iColumn]-primalTolerance) {
870              saveLower[i]=lower[iColumn];
871              solver->setColLower(iColumn,upper[iColumn]);
872            }
873          }
874          solver->resolve();
875          if (!solver->isProvenOptimal()) {
876            // presumably max time or some such
877            exitAll=true;
878            break;
879          }
880          for (i=0;i<numberIntegers;i++) {
881            int iColumn = integerVariable[i];
882            if (saveLower[i]!=COIN_DBL_MAX)
883              solver->setColLower(iColumn,saveLower[i]);
884            if (saveUpper[i]!=-COIN_DBL_MAX)
885              solver->setColUpper(iColumn,saveUpper[i]);
886            saveUpper[i]=-COIN_DBL_MAX;
887          }
888          memcpy(newSolution,solution,numberColumns*sizeof(double));
889          int flip;
890          returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
891                              pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
892          numberPasses++;
893          if (returnCode) {
894            // solution - but may not be better
895            // Compute using dot product
896            double newSolutionValue = -saveOffset;
897            for (  i=0 ; i<numberColumns ; i++ )
898              newSolutionValue += saveObjective[i]*newSolution[i];
899            newSolutionValue *= direction;
900            sprintf(pumpPrint,"Intermediate solution found of %g",newSolutionValue);
901            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
902              << pumpPrint
903              <<CoinMessageEol;
904            if (newSolutionValue<solutionValue) {
905              memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
906              CoinWarmStartBasis * basis =
907                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
908              solutionFound=true;
909              if (exitNow(newSolutionValue))
910                exitAll=true;
911              if (basis) {
912                bestBasis = * basis;
913                delete basis;
914                CbcEventHandler * handler = model_->getEventHandler();
915                if (handler) {
916                  double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
917                  double saveObjectiveValue = model_->getMinimizationObjValue();
918                  model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
919                  int action = handler->event(CbcEventHandler::heuristicSolution);
920                  if (saveOldSolution&&saveObjectiveValue<model_->getMinimizationObjValue())
921                    model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
922                  delete [] saveOldSolution;
923                  if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
924                    exitAll=true; // exit
925                    break;
926                  }
927                }
928              }
929              if ((accumulate_&1)!=0) {
930                model_->incrementUsed(betterSolution); // for local search
931                numberSolutions++;
932              }
933              solutionValue=newSolutionValue;
934              solutionFound=true;
935              if (exitNow(newSolutionValue))
936                exitAll=true;
937            } else {
938              returnCode=0;
939            }
940          }     
941        }
942        int numberIterations=0;
943        if (!doGeneral) {
944          // faster to do from all slack!!!!
945          if (allSlack) {
946            CoinWarmStartBasis dummy;
947            solver->setWarmStart(&dummy);
948          }
949#ifdef COIN_DEVELOP
950          printf("%d perturbed out of %d columns (%d changed)\n",numberPerturbed,numberColumns,numberChanged);
951#endif
952          bool takeHint;
953          OsiHintStrength strength;
954          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
955          if (dualPass) {
956            solver->setHintParam(OsiDoDualInResolve,true); // dual may be better
957            if (dualPass==1) {
958              // but we need to make infeasible
959              CoinWarmStartBasis * basis =
960                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
961              if (basis) {
962                // modify
963                const double * lower = solver->getColLower();
964                const double * upper = solver->getColUpper();
965                double * solution = CoinCopyOfArray(solver->getColSolution(),
966                                                    numberColumns);
967                const double * objective = solver->getObjCoefficients();
968                int nChanged=0;
969                for (i=0;i<numberIntegersOrig;i++) {
970                  int iColumn=integerVariableOrig[i];
971                  if (objective[iColumn]>0.0) {
972                    if (basis->getStructStatus(iColumn)==
973                        CoinWarmStartBasis::atUpperBound) {
974                      solution[iColumn]=lower[iColumn];
975                      basis->setStructStatus(iColumn,CoinWarmStartBasis::atLowerBound);
976                      nChanged++;
977                    }
978                  } else if (objective[iColumn]<0.0) {
979                    if (basis->getStructStatus(iColumn)==
980                        CoinWarmStartBasis::atUpperBound) {
981                      solution[iColumn]=upper[iColumn];
982                      basis->setStructStatus(iColumn,CoinWarmStartBasis::atUpperBound);
983                      nChanged++;
984                    }
985                  }
986                }
987                if (!nChanged) {
988                  for (i=0;i<numberIntegersOrig;i++) {
989                    int iColumn=integerVariableOrig[i];
990                    if (objective[iColumn]>0.0) {
991                      if (basis->getStructStatus(iColumn)==
992                          CoinWarmStartBasis::basic) {
993                        solution[iColumn]=lower[iColumn];
994                        basis->setStructStatus(iColumn,CoinWarmStartBasis::atLowerBound);
995                        break;
996                      }
997                    } else if (objective[iColumn]<0.0) {
998                      if (basis->getStructStatus(iColumn)==
999                          CoinWarmStartBasis::basic) {
1000                        solution[iColumn]=upper[iColumn];
1001                        basis->setStructStatus(iColumn,CoinWarmStartBasis::atUpperBound);
1002                        break;
1003                      }
1004                    }
1005                  }
1006                }
1007                solver->setColSolution(solution);
1008                delete [] solution;
1009                solver->setWarmStart(basis);
1010                delete basis;
1011              }
1012            } else {
1013              // faster to do from all slack!!!! ???
1014              CoinWarmStartBasis dummy;
1015              solver->setWarmStart(&dummy);
1016            }
1017          }
1018          if (numberTries>1&&numberPasses==1&&firstPerturbedObjective) {
1019            // Modify to use convex combination
1020            // use basis from first time
1021            solver->setWarmStart(&saveBasis);
1022            // and objective
1023            if (secondPassOpt<3||(secondPassOpt>=4&&secondPassOpt<6))
1024              solver->setObjective(firstPerturbedObjective);
1025            // and solution
1026            solver->setColSolution(firstPerturbedSolution);
1027            if (secondPassOpt==2||secondPassOpt==5)
1028              solver->setHintParam(OsiDoDualInResolve,true); // dual may be better
1029          }
1030          solver->resolve();
1031          if (!solver->isProvenOptimal()) {
1032            // presumably max time or some such
1033            exitAll=true;
1034            break;
1035          }
1036          if (numberPasses==1&&secondPassOpt) {
1037            if (numberTries==1||secondPassOpt>3) {
1038              // save basis
1039              CoinWarmStartBasis * basis =
1040                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
1041              if (basis) {
1042                saveBasis = * basis;
1043                delete basis;
1044              }
1045              delete [] firstPerturbedObjective;
1046              delete [] firstPerturbedSolution;
1047              firstPerturbedObjective = CoinCopyOfArray(solver->getObjCoefficients(),numberColumns);
1048              firstPerturbedSolution = CoinCopyOfArray(solver->getColSolution(),numberColumns);
1049            }
1050          }
1051          solver->setHintParam(OsiDoDualInResolve,takeHint);
1052          newTrueSolutionValue = -saveOffset;
1053          newSumInfeas=0.0;
1054          newNumberInfeas=0;
1055          {
1056            const double * newSolution = solver->getColSolution();
1057            for (  i=0 ; i<numberColumns ; i++ ) {
1058              if (solver->isInteger(i)) {
1059                double value = newSolution[i];
1060                double nearest = floor(value+0.5);
1061                newSumInfeas += fabs(value-nearest);
1062                if (fabs(value-nearest)>1.0e-6)
1063                  newNumberInfeas++;
1064              }
1065              newTrueSolutionValue += saveObjective[i]*newSolution[i];
1066            }
1067            newTrueSolutionValue *= direction;
1068            if (newNumberInfeas&&newNumberInfeas<15) {
1069#if 0
1070              roundingObjective=solutionValue;
1071              OsiSolverInterface * saveSolver = model_->swapSolver(solver);
1072              double * currentObjective =
1073                CoinCopyOfArray(solver->getObjCoefficients(),numberColumns);
1074              solver->setObjective(saveObjective);
1075              double saveOffset2;
1076              solver->getDblParam(OsiObjOffset,saveOffset2);
1077              solver->setDblParam(OsiObjOffset,saveOffset);
1078              int ifSol = roundingHeuristic.solution(roundingObjective,roundingSolution);
1079              solver->setObjective(currentObjective);
1080              solver->setDblParam(OsiObjOffset,saveOffset2);
1081              delete [] currentObjective;
1082              model_->swapSolver(saveSolver);
1083              if (ifSol>0)
1084                abort();
1085#endif
1086              int numberRows=solver->getNumRows();
1087              double * rowActivity = new double[numberRows];
1088              memset(rowActivity,0,numberRows*sizeof(double));
1089              int * which = new int[newNumberInfeas];
1090              int * stack = new int[newNumberInfeas+1];
1091              double * baseValue = new double[newNumberInfeas];
1092              int * whichRow = new int[numberRows];
1093              double * rowValue = new double[numberRows];
1094              memset(rowValue,0,numberRows*sizeof(double));
1095              int nRow=0;
1096              // Column copy
1097              const double * element = solver->getMatrixByCol()->getElements();
1098              const int * row = solver->getMatrixByCol()->getIndices();
1099              const CoinBigIndex * columnStart = solver->getMatrixByCol()->getVectorStarts();
1100              const int * columnLength = solver->getMatrixByCol()->getVectorLengths();
1101              int n=0;
1102              double contrib = 0.0;
1103              for (  i=0 ; i<numberColumns ; i++ ) {
1104                double value = newSolution[i];
1105                if (solver->isInteger(i)) {
1106                  double nearest = floor(value+0.5);
1107                  if (fabs(value-nearest)>1.0e-6) {
1108                    //printf("Column %d value %g\n",i,value);
1109                    for (CoinBigIndex j=columnStart[i];
1110                         j<columnStart[i]+columnLength[i];j++) {
1111                      int iRow=row[j];
1112                      //printf("row %d element %g\n",iRow,element[j]);
1113                      if (!rowValue[iRow]) {
1114                        rowValue[iRow]=1.0;
1115                        whichRow[nRow++]=iRow;
1116                      }
1117                    }
1118                    baseValue[n]=floor(value);
1119                    contrib += saveObjective[i]*value;
1120                    value=0.0;
1121                    stack[n]=0;
1122                    which[n++]=i;
1123                  }
1124                }
1125                for (CoinBigIndex j=columnStart[i];
1126                     j<columnStart[i]+columnLength[i];j++) {
1127                  int iRow=row[j];
1128                  rowActivity[iRow] += value*element[j];
1129                }
1130              }
1131              if (newNumberInfeas<15) {
1132                stack[n]=newNumberInfeas+100;
1133                int iStack=n;
1134                memset(rowValue,0,numberRows*sizeof(double));
1135                const double * rowLower = solver->getRowLower();
1136                const double * rowUpper = solver->getRowUpper();
1137                while (iStack>=0) {
1138                  double contrib2=0.0;
1139                  // Could do faster
1140                  for (int k=0 ; k<n ; k++ ) {
1141                    i=which[k];
1142                    double value = baseValue[k]+stack[k];
1143                    contrib2 += saveObjective[i]*value;
1144                    for (CoinBigIndex j=columnStart[i];
1145                         j<columnStart[i]+columnLength[i];j++) {
1146                      int iRow=row[j];
1147                      rowValue[iRow] += value*element[j];
1148                    }
1149                  }
1150                  // check if feasible
1151                  bool feasible=true;
1152                  for (int k=0;k<nRow;k++) {
1153                    i=whichRow[k];
1154                    double value = rowValue[i]+rowActivity[i];
1155                    rowValue[i]=0.0;
1156                    if(value<rowLower[i]-1.0e-7||
1157                       value>rowUpper[i]+1.0e-7)
1158                      feasible=false;
1159                  }
1160                  if (feasible) {
1161                    double newObj = newTrueSolutionValue * direction;
1162                    newObj += contrib2-contrib;
1163                    newObj *= direction;
1164#ifdef COIN_DEVELOP
1165                    printf("FFFeasible! - obj %g\n",newObj);
1166#endif
1167                    if (newObj<roundingObjective-1.0e-6) {
1168#ifdef COIN_DEVELOP
1169                      printf("FBetter\n");
1170#endif
1171                      roundingObjective = newObj;
1172                      memcpy(roundingSolution,newSolution,numberColumns*sizeof(double));
1173                      for (int k=0 ; k<n ; k++ ) {
1174                        i=which[k];
1175                        double value = baseValue[k]+stack[k];
1176                        roundingSolution[i]=value;
1177                      }
1178                    }
1179                  }
1180                  while (iStack>=0&&stack[iStack]) {
1181                    stack[iStack]--;
1182                    iStack--;
1183                  }
1184                  if (iStack>=0) {
1185                    stack[iStack]=1;
1186                    iStack=n;
1187                    stack[n]=1;
1188                  }
1189                }
1190              }
1191              delete [] rowActivity;
1192              delete [] which;
1193              delete [] stack;
1194              delete [] baseValue;
1195              delete [] whichRow;
1196              delete [] rowValue;
1197            }
1198          }
1199          if (true) {
1200            OsiSolverInterface * saveSolver = model_->swapSolver(clonedSolver);
1201            clonedSolver->setColSolution(solver->getColSolution());
1202            CbcRounding heuristic1(*model_);
1203            heuristic1.setHeuristicName("rounding in feaspump!");
1204            heuristic1.setWhen(1);
1205            roundingObjective = CoinMin(roundingObjective,solutionValue);
1206            double testSolutionValue=newTrueSolutionValue;
1207            int returnCode = heuristic1.solution(roundingObjective,
1208                                                 roundingSolution,
1209                                                 testSolutionValue) ;
1210            if (returnCode==1) {
1211#ifdef COIN_DEVELOP
1212              printf("rounding obj of %g?\n",roundingObjective);
1213#endif
1214              //roundingObjective = newSolutionValue;
1215            } else {
1216              //  roundingObjective = COIN_DBL_MAX;
1217            }
1218            model_->swapSolver(saveSolver);
1219          }
1220          if (!solver->isProvenOptimal()) {
1221            // presumably max time or some such
1222            exitAll=true;
1223            break;
1224          }
1225          // in case very dubious solver
1226          lower = solver->getColLower();
1227          upper = solver->getColUpper();
1228          solution = solver->getColSolution();
1229          numberIterations = solver->getIterationCount();
1230        } else {
1231          int * addStart = new int[2*general+1];
1232          int * addIndex = new int[4*general];
1233          double * addElement = new double[4*general];
1234          double * addLower = new double[2*general];
1235          double * addUpper = new double[2*general];
1236          double * obj = new double[general];
1237          int nAdd=0;
1238          for (i=0;i<numberIntegers;i++) {
1239            int iColumn = integerVariable[i];
1240            if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
1241                newSolution[iColumn]<upper[iColumn]-primalTolerance) {
1242              assert (upper[iColumn]-lower[iColumn]>1.00001);
1243              obj[nAdd]=1.0;
1244              addLower[nAdd]=0.0;
1245              addUpper[nAdd]=COIN_DBL_MAX;
1246              nAdd++;
1247            }
1248          }
1249          OsiSolverInterface * solver2 = solver;
1250          if (nAdd) {
1251            CoinZeroN(addStart,nAdd+1);
1252            solver2 = solver->clone();
1253            solver2->addCols(nAdd,addStart,NULL,NULL,addLower,addUpper,obj);
1254            // feasible solution
1255            double * sol = new double[nAdd+numberColumns];
1256            memcpy(sol,solution,numberColumns*sizeof(double));
1257            // now rows
1258            int nAdd=0;
1259            int nEl=0;
1260            int nAddRow=0;
1261            for (i=0;i<numberIntegers;i++) {
1262              int iColumn = integerVariable[i];
1263              if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
1264                  newSolution[iColumn]<upper[iColumn]-primalTolerance) {
1265                addLower[nAddRow]=-newSolution[iColumn];;
1266                addUpper[nAddRow]=COIN_DBL_MAX;
1267                addIndex[nEl] = iColumn;
1268                addElement[nEl++]=-1.0;
1269                addIndex[nEl] = numberColumns+nAdd;
1270                addElement[nEl++]=1.0;
1271                nAddRow++;
1272                addStart[nAddRow]=nEl;
1273                addLower[nAddRow]=newSolution[iColumn];;
1274                addUpper[nAddRow]=COIN_DBL_MAX;
1275                addIndex[nEl] = iColumn;
1276                addElement[nEl++]=1.0;
1277                addIndex[nEl] = numberColumns+nAdd;
1278                addElement[nEl++]=1.0;
1279                nAddRow++;
1280                addStart[nAddRow]=nEl;
1281                sol[nAdd+numberColumns] = fabs(sol[iColumn]-newSolution[iColumn]);
1282                nAdd++;
1283              }
1284            }
1285            solver2->setColSolution(sol);
1286            delete [] sol;
1287            solver2->addRows(nAddRow,addStart,addIndex,addElement,addLower,addUpper);
1288          }
1289          delete [] addStart;
1290          delete [] addIndex;
1291          delete [] addElement;
1292          delete [] addLower;
1293          delete [] addUpper;
1294          delete [] obj;
1295          solver2->resolve();
1296          if (!solver2->isProvenOptimal()) {
1297            // presumably max time or some such
1298            exitAll=true;
1299            break;
1300          }
1301          //assert (solver2->isProvenOptimal());
1302          if (nAdd) {
1303            solver->setColSolution(solver2->getColSolution());
1304            numberIterations = solver2->getIterationCount();
1305            delete solver2;
1306          } else {
1307            numberIterations = solver->getIterationCount();
1308          }
1309          lower = solver->getColLower();
1310          upper = solver->getColUpper();
1311          solution = solver->getColSolution();
1312          newTrueSolutionValue = -saveOffset;
1313          newSumInfeas=0.0;
1314          newNumberInfeas=0;
1315          {
1316            const double * newSolution = solver->getColSolution();
1317            for (  i=0 ; i<numberColumns ; i++ ) {
1318              if (solver->isInteger(i)) {
1319                double value = newSolution[i];
1320                double nearest = floor(value+0.5);
1321                newSumInfeas += fabs(value-nearest);
1322                if (fabs(value-nearest)>1.0e-6)
1323                  newNumberInfeas++;
1324              }
1325              newTrueSolutionValue += saveObjective[i]*newSolution[i];
1326            }
1327            newTrueSolutionValue *= direction;
1328          }
1329        }
1330        if (lastMove!=1000000) {
1331          if (newSumInfeas<lastSumInfeas) {
1332            lastMove=numberPasses;
1333            lastSumInfeas=newSumInfeas;
1334          } else if (newSumInfeas>lastSumInfeas+1.0e-5) {
1335            lastMove=1000000; // going up
1336          }
1337        }
1338        totalNumberIterations += numberIterations;
1339        if (solver->getNumRows()<3000)
1340          sprintf(pumpPrint,"Pass %3d: suminf. %10.5f (%d) obj. %g iterations %d", 
1341                  numberPasses+totalNumberPasses,
1342                  newSumInfeas,newNumberInfeas,
1343                  newTrueSolutionValue,numberIterations);
1344        else
1345          sprintf(pumpPrint,"Pass %3d: (%.2f seconds) suminf. %10.5f (%d) obj. %g iterations %d", numberPasses+totalNumberPasses,
1346                  model_->getCurrentSeconds(),newSumInfeas,newNumberInfeas,
1347                  newTrueSolutionValue,numberIterations);
1348        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1349          << pumpPrint
1350          <<CoinMessageEol;
1351        if (closestSolution&&solver->getObjValue()<closestObjectiveValue) {
1352          int i;
1353          const double * objective = solver->getObjCoefficients();
1354          for (i=0;i<numberIntegersOrig;i++) {
1355            int iColumn=integerVariableOrig[i];
1356            if (objective[iColumn]>0.0)
1357              closestSolution[i]=0;
1358            else
1359              closestSolution[i]=1;
1360          }
1361          closestObjectiveValue = solver->getObjValue();
1362        }
1363        // See if we need to think about changing rhs
1364        if ((switches_&12)!=0&&useRhs<1.0e50) {
1365          double oldRhs = useRhs;
1366          bool trying=false;
1367          if ((switches_&4)!=0&&numberPasses&&(numberPasses%50)==0) {
1368            if (solutionValue>1.0e20) {
1369              // only if no genuine solution
1370              double gap = useRhs-continuousObjectiveValue;
1371              useRhs += 0.1*gap;
1372              if (exactMultiple) {
1373                useRhs = exactMultiple*ceil(useRhs/exactMultiple);
1374                useRhs = CoinMax(useRhs,oldRhs+exactMultiple);
1375              }
1376              trying=true;
1377            }
1378          }
1379          if ((switches_&8)!=0) {
1380            // Put in new suminf and check
1381            double largest=newSumInfeas;
1382            double smallest=newSumInfeas;
1383            for (int i=0;i<SIZE_BOBBLE-1;i++) {
1384              double value = saveSumInf[i+1];
1385              saveSumInf[i]=value;
1386              largest = CoinMax(largest,value);
1387              smallest = CoinMin(smallest,value);
1388            }
1389            saveSumInf[SIZE_BOBBLE-1]=newSumInfeas;
1390            if (smallest*1.5>largest&&smallest>2.0) {
1391              if (bobbleMode==0) {
1392                // go closer
1393                double gap = oldRhs-continuousObjectiveValue;
1394                useRhs -= 0.4*gap;
1395                if (exactMultiple) {
1396                  double value = floor(useRhs/exactMultiple);
1397                  useRhs = CoinMin(value*exactMultiple,oldRhs-exactMultiple);
1398                }
1399                if (useRhs<continuousObjectiveValue) {
1400                  // skip decrease
1401                  bobbleMode=1;
1402                  useRhs=oldRhs;
1403                }
1404              } 
1405              if (bobbleMode) {
1406                trying=true;
1407                // weaken
1408                if (solutionValue<1.0e20) {
1409                  double gap = solutionValue-oldRhs;
1410                  useRhs += 0.3*gap;
1411                } else {
1412                  double gap = oldRhs-continuousObjectiveValue;
1413                  useRhs += 0.05*gap;
1414                }
1415                if (exactMultiple) {
1416                  double value = ceil(useRhs/exactMultiple);
1417                  useRhs = CoinMin(value*exactMultiple,
1418                                   solutionValue-exactMultiple);
1419                }
1420              }
1421              bobbleMode++;
1422              // reset
1423              CoinFillN(saveSumInf,SIZE_BOBBLE,COIN_DBL_MAX);
1424            }
1425          }
1426          if (useRhs!=oldRhs) {
1427            // tidy up
1428            if (exactMultiple) {
1429              double value = floor(useRhs/exactMultiple);
1430              double bestPossible = ceil(continuousObjectiveValue/exactMultiple);
1431              useRhs = CoinMax(value,bestPossible)*exactMultiple;
1432            } else {
1433              useRhs = CoinMax(useRhs,continuousObjectiveValue);
1434            }
1435            int k = solver->getNumRows()-1;
1436            solver->setRowUpper(k,useRhs+useOffset);
1437            bool takeHint;
1438            OsiHintStrength strength;
1439            solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
1440            if (useRhs<oldRhs) {
1441              solver->setHintParam(OsiDoDualInResolve,true);
1442              solver->resolve();
1443            } else if (useRhs>oldRhs) {
1444              solver->setHintParam(OsiDoDualInResolve,false);
1445              solver->resolve();
1446            }
1447            solver->setHintParam(OsiDoDualInResolve,takeHint);
1448            if (!solver->isProvenOptimal()) {
1449              // presumably max time or some such
1450              exitAll=true;
1451              break;
1452            }
1453          } else if (trying) {
1454            // doesn't look good
1455            break;
1456          }
1457        }
1458      }
1459      // reduce scale factor
1460      scaleFactor *= weightFactor_;
1461    } // END WHILE
1462    // see if rounding worked!
1463    if (roundingObjective<solutionValue) {
1464      sprintf(pumpPrint,"Rounding solution of %g is better than previous of %g !\n",
1465              roundingObjective,solutionValue);
1466      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1467        << pumpPrint
1468        <<CoinMessageEol;
1469      solutionValue=roundingObjective;
1470      newSolutionValue = solutionValue;
1471      memcpy(betterSolution,roundingSolution,numberColumns*sizeof(double));
1472      solutionFound=true;
1473      if (exitNow(roundingObjective))
1474        exitAll=true;
1475    }
1476    if (!solutionFound) { 
1477      sprintf(pumpPrint,"No solution found this major pass");
1478      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1479        << pumpPrint
1480        <<CoinMessageEol;
1481    }
1482    //}
1483    delete solver;
1484    solver=NULL;
1485    for ( j=0;j<NUMBER_OLD;j++) 
1486      delete [] oldSolution[j];
1487    delete [] oldSolution;
1488    delete [] saveObjective;
1489    if (usedColumn&&!exitAll) {
1490      OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
1491      const double * colLower = newSolver->getColLower();
1492      const double * colUpper = newSolver->getColUpper();
1493      bool stopBAB=false;
1494      int allowedPass=-1;
1495      while (!stopBAB) {
1496        stopBAB=true;
1497        int i;
1498        int nFix=0;
1499        int nFixI=0;
1500        int nFixC=0;
1501        int nFixC2=0;
1502        for (i=0;i<numberIntegersOrig;i++) {
1503          int iColumn=integerVariableOrig[i];
1504          //const OsiObject * object = model_->object(i);
1505          //double originalLower;
1506          //double originalUpper;
1507          //getIntegerInformation( object,originalLower, originalUpper);
1508          //assert(colLower[iColumn]==originalLower);
1509          //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
1510          newSolver->setColLower(iColumn,colLower[iColumn]);
1511          //assert(colUpper[iColumn]==originalUpper);
1512          //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
1513          newSolver->setColUpper(iColumn,colUpper[iColumn]);
1514          if (usedColumn[iColumn]<=allowedPass) {
1515            double value=lastSolution[iColumn];
1516            double nearest = floor(value+0.5);
1517            if (fabs(value-nearest)<1.0e-7) {
1518              if (nearest==colLower[iColumn]) {
1519                newSolver->setColUpper(iColumn,colLower[iColumn]);
1520                nFix++;
1521              } else if (nearest==colUpper[iColumn]) {
1522                newSolver->setColLower(iColumn,colUpper[iColumn]);
1523                nFix++;
1524              } else if (fixInternal) {
1525                newSolver->setColLower(iColumn,nearest);
1526                newSolver->setColUpper(iColumn,nearest);
1527                nFix++;
1528                nFixI++;
1529              }
1530            }
1531          }
1532        }
1533        if (fixContinuous) {
1534          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1535            if (!newSolver->isInteger(iColumn)&&usedColumn[iColumn]<=allowedPass) {
1536              double value=lastSolution[iColumn];
1537              if (value<colLower[iColumn]+1.0e-8) {
1538                newSolver->setColUpper(iColumn,colLower[iColumn]);
1539                nFixC++;
1540              } else if (value>colUpper[iColumn]-1.0e-8) {
1541                newSolver->setColLower(iColumn,colUpper[iColumn]);
1542                nFixC++;
1543              } else if (fixContinuous==2) {
1544                newSolver->setColLower(iColumn,value);
1545                newSolver->setColUpper(iColumn,value);
1546                nFixC++;
1547                nFixC2++;
1548              }
1549            }
1550          }
1551        }
1552        newSolver->initialSolve();
1553        if (!newSolver->isProvenOptimal()) {
1554          //newSolver->writeMps("bad.mps");
1555          //assert (newSolver->isProvenOptimal());
1556          exitAll=true;
1557          break;
1558        }
1559        sprintf(pumpPrint,"Before mini branch and bound, %d integers at bound fixed and %d continuous",
1560                nFix,nFixC);
1561        if (nFixC2+nFixI!=0)
1562          sprintf(pumpPrint+strlen(pumpPrint)," of which %d were internal integer and %d internal continuous",
1563                  nFixI,nFixC2);
1564        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1565          << pumpPrint
1566          <<CoinMessageEol;
1567        double saveValue = newSolutionValue;
1568        if (newSolutionValue-model_->getCutoffIncrement()
1569            >continuousObjectiveValue-1.0e-7) {
1570          double saveFraction = fractionSmall_;
1571          if (numberTries>1&&!numberBandBsolutions)
1572            fractionSmall_ *= 0.5;
1573          // Give branch and bound a bit more freedom
1574          double cutoff2=newSolutionValue+
1575            CoinMax(model_->getCutoffIncrement(),1.0e-3);
1576          int returnCode2 = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
1577                                                cutoff2,"CbcHeuristicLocalAfterFPump");
1578          fractionSmall_ = saveFraction;
1579          if (returnCode2<0) {
1580            if (returnCode2==-2) {
1581              exitAll=true;
1582              returnCode=0;
1583            } else {
1584              returnCode2=0; // returned on size - try changing
1585              //#define ROUND_AGAIN
1586#ifdef ROUND_AGAIN
1587              if (numberTries==1&&numberPasses>20&&allowedPass<numberPasses-1) {
1588                allowedPass=(numberPasses+allowedPass)>>1;
1589                sprintf(pumpPrint,
1590                        "Fixing all variables which were last changed on pass %d and trying again",
1591                       allowedPass);
1592                model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1593                  << pumpPrint
1594                  <<CoinMessageEol;
1595                stopBAB=false;
1596                continue;
1597              }
1598#endif
1599            }
1600          }
1601          if ((returnCode2&2)!=0) {
1602            // could add cut
1603            returnCode2 &= ~2;
1604          }
1605          if (returnCode2)
1606            numberBandBsolutions++;
1607        } else {
1608          // no need
1609          exitAll=true;
1610          //returnCode=0;
1611        }
1612        // recompute solution value
1613        if (returnCode&&true) {
1614          delete newSolver;
1615          newSolver = model_->continuousSolver()->clone();
1616          newSolutionValue = -saveOffset;
1617          double newSumInfeas=0.0;
1618          const double * obj = newSolver->getObjCoefficients();
1619          for (int i=0 ; i<numberColumns ; i++ ) {
1620            if (newSolver->isInteger(i)) {
1621              double value = newSolution[i];
1622              double nearest = floor(value+0.5);
1623              newSumInfeas += fabs(value-nearest);
1624            }
1625            newSolutionValue += obj[i]*newSolution[i];
1626          }
1627          newSolutionValue *= direction;
1628        }
1629        bool gotSolution = false;
1630        if (returnCode&&newSolutionValue<saveValue) {
1631          sprintf(pumpPrint,"Mini branch and bound improved solution from %g to %g (%.2f seconds)",
1632                  saveValue,newSolutionValue,model_->getCurrentSeconds());
1633          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1634            << pumpPrint
1635            <<CoinMessageEol;
1636          memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
1637          gotSolution=true;
1638          if (fixContinuous&&nFixC+nFixC2>0) {
1639            // may be able to do even better
1640            const double * lower = model_->solver()->getColLower();
1641            const double * upper = model_->solver()->getColUpper();
1642            for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1643              if (newSolver->isInteger(iColumn)) {
1644                double value=floor(newSolution[iColumn]+0.5);
1645                newSolver->setColLower(iColumn,value);
1646                newSolver->setColUpper(iColumn,value);
1647              } else {
1648                newSolver->setColLower(iColumn,lower[iColumn]);
1649                newSolver->setColUpper(iColumn,upper[iColumn]);
1650              }
1651            }
1652            newSolver->initialSolve();
1653            if (newSolver->isProvenOptimal()) {
1654              double value = newSolver->getObjValue()*newSolver->getObjSense();
1655              if (value<newSolutionValue) {
1656                //newSolver->writeMps("query","mps");
1657#if 0
1658                {
1659                  double saveOffset;
1660                  newSolver->getDblParam(OsiObjOffset,saveOffset);
1661                  const double * obj = newSolver->getObjCoefficients();
1662                  double newTrueSolutionValue = -saveOffset;
1663                  double newSumInfeas=0.0;
1664                  int numberColumns = newSolver->getNumCols();
1665                  const double * solution = newSolver->getColSolution();
1666                  for (int  i=0 ; i<numberColumns ; i++ ) {
1667                    if (newSolver->isInteger(i)) {
1668                      double value = solution[i];
1669                      double nearest = floor(value+0.5);
1670                      newSumInfeas += fabs(value-nearest);
1671                    }
1672                    if (solution[i])
1673                      printf("%d obj %g val %g - total %g\n",i,obj[i],solution[i],
1674                             newTrueSolutionValue);
1675                    newTrueSolutionValue += obj[i]*solution[i];
1676                  }
1677                  printf("obj %g - inf %g\n",newTrueSolutionValue,
1678                         newSumInfeas);
1679                }
1680#endif
1681                sprintf(pumpPrint,"Freeing continuous variables gives a solution of %g", value);
1682                model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1683                  << pumpPrint
1684                  <<CoinMessageEol;
1685                newSolutionValue=value;
1686                memcpy(betterSolution,newSolver->getColSolution(),numberColumns*sizeof(double));
1687              }
1688            } else {
1689              //newSolver->writeMps("bad3.mps");
1690              exitAll=true;
1691              break;
1692            }
1693          } 
1694        } else {
1695          sprintf(pumpPrint,"Mini branch and bound did not improve solution (%.2f seconds)",
1696                  model_->getCurrentSeconds());
1697          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1698            << pumpPrint
1699            <<CoinMessageEol;
1700          if (returnCode&&newSolutionValue<saveValue+1.0e-3&&nFixC+nFixC2) {
1701            // may be able to do better
1702            const double * lower = model_->solver()->getColLower();
1703            const double * upper = model_->solver()->getColUpper();
1704            for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1705              if (newSolver->isInteger(iColumn)) {
1706                double value=floor(newSolution[iColumn]+0.5);
1707                newSolver->setColLower(iColumn,value);
1708                newSolver->setColUpper(iColumn,value);
1709              } else {
1710                newSolver->setColLower(iColumn,lower[iColumn]);
1711                newSolver->setColUpper(iColumn,upper[iColumn]);
1712              }
1713            }
1714            newSolver->initialSolve();
1715            if (newSolver->isProvenOptimal()) {
1716              double value = newSolver->getObjValue()*newSolver->getObjSense();
1717              if (value<saveValue) {
1718                sprintf(pumpPrint,"Freeing continuous variables gives a solution of %g", value);
1719                model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1720                  << pumpPrint
1721                  <<CoinMessageEol;
1722                newSolutionValue=value;
1723                memcpy(betterSolution,newSolver->getColSolution(),numberColumns*sizeof(double));
1724                gotSolution=true;
1725              }
1726            }
1727          }
1728        }
1729        if (gotSolution) {
1730          if ((accumulate_&1)!=0) {
1731            model_->incrementUsed(betterSolution); // for local search
1732            numberSolutions++;
1733          }
1734          solutionValue=newSolutionValue;
1735          solutionFound=true;
1736          if (exitNow(newSolutionValue))
1737            exitAll=true;
1738          CoinWarmStartBasis * basis =
1739            dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
1740          if (basis) {
1741            bestBasis = * basis;
1742            delete basis;
1743            CbcEventHandler * handler = model_->getEventHandler();
1744            if (handler) {
1745              double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
1746              double saveObjectiveValue = model_->getMinimizationObjValue();
1747              model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
1748              int action = handler->event(CbcEventHandler::heuristicSolution);
1749              //printf("cutoff %g\n",model_->getCutoff());
1750              if (saveOldSolution&&saveObjectiveValue<model_->getMinimizationObjValue())
1751                model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
1752              delete [] saveOldSolution;
1753              if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
1754                exitAll=true; // exit
1755                break;
1756              }
1757            }
1758          }
1759        }
1760      } // end stopBAB while
1761      delete newSolver;
1762    }
1763    if (solutionFound) finalReturnCode=1;
1764    cutoff = CoinMin(cutoff,solutionValue-model_->getCutoffIncrement());
1765    if (numberTries>=maximumRetries_||!solutionFound||exitAll||cutoff<continuousObjectiveValue+1.0e-7) {
1766      break;
1767    } else {
1768      solutionFound=false;
1769      if (absoluteIncrement_>0.0||relativeIncrement_>0.0) {
1770        double gap = relativeIncrement_*fabs(solutionValue);
1771        cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
1772      } else {
1773        //double weights[10]={0.1,0.1,0.2,0.2,0.2,0.3,0.3,0.3,0.4,0.5};
1774        double weights[10]={0.1,0.2,0.3,0.3,0.4,0.4,0.4,0.5,0.5,0.6};
1775        cutoff -= weights[CoinMin(numberTries-1,9)]*(cutoff-continuousObjectiveValue);
1776      }
1777      // But round down
1778      if (exactMultiple) 
1779        cutoff = exactMultiple*floor(cutoff/exactMultiple);
1780      if (cutoff<continuousObjectiveValue)
1781        break;
1782      sprintf(pumpPrint,"Round again with cutoff of %g",cutoff);
1783      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1784        << pumpPrint
1785        <<CoinMessageEol;
1786      if ((accumulate_&3)<2&&usedColumn) {
1787        for (int i=0;i<numberColumns;i++)
1788          usedColumn[i]=-1;
1789      }
1790      if (!totalNumberPasses)
1791        numberIterationsPass1 = totalNumberIterations;
1792      numberIterationsLastPass =  totalNumberIterations;
1793      totalNumberPasses += numberPasses-1;
1794    }
1795  }
1796  delete solver; // probably NULL but do anyway
1797  if (!finalReturnCode&&closestSolution&&closestObjectiveValue <= 10.0&&usedColumn) {
1798    // try a bit of branch and bound
1799    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
1800    const double * colLower = newSolver->getColLower();
1801    const double * colUpper = newSolver->getColUpper();
1802    int i;
1803    double rhs = 0.0;
1804    for (i=0;i<numberIntegersOrig;i++) {
1805      int iColumn=integerVariableOrig[i];
1806      int direction = closestSolution[i];
1807      closestSolution[i]=iColumn;
1808      if (direction==0) {
1809        // keep close to LB
1810        rhs += colLower[iColumn];
1811        lastSolution[i]=1.0;
1812      } else {
1813        // keep close to UB
1814        rhs -= colUpper[iColumn];
1815        lastSolution[i]=-1.0;
1816      }
1817    }
1818    newSolver->addRow(numberIntegersOrig,closestSolution,
1819                      lastSolution,-COIN_DBL_MAX,rhs+10.0);
1820    //double saveValue = newSolutionValue;
1821    //newSolver->writeMps("sub");
1822    int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
1823                                     newSolutionValue,"CbcHeuristicLocalAfterFPump");
1824    if (returnCode<0)
1825      returnCode=0; // returned on size
1826    if ((returnCode&2)!=0) {
1827      // could add cut
1828      returnCode &= ~2;
1829    }
1830    if (returnCode) {
1831      //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
1832      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
1833      //abort();
1834      solutionValue=newSolutionValue;
1835      solutionFound=true;
1836      if (exitNow(newSolutionValue))
1837        exitAll=true;
1838    }
1839    delete newSolver;
1840  }
1841  delete clonedSolver;
1842  delete [] roundingSolution;
1843  delete [] usedColumn;
1844  delete [] lastSolution;
1845  delete [] newSolution;
1846  delete [] closestSolution;
1847  delete [] integerVariable;
1848  delete [] firstPerturbedObjective;
1849  delete [] firstPerturbedSolution;
1850  if (solutionValue==incomingObjective) 
1851    sprintf(pumpPrint,"After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
1852            model_->getCurrentSeconds(),CoinCpuTime()-time1);
1853  else
1854    sprintf(pumpPrint,"After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds",
1855            model_->getCurrentSeconds(),solutionValue,CoinCpuTime()-time1);
1856  model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1857    << pumpPrint
1858    <<CoinMessageEol;
1859  if (bestBasis.getNumStructural())
1860    model_->setBestSolutionBasis(bestBasis);
1861  model_->setMinimizationObjValue(saveBestObjective);
1862  if ((accumulate_&1)!=0&&numberSolutions>1&&!model_->getSolutionCount()) {
1863    model_->setSolutionCount(1); // for local search
1864    model_->setNumberHeuristicSolutions(1); 
1865  }
1866#ifdef COIN_DEVELOP
1867  {
1868    double ncol = model_->solver()->getNumCols();
1869    double nrow = model_->solver()->getNumRows();
1870    printf("XXX total iterations %d ratios - %g %g %g\n",
1871           totalNumberIterations,
1872           static_cast<double> (totalNumberIterations)/nrow,
1873           static_cast<double> (totalNumberIterations)/ncol,
1874           static_cast<double> (totalNumberIterations)/(2*nrow+2*ncol));
1875  }
1876#endif
1877  return finalReturnCode;
1878}
1879
1880/**************************END MAIN PROCEDURE ***********************************/
1881
1882// update model
1883void CbcHeuristicFPump::setModel(CbcModel * model)
1884{
1885  model_ = model;
1886}
1887
1888/* Rounds solution - down if < downValue
1889   returns 1 if current is a feasible solution
1890*/
1891int 
1892CbcHeuristicFPump::rounds(OsiSolverInterface * solver,double * solution,
1893                          const double * objective,
1894                          int numberIntegers, const int * integerVariable,
1895                          char * pumpPrint, int iter,
1896                          bool roundExpensive, double downValue, int *flip)
1897{
1898  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1899  double primalTolerance ;
1900  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
1901
1902  int i;
1903
1904  const double * cost = solver->getObjCoefficients();
1905  int flip_up = 0;
1906  int flip_down  = 0;
1907  double  v = randomNumberGenerator_.randomDouble() * 20.0;
1908  int nn = 10 + static_cast<int> (v);
1909  int nnv = 0;
1910  int * list = new int [nn];
1911  double * val = new double [nn];
1912  for (i = 0; i < nn; i++) val[i] = .001;
1913
1914  const double * rowLower = solver->getRowLower();
1915  const double * rowUpper = solver->getRowUpper();
1916  int numberRows = solver->getNumRows();
1917  if ((iter&1)!=0) {
1918    // Do set covering variables
1919    const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
1920    const double * elementByRow = matrixByRow->getElements();
1921    const int * column = matrixByRow->getIndices();
1922    const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
1923    const int * rowLength = matrixByRow->getVectorLengths();
1924    for (i=0;i<numberRows;i++) {
1925      if (rowLower[i]==1.0&&rowUpper[i]==1.0) {
1926        bool cover=true;
1927        double largest=0.0;
1928        int jColumn=-1;
1929        for (CoinBigIndex k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
1930          int iColumn = column[k];
1931          if (elementByRow[k]!=1.0||!solver->isInteger(iColumn)) {
1932            cover=false;
1933            break;
1934          } else {
1935            if (solution[iColumn]) {
1936              double value = solution[iColumn]*
1937                (randomNumberGenerator_.randomDouble()+5.0);
1938              if (value>largest) {
1939                largest=value;
1940                jColumn=iColumn;
1941              }
1942            }
1943          }
1944        }
1945        if (cover) {
1946          for (CoinBigIndex k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
1947            int iColumn = column[k];
1948            if (iColumn==jColumn) 
1949              solution[iColumn]=1.0;
1950            else 
1951              solution[iColumn]=0.0;
1952          }
1953        }
1954      }
1955    }
1956  }
1957  int numberColumns = solver->getNumCols();
1958#if 0
1959  // Do set covering variables
1960  const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
1961  const double * elementByRow = matrixByRow->getElements();
1962  const int * column = matrixByRow->getIndices();
1963  const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
1964  const int * rowLength = matrixByRow->getVectorLengths();
1965  double * sortTemp = new double[numberColumns];
1966  int * whichTemp = new int [numberColumns];
1967  char * rowTemp = new char [numberRows];
1968  memset(rowTemp,0,numberRows);
1969  for (i=0;i<numberColumns;i++)
1970    whichTemp[i]=-1;
1971  int nSOS=0;
1972  for (i=0;i<numberRows;i++) {
1973    if (rowLower[i]==1.0&&rowUpper[i]==1.0) {
1974      bool cover=true;
1975      for (CoinBigIndex k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
1976        int iColumn = column[k];
1977        if (elementByRow[k]!=1.0||!solver->isInteger(iColumn)) {
1978          cover=false;
1979          break;
1980        }
1981      }
1982      if (cover) {
1983        rowTemp[i]=1;
1984        nSOS++;
1985        for (CoinBigIndex k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
1986          int iColumn = column[k];
1987          double value=solution[iColumn];
1988          whichTemp[iColumn]=iColumn;
1989        }
1990      }
1991    }
1992  }
1993  if (nSOS) {
1994    // Column copy
1995    const CoinPackedMatrix * matrixByColumn = solver->getMatrixByCol();
1996    //const double * element = matrixByColumn->getElements();
1997    const int * row = matrixByColumn->getIndices();
1998    const CoinBigIndex * columnStart = matrixByColumn->getVectorStarts();
1999    const int * columnLength = matrixByColumn->getVectorLengths();
2000    int nLook = 0;
2001    for (i=0;i<numberColumns;i++) {
2002      if (whichTemp[i]>=0) {
2003        whichTemp[nLook]=i;
2004        double value=solution[i];
2005        if (value<0.5)
2006          value *= (0.1*randomNumberGenerator_.randomDouble()+0.3);
2007        sortTemp[nLook++]=-value;
2008      }
2009    }
2010    CoinSort_2(sortTemp,sortTemp+nLook,whichTemp);
2011    double smallest=1.0;
2012    int nFix=0;
2013    int nOne=0;
2014    for (int j=0;j<nLook;j++) {
2015      int jColumn = whichTemp[j];
2016      double thisValue = solution[jColumn];
2017      if (!thisValue)
2018        continue;
2019      if (thisValue==1.0)
2020        nOne++;
2021      smallest = CoinMin(smallest,thisValue);
2022      solution[jColumn]=1.0;
2023      double largest=0.0;
2024      for (CoinBigIndex jEl=columnStart[jColumn];
2025           jEl<columnStart[jColumn]+columnLength[jColumn];jEl++) {
2026        int jRow = row[jEl];
2027        if (rowTemp[jRow]) {
2028          for (CoinBigIndex k=rowStart[jRow];k<rowStart[jRow]+rowLength[jRow];k++) {
2029            int iColumn = column[k];
2030            if (solution[iColumn]) {
2031              if (iColumn!=jColumn) {
2032                double value = solution[iColumn];
2033                if (value>largest)
2034                  largest=value;
2035                solution[iColumn]=0.0;
2036              }
2037            }
2038          }
2039        }
2040      }
2041      if (largest>thisValue)
2042        printf("%d was at %g - chosen over a value of %g\n",
2043               jColumn,thisValue,largest);
2044      nFix++;
2045    }
2046    printf("%d fixed out of %d (%d at one already)\n",
2047           nFix,nLook,nOne);
2048  }
2049  delete [] sortTemp;
2050  delete [] whichTemp;
2051  delete [] rowTemp;
2052#endif
2053  const double * columnLower = solver->getColLower();
2054  const double * columnUpper = solver->getColUpper();
2055  // Check if valid with current solution (allow for 0.99999999s)
2056  for (i=0;i<numberIntegers;i++) {
2057    int iColumn = integerVariable[i];
2058    double value=solution[iColumn];
2059    double round = floor(value+0.5);
2060    if (fabs(value-round)>primalTolerance) 
2061      break;
2062  }
2063  if (i==numberIntegers) {
2064    // may be able to use solution even if 0.99999's
2065    double * saveLower = CoinCopyOfArray(columnLower,numberColumns);
2066    double * saveUpper = CoinCopyOfArray(columnUpper,numberColumns);
2067    double * saveSolution = CoinCopyOfArray(solution,numberColumns);
2068    double * tempSolution = CoinCopyOfArray(solution,numberColumns);
2069    CoinWarmStartBasis  * saveBasis =
2070      dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
2071    for (i=0;i<numberIntegers;i++) {
2072      int iColumn = integerVariable[i];
2073      double value=solution[iColumn];
2074      double round = floor(value+0.5);
2075      solver->setColLower(iColumn,round);
2076      solver->setColUpper(iColumn,round);
2077      tempSolution[iColumn] = round;
2078    }
2079    solver->setColSolution(tempSolution);
2080    delete [] tempSolution;
2081    solver->resolve();
2082    solver->setColLower(saveLower);
2083    solver->setColUpper(saveUpper);
2084    solver->setWarmStart(saveBasis);
2085    delete [] saveLower;
2086    delete [] saveUpper;
2087    delete saveBasis;
2088    if (!solver->isProvenOptimal()) {
2089      solver->setColSolution(saveSolution);
2090    }
2091    delete [] saveSolution;
2092    if (solver->isProvenOptimal()) {
2093      // feasible
2094      delete [] list; 
2095      delete [] val;
2096      return 1;
2097    }
2098  }
2099  //double * saveSolution = CoinCopyOfArray(solution,numberColumns);
2100  // return rounded solution
2101  for (i=0;i<numberIntegers;i++) {
2102    int iColumn = integerVariable[i];
2103    double value=solution[iColumn];
2104    double round = floor(value+primalTolerance);
2105    if (value-round > downValue) round += 1.;
2106#if 1
2107    if (round < integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2108    if (round > 1. - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2109#else
2110    if (round < columnLower[iColumn]+integerTolerance && cost[iColumn] < -1. + integerTolerance) flip_down++;
2111    if (round > columnUpper[iColumn] - integerTolerance && cost[iColumn] > 1. - integerTolerance) flip_up++;
2112#endif
2113    if (flip_up + flip_down == 0) { 
2114       for (int k = 0; k < nn; k++) {
2115           if (fabs(value-round) > val[k]) {
2116              nnv++;
2117              for (int j = nn-2; j >= k; j--) {
2118                  val[j+1] = val[j];
2119                  list[j+1] = list[j];
2120              } 
2121              val[k] = fabs(value-round);
2122              list[k] = iColumn;
2123              break;
2124           }
2125       }
2126    }
2127    solution[iColumn] = round;
2128  }
2129
2130  if (nnv > nn) nnv = nn;
2131  //if (iter != 0)
2132  //sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
2133  *flip = flip_up + flip_down;
2134
2135  if (*flip == 0 && iter != 0) {
2136    //sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
2137     for (i = 0; i < nnv; i++) {
2138       // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
2139       int index = list[i];
2140       double value = solution[index];
2141       if (value<=1.0) {
2142         solution[index] = 1.0-value;
2143       } else if (value<columnLower[index]+integerTolerance) {
2144         solution[index] = value+1.0;
2145       } else if (value>columnUpper[index]-integerTolerance) {
2146         solution[index] = value-1.0;
2147       } else {
2148         solution[index] = value-1.0;
2149       }
2150     }
2151     *flip = nnv;
2152  } else {
2153    //sprintf(pumpPrint+strlen(pumpPrint)," ");
2154  }
2155  delete [] list; delete [] val;
2156  //iter++;
2157   
2158  // get row activities
2159  double * rowActivity = new double[numberRows];
2160  memset(rowActivity,0,numberRows*sizeof(double));
2161  solver->getMatrixByCol()->times(solution,rowActivity) ;
2162  double largestInfeasibility =primalTolerance;
2163  double sumInfeasibility=0.0;
2164  int numberBadRows=0;
2165  for (i=0 ; i < numberRows ; i++) {
2166    double value;
2167    value = rowLower[i]-rowActivity[i];
2168    if (value>primalTolerance) {
2169      numberBadRows++;
2170      largestInfeasibility = CoinMax(largestInfeasibility,value);
2171      sumInfeasibility += value;
2172    }
2173    value = rowActivity[i]-rowUpper[i];
2174    if (value>primalTolerance) {
2175      numberBadRows++;
2176      largestInfeasibility = CoinMax(largestInfeasibility,value);
2177      sumInfeasibility += value;
2178    }
2179  }
2180#if 0
2181  if (largestInfeasibility>primalTolerance&&numberBadRows*10<numberRows) {
2182    // Can we improve by flipping
2183    for (int iPass=0;iPass<10;iPass++) {
2184      int numberColumns = solver->getNumCols();
2185      const CoinPackedMatrix * matrixByCol = solver->getMatrixByCol();
2186      const double * element = matrixByCol->getElements();
2187      const int * row = matrixByCol->getIndices();
2188      const CoinBigIndex * columnStart = matrixByCol->getVectorStarts();
2189      const int * columnLength = matrixByCol->getVectorLengths();
2190      double oldSum = sumInfeasibility;
2191      // First improve by moving continuous ones
2192      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
2193        if (!solver->isInteger(iColumn)) {
2194          double solValue = solution[iColumn];
2195          double thetaUp = columnUpper[iColumn]-solValue;
2196          double improvementUp=0.0;
2197          if (thetaUp>primalTolerance) {
2198            // can go up
2199            for (CoinBigIndex j=columnStart[iColumn];
2200                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
2201              int iRow = row[j];
2202              double distanceUp = rowUpper[iRow]-rowActivity[iRow];
2203              double distanceDown = rowLower[iRow]-rowActivity[iRow];
2204              double el = element[j];
2205              if (el>0.0) {
2206                // positive element
2207                if (distanceUp>0.0) {
2208                  if (thetaUp*el>distanceUp)
2209                    thetaUp=distanceUp/el;
2210                } else {
2211                  improvementUp -= el;
2212                }
2213                if (distanceDown>0.0) {
2214                  if (thetaUp*el>distanceDown)
2215                    thetaUp=distanceDown/el;
2216                  improvementUp += el;
2217                }
2218              } else {
2219                // negative element
2220                if (distanceDown<0.0) {
2221                  if (thetaUp*el<distanceDown)
2222                    thetaUp=distanceDown/el;
2223                } else {
2224                  improvementUp += el;
2225                }
2226                if (distanceUp<0.0) {
2227                  if (thetaUp*el<distanceUp)
2228                    thetaUp=distanceUp/el;
2229                  improvementUp -= el;
2230                }
2231              }
2232            }
2233          }
2234          double thetaDown = solValue-columnLower[iColumn];
2235          double improvementDown=0.0;
2236          if (thetaDown>primalTolerance) {
2237            // can go down
2238            for (CoinBigIndex j=columnStart[iColumn];
2239                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
2240              int iRow = row[j];
2241              double distanceUp = rowUpper[iRow]-rowActivity[iRow];
2242              double distanceDown = rowLower[iRow]-rowActivity[iRow];
2243              double el = -element[j]; // not change in sign form up
2244              if (el>0.0) {
2245                // positive element
2246                if (distanceUp>0.0) {
2247                  if (thetaDown*el>distanceUp)
2248                    thetaDown=distanceUp/el;
2249                } else {
2250                  improvementDown -= el;
2251                }
2252                if (distanceDown>0.0) {
2253                  if (thetaDown*el>distanceDown)
2254                    thetaDown=distanceDown/el;
2255                  improvementDown += el;
2256                }
2257              } else {
2258                // negative element
2259                if (distanceDown<0.0) {
2260                  if (thetaDown*el<distanceDown)
2261                    thetaDown=distanceDown/el;
2262                } else {
2263                  improvementDown += el;
2264                }
2265                if (distanceUp<0.0) {
2266                  if (thetaDown*el<distanceUp)
2267                    thetaDown=distanceUp/el;
2268                  improvementDown -= el;
2269                }
2270              }
2271            }
2272            if (thetaUp<1.0e-8)
2273              improvementUp=0.0;
2274            if (thetaDown<1.0e-8)
2275              improvementDown=0.0;
2276            double theta;
2277            if (improvementUp>=improvementDown) {
2278              theta=thetaUp;
2279            } else {
2280              improvementUp=improvementDown;
2281              theta=-thetaDown;
2282            }
2283            if (improvementUp>1.0e-8&&fabs(theta)>1.0e-8) {
2284              // Could move
2285              double oldSum=0.0;
2286              double newSum=0.0;
2287              solution[iColumn] += theta;
2288              for (CoinBigIndex j=columnStart[iColumn];
2289                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
2290                int iRow = row[j];
2291                double lower = rowLower[iRow];
2292                double upper = rowUpper[iRow];
2293                double value = rowActivity[iRow];
2294                if (value>upper)
2295                  oldSum += value-upper;
2296                else if (value<lower)
2297                  oldSum += lower-value;
2298                value += theta*element[j];
2299                rowActivity[iRow]=value;
2300                if (value>upper)
2301                  newSum += value-upper;
2302                else if (value<lower)
2303                  newSum += lower-value;
2304              }
2305              assert (newSum<=oldSum);
2306              sumInfeasibility += newSum-oldSum;
2307            }
2308          }
2309        }
2310      }
2311      // Now flip some integers?
2312#if 0
2313      for (i=0;i<numberIntegers;i++) {
2314        int iColumn = integerVariable[i];
2315        double solValue = solution[iColumn];
2316        assert (fabs(solValue-floor(solValue+0.5))<1.0e-8);
2317        double improvementUp=0.0;
2318        if (columnUpper[iColumn]>=solValue+1.0) {
2319          // can go up
2320          double oldSum=0.0;
2321          double newSum=0.0;
2322          for (CoinBigIndex j=columnStart[iColumn];
2323               j<columnStart[iColumn]+columnLength[iColumn];j++) {
2324            int iRow = row[j];
2325            double lower = rowLower[iRow];
2326            double upper = rowUpper[iRow];
2327            double value = rowActivity[iRow];
2328            if (value>upper)
2329              oldSum += value-upper;
2330            else if (value<lower)
2331              oldSum += lower-value;
2332            value += element[j];
2333            if (value>upper)
2334              newSum += value-upper;
2335            else if (value<lower)
2336              newSum += lower-value;
2337          }
2338          improvementUp = oldSum-newSum;
2339        }
2340        double improvementDown=0.0;
2341        if (columnLower[iColumn]<=solValue-1.0) {
2342          // can go down
2343          double oldSum=0.0;
2344          double newSum=0.0;
2345          for (CoinBigIndex j=columnStart[iColumn];
2346               j<columnStart[iColumn]+columnLength[iColumn];j++) {
2347            int iRow = row[j];
2348            double lower = rowLower[iRow];
2349            double upper = rowUpper[iRow];
2350            double value = rowActivity[iRow];
2351            if (value>upper)
2352              oldSum += value-upper;
2353            else if (value<lower)
2354              oldSum += lower-value;
2355            value -= element[j];
2356            if (value>upper)
2357              newSum += value-upper;
2358            else if (value<lower)
2359              newSum += lower-value;
2360          }
2361          improvementDown = oldSum-newSum;
2362        }
2363        double theta;
2364        if (improvementUp>=improvementDown) {
2365          theta=1.0;
2366        } else {
2367          improvementUp=improvementDown;
2368          theta=-1.0;
2369        }
2370        if (improvementUp>1.0e-8&&fabs(theta)>1.0e-8) {
2371          // Could move
2372          double oldSum=0.0;
2373          double newSum=0.0;
2374          solution[iColumn] += theta;
2375          for (CoinBigIndex j=columnStart[iColumn];
2376               j<columnStart[iColumn]+columnLength[iColumn];j++) {
2377            int iRow = row[j];
2378            double lower = rowLower[iRow];
2379            double upper = rowUpper[iRow];
2380            double value = rowActivity[iRow];
2381            if (value>upper)
2382              oldSum += value-upper;
2383            else if (value<lower)
2384              oldSum += lower-value;
2385            value += theta*element[j];
2386            rowActivity[iRow]=value;
2387            if (value>upper)
2388              newSum += value-upper;
2389            else if (value<lower)
2390              newSum += lower-value;
2391          }
2392          assert (newSum<=oldSum);
2393          sumInfeasibility += newSum-oldSum;
2394        }
2395      }
2396#else
2397      int bestColumn=-1;
2398      double bestImprovement=primalTolerance;
2399      double theta=0.0;
2400      for (i=0;i<numberIntegers;i++) {
2401        int iColumn = integerVariable[i];
2402        double solValue = solution[iColumn];
2403        assert (fabs(solValue-floor(solValue+0.5))<1.0e-8);
2404        double improvementUp=0.0;
2405        if (columnUpper[iColumn]>=solValue+1.0) {
2406          // can go up
2407          double oldSum=0.0;
2408          double newSum=0.0;
2409          for (CoinBigIndex j=columnStart[iColumn];
2410               j<columnStart[iColumn]+columnLength[iColumn];j++) {
2411            int iRow = row[j];
2412            double lower = rowLower[iRow];
2413            double upper = rowUpper[iRow];
2414            double value = rowActivity[iRow];
2415            if (value>upper)
2416              oldSum += value-upper;
2417            else if (value<lower)
2418              oldSum += lower-value;
2419            value += element[j];
2420            if (value>upper)
2421              newSum += value-upper;
2422            else if (value<lower)
2423              newSum += lower-value;
2424          }
2425          improvementUp = oldSum-newSum;
2426        }
2427        double improvementDown=0.0;
2428        if (columnLower[iColumn]<=solValue-1.0) {
2429          // can go down
2430          double oldSum=0.0;
2431          double newSum=0.0;
2432          for (CoinBigIndex j=columnStart[iColumn];
2433               j<columnStart[iColumn]+columnLength[iColumn];j++) {
2434            int iRow = row[j];
2435            double lower = rowLower[iRow];
2436            double upper = rowUpper[iRow];
2437            double value = rowActivity[iRow];
2438            if (value>upper)
2439              oldSum += value-upper;
2440            else if (value<lower)
2441              oldSum += lower-value;
2442            value -= element[j];
2443            if (value>upper)
2444              newSum += value-upper;
2445            else if (value<lower)
2446              newSum += lower-value;
2447          }
2448          improvementDown = oldSum-newSum;
2449        }
2450        double improvement = CoinMax(improvementUp,improvementDown);
2451        if (improvement>bestImprovement) {
2452          bestImprovement=improvement;
2453          bestColumn = iColumn;
2454          if (improvementUp>improvementDown)
2455            theta=1.0;
2456          else
2457            theta=-1.0;
2458        }
2459      }
2460      if (bestColumn>=0) {
2461        // Could move
2462        int iColumn = bestColumn;
2463        double oldSum=0.0;
2464        double newSum=0.0;
2465        solution[iColumn] += theta;
2466        for (CoinBigIndex j=columnStart[iColumn];
2467             j<columnStart[iColumn]+columnLength[iColumn];j++) {
2468          int iRow = row[j];
2469          double lower = rowLower[iRow];
2470          double upper = rowUpper[iRow];
2471          double value = rowActivity[iRow];
2472          if (value>upper)
2473            oldSum += value-upper;
2474          else if (value<lower)
2475            oldSum += lower-value;
2476          value += theta*element[j];
2477          rowActivity[iRow]=value;
2478          if (value>upper)
2479            newSum += value-upper;
2480          else if (value<lower)
2481            newSum += lower-value;
2482        }
2483        assert (newSum<=oldSum);
2484        sumInfeasibility += newSum-oldSum;
2485      }
2486#endif
2487      if(oldSum <= sumInfeasibility + primalTolerance)
2488        break; // no good
2489    }
2490  }
2491  //delete [] saveSolution;
2492#endif
2493  delete [] rowActivity;
2494  return (largestInfeasibility>primalTolerance) ? 0 : 1;
2495}
2496// Set maximum Time (default off) - also sets starttime to current
2497void 
2498CbcHeuristicFPump::setMaximumTime(double value)
2499{
2500  startTime_=CoinCpuTime();
2501  maximumTime_=value;
2502}
2503
2504# ifdef COIN_HAS_CLP
2505 
2506//#############################################################################
2507// Constructors / Destructor / Assignment
2508//#############################################################################
2509
2510//-------------------------------------------------------------------
2511// Default Constructor
2512//-------------------------------------------------------------------
2513CbcDisasterHandler::CbcDisasterHandler (CbcModel * model) 
2514  : OsiClpDisasterHandler(),
2515    cbcModel_(model)
2516{
2517  if (model) {
2518    osiModel_
2519      = dynamic_cast<OsiClpSolverInterface *> (model->solver());
2520    if (osiModel_)
2521      setSimplex(osiModel_->getModelPtr());
2522  }
2523}
2524
2525//-------------------------------------------------------------------
2526// Copy constructor
2527//-------------------------------------------------------------------
2528CbcDisasterHandler::CbcDisasterHandler (const CbcDisasterHandler & rhs) 
2529  : OsiClpDisasterHandler(rhs),
2530    cbcModel_(rhs.cbcModel_)
2531{ 
2532}
2533
2534
2535//-------------------------------------------------------------------
2536// Destructor
2537//-------------------------------------------------------------------
2538CbcDisasterHandler::~CbcDisasterHandler ()
2539{
2540}
2541
2542//----------------------------------------------------------------
2543// Assignment operator
2544//-------------------------------------------------------------------
2545CbcDisasterHandler &
2546CbcDisasterHandler::operator=(const CbcDisasterHandler& rhs)
2547{
2548  if (this != &rhs) {
2549    OsiClpDisasterHandler::operator=(rhs);
2550    cbcModel_ = rhs.cbcModel_;
2551  }
2552  return *this;
2553}
2554//-------------------------------------------------------------------
2555// Clone
2556//-------------------------------------------------------------------
2557ClpDisasterHandler * CbcDisasterHandler::clone() const
2558{
2559  return new CbcDisasterHandler(*this);
2560}
2561// Type of disaster 0 can fix, 1 abort
2562int 
2563CbcDisasterHandler::typeOfDisaster()
2564{
2565  if (!cbcModel_->parentModel()&&
2566      (cbcModel_->specialOptions()&2048)==0) {
2567    return 0;
2568  } else {
2569    if (cbcModel_->parentModel())
2570      cbcModel_->setMaximumNodes(0);
2571    return 1;
2572  }
2573}
2574/* set model. */
2575void 
2576CbcDisasterHandler::setCbcModel(CbcModel * model)
2577{
2578  cbcModel_ = model;
2579  if (model) {
2580    osiModel_
2581      = dynamic_cast<OsiClpSolverInterface *> (model->solver());
2582    if (osiModel_)
2583      setSimplex(osiModel_->getModelPtr());
2584    else
2585      setSimplex(NULL);
2586  }
2587}
2588#endif
2589 
Note: See TracBrowser for help on using the repository browser.