source: stable/0.1/Couenne/src/main/BonNlpHeuristic.cpp @ 134

Last change on this file since 134 was 134, checked in by pbelotti, 13 years ago

merged trunk revision with Ipopt exception handling

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