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

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

adding diving heuristics from Joao

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