source: branches/devel/Bonmin/src/OaInterface/IpCbcOACutGenerator.cpp @ 44

Last change on this file since 44 was 44, checked in by pbonami, 13 years ago

Can add linear cut to nonlinear formulation in Ipopt

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 5.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 "IpCbcOACutGenerator.hpp"
11#include "OsiAuxInfo.hpp"
12/// Default constructor
13IpCbcOACutGenerator::IpCbcOACutGenerator(IpoptInterface * si,
14    int maxDepth)
15    :
16    CglCutGenerator(),
17    nlp_(si),
18    maxDepth_(maxDepth),
19    nSolve_(0)
20{
21  handler_ = new CoinMessageHandler();
22  handler_ -> setLogLevel(1);
23  messages_ = OaMessages();
24}
25/// Assign an IpoptInterface
26void
27IpCbcOACutGenerator::assignInterface(IpoptInterface * si)
28
29{
30  nlp_ = si;
31}
32static int nCalls = 0;
33/// cut generation method
34void
35IpCbcOACutGenerator::generateCuts( const OsiSolverInterface & si, OsiCuts & cs,
36    const CglTreeInfo info) const
37{
38  if(nlp_ == NULL) {
39    std::cerr<<"Error in cut generator for outer approximation no ipopt NLP assigned"<<std::endl;
40    throw -1;
41  }
42
43  int numcols = nlp_->getNumCols();
44
45  //Get the continuous solution
46  //const double *colsol = si.getColSolution();
47  //Check for integer feasibility
48  nCalls++;
49  if(info.level > maxDepth_ || info.pass > 0 || !info.inTree  )
50    return;
51  //Fix the variable which have to be fixed, after having saved the bounds
52  double * saveColLb = new double[numcols];
53  double * saveColUb = new double[numcols];
54  CoinCopyN(nlp_->getColLower(), numcols , saveColLb);
55  CoinCopyN(nlp_->getColUpper(), numcols , saveColUb);
56  for(int i = 0; i < numcols ; i++) {
57    if(nlp_->isInteger(i)) {
58      nlp_->setColBounds(i,si.getColLower()[i],si.getColUpper()[i]);
59    }
60  }
61
62#ifdef 0
63 //Add tight cuts at LP optimum
64 int numberCuts = si.getNumRows() - nlp_->getNumRows();
65 const OsiRowCut ** cuts = new const OsiRowCut*[numberCuts];
66 int begin = nlp_->getNumRows();
67 numberCuts = 0;
68 int end = si.getNumRows();
69 const double * rowLower = si.getRowLower();
70 const double * rowUpper = si.getRowUpper();
71 const CoinPackedMatrix * mat = si.getMatrixByRow();     
72 const CoinBigIndex * starts = mat->getVectorStarts();
73 const int * lengths = mat->getVectorLengths();
74 const double * elements = mat->getElements();
75 const int * indices = mat->getIndices();
76 
77 for(int i = begin ; i < end ; i++, numberCuts++)
78 {
79   bool nnzExists=false;
80   for(int k = starts[i] ; k < starts[i]+lengths[i] ; k++)
81   {
82     if(indices[k] == nlp_->getNumCols())
83     {
84       nnzExists = true;
85       char sign = (elements[k]>0.)?'+':'-';
86       char type='<';
87       if(rowLower[i]>-1e20) type='>'; 
88       std::cout<<"Non zero with sign: "<<sign<<", type: "<<type<<std::endl;
89     }
90   }
91   if(nnzExists)
92   {
93      numberCuts--;
94      continue;
95   }
96   else
97      std::cout<<"No nonzero element"<<std::endl;
98   int * indsCopy = CoinCopyOfArray(&indices[starts[i]], lengths[i]);
99   double * elemsCopy = CoinCopyOfArray(&elements[starts[i]], lengths[i]);
100   cuts[numberCuts] = new OsiRowCut(rowLower[i], rowUpper[i], lengths[i], lengths[i],
101                           indsCopy, elemsCopy);
102 }
103 nlp_->applyRowCuts(numberCuts,cuts);
104#endif
105  //Now solve the NLP get the cuts, reset bounds and get out
106
107  //  nlp_->turnOnIpoptOutput();
108  nSolve_++;
109  nlp_->resolve();
110  nlp_->getOuterApproximation(cs);
111  if(nlp_->isProvenOptimal()) {
112    handler_->message(LP_ERROR,messages_)
113    <<nlp_->getObjValue()-si.getObjValue()<<CoinMessageEol;
114    bool feasible = 1;
115    const double * colsol2 = nlp_->getColSolution();
116    for(int i = 0 ; i < numcols && feasible; i++) {
117      if(nlp_->isInteger(i)) {
118        if(fabs(colsol2[i] - floor(colsol2[i] + 0.5) ) > 1e-07)
119          feasible = 0;
120      }
121    }
122    if(feasible ) {
123#if 1
124      // Also store into solver
125      OsiAuxInfo * auxInfo = si.getAuxiliaryInfo();
126      OsiBabSolver * auxiliaryInfo = dynamic_cast<OsiBabSolver *> (auxInfo);
127      if (auxiliaryInfo) {
128        double * lpSolution = new double[numcols + 1];
129        CoinCopyN(colsol2, numcols, lpSolution);
130        lpSolution[numcols] = nlp_->getObjValue();
131        auxiliaryInfo->setSolution(lpSolution, numcols + 1, lpSolution[numcols]);
132        delete [] lpSolution;
133      }
134      else
135        printf("No auxiliary info in nlp solve!\n");
136#endif
137
138    }
139  }
140  else if(nlp_->isAbandoned() || nlp_->isIterationLimitReached()) {
141    std::cerr<<"Unsolved NLP ... exit"<<std::endl;
142    nlp_->turnOnIpoptOutput();
143    nlp_->resolve();
144    throw -1;
145  }
146  else {
147    //   std::cout<<"Infeasible NLP => Infeasible node"<<std::endl;
148    //       //Add an infeasibility local constraint
149    //       CoinPackedVector v;
150    //       double rhs = 1.;
151    //       for(int i = 0; i < numcols ; i++)
152    //  {
153    //    if(nlp_->isInteger(i) && (si.getColUpper()[i] - si.getColLower()[i] < 0.9))
154    //      {
155    //        double value = floor(colsol[i] + 0.5);
156    //        assert(fabs(colsol[i]-value)<1e-8 && value >= -1e-08 && value <= 1+ 1e08);
157    //        v.insert(i, -(2*value - 1));
158    //        rhs -= value;
159    //      }
160    //  }
161    //       OsiRowCut cut;
162    //       cut.setRow(v);
163    //       cut.setLb(rhs);
164    //       cut.setUb(1e300);
165    //       cut.setGloballyValid();
166    //       cs.insert(cut);
167  }
168  for(int i = 0; i < numcols ; i++) {
169    if(nlp_->isInteger(i)) {
170      nlp_->setColBounds(i,saveColLb[i],saveColUb[i]);
171    }
172  }
173  nlp_->deleteLastRows(numberCuts);
174  delete [] saveColLb;
175  delete [] saveColUb;
176}
Note: See TracBrowser for help on using the repository browser.