00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "SepaTMINLP2OsiLP.hpp"
00010 #include "BonTypes.hpp"
00011 #include "OsiSolverInterface.hpp"
00012 #include "BonTMINLP2TNLP.hpp"
00013 #include "CoinPackedMatrix.hpp"
00014
00015 #include <vector>
00016 #include <sstream>
00017 #include <climits>
00018
00019 using namespace Ipopt;
00020
00021 namespace Sepa {
00022
00023 void
00024 SepaTMINLP2OsiLP::extract(OsiSolverInterface *si,
00025 const double * x, bool getObj)
00026 {
00027 assert(IsValid(model_));
00028 int n;
00029 int m;
00030 int nnz_jac_g;
00031 int nnz_h_lag;
00032 TNLP::IndexStyleEnum index_style;
00033
00034 model_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
00035
00036
00037 model_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, value_());
00038
00039
00040 Bonmin::vector<double> g(m);
00041 model_->eval_g(n, x, 1, m, g());
00042
00043 Bonmin::vector<double> rowLow(m);
00044 Bonmin::vector<double> rowUp(m);
00045
00046
00047
00048 const double * rowLower = model_->g_l();
00049 const double * rowUpper = model_->g_u();
00050 const double * colLower = model_->x_l();
00051 const double * colUpper = model_->x_u();
00052
00053 double nlp_infty = si->getInfinity();
00054 double infty = DBL_MAX;
00055
00056 for(int i = 0 ; i < m ; i++) {
00057 if(const_types_[i] == Ipopt::TNLP::NON_LINEAR) {
00058 if(rowLower[i] > - nlp_infty){
00059 rowLow[i] = (rowLower[i] - g[i]) - 1e-07;
00060 }
00061 else
00062 rowLow[i] = - infty;
00063 if(rowUpper[i] < nlp_infty)
00064 rowUp[i] = (rowUpper[i] - g[i]) + 1e-07;
00065 else
00066 rowUp[i] = infty;
00067 }
00068 else {
00069 if(rowLower[i] > -nlp_infty){
00070 rowLow[i] = (rowLower[i]);
00071 }
00072 else
00073 rowLow[i] = - infty;
00074 if(rowUpper[i] < nlp_infty){
00075 rowUp[i] = (rowUpper[i]);
00076 }
00077 else
00078 rowUp[i] = infty;
00079 }
00080 }
00081
00082
00083
00084
00085
00086 for(int i = 0 ; i < nnz_jac_g ; i++) {
00087 if(const_types_[iRow_[i]] != TNLP::LINEAR){
00088 if(
00089 cleanNnz(value_[i],colLower[jCol_[i]], colUpper[jCol_[i]],
00090 rowLower[iRow_[i]], rowUpper[iRow_[i]],
00091 x[jCol_[i]],
00092 rowLow[iRow_[i]],
00093 rowUp[iRow_[i]], tiny_, very_tiny_)) {
00094 if(rowLow[iRow_[i]] > -infty)
00095 rowLow[iRow_[i]] += value_[i] * x[jCol_[i]];
00096 if(rowUp[iRow_[i]] < infty)
00097 rowUp[iRow_[i]] += value_[i] *x[jCol_[i]];
00098 }
00099 }
00100 }
00101 CoinPackedMatrix mat(true, iRow_(), jCol_(), value_(), nnz_jac_g);
00102 mat.setDimensions(m,n);
00103
00104 #if 0
00105 vector<double> act(m);
00106 mat.times(x, act());
00107 for(int j = 0 ; j < m ; j++){
00108 if(j==10 && fabs(x[0] - 4.73032) < 1e-4)
00109 assert(act[j] + 1e-5 > rowLow[j] && act[j] - 1e-5 < rowUp[j] + g[j]);
00110 }
00111 #endif
00112
00113 Bonmin::vector<double> obj(n);
00114 for(int i = 0 ; i < n; i++)
00115 obj[i] = 0.;
00116
00117
00118 si->loadProblem(mat, colLower, colUpper, obj(), rowLow(), rowUp());
00119 for(int i = 0 ; i < n ; i++) {
00120 if(model_->var_types()[i] == Bonmin::TMINLP::BINARY || model_->var_types()[i] == Bonmin::TMINLP::INTEGER )
00121 si->setInteger(i);
00122 }
00123 if(getObj) {
00124 if(model_->hasLinearObjective()){
00125 double zero;
00126 Bonmin::vector<double> x0(n,0.);
00127 model_->eval_f(n, x0(), 1, zero);
00128 si->setDblParam(OsiObjOffset, -zero);
00129
00130 model_->eval_grad_f(n, x, 1,obj());
00131 si->setObjective(obj());
00132 }
00133 else {
00134 throw -1;
00135 }
00136
00137 }
00138
00139 OsiCuts cs;
00140 get_oas(cs, x, 0, 1);
00141 si->applyCuts(cs);
00142
00143
00144 }
00145
00146 void
00147 SepaTMINLP2OsiLP::get_oas(OsiCuts &cs, const double *x, bool getObj, bool global) const {
00148
00149 int n,m, nnz_jac_g, nnz_h_lag;
00150 TNLP::IndexStyleEnum index_style;
00151 model_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
00152
00153 Bonmin::vector<double> g(m);
00154
00155 model_->eval_jac_g(n, x, 1, m, nnz_jac_g, NULL, NULL, value_());
00156 model_->eval_g(n,x,0,m,g());
00157
00158
00159 Bonmin::vector<double> lb(m + 1);
00160 Bonmin::vector<double> ub(m + 1);
00161
00162 Bonmin::vector<int> row2cutIdx(m,-1);
00163
00164 std::vector<int> cut2rowIdx;
00165
00166 int numCuts = 0;
00167
00168 const double * rowLower = model_->g_l();
00169 const double * rowUpper = model_->g_u();
00170 const double * colLower = model_->x_l();
00171 const double * colUpper = model_->x_u();
00172 double nlp_infty = infty_;
00173 double infty = DBL_MAX;
00174
00175 for(int rowIdx = 0; rowIdx < m ; rowIdx++) {
00176 if(const_types_[rowIdx] == TNLP::NON_LINEAR) {
00177 row2cutIdx[rowIdx] = numCuts;
00178 cut2rowIdx.push_back(rowIdx);
00179 if(rowLower[rowIdx] > - nlp_infty)
00180 lb[numCuts] = rowLower[rowIdx] - g[rowIdx];
00181 else
00182 lb[numCuts] = - infty;
00183 if(rowUpper[rowIdx] < nlp_infty)
00184 ub[numCuts] = rowUpper[rowIdx] - g[rowIdx];
00185 else
00186 ub[numCuts] = infty;
00187 numCuts++;
00188 }
00189 }
00190
00191 lb.resize(numCuts);
00192 ub.resize(numCuts);
00193 Bonmin::vector<CoinPackedVector> cuts(numCuts);
00194
00195
00196 for(int i = 0 ; i < nnz_jac_g ; i++) {
00197 const int &rowIdx = iRow_[i];
00198 const int & cutIdx = row2cutIdx[ rowIdx ];
00199 if(cutIdx != -1) {
00200 const int &colIdx = jCol_[i];
00201
00202 if(cleanNnz(value_[i],colLower[colIdx], colUpper[colIdx],
00203 rowLower[rowIdx], rowUpper[rowIdx],
00204 x[colIdx],
00205 lb[cutIdx],
00206 ub[cutIdx], tiny_, very_tiny_)) {
00207 if(fabs(value_[i]) > 1e15) {
00208 printf("Not generating cut because of big coefficient %g col %i x[%i] = %g\n", value_[i], colIdx, colIdx, x[colIdx]);
00209 return;
00210 }
00211 cuts[cutIdx].insert(colIdx,value_[i]);
00212 if(lb[cutIdx] > - infty)
00213 lb[cutIdx] += value_[i] * x[colIdx];
00214 if(ub[cutIdx] < infty)
00215 ub[cutIdx] += value_[i] * x[colIdx];
00216 }
00217 }
00218 }
00219
00220 for(int cutIdx = 0; cutIdx < numCuts; cutIdx++) {
00221 OsiRowCut newCut;
00222 if(global)
00223 newCut.setGloballyValidAsInteger(1);
00224
00225 const int* ids = model_->get_const_xtra_id();
00226
00227 int binary_id = (ids == NULL) ? -1 : ids[cut2rowIdx[cutIdx]];
00228 if(binary_id>0) {
00229
00230 if (lb[cutIdx] > -infty) {
00231 cuts[cutIdx].insert(binary_id, -lb[cutIdx]);
00232 newCut.setLb(0);
00233 newCut.setUb(ub[cutIdx]);
00234
00235 }
00236 if (ub[cutIdx] < infty) {
00237 cuts[cutIdx].insert(binary_id, -ub[cutIdx]);
00238 newCut.setLb(lb[cutIdx]);
00239 newCut.setUb(0);
00240
00241 }
00242 }
00243 else {
00244 newCut.setLb(lb[cutIdx]);
00245 newCut.setUb(ub[cutIdx]);
00246 }
00247
00248 newCut.setRow(cuts[cutIdx]);
00249 cs.insert(newCut);
00250 }
00251 printf("++++++++ I have generated %i cuts +++++++++\n", numCuts);
00252
00253 return;
00254
00255
00256 }
00257
00258 void
00259 SepaTMINLP2OsiLP::get_oa(int rowIdx, OsiCuts &cs, const double *x, bool getObj, bool global) const {
00260
00261 int n,m, nnz_jac_g, nnz_h_lag;
00262 TNLP::IndexStyleEnum index_style;
00263 model_->get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
00264
00265 double gi;
00266 model_->eval_gi(n,x,1, rowIdx,gi);
00267 Bonmin::vector<int> jCol(n);
00268 int nnz;
00269 model_->eval_grad_gi(n, x, 0, rowIdx, nnz, jCol(), NULL);
00270 model_->eval_grad_gi(n, x, 0, rowIdx, nnz, NULL, value_());
00271
00272
00273
00274 double lb;
00275 double ub;
00276
00277 const double * rowLower = model_->g_l();
00278 const double * rowUpper = model_->g_u();
00279 const double * colLower = model_->x_l();
00280 const double * colUpper = model_->x_u();
00281 double nlp_infty = infty_;
00282 double infty = DBL_MAX;
00283
00284 if (rowLower[rowIdx] > -nlp_infty)
00285 lb = rowLower[rowIdx] - gi;
00286 else
00287 lb = -infty;
00288 if (rowUpper[rowIdx] < nlp_infty)
00289 ub = rowUpper[rowIdx] - gi;
00290 else
00291 ub = infty;
00292
00293 CoinPackedVector cut;
00294
00295
00296 for(int i = 0 ; i < nnz ; i++) {
00297 if(index_style == Ipopt::TNLP::FORTRAN_STYLE) jCol[i]--;
00298 const int &colIdx = jCol[i];
00299
00300 if(cleanNnz(value_[i],colLower[colIdx], colUpper[colIdx],
00301 rowLower[rowIdx], rowUpper[rowIdx],
00302 x[colIdx],
00303 lb,
00304 ub, tiny_, very_tiny_)) {
00305 if(fabs(value_[i]) > 1e15) {
00306 printf("Not generating cut because of big coefficient %g col %i x[%i] = %g\n", value_[i], colIdx, colIdx, x[colIdx]);
00307 return;
00308 }
00309 cut.insert(colIdx,value_[i]);
00310 if(lb > - infty)
00311 lb += value_[i] * x[colIdx];
00312 if(ub < infty)
00313 ub += value_[i] * x[colIdx];
00314 }
00315 }
00316
00317 OsiRowCut newCut;
00318 if(global)
00319 newCut.setGloballyValidAsInteger(1);
00320
00321 const int* ids = model_->get_const_xtra_id();
00322
00323 int binary_id = (ids == NULL) ? -1 : ids[rowIdx];
00324 if(binary_id>0) {
00325
00326 if (lb > -infty) {
00327 cut.insert(binary_id, -lb);
00328 newCut.setLb(0);
00329 newCut.setUb(ub);
00330
00331 }
00332 if (ub < infty) {
00333 cut.insert(binary_id, -ub);
00334 newCut.setLb(lb);
00335 newCut.setUb(0);
00336 }
00337 }
00338 else {
00339 newCut.setLb(lb);
00340 newCut.setUb(ub);
00341 }
00342
00343 newCut.setRow(cut);
00344 cs.insert(newCut);
00345 return;
00346 }
00347
00348 void
00349 SepaTMINLP2OsiLP::get_refined_oa(OsiCuts & cs) const{
00350 if(num_approx_ <= 0)
00351 return;
00352 int n;
00353 int m;
00354 int nnz_jac_g;
00355 int nnz_h_lag;
00356 Ipopt::TNLP::IndexStyleEnum index_style;
00357
00358 model_->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00359
00360
00361
00362 const double * colLower = model_->x_l();
00363 const double * colUpper = model_->x_u();
00364 Bonmin::vector<Ipopt::TNLP::LinearityType> varTypes(n);
00365
00366 model_->get_variables_linearity(n, varTypes());
00367
00368
00369
00370 double * p = CoinCopyOfArray(colLower, n);
00371 double * pp = CoinCopyOfArray(colLower, n);
00372 double * up = CoinCopyOfArray(colUpper, n);
00373
00374 std::vector<double> step(n);
00375
00376
00377 for (int i = 0; i < n; i++) {
00378 if (p[i] < -1e4){
00379 p[i] = pp[i] = -1e4;
00380 }
00381 if (up[i] > 1e4){
00382 up[i] = 1e4;
00383 }
00384 }
00385
00386
00387
00388 for (int i = 0; i < n; i++) {
00389
00390 if (varTypes[i] == Ipopt::TNLP::LINEAR) {
00391 step[i] = 0;
00392 p[i] = pp[i] = up[i] = 0;
00393 }
00394 else
00395 step[i] = (up[i] - p[i]) / (num_approx_);
00396
00397 }
00398 get_oas(cs, p, 0, true);
00399 for (int j = 1; j <= num_approx_; j++) {
00400
00401 for (int i = 0; i < n; i++) {
00402 pp[i] += step[i];
00403 }
00404
00405 get_oas(cs, pp, 0, true);
00406
00407 }
00408
00409 get_oas(cs, up, 0, true);
00410
00411 delete [] p;
00412 delete [] pp;
00413 delete [] up;
00414 }
00415
00416 void
00417 SepaTMINLP2OsiLP::add_outer_description_function_values(OsiSolverInterface &si) {
00418 int n;
00419 int m;
00420 int nnz_jac_g;
00421 int nnz_h_lag;
00422 Ipopt::TNLP::IndexStyleEnum index_style;
00423
00424 model_->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
00425
00426
00427
00428 const double * colLower = model_->x_l();
00429 const double * colUpper = model_->x_u();
00430 Bonmin::vector<Ipopt::TNLP::LinearityType> varTypes(n);
00431
00432 model_->get_variables_linearity(n, varTypes());
00433 const Bonmin::TMINLP::VariableType* variableType = model_->var_types();
00434
00435 OsiCuts cs;
00436
00437 double * p = CoinCopyOfArray(colLower, n);
00438 double * pp = CoinCopyOfArray(colLower, n);
00439 double * up = CoinCopyOfArray(colUpper, n);
00440 double * low = CoinCopyOfArray(colLower, n);
00441
00442 Bonmin::vector<double> step(n);
00443 Bonmin::vector<int> nbG(m,2);
00444
00445
00446 for (int i = 0; i < n; i++) {
00447 if (low[i] < -1e4){
00448 low[i] = -1e4;
00449 }
00450 if (up[i] > 1e4){
00451 up[i] = 1e4;
00452 }
00453 }
00454
00455
00456 for (int i = 0; i < n; i++) {
00457
00458 if (varTypes[i] == Ipopt::TNLP::LINEAR) {
00459 step[i] = 0;
00460 low[i] = p[i] = pp[i] = up[i] = 0;
00461 }
00462 else
00463 step[i] = (up[i] - low[i]) / (num_approx_*20);
00464
00465 }
00466 get_oas(cs, low, 0, true);
00467 get_oas(cs, up, 0, true);
00468 Bonmin::vector<double> g_p(m);
00469 Bonmin::vector<double> g_pp(m);
00470 model_->eval_g(n, low, 1, m, g_p());
00471 model_->eval_g(n, up, 1, m, g_pp());
00472
00473 for (int i = 0; (i < m); i++) {
00474 if(const_types_[i] != Ipopt::TNLP::NON_LINEAR) continue;
00475 double thresh = std::abs(g_pp[i] - g_p[i])/(num_approx_-1);
00476
00477 printf("Constraint %i threshold is %g lower val %g upper %g\n",i, thresh,
00478 g_p[i], g_pp[i]);
00479 std::copy(low, low + n, p);
00480 std::copy(low, low + n, pp);
00481 double g_p_i, g_pp_i;
00482 g_p_i = g_p[i];
00483 int n_steps = 0;
00484 while (nbG[i] < num_approx_) {
00485 n_steps++;
00486
00487 for (int j = 0; j < n; j++) {
00488 pp[j] += step[j];
00489 }
00490 model_->eval_gi(n, pp, 1, i, g_pp_i);
00491
00492 if (std::abs(g_p_i - g_pp_i)>=thresh ) {
00493 printf("Function value after %i steps %g\n", n_steps, g_pp_i);
00494 get_oa(i, cs, pp, 0, true);
00495 for (int j = 0; j < n; j++) {
00496 p[j] = pp[j];
00497 }
00498 g_p_i = g_pp_i;
00499 nbG[i]++;
00500 }
00501 }
00502
00503 }
00504
00505
00506
00507 si.applyCuts(cs);
00508 delete [] p;
00509 delete [] pp;
00510 delete [] up;
00511 }
00512 }