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

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

Cleaned up a little bit.

File size: 15.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 "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  downLocks_ =NULL;
29  upLocks_ = NULL;
30  // Get a copy of original matrix
31  assert(model.solver());
32  // model may have empty matrix - wait until setModel
33  const CoinPackedMatrix * matrix = model.solver()->getMatrixByCol();
34  if (matrix) {
35    matrix_ = *matrix;
36    validate();
37  }
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_)
66    fprintf(fp,"3  heuristicDiveCoefficient.setPercentageToFix(%.2f);\n",percentageToFix_);
67  else
68    fprintf(fp,"4  heuristicDiveCoefficient.setPercentageToFix(%.2f);\n",percentageToFix_);
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_)
74    fprintf(fp,"3  heuristicDiveCoefficient.setMaxTime(%.2f);\n",maxTime_);
75  else
76    fprintf(fp,"4  heuristicDiveCoefficient.setMaxTime(%.2f);\n",maxTime_);
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();
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  }
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_;
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    }
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());
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  }
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{
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    }
156
157    // Get where we are and create the appropriate CbcHeuristicNode object
158    nodeDesc = new CbcHeuristicNode(*model_);
159    if (!runNodes_->farFrom(nodeDesc)) {
160      delete nodeDesc;
161      return NULL;
162    }
163  }
164
165  if(nodeDesc != NULL)
166    runNodes_->append(nodeDesc);
167 
168  lastRun_ = nodeCount;
169
170#if 0
171  // See if to do
172  if (!when()||(when()%10==1&&model_->phase()!=1)||
173      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
174    return 0; // switched off
175#endif
176
177  double time1 = CoinCpuTime();
178
179  OsiSolverInterface * solver = model_->solver()->clone();
180  const double * lower = solver->getColLower();
181  const double * upper = solver->getColUpper();
182  const double * rowLower = solver->getRowLower();
183  const double * rowUpper = solver->getRowUpper();
184  const double * solution = solver->getColSolution();
185  const double * objective = solver->getObjCoefficients();
186  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
187  double primalTolerance;
188  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
189
190  int numberRows = matrix_.getNumRows();
191  assert (numberRows<=solver->getNumRows());
192  int numberIntegers = model_->numberIntegers();
193  const int * integerVariable = model_->integerVariable();
194  double direction = solver->getObjSense(); // 1 for min, -1 for max
195  double newSolutionValue = direction*solver->getObjValue();
196  int returnCode = 0;
197  // Column copy
198  const double * element = matrix_.getElements();
199  const int * row = matrix_.getIndices();
200  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
201  const int * columnLength = matrix_.getVectorLengths();
202
203  // Get solution array for heuristic solution
204  int numberColumns = solver->getNumCols();
205  double * newSolution = new double [numberColumns];
206  memcpy(newSolution,solution,numberColumns*sizeof(double));
207
208  // vectors to store the latest variables fixed at their bounds
209  int* columnFixed = new int [numberIntegers];
210  double* originalBound = new double [numberIntegers];
211  bool * fixedAtLowerBound = new bool [numberIntegers];
212
213  const int maxNumberAtBoundToFix = (int) floor(percentageToFix_ * numberIntegers);
214
215  // count how many fractional variables
216  int numberFractionalVariables = 0;
217  for (int i=0; i<numberIntegers; i++) {
218    int iColumn = integerVariable[i];
219    double value=newSolution[iColumn];
220    if (fabs(floor(value+0.5)-value)>integerTolerance) {
221      numberFractionalVariables++;
222    }
223  }
224
225  int iteration = 0;
226  while(numberFractionalVariables) {
227    iteration++;
228
229    // select a fractional variable to bound
230    double bestFraction = DBL_MAX;
231    int bestColumn = -1;
232    int bestRound = -1; // -1 rounds down, +1 rounds up
233    int bestLocks = COIN_INT_MAX;
234    int numberAtBoundFixed = 0;
235    bool canRoundSolution = true;
236    for (int i=0; i<numberIntegers; i++) {
237      int iColumn = integerVariable[i];
238      double value=newSolution[iColumn];
239      double fraction=value-floor(value);
240      int round = 0;
241      if (fabs(floor(value+0.5)-value)>integerTolerance) {
242        int nDownLocks = downLocks_[i];
243        int nUpLocks = upLocks_[i];
244        if(nDownLocks>0 && nUpLocks>0) {
245          // the variable cannot be rounded
246          canRoundSolution = false;
247          int nLocks = nDownLocks;
248          if(nDownLocks<nUpLocks)
249            round = -1;
250          else if(nDownLocks>nUpLocks) {
251            round = 1;
252            fraction = 1.0 - fraction;
253            nLocks = nUpLocks;
254          }
255          else if(fraction < 0.5)
256            round = -1;
257          else {
258            round = 1;
259            fraction = 1.0 - fraction;
260            nLocks = nUpLocks;
261          }
262
263          // if variable is not binary, penalize it
264          if(!solver->isBinary(iColumn))
265            fraction *= 1000.0;
266
267          if(nLocks < bestLocks || (nLocks == bestLocks &&
268                                    fraction < bestFraction)) {
269            bestColumn = iColumn;
270            bestLocks = nLocks;
271            bestFraction = fraction;
272            bestRound = round;
273          }
274        }
275      }
276      else if(numberAtBoundFixed < maxNumberAtBoundToFix) {
277        // fix the variable at one of its bounds
278        if (fabs(lower[iColumn]-value)<=integerTolerance &&
279            lower[iColumn] != upper[iColumn]) {
280          columnFixed[numberAtBoundFixed] = iColumn;
281          originalBound[numberAtBoundFixed] = upper[iColumn];
282          fixedAtLowerBound[numberAtBoundFixed] = true;
283          solver->setColUpper(iColumn, lower[iColumn]);
284          numberAtBoundFixed++;
285        }
286        else if(fabs(upper[iColumn]-value)<=integerTolerance &&
287            lower[iColumn] != upper[iColumn]) {
288          columnFixed[numberAtBoundFixed] = iColumn;
289          originalBound[numberAtBoundFixed] = lower[iColumn];
290          fixedAtLowerBound[numberAtBoundFixed] = false;
291          solver->setColLower(iColumn, upper[iColumn]);
292          numberAtBoundFixed++;
293        }
294      }
295    }
296
297    if(canRoundSolution) {
298      // Round all the fractional variables
299      for (int i=0; i<numberIntegers; i++) {
300        int iColumn = integerVariable[i];
301        double value=newSolution[iColumn];
302        if (fabs(floor(value+0.5)-value)>integerTolerance) {
303          if(downLocks_[i]==0 || upLocks_[i]==0) {
304            if(downLocks_[i]==0 && upLocks_[i]==0) {
305              if(direction * objective[iColumn] >= 0.0)
306                newSolution[iColumn] = floor(value);
307              else
308                newSolution[iColumn] = ceil(value);
309            }
310            else if(downLocks_[i]==0)
311              newSolution[iColumn] = floor(value);
312            else
313              newSolution[iColumn] = ceil(value);
314          }
315          else
316            break;
317        }
318      }
319      break;
320    }
321
322
323    double originalBoundBestColumn;
324    if(bestColumn >= 0) {
325      if(bestRound < 0) {
326        originalBoundBestColumn = upper[bestColumn];
327        solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
328      }
329      else {
330        originalBoundBestColumn = lower[bestColumn];
331        solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
332      }
333    } else {
334      break;
335    }
336    int originalBestRound = bestRound;
337    while (1) {
338
339      solver->resolve();
340
341      if(!solver->isProvenOptimal()) {
342        if(numberAtBoundFixed > 0) {
343          // Remove the bound fix for variables that were at bounds
344          for(int i=0; i<numberAtBoundFixed; i++) {
345            int iColFixed = columnFixed[i];
346            if(fixedAtLowerBound[i])
347              solver->setColUpper(iColFixed, originalBound[i]);
348            else
349              solver->setColLower(iColFixed, originalBound[i]);
350          }
351          numberAtBoundFixed = 0;
352        }
353        else if(bestRound == originalBestRound) {
354          bestRound *= (-1);
355          if(bestRound < 0) {
356            solver->setColLower(bestColumn, originalBoundBestColumn);
357            solver->setColUpper(bestColumn, floor(newSolution[bestColumn]));
358          }
359          else {
360            solver->setColLower(bestColumn, ceil(newSolution[bestColumn]));
361            solver->setColUpper(bestColumn, originalBoundBestColumn);
362          }
363        }
364        else
365          break;
366      }
367      else
368        break;
369    }
370
371    if(!solver->isProvenOptimal())
372      break;
373
374    if(iteration > maxIterations_) {
375      break;
376    }
377
378    if(CoinCpuTime()-time1 > maxTime_) {
379      break;
380    }
381
382    memcpy(newSolution,solution,numberColumns*sizeof(double));
383    numberFractionalVariables = 0;
384    for (int i=0; i<numberIntegers; i++) {
385      int iColumn = integerVariable[i];
386      double value=newSolution[iColumn];
387      if (fabs(floor(value+0.5)-value)>integerTolerance) {
388        numberFractionalVariables++;
389      }
390    }
391
392  }
393
394  bool updateHowOften = false;
395  // Good question when... every time the heur was run? every time it
396  // has found a solution? every time it has found a better solution?
397  // for now update every time
398  updateHowOften = true;
399
400  double * rowActivity = new double[numberRows];
401  memset(rowActivity,0,numberRows*sizeof(double));
402
403  // re-compute new solution value
404  double objOffset=0.0;
405  solver->getDblParam(OsiObjOffset,objOffset);
406  newSolutionValue = -objOffset;
407  for (int i=0 ; i<numberColumns ; i++ )
408    newSolutionValue += objective[i]*newSolution[i];
409  newSolutionValue *= direction;
410    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
411  if (newSolutionValue<solutionValue) {
412    // paranoid check
413    memset(rowActivity,0,numberRows*sizeof(double));
414    for (int i=0;i<numberColumns;i++) {
415      int j;
416      double value = newSolution[i];
417      if (value) {
418        for (j=columnStart[i];
419             j<columnStart[i]+columnLength[i];j++) {
420          int iRow=row[j];
421          rowActivity[iRow] += value*element[j];
422        }
423      }
424    }
425    // check was approximately feasible
426    bool feasible=true;
427    for (int i=0;i<numberRows;i++) {
428      if(rowActivity[i]<rowLower[i]) {
429        if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
430          feasible = false;
431      } else if(rowActivity[i]>rowUpper[i]) {
432        if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
433          feasible = false;
434      }
435    }
436    for (int i=0; i<numberIntegers; i++) {
437      int iColumn = integerVariable[i];
438      double value=newSolution[iColumn];
439      if (fabs(floor(value+0.5)-value)>integerTolerance) {
440        feasible = false;
441        break;
442      }
443    }
444    if (feasible) {
445      // new solution
446      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
447      solutionValue = newSolutionValue;
448      //printf("** Solution of %g found by CbcHeuristicDiveCoefficient\n",newSolutionValue);
449      returnCode=1;
450
451      //      setCurrentNode(model_->solver());
452    }
453  }
454
455  if (updateHowOften) {
456    howOften_ += (int) (howOften_*decayFactor_);
457  }
458   
459
460  delete [] newSolution;
461  delete [] columnFixed;
462  delete [] originalBound;
463  delete [] fixedAtLowerBound;
464  delete [] rowActivity;
465  delete solver;
466  delete nodeDesc;
467  return returnCode;
468}
469
470// update model
471void CbcHeuristicDiveCoefficient::setModel(CbcModel * model)
472{
473  model_ = model;
474  // Get a copy of original matrix (and by row for rounding);
475  assert(model_->solver());
476  const CoinPackedMatrix * matrix = model_->solver()->getMatrixByCol();
477  if (matrix) {
478    matrix_ = *matrix;
479    // make sure model okay for heuristic
480    validate();
481  }
482}
483
484// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
485void 
486CbcHeuristicDiveCoefficient::validate() 
487{
488  if (model_&&when()<10) {
489    if (model_->numberIntegers()!=
490        model_->numberObjects())
491      setWhen(0);
492  }
493
494  int numberIntegers = model_->numberIntegers();
495  const int * integerVariable = model_->integerVariable();
496  delete [] downLocks_;
497  delete [] upLocks_;
498  downLocks_ = new unsigned short [numberIntegers];
499  upLocks_ = new unsigned short [numberIntegers];
500  // Column copy
501  const double * element = matrix_.getElements();
502  const int * row = matrix_.getIndices();
503  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
504  const int * columnLength = matrix_.getVectorLengths();
505  const double * rowLower = model_->solver()->getRowLower();
506  const double * rowUpper = model_->solver()->getRowUpper();
507  for (int i=0;i<numberIntegers;i++) {
508    int iColumn = integerVariable[i];
509    int down=0;
510    int up=0;
511    if (columnLength[iColumn]>65535) {
512      setWhen(0);
513      break; // unlikely to work
514    }
515    for (CoinBigIndex j=columnStart[iColumn];
516         j<columnStart[iColumn]+columnLength[iColumn];j++) {
517      int iRow=row[j];
518      if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
519        up++;
520        down++;
521      } else if (element[j]>0.0) {
522        if (rowUpper[iRow]<1.0e20)
523          up++;
524        else
525          down++;
526      } else {
527        if (rowLower[iRow]>-1.0e20)
528          up++;
529        else
530          down++;
531      }
532    }
533    downLocks_[i] = (unsigned short) down;
534    upLocks_[i] = (unsigned short) up;
535  }
536}
Note: See TracBrowser for help on using the repository browser.