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

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

fix tiny bug

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