/home/coin/SVN-release/OS-2.4.2/Bonmin/src/Algorithms/QuadCuts/BonOuterApprox.cpp

Go to the documentation of this file.
00001 // (C) Copyright International Business Machines Corporation 2007
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // Authors :
00006 // Pierre Bonami, International Business Machines Corporation
00007 //
00008 // Date : 10/16/2007
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     Ipopt::TNLP::IndexStyleEnum index_style;
00035 
00036     Bonmin::TMINLP2TNLP * model = minlp.problem();
00037     
00038     //Get problem information
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     //get Jacobian
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     //Put jacobian arrays in c style 
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] != Ipopt::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] != Ipopt::TNLP::LINEAR) rowUp[i] += 1e-07;
00079         }
00080     }
00081     
00082     //Then convert everything to a CoinPackedMatrix
00083     //Go through values, clean coefficients and fix bounds
00084     for(int i = 0 ; i < nnz_jac_g ; i++) {
00085       if(const_types[iRow[i]] != Ipopt::TNLP::LINEAR || //For linear just copy is fine.
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); // In case matrix was empty, this should be enough
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          //if(fabs(zero - 0) > 1e-10)
00113            //addObjVar = true;
00114          //else { 
00115            //Copy the linear objective and don't create a dummy variable.
00116            model->eval_grad_f(n, x, 1,&obj[0]);
00117            si->setObjective(&obj[0]);
00118          //}
00119       }
00120      else {
00121         //add variable alpha
00122         //(alpha should be empty in the matrix with a coefficient of -1 and unbounded)
00123         CoinPackedVector a;
00124         si->addCol(a,-si->getInfinity(), si->getInfinity(), 1.);
00125     
00126         // Now get the objective cuts
00127         // get the gradient, pack it and add the cut
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 //Unconstrained problem can not put clean coefficient
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 

Generated on Wed Nov 30 03:03:53 2011 by  doxygen 1.4.7