/home/coin/SVN-release/OS-2.0.0/Bonmin/src/Interfaces/BonBranchingTQP.cpp

Go to the documentation of this file.
00001 // (C) Copyright International Business Machines Corporation and Carnegie Mellon University 2006, 2008
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // Authors :
00006 // Andreas Waechter, International Business Machines Corporation
00007 //                   (derived from BonTMINLP2TNLP.cpp)            12/22/2006
00008 // Authors :
00009 
00010 
00011 #include "BonBranchingTQP.hpp"
00012 #include "IpBlas.hpp"
00013 #include "IpAlgTypes.hpp"
00014 #include <string>
00015 #include <fstream>
00016 #include <sstream>
00017 namespace Bonmin
00018 {
00019   BranchingTQP::BranchingTQP(SmartPtr<TMINLP2TNLP> tminlp2tnlp)
00020     :
00021     tminlp2tnlp_(tminlp2tnlp)
00022   {
00023     bool retval = tminlp2tnlp_->get_nlp_info(n_, m_, nnz_jac_g_,
00024                                              nnz_h_lag_, index_style_);
00025     ASSERT_EXCEPTION(retval, TMINLP_INVALID,
00026                      "Can't get NLP infor in BranchingTQP");
00027     //DBG_ASSERT(x_sol_);
00028     //DBG_ASSERT(duals_sol_);
00029 
00030     obj_grad_ = new Number[n_];
00031     obj_hess_ = new Number[nnz_h_lag_];
00032     obj_hess_irow_ = new Index[nnz_h_lag_];
00033     obj_hess_jcol_ = new Index[nnz_h_lag_];
00034     g_vals_ = new Number[m_];
00035     g_jac_ = new Number[nnz_jac_g_];
00036     g_jac_irow_ = new Index[nnz_jac_g_];
00037     g_jac_jcol_ = new Index[nnz_jac_g_];
00038 
00039     const Number* x_sol = tminlp2tnlp_->x_sol();
00040     const Number* duals_sol = tminlp2tnlp_->duals_sol();
00041 
00042     // Compute all nonlinear values at the starting point so that we
00043     // have all the information for the QP
00044     bool new_x = true;   // ToDo: maybe NOT new?
00045     retval = tminlp2tnlp_->eval_f(n_, x_sol, new_x, obj_val_);
00046     ASSERT_EXCEPTION(retval, TMINLP_INVALID,
00047                      "Can't evaluate objective function in BranchingTQP");
00048     new_x = false;
00049     retval = tminlp2tnlp_->eval_grad_f(n_, x_sol, new_x, obj_grad_);
00050     ASSERT_EXCEPTION(retval, TMINLP_INVALID,
00051                      "Can't evaluate objective gradient in BranchingTQP");
00052     bool new_lambda = true; // ToDo: maybe NOT new?
00053     retval = tminlp2tnlp_->eval_h(n_, x_sol, new_x, 1., m_, duals_sol + 2 * n_,
00054                              new_lambda, nnz_h_lag_, obj_hess_irow_,
00055                              obj_hess_jcol_, NULL);
00056     ASSERT_EXCEPTION(retval, TMINLP_INVALID,
00057                      "Can't evaluate objective Hessian structure in BranchingTQP");
00058     if (index_style_ == TNLP::FORTRAN_STYLE) {
00059       for (Index i=0; i<nnz_h_lag_; i++) {
00060         obj_hess_irow_[i]--;
00061         obj_hess_jcol_[i]--;
00062       }
00063     }
00064     retval = tminlp2tnlp_->eval_h(n_, x_sol, new_x, 1., m_, duals_sol + 2*n_,
00065                              new_lambda, nnz_h_lag_, NULL, NULL, obj_hess_);
00066     ASSERT_EXCEPTION(retval, TMINLP_INVALID,
00067                      "Can't evaluate objective Hessian values in BranchingTQP");
00068     retval = tminlp2tnlp_->eval_g(n_, x_sol, new_x, m_, g_vals_);
00069     ASSERT_EXCEPTION(retval, TMINLP_INVALID,
00070                      "Can't evaluate constraint values in BranchingTQP");
00071     retval = tminlp2tnlp_->eval_jac_g(n_, x_sol, new_x, m_, nnz_jac_g_,
00072                                  g_jac_irow_, g_jac_jcol_, NULL);
00073     ASSERT_EXCEPTION(retval, TMINLP_INVALID,
00074                      "Can't evaluate constraint Jacobian structure in BranchingTQP");
00075     if (index_style_ == TNLP::FORTRAN_STYLE) {
00076       for (Index i=0; i<nnz_jac_g_; i++) {
00077         g_jac_irow_[i]--;
00078         g_jac_jcol_[i]--;
00079       }
00080     }
00081     retval = tminlp2tnlp_->eval_jac_g(n_, x_sol, new_x, m_, nnz_jac_g_,
00082                                  NULL, NULL, g_jac_);
00083     ASSERT_EXCEPTION(retval, TMINLP_INVALID,
00084                      "Can't evaluate constraint Jacobian values in BranchingTQP");
00085 
00086     // Keep copy of original x_sol and duals_sol values
00087     x_sol_copy_ = new Number[n_];
00088     IpBlasDcopy(n_, x_sol, 1, x_sol_copy_, 1);
00089     duals_sol_copy_ = new Number[m_ + 2*n_];
00090     IpBlasDcopy(m_+2*n_, duals_sol, 1, duals_sol_copy_, 1);
00091 #if 0
00092     for (int i=0; i<n_; i++) {
00093       printf("x_sol_copy_[%3d] = %15.8e duals_sol_copy_[%3d] = %15.8e obj_grad_[%3d] = %15.8e\n", i,x_sol_copy_[i],i,duals_sol_copy_[i],i,obj_grad_[i]);
00094     }
00095     for (int i=0; i<m_; i++) {
00096       printf("g_vals_[%3d] = %15.8e\n", i, g_vals_[i]);
00097     }
00098     for (int i=0; i<nnz_h_lag_; i++) {
00099       printf("Hess[%3d,%3d] = %15.8e\n",obj_hess_irow_[i],obj_hess_jcol_[i],obj_hess_[i]);
00100     }
00101     for (int i=0; i<nnz_jac_g_; i++) {
00102       printf("Jac[%3d,%3d] = %15.8e\n",g_jac_irow_[i],g_jac_jcol_[i],g_jac_[i]);
00103     }
00104 #endif
00105   }
00106 
00107   BranchingTQP::~BranchingTQP()
00108   {
00109     delete [] obj_grad_;
00110     delete [] obj_hess_;
00111     delete [] obj_hess_irow_;
00112     delete [] obj_hess_jcol_;
00113     delete [] g_vals_;
00114     delete [] g_jac_;
00115     delete [] g_jac_irow_;
00116     delete [] g_jac_jcol_;
00117     delete [] x_sol_copy_;
00118     delete [] duals_sol_copy_;
00119   }
00120 
00121   bool BranchingTQP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00122                                   Index& nnz_h_lag,
00123                                   IndexStyleEnum& index_style)
00124   {
00125     n = n_;
00126     m = m_;
00127     nnz_jac_g = nnz_jac_g_;
00128     nnz_h_lag = nnz_h_lag_;
00129     index_style = index_style_;
00130     return true;
00131   }
00132 
00133   bool BranchingTQP::get_bounds_info(Index n, Number* x_l, Number* x_u,
00134                                      Index m, Number* g_l, Number* g_u)
00135   {
00136     DBG_ASSERT(n == n_);
00137     DBG_ASSERT(m == m_);
00138     bool retval = tminlp2tnlp_->get_bounds_info(n, x_l, x_u, m, g_l, g_u);
00139     // correct for displacement
00140     for (int i=0; i<n; i++) {
00141       x_l[i] -= x_sol_copy_[i];
00142       x_u[i] -= x_sol_copy_[i];
00143     }
00144     // include the right hand side of the constraint
00145     for (int i=0; i<m; i++) {
00146       g_l[i] -= g_vals_[i];
00147       g_u[i] -= g_vals_[i];
00148     }
00149     return retval;
00150   }
00151 
00152   bool BranchingTQP::get_starting_point(Index n, bool init_x, Number* x,
00153       bool init_z, Number* z_L, Number* z_U,
00154       Index m, bool init_lambda,
00155       Number* lambda)
00156   {
00157     DBG_ASSERT(n==n_);
00158     if (init_x == true) {
00159       const double zero = 0.;
00160       IpBlasDcopy(n, &zero, 0, x, 1);
00161     }
00162     if (init_z == true) {
00163       DBG_ASSERT(duals_sol_copy_);
00164       IpBlasDcopy(n, duals_sol_copy_, 1, z_L, 1);
00165       IpBlasDcopy(n, duals_sol_copy_ + n, 1, z_U, 1);
00166 
00167     }
00168     if(init_lambda == true) {
00169       DBG_ASSERT(duals_sol_copy_);
00170       IpBlasDcopy(m_, duals_sol_copy_ + 2*n_, 1, lambda, 1);
00171       for(int i = m_ ; i < m; i++)
00172       {
00173         lambda[i] = 0.;
00174       }
00175     }
00176 
00177     return true;
00178   }
00179 
00180   bool
00181   BranchingTQP::get_constraints_linearity(Index m, LinearityType* const_types)
00182   {
00183     for (int i=0; i<m; i++) {
00184       const_types[i] = LINEAR;
00185     }
00186     return true;
00187   }
00188 
00189   bool BranchingTQP::eval_f(Index n, const Number* x, bool new_x,
00190       Number& obj_value)
00191   {
00192     DBG_ASSERT(n == n_);
00193 
00194     obj_value = IpBlasDdot(n, x, 1, obj_grad_, 1);
00195     for (int i=0; i<nnz_h_lag_; i++) {
00196       Index& irow = obj_hess_irow_[i];
00197       Index& jcol = obj_hess_jcol_[i];
00198       if (irow!=jcol) {
00199         obj_value += obj_hess_[i]*x[irow]*x[jcol];
00200       }
00201       else {
00202         obj_value += 0.5*obj_hess_[i]*x[irow]*x[irow];
00203       }
00204     }
00205 
00206     return true;
00207   }
00208 
00209   bool BranchingTQP::eval_grad_f(Index n, const Number* x, bool new_x,
00210                                  Number* grad_f)
00211   {
00212     DBG_ASSERT(n == n_);
00213 
00214     IpBlasDcopy(n_, obj_grad_, 1, grad_f, 1);
00215     for (int i=0; i<nnz_h_lag_; i++) {
00216       Index& irow = obj_hess_irow_[i];
00217       Index& jcol = obj_hess_jcol_[i];
00218       grad_f[irow] += obj_hess_[i]*x[jcol];
00219       if (irow!=jcol) {
00220         grad_f[jcol] += obj_hess_[i]*x[irow];
00221       }
00222     }
00223 
00224     return true;
00225   }
00226 
00227   bool BranchingTQP::eval_g(Index n, const Number* x, bool new_x,
00228                             Index m, Number* g)
00229   {
00230     DBG_ASSERT(n == n_);
00231 
00232     const double zero = 0.;
00233     IpBlasDcopy(m_, &zero, 0, g, 1);
00234     for (Index i=0; i<nnz_jac_g_; i++) {
00235       Index& irow = g_jac_irow_[i];
00236       Index& jcol = g_jac_jcol_[i];
00237       g[irow] += g_jac_[i]*x[jcol];
00238     }
00239 
00240     return true;
00241   }
00242 
00243   bool BranchingTQP::eval_jac_g(Index n, const Number* x, bool new_x,
00244                                 Index m, Index nele_jac, Index* iRow,
00245                                 Index *jCol, Number* values)
00246   {
00247     if (iRow != NULL) {
00248       DBG_ASSERT(jCol != NULL);
00249       DBG_ASSERT(values == NULL);
00250       if (index_style_ == TNLP::FORTRAN_STYLE) {
00251         for (Index i=0; i<nnz_jac_g_; i++) {
00252           iRow[i] = g_jac_irow_[i] + 1;
00253           jCol[i] = g_jac_jcol_[i] + 1;
00254         }
00255       }
00256       else {
00257         for (Index i=0; i<nnz_jac_g_; i++) {
00258           iRow[i] = g_jac_irow_[i];
00259           jCol[i] = g_jac_jcol_[i];
00260         }
00261       }
00262     }
00263     else {
00264       IpBlasDcopy(nnz_jac_g_, g_jac_, 1, values, 1);
00265     }
00266 
00267     return true;
00268   }
00269 
00270   bool BranchingTQP::eval_h(Index n, const Number* x, bool new_x,
00271       Number obj_factor, Index m, const Number* lambda,
00272       bool new_lambda, Index nele_hess,
00273       Index* iRow, Index* jCol, Number* values)
00274   {
00275     DBG_ASSERT(nele_hess == nnz_h_lag_);
00276 
00277     if (iRow != NULL) {
00278       DBG_ASSERT(jCol != NULL);
00279       DBG_ASSERT(values == NULL);
00280       if (index_style_ == TNLP::FORTRAN_STYLE) {
00281         for (Index i=0; i<nele_hess; i++) {
00282           iRow[i] = obj_hess_irow_[i] + 1;
00283           jCol[i] = obj_hess_jcol_[i] + 1;
00284         }
00285       }
00286       else {
00287         for (Index i=0; i<nele_hess; i++) {
00288           iRow[i] = obj_hess_irow_[i];
00289           jCol[i] = obj_hess_jcol_[i];
00290         }
00291       }
00292     }
00293     else {
00294       if (obj_factor==0.) {
00295         const Number zero = 0.;
00296         IpBlasDcopy(nele_hess, &zero, 0, values, 1);
00297       }
00298       else {
00299         IpBlasDcopy(nele_hess, obj_hess_, 1, values, 1);
00300         if (obj_factor != 1.) {
00301           IpBlasDscal(nele_hess, obj_factor, values, 1);
00302         }
00303       }
00304     }
00305 
00306     return true;
00307   }
00308 
00309   void BranchingTQP::finalize_solution(SolverReturn status,
00310                                        Index n, const Number* x,
00311                                        const Number* z_L, const Number* z_U,
00312                                        Index m, const Number* g,
00313                                        const Number* lambda, Number obj_value,
00314                                        const IpoptData* ip_data,
00315                                        IpoptCalculatedQuantities* ip_cq)
00316   {
00317     // Translate displacement solution back to a solution in the real
00318     // variables
00319     double* xx = new double[n];
00320     for (int i=0; i<n; i++) {
00321       xx[i] = x_sol_copy_[i] + x[i];
00322     }
00323 
00324     tminlp2tnlp_->finalize_solution(status, n, xx, z_L, z_U, m, g, lambda,
00325                                     obj_value+obj_val_, ip_data, ip_cq);
00326     delete [] xx;
00327   }
00328 
00329 }
00330 // namespace Ipopt
00331 

Generated on Mon Aug 3 03:02:18 2009 by  doxygen 1.4.7