source: trunk/Cbc/src/CbcHeuristicDiveGuided.cpp @ 871

Last change on this file since 871 was 871, checked in by forrest, 11 years ago

add diving heuristics

File size: 14.0 KB
Line 
1// Copyright (C) 2008, 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
8#include "CbcHeuristicDiveGuided.hpp"
9#include "CbcStrategy.hpp"
10#include  "CoinTime.hpp"
11
12// Default Constructor
13CbcHeuristicDiveGuided::CbcHeuristicDiveGuided() 
14  :CbcHeuristic()
15{
16  // matrix and row copy will automatically be empty
17  downLocks_ =NULL;
18  upLocks_ = NULL;
19  percentageToFix_ = 0.2;
20  maxIterations_ = 100;
21  maxTime_ = 60;
22}
23
24// Constructor from model
25CbcHeuristicDiveGuided::CbcHeuristicDiveGuided(CbcModel & model)
26  :CbcHeuristic(model)
27{
28  downLocks_ =NULL;
29  upLocks_ = NULL;
30  // Get a copy of original matrix
31  assert(model.solver());
32  // model may have empty matrix - wait until setModel
33  const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
34  if (matrix) {
35    matrix_ = *matrix;
36  }
37  percentageToFix_ = 0.2;
38  maxIterations_ = 100;
39  maxTime_ = 60;
40}
41
42// Destructor
43CbcHeuristicDiveGuided::~CbcHeuristicDiveGuided ()
44{
45  delete [] downLocks_;
46  delete [] upLocks_;
47}
48
49// Clone
50CbcHeuristicDiveGuided *
51CbcHeuristicDiveGuided::clone() const
52{
53  return new CbcHeuristicDiveGuided(*this);
54}
55
56// Create C++ lines to get to current state
57void 
58CbcHeuristicDiveGuided::generateCpp( FILE * fp) 
59{
60  CbcHeuristicDiveGuided other;
61  fprintf(fp,"0#include \"CbcHeuristicDiveGuided.hpp\"\n");
62  fprintf(fp,"3  CbcHeuristicDiveGuided heuristicDiveGuided(*cbcModel);\n");
63  CbcHeuristic::generateCpp(fp,"heuristicDiveGuided");
64  if (percentageToFix_!=other.percentageToFix_)
65    fprintf(fp,"3  heuristicDiveGuided.setPercentageToFix(%.2f);\n",percentageToFix_);
66  else
67    fprintf(fp,"4  heuristicDiveGuided.setPercentageToFix(%.2f);\n",percentageToFix_);
68  if (maxIterations_!=other.maxIterations_)
69    fprintf(fp,"3  heuristicDiveGuided.setMaxIterations(%d);\n",maxIterations_);
70  else
71    fprintf(fp,"4  heuristicDiveGuided.setMaxIterations(%d);\n",maxIterations_);
72  if (maxTime_!=other.maxTime_)
73    fprintf(fp,"3  heuristicDiveGuided.setMaxTime(%.2f);\n",maxTime_);
74  else
75    fprintf(fp,"4  heuristicDiveGuided.setMaxTime(%.2f);\n",maxTime_);
76  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicDiveGuided);\n");
77}
78
79// Copy constructor
80CbcHeuristicDiveGuided::CbcHeuristicDiveGuided(const CbcHeuristicDiveGuided & rhs)
81:
82  CbcHeuristic(rhs),
83  matrix_(rhs.matrix_),
84  percentageToFix_(rhs.percentageToFix_),
85  maxIterations_(rhs.maxIterations_),
86  maxTime_(rhs.maxTime_)
87{
88  int numberIntegers = model_->numberIntegers();
89  if (rhs.downLocks_) {
90    int numberIntegers = model_->numberIntegers();
91    downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
92    upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
93  } else {
94    downLocks_ = NULL;
95    upLocks_ = NULL;
96  }
97}
98
99// Assignment operator
100CbcHeuristicDiveGuided & 
101CbcHeuristicDiveGuided::operator=( const CbcHeuristicDiveGuided& rhs)
102{
103  if (this!=&rhs) {
104    CbcHeuristic::operator=(rhs);
105    matrix_ = rhs.matrix_;
106    percentageToFix_ = rhs.percentageToFix_;
107    maxIterations_ = rhs.maxIterations_;
108    maxTime_ = rhs.maxTime_;
109    delete [] downLocks_;
110    delete [] upLocks_;
111    if (rhs.downLocks_) {
112      int numberIntegers = model_->numberIntegers();
113      downLocks_ = CoinCopyOfArray(rhs.downLocks_,numberIntegers);
114      upLocks_ = CoinCopyOfArray(rhs.upLocks_,numberIntegers);
115    } else {
116      downLocks_ = NULL;
117      upLocks_ = NULL;
118    }
119  }
120  return *this;
121}
122
123// Resets stuff if model changes
124void 
125CbcHeuristicDiveGuided::resetModel(CbcModel * model)
126{
127  model_=model;
128  // Get a copy of original matrix
129  assert(model_->solver());
130  // model may have empty matrix - wait until setModel
131  const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
132  if (matrix) {
133    matrix_ = *matrix;
134    validate();
135  }
136}
137
138// See if dive guided will give better solution
139// Sets value of solution
140// Returns 1 if solution, 0 if not
141int
142CbcHeuristicDiveGuided::solution(double & solutionValue,
143                                 double * betterSolution)
144{
145
146  // See if to do
147  if (!when()||(when()%10==1&&model_->phase()!=1)||
148      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
149    return 0; // switched off
150
151  double* bestIntegerSolution = model_->bestSolution();
152  if(bestIntegerSolution == NULL)
153    return 0; // no integer solution available. Switch off heuristic
154
155  double time1 = CoinCpuTime();
156
157  OsiSolverInterface * solver = model_->solver()->clone();
158  const double * lower = solver->getColLower();
159  const double * upper = solver->getColUpper();
160  const double * rowLower = solver->getRowLower();
161  const double * rowUpper = solver->getRowUpper();
162  const double * solution = solver->getColSolution();
163  const double * objective = solver->getObjCoefficients();
164  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
165  double primalTolerance;
166  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
167
168  int numberRows = matrix_.getNumRows();
169  assert (numberRows<=solver->getNumRows());
170  int numberIntegers = model_->numberIntegers();
171  const int * integerVariable = model_->integerVariable();
172  double direction = solver->getObjSense(); // 1 for min, -1 for max
173  double newSolutionValue = direction*solver->getObjValue();
174  int returnCode = 0;
175  // Column copy
176  const double * element = matrix_.getElements();
177  const int * row = matrix_.getIndices();
178  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
179  const int * columnLength = matrix_.getVectorLengths();
180
181  // Get solution array for heuristic solution
182  int numberColumns = solver->getNumCols();
183  double * newSolution = new double [numberColumns];
184  memcpy(newSolution,solution,numberColumns*sizeof(double));
185
186  // vectors to store the latest variables fixed at their bounds
187  int* columnFixed = new int [numberIntegers];
188  double* originalBound = new double [numberIntegers];
189  bool * fixedAtLowerBound = new bool [numberIntegers];
190
191  const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers);
192
193  // count how many fractional variables
194  int numberFractionalVariables = 0;
195  for (int i=0; i<numberIntegers; i++) {
196    int iColumn = integerVariable[i];
197    double value=newSolution[iColumn];
198    if (fabs(floor(value+0.5)-value)>integerTolerance) {
199      numberFractionalVariables++;
200    }
201  }
202
203  int iteration = 0;
204  while(numberFractionalVariables) {
205    iteration++;
206
207    // select a fractional variable to bound
208    double bestFraction = DBL_MAX;
209    int bestColumn = -1;
210    int bestRound = -1; // -1 rounds down, +1 rounds up
211    int numberAtBoundFixed = 0;
212    bool canRoundSolution = true;
213    for (int i=0; i<numberIntegers; i++) {
214      int iColumn = integerVariable[i];
215      double value=newSolution[iColumn];
216      double fraction=value-floor(value);
217      int round = 0;
218      if (fabs(floor(value+0.5)-value)>integerTolerance) {
219        if(downLocks_[i]>0&&upLocks_[i]>0) {
220          // the variable cannot be rounded
221          canRoundSolution = false;
222          if(value >= bestIntegerSolution[iColumn])
223            round = -1;
224          else {
225            round = 1;
226            fraction = 1.0 - fraction;
227          }
228
229          // if variable is not binary, penalize it
230          if(!solver->isBinary(iColumn))
231            fraction *= 1000.0;
232
233          if(fraction < bestFraction) {
234            bestColumn = iColumn;
235            bestFraction = fraction;
236            bestRound = round;
237          }
238        }
239      }
240      else if(numberAtBoundFixed < maxNumberAtBoundToFix) {
241        // fix the variable at one of its bounds
242        if (fabs(lower[iColumn]-value)<=integerTolerance &&
243            lower[iColumn] != upper[iColumn]) {
244          columnFixed[numberAtBoundFixed] = iColumn;
245          originalBound[numberAtBoundFixed] = upper[iColumn];
246          fixedAtLowerBound[numberAtBoundFixed] = true;
247          solver->setColUpper(iColumn, lower[iColumn]);
248          numberAtBoundFixed++;
249        }
250        else if(fabs(upper[iColumn]-value)<=integerTolerance &&
251            lower[iColumn] != upper[iColumn]) {
252          columnFixed[numberAtBoundFixed] = iColumn;
253          originalBound[numberAtBoundFixed] = lower[iColumn];
254          fixedAtLowerBound[numberAtBoundFixed] = false;
255          solver->setColLower(iColumn, upper[iColumn]);
256          numberAtBoundFixed++;
257        }
258      }
259    }
260
261    if(canRoundSolution) {
262      // Round all the fractional variables
263      for (int i=0; i<numberIntegers; i++) {
264        int iColumn = integerVariable[i];
265        double value=newSolution[iColumn];
266        if (fabs(floor(value+0.5)-value)>integerTolerance) {
267          if(downLocks_[i]==0 || upLocks_[i]==0) {
268            if(downLocks_[i]==0 && upLocks_[i]==0) {
269              if(direction * objective[iColumn] >= 0.0)
270                newSolution[iColumn] = floor(value);
271              else
272                newSolution[iColumn] = ceil(value);
273            }
274            else if(downLocks_[i]==0)
275              newSolution[iColumn] = floor(value);
276            else
277              newSolution[iColumn] = ceil(value);
278          }
279          else
280            break;
281        }
282      }
283      break;
284    }
285
286
287    double originalBoundBestColumn;
288    if(bestColumn >= 0) {
289      if(bestRound < 0) {
290        originalBoundBestColumn = upper[bestColumn];
291        solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
292      }
293      else {
294        originalBoundBestColumn = lower[bestColumn];
295        solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
296      }
297    } else {
298      break;
299    }
300    int originalBestRound = bestRound;
301    while (1) {
302
303      solver->resolve();
304
305      if(!solver->isProvenOptimal()) {
306        if(numberAtBoundFixed > 0) {
307          // Remove the bound fix for variables that were at bounds
308          for(int i=0; i<numberAtBoundFixed; i++) {
309            int iColFixed = columnFixed[i];
310            if(fixedAtLowerBound[i])
311              solver->setColUpper(iColFixed, originalBound[i]);
312            else
313              solver->setColLower(iColFixed, originalBound[i]);
314          }
315          numberAtBoundFixed = 0;
316        }
317        else if(bestRound == originalBestRound) {
318          bestRound *= (-1);
319          if(bestRound < 0) {
320            solver->setColLower(bestColumn, originalBoundBestColumn);
321            solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
322          }
323          else {
324            solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
325            solver->setColUpper(bestColumn, originalBoundBestColumn);
326          }
327        }
328        else
329          break;
330      }
331      else
332        break;
333    }
334
335    if(!solver->isProvenOptimal())
336      break;
337
338    if(iteration > maxIterations_) {
339      break;
340    }
341
342    if(CoinCpuTime()-time1 > maxTime_) {
343      break;
344    }
345
346    memcpy(newSolution,solution,numberColumns*sizeof(double));
347    numberFractionalVariables = 0;
348    for (int i=0; i<numberIntegers; i++) {
349      int iColumn = integerVariable[i];
350      double value=newSolution[iColumn];
351      if (fabs(floor(value+0.5)-value)>integerTolerance) {
352        numberFractionalVariables++;
353      }
354    }
355
356  }
357
358
359  double * rowActivity = new double[numberRows];
360  memset(rowActivity,0,numberRows*sizeof(double));
361
362  // re-compute new solution value
363  double objOffset=0.0;
364  solver->getDblParam(OsiObjOffset,objOffset);
365  newSolutionValue = -objOffset;
366  for (int i=0 ; i<numberColumns ; i++ )
367    newSolutionValue += objective[i]*newSolution[i];
368  newSolutionValue *= direction;
369    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
370  if (newSolutionValue<solutionValue) {
371    // paranoid check
372    memset(rowActivity,0,numberRows*sizeof(double));
373    for (int i=0;i<numberColumns;i++) {
374      int j;
375      double value = newSolution[i];
376      if (value) {
377        for (j=columnStart[i];
378             j<columnStart[i]+columnLength[i];j++) {
379          int iRow=row[j];
380          rowActivity[iRow] += value*element[j];
381        }
382      }
383    }
384    // check was approximately feasible
385    bool feasible=true;
386    for (int i=0;i<numberRows;i++) {
387      if(rowActivity[i]<rowLower[i]) {
388        if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
389          feasible = false;
390      } else if(rowActivity[i]>rowUpper[i]) {
391        if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
392          feasible = false;
393      }
394    }
395    for (int i=0; i<numberIntegers; i++) {
396      int iColumn = integerVariable[i];
397      double value=newSolution[iColumn];
398      if (fabs(floor(value+0.5)-value)>integerTolerance) {
399        feasible = false;
400        break;
401      }
402    }
403    if (feasible) {
404      // new solution
405      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
406      solutionValue = newSolutionValue;
407      //printf("** Solution of %g found by CbcHeuristicDiveGuided\n",newSolutionValue);
408      returnCode=1;
409    } else {
410      // Can easily happen
411      //printf("Debug CbcHeuristicDiveGuided giving bad solution\n");
412    }
413  }
414
415  delete [] newSolution;
416  delete [] columnFixed;
417  delete [] originalBound;
418  delete [] fixedAtLowerBound;
419  delete [] rowActivity;
420  delete solver;
421  return returnCode;
422}
423
424// update model
425void CbcHeuristicDiveGuided::setModel(CbcModel * model)
426{
427  model_ = model;
428  // Get a copy of original matrix (and by row for rounding);
429  assert(model_->solver());
430  const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
431  if (matrix) {
432    matrix_ = *matrix;
433    // make sure model okay for heuristic
434    validate();
435  }
436}
437
438// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
439void 
440CbcHeuristicDiveGuided::validate() 
441{
442  if (model_&&when()<10) {
443    if (model_->numberIntegers()!=
444        model_->numberObjects())
445      setWhen(0);
446  }
447
448  int numberIntegers = model_->numberIntegers();
449  const int * integerVariable = model_->integerVariable();
450  delete [] downLocks_;
451  delete [] upLocks_;
452  downLocks_ = new unsigned short [numberIntegers];
453  upLocks_ = new unsigned short [numberIntegers];
454  // Column copy
455  const double * element = matrix_.getElements();
456  const int * row = matrix_.getIndices();
457  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
458  const int * columnLength = matrix_.getVectorLengths();
459  const double * rowLower = model_->solver()->getRowLower();
460  const double * rowUpper = model_->solver()->getRowUpper();
461  for (int i=0;i<numberIntegers;i++) {
462    int iColumn = integerVariable[i];
463    int down=0;
464    int up=0;
465    if (columnLength[iColumn]>65535) {
466      setWhen(0);
467      break; // unlikely to work
468    }
469    for (CoinBigIndex j=columnStart[iColumn];
470         j<columnStart[iColumn]+columnLength[iColumn];j++) {
471      int iRow=row[j];
472      if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
473        up++;
474        down++;
475      } else if (element[j]>0.0) {
476        if (rowUpper[iRow]<1.0e20)
477          up++;
478        else
479          down++;
480      } else {
481        if (rowLower[iRow]>-1.0e20)
482          up++;
483        else
484          down++;
485      }
486    }
487    downLocks_[i] = (unsigned short) down;
488    upLocks_[i] = (unsigned short) up;
489  }
490}
Note: See TracBrowser for help on using the repository browser.