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