source: stable/2.0/Cbc/src/CbcHeuristicFPump.cpp @ 905

Last change on this file since 905 was 905, checked in by ladanyi, 13 years ago

include cstdlib before cmath to get things to compile on AIX with xlC (same as changeset 904 in trunk)

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