source: stable/2.1/Cbc/src/CbcHeuristicFPump.cpp @ 942

Last change on this file since 942 was 942, checked in by forrest, 12 years ago

deal with cutoffs after solutions

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 45.5 KB
Line 
1// Copyright (C) 2004, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <cassert>
8#include <cstdlib>
9#include <cmath>
10#include <cfloat>
11
12#include "OsiSolverInterface.hpp"
13#include "CbcModel.hpp"
14#include "CbcMessage.hpp"
15#include "CbcHeuristicFPump.hpp"
16#include "CbcBranchActual.hpp"
17#include "CoinHelperFunctions.hpp"
18#include "CoinWarmStartBasis.hpp"
19#include "CoinTime.hpp"
20#include "CbcEventHandler.hpp"
21
22
23// Default Constructor
24CbcHeuristicFPump::CbcHeuristicFPump() 
25  :CbcHeuristic(),
26   startTime_(0.0),
27   maximumTime_(0.0),
28   fakeCutoff_(COIN_DBL_MAX),
29   absoluteIncrement_(0.0),
30   relativeIncrement_(0.0),
31   defaultRounding_(0.49999),
32   initialWeight_(0.0),
33   weightFactor_(0.1),
34   artificialCost_(COIN_DBL_MAX),
35   maximumPasses_(100),
36   maximumRetries_(1),
37   accumulate_(0),
38   fixOnReducedCosts_(1),
39   roundExpensive_(false)
40{
41  setWhen(1);
42}
43
44// Constructor from model
45CbcHeuristicFPump::CbcHeuristicFPump(CbcModel & model,
46                                     double downValue,bool roundExpensive)
47  :CbcHeuristic(model),
48   startTime_(0.0),
49   maximumTime_(0.0),
50   fakeCutoff_(COIN_DBL_MAX),
51   absoluteIncrement_(0.0),
52   relativeIncrement_(0.0),
53   defaultRounding_(downValue),
54   initialWeight_(0.0),
55   weightFactor_(0.1),
56   artificialCost_(COIN_DBL_MAX),
57   maximumPasses_(100),
58   maximumRetries_(1),
59   accumulate_(0),
60   fixOnReducedCosts_(1),
61   roundExpensive_(roundExpensive)
62{
63  setWhen(1);
64}
65
66// Destructor
67CbcHeuristicFPump::~CbcHeuristicFPump ()
68{
69}
70
71// Clone
72CbcHeuristic *
73CbcHeuristicFPump::clone() const
74{
75  return new CbcHeuristicFPump(*this);
76}
77// Create C++ lines to get to current state
78void 
79CbcHeuristicFPump::generateCpp( FILE * fp) 
80{
81  CbcHeuristicFPump other;
82  fprintf(fp,"0#include \"CbcHeuristicFPump.hpp\"\n");
83  fprintf(fp,"3  CbcHeuristicFPump heuristicFPump(*cbcModel);\n");
84  CbcHeuristic::generateCpp(fp,"heuristicFPump");
85  if (maximumPasses_!=other.maximumPasses_)
86    fprintf(fp,"3  heuristicFPump.setMaximumPasses(%d);\n",maximumPasses_);
87  else
88    fprintf(fp,"4  heuristicFPump.setMaximumPasses(%d);\n",maximumPasses_);
89  if (maximumRetries_!=other.maximumRetries_)
90    fprintf(fp,"3  heuristicFPump.setMaximumRetries(%d);\n",maximumRetries_);
91  else
92    fprintf(fp,"4  heuristicFPump.setMaximumRetries(%d);\n",maximumRetries_);
93  if (accumulate_!=other.accumulate_)
94    fprintf(fp,"3  heuristicFPump.setAccumulate(%d);\n",accumulate_);
95  else
96    fprintf(fp,"4  heuristicFPump.setAccumulate(%d);\n",accumulate_);
97  if (fixOnReducedCosts_!=other.fixOnReducedCosts_)
98    fprintf(fp,"3  heuristicFPump.setFixOnReducedCosts(%d);\n",fixOnReducedCosts_);
99  else
100    fprintf(fp,"4  heuristicFPump.setFixOnReducedCosts(%d);\n",fixOnReducedCosts_);
101  if (maximumTime_!=other.maximumTime_)
102    fprintf(fp,"3  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
103  else
104    fprintf(fp,"4  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
105  if (fakeCutoff_!=other.fakeCutoff_)
106    fprintf(fp,"3  heuristicFPump.setFakeCutoff(%g);\n",fakeCutoff_);
107  else
108    fprintf(fp,"4  heuristicFPump.setFakeCutoff(%g);\n",fakeCutoff_);
109  if (absoluteIncrement_!=other.absoluteIncrement_)
110    fprintf(fp,"3  heuristicFPump.setAbsoluteIncrement(%g);\n",absoluteIncrement_);
111  else
112    fprintf(fp,"4  heuristicFPump.setAbsoluteIncrement(%g);\n",absoluteIncrement_);
113  if (relativeIncrement_!=other.relativeIncrement_)
114    fprintf(fp,"3  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
115  else
116    fprintf(fp,"4  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
117  if (defaultRounding_!=other.defaultRounding_)
118    fprintf(fp,"3  heuristicFPump.setDefaultRounding(%g);\n",defaultRounding_);
119  else
120    fprintf(fp,"4  heuristicFPump.setDefaultRounding(%g);\n",defaultRounding_);
121  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicFPump);\n");
122  if (initialWeight_!=other.initialWeight_)
123    fprintf(fp,"3  heuristicFPump.setInitialWeight(%g);\n",initialWeight_);
124  else
125    fprintf(fp,"4  heuristicFPump.setInitialWeight(%g);\n",initialWeight_);
126  if (weightFactor_!=other.weightFactor_)
127    fprintf(fp,"3  heuristicFPump.setWeightFactor(%g);\n",weightFactor_);
128  else
129    fprintf(fp,"4  heuristicFPump.setWeightFactor(%g);\n",weightFactor_);
130}
131
132// Copy constructor
133CbcHeuristicFPump::CbcHeuristicFPump(const CbcHeuristicFPump & rhs)
134:
135  CbcHeuristic(rhs),
136  startTime_(rhs.startTime_),
137  maximumTime_(rhs.maximumTime_),
138  fakeCutoff_(rhs.fakeCutoff_),
139  absoluteIncrement_(rhs.absoluteIncrement_),
140  relativeIncrement_(rhs.relativeIncrement_),
141  defaultRounding_(rhs.defaultRounding_),
142  initialWeight_(rhs.initialWeight_),
143  weightFactor_(rhs.weightFactor_),
144  artificialCost_(rhs.artificialCost_),
145  maximumPasses_(rhs.maximumPasses_),
146  maximumRetries_(rhs.maximumRetries_),
147  accumulate_(rhs.accumulate_),
148  fixOnReducedCosts_(rhs.fixOnReducedCosts_),
149  roundExpensive_(rhs.roundExpensive_)
150{
151}
152
153// Assignment operator
154CbcHeuristicFPump & 
155CbcHeuristicFPump::operator=( const CbcHeuristicFPump& rhs)
156{
157  if (this!=&rhs) {
158    CbcHeuristic::operator=(rhs);
159    startTime_ = rhs.startTime_;
160    maximumTime_ = rhs.maximumTime_;
161    fakeCutoff_ = rhs.fakeCutoff_;
162    absoluteIncrement_ = rhs.absoluteIncrement_;
163    relativeIncrement_ = rhs.relativeIncrement_;
164    defaultRounding_ = rhs.defaultRounding_;
165    initialWeight_ = rhs.initialWeight_;
166    weightFactor_ = rhs.weightFactor_;
167    artificialCost_ = rhs.artificialCost_;
168    maximumPasses_ = rhs.maximumPasses_;
169    maximumRetries_ = rhs.maximumRetries_;
170    accumulate_ = rhs.accumulate_;
171    fixOnReducedCosts_ = rhs.fixOnReducedCosts_;
172    roundExpensive_ = rhs.roundExpensive_;
173  }
174  return *this;
175}
176
177// Resets stuff if model changes
178void 
179CbcHeuristicFPump::resetModel(CbcModel * model)
180{
181}
182
183/**************************BEGIN MAIN PROCEDURE ***********************************/
184
185// See if feasibility pump will give better solution
186// Sets value of solution
187// Returns 1 if solution, 0 if not
188int
189CbcHeuristicFPump::solution(double & solutionValue,
190                         double * betterSolution)
191{
192  double incomingObjective = solutionValue;
193#define LEN_PRINT 200
194  char pumpPrint[LEN_PRINT];
195  pumpPrint[0]='\0';
196  if (!when()||(when()==1&&model_->phase()!=1))
197    return 0; // switched off
198  // See if at root node
199  bool atRoot = model_->getNodeCount()==0;
200  int passNumber = model_->getCurrentPassNumber();
201  // just do once
202  if (!atRoot||passNumber!=1)
203    return 0;
204  // probably a good idea
205  //if (model_->getSolutionCount()) return 0;
206  // loop round doing repeated pumps
207  double cutoff;
208  model_->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
209  double direction = model_->solver()->getObjSense();
210  cutoff *= direction;
211  double firstCutoff = fabs(cutoff);
212  cutoff = CoinMin(cutoff,solutionValue);
213  // check plausible and space for rounded solution
214  int numberColumns = model_->getNumCols();
215  int numberIntegers = model_->numberIntegers();
216  const int * integerVariableOrig = model_->integerVariable();
217
218  // 1. initially check 0-1
219  int i,j;
220  int general=0;
221  int * integerVariable = new int[numberIntegers];
222  const double * lower = model_->solver()->getColLower();
223  const double * upper = model_->solver()->getColUpper();
224  bool doGeneral = (accumulate_&(4+8))!=0;
225  j=0;
226  for (i=0;i<numberIntegers;i++) {
227    int iColumn = integerVariableOrig[i];
228#ifndef NDEBUG
229    const OsiObject * object = model_->object(i);
230    const CbcSimpleInteger * integerObject = 
231      dynamic_cast<const  CbcSimpleInteger *> (object);
232    const OsiSimpleInteger * integerObject2 = 
233      dynamic_cast<const  OsiSimpleInteger *> (object);
234    assert(integerObject||integerObject2);
235#endif
236    if (upper[iColumn]-lower[iColumn]>1.000001) {
237      general++;
238      if (doGeneral)
239        integerVariable[j++]=iColumn;
240    } else {
241      integerVariable[j++]=iColumn;
242    }
243  }
244  if (general*3>2*numberIntegers&&!doGeneral) {
245    delete [] integerVariable;
246    return 0;
247  } else if ((accumulate_&8)==0) {
248    doGeneral=false;
249  }
250  if (!general)
251    doGeneral=false;
252  // For solution closest to feasible if none found
253  int * closestSolution = general ? NULL : new int[numberIntegers];
254  double closestObjectiveValue = COIN_DBL_MAX;
255 
256  int numberIntegersOrig = numberIntegers;
257  numberIntegers = j;
258  double * newSolution = new double [numberColumns];
259  double newSolutionValue=COIN_DBL_MAX;
260  bool solutionFound=false;
261  char * usedColumn = NULL;
262  double * lastSolution=NULL;
263  int fixContinuous=0;
264  bool fixInternal=false;
265  bool allSlack=false;
266  if (when_>=21&&when_<=25) {
267    when_ -= 10;
268    allSlack=true;
269  }
270  double time1 = CoinCpuTime();
271  model_->solver()->resolve();
272  if (!model_->solver()->isProvenOptimal()) {
273    // presumably max time or some such
274    return 0;
275  }
276  if (cutoff<1.0e50&&false) {
277    // Fix on djs
278    double direction = model_->solver()->getObjSense() ;
279    double gap = cutoff - model_->solver()->getObjValue()*direction ;
280    double tolerance;
281    model_->solver()->getDblParam(OsiDualTolerance,tolerance) ;
282    if (gap>0.0) {
283      gap += 100.0*tolerance;
284      int nFix=model_->solver()->reducedCostFix(gap);
285      printf("dj fixing fixed %d variables\n",nFix);
286    }
287  }
288  CoinWarmStartBasis saveBasis;
289  CoinWarmStartBasis * basis =
290    dynamic_cast<CoinWarmStartBasis *>(model_->solver()->getWarmStart()) ;
291  if (basis) {
292    saveBasis = * basis;
293    delete basis;
294  }
295  double continuousObjectiveValue = model_->solver()->getObjValue()*model_->solver()->getObjSense();
296  double * firstPerturbedObjective = NULL;
297  double * firstPerturbedSolution = NULL;
298  if (when_>=11&&when_<=15) {
299    fixInternal = when_ >11&&when_<15;
300    if (when_<13)
301      fixContinuous = 0;
302    else if (when_!=14)
303      fixContinuous=1;
304    else
305      fixContinuous=2;
306    when_=1;
307    usedColumn = new char [numberColumns];
308    memset(usedColumn,0,numberColumns);
309    lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(),numberColumns);
310  }
311  int finalReturnCode=0;
312  int totalNumberPasses=0;
313  int numberTries=0;
314  CoinWarmStartBasis bestBasis;
315  bool exitAll=false;
316  double saveBestObjective = model_->getMinimizationObjValue();
317  int numberSolutions=0;
318  OsiSolverInterface * solver = NULL;
319  double artificialFactor = 0.00001;
320  // also try rounding!
321  double * roundingSolution = new double[numberColumns];
322  double roundingObjective = solutionValue;
323  while (!exitAll) {
324    int numberPasses=0;
325    artificialFactor *= 10.0;
326    numberTries++;
327    // Clone solver - otherwise annoys root node computations
328    solver = model_->solver()->clone();
329    if (CoinMin(fakeCutoff_,cutoff)<1.0e50) {
330      // Fix on djs
331      double direction = solver->getObjSense() ;
332      double gap = CoinMin(fakeCutoff_,cutoff) - solver->getObjValue()*direction ;
333      double tolerance;
334      solver->getDblParam(OsiDualTolerance,tolerance) ;
335      if (gap>0.0&&(fixOnReducedCosts_==1||(numberTries==1&&fixOnReducedCosts_==2))) {
336        gap += 100.0*tolerance;
337        int nFix=solver->reducedCostFix(gap);
338        if (nFix) {
339          sprintf(pumpPrint,"Reduced cost fixing fixed %d variables on major pass %d",nFix,numberTries);
340          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
341            << pumpPrint
342            <<CoinMessageEol;
343          //pumpPrint[0]='\0';
344        }
345      }
346    }
347    // if cutoff exists then add constraint
348    bool useCutoff = (fabs(cutoff)<1.0e20&&(fakeCutoff_!=COIN_DBL_MAX||numberTries>1));
349    // but there may be a close one
350    if (firstCutoff<2.0*solutionValue&&numberTries==1) 
351      useCutoff=true;
352    if (useCutoff) {
353      cutoff = CoinMin(cutoff,fakeCutoff_);
354      const double * objective = solver->getObjCoefficients();
355      int numberColumns = solver->getNumCols();
356      int * which = new int[numberColumns];
357      double * els = new double[numberColumns];
358      int nel=0;
359      for (int i=0;i<numberColumns;i++) {
360        double value = objective[i];
361        if (value) {
362          which[nel]=i;
363          els[nel++] = direction*value;
364        }
365      }
366      double offset;
367      solver->getDblParam(OsiObjOffset,offset);
368#ifdef COIN_DEVELOP
369      if (offset)
370        printf("CbcHeuristicFPump obj offset %g\n",offset);
371#endif
372      solver->addRow(nel,which,els,-COIN_DBL_MAX,cutoff+offset*direction);
373      delete [] which;
374      delete [] els;
375      bool takeHint;
376      OsiHintStrength strength;
377      solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
378      solver->setHintParam(OsiDoDualInResolve,true,OsiHintDo);
379      solver->resolve();
380      solver->setHintParam(OsiDoDualInResolve,takeHint,strength);
381      if (!solver->isProvenOptimal()) {
382        // presumably max time or some such
383        exitAll=true;
384        break;
385      }
386    }
387    solver->setDblParam(OsiDualObjectiveLimit,1.0e50);
388    solver->resolve();
389    // Solver may not be feasible
390    if (!solver->isProvenOptimal()) {
391      exitAll=true;
392      break;
393    }
394    const double * lower = solver->getColLower();
395    const double * upper = solver->getColUpper();
396    const double * solution = solver->getColSolution();
397    if (lastSolution)
398      memcpy(lastSolution,solution,numberColumns*sizeof(double));
399    double primalTolerance;
400    solver->getDblParam(OsiPrimalTolerance,primalTolerance);
401   
402   
403    // 2 space for last rounded solutions
404#define NUMBER_OLD 4
405    double ** oldSolution = new double * [NUMBER_OLD];
406    for (j=0;j<NUMBER_OLD;j++) {
407      oldSolution[j]= new double[numberColumns];
408      for (i=0;i<numberColumns;i++) oldSolution[j][i]=-COIN_DBL_MAX;
409    }
410   
411    // 3. Replace objective with an initial 0-valued objective
412    double * saveObjective = new double [numberColumns];
413    memcpy(saveObjective,solver->getObjCoefficients(),numberColumns*sizeof(double));
414    for (i=0;i<numberColumns;i++) {
415      solver->setObjCoeff(i,0.0);
416    }
417    bool finished=false;
418    double direction = solver->getObjSense();
419    int returnCode=0;
420    bool takeHint;
421    OsiHintStrength strength;
422    solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
423    solver->setHintParam(OsiDoDualInResolve,false);
424    //solver->messageHandler()->setLogLevel(0);
425   
426    // 4. Save objective offset so we can see progress
427    double saveOffset;
428    solver->getDblParam(OsiObjOffset,saveOffset);
429    // Get amount for original objective
430    double scaleFactor = 0.0;
431#ifdef COIN_DEVELOP
432    double largestCost=0.0;
433    int nArtificial=0;
434#endif
435    for (i=0;i<numberColumns;i++) {
436      double value = saveObjective[i];
437      scaleFactor += value*value;
438#ifdef COIN_DEVELOP
439      largestCost=CoinMax(largestCost,fabs(value));
440      if (value*direction>=artificialCost_)
441        nArtificial++;
442#endif
443    }
444    if (scaleFactor)
445      scaleFactor = (initialWeight_*sqrt((double) numberIntegers))/sqrt(scaleFactor);
446#ifdef COIN_DEVELOP
447    if (scaleFactor)
448      printf("Using %g fraction of original objective - largest %g - %d artificials\n",scaleFactor,
449             largestCost,nArtificial);
450#endif
451    // 5. MAIN WHILE LOOP
452    //bool newLineNeeded=false;
453    while (!finished) {
454      double newTrueSolutionValue=0.0;
455      double newSumInfeas=0.0;
456      returnCode=0;
457      if (model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
458        exitAll=true;
459        break;
460      }
461      // see what changed
462      if (usedColumn) {
463        for (i=0;i<numberColumns;i++) {
464          if (fabs(solution[i]-lastSolution[i])>1.0e-8) 
465            usedColumn[i]=1;
466          lastSolution[i]=solution[i];
467        }
468      }
469      if (numberPasses>=maximumPasses_) {
470        break;
471      }
472      if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
473      memcpy(newSolution,solution,numberColumns*sizeof(double));
474      int flip;
475      if (numberPasses==0&&false) {
476        // always use same seed
477        randomNumberGenerator_.setSeed(987654321);
478      }
479      returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
480                          pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
481      if (numberPasses==0&&false) {
482        // Make sure random will be different
483        for (i=1;i<numberTries;i++)
484          randomNumberGenerator_.randomDouble();
485      }
486      numberPasses++;
487      if (returnCode) {
488        // SOLUTION IS INTEGER
489        // Put back correct objective
490        for (i=0;i<numberColumns;i++)
491          solver->setObjCoeff(i,saveObjective[i]);
492
493        // solution - but may not be better
494        // Compute using dot product
495        solver->setDblParam(OsiObjOffset,saveOffset);
496        newSolutionValue = -saveOffset;
497        for (  i=0 ; i<numberColumns ; i++ )
498          newSolutionValue += saveObjective[i]*newSolution[i];
499        newSolutionValue *= direction;
500        sprintf(pumpPrint,"Solution found of %g",newSolutionValue);
501        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
502          << pumpPrint
503          <<CoinMessageEol;
504        //newLineNeeded=false;
505        if (newSolutionValue<solutionValue) {
506          double saveValue = solutionValue;
507          if (!doGeneral) {
508            int numberLeft=0;
509            for (i=0;i<numberIntegersOrig;i++) {
510              int iColumn = integerVariableOrig[i];
511              double value = floor(newSolution[iColumn]+0.5);
512              if(solver->isBinary(iColumn)) {
513                solver->setColLower(iColumn,value);
514                solver->setColUpper(iColumn,value);
515              } else {
516                if (fabs(value-newSolution[iColumn])>1.0e-7) 
517                  numberLeft++;
518              }
519            }
520            if (numberLeft) {
521              returnCode = smallBranchAndBound(solver,numberNodes_,newSolution,newSolutionValue,
522                                               solutionValue,"CbcHeuristicFpump");
523              if (returnCode<0) {
524                if (returnCode==-2)
525                  exitAll=true;
526                returnCode=0; // returned on size or event
527              }
528              if ((returnCode&2)!=0) {
529                // could add cut
530                returnCode &= ~2;
531              }
532              if (returnCode!=1)
533                newSolutionValue=saveValue;
534            }
535          }
536          if (returnCode&&newSolutionValue<saveValue) {
537            memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
538            solutionFound=true;
539            CoinWarmStartBasis * basis =
540              dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
541            if (basis) {
542              bestBasis = * basis;
543              delete basis;
544              CbcEventHandler * handler = model_->getEventHandler();
545              if (handler) {
546                double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
547                double saveObjectiveValue = model_->getMinimizationObjValue();
548                model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
549                int action = handler->event(CbcEventHandler::heuristicSolution);
550                if (saveOldSolution) {
551                  model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
552                  delete [] saveOldSolution;
553                }
554                if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
555                  exitAll=true; // exit
556                  break;
557                }
558              }
559            }
560            if ((accumulate_&1)!=0) {
561              model_->incrementUsed(betterSolution); // for local search
562              numberSolutions++;
563            }
564            solutionValue=newSolutionValue;
565            solutionFound=true;
566            if (general&&saveValue!=newSolutionValue) {
567              sprintf(pumpPrint,"Cleaned solution of %g",solutionValue);
568              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
569                << pumpPrint
570                <<CoinMessageEol;
571            }
572          } else {
573            sprintf(pumpPrint,"Mini branch and bound could not fix general integers");
574            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
575              << pumpPrint
576              <<CoinMessageEol;
577          }
578        } else {
579          sprintf(pumpPrint,"After further testing solution no better than previous of %g",solutionValue);
580
581          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
582            << pumpPrint
583            <<CoinMessageEol;
584          //newLineNeeded=false;
585          returnCode=0;
586        }     
587        break;
588      } else {
589        // SOLUTION IS not INTEGER
590        // 1. check for loop
591        bool matched;
592        for (int k = NUMBER_OLD-1; k > 0; k--) {
593          double * b = oldSolution[k];
594          matched = true;
595          for (i = 0; i <numberIntegers; i++) {
596            int iColumn = integerVariable[i];
597            if (newSolution[iColumn]!=b[iColumn]) {
598              matched=false;
599              break;
600            }
601          }
602          if (matched) break;
603        }
604        int numberPerturbed=0;
605        if (matched || numberPasses%100 == 0) {
606          // perturbation
607          //sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
608          //newLineNeeded=true;
609          double factorX[10]={0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0};
610          double factor=1.0;
611          double target=-1.0;
612          double * randomX = new double [numberIntegers];
613          for (i=0;i<numberIntegers;i++) 
614            randomX[i] = max(0.0,randomNumberGenerator_.randomDouble()-0.3);
615          for (int k=0;k<10;k++) {
616#ifdef COIN_DEVELOP_x
617            printf("kpass %d\n",k);
618#endif
619            int numberX[10]={0,0,0,0,0,0,0,0,0,0};
620            for (i=0;i<numberIntegers;i++) {
621              int iColumn = integerVariable[i];
622              double value = randomX[i];
623              double difference = fabs(solution[iColumn]-newSolution[iColumn]);
624              for (int j=0;j<10;j++) {
625                if (difference+value*factorX[j]>0.5) 
626                  numberX[j]++;
627              }
628            }
629            if (target<0.0) {
630              if (numberX[9]<=200)
631                break; // not very many changes
632              target=CoinMax(200.0,CoinMin(0.05*numberX[9],1000.0));
633            }
634            int iX=-1;
635            int iBand=-1;
636            for (i=0;i<10;i++) {
637#ifdef COIN_DEVELOP_x
638              printf("** %d changed at %g\n",numberX[i],factorX[i]);
639#endif
640              if (numberX[i]>=target&&numberX[i]<2.0*target&&iX<0)
641                iX=i;
642              if (iBand<0&&numberX[i]>target) {
643                iBand=i;
644                factor=factorX[i];
645              }
646            }
647            if (iX>=0) {
648              factor=factorX[iX];
649              break;
650            } else {
651              assert (iBand>=0);
652              double hi = factor;
653              double lo = (iBand>0) ? factorX[iBand-1] : 0.0;
654              double diff = (hi-lo)/9.0;
655              for (i=0;i<10;i++) {
656                factorX[i]=lo;
657                lo += diff;
658              }
659            }
660          }
661          for (i=0;i<numberIntegers;i++) {
662            int iColumn = integerVariable[i];
663            double value = randomX[i];
664            double difference = fabs(solution[iColumn]-newSolution[iColumn]);
665            if (difference+value*factor>0.5) {
666              numberPerturbed++;
667              if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
668                newSolution[iColumn] += 1.0;
669              } else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
670                newSolution[iColumn] -= 1.0;
671              } else {
672                // general integer
673                if (difference+value>0.75)
674                  newSolution[iColumn] += 1.0;
675                else
676                  newSolution[iColumn] -= 1.0;
677              }
678            }
679          }
680          delete [] randomX;
681        } else {
682          for (j=NUMBER_OLD-1;j>0;j--) {
683            for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
684          }
685          for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
686        }
687       
688        // 2. update the objective function based on the new rounded solution
689        double offset=0.0;
690        double costValue = (1.0-scaleFactor)*solver->getObjSense();
691        int numberChanged=0;
692        const double * oldObjective = solver->getObjCoefficients();
693        for (i=0;i<numberColumns;i++) {
694          // below so we can keep original code and allow for objective
695          int iColumn = i;
696          // Special code for "artificials"
697          if (direction*saveObjective[iColumn]>=artificialCost_) {
698            //solver->setObjCoeff(iColumn,scaleFactor*saveObjective[iColumn]);
699            solver->setObjCoeff(iColumn,(artificialFactor*saveObjective[iColumn])/artificialCost_);
700          }
701          if(!solver->isBinary(iColumn)&&!doGeneral)
702            continue;
703          // deal with fixed variables (i.e., upper=lower)
704          if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance||!solver->isInteger(iColumn)) {
705            //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
706            //else                                       solver->setObjCoeff(iColumn,costValue);
707            continue;
708          }
709          double newValue=0.0;
710          if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
711            newValue = costValue+scaleFactor*saveObjective[iColumn];
712          } else {
713            if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
714              newValue = -costValue+scaleFactor*saveObjective[iColumn];
715            }
716          }
717          if (newValue!=oldObjective[iColumn])
718            numberChanged++;
719          solver->setObjCoeff(iColumn,newValue);
720          offset += costValue*newSolution[iColumn];
721        }
722        solver->setDblParam(OsiObjOffset,-offset);
723        if (!general&&false) {
724          // Solve in two goes - first keep satisfied ones fixed
725          double * saveLower = new double [numberIntegers];
726          double * saveUpper = new double [numberIntegers];
727          for (i=0;i<numberIntegers;i++) {
728            int iColumn = integerVariable[i];
729            saveLower[i]=COIN_DBL_MAX;
730            saveUpper[i]=-COIN_DBL_MAX;
731            if (solution[iColumn]<lower[iColumn]+primalTolerance) {
732              saveUpper[i]=upper[iColumn];
733              solver->setColUpper(iColumn,lower[iColumn]);
734            } else if (solution[iColumn]>upper[iColumn]-primalTolerance) {
735              saveLower[i]=lower[iColumn];
736              solver->setColLower(iColumn,upper[iColumn]);
737            }
738          }
739          solver->resolve();
740          if (!solver->isProvenOptimal()) {
741            // presumably max time or some such
742            exitAll=true;
743            break;
744          }
745          for (i=0;i<numberIntegers;i++) {
746            int iColumn = integerVariable[i];
747            if (saveLower[i]!=COIN_DBL_MAX)
748              solver->setColLower(iColumn,saveLower[i]);
749            if (saveUpper[i]!=-COIN_DBL_MAX)
750              solver->setColUpper(iColumn,saveUpper[i]);
751            saveUpper[i]=-COIN_DBL_MAX;
752          }
753          memcpy(newSolution,solution,numberColumns*sizeof(double));
754          int flip;
755          returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
756                              pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
757          numberPasses++;
758          if (returnCode) {
759            // solution - but may not be better
760            // Compute using dot product
761            double newSolutionValue = -saveOffset;
762            for (  i=0 ; i<numberColumns ; i++ )
763              newSolutionValue += saveObjective[i]*newSolution[i];
764            newSolutionValue *= direction;
765            sprintf(pumpPrint,"Intermediate solution found of %g",newSolutionValue);
766            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
767              << pumpPrint
768              <<CoinMessageEol;
769            if (newSolutionValue<solutionValue) {
770              memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
771              CoinWarmStartBasis * basis =
772                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
773              solutionFound=true;
774              if (basis) {
775                bestBasis = * basis;
776                delete basis;
777                CbcEventHandler * handler = model_->getEventHandler();
778                if (handler) {
779                  double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
780                  double saveObjectiveValue = model_->getMinimizationObjValue();
781                  model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
782                  int action = handler->event(CbcEventHandler::heuristicSolution);
783                  if (saveOldSolution) {
784                    model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
785                    delete [] saveOldSolution;
786                  }
787                  if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
788                    exitAll=true; // exit
789                    break;
790                  }
791                }
792              }
793              if ((accumulate_&1)!=0) {
794                model_->incrementUsed(betterSolution); // for local search
795                numberSolutions++;
796              }
797              solutionValue=newSolutionValue;
798              solutionFound=true;
799            } else {
800              returnCode=0;
801            }
802          }     
803        }
804        if (!doGeneral) {
805          // faster to do from all slack!!!!
806          if (allSlack) {
807            CoinWarmStartBasis dummy;
808            solver->setWarmStart(&dummy);
809          }
810#ifdef COIN_DEVELOP
811          printf("%d perturbed out of %d columns (%d changed)\n",numberPerturbed,numberColumns,numberChanged);
812#endif
813          bool takeHint;
814          OsiHintStrength strength;
815          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
816          if (numberPerturbed>numberColumns)
817            solver->setHintParam(OsiDoDualInResolve,true); // dual may be better
818          if (numberTries>1&&numberPasses==1&&false) {
819            // use basis from first time
820            solver->setWarmStart(&saveBasis);
821            // and objective
822            solver->setObjective(firstPerturbedObjective);
823            // and solution
824            solver->setColSolution(firstPerturbedSolution);
825          }
826          solver->resolve();
827          if (!solver->isProvenOptimal()) {
828            // presumably max time or some such
829            exitAll=true;
830            break;
831          }
832          if (numberTries==1&&numberPasses==1&&false) {
833            // save basis
834            CoinWarmStartBasis * basis =
835              dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
836            if (basis) {
837              saveBasis = * basis;
838              delete basis;
839            }
840            firstPerturbedObjective = CoinCopyOfArray(solver->getObjCoefficients(),numberColumns);
841            firstPerturbedSolution = CoinCopyOfArray(solver->getColSolution(),numberColumns);
842          }
843          solver->setHintParam(OsiDoDualInResolve,takeHint);
844          newTrueSolutionValue = -saveOffset;
845          newSumInfeas=0.0;
846          {
847            const double * newSolution = solver->getColSolution();
848            for (  i=0 ; i<numberColumns ; i++ ) {
849              if (solver->isInteger(i)) {
850                double value = newSolution[i];
851                double nearest = floor(value+0.5);
852                newSumInfeas += fabs(value-nearest);
853              }
854              newTrueSolutionValue += saveObjective[i]*newSolution[i];
855            }
856            newTrueSolutionValue *= direction;
857            OsiSolverInterface * saveSolver = model_->swapSolver(solver);
858            CbcRounding heuristic1(*model_);
859            heuristic1.setHeuristicName("rounding in feaspump!");
860            heuristic1.setWhen(1);
861            double testObjectiveValue = CoinMin(solutionValue,roundingObjective);
862            int returnCode = heuristic1.solution(testObjectiveValue,roundingSolution,newTrueSolutionValue) ;
863            if (returnCode==1) {
864              assert(testObjectiveValue < CoinMin(solutionValue,roundingObjective));
865              roundingObjective = testObjectiveValue;
866            }
867            model_->swapSolver(saveSolver);
868          }
869          if (!solver->isProvenOptimal()) {
870            // presumably max time or some such
871            exitAll=true;
872            break;
873          }
874          // in case very dubious solver
875          lower = solver->getColLower();
876          upper = solver->getColUpper();
877          solution = solver->getColSolution();
878        } else {
879          int * addStart = new int[2*general+1];
880          int * addIndex = new int[4*general];
881          double * addElement = new double[4*general];
882          double * addLower = new double[2*general];
883          double * addUpper = new double[2*general];
884          double * obj = new double[general];
885          int nAdd=0;
886          for (i=0;i<numberIntegers;i++) {
887            int iColumn = integerVariable[i];
888            if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
889                newSolution[iColumn]<upper[iColumn]-primalTolerance) {
890              obj[nAdd]=1.0;
891              addLower[nAdd]=0.0;
892              addUpper[nAdd]=COIN_DBL_MAX;
893              nAdd++;
894            }
895          }
896          OsiSolverInterface * solver2 = solver;
897          if (nAdd) {
898            CoinZeroN(addStart,nAdd+1);
899            solver2 = solver->clone();
900            solver2->addCols(nAdd,addStart,NULL,NULL,addLower,addUpper,obj);
901            // feasible solution
902            double * sol = new double[nAdd+numberColumns];
903            memcpy(sol,solution,numberColumns*sizeof(double));
904            // now rows
905            int nAdd=0;
906            int nEl=0;
907            int nAddRow=0;
908            for (i=0;i<numberIntegers;i++) {
909              int iColumn = integerVariable[i];
910              if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
911                  newSolution[iColumn]<upper[iColumn]-primalTolerance) {
912                addLower[nAddRow]=-newSolution[iColumn];;
913                addUpper[nAddRow]=COIN_DBL_MAX;
914                addIndex[nEl] = iColumn;
915                addElement[nEl++]=-1.0;
916                addIndex[nEl] = numberColumns+nAdd;
917                addElement[nEl++]=1.0;
918                nAddRow++;
919                addStart[nAddRow]=nEl;
920                addLower[nAddRow]=newSolution[iColumn];;
921                addUpper[nAddRow]=COIN_DBL_MAX;
922                addIndex[nEl] = iColumn;
923                addElement[nEl++]=1.0;
924                addIndex[nEl] = numberColumns+nAdd;
925                addElement[nEl++]=1.0;
926                nAddRow++;
927                addStart[nAddRow]=nEl;
928                sol[nAdd+numberColumns] = fabs(sol[iColumn]-newSolution[iColumn]);
929                nAdd++;
930              }
931            }
932            solver2->setColSolution(sol);
933            delete [] sol;
934            solver2->addRows(nAddRow,addStart,addIndex,addElement,addLower,addUpper);
935          }
936          delete [] addStart;
937          delete [] addIndex;
938          delete [] addElement;
939          delete [] addLower;
940          delete [] addUpper;
941          delete [] obj;
942          solver2->resolve();
943          if (!solver2->isProvenOptimal()) {
944            // presumably max time or some such
945            exitAll=true;
946            break;
947          }
948          //assert (solver2->isProvenOptimal());
949          if (nAdd) {
950            solver->setColSolution(solver2->getColSolution());
951            delete solver2;
952          }
953        }
954        if (solver->getNumRows()<3000)
955          sprintf(pumpPrint,"Pass %3d: suminf. %10.5f obj. %g iterations %d", numberPasses+totalNumberPasses,
956                  newSumInfeas,newTrueSolutionValue,solver->getIterationCount());
957        else
958          sprintf(pumpPrint,"Pass %3d: (%.2f seconds) suminf. %10.5f obj. %g iterations %d", numberPasses+totalNumberPasses,
959                  model_->getCurrentSeconds(),newSumInfeas,newTrueSolutionValue,solver->getIterationCount());
960        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
961          << pumpPrint
962          <<CoinMessageEol;
963        if (closestSolution&&solver->getObjValue()<closestObjectiveValue) {
964          int i;
965          const double * objective = solver->getObjCoefficients();
966          for (i=0;i<numberIntegersOrig;i++) {
967            int iColumn=integerVariableOrig[i];
968            if (objective[iColumn]>0.0)
969              closestSolution[i]=0;
970            else
971              closestSolution[i]=1;
972          }
973          closestObjectiveValue = solver->getObjValue();
974        }
975        //newLineNeeded=true;
976       
977      }
978      // reduce scale factor
979      scaleFactor *= weightFactor_;
980    } // END WHILE
981    // see if rounding worked!
982    if (roundingObjective<solutionValue) {
983      sprintf(pumpPrint,"Rounding solution of %g is better than previous of %g !\n",
984              roundingObjective,solutionValue);
985      solutionValue=roundingObjective;
986      memcpy(betterSolution,roundingSolution,numberColumns*sizeof(double));
987      solutionFound=true;
988    }
989    if (!solutionFound) { 
990      sprintf(pumpPrint,"No solution found this major pass");
991      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
992        << pumpPrint
993        <<CoinMessageEol;
994    }
995    //}
996    delete solver;
997    solver=NULL;
998    for ( j=0;j<NUMBER_OLD;j++) 
999      delete [] oldSolution[j];
1000    delete [] oldSolution;
1001    delete [] saveObjective;
1002    if (usedColumn&&!exitAll) {
1003      OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
1004      const double * colLower = newSolver->getColLower();
1005      const double * colUpper = newSolver->getColUpper();
1006      int i;
1007      int nFix=0;
1008      int nFixI=0;
1009      int nFixC=0;
1010      int nFixC2=0;
1011      for (i=0;i<numberIntegersOrig;i++) {
1012        int iColumn=integerVariableOrig[i];
1013        //const OsiObject * object = model_->object(i);
1014        //double originalLower;
1015        //double originalUpper;
1016        //getIntegerInformation( object,originalLower, originalUpper);
1017        //assert(colLower[iColumn]==originalLower);
1018        //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
1019        newSolver->setColLower(iColumn,colLower[iColumn]);
1020        //assert(colUpper[iColumn]==originalUpper);
1021        //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
1022        newSolver->setColUpper(iColumn,colUpper[iColumn]);
1023        if (!usedColumn[iColumn]) {
1024          double value=lastSolution[iColumn];
1025          double nearest = floor(value+0.5);
1026          if (fabs(value-nearest)<1.0e-7) {
1027            if (nearest==colLower[iColumn]) {
1028              newSolver->setColUpper(iColumn,colLower[iColumn]);
1029              nFix++;
1030            } else if (nearest==colUpper[iColumn]) {
1031              newSolver->setColLower(iColumn,colUpper[iColumn]);
1032              nFix++;
1033            } else if (fixInternal) {
1034              newSolver->setColLower(iColumn,nearest);
1035              newSolver->setColUpper(iColumn,nearest);
1036              nFix++;
1037              nFixI++;
1038            }
1039          }
1040        }
1041      }
1042      if (fixContinuous) {
1043        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1044          if (!newSolver->isInteger(iColumn)&&!usedColumn[iColumn]) {
1045            double value=lastSolution[iColumn];
1046            if (value<colLower[iColumn]+1.0e-8) {
1047              newSolver->setColUpper(iColumn,colLower[iColumn]);
1048              nFixC++;
1049            } else if (value>colUpper[iColumn]-1.0e-8) {
1050              newSolver->setColLower(iColumn,colUpper[iColumn]);
1051              nFixC++;
1052            } else if (fixContinuous==2) {
1053              newSolver->setColLower(iColumn,value);
1054              newSolver->setColUpper(iColumn,value);
1055              nFixC++;
1056              nFixC2++;
1057            }
1058          }
1059        }
1060      }
1061      newSolver->initialSolve();
1062      if (!newSolver->isProvenOptimal()) {
1063        //newSolver->writeMps("bad.mps");
1064        //assert (newSolver->isProvenOptimal());
1065        exitAll=true;
1066        break;
1067      }
1068      sprintf(pumpPrint,"Before mini branch and bound, %d integers at bound fixed and %d continuous",
1069             nFix,nFixC);
1070      if (nFixC2+nFixI!=0)
1071        sprintf(pumpPrint+strlen(pumpPrint)," of which %d were internal integer and %d internal continuous",
1072                nFixI,nFixC2);
1073      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1074        << pumpPrint
1075        <<CoinMessageEol;
1076      double saveValue = newSolutionValue;
1077      returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
1078                                       cutoff,"CbcHeuristicLocalAfterFPump");
1079      if (returnCode<0) {
1080        if (returnCode==-2)
1081          exitAll=true;
1082        returnCode=0; // returned on size (or event) - could try changing
1083      }
1084      if ((returnCode&2)!=0) {
1085        // could add cut
1086        returnCode &= ~2;
1087      }
1088      if (returnCode&&newSolutionValue<saveValue) {
1089        sprintf(pumpPrint,"Mini branch and bound improved solution from %g to %g (%.2f seconds)",
1090                saveValue,newSolutionValue,model_->getCurrentSeconds());
1091        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1092          << pumpPrint
1093          <<CoinMessageEol;
1094        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
1095        if (fixContinuous) {
1096          // may be able to do even better
1097          const double * lower = model_->solver()->getColLower();
1098          const double * upper = model_->solver()->getColUpper();
1099          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1100            if (newSolver->isInteger(iColumn)) {
1101              double value=floor(newSolution[iColumn]+0.5);
1102              newSolver->setColLower(iColumn,value);
1103              newSolver->setColUpper(iColumn,value);
1104            } else {
1105              newSolver->setColLower(iColumn,lower[iColumn]);
1106              newSolver->setColUpper(iColumn,upper[iColumn]);
1107            }
1108          }
1109          newSolver->initialSolve();
1110          if (newSolver->isProvenOptimal()) {
1111            double value = newSolver->getObjValue()*newSolver->getObjSense();
1112            if (value<newSolutionValue) {
1113              sprintf(pumpPrint,"Freeing continuous variables gives a solution of %g", value);
1114              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1115                << pumpPrint
1116                <<CoinMessageEol;
1117              newSolutionValue=value;
1118              memcpy(betterSolution,newSolver->getColSolution(),numberColumns*sizeof(double));
1119            }
1120          } else {
1121            //newSolver->writeMps("bad3.mps");
1122            exitAll=true;
1123            break;
1124          }
1125        } 
1126        if ((accumulate_&1)!=0) {
1127          model_->incrementUsed(betterSolution); // for local search
1128          numberSolutions++;
1129        }
1130        solutionValue=newSolutionValue;
1131        solutionFound=true;
1132        CoinWarmStartBasis * basis =
1133          dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
1134        if (basis) {
1135          bestBasis = * basis;
1136          delete basis;
1137          CbcEventHandler * handler = model_->getEventHandler();
1138          if (handler) {
1139            double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
1140            double saveObjectiveValue = model_->getMinimizationObjValue();
1141            model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
1142            int action = handler->event(CbcEventHandler::heuristicSolution);
1143            if (saveOldSolution) {
1144              model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
1145              delete [] saveOldSolution;
1146            }
1147            if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
1148              exitAll=true; // exit
1149              break;
1150            }
1151          }
1152        }
1153      } else {
1154        sprintf(pumpPrint,"Mini branch and bound did not improve solution (%.2f seconds)",
1155                model_->getCurrentSeconds());
1156        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1157          << pumpPrint
1158          <<CoinMessageEol;
1159      }
1160      delete newSolver;
1161    }
1162    if (solutionFound) finalReturnCode=1;
1163    cutoff = CoinMin(cutoff,solutionValue);
1164    if (numberTries>=maximumRetries_||!solutionFound||exitAll||cutoff<continuousObjectiveValue+1.0e-7) {
1165      break;
1166    } else {
1167      solutionFound=false;
1168      if (absoluteIncrement_>0.0||relativeIncrement_>0.0) {
1169        double gap = relativeIncrement_*fabs(solutionValue);
1170        cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
1171      } else {
1172        cutoff -= 0.1*(cutoff-continuousObjectiveValue);
1173      }
1174      if (cutoff<continuousObjectiveValue)
1175        break;
1176      sprintf(pumpPrint,"Round again with cutoff of %g",cutoff);
1177      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1178        << pumpPrint
1179        <<CoinMessageEol;
1180      if ((accumulate_&3)<2&&usedColumn)
1181        memset(usedColumn,0,numberColumns);
1182      totalNumberPasses += numberPasses-1;
1183    }
1184  }
1185  delete solver; // probably NULL but do anyway
1186  if (!finalReturnCode&&closestSolution&&closestObjectiveValue <= 10.0&&usedColumn) {
1187    // try a bit of branch and bound
1188    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
1189    const double * colLower = newSolver->getColLower();
1190    const double * colUpper = newSolver->getColUpper();
1191    int i;
1192    double rhs = 0.0;
1193    for (i=0;i<numberIntegersOrig;i++) {
1194      int iColumn=integerVariableOrig[i];
1195      int direction = closestSolution[i];
1196      closestSolution[i]=iColumn;
1197      if (direction==0) {
1198        // keep close to LB
1199        rhs += colLower[iColumn];
1200        lastSolution[i]=1.0;
1201      } else {
1202        // keep close to UB
1203        rhs -= colUpper[iColumn];
1204        lastSolution[i]=-1.0;
1205      }
1206    }
1207    newSolver->addRow(numberIntegersOrig,closestSolution,
1208                      lastSolution,-COIN_DBL_MAX,rhs+10.0);
1209    //double saveValue = newSolutionValue;
1210    //newSolver->writeMps("sub");
1211    int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
1212                                     newSolutionValue,"CbcHeuristicLocalAfterFPump");
1213    if (returnCode<0)
1214      returnCode=0; // returned on size
1215    if ((returnCode&2)!=0) {
1216      // could add cut
1217      returnCode &= ~2;
1218    }
1219    if (returnCode) {
1220      //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
1221      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
1222      //abort();
1223      solutionValue=newSolutionValue;
1224      solutionFound=true;
1225    }
1226    delete newSolver;
1227  }
1228  delete [] roundingSolution;
1229  delete [] usedColumn;
1230  delete [] lastSolution;
1231  delete [] newSolution;
1232  delete [] closestSolution;
1233  delete [] integerVariable;
1234  delete [] firstPerturbedObjective;
1235  delete [] firstPerturbedSolution;
1236  if (solutionValue==incomingObjective) 
1237    sprintf(pumpPrint,"After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
1238            model_->getCurrentSeconds(),CoinCpuTime()-time1);
1239  else
1240    sprintf(pumpPrint,"After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds",
1241            model_->getCurrentSeconds(),solutionValue,CoinCpuTime()-time1);
1242  model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1243    << pumpPrint
1244    <<CoinMessageEol;
1245  if (bestBasis.getNumStructural())
1246    model_->setBestSolutionBasis(bestBasis);
1247  model_->setMinimizationObjValue(saveBestObjective);
1248  if ((accumulate_&1)!=0&&numberSolutions>1&&!model_->getSolutionCount()) {
1249    model_->setSolutionCount(1); // for local search
1250    model_->setNumberHeuristicSolutions(1); 
1251  }
1252  return finalReturnCode;
1253}
1254
1255/**************************END MAIN PROCEDURE ***********************************/
1256
1257// update model
1258void CbcHeuristicFPump::setModel(CbcModel * model)
1259{
1260  model_ = model;
1261}
1262
1263/* Rounds solution - down if < downValue
1264   returns 1 if current is a feasible solution
1265*/
1266int 
1267CbcHeuristicFPump::rounds(OsiSolverInterface * solver,double * solution,
1268                          const double * objective,
1269                          int numberIntegers, const int * integerVariable,
1270                          char * pumpPrint, int & iter,
1271                          bool roundExpensive, double downValue, int *flip)
1272{
1273  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1274  double primalTolerance ;
1275  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
1276
1277  int i;
1278
1279  int numberColumns = model_->getNumCols();
1280  // tmp contains the current obj coefficients
1281  double * tmp = new double [numberColumns];
1282  memcpy(tmp,solver->getObjCoefficients(),numberColumns*sizeof(double));
1283  int flip_up = 0;
1284  int flip_down  = 0;
1285  double  v = randomNumberGenerator_.randomDouble() * 20.0;
1286  int nn = 10 + (int) v;
1287  int nnv = 0;
1288  int * list = new int [nn];
1289  double * val = new double [nn];
1290  for (i = 0; i < nn; i++) val[i] = .001;
1291
1292  // return rounded solution
1293  for (i=0;i<numberIntegers;i++) {
1294    int iColumn = integerVariable[i];
1295    double value=solution[iColumn];
1296    double round = floor(value+primalTolerance);
1297    if (value-round > downValue) round += 1.;
1298    if (round < integerTolerance && tmp[iColumn] < -1. + integerTolerance) flip_down++;
1299    if (round > 1. - integerTolerance && tmp[iColumn] > 1. - integerTolerance) flip_up++;
1300    if (flip_up + flip_down == 0) { 
1301       for (int k = 0; k < nn; k++) {
1302           if (fabs(value-round) > val[k]) {
1303              nnv++;
1304              for (int j = nn-2; j >= k; j--) {
1305                  val[j+1] = val[j];
1306                  list[j+1] = list[j];
1307              } 
1308              val[k] = fabs(value-round);
1309              list[k] = iColumn;
1310              break;
1311           }
1312       }
1313    }
1314    solution[iColumn] = round;
1315  }
1316
1317  if (nnv > nn) nnv = nn;
1318  //if (iter != 0)
1319  //sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
1320  *flip = flip_up + flip_down;
1321  delete [] tmp;
1322
1323  const double * columnLower = solver->getColLower();
1324  const double * columnUpper = solver->getColUpper();
1325  if (*flip == 0 && iter != 0) {
1326    //sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
1327     for (i = 0; i < nnv; i++) {
1328       // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
1329       int index = list[i];
1330       double value = solution[index];
1331       if (value<=1.0) {
1332         solution[index] = 1.0-value;
1333       } else if (value<columnLower[index]+integerTolerance) {
1334         solution[index] = value+1.0;
1335       } else if (value>columnUpper[index]-integerTolerance) {
1336         solution[index] = value-1.0;
1337       } else {
1338         solution[index] = value-1.0;
1339       }
1340     }
1341     *flip = nnv;
1342  } else {
1343    //sprintf(pumpPrint+strlen(pumpPrint)," ");
1344  }
1345  delete [] list; delete [] val;
1346  //iter++;
1347   
1348  const double * rowLower = solver->getRowLower();
1349  const double * rowUpper = solver->getRowUpper();
1350
1351  int numberRows = solver->getNumRows();
1352  // get row activities
1353  double * rowActivity = new double[numberRows];
1354  memset(rowActivity,0,numberRows*sizeof(double));
1355  solver->getMatrixByCol()->times(solution,rowActivity) ;
1356  double largestInfeasibility =0.0;
1357  for (i=0 ; i < numberRows ; i++) {
1358    largestInfeasibility = max(largestInfeasibility,
1359                               rowLower[i]-rowActivity[i]);
1360    largestInfeasibility = max(largestInfeasibility,
1361                               rowActivity[i]-rowUpper[i]);
1362  }
1363  delete [] rowActivity;
1364  return (largestInfeasibility>primalTolerance) ? 0 : 1;
1365}
1366// Set maximum Time (default off) - also sets starttime to current
1367void 
1368CbcHeuristicFPump::setMaximumTime(double value)
1369{
1370  startTime_=CoinCpuTime();
1371  maximumTime_=value;
1372}
1373
1374 
Note: See TracBrowser for help on using the repository browser.