source: trunk/Cbc/src/CbcHeuristicDive.cpp @ 912

Last change on this file since 912 was 912, checked in by ladanyi, 13 years ago

Incorporated changes from branches/heur

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