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

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

to make easier to use as function

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.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 "CoinWarmStartBasis.hpp"
18#include "CoinTime.hpp"
19
20
21// Default Constructor
22CbcHeuristicFPump::CbcHeuristicFPump() 
23  :CbcHeuristic(),
24   startTime_(0.0),
25   maximumTime_(0.0),
26   fakeCutoff_(COIN_DBL_MAX),
27   absoluteIncrement_(0.0),
28   relativeIncrement_(0.0),
29   defaultRounding_(0.49999),
30   initialWeight_(0.0),
31   weightFactor_(0.1),
32   maximumPasses_(100),
33   maximumRetries_(1),
34   accumulate_(0),
35   roundExpensive_(false)
36{
37  setWhen(1);
38}
39
40// Constructor from model
41CbcHeuristicFPump::CbcHeuristicFPump(CbcModel & model,
42                                     double downValue,bool roundExpensive)
43  :CbcHeuristic(model),
44   startTime_(0.0),
45   maximumTime_(0.0),
46   fakeCutoff_(COIN_DBL_MAX),
47   absoluteIncrement_(0.0),
48   relativeIncrement_(0.0),
49   defaultRounding_(downValue),
50   initialWeight_(0.0),
51   weightFactor_(0.1),
52   maximumPasses_(100),
53   maximumRetries_(1),
54   accumulate_(0),
55   roundExpensive_(roundExpensive)
56{
57  setWhen(1);
58}
59
60// Destructor
61CbcHeuristicFPump::~CbcHeuristicFPump ()
62{
63}
64
65// Clone
66CbcHeuristic *
67CbcHeuristicFPump::clone() const
68{
69  return new CbcHeuristicFPump(*this);
70}
71// Create C++ lines to get to current state
72void 
73CbcHeuristicFPump::generateCpp( FILE * fp) 
74{
75  CbcHeuristicFPump other;
76  fprintf(fp,"0#include \"CbcHeuristicFPump.hpp\"\n");
77  fprintf(fp,"3  CbcHeuristicFPump heuristicFPump(*cbcModel);\n");
78  CbcHeuristic::generateCpp(fp,"heuristicFPump");
79  if (maximumPasses_!=other.maximumPasses_)
80    fprintf(fp,"3  heuristicFPump.setMaximumPasses(%d);\n",maximumPasses_);
81  else
82    fprintf(fp,"4  heuristicFPump.setMaximumPasses(%d);\n",maximumPasses_);
83  if (maximumRetries_!=other.maximumRetries_)
84    fprintf(fp,"3  heuristicFPump.setMaximumRetries(%d);\n",maximumRetries_);
85  else
86    fprintf(fp,"4  heuristicFPump.setMaximumRetries(%d);\n",maximumRetries_);
87  if (accumulate_!=other.accumulate_)
88    fprintf(fp,"3  heuristicFPump.setAccumulate(%d);\n",accumulate_);
89  else
90    fprintf(fp,"4  heuristicFPump.setAccumulate(%d);\n",accumulate_);
91  if (maximumTime_!=other.maximumTime_)
92    fprintf(fp,"3  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
93  else
94    fprintf(fp,"4  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
95  if (fakeCutoff_!=other.fakeCutoff_)
96    fprintf(fp,"3  heuristicFPump.setFakeCutoff(%g);\n",fakeCutoff_);
97  else
98    fprintf(fp,"4  heuristicFPump.setFakeCutoff(%g);\n",fakeCutoff_);
99  if (absoluteIncrement_!=other.absoluteIncrement_)
100    fprintf(fp,"3  heuristicFPump.setAbsoluteIncrement(%g);\n",absoluteIncrement_);
101  else
102    fprintf(fp,"4  heuristicFPump.setAbsoluteIncrement(%g);\n",absoluteIncrement_);
103  if (relativeIncrement_!=other.relativeIncrement_)
104    fprintf(fp,"3  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
105  else
106    fprintf(fp,"4  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
107  if (defaultRounding_!=other.defaultRounding_)
108    fprintf(fp,"3  heuristicFPump.setDefaultRounding(%g);\n",defaultRounding_);
109  else
110    fprintf(fp,"4  heuristicFPump.setDefaultRounding(%g);\n",defaultRounding_);
111  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicFPump);\n");
112  if (initialWeight_!=other.initialWeight_)
113    fprintf(fp,"3  heuristicFPump.setInitialWeight(%g);\n",initialWeight_);
114  else
115    fprintf(fp,"4  heuristicFPump.setInitialWeight(%g);\n",initialWeight_);
116  if (weightFactor_!=other.weightFactor_)
117    fprintf(fp,"3  heuristicFPump.setWeightFactor(%g);\n",weightFactor_);
118  else
119    fprintf(fp,"4  heuristicFPump.setWeightFactor(%g);\n",weightFactor_);
120}
121
122// Copy constructor
123CbcHeuristicFPump::CbcHeuristicFPump(const CbcHeuristicFPump & rhs)
124:
125  CbcHeuristic(rhs),
126  startTime_(rhs.startTime_),
127  maximumTime_(rhs.maximumTime_),
128  fakeCutoff_(rhs.fakeCutoff_),
129  absoluteIncrement_(rhs.absoluteIncrement_),
130  relativeIncrement_(rhs.relativeIncrement_),
131  defaultRounding_(rhs.defaultRounding_),
132  initialWeight_(rhs.initialWeight_),
133  weightFactor_(rhs.weightFactor_),
134  maximumPasses_(rhs.maximumPasses_),
135  maximumRetries_(rhs.maximumRetries_),
136  accumulate_(rhs.accumulate_),
137  roundExpensive_(rhs.roundExpensive_)
138{
139}
140
141// Assignment operator
142CbcHeuristicFPump & 
143CbcHeuristicFPump::operator=( const CbcHeuristicFPump& rhs)
144{
145  if (this!=&rhs) {
146    CbcHeuristic::operator=(rhs);
147    startTime_ = rhs.startTime_;
148    maximumTime_ = rhs.maximumTime_;
149    fakeCutoff_ = rhs.fakeCutoff_;
150    absoluteIncrement_ = rhs.absoluteIncrement_;
151    relativeIncrement_ = rhs.relativeIncrement_;
152    defaultRounding_ = rhs.defaultRounding_;
153    initialWeight_ = rhs.initialWeight_;
154    weightFactor_ = rhs.weightFactor_;
155    maximumPasses_ = rhs.maximumPasses_;
156    maximumRetries_ = rhs.maximumRetries_;
157    accumulate_ = rhs.accumulate_;
158    roundExpensive_ = rhs.roundExpensive_;
159  }
160  return *this;
161}
162
163// Resets stuff if model changes
164void 
165CbcHeuristicFPump::resetModel(CbcModel * model)
166{
167}
168
169/**************************BEGIN MAIN PROCEDURE ***********************************/
170
171// See if feasibility pump will give better solution
172// Sets value of solution
173// Returns 1 if solution, 0 if not
174int
175CbcHeuristicFPump::solution(double & solutionValue,
176                         double * betterSolution)
177{
178#define LEN_PRINT 200
179  char pumpPrint[LEN_PRINT];
180  pumpPrint[0]='\0';
181  if (!when()||(when()==1&&model_->phase()!=1))
182    return 0; // switched off
183  // See if at root node
184  bool atRoot = model_->getNodeCount()==0;
185  int passNumber = model_->getCurrentPassNumber();
186  // just do once
187  if (!atRoot||passNumber!=1)
188    return 0;
189  // probably a good idea
190  if (model_->getSolutionCount()) return 0;
191  // loop round doing repeated pumps
192  double cutoff;
193  model_->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
194  double direction = model_->solver()->getObjSense();
195  cutoff *= direction;
196  double firstCutoff = fabs(cutoff);
197  cutoff = CoinMin(cutoff,solutionValue);
198  // check plausible and space for rounded solution
199  int numberColumns = model_->getNumCols();
200  int numberIntegers = model_->numberIntegers();
201  const int * integerVariableOrig = model_->integerVariable();
202
203  // 1. initially check 0-1
204  int i,j;
205  int general=0;
206  int * integerVariable = new int[numberIntegers];
207  const double * lower = model_->solver()->getColLower();
208  const double * upper = model_->solver()->getColUpper();
209  bool doGeneral = (accumulate_&(4+8))!=0;
210  j=0;
211  for (i=0;i<numberIntegers;i++) {
212    int iColumn = integerVariableOrig[i];
213#ifndef NDEBUG
214    const OsiObject * object = model_->object(i);
215    const CbcSimpleInteger * integerObject = 
216      dynamic_cast<const  CbcSimpleInteger *> (object);
217    const OsiSimpleInteger * integerObject2 = 
218      dynamic_cast<const  OsiSimpleInteger *> (object);
219    assert(integerObject||integerObject2);
220#endif
221    if (upper[iColumn]-lower[iColumn]>1.000001) {
222      general++;
223      if (doGeneral)
224        integerVariable[j++]=iColumn;
225    } else {
226      integerVariable[j++]=iColumn;
227    }
228  }
229  if (general*3>2*numberIntegers&&!doGeneral) {
230    delete [] integerVariable;
231    return 0;
232  } else if ((accumulate_&8)==0) {
233    doGeneral=false;
234  }
235  if (!general)
236    doGeneral=false;
237  // For solution closest to feasible if none found
238  int * closestSolution = general ? NULL : new int[numberIntegers];
239  double closestObjectiveValue = COIN_DBL_MAX;
240 
241  int numberIntegersOrig = numberIntegers;
242  numberIntegers = j;
243  double * newSolution = new double [numberColumns];
244  double newSolutionValue=COIN_DBL_MAX;
245  bool solutionFound=false;
246  char * usedColumn = NULL;
247  double * lastSolution=NULL;
248  int fixContinuous=0;
249  bool fixInternal=false;
250  bool allSlack=false;
251  if (when_>=21&&when_<=25) {
252    when_ -= 10;
253    allSlack=true;
254  }
255  double time1 = CoinCpuTime();
256  if (when_>=11&&when_<=15) {
257    fixInternal = when_ >11&&when_<15;
258    if (when_<13)
259      fixContinuous = 0;
260    else if (when_!=14)
261      fixContinuous=1;
262    else
263      fixContinuous=2;
264    when_=1;
265    usedColumn = new char [numberColumns];
266    memset(usedColumn,0,numberColumns);
267    model_->solver()->resolve();
268    assert (model_->solver()->isProvenOptimal());
269    lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(),numberColumns);
270  }
271  int finalReturnCode=0;
272  int totalNumberPasses=0;
273  int numberTries=0;
274  CoinWarmStartBasis bestBasis;
275  while (true) {
276    int numberPasses=0;
277    numberTries++;
278    // Clone solver - otherwise annoys root node computations
279    OsiSolverInterface * solver = model_->solver()->clone();
280    // if cutoff exists then add constraint
281    bool useCutoff = (fabs(cutoff)<1.0e20&&(fakeCutoff_!=COIN_DBL_MAX||numberTries>1));
282    // but there may be a close one
283    if (firstCutoff<2.0*solutionValue&&numberTries==1) 
284      useCutoff=true;
285    if (useCutoff) {
286      cutoff = CoinMin(cutoff,fakeCutoff_);
287      const double * objective = solver->getObjCoefficients();
288      int numberColumns = solver->getNumCols();
289      int * which = new int[numberColumns];
290      double * els = new double[numberColumns];
291      int nel=0;
292      for (int i=0;i<numberColumns;i++) {
293        double value = objective[i];
294        if (value) {
295          which[nel]=i;
296          els[nel++] = direction*value;
297        }
298      }
299      solver->addRow(nel,which,els,-COIN_DBL_MAX,cutoff);
300      delete [] which;
301      delete [] els;
302      bool takeHint;
303      OsiHintStrength strength;
304      solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
305      solver->setHintParam(OsiDoDualInResolve,true,OsiHintDo);
306      solver->resolve();
307      solver->setHintParam(OsiDoDualInResolve,takeHint,strength);
308    }
309    solver->setDblParam(OsiDualObjectiveLimit,1.0e50);
310    solver->resolve();
311    // Solver may not be feasible
312    if (!solver->isProvenOptimal()) {
313      break;
314    }
315    const double * lower = solver->getColLower();
316    const double * upper = solver->getColUpper();
317    const double * solution = solver->getColSolution();
318    if (lastSolution)
319      memcpy(lastSolution,solution,numberColumns*sizeof(double));
320    double primalTolerance;
321    solver->getDblParam(OsiPrimalTolerance,primalTolerance);
322   
323   
324    // 2 space for last rounded solutions
325#define NUMBER_OLD 4
326    double ** oldSolution = new double * [NUMBER_OLD];
327    for (j=0;j<NUMBER_OLD;j++) {
328      oldSolution[j]= new double[numberColumns];
329      for (i=0;i<numberColumns;i++) oldSolution[j][i]=-COIN_DBL_MAX;
330    }
331   
332    // 3. Replace objective with an initial 0-valued objective
333    double * saveObjective = new double [numberColumns];
334    memcpy(saveObjective,solver->getObjCoefficients(),numberColumns*sizeof(double));
335    for (i=0;i<numberColumns;i++) {
336      solver->setObjCoeff(i,0.0);
337    }
338    bool finished=false;
339    double direction = solver->getObjSense();
340    int returnCode=0;
341    bool takeHint;
342    OsiHintStrength strength;
343    solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
344    solver->setHintParam(OsiDoDualInResolve,false);
345    //solver->messageHandler()->setLogLevel(0);
346   
347    // 4. Save objective offset so we can see progress
348    double saveOffset;
349    solver->getDblParam(OsiObjOffset,saveOffset);
350    // Get amount for original objective
351    double scaleFactor = 0.0;
352    for (i=0;i<numberColumns;i++)
353      scaleFactor += saveObjective[i]*saveObjective[i];
354    if (scaleFactor)
355      scaleFactor = (initialWeight_*sqrt((double) numberIntegers))/sqrt(scaleFactor);
356    // 5. MAIN WHILE LOOP
357    bool newLineNeeded=false;
358    while (!finished) {
359      returnCode=0;
360      // see what changed
361      if (usedColumn) {
362        for (i=0;i<numberColumns;i++) {
363          if (fabs(solution[i]-lastSolution[i])>1.0e-8) 
364            usedColumn[i]=1;
365          lastSolution[i]=solution[i];
366        }
367      }
368      if (numberPasses>=maximumPasses_) {
369        break;
370      }
371      if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
372      memcpy(newSolution,solution,numberColumns*sizeof(double));
373      int flip;
374      returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
375                          pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
376      numberPasses++;
377      if (returnCode) {
378        // SOLUTION IS INTEGER
379        // Put back correct objective
380        for (i=0;i<numberColumns;i++)
381          solver->setObjCoeff(i,saveObjective[i]);
382
383        // solution - but may not be better
384        // Compute using dot product
385        solver->setDblParam(OsiObjOffset,saveOffset);
386        newSolutionValue = -saveOffset;
387        for (  i=0 ; i<numberColumns ; i++ )
388          newSolutionValue += saveObjective[i]*newSolution[i];
389        newSolutionValue *= direction;
390        sprintf(pumpPrint+strlen(pumpPrint)," - solution found of %g",newSolutionValue);
391        newLineNeeded=false;
392        if (newSolutionValue<solutionValue) {
393          double saveValue = newSolutionValue;
394          if (!doGeneral) {
395            int numberLeft=0;
396            for (i=0;i<numberIntegersOrig;i++) {
397              int iColumn = integerVariableOrig[i];
398              double value = floor(newSolution[iColumn]+0.5);
399              if(solver->isBinary(iColumn)) {
400                solver->setColLower(iColumn,value);
401                solver->setColUpper(iColumn,value);
402              } else {
403                if (fabs(value-newSolution[iColumn])>1.0e-7) 
404                  numberLeft++;
405              }
406            }
407            if (numberLeft) {
408              returnCode = smallBranchAndBound(solver,numberNodes_,newSolution,newSolutionValue,
409                                               solutionValue,"CbcHeuristicFpump");
410              if ((returnCode&2)!=0) {
411                // could add cut
412                returnCode &= ~2;
413              }
414            }
415          }
416          if (returnCode&&newSolutionValue<saveValue) {
417            memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
418            CoinWarmStartBasis * basis =
419              dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
420            if (basis) {
421              bestBasis = * basis;
422              delete basis;
423            }
424            if ((accumulate_&1)!=0)
425              model_->incrementUsed(betterSolution); // for local search
426            solutionValue=newSolutionValue;
427            solutionFound=true;
428            if (general&&saveValue!=newSolutionValue)
429              sprintf(pumpPrint+strlen(pumpPrint)," - cleaned solution of %g",solutionValue);
430            if (pumpPrint[0]!='\0')
431              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
432                << pumpPrint
433                <<CoinMessageEol;
434            pumpPrint[0]='\0';
435          } else {
436            sprintf(pumpPrint+strlen(pumpPrint)," - not improved by mini branch and bound");
437            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
438              << pumpPrint
439              <<CoinMessageEol;
440            pumpPrint[0]='\0';
441          }
442        } else {
443          sprintf(pumpPrint+strlen(pumpPrint)," - not good enough");
444          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
445            << pumpPrint
446            <<CoinMessageEol;
447          pumpPrint[0]='\0';
448          newLineNeeded=false;
449          returnCode=0;
450        }     
451        break;
452      } else {
453        // SOLUTION IS not INTEGER
454        // 1. check for loop
455        bool matched;
456        for (int k = NUMBER_OLD-1; k > 0; k--) {
457          double * b = oldSolution[k];
458          matched = true;
459          for (i = 0; i <numberIntegers; i++) {
460            int iColumn = integerVariable[i];
461            if (newSolution[iColumn]!=b[iColumn]) {
462              matched=false;
463              break;
464            }
465          }
466          if (matched) break;
467        }
468        if (matched || numberPasses%100 == 0) {
469          // perturbation
470          sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
471          newLineNeeded=true;
472          for (i=0;i<numberIntegers;i++) {
473            int iColumn = integerVariable[i];
474            double value = max(0.0,CoinDrand48()-0.3);
475            double difference = fabs(solution[iColumn]-newSolution[iColumn]);
476            if (difference+value>0.5) {
477              if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
478                newSolution[iColumn] += 1.0;
479              } else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
480                newSolution[iColumn] -= 1.0;
481              } else {
482                // general integer
483                if (difference+value>0.75)
484                  newSolution[iColumn] += 1.0;
485                else
486                  newSolution[iColumn] -= 1.0;
487              }
488            }
489          }
490        } else {
491          for (j=NUMBER_OLD-1;j>0;j--) {
492            for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
493          }
494          for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
495        }
496       
497        // 2. update the objective function based on the new rounded solution
498        double offset=0.0;
499        double costValue = (1.0-scaleFactor)*solver->getObjSense();
500       
501        for (i=0;i<numberColumns;i++) {
502          // below so we can keep original code and allow for objective
503          int iColumn = i;
504          if(!solver->isBinary(iColumn)&&!doGeneral)
505            continue;
506          // deal with fixed variables (i.e., upper=lower)
507          if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) {
508            //if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
509            //else                                       solver->setObjCoeff(iColumn,costValue);
510            continue;
511          }
512          if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
513            solver->setObjCoeff(iColumn,costValue+scaleFactor*saveObjective[iColumn]);
514          } else {
515            if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
516              solver->setObjCoeff(iColumn,-costValue+scaleFactor*saveObjective[iColumn]);
517            } else {
518              solver->setObjCoeff(iColumn,0.0);
519            }
520          }
521          offset += costValue*newSolution[iColumn];
522        }
523        solver->setDblParam(OsiObjOffset,-offset);
524        if (!general&&false) {
525          // Solve in two goes - first keep satisfied ones fixed
526          double * saveLower = new double [numberIntegers];
527          double * saveUpper = new double [numberIntegers];
528          for (i=0;i<numberIntegers;i++) {
529            int iColumn = integerVariable[i];
530            saveLower[i]=COIN_DBL_MAX;
531            saveUpper[i]=-COIN_DBL_MAX;
532            if (solution[iColumn]<lower[iColumn]+primalTolerance) {
533              saveUpper[i]=upper[iColumn];
534              solver->setColUpper(iColumn,lower[iColumn]);
535            } else if (solution[iColumn]>upper[iColumn]-primalTolerance) {
536              saveLower[i]=lower[iColumn];
537              solver->setColLower(iColumn,upper[iColumn]);
538            }
539          }
540          solver->resolve();
541          assert (solver->isProvenOptimal());
542          for (i=0;i<numberIntegers;i++) {
543            int iColumn = integerVariable[i];
544            if (saveLower[i]!=COIN_DBL_MAX)
545              solver->setColLower(iColumn,saveLower[i]);
546            if (saveUpper[i]!=-COIN_DBL_MAX)
547              solver->setColUpper(iColumn,saveUpper[i]);
548            saveUpper[i]=-COIN_DBL_MAX;
549          }
550          memcpy(newSolution,solution,numberColumns*sizeof(double));
551          int flip;
552          returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
553                              pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
554          numberPasses++;
555          if (returnCode) {
556            // solution - but may not be better
557            // Compute using dot product
558            double newSolutionValue = -saveOffset;
559            for (  i=0 ; i<numberColumns ; i++ )
560              newSolutionValue += saveObjective[i]*newSolution[i];
561            newSolutionValue *= direction;
562            sprintf(pumpPrint+strlen(pumpPrint)," - intermediate solution found of %g",newSolutionValue);
563            if (newSolutionValue<solutionValue) {
564              memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
565              CoinWarmStartBasis * basis =
566                dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
567              if (basis) {
568                bestBasis = * basis;
569                delete basis;
570              }
571              if ((accumulate_&1)!=0)
572                model_->incrementUsed(betterSolution); // for local search
573              solutionValue=newSolutionValue;
574              solutionFound=true;
575            } else {
576              returnCode=0;
577            }
578          }     
579        }
580        if (!doGeneral) {
581          // faster to do from all slack!!!!
582          if (allSlack) {
583            CoinWarmStartBasis dummy;
584            solver->setWarmStart(&dummy);
585          }
586          solver->resolve();
587          assert (solver->isProvenOptimal());
588          // in case very dubious solver
589          lower = solver->getColLower();
590          upper = solver->getColUpper();
591          solution = solver->getColSolution();
592        } else {
593          int * addStart = new int[2*general+1];
594          int * addIndex = new int[4*general];
595          double * addElement = new double[4*general];
596          double * addLower = new double[2*general];
597          double * addUpper = new double[2*general];
598          double * obj = new double[general];
599          int nAdd=0;
600          for (i=0;i<numberIntegers;i++) {
601            int iColumn = integerVariable[i];
602            if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
603                newSolution[iColumn]<upper[iColumn]-primalTolerance) {
604              obj[nAdd]=1.0;
605              addLower[nAdd]=0.0;
606              addUpper[nAdd]=COIN_DBL_MAX;
607              nAdd++;
608            }
609          }
610          OsiSolverInterface * solver2 = solver;
611          if (nAdd) {
612            CoinZeroN(addStart,nAdd+1);
613            solver2 = solver->clone();
614            solver2->addCols(nAdd,addStart,NULL,NULL,addLower,addUpper,obj);
615            // feasible solution
616            double * sol = new double[nAdd+numberColumns];
617            memcpy(sol,solution,numberColumns*sizeof(double));
618            // now rows
619            int nAdd=0;
620            int nEl=0;
621            int nAddRow=0;
622            for (i=0;i<numberIntegers;i++) {
623              int iColumn = integerVariable[i];
624              if (newSolution[iColumn]>lower[iColumn]+primalTolerance&&
625                  newSolution[iColumn]<upper[iColumn]-primalTolerance) {
626                addLower[nAddRow]=-newSolution[iColumn];;
627                addUpper[nAddRow]=COIN_DBL_MAX;
628                addIndex[nEl] = iColumn;
629                addElement[nEl++]=-1.0;
630                addIndex[nEl] = numberColumns+nAdd;
631                addElement[nEl++]=1.0;
632                nAddRow++;
633                addStart[nAddRow]=nEl;
634                addLower[nAddRow]=newSolution[iColumn];;
635                addUpper[nAddRow]=COIN_DBL_MAX;
636                addIndex[nEl] = iColumn;
637                addElement[nEl++]=1.0;
638                addIndex[nEl] = numberColumns+nAdd;
639                addElement[nEl++]=1.0;
640                nAddRow++;
641                addStart[nAddRow]=nEl;
642                sol[nAdd+numberColumns] = fabs(sol[iColumn]-newSolution[iColumn]);
643                nAdd++;
644              }
645            }
646            solver2->setColSolution(sol);
647            delete [] sol;
648            solver2->addRows(nAddRow,addStart,addIndex,addElement,addLower,addUpper);
649          }
650          delete [] addStart;
651          delete [] addIndex;
652          delete [] addElement;
653          delete [] addLower;
654          delete [] addUpper;
655          delete [] obj;
656          solver2->resolve();
657          assert (solver2->isProvenOptimal());
658          if (nAdd) {
659            solver->setColSolution(solver2->getColSolution());
660            delete solver2;
661          }
662        }
663        if (pumpPrint[0]!='\0')
664          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
665            << pumpPrint
666            <<CoinMessageEol;
667        pumpPrint[0]='\0';
668        if (solver->getNumRows()<3000)
669          sprintf(pumpPrint+strlen(pumpPrint),"Pass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue());
670        else
671          sprintf(pumpPrint+strlen(pumpPrint),"Pass %3d: (%.2f seconds) obj. %10.5f --> ", numberPasses+totalNumberPasses,
672                  model_->getCurrentSeconds(),solver->getObjValue());
673        if (closestSolution&&solver->getObjValue()<closestObjectiveValue) {
674          int i;
675          const double * objective = solver->getObjCoefficients();
676          for (i=0;i<numberIntegersOrig;i++) {
677            int iColumn=integerVariableOrig[i];
678            if (objective[iColumn]>0.0)
679              closestSolution[i]=0;
680            else
681              closestSolution[i]=1;
682          }
683          closestObjectiveValue = solver->getObjValue();
684        }
685        newLineNeeded=true;
686       
687      }
688      // reduce scale factor
689      scaleFactor *= weightFactor_;
690    } // END WHILE
691    if (!solutionFound) {
692      sprintf(pumpPrint+strlen(pumpPrint),"No solution found this major pass");
693      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
694        << pumpPrint
695        <<CoinMessageEol;
696      pumpPrint[0]='\0';
697    }
698    delete solver;
699    for ( j=0;j<NUMBER_OLD;j++) 
700      delete [] oldSolution[j];
701    delete [] oldSolution;
702    delete [] saveObjective;
703    if (usedColumn) {
704      OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
705      const double * colLower = newSolver->getColLower();
706      const double * colUpper = newSolver->getColUpper();
707      int i;
708      int nFix=0;
709      int nFixI=0;
710      int nFixC=0;
711      int nFixC2=0;
712      for (i=0;i<numberIntegersOrig;i++) {
713        int iColumn=integerVariableOrig[i];
714        //const OsiObject * object = model_->object(i);
715        //double originalLower;
716        //double originalUpper;
717        //getIntegerInformation( object,originalLower, originalUpper);
718        //assert(colLower[iColumn]==originalLower);
719        //newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
720        newSolver->setColLower(iColumn,colLower[iColumn]);
721        //assert(colUpper[iColumn]==originalUpper);
722        //newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
723        newSolver->setColUpper(iColumn,colUpper[iColumn]);
724        if (!usedColumn[iColumn]) {
725          double value=lastSolution[iColumn];
726          double nearest = floor(value+0.5);
727          if (fabs(value-nearest)<1.0e-7) {
728            if (nearest==colLower[iColumn]) {
729              newSolver->setColUpper(iColumn,colLower[iColumn]);
730              nFix++;
731            } else if (nearest==colUpper[iColumn]) {
732              newSolver->setColLower(iColumn,colUpper[iColumn]);
733              nFix++;
734            } else if (fixInternal) {
735              newSolver->setColLower(iColumn,nearest);
736              newSolver->setColUpper(iColumn,nearest);
737              nFix++;
738              nFixI++;
739            }
740          }
741        }
742      }
743      if (fixContinuous) {
744        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
745          if (!newSolver->isInteger(iColumn)&&!usedColumn[iColumn]) {
746            double value=lastSolution[iColumn];
747            if (value<colLower[iColumn]+1.0e-8) {
748              newSolver->setColUpper(iColumn,colLower[iColumn]);
749              nFixC++;
750            } else if (value>colUpper[iColumn]-1.0e-8) {
751              newSolver->setColLower(iColumn,colUpper[iColumn]);
752              nFixC++;
753            } else if (fixContinuous==2) {
754              newSolver->setColLower(iColumn,value);
755              newSolver->setColUpper(iColumn,value);
756              nFixC++;
757              nFixC2++;
758            }
759          }
760        }
761      }
762      newSolver->initialSolve();
763      if (!newSolver->isProvenOptimal()) {
764        //newSolver->writeMps("bad.mps");
765        assert (newSolver->isProvenOptimal());
766        break;
767      }
768      sprintf(pumpPrint+strlen(pumpPrint),"Before mini branch and bound, %d integers at bound fixed and %d continuous",
769             nFix,nFixC);
770      if (nFixC2+nFixI!=0)
771        sprintf(pumpPrint+strlen(pumpPrint)," of which %d were internal integer and %d internal continuous",
772                nFixI,nFixC2);
773      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
774        << pumpPrint
775        <<CoinMessageEol;
776      pumpPrint[0]='\0';
777      double saveValue = newSolutionValue;
778      returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
779                                       cutoff,"CbcHeuristicLocalAfterFPump");
780      if ((returnCode&2)!=0) {
781        // could add cut
782        returnCode &= ~2;
783      }
784      if (returnCode&&newSolutionValue<saveValue) {
785        sprintf(pumpPrint+strlen(pumpPrint),"Mini branch and bound improved solution from %g to %g (%.2f seconds)",
786                saveValue,newSolutionValue,model_->getCurrentSeconds());
787        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
788          << pumpPrint
789          <<CoinMessageEol;
790        pumpPrint[0]='\0';
791        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
792        if (fixContinuous) {
793          // may be able to do even better
794          const double * lower = model_->solver()->getColLower();
795          const double * upper = model_->solver()->getColUpper();
796          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
797            if (newSolver->isInteger(iColumn)) {
798              double value=floor(newSolution[iColumn]+0.5);
799              newSolver->setColLower(iColumn,value);
800              newSolver->setColUpper(iColumn,value);
801            } else {
802              newSolver->setColLower(iColumn,lower[iColumn]);
803              newSolver->setColUpper(iColumn,upper[iColumn]);
804            }
805          }
806          newSolver->initialSolve();
807          if (newSolver->isProvenOptimal()) {
808            double value = newSolver->getObjValue()*newSolver->getObjSense();
809            if (value<newSolutionValue) {
810              sprintf(pumpPrint+strlen(pumpPrint),"Freeing continuous variables gives a solution of %g", value);
811              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
812                << pumpPrint
813                <<CoinMessageEol;
814              pumpPrint[0]='\0';
815              newSolutionValue=value;
816              memcpy(betterSolution,newSolver->getColSolution(),numberColumns*sizeof(double));
817            }
818          } else {
819            //newSolver->writeMps("bad3.mps");
820          }
821        } 
822        if ((accumulate_&1)!=0)
823          model_->incrementUsed(betterSolution); // for local search
824        solutionValue=newSolutionValue;
825        solutionFound=true;
826        CoinWarmStartBasis * basis =
827          dynamic_cast<CoinWarmStartBasis *>(newSolver->getWarmStart()) ;
828        if (basis) {
829          bestBasis = * basis;
830          delete basis;
831        }
832      } else {
833        sprintf(pumpPrint+strlen(pumpPrint),"Mini branch and bound did not improve solution (%.2f seconds)",
834                model_->getCurrentSeconds());
835        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
836          << pumpPrint
837          <<CoinMessageEol;
838        pumpPrint[0]='\0';
839      }
840      delete newSolver;
841    }
842    if (solutionFound) finalReturnCode=1;
843    cutoff = CoinMin(cutoff,solutionValue);
844    if (numberTries>=maximumRetries_||!solutionFound) {
845      break;
846    } else if (absoluteIncrement_>0.0||relativeIncrement_>0.0) {
847      solutionFound=false;
848      double gap = relativeIncrement_*fabs(solutionValue);
849      cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
850      sprintf(pumpPrint+strlen(pumpPrint),"Round again with cutoff of %g",cutoff);
851      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
852        << pumpPrint
853        <<CoinMessageEol;
854      pumpPrint[0]='\0';
855      if ((accumulate_&3)<2&&usedColumn)
856        memset(usedColumn,0,numberColumns);
857      totalNumberPasses += numberPasses-1;
858    } else {
859      break;
860    }
861  }
862  if (!finalReturnCode&&closestSolution&&closestObjectiveValue <= 10.0&&usedColumn) {
863    // try a bit of branch and bound
864    OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
865    const double * colLower = newSolver->getColLower();
866    const double * colUpper = newSolver->getColUpper();
867    int i;
868    double rhs = 0.0;
869    for (i=0;i<numberIntegersOrig;i++) {
870      int iColumn=integerVariableOrig[i];
871      int direction = closestSolution[i];
872      closestSolution[i]=iColumn;
873      if (direction==0) {
874        // keep close to LB
875        rhs += colLower[iColumn];
876        lastSolution[i]=1.0;
877      } else {
878        // keep close to UB
879        rhs -= colUpper[iColumn];
880        lastSolution[i]=-1.0;
881      }
882    }
883    newSolver->addRow(numberIntegersOrig,closestSolution,
884                      lastSolution,-COIN_DBL_MAX,rhs+10.0);
885    //double saveValue = newSolutionValue;
886    //newSolver->writeMps("sub");
887    int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
888                                     newSolutionValue,"CbcHeuristicLocalAfterFPump");
889    if ((returnCode&2)!=0) {
890      // could add cut
891      returnCode &= ~2;
892    }
893    if (returnCode) {
894      //printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
895      memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
896      //abort();
897      solutionValue=newSolutionValue;
898      solutionFound=true;
899    }
900    delete newSolver;
901  }
902  delete [] usedColumn;
903  delete [] lastSolution;
904  delete [] newSolution;
905  delete [] closestSolution;
906  delete [] integerVariable;
907  sprintf(pumpPrint,"After %.2f seconds - Feasibility pump exiting - took %.2f seconds",
908          model_->getCurrentSeconds(),CoinCpuTime()-time1);
909  model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
910    << pumpPrint
911    <<CoinMessageEol;
912  if (bestBasis.getNumStructural())
913    model_->setBestSolutionBasis(bestBasis);
914  return finalReturnCode;
915}
916
917/**************************END MAIN PROCEDURE ***********************************/
918
919// update model
920void CbcHeuristicFPump::setModel(CbcModel * model)
921{
922  model_ = model;
923}
924
925/* Rounds solution - down if < downValue
926   returns 1 if current is a feasible solution
927*/
928int 
929CbcHeuristicFPump::rounds(OsiSolverInterface * solver,double * solution,
930                          const double * objective,
931                          int numberIntegers, const int * integerVariable,
932                          char * pumpPrint, int & iter,
933                          bool roundExpensive, double downValue, int *flip)
934{
935  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
936  double primalTolerance ;
937  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
938
939  int i;
940
941  int numberColumns = model_->getNumCols();
942  // tmp contains the current obj coefficients
943  double * tmp = new double [numberColumns];
944  memcpy(tmp,solver->getObjCoefficients(),numberColumns*sizeof(double));
945  int flip_up = 0;
946  int flip_down  = 0;
947  double  v = CoinDrand48() * 20;
948  int nn = 10 + (int) v;
949  int nnv = 0;
950  int * list = new int [nn];
951  double * val = new double [nn];
952  for (i = 0; i < nn; i++) val[i] = .001;
953
954  // return rounded solution
955  for (i=0;i<numberIntegers;i++) {
956    int iColumn = integerVariable[i];
957    double value=solution[iColumn];
958    double round = floor(value+primalTolerance);
959    if (value-round > downValue) round += 1.;
960    if (round < integerTolerance && tmp[iColumn] < -1. + integerTolerance) flip_down++;
961    if (round > 1. - integerTolerance && tmp[iColumn] > 1. - integerTolerance) flip_up++;
962    if (flip_up + flip_down == 0) { 
963       for (int k = 0; k < nn; k++) {
964           if (fabs(value-round) > val[k]) {
965              nnv++;
966              for (int j = nn-2; j >= k; j--) {
967                  val[j+1] = val[j];
968                  list[j+1] = list[j];
969              } 
970              val[k] = fabs(value-round);
971              list[k] = iColumn;
972              break;
973           }
974       }
975    }
976    solution[iColumn] = round;
977  }
978
979  if (nnv > nn) nnv = nn;
980  if (iter != 0)
981    sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
982  *flip = flip_up + flip_down;
983  delete [] tmp;
984
985  const double * columnLower = solver->getColLower();
986  const double * columnUpper = solver->getColUpper();
987  if (*flip == 0 && iter != 0) {
988    sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
989     for (i = 0; i < nnv; i++) {
990       // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
991       int index = list[i];
992       double value = solution[index];
993       if (value<=1.0) {
994         solution[index] = 1.0-value;
995       } else if (value<columnLower[index]+integerTolerance) {
996         solution[index] = value+1.0;
997       } else if (value>columnUpper[index]-integerTolerance) {
998         solution[index] = value-1.0;
999       } else {
1000         solution[index] = value-1.0;
1001       }
1002     }
1003     *flip = nnv;
1004  } else {
1005    //sprintf(pumpPrint+strlen(pumpPrint)," ");
1006  }
1007  delete [] list; delete [] val;
1008  //iter++;
1009   
1010  const double * rowLower = solver->getRowLower();
1011  const double * rowUpper = solver->getRowUpper();
1012
1013  int numberRows = solver->getNumRows();
1014  // get row activities
1015  double * rowActivity = new double[numberRows];
1016  memset(rowActivity,0,numberRows*sizeof(double));
1017  solver->getMatrixByCol()->times(solution,rowActivity) ;
1018  double largestInfeasibility =0.0;
1019  for (i=0 ; i < numberRows ; i++) {
1020    largestInfeasibility = max(largestInfeasibility,
1021                               rowLower[i]-rowActivity[i]);
1022    largestInfeasibility = max(largestInfeasibility,
1023                               rowActivity[i]-rowUpper[i]);
1024  }
1025  delete [] rowActivity;
1026  return (largestInfeasibility>primalTolerance) ? 0 : 1;
1027}
1028// Set maximum Time (default off) - also sets starttime to current
1029void 
1030CbcHeuristicFPump::setMaximumTime(double value)
1031{
1032  startTime_=CoinCpuTime();
1033  maximumTime_=value;
1034}
1035
1036 
Note: See TracBrowser for help on using the repository browser.