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

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

for deterministic parallel

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