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

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

update feasibility pump

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