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

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

probably breaking everything

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