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

Last change on this file since 625 was 625, checked in by forrest, 12 years ago

modify feasibility pump

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