source: branches/heur/Cbc/src/CbcHeuristicDiveCoefficient.cpp @ 884

Last change on this file since 884 was 884, checked in by jpgoncal, 11 years ago

Draft code for heuristic to control heuristics.

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