00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "BonOuterApprox.hpp"
00010 #include "BonBabSetupBase.hpp"
00011 #include "BonOsiTMINLPInterface.hpp"
00012 #include "BonTypes.hpp"
00013 #include <vector>
00014 #include <sstream>
00015 namespace Bonmin {
00016
00017 int OuterApprox::nTimesCalled = 0;
00018
00019 void
00020 OuterApprox::initialize(Bonmin::BabSetupBase & b){
00021 b.options()->GetNumericValue("tiny_element", tiny_, "bonmin.");
00022 b.options()->GetNumericValue("very_tiny_element", veryTiny_, "bonmin.");
00023 }
00024 void
00025 OuterApprox::extractLinearRelaxation(Bonmin::OsiTMINLPInterface &minlp,
00026 OsiSolverInterface *si,
00027 const double * x, bool getObj)
00028 {
00029
00030 int n;
00031 int m;
00032 int nnz_jac_g;
00033 int nnz_h_lag;
00034 TNLP::IndexStyleEnum index_style;
00035
00036 Bonmin::TMINLP2TNLP * model = minlp.problem();
00037
00038
00039 model->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
00040
00041 vector<int> iRow(nnz_jac_g);
00042 vector<int> jCol(nnz_jac_g);
00043 vector<double> vals(nnz_jac_g);
00044
00045
00046 model->eval_jac_g(n, x, 1, m, nnz_jac_g, &iRow[0], &jCol[0], NULL);
00047 model->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, &vals[0]);
00048
00049
00050 if(index_style == Ipopt::TNLP::FORTRAN_STYLE){
00051 for(int k = 0 ; k < nnz_jac_g ; k++){
00052 iRow[k]--;
00053 jCol[k]--;
00054 }
00055 }
00056 vector<double> g(m);
00057 model->eval_g(n, x, 1, m, &g[0]);
00058
00059 vector<double> rowLow(m);
00060 vector<double> rowUp(m);
00061 vector<double> colUp(n);
00062 vector<double> colLow(n);
00063
00064 model->get_bounds_info(n, &colLow[0], &colUp[0], m, &rowLow[0], &rowUp[0]);
00065
00066 double infty = si->getInfinity();
00067 double nlp_infty = infty;
00068
00069 vector<Ipopt::TNLP::LinearityType> const_types(m);
00070 model->get_constraints_linearity(m, &const_types[0]);
00071 for(int i = 0 ; i < m ; i++) {
00072 if(rowLow[i] > - nlp_infty){
00073 rowLow[i] = (rowLow[i] - g[i]);
00074 if(1 || const_types[i] != TNLP::LINEAR) rowLow[i] -= 1e-07;
00075 }
00076 if(rowUp[i] < nlp_infty){
00077 rowUp[i] = (rowUp[i] - g[i]);
00078 if(1 || const_types[i] != TNLP::LINEAR) rowUp[i] += 1e-07;
00079 }
00080 }
00081
00082
00083
00084 for(int i = 0 ; i < nnz_jac_g ; i++) {
00085 if(const_types[iRow[i]] != TNLP::LINEAR ||
00086 cleanNnz(vals[i],colLow[jCol[i]], colUp[jCol[i]],
00087 rowLow[iRow[i]], rowUp[iRow[i]],
00088 x[jCol[i]],
00089 rowLow[iRow[i]],
00090 rowUp[iRow[i]], tiny_, veryTiny_)) {
00091 rowLow[iRow[i]] += vals[i] * x[jCol[i]];
00092 rowUp[iRow[i]] += vals[i] *x[jCol[i]];
00093 }
00094 }
00095 CoinPackedMatrix mat(true, &iRow[0], &jCol[0], &vals[0], nnz_jac_g);
00096 mat.setDimensions(m,n);
00097
00098 vector<double> obj(n,0.);
00099
00100 si->loadProblem(mat, &colLow[0], &colUp[0], &obj[0], &rowLow[0], &rowUp[0]);
00101
00102 for(int i = 0 ; i < n ; i++) {
00103 if(minlp.isInteger(i))
00104 si->setInteger(i);
00105 }
00106 if(getObj) {
00107 if(model->hasLinearObjective()){
00108 std::cout<<"Linear stuff"<<std::endl;
00109 double zero;
00110 model->eval_f(n, &obj[0], 1, zero);
00111 si->setDblParam(OsiObjOffset, -zero);
00112
00113
00114
00115
00116 model->eval_grad_f(n, x, 1,&obj[0]);
00117 si->setObjective(&obj[0]);
00118
00119 }
00120 else {
00121
00122
00123 CoinPackedVector a;
00124 si->addCol(a,-si->getInfinity(), si->getInfinity(), 1.);
00125
00126
00127
00128 model->eval_grad_f(n, x, 1,&obj[0]);
00129 double ub;
00130 model->eval_f(n, x, 1, ub);
00131 ub*=-1;
00132 double lb = -1e300;
00133 CoinPackedVector objCut;
00134 CoinPackedVector * v = &objCut;
00135 v->reserve(n);
00136 for(int i = 0; i<n ; i++) {
00137 if(nnz_jac_g)
00138 {
00139 if(cleanNnz(obj[i],colLow[i], colUp[i],
00140 -si->getInfinity(), 0,
00141 x[i],
00142 lb,
00143 ub, tiny_, veryTiny_)) {
00144 v->insert(i,obj[i]);
00145 lb += obj[i] * x[i];
00146 ub += obj[i] * x[i];
00147 }
00148 }
00149 else
00150 {
00151 if(cleanNnz(obj[i],colLow[i], colUp[i],
00152 -si->getInfinity(), 0,
00153 x[i],
00154 lb,
00155 ub, 1e-03, 1e-08)) {
00156 v->insert(i,obj[i]);
00157 lb += obj[i] * x[i];
00158 ub += obj[i] * x[i];
00159 }
00160 }
00161 }
00162 v->insert(n,-1);
00163 si->addRow(objCut, lb, ub);
00164 }
00165 }
00166
00167 #if 0
00168 std::ostringstream os;
00169 os<<"OA_"<<nTimesCalled;
00170 nTimesCalled++;
00171 std::string f_name = os.str();
00172 si->writeMps(f_name.c_str());
00173 #endif
00174 }
00175
00176 }
00177