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

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

dj fixing and fpump changes

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