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