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

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

exit on cutoff a bit quicker

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