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

Last change on this file since 827 was 827, checked in by forrest, 12 years ago

fix bug if maximize

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