source: trunk/Bonmin/src/Algorithms/OaGenerators/BonOaNlpOptim.cpp @ 1360

Last change on this file since 1360 was 1360, checked in by pbonami, 12 years ago

Make "bonmin." prefix a variable

File size: 7.3 KB
Line 
1// (C) Copyright Carnegie Mellon University 2005
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// P. Bonami, Carnegie Mellon University
7//
8// Date :  05/26/2005
9
10#include "BonOaNlpOptim.hpp"
11#include "OsiAuxInfo.hpp"
12
13namespace Bonmin
14{
15/// Default constructor
16  OaNlpOptim::OaNlpOptim(OsiTMINLPInterface * si,
17      int maxDepth, bool addOnlyViolated, bool global)
18      :
19      CglCutGenerator(),
20      nlp_(si),
21      maxDepth_(maxDepth),
22      nSolve_(0),
23      addOnlyViolated_(addOnlyViolated),
24      global_(global)
25  {
26    handler_ = new CoinMessageHandler();
27    handler_ -> setLogLevel(1);
28    messages_ = OaMessages();
29  }
30
31  OaNlpOptim::OaNlpOptim(BabSetupBase &b):
32      CglCutGenerator(),
33      nlp_(b.nonlinearSolver()),
34      maxDepth_(1000),
35      nSolve_(0)
36  {
37    int ivalue;
38    b.options()->GetEnumValue("add_only_violated_oa", ivalue,b.prefix());
39    addOnlyViolated_ = ivalue;
40    b.options()->GetEnumValue("oa_cuts_scope", ivalue,b.prefix());
41    global_ = ivalue;
42
43    b.options()->GetIntegerValue("nlp_solve_max_depth", maxDepth_,b.prefix());
44    b.options()->GetNumericValue("nlp_solves_per_depth", solves_per_level_,b.prefix());
45    handler_ = new CoinMessageHandler();
46    handler_ -> setLogLevel(1);
47    messages_ = OaMessages();
48  }
49/// Assign an OsiTMINLPInterface
50  void
51  OaNlpOptim::assignInterface(OsiTMINLPInterface * si)
52
53  {
54    nlp_ = si;
55  }
56  static int nCalls = 0;
57/// cut generation method
58  void
59  OaNlpOptim::generateCuts( const OsiSolverInterface & si, OsiCuts & cs,
60      const CglTreeInfo info) const
61  {
62    if (nlp_ == NULL) {
63      CoinError("Error in cut generator for outer approximation no ipopt NLP assigned", "generateCuts", "OaNlpOptim");
64    }
65
66    int numcols = nlp_->getNumCols();
67
68    //Get the continuous solution
69    //const double *colsol = si.getColSolution();
70    //Check for integer feasibility
71    nCalls++;
72    double rand = CoinDrand48();
73    if (info.level > maxDepth_ || info.pass > 0 || !info.inTree  ||
74        pow(2.,-info.level)*solves_per_level_ <= rand)
75      return;
76    //Fix the variable which have to be fixed, after having saved the bounds
77    double * saveColLb = new double[numcols];
78    double * saveColUb = new double[numcols];
79    CoinCopyN(nlp_->getColLower(), numcols , saveColLb);
80    CoinCopyN(nlp_->getColUpper(), numcols , saveColUb);
81    for (int i = 0; i < numcols ; i++) {
82      if (nlp_->isInteger(i)) {
83        nlp_->setColBounds(i,si.getColLower()[i],si.getColUpper()[i]);
84      }
85    }
86
87#if 0
88//Add tight cuts at LP optimum
89    int numberCuts = si.getNumRows() - nlp_->getNumRows();
90    const OsiRowCut ** cuts = new const OsiRowCut*[numberCuts];
91    int begin = nlp_->getNumRows();
92    numberCuts = 0;
93    int end = si.getNumRows();
94    const double * rowLower = si.getRowLower();
95    const double * rowUpper = si.getRowUpper();
96    const CoinPackedMatrix * mat = si.getMatrixByRow();
97    const CoinBigIndex * starts = mat->getVectorStarts();
98    const int * lengths = mat->getVectorLengths();
99    const double * elements = mat->getElements();
100    const int * indices = mat->getIndices();
101
102    for (int i = begin ; i < end ; i++, numberCuts++) {
103      bool nnzExists=false;
104      for (int k = starts[i] ; k < starts[i]+lengths[i] ; k++) {
105        if (indices[k] == nlp_->getNumCols()) {
106          nnzExists = true;
107          char sign = (elements[k]>0.)?'+':'-';
108          char type='<';
109          if (rowLower[i]>-1e20) type='>';
110        }
111      }
112      if (nnzExists) {
113        numberCuts--;
114        continue;
115      }
116      int * indsCopy = CoinCopyOfArray(&indices[starts[i]], lengths[i]);
117      double * elemsCopy = CoinCopyOfArray(&elements[starts[i]], lengths[i]);
118      cuts[numberCuts] = new OsiRowCut(rowLower[i], rowUpper[i], lengths[i], lengths[i],
119          indsCopy, elemsCopy);
120    }
121    nlp_->applyRowCuts(numberCuts,cuts);
122#endif
123    //Now solve the NLP get the cuts, reset bounds and get out
124
125    //  nlp_->turnOnIpoptOutput();
126    nSolve_++;
127    nlp_->resolve();
128    const double * violatedPoint = (addOnlyViolated_)? si.getColSolution():
129        NULL;
130    nlp_->getOuterApproximation(cs, 1, violatedPoint,global_);
131    if (nlp_->isProvenOptimal()) {
132      handler_->message(LP_ERROR,messages_)
133      <<nlp_->getObjValue()-si.getObjValue()<<CoinMessageEol;
134      bool feasible = 1;
135      const double * colsol2 = nlp_->getColSolution();
136      for (int i = 0 ; i < numcols && feasible; i++) {
137        if (nlp_->isInteger(i)) {
138          if (fabs(colsol2[i] - floor(colsol2[i] + 0.5) ) > 1e-07)
139            feasible = 0;
140        }
141      }
142      if (feasible ) {
143#if 1
144        // Also store into solver
145        OsiAuxInfo * auxInfo = si.getAuxiliaryInfo();
146        OsiBabSolver * auxiliaryInfo = dynamic_cast<OsiBabSolver *> (auxInfo);
147        if (auxiliaryInfo) {
148          double * lpSolution = new double[numcols + 1];
149          CoinCopyN(colsol2, numcols, lpSolution);
150          lpSolution[numcols] = nlp_->getObjValue();
151          auxiliaryInfo->setSolution(lpSolution, numcols + 1, lpSolution[numcols]);
152          delete [] lpSolution;
153        }
154        else
155          printf("No auxiliary info in nlp solve!\n");
156#endif
157
158      }
159    }
160    else if (nlp_->isAbandoned() || nlp_->isIterationLimitReached()) {
161      throw CoinError("Unsolved NLP ... exit", "generateCuts", "OaNlpOptim");
162    }
163    else {
164      //       //Add an infeasibility local constraint
165      //       CoinPackedVector v;
166      //       double rhs = 1.;
167      //       for(int i = 0; i < numcols ; i++)
168      //        {
169      //          if(nlp_->isInteger(i) && (si.getColUpper()[i] - si.getColLower()[i] < 0.9))
170      //            {
171      //              double value = floor(colsol[i] + 0.5);
172      //              assert(fabs(colsol[i]-value)<1e-8 && value >= -1e-08 && value <= 1+ 1e08);
173      //              v.insert(i, -(2*value - 1));
174      //              rhs -= value;
175      //            }
176      //        }
177      //       OsiRowCut cut;
178      //       cut.setRow(v);
179      //       cut.setLb(rhs);
180      //       cut.setUb(1e300);
181      //       cut.setGloballyValid();
182      //       cs.insert(cut);
183    }
184    for (int i = 0; i < numcols ; i++) {
185      if (nlp_->isInteger(i)) {
186        nlp_->setColBounds(i,saveColLb[i],saveColUb[i]);
187      }
188    }
189#if 0
190    nlp_->deleteLastRows(numberCuts);
191#endif
192    delete [] saveColLb;
193    delete [] saveColUb;
194  }
195
196  void
197  OaNlpOptim::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)
198  {
199    roptions->SetRegisteringCategory("Nlp solve options in B-Hyb", RegisteredOptions::BonminCategory);
200    roptions->AddLowerBoundedIntegerOption("nlp_solve_frequency",
201        "Specify the frequency (in terms of nodes) at which NLP relaxations are solved in B-Hyb.",
202        0,10,
203        "A frequency of 0 amounts to to never solve the NLP relaxation.");
204    roptions->setOptionExtraInfo("nlp_solve_frequency",1);
205    roptions->AddLowerBoundedIntegerOption("nlp_solve_max_depth",
206        "Set maximum depth in the tree at which NLP relaxations are solved in B-Hyb.",
207        0,10,
208        "A depth of 0 amounts to to never solve the NLP relaxation.");
209    roptions->setOptionExtraInfo("nlp_solve_max_depth",1);
210    roptions->AddLowerBoundedNumberOption("nlp_solves_per_depth",
211        "Set average number of nodes in the tree at which NLP relaxations are solved in B-Hyb for each depth.",
212        0.,false,1e30);
213    roptions->setOptionExtraInfo("nlp_solves_per_depth",1);
214  }
215
216
217}/* End namespace Bonmin. */
Note: See TracBrowser for help on using the repository browser.