source: stable/2.4/Cbc/src/CbcHeuristicRINS.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

File size: 33.9 KB
Line 
1/* $Id: CbcHeuristicRINS.cpp 1261 2009-10-30 12:45:20Z forrest $ */
2// Copyright (C) 2006, International Business Machines
3// Corporation and others.  All Rights Reserved.
4#if defined(_MSC_VER)
5// Turn off compiler warning about long names
6#  pragma warning(disable:4786)
7#endif
8#include <cassert>
9#include <cstdlib>
10#include <cmath>
11#include <cfloat>
12
13#include "OsiSolverInterface.hpp"
14#include "CbcModel.hpp"
15#include "CbcMessage.hpp"
16#include "CbcHeuristicRINS.hpp"
17#include "CbcBranchActual.hpp"
18#include "CbcStrategy.hpp"
19#include "CglPreProcess.hpp"
20
21// Default Constructor
22CbcHeuristicRINS::CbcHeuristicRINS() 
23  :CbcHeuristic()
24{
25  numberSolutions_=0;
26  numberSuccesses_=0;
27  numberTries_=0;
28  stateOfFixing_ =0;
29  lastNode_=-999999;
30  howOften_=100;
31  decayFactor_ = 0.5; 
32  used_=NULL;
33  whereFrom_=1+8+255*256;
34}
35
36// Constructor with model - assumed before cuts
37
38CbcHeuristicRINS::CbcHeuristicRINS(CbcModel & model)
39  :CbcHeuristic(model)
40{
41  numberSolutions_=0;
42  numberSuccesses_=0;
43  numberTries_=0;
44  stateOfFixing_ =0;
45  lastNode_=-999999;
46  howOften_=100;
47  decayFactor_ = 0.5; 
48  assert(model.solver());
49  int numberColumns = model.solver()->getNumCols();
50  used_ = new char[numberColumns];
51  memset(used_,0,numberColumns);
52  whereFrom_=1+8+255*256;
53}
54
55// Destructor
56CbcHeuristicRINS::~CbcHeuristicRINS ()
57{
58  delete [] used_;
59}
60
61// Clone
62CbcHeuristic *
63CbcHeuristicRINS::clone() const
64{
65  return new CbcHeuristicRINS(*this);
66}
67
68// Assignment operator
69CbcHeuristicRINS & 
70CbcHeuristicRINS::operator=( const CbcHeuristicRINS& rhs)
71{
72  if (this!=&rhs) {
73    CbcHeuristic::operator=(rhs);
74    numberSolutions_ = rhs.numberSolutions_;
75    howOften_ = rhs.howOften_;
76    numberSuccesses_ = rhs.numberSuccesses_;
77    numberTries_ = rhs.numberTries_;
78    stateOfFixing_ = rhs.stateOfFixing_;
79    lastNode_=rhs.lastNode_;
80    delete [] used_;
81    if (model_&&rhs.used_) {
82      int numberColumns = model_->solver()->getNumCols();
83      used_ = new char[numberColumns];
84      memcpy(used_,rhs.used_,numberColumns);
85    } else {
86      used_=NULL;
87    }
88  }
89  return *this;
90}
91
92// Create C++ lines to get to current state
93void 
94CbcHeuristicRINS::generateCpp( FILE * fp) 
95{
96  CbcHeuristicRINS other;
97  fprintf(fp,"0#include \"CbcHeuristicRINS.hpp\"\n");
98  fprintf(fp,"3  CbcHeuristicRINS heuristicRINS(*cbcModel);\n");
99  CbcHeuristic::generateCpp(fp,"heuristicRINS");
100  if (howOften_!=other.howOften_)
101    fprintf(fp,"3  heuristicRINS.setHowOften(%d);\n",howOften_);
102  else
103    fprintf(fp,"4  heuristicRINS.setHowOften(%d);\n",howOften_);
104  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicRINS);\n");
105}
106
107// Copy constructor
108CbcHeuristicRINS::CbcHeuristicRINS(const CbcHeuristicRINS & rhs)
109:
110  CbcHeuristic(rhs),
111  numberSolutions_(rhs.numberSolutions_),
112  howOften_(rhs.howOften_),
113  numberSuccesses_(rhs.numberSuccesses_),
114  numberTries_(rhs.numberTries_),
115  stateOfFixing_(rhs.stateOfFixing_),
116  lastNode_(rhs.lastNode_)
117{
118  if (model_&&rhs.used_) {
119    int numberColumns = model_->solver()->getNumCols();
120    used_ = new char[numberColumns];
121    memcpy(used_,rhs.used_,numberColumns);
122  } else {
123    used_=NULL;
124  }
125}
126// Resets stuff if model changes
127void 
128CbcHeuristicRINS::resetModel(CbcModel * /*model*/)
129{
130  //CbcHeuristic::resetModel(model);
131  delete [] used_;
132  stateOfFixing_ = 0;
133  if (model_&&used_) {
134    int numberColumns = model_->solver()->getNumCols();
135    used_ = new char[numberColumns];
136    memset(used_,0,numberColumns);
137  } else {
138    used_=NULL;
139  }
140}
141/*
142  First tries setting a variable to better value.  If feasible then
143  tries setting others.  If not feasible then tries swaps
144  Returns 1 if solution, 0 if not */
145int
146CbcHeuristicRINS::solution(double & solutionValue,
147                         double * betterSolution)
148{
149  numCouldRun_++;
150  int returnCode=0;
151  const double * bestSolution = model_->bestSolution();
152  if (!bestSolution)
153    return 0; // No solution found yet
154  if (numberSolutions_<model_->getSolutionCount()) {
155    // new solution - add info
156    numberSolutions_=model_->getSolutionCount();
157
158    int numberIntegers = model_->numberIntegers();
159    const int * integerVariable = model_->integerVariable();
160 
161    int i;
162    for (i=0;i<numberIntegers;i++) {
163      int iColumn = integerVariable[i];
164      const OsiObject * object = model_->object(i);
165      // get original bounds
166      double originalLower;
167      double originalUpper;
168      getIntegerInformation( object,originalLower, originalUpper); 
169      double value=bestSolution[iColumn];
170      if (value<originalLower) {
171        value=originalLower;
172      } else if (value>originalUpper) {
173        value=originalUpper;
174      }
175      double nearest=floor(value+0.5);
176      // if away from lower bound mark that fact
177      if (nearest>originalLower) {
178        used_[iColumn]=1;
179      }
180    }
181  } 
182  int numberNodes=model_->getNodeCount();
183  if (howOften_==100) {
184    if (numberNodes<lastNode_+12)
185      return 0;
186    // Do at 50 and 100
187    if ((numberNodes>40&&numberNodes<=50)||(numberNodes>90&&numberNodes<100))
188      numberNodes=howOften_;
189  }
190  if ((numberNodes%howOften_)==0&&(model_->getCurrentPassNumber()==1||
191                                   model_->getCurrentPassNumber()==999999)) {
192    lastNode_=model_->getNodeCount();
193    OsiSolverInterface * solver = model_->solver();
194
195    int numberIntegers = model_->numberIntegers();
196    const int * integerVariable = model_->integerVariable();
197 
198    const double * currentSolution = solver->getColSolution();
199    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
200    int numberColumns=newSolver->getNumCols();
201    int numberContinuous = numberColumns-numberIntegers;
202
203    double primalTolerance;
204    solver->getDblParam(OsiPrimalTolerance,primalTolerance);
205   
206    int i;
207    int nFix=0;
208    for (i=0;i<numberIntegers;i++) {
209      int iColumn=integerVariable[i];
210      const OsiObject * object = model_->object(i);
211      // get original bounds
212      double originalLower;
213      double originalUpper;
214      getIntegerInformation( object,originalLower, originalUpper); 
215      double valueInt=bestSolution[iColumn];
216      if (valueInt<originalLower) {
217        valueInt=originalLower;
218      } else if (valueInt>originalUpper) {
219        valueInt=originalUpper;
220      }
221      if (fabs(currentSolution[iColumn]-valueInt)<10.0*primalTolerance) {
222        double nearest=floor(valueInt+0.5);
223        newSolver->setColLower(iColumn,nearest);
224        newSolver->setColUpper(iColumn,nearest);
225        nFix++;
226      }
227    }
228    int divisor=0;
229    if (5*nFix>numberIntegers) {
230      if (numberContinuous>2*numberIntegers&&nFix*10<numberColumns&&
231          ((!numRuns_&&numberTries_>2)||stateOfFixing_)) {
232#define RINS_FIX_CONTINUOUS
233#ifdef RINS_FIX_CONTINUOUS
234        const double * colLower = newSolver->getColLower();
235        //const double * colUpper = newSolver->getColUpper();
236        int nAtLb=0;
237        //double sumDj=0.0;
238        const double * dj = newSolver->getReducedCost();
239        double direction=newSolver->getObjSense();
240        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
241          if (!newSolver->isInteger(iColumn)) {
242            double value=bestSolution[iColumn];
243            if (value<colLower[iColumn]+1.0e-8) {
244              //double djValue = dj[iColumn]*direction;
245              nAtLb++;
246              //sumDj += djValue;
247            }
248          }
249        }
250        if (nAtLb) {
251          // fix some continuous
252          double * sort = new double[nAtLb];
253          int * which = new int [nAtLb];
254          //double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
255          int nFix2=0;
256          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
257            if (!newSolver->isInteger(iColumn)) {
258              double value=bestSolution[iColumn];
259              if (value<colLower[iColumn]+1.0e-8) {
260                double djValue = dj[iColumn]*direction;
261                if (djValue>1.0e-6) {
262                  sort[nFix2]=-djValue;
263                  which[nFix2++]=iColumn;
264                }
265              }
266            }
267          }
268          CoinSort_2(sort,sort+nFix2,which);
269          divisor=4;
270          if (stateOfFixing_>0)
271            divisor=stateOfFixing_;
272          else if (stateOfFixing_<-1)
273            divisor=(-stateOfFixing_)-1;
274          nFix2=CoinMin(nFix2,(numberColumns-nFix)/divisor);
275          for (int i=0;i<nFix2;i++) {
276            int iColumn = which[i];
277            newSolver->setColUpper(iColumn,colLower[iColumn]);
278          }
279          delete [] sort;
280          delete [] which;
281#ifdef CLP_INVESTIGATE2
282          printf("%d integers have same value, and %d continuous fixed at lb\n",
283                 nFix,nFix2);
284#endif
285        }
286#endif
287      }
288      //printf("%d integers have same value\n",nFix);
289      returnCode = smallBranchAndBound(newSolver,numberNodes_,betterSolution,solutionValue,
290                                         model_->getCutoff(),"CbcHeuristicRINS");
291      if (returnCode<0) {
292        returnCode=0; // returned on size
293        if (divisor)
294          stateOfFixing_ = - divisor; // say failed
295      } else {
296        numRuns_++;
297        if (divisor)
298          stateOfFixing_ =  divisor; // say small enough
299      }
300      if ((returnCode&1)!=0)
301        numberSuccesses_++;
302      //printf("return code %d",returnCode);
303      if ((returnCode&2)!=0) {
304        // could add cut
305        returnCode &= ~2;
306        //printf("could add cut with %d elements (if all 0-1)\n",nFix);
307      } else {
308        //printf("\n");
309      }
310    }
311
312    numberTries_++;
313    if ((numberTries_%10)==0&&numberSuccesses_*3<numberTries_)
314      howOften_ += static_cast<int> (howOften_*decayFactor_);
315    delete newSolver;
316  }
317  return returnCode;
318}
319// update model
320void CbcHeuristicRINS::setModel(CbcModel * model)
321{
322  model_ = model;
323  // Get a copy of original matrix
324  assert(model_->solver());
325  delete [] used_;
326  int numberColumns = model->solver()->getNumCols();
327  used_ = new char[numberColumns];
328  memset(used_,0,numberColumns);
329}
330// Default Constructor
331CbcHeuristicRENS::CbcHeuristicRENS() 
332  :CbcHeuristic()
333{
334  numberTries_=0;
335  whereFrom_=256+1;
336}
337
338// Constructor with model - assumed before cuts
339
340CbcHeuristicRENS::CbcHeuristicRENS(CbcModel & model)
341  :CbcHeuristic(model)
342{
343  numberTries_=0;
344  whereFrom_=256+1;
345}
346
347// Destructor
348CbcHeuristicRENS::~CbcHeuristicRENS ()
349{
350}
351
352// Clone
353CbcHeuristic *
354CbcHeuristicRENS::clone() const
355{
356  return new CbcHeuristicRENS(*this);
357}
358
359// Assignment operator
360CbcHeuristicRENS & 
361CbcHeuristicRENS::operator=( const CbcHeuristicRENS& rhs)
362{
363  if (this!=&rhs) {
364    CbcHeuristic::operator=(rhs);
365    numberTries_ = rhs.numberTries_;
366  }
367  return *this;
368}
369
370// Copy constructor
371CbcHeuristicRENS::CbcHeuristicRENS(const CbcHeuristicRENS & rhs)
372:
373  CbcHeuristic(rhs),
374  numberTries_(rhs.numberTries_)
375{
376}
377// Resets stuff if model changes
378void 
379CbcHeuristicRENS::resetModel(CbcModel * )
380{
381}
382int
383CbcHeuristicRENS::solution(double & solutionValue,
384                         double * betterSolution)
385{
386  int returnCode=0;
387  const double * bestSolution = model_->bestSolution();
388  if (numberTries_||(when()<2&&bestSolution))
389    return 0; 
390  numberTries_++;
391  OsiSolverInterface * solver = model_->solver();
392 
393  int numberIntegers = model_->numberIntegers();
394  const int * integerVariable = model_->integerVariable();
395 
396  const double * currentSolution = solver->getColSolution();
397  OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
398  const double * colLower = newSolver->getColLower();
399  const double * colUpper = newSolver->getColUpper();
400
401  double primalTolerance;
402  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
403   
404  int i;
405  int numberFixed=0;
406  int numberTightened=0;
407  int numberAtBound=0;
408  int numberColumns=newSolver->getNumCols();
409  int numberContinuous = numberColumns-numberIntegers;
410
411  for (i=0;i<numberIntegers;i++) {
412    int iColumn=integerVariable[i];
413    double value = currentSolution[iColumn];
414    double lower = colLower[iColumn];
415    double upper = colUpper[iColumn];
416    value = CoinMax(value,lower);
417    value = CoinMin(value,upper);
418#define RENS_FIX_ONLY_LOWER
419#ifndef RENS_FIX_ONLY_LOWER
420    if (fabs(value-floor(value+0.5))<1.0e-8) {
421      value = floor(value+0.5);
422      if (value==lower||value==upper)
423        numberAtBound++;
424      newSolver->setColLower(iColumn,value);
425      newSolver->setColUpper(iColumn,value);
426      numberFixed++;
427    } else if (colUpper[iColumn]-colLower[iColumn]>=2.0) {
428      numberTightened++;
429      newSolver->setColLower(iColumn,floor(value));
430      newSolver->setColUpper(iColumn,ceil(value));
431    }
432#else
433    if (fabs(value-floor(value+0.5))<1.0e-8&&
434        floor(value+0.5)==lower) {
435      value = floor(value+0.5);
436      numberAtBound++;
437      newSolver->setColLower(iColumn,value);
438      newSolver->setColUpper(iColumn,value);
439      numberFixed++;
440    } else if (colUpper[iColumn]-colLower[iColumn]>=2.0) {
441      numberTightened++;
442      if (fabs(value-floor(value+0.5))<1.0e-8) {
443        value = floor(value+0.5);
444        if (value<upper) {
445          newSolver->setColLower(iColumn,CoinMax(value-1.0,lower));
446          newSolver->setColUpper(iColumn,CoinMin(value+1.0,upper));
447        } else {
448          newSolver->setColLower(iColumn,upper-1.0);
449        }
450      } else {
451        newSolver->setColLower(iColumn,floor(value));
452        newSolver->setColUpper(iColumn,ceil(value));
453      }
454    }
455#endif
456  }
457  if (numberFixed>numberIntegers/5) {
458      if (numberContinuous>numberIntegers&&numberFixed<numberColumns/5) {
459#define RENS_FIX_CONTINUOUS
460#ifdef RENS_FIX_CONTINUOUS
461        const double * colLower = newSolver->getColLower();
462        //const double * colUpper = newSolver->getColUpper();
463        int nAtLb=0;
464        double sumDj=0.0;
465        const double * dj = newSolver->getReducedCost();
466        double direction=newSolver->getObjSense();
467        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
468          if (!newSolver->isInteger(iColumn)) {
469            double value=currentSolution[iColumn];
470            if (value<colLower[iColumn]+1.0e-8) {
471              double djValue = dj[iColumn]*direction;
472              nAtLb++;
473              sumDj += djValue;
474            }
475          }
476        }
477        if (nAtLb) {
478          // fix some continuous
479          double * sort = new double[nAtLb];
480          int * which = new int [nAtLb];
481          double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
482          int nFix2=0;
483          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
484            if (!newSolver->isInteger(iColumn)) {
485              double value=currentSolution[iColumn];
486              if (value<colLower[iColumn]+1.0e-8) {
487                double djValue = dj[iColumn]*direction;
488                if (djValue>threshold) {
489                  sort[nFix2]=-djValue;
490                  which[nFix2++]=iColumn;
491                }
492              }
493            }
494          }
495          CoinSort_2(sort,sort+nFix2,which);
496          nFix2=CoinMin(nFix2,(numberColumns-numberFixed)/2);
497          for (int i=0;i<nFix2;i++) {
498            int iColumn = which[i];
499            newSolver->setColUpper(iColumn,colLower[iColumn]);
500          }
501          delete [] sort;
502          delete [] which;
503#ifdef CLP_INVESTIGATE2
504          printf("%d integers fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
505                 numberFixed,numberTightened,numberAtBound,nFix2);
506#endif
507        }
508#endif
509      }
510#ifdef COIN_DEVELOP
511    printf("%d integers fixed and %d tightened\n",numberFixed,numberTightened);
512#endif
513    returnCode = smallBranchAndBound(newSolver,numberNodes_,betterSolution,solutionValue,
514                                     model_->getCutoff(),"CbcHeuristicRENS");
515    if (returnCode<0)
516      returnCode=0; // returned on size
517    //printf("return code %d",returnCode);
518    if ((returnCode&2)!=0) {
519      // could add cut
520      returnCode &= ~2;
521#ifdef COIN_DEVELOP
522      if (!numberTightened&&numberFixed==numberAtBound)
523        printf("could add cut with %d elements\n",numberFixed);
524#endif
525    } else {
526      //printf("\n");
527    }
528  }
529 
530  delete newSolver;
531  return returnCode;
532}
533// update model
534void CbcHeuristicRENS::setModel(CbcModel * model)
535{
536  model_ = model;
537}
538
539// Default Constructor
540CbcHeuristicDINS::CbcHeuristicDINS() 
541  :CbcHeuristic()
542{
543  numberSolutions_=0;
544  numberSuccesses_=0;
545  numberTries_=0;
546  howOften_=100;
547  decayFactor_ = 0.5; 
548  maximumKeepSolutions_=5;
549  numberKeptSolutions_=0;
550  numberIntegers_=-1;
551  localSpace_ = 10;
552  values_=NULL;
553}
554
555// Constructor with model - assumed before cuts
556
557CbcHeuristicDINS::CbcHeuristicDINS(CbcModel & model)
558  :CbcHeuristic(model)
559{
560  numberSolutions_=0;
561  numberSuccesses_=0;
562  numberTries_=0;
563  howOften_=100;
564  decayFactor_ = 0.5; 
565  assert(model.solver());
566  maximumKeepSolutions_=5;
567  numberKeptSolutions_=0;
568  numberIntegers_=-1;
569  localSpace_ = 10;
570  values_=NULL;
571}
572
573// Destructor
574CbcHeuristicDINS::~CbcHeuristicDINS ()
575{
576  for (int i=0;i<numberKeptSolutions_;i++)
577    delete [] values_[i];
578  delete [] values_;
579}
580
581// Clone
582CbcHeuristic *
583CbcHeuristicDINS::clone() const
584{
585  return new CbcHeuristicDINS(*this);
586}
587
588// Assignment operator
589CbcHeuristicDINS & 
590CbcHeuristicDINS::operator=( const CbcHeuristicDINS& rhs)
591{
592  if (this!=&rhs) {
593    CbcHeuristic::operator=(rhs);
594    numberSolutions_ = rhs.numberSolutions_;
595    howOften_ = rhs.howOften_;
596    numberSuccesses_ = rhs.numberSuccesses_;
597    numberTries_ = rhs.numberTries_;
598    for (int i=0;i<numberKeptSolutions_;i++)
599      delete [] values_[i];
600    delete [] values_;
601    maximumKeepSolutions_ = rhs.maximumKeepSolutions_;
602    numberKeptSolutions_ = rhs.numberKeptSolutions_;
603    numberIntegers_ = rhs.numberIntegers_;
604    localSpace_ = rhs.localSpace_;
605    if (model_&&rhs.values_) {
606      assert (numberIntegers_>=0);
607      values_ = new int * [maximumKeepSolutions_];
608      for (int i=0;i<maximumKeepSolutions_;i++)
609        values_[i]=CoinCopyOfArray(rhs.values_[i],numberIntegers_);
610    } else {
611      values_=NULL;
612    }
613  }
614  return *this;
615}
616
617// Create C++ lines to get to current state
618void 
619CbcHeuristicDINS::generateCpp( FILE * fp) 
620{
621  CbcHeuristicDINS other;
622  fprintf(fp,"0#include \"CbcHeuristicDINS.hpp\"\n");
623  fprintf(fp,"3  CbcHeuristicDINS heuristicDINS(*cbcModel);\n");
624  CbcHeuristic::generateCpp(fp,"heuristicDINS");
625  if (howOften_!=other.howOften_)
626    fprintf(fp,"3  heuristicDINS.setHowOften(%d);\n",howOften_);
627  else
628    fprintf(fp,"4  heuristicDINS.setHowOften(%d);\n",howOften_);
629  if (maximumKeepSolutions_!=other.maximumKeepSolutions_)
630    fprintf(fp,"3  heuristicDINS.setMaximumKeep(%d);\n",maximumKeepSolutions_);
631  else
632    fprintf(fp,"4  heuristicDINS.setMaximumKeep(%d);\n",maximumKeepSolutions_);
633  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicDINS);\n");
634}
635
636// Copy constructor
637CbcHeuristicDINS::CbcHeuristicDINS(const CbcHeuristicDINS & rhs)
638:
639  CbcHeuristic(rhs),
640  numberSolutions_(rhs.numberSolutions_),
641  howOften_(rhs.howOften_),
642  numberSuccesses_(rhs.numberSuccesses_),
643  numberTries_(rhs.numberTries_),
644  maximumKeepSolutions_(rhs.maximumKeepSolutions_),
645  numberKeptSolutions_(rhs.numberKeptSolutions_),
646  numberIntegers_(rhs.numberIntegers_),
647  localSpace_(rhs.localSpace_)
648{
649  if (model_&&rhs.values_) {
650    assert (numberIntegers_>=0);
651    values_ = new int * [maximumKeepSolutions_];
652    for (int i=0;i<maximumKeepSolutions_;i++)
653      values_[i]=CoinCopyOfArray(rhs.values_[i],numberIntegers_);
654  } else {
655    values_=NULL;
656  }
657}
658// Resets stuff if model changes
659void 
660CbcHeuristicDINS::resetModel(CbcModel * )
661{
662  //CbcHeuristic::resetModel(model);
663  for (int i=0;i<numberKeptSolutions_;i++)
664    delete [] values_[i];
665  delete [] values_;
666  numberKeptSolutions_=0;
667  numberIntegers_=-1;
668  numberSolutions_=0;
669  values_=NULL;
670}
671/*
672  First tries setting a variable to better value.  If feasible then
673  tries setting others.  If not feasible then tries swaps
674  Returns 1 if solution, 0 if not */
675int
676CbcHeuristicDINS::solution(double & solutionValue,
677                         double * betterSolution)
678{
679  numCouldRun_++;
680  int returnCode=0;
681  const double * bestSolution = model_->bestSolution();
682  if (!bestSolution)
683    return 0; // No solution found yet
684  if (numberSolutions_<model_->getSolutionCount()) {
685    // new solution - add info
686    numberSolutions_=model_->getSolutionCount();
687
688    int numberIntegers = model_->numberIntegers();
689    const int * integerVariable = model_->integerVariable();
690    if (numberIntegers_<0) {
691      numberIntegers_ = numberIntegers;
692      assert (!values_);
693      values_ = new int * [maximumKeepSolutions_];
694      for (int i=0;i<maximumKeepSolutions_;i++)
695        values_[i]=NULL;
696    } else {
697      assert (numberIntegers==numberIntegers_);
698    }
699    // move solutions (0 will be most recent)
700    {
701      int * temp = values_[maximumKeepSolutions_-1];
702      for (int i = maximumKeepSolutions_-1;i>0;i--)
703        values_[i]=values_[i-1];
704      if (!temp)
705        temp = new int [numberIntegers_];
706      values_[0]=temp;
707    }
708    int i;
709    for (i=0;i<numberIntegers;i++) {
710      int iColumn = integerVariable[i];
711      double value=bestSolution[iColumn];
712      double nearest=floor(value+0.5);
713      values_[0][i]=static_cast<int> (nearest);
714    }
715    numberKeptSolutions_ = CoinMin(numberKeptSolutions_+1,maximumKeepSolutions_);
716  } 
717  int finalReturnCode=0;
718  if (((model_->getNodeCount()%howOften_)==howOften_/2||!model_->getNodeCount())&&(model_->getCurrentPassNumber()==1||model_->getCurrentPassNumber()==999999)) {
719    OsiSolverInterface * solver = model_->solver();
720
721    int numberIntegers = model_->numberIntegers();
722    const int * integerVariable = model_->integerVariable();
723 
724    const double * currentSolution = solver->getColSolution();
725    int localSpace = localSpace_;
726    // 0 means finished but no solution, 1 solution, 2 node limit
727    int status=-1;
728    double cutoff = model_->getCutoff();
729    while(status) {
730      status=0;
731      OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
732      const double * colLower = solver->getColLower();
733      const double * colUpper = solver->getColUpper();
734     
735      double primalTolerance;
736      solver->getDblParam(OsiPrimalTolerance,primalTolerance);
737      const double * continuousSolution = newSolver->getColSolution();
738      // Space for added constraint
739      double * element = new double [numberIntegers];
740      int * column = new int [numberIntegers];
741      int i;
742      int nFix=0;
743      int nCouldFix=0;
744      int nCouldFix2=0;
745      int nBound=0;
746      int nEl=0;
747      double bias=localSpace;
748      int okSame = numberKeptSolutions_-1;
749      for (i=0;i<numberIntegers;i++) {
750        int iColumn=integerVariable[i];
751        const OsiObject * object = model_->object(i);
752        // get original bounds
753        double originalLower;
754        double originalUpper;
755        getIntegerInformation( object,originalLower, originalUpper); 
756        double valueInt=bestSolution[iColumn];
757        if (valueInt<originalLower) {
758          valueInt=originalLower;
759        } else if (valueInt>originalUpper) {
760          valueInt=originalUpper;
761        }
762        int intValue = static_cast<int> (floor(valueInt+0.5));
763        double currentValue = currentSolution[iColumn];
764        double currentLower = colLower[iColumn];
765        double currentUpper = colUpper[iColumn];
766        if (fabs(valueInt-currentValue)>=0.5) {
767          // Re-bound
768          nBound++;
769          if (intValue>=currentValue) {
770            currentLower = CoinMax(currentLower,ceil(2*currentValue-intValue));
771            currentUpper = intValue;
772          } else {
773            currentLower = intValue;
774            currentUpper = CoinMin(currentUpper,floor(2*currentValue-intValue));
775          }
776          newSolver->setColLower(iColumn,currentLower);
777          newSolver->setColUpper(iColumn,currentUpper);
778        } else {
779          // See if can fix
780          bool canFix=false;
781          double continuousValue = continuousSolution[iColumn];
782          if (fabs(currentValue-valueInt)<10.0*primalTolerance) {
783            if (currentUpper-currentLower>1.0) {
784              // General integer variable
785              canFix=true;
786            } else if(fabs(continuousValue-valueInt)<10.0*primalTolerance) {
787              int nSame=1;
788              //assert (intValue==values_[0][i]);
789              for (int k=1;k<numberKeptSolutions_;k++) {
790                if (intValue==values_[k][i])
791                  nSame++;
792              }
793              if (nSame>=okSame) {
794                // can fix
795                canFix=true;
796              } else {
797                nCouldFix++;
798              }
799            } else {
800              nCouldFix2++;
801            }
802          }
803          if (canFix) {
804            newSolver->setColLower(iColumn,intValue);
805            newSolver->setColUpper(iColumn,intValue);
806            nFix++;
807          } else {
808            if (currentUpper-currentLower>1.0) {
809              // General integer variable
810              currentLower = floor(currentValue);
811              if (intValue>=currentLower&&intValue<=currentLower+1) {
812                newSolver->setColLower(iColumn,currentLower);
813                newSolver->setColUpper(iColumn,currentLower+1.0);
814              } else {
815                // fix
816                double value;
817                if (intValue<currentLower)
818                  value = currentLower;
819                else
820                  value = currentLower+1;
821                newSolver->setColLower(iColumn,value);
822                newSolver->setColUpper(iColumn,value);
823                nFix++;
824              }
825            } else {
826              // 0-1 (ish)
827              column[nEl]=iColumn;
828              if (intValue==currentLower) {
829                bias += currentLower;
830                element[nEl++]=1.0;
831              } else if (intValue==currentUpper) {
832                bias += currentUpper;
833                element[nEl++]=-1.0;
834              } else {
835                printf("bad DINS logic\n");
836                abort();
837              }
838            }
839          }
840        }
841      }
842      char generalPrint[200];
843      sprintf(generalPrint,
844              "%d fixed, %d same as cont/int, %d same as int - %d bounded %d in cut\n",
845              nFix,nCouldFix,nCouldFix2,nBound,nEl);
846      model_->messageHandler()->message(CBC_FPUMP2,model_->messages())
847        << generalPrint
848        <<CoinMessageEol;
849      if (nFix>numberIntegers/10) {
850#if 0
851        newSolver->initialSolve();
852        printf("obj %g\n",newSolver->getObjValue());
853        for (i=0;i<numberIntegers;i++) {
854          int iColumn=integerVariable[i];
855          printf("%d new bounds %g %g - solutions %g %g\n",
856                 iColumn,newSolver->getColLower()[iColumn],
857                 newSolver->getColUpper()[iColumn],
858                 bestSolution[iColumn],
859                 currentSolution[iColumn]);
860        }
861#endif
862        if (nEl>0)
863          newSolver->addRow(nEl,column,element,-COIN_DBL_MAX,bias);
864        //printf("%d integers have same value\n",nFix);
865        returnCode = smallBranchAndBound(newSolver,numberNodes_,betterSolution,solutionValue,
866                                         cutoff,"CbcHeuristicDINS");
867        if (returnCode<0) {
868          returnCode=0; // returned on size
869          status=0;
870        } else {
871          numRuns_++;
872          if ((returnCode&2)!=0) {
873            // could add cut as complete search
874            returnCode &= ~2;
875            if ((returnCode&1)!=0) {
876              numberSuccesses_++;
877              status=1;
878            } else {
879              // no solution
880              status=0;
881            }
882          } else {
883            if ((returnCode&1)!=0) {
884              numberSuccesses_++;
885              status=1;
886            } else {
887              // no solution but node limit
888              status=2;
889              if (nEl)
890                localSpace -= 5;
891              else
892                localSpace = -1;
893              if (localSpace<0)
894                status=0;
895            }
896          }
897          if ((returnCode&1)!=0) {
898            cutoff = CoinMin(cutoff,solutionValue-model_->getCutoffIncrement());
899            finalReturnCode=1;
900          }
901        }
902      }
903      delete [] element;
904      delete [] column;
905      delete newSolver;
906    }
907    numberTries_++;
908    if ((numberTries_%10)==0&&numberSuccesses_*3<numberTries_)
909      howOften_ += static_cast<int> (howOften_*decayFactor_);
910  }
911  return finalReturnCode;
912}
913// update model
914void CbcHeuristicDINS::setModel(CbcModel * model)
915{
916  model_ = model;
917  // Get a copy of original matrix
918  assert(model_->solver());
919  for (int i=0;i<numberKeptSolutions_;i++)
920    delete [] values_[i];
921  delete [] values_;
922  numberKeptSolutions_=0;
923  numberIntegers_=-1;
924  numberSolutions_=0;
925  values_=NULL;
926}
927
928// Default Constructor
929CbcHeuristicVND::CbcHeuristicVND() 
930  :CbcHeuristic()
931{
932  numberSolutions_=0;
933  numberSuccesses_=0;
934  numberTries_=0;
935  lastNode_=-999999;
936  howOften_=100;
937  decayFactor_ = 0.5; 
938  baseSolution_=NULL;
939  whereFrom_=1+8+255*256;
940  stepSize_=0;
941  k_=0;
942  kmax_=0;
943  nDifferent_=0;
944}
945
946// Constructor with model - assumed before cuts
947
948CbcHeuristicVND::CbcHeuristicVND(CbcModel & model)
949  :CbcHeuristic(model)
950{
951  numberSolutions_=0;
952  numberSuccesses_=0;
953  numberTries_=0;
954  lastNode_=-999999;
955  howOften_=100;
956  decayFactor_ = 0.5; 
957  assert(model.solver());
958  int numberColumns = model.solver()->getNumCols();
959  baseSolution_ = new double [numberColumns];
960  memset(baseSolution_,0,numberColumns*sizeof(double));
961  whereFrom_=1+8+255*256;
962  stepSize_=0;
963  k_=0;
964  kmax_=0;
965  nDifferent_=0;
966}
967
968// Destructor
969CbcHeuristicVND::~CbcHeuristicVND ()
970{
971  delete [] baseSolution_;
972}
973
974// Clone
975CbcHeuristic *
976CbcHeuristicVND::clone() const
977{
978  return new CbcHeuristicVND(*this);
979}
980
981// Assignment operator
982CbcHeuristicVND & 
983CbcHeuristicVND::operator=( const CbcHeuristicVND& rhs)
984{
985  if (this!=&rhs) {
986    CbcHeuristic::operator=(rhs);
987    numberSolutions_ = rhs.numberSolutions_;
988    howOften_ = rhs.howOften_;
989    numberSuccesses_ = rhs.numberSuccesses_;
990    numberTries_ = rhs.numberTries_;
991    lastNode_=rhs.lastNode_;
992    delete [] baseSolution_;
993    if (model_&&rhs.baseSolution_) {
994      int numberColumns = model_->solver()->getNumCols();
995      baseSolution_ = new double [numberColumns];
996      memcpy(baseSolution_,rhs.baseSolution_,numberColumns*sizeof(double));
997    } else {
998      baseSolution_=NULL;
999    }
1000    stepSize_=rhs.stepSize_;
1001    k_=rhs.k_;
1002    kmax_=rhs.kmax_;
1003    nDifferent_=rhs.nDifferent_;
1004  }
1005  return *this;
1006}
1007
1008// Create C++ lines to get to current state
1009void 
1010CbcHeuristicVND::generateCpp( FILE * fp) 
1011{
1012  CbcHeuristicVND other;
1013  fprintf(fp,"0#include \"CbcHeuristicVND.hpp\"\n");
1014  fprintf(fp,"3  CbcHeuristicVND heuristicVND(*cbcModel);\n");
1015  CbcHeuristic::generateCpp(fp,"heuristicVND");
1016  if (howOften_!=other.howOften_)
1017    fprintf(fp,"3  heuristicVND.setHowOften(%d);\n",howOften_);
1018  else
1019    fprintf(fp,"4  heuristicVND.setHowOften(%d);\n",howOften_);
1020  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicVND);\n");
1021}
1022
1023// Copy constructor
1024CbcHeuristicVND::CbcHeuristicVND(const CbcHeuristicVND & rhs)
1025:
1026  CbcHeuristic(rhs),
1027  numberSolutions_(rhs.numberSolutions_),
1028  howOften_(rhs.howOften_),
1029  numberSuccesses_(rhs.numberSuccesses_),
1030  numberTries_(rhs.numberTries_),
1031  lastNode_(rhs.lastNode_)
1032{
1033  if (model_&&rhs.baseSolution_) {
1034    int numberColumns = model_->solver()->getNumCols();
1035    baseSolution_ = new double [numberColumns];
1036    memcpy(baseSolution_,rhs.baseSolution_,numberColumns*sizeof(double));
1037  } else {
1038    baseSolution_=NULL;
1039  }
1040  stepSize_=rhs.stepSize_;
1041  k_=rhs.k_;
1042  kmax_=rhs.kmax_;
1043  nDifferent_=rhs.nDifferent_;
1044}
1045// Resets stuff if model changes
1046void 
1047CbcHeuristicVND::resetModel(CbcModel * /*model*/)
1048{
1049  //CbcHeuristic::resetModel(model);
1050  delete [] baseSolution_;
1051  if (model_&&baseSolution_) {
1052    int numberColumns = model_->solver()->getNumCols();
1053    baseSolution_ = new double [numberColumns];
1054    memset(baseSolution_,0,numberColumns*sizeof(double));
1055  } else {
1056    baseSolution_=NULL;
1057  }
1058}
1059/*
1060  First tries setting a variable to better value.  If feasible then
1061  tries setting others.  If not feasible then tries swaps
1062  Returns 1 if solution, 0 if not */
1063int
1064CbcHeuristicVND::solution(double & solutionValue,
1065                         double * betterSolution)
1066{
1067  numCouldRun_++;
1068  int returnCode=0;
1069  const double * bestSolution = model_->bestSolution();
1070  if (!bestSolution)
1071    return 0; // No solution found yet
1072  if (numberSolutions_<model_->getSolutionCount()) {
1073    // new solution - add info
1074    numberSolutions_=model_->getSolutionCount();
1075
1076    int numberIntegers = model_->numberIntegers();
1077    const int * integerVariable = model_->integerVariable();
1078 
1079    int i;
1080    for (i=0;i<numberIntegers;i++) {
1081      int iColumn = integerVariable[i];
1082      const OsiObject * object = model_->object(i);
1083      // get original bounds
1084      double originalLower;
1085      double originalUpper;
1086      getIntegerInformation( object,originalLower, originalUpper); 
1087      double value=bestSolution[iColumn];
1088      if (value<originalLower) {
1089        value=originalLower;
1090      } else if (value>originalUpper) {
1091        value=originalUpper;
1092      }
1093    }
1094  } 
1095  int numberNodes=model_->getNodeCount();
1096  if (howOften_==100) {
1097    if (numberNodes<lastNode_+12)
1098      return 0;
1099    // Do at 50 and 100
1100    if ((numberNodes>40&&numberNodes<=50)||(numberNodes>90&&numberNodes<100))
1101      numberNodes=howOften_;
1102  }
1103  if ((numberNodes%howOften_)==0&&(model_->getCurrentPassNumber()==1||
1104                                   model_->getCurrentPassNumber()==999999)) {
1105    lastNode_=model_->getNodeCount();
1106    OsiSolverInterface * solver = model_->solver();
1107
1108    int numberIntegers = model_->numberIntegers();
1109    const int * integerVariable = model_->integerVariable();
1110 
1111    const double * currentSolution = solver->getColSolution();
1112    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
1113    //const double * colLower = newSolver->getColLower();
1114    //const double * colUpper = newSolver->getColUpper();
1115
1116    double primalTolerance;
1117    solver->getDblParam(OsiPrimalTolerance,primalTolerance);
1118
1119    // Sort on distance
1120    double * distance = new double [numberIntegers];
1121    int * which = new int [numberIntegers];
1122   
1123    int i;
1124    int nFix=0;
1125    double tolerance = 10.0*primalTolerance;
1126    for (i=0;i<numberIntegers;i++) {
1127      int iColumn=integerVariable[i];
1128      const OsiObject * object = model_->object(i);
1129      // get original bounds
1130      double originalLower;
1131      double originalUpper;
1132      getIntegerInformation( object,originalLower, originalUpper); 
1133      double valueInt=bestSolution[iColumn];
1134      if (valueInt<originalLower) {
1135        valueInt=originalLower;
1136      } else if (valueInt>originalUpper) {
1137        valueInt=originalUpper;
1138      }
1139      baseSolution_[iColumn]=currentSolution[iColumn];
1140      distance[i]=fabs(currentSolution[iColumn]-valueInt);
1141      which[i]=i;
1142      if (fabs(currentSolution[iColumn]-valueInt)<tolerance) 
1143        nFix++;
1144    }
1145    CoinSort_2(distance,distance+numberIntegers,which);
1146    nDifferent_=numberIntegers-nFix;
1147    stepSize_=nDifferent_/10;
1148    k_ = stepSize_;
1149    //nFix = numberIntegers-stepSize_;
1150    for (i=0;i<nFix;i++) {
1151      int j=which[i];
1152      int iColumn=integerVariable[j];
1153      const OsiObject * object = model_->object(i);
1154      // get original bounds
1155      double originalLower;
1156      double originalUpper;
1157      getIntegerInformation( object,originalLower, originalUpper); 
1158      double valueInt=bestSolution[iColumn];
1159      if (valueInt<originalLower) {
1160        valueInt=originalLower;
1161      } else if (valueInt>originalUpper) {
1162        valueInt=originalUpper;
1163      }
1164      double nearest=floor(valueInt+0.5);
1165      newSolver->setColLower(iColumn,nearest);
1166      newSolver->setColUpper(iColumn,nearest);
1167    }
1168    delete [] distance;
1169    delete [] which;
1170    if (nFix>numberIntegers/5) {
1171      //printf("%d integers have samish value\n",nFix);
1172      returnCode = smallBranchAndBound(newSolver,numberNodes_,betterSolution,solutionValue,
1173                                         model_->getCutoff(),"CbcHeuristicVND");
1174      if (returnCode<0)
1175        returnCode=0; // returned on size
1176      else
1177        numRuns_++;
1178      if ((returnCode&1)!=0)
1179        numberSuccesses_++;
1180      //printf("return code %d",returnCode);
1181      if ((returnCode&2)!=0) {
1182        // could add cut
1183        returnCode &= ~2;
1184        //printf("could add cut with %d elements (if all 0-1)\n",nFix);
1185      } else {
1186        //printf("\n");
1187      }
1188      numberTries_++;
1189      if ((numberTries_%10)==0&&numberSuccesses_*3<numberTries_)
1190        howOften_ += static_cast<int> (howOften_*decayFactor_);
1191    }
1192
1193    delete newSolver;
1194  }
1195  return returnCode;
1196}
1197// update model
1198void CbcHeuristicVND::setModel(CbcModel * model)
1199{
1200  model_ = model;
1201  // Get a copy of original matrix
1202  assert(model_->solver());
1203  delete [] baseSolution_;
1204  int numberColumns = model->solver()->getNumCols();
1205  baseSolution_ = new double [numberColumns];
1206  memset(baseSolution_,0,numberColumns*sizeof(double));
1207}
1208
1209 
Note: See TracBrowser for help on using the repository browser.