source: branches/heur/Cbc/src/CbcHeuristicDiveVectorLength.cpp @ 885

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

add diving heuristics

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