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

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