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

Last change on this file since 759 was 759, checked in by forrest, 14 years ago

start of work on new vub heuristic

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