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

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

fixes

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