source: stable/0.2/Couenne/src/main/BonNlpHeuristic.cpp @ 159

Last change on this file since 159 was 159, checked in by pbelotti, 11 years ago

created new stable branch 0.2 from trunk (rev. 157)

File size: 9.8 KB
Line 
1/* $Id: BonNlpHeuristic.cpp 154 2009-06-16 18:52:53Z pbelotti $ */
2// (C) Copyright International Business Machines Corporation 2007
3// All Rights Reserved.
4// This code is published under the Common Public License.
5//
6// Authors :
7// Pierre Bonami, International Business Machines Corporation
8//
9// Date : 04/09/2007
10
11#include "BonNlpHeuristic.hpp"
12#include "BonCouenneInterface.hpp"
13#include "CouenneObject.hpp"
14#include "CouenneProblem.hpp"
15#include "CbcCutGenerator.hpp"
16#include "CbcBranchActual.hpp"
17#include "BonAuxInfos.hpp"
18#include "CoinHelperFunctions.hpp"
19
20#include "CouenneCutGenerator.hpp"
21#include "CouenneProblem.hpp"
22
23namespace Bonmin{
24  NlpSolveHeuristic::NlpSolveHeuristic():
25    CbcHeuristic(),
26    nlp_(NULL),
27    hasCloned_(false),
28    maxNlpInf_(maxNlpInf_0),
29    numberSolvePerLevel_(-1),
30    couenne_(NULL){
31    setHeuristicName("NlpSolveHeuristic");
32  }
33 
34  NlpSolveHeuristic::NlpSolveHeuristic(CbcModel & model, OsiSolverInterface &nlp, bool cloneNlp, CouenneProblem * couenne):
35  CbcHeuristic(model), nlp_(&nlp), hasCloned_(cloneNlp),maxNlpInf_(maxNlpInf_0),
36  numberSolvePerLevel_(-1),
37  couenne_(couenne){
38    setHeuristicName("NlpSolveHeuristic");
39    if(cloneNlp)
40      nlp_ = nlp.clone();
41  }
42 
43  NlpSolveHeuristic::NlpSolveHeuristic(const NlpSolveHeuristic & other):
44  CbcHeuristic(other), nlp_(other.nlp_), 
45  hasCloned_(other.hasCloned_),
46  maxNlpInf_(other.maxNlpInf_),
47  numberSolvePerLevel_(other.numberSolvePerLevel_),
48  couenne_(other.couenne_){
49    if(hasCloned_ && nlp_ != NULL)
50      nlp_ = other.nlp_->clone();
51  }
52 
53  CbcHeuristic * 
54  NlpSolveHeuristic::clone() const{
55    return new NlpSolveHeuristic(*this);
56  }
57 
58  NlpSolveHeuristic &
59  NlpSolveHeuristic::operator=(const NlpSolveHeuristic & rhs){
60    if(this != &rhs){
61      CbcHeuristic::operator=(rhs);
62      if(hasCloned_ && nlp_)
63        delete nlp_;
64     
65      hasCloned_ = rhs.hasCloned_;
66      if(nlp_ != NULL){
67        if(hasCloned_)
68          nlp_ = rhs.nlp_->clone();
69        else
70          nlp_ = rhs.nlp_;
71      }
72    }
73    maxNlpInf_ = rhs.maxNlpInf_;
74    numberSolvePerLevel_ = rhs.numberSolvePerLevel_;
75    couenne_ = rhs.couenne_;
76    return *this;
77  }
78 
79  NlpSolveHeuristic::~NlpSolveHeuristic(){
80    if(hasCloned_)
81      delete nlp_;
82    nlp_ = NULL;
83  }
84 
85  void
86  NlpSolveHeuristic::setNlp(OsiSolverInterface &nlp, bool cloneNlp){
87    if(hasCloned_ && nlp_ != NULL)
88      delete nlp_;
89    hasCloned_ = cloneNlp;
90    if(cloneNlp)
91      nlp_ = nlp.clone();
92    else
93      nlp_ = &nlp;
94  }
95 
96  void
97  NlpSolveHeuristic::setCouenneProblem(CouenneProblem * couenne){
98    couenne_ = couenne;}
99
100
101  int
102  NlpSolveHeuristic::solution( double & objectiveValue, double * newSolution){
103    OsiSolverInterface * solver = model_->solver();
104
105    OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
106    BabInfo * babInfo = dynamic_cast<BabInfo *> (auxInfo);
107
108    if(babInfo){
109      babInfo->setHasNlpSolution(false);
110      if(babInfo->infeasibleNode()){
111        return 0;
112      }
113    }
114
115    // if too deep in the BB tree, only run NLP heuristic if
116    // feasibility is low
117    bool too_deep = false;
118
119    // check depth
120    if (numberSolvePerLevel_ > -1){
121      if (numberSolvePerLevel_ == 0) 
122        return 0;
123
124      const int depth = (model_ -> currentNode ()) ? model_ -> currentNode () -> depth () : 0;
125
126      //if (CoinDrand48 () > pow (2., numberSolvePerLevel_ - depth))
127      if (CoinDrand48 () > 1. / CoinMax
128          (1., (double) ((depth - numberSolvePerLevel_) * (depth - numberSolvePerLevel_))))
129        too_deep = true;
130    }
131
132    if (too_deep)
133      return 0;
134
135    double *lower = new double [couenne_ -> nVars ()];
136    double *upper = new double [couenne_ -> nVars ()];
137
138    CoinFillN (lower, couenne_ -> nVars (), -COUENNE_INFINITY);
139    CoinFillN (upper, couenne_ -> nVars (),  COUENNE_INFINITY);
140
141    CoinCopyN (solver->getColLower(), nlp_ -> getNumCols (), lower);
142    CoinCopyN (solver->getColUpper(), nlp_ -> getNumCols (), upper);
143
144    /*printf ("-- int candidate, before: ");
145    for (int i=0; i<couenne_ -> nOrig (); i++)
146      printf ("[%g %g] ", lower [i], upper [i]);
147      printf ("\n");*/
148
149    const double * solution = solver->getColSolution();
150    OsiBranchingInformation info (solver, true);
151    const int & numberObjects = model_->numberObjects();
152    OsiObject ** objects = model_->objects();
153    double maxInfeasibility = 0;
154
155    bool haveRoundedIntVars = false;
156
157    for(int i = 0 ; i < numberObjects ; i++){
158      CouenneObject * couObj = dynamic_cast <CouenneObject *> (objects [i]);
159      if (couObj) {
160        if (too_deep) { // only test infeasibility if BB level is high
161          int dummy;
162          double infeas;
163          maxInfeasibility = CoinMax ( maxInfeasibility, infeas = couObj->infeasibility(&info, dummy));
164          if(maxInfeasibility > maxNlpInf_){
165            delete [] lower;
166            delete [] upper;
167            return 0;
168          }
169        }
170      } else {
171
172        OsiSimpleInteger * intObj = dynamic_cast<OsiSimpleInteger *>(objects[i]);
173
174        if (intObj) {
175          const int & i = intObj -> columnNumber ();
176          // Round the variable in the solver
177          double value = solution [i];
178          if (value < lower[i])
179            value = lower[i];
180          else if (value > upper[i])
181            value = upper[i];
182
183          double rounded = floor (value + 0.5);
184
185          if (fabs (value - rounded) > COUENNE_EPS) {
186            haveRoundedIntVars = true;
187            //value = rounded;
188          }
189
190          // fix bounds anyway, if a better candidate is not found
191          // below at least we have an integer point
192          //lower[i] = upper[i] = value;
193        }
194        else{
195           throw CoinError("Bonmin::NlpSolveHeuristic","solution",
196                           "Unknown object.");
197        }
198      }
199    }
200
201    // if here, it means the infeasibility is not too high. Generate a
202    // better integer point as there are rounded integer variables
203
204    bool skipOnInfeasibility = false;
205
206    double *Y = new double [couenne_ -> nVars ()];
207    CoinFillN (Y, couenne_ -> nVars (), 0.);
208    CoinCopyN (solution, nlp_ -> getNumCols (), Y);
209
210    /*printf ("-- int candidate, upon call: ");
211    for (int i=0; i<couenne_ -> nOrig (); i++)
212      if (couenne_ -> Var (i) -> isInteger ())
213        printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
214      else printf ("%g ", Y [i]);
215      printf ("\n");*/
216
217    if (haveRoundedIntVars) // create "good" integer candidate for Ipopt
218      skipOnInfeasibility = (couenne_ -> getIntegerCandidate (solution, Y, lower, upper) < 0);
219
220    /*printf ("-- int candidate, after: ");
221    for (int i=0; i<couenne_ -> nOrig (); i++)
222      if (couenne_ -> Var (i) -> isInteger ())
223        printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
224      else printf ("%g ", Y [i]);
225      printf ("\n");*/
226
227    bool foundSolution = false;
228
229    if (haveRoundedIntVars && skipOnInfeasibility) 
230      // no integer initial point could be found, make up some random rounding
231
232      for (int i = couenne_ -> nOrigVars (); i--;) 
233
234        if (couenne_ -> Var (i) -> isDefinedInteger ())
235          lower [i] = upper [i] = Y [i] = 
236            (CoinDrand48 () < 0.5) ? 
237            floor (Y [i] + COUENNE_EPS) : 
238            ceil  (Y [i] - COUENNE_EPS);
239
240        else if (lower [i] > upper [i]) { 
241
242          // sanity check (should avoid problems in ex1263 with
243          // couenne.opt.obbt)
244
245          double swap = lower [i];
246          lower [i] = upper [i];
247          upper [i] = swap;
248        }
249
250
251    {
252      //        printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
253
254      /*printf ("int candidate: ");
255        for (int i=0; i<couenne_ -> nOrig (); i++)
256        if (couenne_ -> Var (i) -> isInteger ())
257        printf ("[%g <%g> %g] ", lower [i], Y [i], upper [i]);
258        else printf ("%g ", Y [i]);
259        printf ("\n");*/
260
261      // Now set column bounds and solve the NLP with starting point
262      double * saveColLower = CoinCopyOfArray (nlp_ -> getColLower (), nlp_ -> getNumCols ());
263      double * saveColUpper = CoinCopyOfArray (nlp_ -> getColUpper (), nlp_ -> getNumCols ());
264
265      for (int i = nlp_ -> getNumCols (); i--;) {
266
267        if (lower [i] > upper [i]) {
268          double swap = lower [i];
269          lower [i] = upper [i];
270          upper [i] = swap;
271        }
272
273        if      (Y [i] < lower [i]) Y [i] = lower [i];
274        else if (Y [i] > upper [i]) Y [i] = upper [i];
275      }
276
277
278      nlp_ -> setColLower    (lower);
279      nlp_ -> setColUpper    (upper);
280      nlp_ -> setColSolution (Y);
281
282      // apply NLP solver /////////////////////////////////
283      try {
284        nlp_ -> initialSolve ();
285      }
286      catch (TNLPSolver::UnsolvedError *E) {}
287
288      double obj = (nlp_ -> isProvenOptimal()) ? nlp_ -> getObjValue (): COIN_DBL_MAX;
289
290      if (nlp_ -> isProvenOptimal () &&
291          couenne_ -> checkNLP (nlp_ -> getColSolution (), obj, true) && // true for recomputing obj
292          (obj < couenne_ -> getCutOff ())) {
293
294        // store solution in Aux info
295
296        const int nVars = solver->getNumCols();
297        double* tmpSolution = new double [nVars];
298        CoinCopyN (nlp_ -> getColSolution(), nlp_ -> getNumCols(), tmpSolution);
299
300        //Get correct values for all auxiliary variables
301        CouenneInterface * couenne = dynamic_cast <CouenneInterface *> (nlp_);
302
303        if (couenne)
304          couenne_ -> getAuxs (tmpSolution);
305
306        if (babInfo){
307          babInfo->setNlpSolution (tmpSolution, nVars, obj);
308          babInfo->setHasNlpSolution (true);
309        }
310
311        if (obj < objectiveValue) { // found better solution?
312
313          const CouNumber
314            *lb = solver -> getColLower (),
315            *ub = solver -> getColUpper ();
316
317          // check bounds once more after getAux. This avoids false
318          // asserts in CbcModel.cpp:8305 on integerTolerance violated
319          for (int i=0; i < nVars; i++, lb++, ub++) {
320
321            CouNumber &t = tmpSolution [i];
322            if      (t < *lb) t = *lb;
323            else if (t > *ub) t = *ub;
324          }
325         
326          //printf ("new cutoff %g from BonNlpHeuristic\n", obj);
327          couenne_ -> setCutOff (obj);
328          foundSolution = true;
329          objectiveValue = obj;
330          CoinCopyN (tmpSolution, nVars, newSolution);
331        }
332        delete [] tmpSolution;
333      }
334
335      nlp_->setColLower (saveColLower);
336      nlp_->setColUpper (saveColUpper);
337
338      delete [] saveColLower;
339      delete [] saveColUpper;
340    }
341
342    delete [] Y;
343
344    delete [] lower;
345    delete [] upper;
346
347    return foundSolution;
348  }
349}/** Ends namespace Bonmin.*/
Note: See TracBrowser for help on using the repository browser.