source: branches/devel/Cbc/src/CbcHeuristicFPump.cpp @ 435

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

extend fpump heuristic and allow more probing

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.1 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  char * usedColumn = NULL;
123  double * lastSolution=NULL;
124  bool fixContinuous=false;
125  bool fixInternal=false;
126  if (when_>=11&&when_<=13) {
127    fixInternal = when_ >11;
128    fixContinuous = when_==13;
129    when_=1;
130    usedColumn = new char [numberColumns];
131    memset(usedColumn,0,numberColumns);
132    lastSolution = CoinCopyOfArray(solver->getColSolution(),numberColumns);
133  }
134  int numberIntegers = model_->numberIntegers();
135  const int * integerVariable = model_->integerVariable();
136
137// 1. initially check 0-1
138  int i,j;
139  int general=0;
140  for (i=0;i<numberIntegers;i++) {
141    int iColumn = integerVariable[i];
142#ifndef NDEBUG
143    const CbcObject * object = model_->object(i);
144    const CbcSimpleInteger * integerObject = 
145      dynamic_cast<const  CbcSimpleInteger *> (object);
146    assert(integerObject);
147#endif
148    if (upper[iColumn]-lower[iColumn]>1.000001) {
149      general++;
150      break;
151    }
152  }
153  if (general*3>numberIntegers) {
154    delete solver;
155    return 0;
156  }
157
158// 2. space for rounded solution
159  double * newSolution = new double [numberColumns];
160  // space for last rounded solutions
161#define NUMBER_OLD 4
162  double ** oldSolution = new double * [NUMBER_OLD];
163  for (j=0;j<NUMBER_OLD;j++) {
164    oldSolution[j]= new double[numberColumns];
165    for (i=0;i<numberColumns;i++) oldSolution[j][i]=-COIN_DBL_MAX;
166  }
167
168// 3. Replace objective with an initial 0-valued objective
169  double * saveObjective = new double [numberColumns];
170  memcpy(saveObjective,solver->getObjCoefficients(),numberColumns*sizeof(double));
171  for (i=0;i<numberColumns;i++) {
172    solver->setObjCoeff(i,0.0);
173  }
174  bool finished=false;
175  double direction = solver->getObjSense();
176  int returnCode=0;
177  bool takeHint;
178  OsiHintStrength strength;
179  solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
180  solver->setHintParam(OsiDoDualInResolve,false);
181  solver->messageHandler()->setLogLevel(0);
182
183// 4. Save objective offset so we can see progress
184  double saveOffset;
185  solver->getDblParam(OsiObjOffset,saveOffset);
186
187// 5. MAIN WHILE LOOP
188  int numberPasses=0;
189  double newSolutionValue=COIN_DBL_MAX;
190  bool newLineNeeded=false;
191  while (!finished) {
192    returnCode=0;
193    // see what changed
194    if (usedColumn) {
195      const double * solution = solver->getColSolution();
196      for (i=0;i<numberColumns;i++) {
197        if (fabs(solution[i]-lastSolution[i])>1.0e-8) 
198          usedColumn[i]=1;
199        lastSolution[i]=solution[i];
200      }
201    }
202    if (numberPasses>=maximumPasses_) {
203      break;
204    }
205    if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
206    numberPasses++;
207    memcpy(newSolution,solution,numberColumns*sizeof(double));
208    int flip;
209    returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip);
210    if (returnCode) {
211      // SOLUTION IS INTEGER
212      // Put back correct objective
213      for (i=0;i<numberColumns;i++)
214        solver->setObjCoeff(i,saveObjective[i]);
215      // solution - but may not be better
216      // Compute using dot product
217      solver->setDblParam(OsiObjOffset,saveOffset);
218      newSolutionValue = direction*solver->OsiSolverInterface::getObjValue();
219      if (model_->logLevel())
220        printf(" - solution found\n");
221      newLineNeeded=false;
222      if (newSolutionValue<solutionValue) {
223        if (general) {
224          int numberLeft=0;
225          for (i=0;i<numberIntegers;i++) {
226            int iColumn = integerVariable[i];
227            double value = floor(newSolution[iColumn]+0.5);
228            if(solver->isBinary(iColumn)) {
229              solver->setColLower(iColumn,value);
230              solver->setColUpper(iColumn,value);
231            } else {
232              if (fabs(value-newSolution[iColumn])>1.0e-7) 
233                numberLeft++;
234            }
235          }
236          if (numberLeft) {
237            returnCode = smallBranchAndBound(solver,200,newSolution,newSolutionValue,
238                                         solutionValue,"CbcHeuristicFpump");
239          }
240        }
241        if (returnCode) {
242          memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
243          solutionValue=newSolutionValue;
244        }
245      } else {
246        returnCode=0;
247      }     
248      break;
249    } else {
250      // SOLUTION IS not INTEGER
251      // 1. check for loop
252      bool matched;
253      for (int k = NUMBER_OLD-1; k > 0; k--) {
254          double * b = oldSolution[k];
255          matched = true;
256          for (i = 0; i <numberIntegers; i++) {
257              int iColumn = integerVariable[i];
258              if(!solver->isBinary(iColumn))
259                continue;
260              if (newSolution[iColumn]!=b[iColumn]) {
261                matched=false;
262                break;
263              }
264          }
265          if (matched) break;
266      }
267      if (matched || numberPasses%100 == 0) {
268         // perturbation
269        if (model_->logLevel())
270          printf("Perturbation applied");
271         newLineNeeded=true;
272         for (i=0;i<numberIntegers;i++) {
273             int iColumn = integerVariable[i];
274             if(!solver->isBinary(iColumn))
275               continue;
276             double value = max(0.0,CoinDrand48()-0.3);
277             double difference = fabs(solution[iColumn]-newSolution[iColumn]);
278             if (difference+value>0.5) {
279                if (newSolution[iColumn]<lower[iColumn]+primalTolerance) newSolution[iColumn] += 1.0;
280             else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) newSolution[iColumn] -= 1.0;
281                  else abort();
282             }
283         }
284      } else {
285         for (j=NUMBER_OLD-1;j>0;j--) {
286             for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
287         }
288         for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
289      }
290
291      // 2. update the objective function based on the new rounded solution
292      double offset=0.0;
293      for (i=0;i<numberIntegers;i++) {
294        int iColumn = integerVariable[i];
295        if(!solver->isBinary(iColumn))
296          continue;
297        double costValue = 1.0;
298        // deal with fixed variables (i.e., upper=lower)
299        if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) {
300           if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
301           else                                       solver->setObjCoeff(iColumn,costValue);
302           continue;
303        }
304        if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
305          solver->setObjCoeff(iColumn,costValue);
306        } else {
307          if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
308            solver->setObjCoeff(iColumn,-costValue);
309          } else {
310            abort();
311          }
312        }
313        offset += costValue*newSolution[iColumn];
314      }
315      solver->setDblParam(OsiObjOffset,-offset);
316      solver->resolve();
317      if (model_->logLevel())
318        printf("\npass %3d: obj. %10.5f --> ", numberPasses,solver->getObjValue());
319      newLineNeeded=true;
320
321    }
322  } // END WHILE
323  if (newLineNeeded&&model_->logLevel())
324    printf(" - no solution found\n");
325  delete solver;
326  for ( j=0;j<NUMBER_OLD;j++) 
327    delete [] oldSolution[j];
328  delete [] oldSolution;
329  delete [] saveObjective;
330  if (usedColumn) {
331    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
332    const double * colLower = newSolver->getColLower();
333    const double * colUpper = newSolver->getColUpper();
334    int i;
335    int nFix=0;
336    int nFixI=0;
337    int nFixC=0;
338    for (i=0;i<numberIntegers;i++) {
339      int iColumn=integerVariable[i];
340      const CbcObject * object = model_->object(i);
341      const CbcSimpleInteger * integerObject = 
342        dynamic_cast<const  CbcSimpleInteger *> (object);
343      assert(integerObject);
344      // get original bounds
345      double originalLower = integerObject->originalLowerBound();
346      assert(colLower[iColumn]==originalLower);
347      newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
348      double originalUpper = integerObject->originalUpperBound();
349      assert(colUpper[iColumn]==originalUpper);
350      newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
351      if (!usedColumn[iColumn]) {
352        double value=lastSolution[iColumn];
353        double nearest = floor(value+0.5);
354        if (fabs(value-nearest)<1.0e-7) {
355          if (nearest==colLower[iColumn]) {
356            newSolver->setColUpper(iColumn,colLower[iColumn]);
357            nFix++;
358          } else if (nearest==colUpper[iColumn]) {
359            newSolver->setColLower(iColumn,colUpper[iColumn]);
360            nFix++;
361          } else if (fixInternal) {
362            newSolver->setColLower(iColumn,nearest);
363            newSolver->setColUpper(iColumn,nearest);
364            nFix++;
365            nFixI++;
366          }
367        }
368      }
369    }
370    if (fixContinuous) {
371      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
372        if (!newSolver->isInteger(iColumn)&&!usedColumn[iColumn]) {
373          double value=lastSolution[iColumn];
374          if (value<colLower[iColumn]+1.0e-8) {
375            newSolver->setColUpper(iColumn,colLower[iColumn]);
376            nFix++;
377            nFixC++;
378          } else if (value>colUpper[iColumn]-1.0e-8) {
379            newSolver->setColLower(iColumn,colUpper[iColumn]);
380            nFix++;
381            nFixC++;
382          }
383        }
384      }
385    }
386    newSolver->initialSolve();
387    assert (newSolver->isProvenOptimal());
388    printf("%d integers at bound fixed, %d internal and %d continuous\n",
389           nFix,nFixI,nFixC);
390    double saveValue = newSolutionValue;
391    returnCode = smallBranchAndBound(newSolver,200,newSolution,newSolutionValue,
392                                     newSolutionValue,"CbcHeuristicLocalAfterFPump");
393    if (returnCode) {
394      printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
395      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
396      solutionValue=newSolutionValue;
397    }
398    delete newSolver;
399    delete [] usedColumn;
400    delete [] lastSolution;
401  }
402  delete [] newSolution;
403  return returnCode;
404}
405
406/**************************END MAIN PROCEDURE ***********************************/
407
408// update model
409void CbcHeuristicFPump::setModel(CbcModel * model)
410{
411  model_ = model;
412}
413
414/* Rounds solution - down if < downValue
415   returns 1 if current is a feasible solution
416*/
417int 
418CbcHeuristicFPump::rounds(double * solution,
419                          const double * objective,
420                          bool roundExpensive, double downValue, int *flip)
421{
422  OsiSolverInterface * solver = model_->solver();
423  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
424  double primalTolerance ;
425  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
426  int numberIntegers = model_->numberIntegers();
427  const int * integerVariable = model_->integerVariable();
428
429  int i;
430
431  static int iter = 0;
432  int numberColumns = model_->getNumCols();
433  // tmp contains the current obj coefficients
434  double * tmp = new double [numberColumns];
435  memcpy(tmp,solver->getObjCoefficients(),numberColumns*sizeof(double));
436  int flip_up = 0;
437  int flip_down  = 0;
438  double  v = CoinDrand48() * 20;
439  int nn = 10 + (int) v;
440  int nnv = 0;
441  int * list = new int [nn];
442  double * val = new double [nn];
443  for (i = 0; i < nn; i++) val[i] = .001;
444
445  // return rounded solution
446  for (i=0;i<numberIntegers;i++) {
447    int iColumn = integerVariable[i];
448    if(!solver->isBinary(iColumn))
449      continue;
450#ifndef NDEBUG
451    const CbcObject * object = model_->object(i);
452    const CbcSimpleInteger * integerObject = 
453      dynamic_cast<const  CbcSimpleInteger *> (object);
454    assert(integerObject);
455#endif
456    double value=solution[iColumn];
457    double round = floor(value+primalTolerance);
458    if (value-round > .5) round += 1.;
459    if (round < integerTolerance && tmp[iColumn] < -1. + integerTolerance) flip_down++;
460    if (round > 1. - integerTolerance && tmp[iColumn] > 1. - integerTolerance) flip_up++;
461    if (flip_up + flip_down == 0) { 
462       for (int k = 0; k < nn; k++) {
463           if (fabs(value-round) > val[k]) {
464              nnv++;
465              for (int j = nn-2; j >= k; j--) {
466                  val[j+1] = val[j];
467                  list[j+1] = list[j];
468              } 
469              val[k] = fabs(value-round);
470              list[k] = iColumn;
471              break;
472           }
473       }
474    }
475    solution[iColumn] = round;
476  }
477
478  if (nnv > nn) nnv = nn;
479  if (iter != 0&&model_->logLevel())
480    printf("up = %5d , down = %5d", flip_up, flip_down); fflush(stdout);
481  *flip = flip_up + flip_down;
482  delete [] tmp;
483
484  if (*flip == 0 && iter != 0) {
485    if(model_->logLevel())
486      printf(" -- rand = %4d (%4d) ", nnv, nn);
487     for (i = 0; i < nnv; i++) solution[list[i]] = 1. - solution[list[i]];
488     *flip = nnv;
489  } else if (model_->logLevel()) {
490    printf(" ");
491  }
492  delete [] list; delete [] val;
493  iter++;
494   
495  const double * rowLower = solver->getRowLower();
496  const double * rowUpper = solver->getRowUpper();
497
498  int numberRows = solver->getNumRows();
499  // get row activities
500  double * rowActivity = new double[numberRows];
501  memset(rowActivity,0,numberRows*sizeof(double));
502  solver->getMatrixByCol()->times(solution,rowActivity) ;
503  double largestInfeasibility =0.0;
504  for (i=0 ; i < numberRows ; i++) {
505    largestInfeasibility = max(largestInfeasibility,
506                               rowLower[i]-rowActivity[i]);
507    largestInfeasibility = max(largestInfeasibility,
508                               rowActivity[i]-rowUpper[i]);
509  }
510  delete [] rowActivity;
511  return (largestInfeasibility>primalTolerance) ? 0 : 1;
512}
513// Set maximum Time (default off) - also sets starttime to current
514void 
515CbcHeuristicFPump::setMaximumTime(double value)
516{
517  startTime_=CoinCpuTime();
518  maximumTime_=value;
519}
520
521 
Note: See TracBrowser for help on using the repository browser.