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

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

changes for heuristic clone

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