source: trunk/Cbc/src/CbcHeuristicFPump.cpp @ 310

Last change on this file since 310 was 310, checked in by andreasw, 13 years ago

first commit for autotools conversion to be able to move more files

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 12.2 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->setDblParam(OsiDualObjectiveLimit,1.0e50);
97  solver->resolve();
98  const double * lower = solver->getColLower();
99  const double * upper = solver->getColUpper();
100  const double * solution = solver->getColSolution();
101  double primalTolerance;
102  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
103 
104  int numberColumns = model_->getNumCols();
105  int numberIntegers = model_->numberIntegers();
106  const int * integerVariable = model_->integerVariable();
107
108// 1. initially check 0-1
109  int i,j;
110  int general=0;
111  for (i=0;i<numberIntegers;i++) {
112    int iColumn = integerVariable[i];
113    const CbcObject * object = model_->object(i);
114#ifndef NDEBUG
115    const CbcSimpleInteger * integerObject = 
116      dynamic_cast<const  CbcSimpleInteger *> (object);
117    assert(integerObject);
118#endif
119    if (upper[iColumn]-lower[iColumn]>1.000001) {
120      general++;
121      break;
122    }
123  }
124  if (general*3>numberIntegers) {
125    delete solver;
126    return 0;
127  }
128
129// 2. space for rounded solution
130  double * newSolution = new double [numberColumns];
131  // space for last rounded solutions
132#define NUMBER_OLD 4
133  double ** oldSolution = new double * [NUMBER_OLD];
134  for (j=0;j<NUMBER_OLD;j++) {
135    oldSolution[j]= new double[numberColumns];
136    for (i=0;i<numberColumns;i++) oldSolution[j][i]=-COIN_DBL_MAX;
137  }
138
139// 3. Replace objective with an initial 0-valued objective
140  double * saveObjective = new double [numberColumns];
141  memcpy(saveObjective,solver->getObjCoefficients(),numberColumns*sizeof(double));
142  for (i=0;i<numberColumns;i++) {
143    solver->setObjCoeff(i,0.0);
144  }
145  bool finished=false;
146  double direction = solver->getObjSense();
147  int returnCode=0;
148  bool takeHint;
149  OsiHintStrength strength;
150  solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
151  solver->setHintParam(OsiDoDualInResolve,false);
152  solver->messageHandler()->setLogLevel(0);
153
154// 4. Save objective offset so we can see progress
155  double saveOffset;
156  solver->getDblParam(OsiObjOffset,saveOffset);
157
158// 5. MAIN WHILE LOOP
159  int numberPasses=0;
160  bool newLineNeeded=false;
161  while (!finished) {
162    returnCode=0;
163    if (numberPasses>=maximumPasses_) {
164      break;
165    }
166    if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
167    numberPasses++;
168    memcpy(newSolution,solution,numberColumns*sizeof(double));
169    int flip;
170    returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip);
171    if (returnCode) {
172      // SOLUTION IS INTEGER
173      // Put back correct objective
174      for (i=0;i<numberColumns;i++)
175        solver->setObjCoeff(i,saveObjective[i]);
176      // solution - but may not be better
177      // Compute using dot product
178      solver->setDblParam(OsiObjOffset,saveOffset);
179      double newSolutionValue = direction*solver->OsiSolverInterface::getObjValue();
180      if (model_->logLevel())
181        printf(" - solution found\n");
182      newLineNeeded=false;
183      if (newSolutionValue<solutionValue) {
184        if (general) {
185          int numberLeft=0;
186          for (i=0;i<numberIntegers;i++) {
187            int iColumn = integerVariable[i];
188            double value = floor(newSolution[iColumn]+0.5);
189            if(solver->isBinary(iColumn)) {
190              solver->setColLower(iColumn,value);
191              solver->setColUpper(iColumn,value);
192            } else {
193              if (fabs(value-newSolution[iColumn])>1.0e-7) 
194                numberLeft++;
195            }
196          }
197          if (numberLeft) {
198            returnCode = smallBranchAndBound(solver,200,newSolution,newSolutionValue,
199                                         solutionValue,"CbcHeuristicFpump");
200          }
201        }
202        if (returnCode) {
203          memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
204          solutionValue=newSolutionValue;
205        }
206      } else {
207        returnCode=0;
208      }     
209      break;
210    } else {
211      // SOLUTION IS not INTEGER
212      // 1. check for loop
213      bool matched;
214      for (int k = NUMBER_OLD-1; k > 0; k--) {
215          double * b = oldSolution[k];
216          matched = true;
217          for (i = 0; i <numberIntegers; i++) {
218              int iColumn = integerVariable[i];
219              if(!solver->isBinary(iColumn))
220                continue;
221              if (newSolution[iColumn]!=b[iColumn]) {
222                matched=false;
223                break;
224              }
225          }
226          if (matched) break;
227      }
228      if (matched || numberPasses%100 == 0) {
229         // perturbation
230        if (model_->logLevel())
231          printf("Perturbation applied");
232         newLineNeeded=true;
233         for (i=0;i<numberIntegers;i++) {
234             int iColumn = integerVariable[i];
235             if(!solver->isBinary(iColumn))
236               continue;
237             double value = max(0.0,CoinDrand48()-0.3);
238             double difference = fabs(solution[iColumn]-newSolution[iColumn]);
239             if (difference+value>0.5) {
240                if (newSolution[iColumn]<lower[iColumn]+primalTolerance) newSolution[iColumn] += 1.0;
241             else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) newSolution[iColumn] -= 1.0;
242                  else abort();
243             }
244         }
245      } else {
246         for (j=NUMBER_OLD-1;j>0;j--) {
247             for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
248         }
249         for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
250      }
251
252      // 2. update the objective function based on the new rounded solution
253      double offset=0.0;
254      for (i=0;i<numberIntegers;i++) {
255        int iColumn = integerVariable[i];
256        if(!solver->isBinary(iColumn))
257          continue;
258        double costValue = 1.0;
259        // deal with fixed variables (i.e., upper=lower)
260        if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) {
261           if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
262           else                                       solver->setObjCoeff(iColumn,costValue);
263           continue;
264        }
265        if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
266          solver->setObjCoeff(iColumn,costValue);
267        } else {
268          if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
269            solver->setObjCoeff(iColumn,-costValue);
270          } else {
271            abort();
272          }
273        }
274        offset += costValue*newSolution[iColumn];
275      }
276      solver->setDblParam(OsiObjOffset,-offset);
277      solver->resolve();
278      if (model_->logLevel())
279        printf("\npass %3d: obj. %10.5lf --> ", numberPasses,solver->getObjValue());
280      newLineNeeded=true;
281
282    }
283  } // END WHILE
284  if (newLineNeeded&&model_->logLevel())
285    printf(" - no solution found\n");
286  delete solver;
287  delete [] newSolution;
288  for ( j=0;j<NUMBER_OLD;j++) 
289    delete [] oldSolution[j];
290  delete [] oldSolution;
291  delete [] saveObjective;
292  return returnCode;
293}
294
295/**************************END MAIN PROCEDURE ***********************************/
296
297// update model
298void CbcHeuristicFPump::setModel(CbcModel * model)
299{
300  model_ = model;
301}
302
303/* Rounds solution - down if < downValue
304   returns 1 if current is a feasible solution
305*/
306int 
307CbcHeuristicFPump::rounds(double * solution,
308                          const double * objective,
309                          bool roundExpensive, double downValue, int *flip)
310{
311  OsiSolverInterface * solver = model_->solver();
312  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
313  double primalTolerance ;
314  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
315  int numberIntegers = model_->numberIntegers();
316  const int * integerVariable = model_->integerVariable();
317
318  int i;
319
320  static int iter = 0;
321  int numberColumns = model_->getNumCols();
322  // tmp contains the current obj coefficients
323  double * tmp = new double [numberColumns];
324  memcpy(tmp,solver->getObjCoefficients(),numberColumns*sizeof(double));
325  int flip_up = 0;
326  int flip_down  = 0;
327  double  v = CoinDrand48() * 20;
328  int nn = 10 + (int) v;
329  int nnv = 0;
330  int * list = new int [nn];
331  double * val = new double [nn];
332  for (i = 0; i < nn; i++) val[i] = .001;
333
334  // return rounded solution
335  for (i=0;i<numberIntegers;i++) {
336    int iColumn = integerVariable[i];
337    if(!solver->isBinary(iColumn))
338      continue;
339    const CbcObject * object = model_->object(i);
340#ifndef NDEBUG
341    const CbcSimpleInteger * integerObject = 
342      dynamic_cast<const  CbcSimpleInteger *> (object);
343    assert(integerObject);
344#endif
345    double value=solution[iColumn];
346    double round = floor(value+primalTolerance);
347    if (value-round > .5) round += 1.;
348    if (round < integerTolerance && tmp[iColumn] < -1. + integerTolerance) flip_down++;
349    if (round > 1. - integerTolerance && tmp[iColumn] > 1. - integerTolerance) flip_up++;
350    if (flip_up + flip_down == 0) { 
351       for (int k = 0; k < nn; k++) {
352           if (fabs(value-round) > val[k]) {
353              nnv++;
354              for (int j = nn-2; j >= k; j--) {
355                  val[j+1] = val[j];
356                  list[j+1] = list[j];
357              } 
358              val[k] = fabs(value-round);
359              list[k] = iColumn;
360              break;
361           }
362       }
363    }
364    solution[iColumn] = round;
365  }
366
367  if (nnv > nn) nnv = nn;
368  if (iter != 0&&model_->logLevel())
369    printf("up = %5d , down = %5d", flip_up, flip_down); fflush(stdout);
370  *flip = flip_up + flip_down;
371  delete [] tmp;
372
373  if (*flip == 0 && iter != 0) {
374    if(model_->logLevel())
375      printf(" -- rand = %4d (%4d) ", nnv, nn);
376     for (i = 0; i < nnv; i++) solution[list[i]] = 1. - solution[list[i]];
377     *flip = nnv;
378  } else if (model_->logLevel()) {
379    printf(" ");
380  }
381  delete [] list; delete [] val;
382  iter++;
383   
384  const double * rowLower = solver->getRowLower();
385  const double * rowUpper = solver->getRowUpper();
386
387  int numberRows = solver->getNumRows();
388  // get row activities
389  double * rowActivity = new double[numberRows];
390  memset(rowActivity,0,numberRows*sizeof(double));
391  solver->getMatrixByCol()->times(solution,rowActivity) ;
392  double largestInfeasibility =0.0;
393  for (i=0 ; i < numberRows ; i++) {
394    largestInfeasibility = max(largestInfeasibility,
395                               rowLower[i]-rowActivity[i]);
396    largestInfeasibility = max(largestInfeasibility,
397                               rowActivity[i]-rowUpper[i]);
398  }
399  delete [] rowActivity;
400  return (largestInfeasibility>primalTolerance) ? 0 : 1;
401}
402// Set maximum Time (default off) - also sets starttime to current
403void 
404CbcHeuristicFPump::setMaximumTime(double value)
405{
406  startTime_=CoinCpuTime();
407  maximumTime_=value;
408}
409
410 
Note: See TracBrowser for help on using the repository browser.