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

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

fix bug I just put in fpump

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