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

Last change on this file since 491 was 482, checked in by forrest, 14 years ago

for nonlinear and sprint

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