source: trunk/Samples/CbcHeuristicFPump.cpp @ 144

Last change on this file since 144 was 144, checked in by forrest, 16 years ago

stuffCVS: ----------------------------------------------------------------------

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 10.7 KB
Line 
1// Copyright (C) 2004, 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 "CbcHeuristicFPump.hpp"
15#include "CbcBranchActual.hpp"
16#include "CoinHelperFunctions.hpp"
17#include "CoinTime.hpp"
18
19
20// Default Constructor
21CbcHeuristicFPump::CbcHeuristicFPump() 
22  :CbcHeuristic(),
23   startTime_(0.0),
24   maximumTime_(0.0),
25   maximumPasses_(100),
26   downValue_(0.5),
27   roundExpensive_(false)
28{
29  setWhen(1);
30}
31
32// Constructor from model
33CbcHeuristicFPump::CbcHeuristicFPump(CbcModel & model,
34                                     double downValue,bool roundExpensive)
35  :CbcHeuristic(model),
36   startTime_(0.0),
37   maximumTime_(0.0),
38   maximumPasses_(100),
39   downValue_(downValue),
40   roundExpensive_(roundExpensive)
41{
42  setWhen(1);
43}
44
45// Destructor
46CbcHeuristicFPump::~CbcHeuristicFPump ()
47{
48}
49
50// Clone
51CbcHeuristic *
52CbcHeuristicFPump::clone() const
53{
54  return new CbcHeuristicFPump(*this);
55}
56
57// Copy constructor
58CbcHeuristicFPump::CbcHeuristicFPump(const CbcHeuristicFPump & rhs)
59:
60  CbcHeuristic(rhs),
61  startTime_(rhs.startTime_),
62  maximumTime_(rhs.maximumTime_),
63  maximumPasses_(rhs.maximumPasses_),
64  downValue_(rhs.downValue_),
65  roundExpensive_(rhs.roundExpensive_)
66{
67  setWhen(rhs.when());
68}
69// Resets stuff if model changes
70void 
71CbcHeuristicFPump::resetModel(CbcModel * model)
72{
73}
74
75/**************************BEGIN MAIN PROCEDURE ***********************************/
76
77// See if feasibility pump will give better solution
78// Sets value of solution
79// Returns 1 if solution, 0 if not
80int
81CbcHeuristicFPump::solution(double & solutionValue,
82                         double * betterSolution)
83{
84  if (!when()||(when()==1&&model_->phase()!=1))
85    return 0; // switched off
86  // See if at root node
87  bool atRoot = model_->getNodeCount()==0;
88  int passNumber = model_->getCurrentPassNumber();
89  // just do once
90  if (!atRoot||passNumber!=1)
91    return 0;
92  // probably a good idea
93  if (model_->getSolutionCount()) return 0;
94  // Clone solver - otherwise annoys root node computations
95  OsiSolverInterface * solver = model_->solver()->clone();
96  solver->resolve();
97  const double * lower = solver->getColLower();
98  const double * upper = solver->getColUpper();
99  const double * solution = solver->getColSolution();
100  double primalTolerance;
101  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
102 
103  int numberColumns = model_->getNumCols();
104  int numberIntegers = model_->numberIntegers();
105  const int * integerVariable = model_->integerVariable();
106
107// 1. initially check 0-1
108  int i,j;
109  bool zeroOne=true;
110  for (i=0;i<numberIntegers;i++) {
111    int iColumn = integerVariable[i];
112    const CbcObject * object = model_->object(i);
113    const CbcSimpleInteger * integerObject = 
114      dynamic_cast<const  CbcSimpleInteger *> (object);
115    assert(integerObject);
116    if (upper[iColumn]-lower[iColumn]>1.000001) {
117      zeroOne=false;
118      break;
119    }
120  }
121  if (!zeroOne) {
122    delete solver;
123    return 0;
124  }
125
126// 2. space for rounded solution
127  double * newSolution = new double [numberColumns];
128  // space for last rounded solutions
129#define NUMBER_OLD 4
130  double ** oldSolution = new double * [NUMBER_OLD];
131  for (j=0;j<NUMBER_OLD;j++) {
132    oldSolution[j]= new double[numberColumns];
133    for (i=0;i<numberColumns;i++) oldSolution[j][i]=-COIN_DBL_MAX;
134  }
135
136// 3. Replace objective with an initial 0-valued objective
137  double * saveObjective = new double [numberColumns];
138  memcpy(saveObjective,solver->getObjCoefficients(),numberColumns*sizeof(double));
139  for (i=0;i<numberColumns;i++) {
140    solver->setObjCoeff(i,0.0);
141  }
142  bool finished=false;
143  double direction = solver->getObjSense();
144  int returnCode=0;
145  bool takeHint;
146  OsiHintStrength strength;
147  solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
148  solver->setHintParam(OsiDoDualInResolve,false);
149  solver->messageHandler()->setLogLevel(1);
150
151// 4. Save objective offset so we can see progress
152  double saveOffset;
153  solver->getDblParam(OsiObjOffset,saveOffset);
154
155// 5. MAIN WHILE LOOP
156  int numberPasses=0;
157  while (!finished) {
158    if (numberPasses>=maximumPasses_) {
159      break;
160    }
161    if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
162    numberPasses++;
163    memcpy(newSolution,solution,numberColumns*sizeof(double));
164    int flip;
165    returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip);
166    if (returnCode) {
167      // SOLUTION IS INTEGER
168      // Put back correct objective
169      printf("\n");
170      for (i=0;i<numberColumns;i++)
171        solver->setObjCoeff(i,saveObjective[i]);
172      // solution - but may not be better
173      // Compute using dot product
174      double newSolutionValue = direction*solver->OsiSolverInterface::getObjValue();
175      if (newSolutionValue<solutionValue) {
176        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
177      } else {
178        returnCode=0;
179      }     
180      break;
181    } else {
182      // SOLUTION IS not INTEGER
183      // 1. check for loop
184      bool matched;
185      for (int k = NUMBER_OLD-1; k > 0; k--) {
186          double * b = oldSolution[k];
187          matched = true;
188          for (i = 0; i <numberIntegers; i++) {
189              int iColumn = integerVariable[i];
190              if (newSolution[iColumn]!=b[iColumn]) {
191                matched=false;
192                break;
193              }
194          }
195          if (matched) break;
196      }
197      if (matched || numberPasses%100 == 0) {
198         // perturbation
199         printf("Perturbation applied");
200         for (i=0;i<numberIntegers;i++) {
201             int iColumn = integerVariable[i];
202             double value = max(0.0,CoinDrand48()-0.3);
203             double difference = fabs(solution[iColumn]-newSolution[iColumn]);
204             if (difference+value>0.5) {
205                if (newSolution[iColumn]<lower[iColumn]+primalTolerance) newSolution[iColumn] += 1.0;
206             else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) newSolution[iColumn] -= 1.0;
207                  else abort();
208             }
209         }
210      } else {
211         for (j=NUMBER_OLD-1;j>0;j--) {
212             for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
213         }
214         for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
215      }
216
217      // 2. update the objective function based on the new rounded solution
218      double offset=0.0;
219      for (i=0;i<numberIntegers;i++) {
220        int iColumn = integerVariable[i];
221        double costValue = 1.0;
222        // deal with fixed variables (i.e., upper=lower)
223        if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) {
224           if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
225           else                                       solver->setObjCoeff(iColumn,costValue);
226           continue;
227        }
228        if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
229          solver->setObjCoeff(iColumn,costValue);
230        } else {
231          if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
232            solver->setObjCoeff(iColumn,-costValue);
233          } else {
234            abort();
235          }
236        }
237        offset += costValue*newSolution[iColumn];
238      }
239      solver->setDblParam(OsiObjOffset,-offset);
240      solver->resolve();
241      printf("\npass %3d: obj. %10.5lf --> ", numberPasses,solver->getObjValue());
242
243
244    }
245  } // END WHILE
246
247  delete solver;
248  delete [] newSolution;
249  for ( j=0;j<NUMBER_OLD;j++) 
250    delete [] oldSolution[j];
251  delete [] oldSolution;
252  delete [] saveObjective;
253  return returnCode;
254}
255
256/**************************END MAIN PROCEDURE ***********************************/
257
258// update model
259void CbcHeuristicFPump::setModel(CbcModel * model)
260{
261  model_ = model;
262}
263
264/* Rounds solution - down if < downValue
265   returns 1 if current is a feasible solution
266*/
267int 
268CbcHeuristicFPump::rounds(double * solution,
269                          const double * objective,
270                          bool roundExpensive, double downValue, int *flip)
271{
272  OsiSolverInterface * solver = model_->solver();
273  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
274  double primalTolerance ;
275  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
276  int numberIntegers = model_->numberIntegers();
277  const int * integerVariable = model_->integerVariable();
278
279  int i;
280
281  static int iter = 0;
282  int numberColumns = model_->getNumCols();
283  // tmp contains the current obj coefficients
284  double * tmp = new double [numberColumns];
285  memcpy(tmp,solver->getObjCoefficients(),numberColumns*sizeof(double));
286  int flip_up = 0;
287  int flip_down  = 0;
288  double  v = CoinDrand48() * 20;
289  int nn = 10 + (int) v;
290  int nnv = 0;
291  int * list = new int [nn];
292  double * val = new double [nn];
293  for (i = 0; i < nn; i++) val[i] = .001;
294
295  // return rounded solution
296  for (i=0;i<numberIntegers;i++) {
297    int iColumn = integerVariable[i];
298    const CbcObject * object = model_->object(i);
299    const CbcSimpleInteger * integerObject = 
300      dynamic_cast<const  CbcSimpleInteger *> (object);
301    assert(integerObject);
302    double value=solution[iColumn];
303    double round = floor(value+primalTolerance);
304    if (value-round > .5) round += 1.;
305    if (round < integerTolerance && tmp[iColumn] < -1. + integerTolerance) flip_down++;
306    if (round > 1. - integerTolerance && tmp[iColumn] > 1. - integerTolerance) flip_up++;
307    if (flip_up + flip_down == 0) { 
308       for (int k = 0; k < nn; k++) {
309           if (fabs(value-round) > val[k]) {
310              nnv++;
311              for (int j = nn-2; j >= k; j--) {
312                  val[j+1] = val[j];
313                  list[j+1] = list[j];
314              } 
315              val[k] = fabs(value-round);
316              list[k] = iColumn;
317              break;
318           }
319       }
320    }
321    solution[iColumn] = round;
322  }
323
324  if (nnv > nn) nnv = nn;
325  if (iter != 0) printf("up = %5d , down = %5d", flip_up, flip_down); fflush(stdout);
326  *flip = flip_up + flip_down;
327  delete [] tmp;
328
329  if (*flip == 0 && iter != 0) {
330     printf(" -- rand = %4d (%4d) ", nnv, nn);
331     for (i = 0; i < nnv; i++) solution[list[i]] = 1. - solution[list[i]];
332     *flip = nnv;
333  } else printf(" ");
334  delete [] list; delete [] val;
335  iter++;
336   
337  const double * rowLower = solver->getRowLower();
338  const double * rowUpper = solver->getRowUpper();
339
340  int numberRows = solver->getNumRows();
341  // get row activities
342  double * rowActivity = new double[numberRows];
343  memset(rowActivity,0,numberRows*sizeof(double));
344  solver->getMatrixByCol()->times(solution,rowActivity) ;
345  double largestInfeasibility =0.0;
346  for (i=0 ; i < numberRows ; i++) {
347    largestInfeasibility = max(largestInfeasibility,
348                               rowLower[i]-rowActivity[i]);
349    largestInfeasibility = max(largestInfeasibility,
350                               rowActivity[i]-rowUpper[i]);
351  }
352  return (largestInfeasibility>primalTolerance) ? 0 : 1;
353}
354// Set maximum Time (default off) - also sets starttime to current
355void 
356CbcHeuristicFPump::setMaximumTime(double value)
357{
358  startTime_=CoinCpuTime();
359  maximumTime_=value;
360}
361
362 
Note: See TracBrowser for help on using the repository browser.