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

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

max times and fix bug in fpump

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