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

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

for osi methods

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