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

Last change on this file since 904 was 904, checked in by ladanyi, 11 years ago

include cstdlib before cmath to get things to compile on AIX with xlC

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 45.4 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,"Not good enough");
580          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
581            << pumpPrint
582            <<CoinMessageEol;
583          //newLineNeeded=false;
584          returnCode=0;
585        }     
586        break;
587      } else {
588        // SOLUTION IS not INTEGER
589        // 1. check for loop
590        bool matched;
591        for (int k = NUMBER_OLD-1; k > 0; k--) {
592          double * b = oldSolution[k];
593          matched = true;
594          for (i = 0; i <numberIntegers; i++) {
595            int iColumn = integerVariable[i];
596            if (newSolution[iColumn]!=b[iColumn]) {
597              matched=false;
598              break;
599            }
600          }
601          if (matched) break;
602        }
603        int numberPerturbed=0;
604        if (matched || numberPasses%100 == 0) {
605          // perturbation
606          //sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
607          //newLineNeeded=true;
608          double factorX[10]={0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0};
609          double factor=1.0;
610          double target=-1.0;
611          double * randomX = new double [numberIntegers];
612          for (i=0;i<numberIntegers;i++) 
613            randomX[i] = max(0.0,randomNumberGenerator_.randomDouble()-0.3);
614          for (int k=0;k<10;k++) {
615#ifdef COIN_DEVELOP_x
616            printf("kpass %d\n",k);
617#endif
618            int numberX[10]={0,0,0,0,0,0,0,0,0,0};
619            for (i=0;i<numberIntegers;i++) {
620              int iColumn = integerVariable[i];
621              double value = randomX[i];
622              double difference = fabs(solution[iColumn]-newSolution[iColumn]);
623              for (int j=0;j<10;j++) {
624                if (difference+value*factorX[j]>0.5) 
625                  numberX[j]++;
626              }
627            }
628            if (target<0.0) {
629              if (numberX[9]<=200)
630                break; // not very many changes
631              target=CoinMax(200.0,CoinMin(0.05*numberX[9],1000.0));
632            }
633            int iX=-1;
634            int iBand=-1;
635            for (i=0;i<10;i++) {
636#ifdef COIN_DEVELOP_x
637              printf("** %d changed at %g\n",numberX[i],factorX[i]);
638#endif
639              if (numberX[i]>=target&&numberX[i]<2.0*target&&iX<0)
640                iX=i;
641              if (iBand<0&&numberX[i]>target) {
642                iBand=i;
643                factor=factorX[i];
644              }
645            }
646            if (iX>=0) {
647              factor=factorX[iX];
648              break;
649            } else {
650              assert (iBand>=0);
651              double hi = factor;
652              double lo = (iBand>0) ? factorX[iBand-1] : 0.0;
653              double diff = (hi-lo)/9.0;
654              for (i=0;i<10;i++) {
655                factorX[i]=lo;
656                lo += diff;
657              }
658            }
659          }
660          for (i=0;i<numberIntegers;i++) {
661            int iColumn = integerVariable[i];
662            double value = randomX[i];
663            double difference = fabs(solution[iColumn]-newSolution[iColumn]);
664            if (difference+value*factor>0.5) {
665              numberPerturbed++;
666              if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
667                newSolution[iColumn] += 1.0;
668              } else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
669                newSolution[iColumn] -= 1.0;
670              } else {
671                // general integer
672                if (difference+value>0.75)
673                  newSolution[iColumn] += 1.0;
674                else
675                  newSolution[iColumn] -= 1.0;
676              }
677            }
678          }
679          delete [] randomX;
680        } else {
681          for (j=NUMBER_OLD-1;j>0;j--) {
682            for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
683          }
684          for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
685        }
686       
687        // 2. update the objective function based on the new rounded solution
688        double offset=0.0;
689        double costValue = (1.0-scaleFactor)*solver->getObjSense();
690        int numberChanged=0;
691        const double * oldObjective = solver->getObjCoefficients();
692        for (i=0;i<numberColumns;i++) {
693          // below so we can keep original code and allow for objective
694          int iColumn = i;
695          // Special code for "artificials"
696          if (direction*saveObjective[iColumn]>=artificialCost_) {
697            //solver->setObjCoeff(iColumn,scaleFactor*saveObjective[iColumn]);
698            solver->setObjCoeff(iColumn,(artificialFactor*saveObjective[iColumn])/artificialCost_);
699          }
700          if(!solver->isBinary(iColumn)&&!doGeneral)
701            continue;
702          // deal with fixed variables (i.e., upper=lower)
703          if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance||!solver->isInteger(iColumn)) {
704            //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
705            //else                                       solver->setObjCoeff(iColumn,costValue);
706            continue;
707          }
708          double newValue=0.0;
709          if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
710            newValue = costValue+scaleFactor*saveObjective[iColumn];
711          } else {
712            if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
713              newValue = -costValue+scaleFactor*saveObjective[iColumn];
714            }
715          }
716          if (newValue!=oldObjective[iColumn])
717            numberChanged++;
718          solver->setObjCoeff(iColumn,newValue);
719          offset += costValue*newSolution[iColumn];
720        }
721        solver->setDblParam(OsiObjOffset,-offset);
722        if (!general&&false) {
723          // Solve in two goes - first keep satisfied ones fixed
724          double * saveLower = new double [numberIntegers];
725          double * saveUpper = new double [numberIntegers];
726          for (i=0;i<numberIntegers;i++) {
727            int iColumn = integerVariable[i];
728            saveLower[i]=COIN_DBL_MAX;
729            saveUpper[i]=-COIN_DBL_MAX;
730            if (solution[iColumn]<lower[iColumn]+primalTolerance) {
731              saveUpper[i]=upper[iColumn];
732              solver->setColUpper(iColumn,lower[iColumn]);
733            } else if (solution[iColumn]>upper[iColumn]-primalTolerance) {
734              saveLower[i]=lower[iColumn];
735              solver->setColLower(iColumn,upper[iColumn]);
736            }
737          }
738          solver->resolve();
739          if (!solver->isProvenOptimal()) {
740            // presumably max time or some such
741            exitAll=true;
742            break;
743          }
744          for (i=0;i<numberIntegers;i++) {
745            int iColumn = integerVariable[i];
746            if (saveLower[i]!=COIN_DBL_MAX)
747              solver->setColLower(iColumn,saveLower[i]);
748            if (saveUpper[i]!=-COIN_DBL_MAX)
749              solver->setColUpper(iColumn,saveUpper[i]);
750            saveUpper[i]=-COIN_DBL_MAX;
751          }
752          memcpy(newSolution,solution,numberColumns*sizeof(double));
753          int flip;
754          returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
755                              pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
756          numberPasses++;
757          if (returnCode) {
758            // solution - but may not be better
759            // Compute using dot product
760            double newSolutionValue = -saveOffset;
761            for (  i=0 ; i<numberColumns ; i++ )
762              newSolutionValue += saveObjective[i]*newSolution[i];
763            newSolutionValue *= direction;
764            sprintf(pumpPrint,"Intermediate solution found of %g",newSolutionValue);
765            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
766              << pumpPrint
767              <<CoinMessageEol;
768            if (newSolutionValue<solutionValue) {
769              memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
770              CoinWarmStartBasis * basis =
771                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
772              solutionFound=true;
773              if (basis) {
774                bestBasis = * basis;
775                delete basis;
776                CbcEventHandler * handler = model_->getEventHandler();
777                if (handler) {
778                  double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
779                  double saveObjectiveValue = model_->getMinimizationObjValue();
780                  model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
781                  int action = handler->event(CbcEventHandler::heuristicSolution);
782                  if (saveOldSolution) {
783                    model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
784                    delete [] saveOldSolution;
785                  }
786                  if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
787                    exitAll=true; // exit
788                    break;
789                  }
790                }
791              }
792              if ((accumulate_&1)!=0) {
793                model_->incrementUsed(betterSolution); // for local search
794                numberSolutions++;
795              }
796              solutionValue=newSolutionValue;
797              solutionFound=true;
798            } else {
799              returnCode=0;
800            }
801          }     
802        }
803        if (!doGeneral) {
804          // faster to do from all slack!!!!
805          if (allSlack) {
806            CoinWarmStartBasis dummy;
807            solver->setWarmStart(&dummy);
808          }
809#ifdef COIN_DEVELOP
810          printf("%d perturbed out of %d columns (%d changed)\n",numberPerturbed,numberColumns,numberChanged);
811#endif
812          bool takeHint;
813          OsiHintStrength strength;
814          solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
815          if (numberPerturbed>numberColumns)
816            solver->setHintParam(OsiDoDualInResolve,true); // dual may be better
817          if (numberTries>1&&numberPasses==1&&false) {
818            // use basis from first time
819            solver->setWarmStart(&saveBasis);
820            // and objective
821            solver->setObjective(firstPerturbedObjective);
822            // and solution
823            solver->setColSolution(firstPerturbedSolution);
824          }
825          solver->resolve();
826          if (!solver->isProvenOptimal()) {
827            // presumably max time or some such
828            exitAll=true;
829            break;
830          }
831          if (numberTries==1&&numberPasses==1&&false) {
832            // save basis
833            CoinWarmStartBasis * basis =
834              dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
835            if (basis) {
836              saveBasis = * basis;
837              delete basis;
838            }
839            firstPerturbedObjective = CoinCopyOfArray(solver->getObjCoefficients(),numberColumns);
840            firstPerturbedSolution = CoinCopyOfArray(solver->getColSolution(),numberColumns);
841          }
842          solver->setHintParam(OsiDoDualInResolve,takeHint);
843          newTrueSolutionValue = -saveOffset;
844          newSumInfeas=0.0;
845          {
846            const double * newSolution = solver->getColSolution();
847            for (  i=0 ; i<numberColumns ; i++ ) {
848              if (solver->isInteger(i)) {
849                double value = newSolution[i];
850                double nearest = floor(value+0.5);
851                newSumInfeas += fabs(value-nearest);
852              }
853              newTrueSolutionValue += saveObjective[i]*newSolution[i];
854            }
855            newTrueSolutionValue *= direction;
856            OsiSolverInterface * saveSolver = model_->swapSolver(solver);
857            CbcRounding heuristic1(*model_);
858            heuristic1.setHeuristicName("rounding in feaspump!");
859            heuristic1.setWhen(1);
860            double testObjectiveValue = CoinMin(solutionValue,roundingObjective);
861            int returnCode = heuristic1.solution(testObjectiveValue,roundingSolution,newTrueSolutionValue) ;
862            if (returnCode==1) {
863              assert(testObjectiveValue < CoinMin(solutionValue,roundingObjective));
864              roundingObjective = testObjectiveValue;
865            }
866            model_->swapSolver(saveSolver);
867          }
868          if (!solver->isProvenOptimal()) {
869            // presumably max time or some such
870            exitAll=true;
871            break;
872          }
873          // in case very dubious solver
874          lower = solver->getColLower();
875          upper = solver->getColUpper();
876          solution = solver->getColSolution();
877        } else {
878          int * addStart = new int[2*general+1];
879          int * addIndex = new int[4*general];
880          double * addElement = new double[4*general];
881          double * addLower = new double[2*general];
882          double * addUpper = new double[2*general];
883          double * obj = new double[general];
884          int nAdd=0;
885          for (i=0;i<numberIntegers;i++) {
886            int iColumn = integerVariable[i];
887            if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
888                newSolution[iColumn]<upper[iColumn]-primalTolerance) {
889              obj[nAdd]=1.0;
890              addLower[nAdd]=0.0;
891              addUpper[nAdd]=COIN_DBL_MAX;
892              nAdd++;
893            }
894          }
895          OsiSolverInterface * solver2 = solver;
896          if (nAdd) {
897            CoinZeroN(addStart,nAdd+1);
898            solver2 = solver->clone();
899            solver2->addCols(nAdd,addStart,NULL,NULL,addLower,addUpper,obj);
900            // feasible solution
901            double * sol = new double[nAdd+numberColumns];
902            memcpy(sol,solution,numberColumns*sizeof(double));
903            // now rows
904            int nAdd=0;
905            int nEl=0;
906            int nAddRow=0;
907            for (i=0;i<numberIntegers;i++) {
908              int iColumn = integerVariable[i];
909              if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
910                  newSolution[iColumn]<upper[iColumn]-primalTolerance) {
911                addLower[nAddRow]=-newSolution[iColumn];;
912                addUpper[nAddRow]=COIN_DBL_MAX;
913                addIndex[nEl] = iColumn;
914                addElement[nEl++]=-1.0;
915                addIndex[nEl] = numberColumns+nAdd;
916                addElement[nEl++]=1.0;
917                nAddRow++;
918                addStart[nAddRow]=nEl;
919                addLower[nAddRow]=newSolution[iColumn];;
920                addUpper[nAddRow]=COIN_DBL_MAX;
921                addIndex[nEl] = iColumn;
922                addElement[nEl++]=1.0;
923                addIndex[nEl] = numberColumns+nAdd;
924                addElement[nEl++]=1.0;
925                nAddRow++;
926                addStart[nAddRow]=nEl;
927                sol[nAdd+numberColumns] = fabs(sol[iColumn]-newSolution[iColumn]);
928                nAdd++;
929              }
930            }
931            solver2->setColSolution(sol);
932            delete [] sol;
933            solver2->addRows(nAddRow,addStart,addIndex,addElement,addLower,addUpper);
934          }
935          delete [] addStart;
936          delete [] addIndex;
937          delete [] addElement;
938          delete [] addLower;
939          delete [] addUpper;
940          delete [] obj;
941          solver2->resolve();
942          if (!solver2->isProvenOptimal()) {
943            // presumably max time or some such
944            exitAll=true;
945            break;
946          }
947          //assert (solver2->isProvenOptimal());
948          if (nAdd) {
949            solver->setColSolution(solver2->getColSolution());
950            delete solver2;
951          }
952        }
953        if (solver->getNumRows()<3000)
954          sprintf(pumpPrint,"Pass %3d: suminf. %10.5f obj. %g iterations %d", numberPasses+totalNumberPasses,
955                  newSumInfeas,newTrueSolutionValue,solver->getIterationCount());
956        else
957          sprintf(pumpPrint,"Pass %3d: (%.2f seconds) suminf. %10.5f obj. %g iterations %d", numberPasses+totalNumberPasses,
958                  model_->getCurrentSeconds(),newSumInfeas,newTrueSolutionValue,solver->getIterationCount());
959        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
960          << pumpPrint
961          <<CoinMessageEol;
962        if (closestSolution&&solver->getObjValue()<closestObjectiveValue) {
963          int i;
964          const double * objective = solver->getObjCoefficients();
965          for (i=0;i<numberIntegersOrig;i++) {
966            int iColumn=integerVariableOrig[i];
967            if (objective[iColumn]>0.0)
968              closestSolution[i]=0;
969            else
970              closestSolution[i]=1;
971          }
972          closestObjectiveValue = solver->getObjValue();
973        }
974        //newLineNeeded=true;
975       
976      }
977      // reduce scale factor
978      scaleFactor *= weightFactor_;
979    } // END WHILE
980    // see if rounding worked!
981    if (roundingObjective<solutionValue) {
982      sprintf(pumpPrint,"Rounding solution of %g is better than previous of %g !\n",
983              roundingObjective,solutionValue);
984      solutionValue=roundingObjective;
985      memcpy(betterSolution,roundingSolution,numberColumns*sizeof(double));
986      solutionFound=true;
987    }
988    if (!solutionFound) { 
989      sprintf(pumpPrint,"No solution found this major pass");
990      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
991        << pumpPrint
992        <<CoinMessageEol;
993    }
994    //}
995    delete solver;
996    solver=NULL;
997    for ( j=0;j<NUMBER_OLD;j++) 
998      delete [] oldSolution[j];
999    delete [] oldSolution;
1000    delete [] saveObjective;
1001    if (usedColumn&&!exitAll) {
1002      OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
1003      const double * colLower = newSolver->getColLower();
1004      const double * colUpper = newSolver->getColUpper();
1005      int i;
1006      int nFix=0;
1007      int nFixI=0;
1008      int nFixC=0;
1009      int nFixC2=0;
1010      for (i=0;i<numberIntegersOrig;i++) {
1011        int iColumn=integerVariableOrig[i];
1012        //const OsiObject * object = model_->object(i);
1013        //double originalLower;
1014        //double originalUpper;
1015        //getIntegerInformation( object,originalLower, originalUpper);
1016        //assert(colLower[iColumn]==originalLower);
1017        //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
1018        newSolver->setColLower(iColumn,colLower[iColumn]);
1019        //assert(colUpper[iColumn]==originalUpper);
1020        //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
1021        newSolver->setColUpper(iColumn,colUpper[iColumn]);
1022        if (!usedColumn[iColumn]) {
1023          double value=lastSolution[iColumn];
1024          double nearest = floor(value+0.5);
1025          if (fabs(value-nearest)<1.0e-7) {
1026            if (nearest==colLower[iColumn]) {
1027              newSolver->setColUpper(iColumn,colLower[iColumn]);
1028              nFix++;
1029            } else if (nearest==colUpper[iColumn]) {
1030              newSolver->setColLower(iColumn,colUpper[iColumn]);
1031              nFix++;
1032            } else if (fixInternal) {
1033              newSolver->setColLower(iColumn,nearest);
1034              newSolver->setColUpper(iColumn,nearest);
1035              nFix++;
1036              nFixI++;
1037            }
1038          }
1039        }
1040      }
1041      if (fixContinuous) {
1042        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1043          if (!newSolver->isInteger(iColumn)&&!usedColumn[iColumn]) {
1044            double value=lastSolution[iColumn];
1045            if (value<colLower[iColumn]+1.0e-8) {
1046              newSolver->setColUpper(iColumn,colLower[iColumn]);
1047              nFixC++;
1048            } else if (value>colUpper[iColumn]-1.0e-8) {
1049              newSolver->setColLower(iColumn,colUpper[iColumn]);
1050              nFixC++;
1051            } else if (fixContinuous==2) {
1052              newSolver->setColLower(iColumn,value);
1053              newSolver->setColUpper(iColumn,value);
1054              nFixC++;
1055              nFixC2++;
1056            }
1057          }
1058        }
1059      }
1060      newSolver->initialSolve();
1061      if (!newSolver->isProvenOptimal()) {
1062        //newSolver->writeMps("bad.mps");
1063        //assert (newSolver->isProvenOptimal());
1064        exitAll=true;
1065        break;
1066      }
1067      sprintf(pumpPrint,"Before mini branch and bound, %d integers at bound fixed and %d continuous",
1068             nFix,nFixC);
1069      if (nFixC2+nFixI!=0)
1070        sprintf(pumpPrint+strlen(pumpPrint)," of which %d were internal integer and %d internal continuous",
1071                nFixI,nFixC2);
1072      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1073        << pumpPrint
1074        <<CoinMessageEol;
1075      double saveValue = newSolutionValue;
1076      returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
1077                                       cutoff,"CbcHeuristicLocalAfterFPump");
1078      if (returnCode<0) {
1079        if (returnCode==-2)
1080          exitAll=true;
1081        returnCode=0; // returned on size (or event) - could try changing
1082      }
1083      if ((returnCode&2)!=0) {
1084        // could add cut
1085        returnCode &= ~2;
1086      }
1087      if (returnCode&&newSolutionValue<saveValue) {
1088        sprintf(pumpPrint,"Mini branch and bound improved solution from %g to %g (%.2f seconds)",
1089                saveValue,newSolutionValue,model_->getCurrentSeconds());
1090        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1091          << pumpPrint
1092          <<CoinMessageEol;
1093        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
1094        if (fixContinuous) {
1095          // may be able to do even better
1096          const double * lower = model_->solver()->getColLower();
1097          const double * upper = model_->solver()->getColUpper();
1098          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1099            if (newSolver->isInteger(iColumn)) {
1100              double value=floor(newSolution[iColumn]+0.5);
1101              newSolver->setColLower(iColumn,value);
1102              newSolver->setColUpper(iColumn,value);
1103            } else {
1104              newSolver->setColLower(iColumn,lower[iColumn]);
1105              newSolver->setColUpper(iColumn,upper[iColumn]);
1106            }
1107          }
1108          newSolver->initialSolve();
1109          if (newSolver->isProvenOptimal()) {
1110            double value = newSolver->getObjValue()*newSolver->getObjSense();
1111            if (value<newSolutionValue) {
1112              sprintf(pumpPrint,"Freeing continuous variables gives a solution of %g", value);
1113              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1114                << pumpPrint
1115                <<CoinMessageEol;
1116              newSolutionValue=value;
1117              memcpy(betterSolution,newSolver->getColSolution(),numberColumns*sizeof(double));
1118            }
1119          } else {
1120            //newSolver->writeMps("bad3.mps");
1121            exitAll=true;
1122            break;
1123          }
1124        } 
1125        if ((accumulate_&1)!=0) {
1126          model_->incrementUsed(betterSolution); // for local search
1127          numberSolutions++;
1128        }
1129        solutionValue=newSolutionValue;
1130        solutionFound=true;
1131        CoinWarmStartBasis * basis =
1132          dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
1133        if (basis) {
1134          bestBasis = * basis;
1135          delete basis;
1136          CbcEventHandler * handler = model_->getEventHandler();
1137          if (handler) {
1138            double * saveOldSolution = CoinCopyOfArray(model_->bestSolution(),numberColumns);
1139            double saveObjectiveValue = model_->getMinimizationObjValue();
1140            model_->setBestSolution(betterSolution,numberColumns,newSolutionValue);
1141            int action = handler->event(CbcEventHandler::heuristicSolution);
1142            if (saveOldSolution) {
1143              model_->setBestSolution(saveOldSolution,numberColumns,saveObjectiveValue);
1144              delete [] saveOldSolution;
1145            }
1146            if (!action||model_->getCurrentSeconds()>model_->getMaximumSeconds()) {
1147              exitAll=true; // exit
1148              break;
1149            }
1150          }
1151        }
1152      } else {
1153        sprintf(pumpPrint,"Mini branch and bound did not improve solution (%.2f seconds)",
1154                model_->getCurrentSeconds());
1155        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1156          << pumpPrint
1157          <<CoinMessageEol;
1158      }
1159      delete newSolver;
1160    }
1161    if (solutionFound) finalReturnCode=1;
1162    cutoff = CoinMin(cutoff,solutionValue);
1163    if (numberTries>=maximumRetries_||!solutionFound||exitAll||cutoff<continuousObjectiveValue+1.0e-7) {
1164      break;
1165    } else {
1166      solutionFound=false;
1167      if (absoluteIncrement_>0.0||relativeIncrement_>0.0) {
1168        double gap = relativeIncrement_*fabs(solutionValue);
1169        cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
1170      } else {
1171        cutoff -= 0.1*(cutoff-continuousObjectiveValue);
1172      }
1173      if (cutoff<continuousObjectiveValue)
1174        break;
1175      sprintf(pumpPrint,"Round again with cutoff of %g",cutoff);
1176      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1177        << pumpPrint
1178        <<CoinMessageEol;
1179      if ((accumulate_&3)<2&&usedColumn)
1180        memset(usedColumn,0,numberColumns);
1181      totalNumberPasses += numberPasses-1;
1182    }
1183  }
1184  delete solver; // probably NULL but do anyway
1185  if (!finalReturnCode&&closestSolution&&closestObjectiveValue <= 10.0&&usedColumn) {
1186    // try a bit of branch and bound
1187    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
1188    const double * colLower = newSolver->getColLower();
1189    const double * colUpper = newSolver->getColUpper();
1190    int i;
1191    double rhs = 0.0;
1192    for (i=0;i<numberIntegersOrig;i++) {
1193      int iColumn=integerVariableOrig[i];
1194      int direction = closestSolution[i];
1195      closestSolution[i]=iColumn;
1196      if (direction==0) {
1197        // keep close to LB
1198        rhs += colLower[iColumn];
1199        lastSolution[i]=1.0;
1200      } else {
1201        // keep close to UB
1202        rhs -= colUpper[iColumn];
1203        lastSolution[i]=-1.0;
1204      }
1205    }
1206    newSolver->addRow(numberIntegersOrig,closestSolution,
1207                      lastSolution,-COIN_DBL_MAX,rhs+10.0);
1208    //double saveValue = newSolutionValue;
1209    //newSolver->writeMps("sub");
1210    int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
1211                                     newSolutionValue,"CbcHeuristicLocalAfterFPump");
1212    if (returnCode<0)
1213      returnCode=0; // returned on size
1214    if ((returnCode&2)!=0) {
1215      // could add cut
1216      returnCode &= ~2;
1217    }
1218    if (returnCode) {
1219      //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
1220      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
1221      //abort();
1222      solutionValue=newSolutionValue;
1223      solutionFound=true;
1224    }
1225    delete newSolver;
1226  }
1227  delete [] roundingSolution;
1228  delete [] usedColumn;
1229  delete [] lastSolution;
1230  delete [] newSolution;
1231  delete [] closestSolution;
1232  delete [] integerVariable;
1233  delete [] firstPerturbedObjective;
1234  delete [] firstPerturbedSolution;
1235  if (solutionValue==incomingObjective) 
1236    sprintf(pumpPrint,"After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
1237            model_->getCurrentSeconds(),CoinCpuTime()-time1);
1238  else
1239    sprintf(pumpPrint,"After %.2f seconds - Feasibility pump exiting with objective of %g - took %.2f seconds",
1240            model_->getCurrentSeconds(),solutionValue,CoinCpuTime()-time1);
1241  model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
1242    << pumpPrint
1243    <<CoinMessageEol;
1244  if (bestBasis.getNumStructural())
1245    model_->setBestSolutionBasis(bestBasis);
1246  model_->setMinimizationObjValue(saveBestObjective);
1247  if ((accumulate_&1)!=0&&numberSolutions>1&&!model_->getSolutionCount()) {
1248    model_->setSolutionCount(1); // for local search
1249    model_->setNumberHeuristicSolutions(1); 
1250  }
1251  return finalReturnCode;
1252}
1253
1254/**************************END MAIN PROCEDURE ***********************************/
1255
1256// update model
1257void CbcHeuristicFPump::setModel(CbcModel * model)
1258{
1259  model_ = model;
1260}
1261
1262/* Rounds solution - down if < downValue
1263   returns 1 if current is a feasible solution
1264*/
1265int 
1266CbcHeuristicFPump::rounds(OsiSolverInterface * solver,double * solution,
1267                          const double * objective,
1268                          int numberIntegers, const int * integerVariable,
1269                          char * pumpPrint, int & iter,
1270                          bool roundExpensive, double downValue, int *flip)
1271{
1272  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1273  double primalTolerance ;
1274  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
1275
1276  int i;
1277
1278  int numberColumns = model_->getNumCols();
1279  // tmp contains the current obj coefficients
1280  double * tmp = new double [numberColumns];
1281  memcpy(tmp,solver->getObjCoefficients(),numberColumns*sizeof(double));
1282  int flip_up = 0;
1283  int flip_down  = 0;
1284  double  v = randomNumberGenerator_.randomDouble() * 20.0;
1285  int nn = 10 + (int) v;
1286  int nnv = 0;
1287  int * list = new int [nn];
1288  double * val = new double [nn];
1289  for (i = 0; i < nn; i++) val[i] = .001;
1290
1291  // return rounded solution
1292  for (i=0;i<numberIntegers;i++) {
1293    int iColumn = integerVariable[i];
1294    double value=solution[iColumn];
1295    double round = floor(value+primalTolerance);
1296    if (value-round > downValue) round += 1.;
1297    if (round < integerTolerance && tmp[iColumn] < -1. + integerTolerance) flip_down++;
1298    if (round > 1. - integerTolerance && tmp[iColumn] > 1. - integerTolerance) flip_up++;
1299    if (flip_up + flip_down == 0) { 
1300       for (int k = 0; k < nn; k++) {
1301           if (fabs(value-round) > val[k]) {
1302              nnv++;
1303              for (int j = nn-2; j >= k; j--) {
1304                  val[j+1] = val[j];
1305                  list[j+1] = list[j];
1306              } 
1307              val[k] = fabs(value-round);
1308              list[k] = iColumn;
1309              break;
1310           }
1311       }
1312    }
1313    solution[iColumn] = round;
1314  }
1315
1316  if (nnv > nn) nnv = nn;
1317  //if (iter != 0)
1318  //sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
1319  *flip = flip_up + flip_down;
1320  delete [] tmp;
1321
1322  const double * columnLower = solver->getColLower();
1323  const double * columnUpper = solver->getColUpper();
1324  if (*flip == 0 && iter != 0) {
1325    //sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
1326     for (i = 0; i < nnv; i++) {
1327       // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
1328       int index = list[i];
1329       double value = solution[index];
1330       if (value<=1.0) {
1331         solution[index] = 1.0-value;
1332       } else if (value<columnLower[index]+integerTolerance) {
1333         solution[index] = value+1.0;
1334       } else if (value>columnUpper[index]-integerTolerance) {
1335         solution[index] = value-1.0;
1336       } else {
1337         solution[index] = value-1.0;
1338       }
1339     }
1340     *flip = nnv;
1341  } else {
1342    //sprintf(pumpPrint+strlen(pumpPrint)," ");
1343  }
1344  delete [] list; delete [] val;
1345  //iter++;
1346   
1347  const double * rowLower = solver->getRowLower();
1348  const double * rowUpper = solver->getRowUpper();
1349
1350  int numberRows = solver->getNumRows();
1351  // get row activities
1352  double * rowActivity = new double[numberRows];
1353  memset(rowActivity,0,numberRows*sizeof(double));
1354  solver->getMatrixByCol()->times(solution,rowActivity) ;
1355  double largestInfeasibility =0.0;
1356  for (i=0 ; i < numberRows ; i++) {
1357    largestInfeasibility = max(largestInfeasibility,
1358                               rowLower[i]-rowActivity[i]);
1359    largestInfeasibility = max(largestInfeasibility,
1360                               rowActivity[i]-rowUpper[i]);
1361  }
1362  delete [] rowActivity;
1363  return (largestInfeasibility>primalTolerance) ? 0 : 1;
1364}
1365// Set maximum Time (default off) - also sets starttime to current
1366void 
1367CbcHeuristicFPump::setMaximumTime(double value)
1368{
1369  startTime_=CoinCpuTime();
1370  maximumTime_=value;
1371}
1372
1373 
Note: See TracBrowser for help on using the repository browser.