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

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

try and speed up feasibility pump

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