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

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

add diving heuristics

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