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