source: branches/devel/Cbc/src/CbcHeuristicFPump.cpp @ 441

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

for local tree search and feasibility pump

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.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 "CoinTime.hpp"
18
19
20// Default Constructor
21CbcHeuristicFPump::CbcHeuristicFPump() 
22  :CbcHeuristic(),
23   startTime_(0.0),
24   maximumTime_(0.0),
25   downValue_(0.5),
26   fakeCutoff_(COIN_DBL_MAX),
27   absoluteIncrement_(0.0),
28   relativeIncrement_(0.0),
29   maximumPasses_(100),
30   maximumRetries_(1),
31   accumulate_(0),
32   roundExpensive_(false)
33{
34  setWhen(1);
35}
36
37// Constructor from model
38CbcHeuristicFPump::CbcHeuristicFPump(CbcModel & model,
39                                     double downValue,bool roundExpensive)
40  :CbcHeuristic(model),
41   startTime_(0.0),
42   maximumTime_(0.0),
43   downValue_(downValue),
44   fakeCutoff_(COIN_DBL_MAX),
45   absoluteIncrement_(0.0),
46   relativeIncrement_(0.0),
47   maximumPasses_(100),
48   maximumRetries_(1),
49   accumulate_(0),
50   roundExpensive_(roundExpensive)
51{
52  setWhen(1);
53}
54
55// Destructor
56CbcHeuristicFPump::~CbcHeuristicFPump ()
57{
58}
59
60// Clone
61CbcHeuristic *
62CbcHeuristicFPump::clone() const
63{
64  return new CbcHeuristicFPump(*this);
65}
66// Create C++ lines to get to current state
67void 
68CbcHeuristicFPump::generateCpp( FILE * fp) 
69{
70  CbcHeuristicFPump other;
71  fprintf(fp,"0#include \"CbcHeuristicFPump.hpp\"\n");
72  fprintf(fp,"3  CbcHeuristicFPump heuristicFPump(*cbcModel);\n");
73  if (maximumPasses_!=other.maximumPasses_)
74    fprintf(fp,"3  heuristicFPump.setMaximumPasses(%d);\n",maximumPasses_);
75  else
76    fprintf(fp,"4  heuristicFPump.setMaximumPasses(%d);\n",maximumPasses_);
77  if (maximumRetries_!=other.maximumRetries_)
78    fprintf(fp,"3  heuristicFPump.setMaximumRetries(%d);\n",maximumRetries_);
79  else
80    fprintf(fp,"4  heuristicFPump.setMaximumRetries(%d);\n",maximumRetries_);
81  if (accumulate_!=other.accumulate_)
82    fprintf(fp,"3  heuristicFPump.setAccumulate(%d);\n",accumulate_);
83  else
84    fprintf(fp,"4  heuristicFPump.setAccumulate(%d);\n",accumulate_);
85  if (maximumTime_!=other.maximumTime_)
86    fprintf(fp,"3  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
87  else
88    fprintf(fp,"4  heuristicFPump.setMaximumTime(%g);\n",maximumTime_);
89  if (fakeCutoff_!=other.fakeCutoff_)
90    fprintf(fp,"3  heuristicFPump.setFakeCutoff(%g);\n",fakeCutoff_);
91  else
92    fprintf(fp,"4  heuristicFPump.setFakeCutoff(%g);\n",fakeCutoff_);
93  if (absoluteIncrement_!=other.absoluteIncrement_)
94    fprintf(fp,"3  heuristicFPump.setAbsoluteIncrement(%g);\n",absoluteIncrement_);
95  else
96    fprintf(fp,"4  heuristicFPump.setAbsoluteIncrement(%g);\n",absoluteIncrement_);
97  if (relativeIncrement_!=other.relativeIncrement_)
98    fprintf(fp,"3  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
99  else
100    fprintf(fp,"4  heuristicFPump.setRelativeIncrement(%g);\n",relativeIncrement_);
101  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicFPump);\n");
102}
103
104// Copy constructor
105CbcHeuristicFPump::CbcHeuristicFPump(const CbcHeuristicFPump & rhs)
106:
107  CbcHeuristic(rhs),
108  startTime_(rhs.startTime_),
109  maximumTime_(rhs.maximumTime_),
110  downValue_(rhs.downValue_),
111  fakeCutoff_(rhs.fakeCutoff_),
112  absoluteIncrement_(rhs.absoluteIncrement_),
113  relativeIncrement_(rhs.relativeIncrement_),
114  maximumPasses_(rhs.maximumPasses_),
115  maximumRetries_(rhs.maximumRetries_),
116  accumulate_(rhs.accumulate_),
117  roundExpensive_(rhs.roundExpensive_)
118{
119  setWhen(rhs.when());
120}
121// Resets stuff if model changes
122void 
123CbcHeuristicFPump::resetModel(CbcModel * model)
124{
125}
126
127/**************************BEGIN MAIN PROCEDURE ***********************************/
128
129// See if feasibility pump will give better solution
130// Sets value of solution
131// Returns 1 if solution, 0 if not
132int
133CbcHeuristicFPump::solution(double & solutionValue,
134                         double * betterSolution)
135{
136  if (!when()||(when()==1&&model_->phase()!=1))
137    return 0; // switched off
138  // See if at root node
139  bool atRoot = model_->getNodeCount()==0;
140  int passNumber = model_->getCurrentPassNumber();
141  // just do once
142  if (!atRoot||passNumber!=1)
143    return 0;
144  // probably a good idea
145  if (model_->getSolutionCount()) return 0;
146  // loop round doing repeated pumps
147  double cutoff;
148  model_->solver()->getDblParam(OsiDualObjectiveLimit,cutoff);
149  double direction = model_->solver()->getObjSense();
150  cutoff *= direction;
151  cutoff = CoinMin(cutoff,solutionValue);
152  // space for rounded solution
153  int numberColumns = model_->getNumCols();
154  double * newSolution = new double [numberColumns];
155  double newSolutionValue=COIN_DBL_MAX;
156  bool solutionFound=false;
157  char * usedColumn = NULL;
158  double * lastSolution=NULL;
159  int fixContinuous=0;
160  bool fixInternal=false;
161  if (when_>=11&&when_<=15) {
162    fixInternal = when_ >11&&when_<15;
163    if (when_<13)
164      fixContinuous = 0;
165    else if (when_!=14)
166      fixContinuous=1;
167    else
168      fixContinuous=2;
169    when_=1;
170    usedColumn = new char [numberColumns];
171    memset(usedColumn,0,numberColumns);
172    model_->solver()->resolve();
173    assert (model_->solver()->isProvenOptimal());
174    lastSolution = CoinCopyOfArray(model_->solver()->getColSolution(),numberColumns);
175  }
176  int finalReturnCode=0;
177  int totalNumberPasses=0;
178  int numberTries=0;
179  while (true) {
180    int numberPasses=0;
181    numberTries++;
182    // Clone solver - otherwise annoys root node computations
183    OsiSolverInterface * solver = model_->solver()->clone();
184    // if cutoff exists then add constraint
185    if (fabs(cutoff)<1.0e20&&(fakeCutoff_!=COIN_DBL_MAX||numberTries>1)) {
186      cutoff = CoinMin(cutoff,fakeCutoff_);
187      const double * objective = solver->getObjCoefficients();
188      int numberColumns = solver->getNumCols();
189      int * which = new int[numberColumns];
190      double * els = new double[numberColumns];
191      int nel=0;
192      for (int i=0;i<numberColumns;i++) {
193        double value = objective[i];
194        if (value) {
195          which[nel]=i;
196          els[nel++] = direction*value;
197        }
198      }
199      solver->addRow(nel,which,els,-COIN_DBL_MAX,cutoff);
200      delete [] which;
201      delete [] els;
202    }
203    solver->setDblParam(OsiDualObjectiveLimit,1.0e50);
204    solver->resolve();
205    // Solver may not be feasible
206    if (!solver->isProvenOptimal()) {
207      break;
208    }
209    const double * lower = solver->getColLower();
210    const double * upper = solver->getColUpper();
211    const double * solution = solver->getColSolution();
212    if (lastSolution)
213      memcpy(lastSolution,solution,numberColumns*sizeof(double));
214    double primalTolerance;
215    solver->getDblParam(OsiPrimalTolerance,primalTolerance);
216   
217    int numberIntegers = model_->numberIntegers();
218    const int * integerVariable = model_->integerVariable();
219
220    // 1. initially check 0-1
221    int i,j;
222    int general=0;
223    for (i=0;i<numberIntegers;i++) {
224      int iColumn = integerVariable[i];
225#ifndef NDEBUG
226      const OsiObject * object = model_->object(i);
227      const CbcSimpleInteger * integerObject = 
228        dynamic_cast<const  CbcSimpleInteger *> (object);
229      const OsiSimpleInteger * integerObject2 = 
230        dynamic_cast<const  OsiSimpleInteger *> (object);
231      assert(integerObject||integerObject2);
232#endif
233      if (upper[iColumn]-lower[iColumn]>1.000001) {
234        general++;
235        break;
236      }
237    }
238    if (general*3>numberIntegers) {
239      delete solver;
240      return 0;
241    }
242   
243    // 2 space for last rounded solutions
244#define NUMBER_OLD 4
245    double ** oldSolution = new double * [NUMBER_OLD];
246    for (j=0;j<NUMBER_OLD;j++) {
247      oldSolution[j]= new double[numberColumns];
248      for (i=0;i<numberColumns;i++) oldSolution[j][i]=-COIN_DBL_MAX;
249    }
250   
251    // 3. Replace objective with an initial 0-valued objective
252    double * saveObjective = new double [numberColumns];
253    memcpy(saveObjective,solver->getObjCoefficients(),numberColumns*sizeof(double));
254    for (i=0;i<numberColumns;i++) {
255      solver->setObjCoeff(i,0.0);
256    }
257    bool finished=false;
258    double direction = solver->getObjSense();
259    int returnCode=0;
260    bool takeHint;
261    OsiHintStrength strength;
262    solver->getHintParam(OsiDoDualInResolve,takeHint,strength);
263    solver->setHintParam(OsiDoDualInResolve,false);
264    solver->messageHandler()->setLogLevel(0);
265   
266    // 4. Save objective offset so we can see progress
267    double saveOffset;
268    solver->getDblParam(OsiObjOffset,saveOffset);
269   
270    // 5. MAIN WHILE LOOP
271    bool newLineNeeded=false;
272    while (!finished) {
273      returnCode=0;
274      // see what changed
275      if (usedColumn) {
276        for (i=0;i<numberColumns;i++) {
277          if (fabs(solution[i]-lastSolution[i])>1.0e-8) 
278            usedColumn[i]=1;
279          lastSolution[i]=solution[i];
280        }
281      }
282      if (numberPasses>=maximumPasses_) {
283        break;
284      }
285      if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
286      numberPasses++;
287      memcpy(newSolution,solution,numberColumns*sizeof(double));
288      int flip;
289      returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip);
290      if (returnCode) {
291        // SOLUTION IS INTEGER
292        // Put back correct objective
293        for (i=0;i<numberColumns;i++)
294          solver->setObjCoeff(i,saveObjective[i]);
295
296        // solution - but may not be better
297        // Compute using dot product
298        solver->setDblParam(OsiObjOffset,saveOffset);
299        newSolutionValue = -saveOffset;
300        for (  i=0 ; i<numberColumns ; i++ )
301          newSolutionValue += saveObjective[i]*newSolution[i];
302        newSolutionValue *= direction;
303        if (model_->logLevel())
304          printf(" - solution found of %g",newSolutionValue);
305        newLineNeeded=false;
306        if (newSolutionValue<solutionValue) {
307          double saveValue = newSolutionValue;
308          if (general) {
309            int numberLeft=0;
310            for (i=0;i<numberIntegers;i++) {
311              int iColumn = integerVariable[i];
312              double value = floor(newSolution[iColumn]+0.5);
313              if(solver->isBinary(iColumn)) {
314                solver->setColLower(iColumn,value);
315                solver->setColUpper(iColumn,value);
316              } else {
317                if (fabs(value-newSolution[iColumn])>1.0e-7) 
318                  numberLeft++;
319              }
320            }
321            if (numberLeft) {
322              returnCode = smallBranchAndBound(solver,200,newSolution,newSolutionValue,
323                                               solutionValue,"CbcHeuristicFpump");
324            }
325          }
326          if (returnCode) {
327            memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
328            if ((accumulate_&1)!=0)
329              model_->incrementUsed(betterSolution); // for local search
330            solutionValue=newSolutionValue;
331            solutionFound=true;
332            if (general&&saveValue!=newSolutionValue)
333              printf(" - cleaned solution of %g\n",solutionValue);
334            else
335              printf("\n");
336          } else {
337          if (model_->logLevel()) 
338            printf(" - not good enough after small branch and bound\n");
339          }
340        } else {
341          if (model_->logLevel()) 
342            printf(" - not good enough\n");
343          newLineNeeded=false;
344          returnCode=0;
345        }     
346        break;
347      } else {
348        // SOLUTION IS not INTEGER
349        // 1. check for loop
350        bool matched;
351        for (int k = NUMBER_OLD-1; k > 0; k--) {
352          double * b = oldSolution[k];
353          matched = true;
354          for (i = 0; i <numberIntegers; i++) {
355            int iColumn = integerVariable[i];
356            if(!solver->isBinary(iColumn))
357              continue;
358            if (newSolution[iColumn]!=b[iColumn]) {
359              matched=false;
360              break;
361            }
362          }
363          if (matched) break;
364        }
365        if (matched || numberPasses%100 == 0) {
366          // perturbation
367          if (model_->logLevel())
368            printf("Perturbation applied");
369          newLineNeeded=true;
370          for (i=0;i<numberIntegers;i++) {
371            int iColumn = integerVariable[i];
372            if(!solver->isBinary(iColumn))
373              continue;
374            double value = max(0.0,CoinDrand48()-0.3);
375            double difference = fabs(solution[iColumn]-newSolution[iColumn]);
376            if (difference+value>0.5) {
377              if (newSolution[iColumn]<lower[iColumn]+primalTolerance) newSolution[iColumn] += 1.0;
378              else if (newSolution[iColumn]>upper[iColumn]-primalTolerance) newSolution[iColumn] -= 1.0;
379              else abort();
380            }
381          }
382        } else {
383          for (j=NUMBER_OLD-1;j>0;j--) {
384            for (i = 0; i < numberColumns; i++) oldSolution[j][i]=oldSolution[j-1][i];
385          }
386          for (j = 0; j < numberColumns; j++) oldSolution[0][j] = newSolution[j];
387        }
388       
389        // 2. update the objective function based on the new rounded solution
390        double offset=0.0;
391        for (i=0;i<numberIntegers;i++) {
392          int iColumn = integerVariable[i];
393          if(!solver->isBinary(iColumn))
394            continue;
395          double costValue = 1.0;
396          // deal with fixed variables (i.e., upper=lower)
397          if (fabs(lower[iColumn]-upper[iColumn]) < primalTolerance) {
398            if (lower[iColumn] > 1. - primalTolerance) solver->setObjCoeff(iColumn,-costValue);
399            else                                       solver->setObjCoeff(iColumn,costValue);
400            continue;
401          }
402          if (newSolution[iColumn]<lower[iColumn]+primalTolerance) {
403            solver->setObjCoeff(iColumn,costValue);
404          } else {
405            if (newSolution[iColumn]>upper[iColumn]-primalTolerance) {
406              solver->setObjCoeff(iColumn,-costValue);
407            } else {
408              abort();
409            }
410          }
411          offset += costValue*newSolution[iColumn];
412        }
413        solver->setDblParam(OsiObjOffset,-offset);
414        if (!general&false) {
415          // Solve in two goes - first keep satisfied ones fixed
416          double * saveLower = new double [numberIntegers];
417          double * saveUpper = new double [numberIntegers];
418          for (i=0;i<numberIntegers;i++) {
419            int iColumn = integerVariable[i];
420            saveLower[i]=COIN_DBL_MAX;
421            saveUpper[i]=-COIN_DBL_MAX;
422            if (solution[iColumn]<lower[iColumn]+primalTolerance) {
423              saveUpper[i]=upper[iColumn];
424              solver->setColUpper(iColumn,lower[iColumn]);
425            } else if (solution[iColumn]>upper[iColumn]-primalTolerance) {
426              saveLower[i]=lower[iColumn];
427              solver->setColLower(iColumn,upper[iColumn]);
428            }
429          }
430          solver->resolve();
431          assert (solver->isProvenOptimal());
432          for (i=0;i<numberIntegers;i++) {
433            int iColumn = integerVariable[i];
434            if (saveLower[i]!=COIN_DBL_MAX)
435              solver->setColLower(iColumn,saveLower[i]);
436            if (saveUpper[i]!=-COIN_DBL_MAX)
437              solver->setColUpper(iColumn,saveUpper[i]);
438            saveUpper[i]=-COIN_DBL_MAX;
439          }
440          memcpy(newSolution,solution,numberColumns*sizeof(double));
441          int flip;
442          returnCode = rounds(newSolution,saveObjective,roundExpensive_,downValue_,&flip);
443          if (returnCode) {
444            // solution - but may not be better
445            // Compute using dot product
446            double newSolutionValue = -saveOffset;
447            for (  i=0 ; i<numberColumns ; i++ )
448              newSolutionValue += saveObjective[i]*newSolution[i];
449            newSolutionValue *= direction;
450            if (model_->logLevel())
451              printf(" - intermediate solution found of %g",newSolutionValue);
452            if (newSolutionValue<solutionValue) {
453              memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
454              if ((accumulate_&1)!=0)
455                model_->incrementUsed(betterSolution); // for local search
456              solutionValue=newSolutionValue;
457              solutionFound=true;
458            } else {
459              returnCode=0;
460            }
461          }     
462        }
463        solver->resolve();
464        assert (solver->isProvenOptimal());
465        if (model_->logLevel())
466          printf("\npass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue());
467        newLineNeeded=true;
468       
469      }
470    } // END WHILE
471    if (newLineNeeded&&model_->logLevel())
472      printf(" - no solution found\n");
473    delete solver;
474    for ( j=0;j<NUMBER_OLD;j++) 
475      delete [] oldSolution[j];
476    delete [] oldSolution;
477    delete [] saveObjective;
478    if (usedColumn) {
479      OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
480      const double * colLower = newSolver->getColLower();
481      const double * colUpper = newSolver->getColUpper();
482      int i;
483      int nFix=0;
484      int nFixI=0;
485      int nFixC=0;
486      int nFixC2=0;
487      for (i=0;i<numberIntegers;i++) {
488        int iColumn=integerVariable[i];
489        const OsiObject * object = model_->object(i);
490        double originalLower;
491        double originalUpper;
492        getIntegerInformation( object,originalLower, originalUpper); 
493        assert(colLower[iColumn]==originalLower);
494        newSolver->setColLower(iColumn,CoinMax(colLower[iColumn],originalLower));
495        assert(colUpper[iColumn]==originalUpper);
496        newSolver->setColUpper(iColumn,CoinMin(colUpper[iColumn],originalUpper));
497        if (!usedColumn[iColumn]) {
498          double value=lastSolution[iColumn];
499          double nearest = floor(value+0.5);
500          if (fabs(value-nearest)<1.0e-7) {
501            if (nearest==colLower[iColumn]) {
502              newSolver->setColUpper(iColumn,colLower[iColumn]);
503              nFix++;
504            } else if (nearest==colUpper[iColumn]) {
505              newSolver->setColLower(iColumn,colUpper[iColumn]);
506              nFix++;
507            } else if (fixInternal) {
508              newSolver->setColLower(iColumn,nearest);
509              newSolver->setColUpper(iColumn,nearest);
510              nFix++;
511              nFixI++;
512            }
513          }
514        }
515      }
516      if (fixContinuous) {
517        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
518          if (!newSolver->isInteger(iColumn)&&!usedColumn[iColumn]) {
519            double value=lastSolution[iColumn];
520            if (value<colLower[iColumn]+1.0e-8) {
521              newSolver->setColUpper(iColumn,colLower[iColumn]);
522              nFixC++;
523            } else if (value>colUpper[iColumn]-1.0e-8) {
524              newSolver->setColLower(iColumn,colUpper[iColumn]);
525              nFixC++;
526            } else if (fixContinuous==2) {
527              newSolver->setColLower(iColumn,value);
528              newSolver->setColUpper(iColumn,value);
529              nFixC++;
530              nFixC2++;
531            }
532          }
533        }
534      }
535      newSolver->initialSolve();
536      if (!newSolver->isProvenOptimal()) {
537        newSolver->writeMps("bad.mps");
538        assert (newSolver->isProvenOptimal());
539        break;
540      }
541      printf("%d integers at bound fixed and %d continuous",
542             nFix,nFixC);
543      if (nFixC2+nFixI==0)
544        printf("\n");
545      else
546        printf("of which %d were internal integer and %d internal continuous\n",
547             nFixI,nFixC2);
548      double saveValue = newSolutionValue;
549      returnCode = smallBranchAndBound(newSolver,200,newSolution,newSolutionValue,
550                                       newSolutionValue,"CbcHeuristicLocalAfterFPump");
551      if (returnCode) {
552        printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
553        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
554        if (fixContinuous) {
555          // may be able to do even better
556          const double * lower = model_->solver()->getColLower();
557          const double * upper = model_->solver()->getColUpper();
558          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
559            if (newSolver->isInteger(iColumn)) {
560              double value=floor(newSolution[iColumn]+0.5);
561              newSolver->setColLower(iColumn,value);
562              newSolver->setColUpper(iColumn,value);
563            } else {
564              newSolver->setColLower(iColumn,lower[iColumn]);
565              newSolver->setColUpper(iColumn,upper[iColumn]);
566            }
567          }
568          newSolver->initialSolve();
569          if (newSolver->isProvenOptimal()) {
570            double value = newSolver->getObjValue()*newSolver->getObjSense();
571            if (value<newSolutionValue) {
572              printf("freeing continuous gives a solution of %g\n", value);
573              newSolutionValue=value;
574              memcpy(betterSolution,newSolver->getColSolution(),numberColumns*sizeof(double));
575            }
576          } else {
577            newSolver->writeMps("bad3.mps");
578          }
579        } 
580        if ((accumulate_&1)!=0)
581          model_->incrementUsed(betterSolution); // for local search
582        solutionValue=newSolutionValue;
583        solutionFound=true;
584      }
585      delete newSolver;
586    }
587    if (solutionFound) finalReturnCode=1;
588    cutoff = CoinMin(cutoff,solutionValue);
589    if (numberTries>=maximumRetries_||!solutionFound) {
590      break;
591    } else if (absoluteIncrement_>0.0||relativeIncrement_>0.0) {
592      solutionFound=false;
593      double gap = relativeIncrement_*fabs(solutionValue);
594      cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
595      printf("round again with cutoff of %g\n",cutoff);
596      if (accumulate_<2)
597        memset(usedColumn,0,numberColumns);
598      totalNumberPasses += numberPasses;
599    } else {
600      break;
601    }
602  }
603  delete [] usedColumn;
604  delete [] lastSolution;
605  delete [] newSolution;
606  return finalReturnCode;
607}
608
609/**************************END MAIN PROCEDURE ***********************************/
610
611// update model
612void CbcHeuristicFPump::setModel(CbcModel * model)
613{
614  model_ = model;
615}
616
617/* Rounds solution - down if < downValue
618   returns 1 if current is a feasible solution
619*/
620int 
621CbcHeuristicFPump::rounds(double * solution,
622                          const double * objective,
623                          bool roundExpensive, double downValue, int *flip)
624{
625  OsiSolverInterface * solver = model_->solver();
626  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
627  double primalTolerance ;
628  solver->getDblParam(OsiPrimalTolerance,primalTolerance) ;
629  int numberIntegers = model_->numberIntegers();
630  const int * integerVariable = model_->integerVariable();
631
632  int i;
633
634  static int iter = 0;
635  int numberColumns = model_->getNumCols();
636  // tmp contains the current obj coefficients
637  double * tmp = new double [numberColumns];
638  memcpy(tmp,solver->getObjCoefficients(),numberColumns*sizeof(double));
639  int flip_up = 0;
640  int flip_down  = 0;
641  double  v = CoinDrand48() * 20;
642  int nn = 10 + (int) v;
643  int nnv = 0;
644  int * list = new int [nn];
645  double * val = new double [nn];
646  for (i = 0; i < nn; i++) val[i] = .001;
647
648  // return rounded solution
649  for (i=0;i<numberIntegers;i++) {
650    int iColumn = integerVariable[i];
651    if(!solver->isBinary(iColumn))
652      continue;
653#ifndef NDEBUG
654      const OsiObject * object = model_->object(i);
655      const CbcSimpleInteger * integerObject = 
656        dynamic_cast<const  CbcSimpleInteger *> (object);
657      const OsiSimpleInteger * integerObject2 = 
658        dynamic_cast<const  OsiSimpleInteger *> (object);
659      assert(integerObject||integerObject2);
660#endif
661    double value=solution[iColumn];
662    double round = floor(value+primalTolerance);
663    if (value-round > .5) round += 1.;
664    if (round < integerTolerance && tmp[iColumn] < -1. + integerTolerance) flip_down++;
665    if (round > 1. - integerTolerance && tmp[iColumn] > 1. - integerTolerance) flip_up++;
666    if (flip_up + flip_down == 0) { 
667       for (int k = 0; k < nn; k++) {
668           if (fabs(value-round) > val[k]) {
669              nnv++;
670              for (int j = nn-2; j >= k; j--) {
671                  val[j+1] = val[j];
672                  list[j+1] = list[j];
673              } 
674              val[k] = fabs(value-round);
675              list[k] = iColumn;
676              break;
677           }
678       }
679    }
680    solution[iColumn] = round;
681  }
682
683  if (nnv > nn) nnv = nn;
684  if (iter != 0&&model_->logLevel())
685    printf("up = %5d , down = %5d", flip_up, flip_down); fflush(stdout);
686  *flip = flip_up + flip_down;
687  delete [] tmp;
688
689  if (*flip == 0 && iter != 0) {
690    if(model_->logLevel())
691      printf(" -- rand = %4d (%4d) ", nnv, nn);
692     for (i = 0; i < nnv; i++) solution[list[i]] = 1. - solution[list[i]];
693     *flip = nnv;
694  } else if (model_->logLevel()) {
695    printf(" ");
696  }
697  delete [] list; delete [] val;
698  iter++;
699   
700  const double * rowLower = solver->getRowLower();
701  const double * rowUpper = solver->getRowUpper();
702
703  int numberRows = solver->getNumRows();
704  // get row activities
705  double * rowActivity = new double[numberRows];
706  memset(rowActivity,0,numberRows*sizeof(double));
707  solver->getMatrixByCol()->times(solution,rowActivity) ;
708  double largestInfeasibility =0.0;
709  for (i=0 ; i < numberRows ; i++) {
710    largestInfeasibility = max(largestInfeasibility,
711                               rowLower[i]-rowActivity[i]);
712    largestInfeasibility = max(largestInfeasibility,
713                               rowActivity[i]-rowUpper[i]);
714  }
715  delete [] rowActivity;
716  return (largestInfeasibility>primalTolerance) ? 0 : 1;
717}
718// Set maximum Time (default off) - also sets starttime to current
719void 
720CbcHeuristicFPump::setMaximumTime(double value)
721{
722  startTime_=CoinCpuTime();
723  maximumTime_=value;
724}
725
726 
Note: See TracBrowser for help on using the repository browser.