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

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

changes to heuristics and allow collection of more statistics

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