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

Last change on this file since 758 was 758, checked in by forrest, 13 years ago

fpump - make dj fixing optional

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