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

Last change on this file since 581 was 581, checked in by forrest, 13 years ago

fix abort in fpump

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