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

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

update branches/devel for threads

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