source: trunk/CbcHeuristic.cpp @ 2

Last change on this file since 2 was 2, checked in by ladanyi, 15 years ago

Import of Coin Branch-and-Cut (formerly known as Sbb)

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