source: trunk/CbcHeuristic.cpp @ 73

Last change on this file since 73 was 73, checked in by forrest, 17 years ago

switch off heuristic

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.6 KB
Line 
1// Copyright (C) 2002, 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#include <cassert>
8#include <cmath>
9#include <cfloat>
10
11#include "OsiSolverInterface.hpp"
12#include "CbcModel.hpp"
13#include "CbcMessage.hpp"
14#include "CbcHeuristic.hpp"
15
16// Default Constructor
17CbcHeuristic::CbcHeuristic() 
18  :model_(NULL),
19   when_(2)
20{
21}
22
23// Constructor from model
24CbcHeuristic::CbcHeuristic(CbcModel & model)
25:
26  model_(&model),
27  when_(2)
28{
29}
30// Resets stuff if model changes
31void 
32CbcHeuristic::resetModel(CbcModel * model)
33{
34  model_=model;
35}
36
37// Destructor
38CbcHeuristic::~CbcHeuristic ()
39{
40}
41
42// update model
43void CbcHeuristic::setModel(CbcModel * model)
44{
45  model_ = model;
46}
47
48// Default Constructor
49CbcRounding::CbcRounding() 
50  :CbcHeuristic()
51{
52  // matrix and row copy will automatically be empty
53}
54
55// Constructor from model
56CbcRounding::CbcRounding(CbcModel & model)
57  :CbcHeuristic(model)
58{
59  // Get a copy of original matrix (and by row for rounding);
60  assert(model.solver());
61  matrix_ = *model.solver()->getMatrixByCol();
62  matrixByRow_ = *model.solver()->getMatrixByRow();
63  seed_=1;
64}
65
66// Destructor
67CbcRounding::~CbcRounding ()
68{
69}
70
71// Clone
72CbcHeuristic *
73CbcRounding::clone() const
74{
75  return new CbcRounding(*this);
76}
77
78// Copy constructor
79CbcRounding::CbcRounding(const CbcRounding & rhs)
80:
81  CbcHeuristic(rhs),
82  matrix_(rhs.matrix_),
83  matrixByRow_(rhs.matrixByRow_),
84  seed_(rhs.seed_)
85{
86  setWhen(rhs.when());
87}
88
89// Resets stuff if model changes
90void 
91CbcRounding::resetModel(CbcModel * model)
92{
93  model_=model;
94  // Get a copy of original matrix (and by row for rounding);
95  assert(model_->solver());
96  matrix_ = *model_->solver()->getMatrixByCol();
97  matrixByRow_ = *model_->solver()->getMatrixByRow();
98}
99// See if rounding will give solution
100// Sets value of solution
101// Assumes rhs for original matrix still okay
102// At present only works with integers
103// Fix values if asked for
104// Returns 1 if solution, 0 if not
105int
106CbcRounding::solution(double & solutionValue,
107                         double * betterSolution)
108{
109
110  // See if to do
111  if (!when()||(when()%10==1&&model_->phase()!=1)||
112      (when()%10==2&&model_->phase()!=2))
113    return 0; // switched off
114
115  OsiSolverInterface * solver = model_->solver();
116  const double * lower = solver->getColLower();
117  const double * upper = solver->getColUpper();
118  const double * rowLower = solver->getRowLower();
119  const double * rowUpper = solver->getRowUpper();
120  const double * solution = solver->getColSolution();
121  const double * objective = solver->getObjCoefficients();
122  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
123  double primalTolerance;
124  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
125
126  int numberRows = matrix_.getNumRows();
127
128  int numberIntegers = model_->numberIntegers();
129  const int * integerVariable = model_->integerVariable();
130  int i;
131  double direction = solver->getObjSense();
132  double newSolutionValue = direction*solver->getObjValue();
133  int returnCode = 0;
134
135  // Column copy
136  const double * element = matrix_.getElements();
137  const int * row = matrix_.getIndices();
138  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
139  const int * columnLength = matrix_.getVectorLengths();
140  // Row copy
141  const double * elementByRow = matrixByRow_.getElements();
142  const int * column = matrixByRow_.getIndices();
143  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
144  const int * rowLength = matrixByRow_.getVectorLengths();
145
146  // Get solution array for heuristic solution
147  int numberColumns = solver->getNumCols();
148  double * newSolution = new double [numberColumns];
149  memcpy(newSolution,solution,numberColumns*sizeof(double));
150
151  double * rowActivity = new double[numberRows];
152  memset(rowActivity,0,numberRows*sizeof(double));
153  for (i=0;i<numberColumns;i++) {
154    int j;
155    double value = newSolution[i];
156    if (value<lower[i]) {
157      value=lower[i];
158      newSolution[i]=value;
159    } else if (value>upper[i]) {
160      value=upper[i];
161      newSolution[i]=value;
162    }
163    if (value) {
164      for (j=columnStart[i];
165           j<columnStart[i]+columnLength[i];j++) {
166        int iRow=row[j];
167        rowActivity[iRow] += value*element[j];
168      }
169    }
170  }
171  // check was feasible - if not adjust (cleaning may move)
172  for (i=0;i<numberRows;i++) {
173    if(rowActivity[i]<rowLower[i]) {
174      //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
175      rowActivity[i]=rowLower[i];
176    } else if(rowActivity[i]>rowUpper[i]) {
177      //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
178      rowActivity[i]=rowUpper[i];
179    }
180  }
181  for (i=0;i<numberIntegers;i++) {
182    int iColumn = integerVariable[i];
183    double value=newSolution[iColumn];
184    if (fabs(floor(value+0.5)-value)>integerTolerance) {
185      double below = floor(value);
186      double newValue=newSolution[iColumn];
187      double cost = direction * objective[iColumn];
188      double move;
189      if (cost>0.0) {
190        // try up
191        move = 1.0 -(value-below);
192      } else if (cost<0.0) {
193        // try down
194        move = below-value;
195      } else {
196        // won't be able to move unless we can grab another variable
197        // just for now go down
198        move = below-value;
199      }
200      newValue += move;
201      newSolution[iColumn] = newValue;
202      newSolutionValue += move*cost;
203      int j;
204      for (j=columnStart[iColumn];
205           j<columnStart[iColumn]+columnLength[iColumn];j++) {
206        int iRow = row[j];
207        rowActivity[iRow] += move*element[j];
208      }
209    }
210  }
211
212  double penalty=0.0;
213 
214  // see if feasible
215  for (i=0;i<numberRows;i++) {
216    double value = rowActivity[i];
217    double thisInfeasibility=0.0;
218    if (value<rowLower[i]-primalTolerance)
219      thisInfeasibility = value-rowLower[i];
220    else if (value>rowUpper[i]+primalTolerance)
221      thisInfeasibility = value-rowUpper[i];
222    if (thisInfeasibility) {
223      // See if there are any slacks I can use to fix up
224      // maybe put in coding for multiple slacks?
225      double bestCost = 1.0e50;
226      int k;
227      int iBest=-1;
228      double addCost=0.0;
229      double newValue=0.0;
230      double changeRowActivity=0.0;
231      double absInfeasibility = fabs(thisInfeasibility);
232      for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
233        int iColumn = column[k];
234        if (columnLength[iColumn]==1) {
235          double currentValue = newSolution[iColumn];
236          double elementValue = elementByRow[k];
237          double lowerValue = lower[iColumn];
238          double upperValue = upper[iColumn];
239          double gap = rowUpper[i]-rowLower[i];
240          double absElement=fabs(elementValue);
241          if (thisInfeasibility*elementValue>0.0) {
242            // we want to reduce
243            if ((currentValue-lowerValue)*absElement>=absInfeasibility) {
244              // possible - check if integer
245              double distance = absInfeasibility/absElement;
246              double thisCost = -direction*objective[iColumn]*distance;
247              if (solver->isInteger(iColumn)) {
248                distance = ceil(distance-primalTolerance);
249                if (currentValue-distance>=lowerValue-primalTolerance) {
250                  if (absInfeasibility-distance*absElement< -gap-primalTolerance)
251                    thisCost=1.0e100; // no good
252                  else
253                    thisCost = -direction*objective[iColumn]*distance;
254                } else {
255                  thisCost=1.0e100; // no good
256                }
257              }
258              if (thisCost<bestCost) {
259                bestCost=thisCost;
260                iBest=iColumn;
261                addCost = thisCost;
262                newValue = currentValue-distance;
263                changeRowActivity = -distance*elementValue;
264              }
265            }
266          } else {
267            // we want to increase
268            if ((upperValue-currentValue)*absElement>=absInfeasibility) {
269              // possible - check if integer
270              double distance = absInfeasibility/absElement;
271              double thisCost = direction*objective[iColumn]*distance;
272              if (solver->isInteger(iColumn)) {
273                distance = ceil(distance-1.0e-7);
274                assert (currentValue-distance<=upperValue+primalTolerance);
275                if (absInfeasibility-distance*absElement< -gap-primalTolerance)
276                  thisCost=1.0e100; // no good
277                else
278                  thisCost = direction*objective[iColumn]*distance;
279              }
280              if (thisCost<bestCost) {
281                bestCost=thisCost;
282                iBest=iColumn;
283                addCost = thisCost;
284                newValue = currentValue+distance;
285                changeRowActivity = distance*elementValue;
286              }
287            }
288          }
289        }
290      }
291      if (iBest>=0) {
292        /*printf("Infeasibility of %g on row %d cost %g\n",
293          thisInfeasibility,i,addCost);*/
294        newSolution[iBest]=newValue;
295        thisInfeasibility=0.0;
296        newSolutionValue += addCost;
297        rowActivity[i] += changeRowActivity;
298      }
299      penalty += fabs(thisInfeasibility);
300    }
301  }
302
303  // Could also set SOS (using random) and repeat
304  if (!penalty) {
305    // See if we can do better
306    //seed_++;
307    //CoinSeedRandom(seed_);
308    // Random number between 0 and 1.
309    double randomNumber = CoinDrand48();
310    int iPass;
311    int start[2];
312    int end[2];
313    int iRandom = (int) (randomNumber*((double) numberIntegers));
314    start[0]=iRandom;
315    end[0]=numberIntegers;
316    start[1]=0;
317    end[1]=iRandom;
318    for (iPass=0;iPass<2;iPass++) {
319      int i;
320      for (i=start[iPass];i<end[iPass];i++) {
321        int iColumn = integerVariable[i];
322        double value=newSolution[iColumn];
323        assert (fabs(floor(value+0.5)-value)<integerTolerance);
324        double cost = direction * objective[iColumn];
325        double move=0.0;
326        if (cost>0.0)
327          move = -1.0;
328        else if (cost<0.0)
329          move=1.0;
330        while (move) {
331          bool good=true;
332          double newValue=newSolution[iColumn]+move;
333          if (newValue<lower[iColumn]-primalTolerance||
334              newValue>upper[iColumn]+primalTolerance) {
335            move=0.0;
336          } else {
337            // see if we can move
338            int j;
339            for (j=columnStart[iColumn];
340                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
341              int iRow = row[j];
342              double newActivity = rowActivity[iRow] + move*element[j];
343              if (newActivity<rowLower[iRow]-primalTolerance||
344                  newActivity>rowUpper[iRow]+primalTolerance) {
345                good=false;
346                break;
347              }
348            }
349            if (good) {
350              newSolution[iColumn] = newValue;
351              newSolutionValue += move*cost;
352              int j;
353              for (j=columnStart[iColumn];
354                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
355                int iRow = row[j];
356                rowActivity[iRow] += move*element[j];
357              }
358            } else {
359              move=0.0;
360            }
361          }
362        }
363      }
364    }
365    if (newSolutionValue<solutionValue) {
366      // paranoid check
367      memset(rowActivity,0,numberRows*sizeof(double));
368      for (i=0;i<numberColumns;i++) {
369        int j;
370        double value = newSolution[i];
371        if (value) {
372          for (j=columnStart[i];
373               j<columnStart[i]+columnLength[i];j++) {
374            int iRow=row[j];
375            rowActivity[iRow] += value*element[j];
376          }
377        }
378      }
379      // check was approximately feasible
380      bool feasible=true;
381      for (i=0;i<numberRows;i++) {
382        if(rowActivity[i]<rowLower[i]) {
383          if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
384            feasible = false;
385        } else if(rowActivity[i]>rowUpper[i]) {
386          if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
387            feasible = false;
388        }
389      }
390      if (feasible) {
391        // new solution
392        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
393        solutionValue = newSolutionValue;
394        //printf("** Solution of %g found by rounding\n",newSolutionValue);
395        returnCode=1;
396      } else {
397        // Can easily happen
398        //printf("Debug CbcRounding giving bad solution\n");
399      }
400    }
401  }
402  delete [] newSolution;
403  delete [] rowActivity;
404  return returnCode;
405}
406// update model
407void CbcRounding::setModel(CbcModel * model)
408{
409  model_ = model;
410  // Get a copy of original matrix (and by row for rounding);
411  assert(model_->solver());
412  matrix_ = *model_->solver()->getMatrixByCol();
413  matrixByRow_ = *model_->solver()->getMatrixByRow();
414  // make sure model okay for heuristic
415  validate();
416}
417// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
418void 
419CbcRounding::validate() 
420{
421  if (model_&&when()<10) {
422    if (model_->numberIntegers()!=
423        model_->numberObjects())
424      setWhen(0);
425  }
426}
427
428 
Note: See TracBrowser for help on using the repository browser.