source: trunk/CbcHeuristic.cpp @ 75

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

add random rounding if zero cost

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.7 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        double randomNumber = CoinDrand48();
198        // which way?
199        if (randomNumber<0.5) 
200          move = below-value;
201        else
202          move = 1.0 -(value-below);
203      }
204      newValue += move;
205      newSolution[iColumn] = newValue;
206      newSolutionValue += move*cost;
207      int j;
208      for (j=columnStart[iColumn];
209           j<columnStart[iColumn]+columnLength[iColumn];j++) {
210        int iRow = row[j];
211        rowActivity[iRow] += move*element[j];
212      }
213    }
214  }
215
216  double penalty=0.0;
217 
218  // see if feasible
219  for (i=0;i<numberRows;i++) {
220    double value = rowActivity[i];
221    double thisInfeasibility=0.0;
222    if (value<rowLower[i]-primalTolerance)
223      thisInfeasibility = value-rowLower[i];
224    else if (value>rowUpper[i]+primalTolerance)
225      thisInfeasibility = value-rowUpper[i];
226    if (thisInfeasibility) {
227      // See if there are any slacks I can use to fix up
228      // maybe put in coding for multiple slacks?
229      double bestCost = 1.0e50;
230      int k;
231      int iBest=-1;
232      double addCost=0.0;
233      double newValue=0.0;
234      double changeRowActivity=0.0;
235      double absInfeasibility = fabs(thisInfeasibility);
236      for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
237        int iColumn = column[k];
238        if (columnLength[iColumn]==1) {
239          double currentValue = newSolution[iColumn];
240          double elementValue = elementByRow[k];
241          double lowerValue = lower[iColumn];
242          double upperValue = upper[iColumn];
243          double gap = rowUpper[i]-rowLower[i];
244          double absElement=fabs(elementValue);
245          if (thisInfeasibility*elementValue>0.0) {
246            // we want to reduce
247            if ((currentValue-lowerValue)*absElement>=absInfeasibility) {
248              // possible - check if integer
249              double distance = absInfeasibility/absElement;
250              double thisCost = -direction*objective[iColumn]*distance;
251              if (solver->isInteger(iColumn)) {
252                distance = ceil(distance-primalTolerance);
253                if (currentValue-distance>=lowerValue-primalTolerance) {
254                  if (absInfeasibility-distance*absElement< -gap-primalTolerance)
255                    thisCost=1.0e100; // no good
256                  else
257                    thisCost = -direction*objective[iColumn]*distance;
258                } else {
259                  thisCost=1.0e100; // no good
260                }
261              }
262              if (thisCost<bestCost) {
263                bestCost=thisCost;
264                iBest=iColumn;
265                addCost = thisCost;
266                newValue = currentValue-distance;
267                changeRowActivity = -distance*elementValue;
268              }
269            }
270          } else {
271            // we want to increase
272            if ((upperValue-currentValue)*absElement>=absInfeasibility) {
273              // possible - check if integer
274              double distance = absInfeasibility/absElement;
275              double thisCost = direction*objective[iColumn]*distance;
276              if (solver->isInteger(iColumn)) {
277                distance = ceil(distance-1.0e-7);
278                assert (currentValue-distance<=upperValue+primalTolerance);
279                if (absInfeasibility-distance*absElement< -gap-primalTolerance)
280                  thisCost=1.0e100; // no good
281                else
282                  thisCost = direction*objective[iColumn]*distance;
283              }
284              if (thisCost<bestCost) {
285                bestCost=thisCost;
286                iBest=iColumn;
287                addCost = thisCost;
288                newValue = currentValue+distance;
289                changeRowActivity = distance*elementValue;
290              }
291            }
292          }
293        }
294      }
295      if (iBest>=0) {
296        /*printf("Infeasibility of %g on row %d cost %g\n",
297          thisInfeasibility,i,addCost);*/
298        newSolution[iBest]=newValue;
299        thisInfeasibility=0.0;
300        newSolutionValue += addCost;
301        rowActivity[i] += changeRowActivity;
302      }
303      penalty += fabs(thisInfeasibility);
304    }
305  }
306
307  // Could also set SOS (using random) and repeat
308  if (!penalty) {
309    // See if we can do better
310    //seed_++;
311    //CoinSeedRandom(seed_);
312    // Random number between 0 and 1.
313    double randomNumber = CoinDrand48();
314    int iPass;
315    int start[2];
316    int end[2];
317    int iRandom = (int) (randomNumber*((double) numberIntegers));
318    start[0]=iRandom;
319    end[0]=numberIntegers;
320    start[1]=0;
321    end[1]=iRandom;
322    for (iPass=0;iPass<2;iPass++) {
323      int i;
324      for (i=start[iPass];i<end[iPass];i++) {
325        int iColumn = integerVariable[i];
326        double value=newSolution[iColumn];
327        assert (fabs(floor(value+0.5)-value)<integerTolerance);
328        double cost = direction * objective[iColumn];
329        double move=0.0;
330        if (cost>0.0)
331          move = -1.0;
332        else if (cost<0.0)
333          move=1.0;
334        while (move) {
335          bool good=true;
336          double newValue=newSolution[iColumn]+move;
337          if (newValue<lower[iColumn]-primalTolerance||
338              newValue>upper[iColumn]+primalTolerance) {
339            move=0.0;
340          } else {
341            // see if we can move
342            int j;
343            for (j=columnStart[iColumn];
344                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
345              int iRow = row[j];
346              double newActivity = rowActivity[iRow] + move*element[j];
347              if (newActivity<rowLower[iRow]-primalTolerance||
348                  newActivity>rowUpper[iRow]+primalTolerance) {
349                good=false;
350                break;
351              }
352            }
353            if (good) {
354              newSolution[iColumn] = newValue;
355              newSolutionValue += move*cost;
356              int j;
357              for (j=columnStart[iColumn];
358                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
359                int iRow = row[j];
360                rowActivity[iRow] += move*element[j];
361              }
362            } else {
363              move=0.0;
364            }
365          }
366        }
367      }
368    }
369    if (newSolutionValue<solutionValue) {
370      // paranoid check
371      memset(rowActivity,0,numberRows*sizeof(double));
372      for (i=0;i<numberColumns;i++) {
373        int j;
374        double value = newSolution[i];
375        if (value) {
376          for (j=columnStart[i];
377               j<columnStart[i]+columnLength[i];j++) {
378            int iRow=row[j];
379            rowActivity[iRow] += value*element[j];
380          }
381        }
382      }
383      // check was approximately feasible
384      bool feasible=true;
385      for (i=0;i<numberRows;i++) {
386        if(rowActivity[i]<rowLower[i]) {
387          if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
388            feasible = false;
389        } else if(rowActivity[i]>rowUpper[i]) {
390          if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
391            feasible = false;
392        }
393      }
394      if (feasible) {
395        // new solution
396        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
397        solutionValue = newSolutionValue;
398        //printf("** Solution of %g found by rounding\n",newSolutionValue);
399        returnCode=1;
400      } else {
401        // Can easily happen
402        //printf("Debug CbcRounding giving bad solution\n");
403      }
404    }
405  }
406  delete [] newSolution;
407  delete [] rowActivity;
408  return returnCode;
409}
410// update model
411void CbcRounding::setModel(CbcModel * model)
412{
413  model_ = model;
414  // Get a copy of original matrix (and by row for rounding);
415  assert(model_->solver());
416  matrix_ = *model_->solver()->getMatrixByCol();
417  matrixByRow_ = *model_->solver()->getMatrixByRow();
418  // make sure model okay for heuristic
419  validate();
420}
421// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
422void 
423CbcRounding::validate() 
424{
425  if (model_&&when()<10) {
426    if (model_->numberIntegers()!=
427        model_->numberObjects())
428      setWhen(0);
429  }
430}
431
432 
Note: See TracBrowser for help on using the repository browser.