1 // 2 // IpoptProgram.h 3 // 4 5 #ifndef __Gravity____IpoptProgram__ 6 #define __Gravity____IpoptProgram__ 7 8 #include <stdio.h> 9 #include <assert.h> 10 #ifdef USE_IPOPT 11 #include <IpIpoptApplication.hpp> 12 #include <IpTNLP.hpp> 13 #endif 14 #include <gravity/model.h> 15 16 using Ipopt::IpoptApplication; 17 18 using namespace Ipopt; 19 using namespace gravity; 20 21 template<typename type = double> 22 class IpoptProgram : public TNLP, public Program<type>{ 23 24 25 public: 26 Model<type>* _model = nullptr; 27 IpoptProgram(Model<type> * m)28 IpoptProgram(Model<type>* m):_model(m){ 29 } 30 31 32 33 update_model()34 void update_model(){ 35 _model->fill_in_maps(); 36 } 37 38 39 40 41 /** Method to return some info about the nlp */ get_nlp_info(Index & n,Index & m,Index & nnz_jac_g,Index & nnz_h_lag,Ipopt::TNLP::IndexStyleEnum & index_style)42 bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, 43 Index& nnz_h_lag, Ipopt::TNLP::IndexStyleEnum& index_style){ 44 index_style = Ipopt::TNLP::C_STYLE; 45 n = (Index)_model->get_nb_vars(); 46 m = (Index)_model->get_nb_cons(); 47 nnz_jac_g = (Index)_model->get_nb_nnz_g(); 48 _model->_jac_vals.resize(nnz_jac_g,0); 49 nnz_h_lag = (Index)_model->get_nb_nnz_h(); 50 _model->_first_call_jac = true; 51 //If quadratic model and we're resolving no need to reset these 52 _model->_first_call_hess = true; 53 _model->_first_call_grad_obj = true; 54 return true; 55 } 56 finalize_solution(Ipopt::SolverReturn status,Index n,const Number * x,const Number * z_L,const Number * z_U,Index m,const Number * g,const Number * lambda,Number obj_value,const Ipopt::IpoptData * ip_data,Ipopt::IpoptCalculatedQuantities * ip_cq)57 void finalize_solution(Ipopt::SolverReturn status , 58 Index n , 59 const Number* x , 60 const Number* z_L , 61 const Number* z_U , 62 Index m , 63 const Number* g , 64 const Number* lambda , 65 Number obj_value , 66 const Ipopt::IpoptData* ip_data , 67 Ipopt::IpoptCalculatedQuantities* ip_cq 68 ) 69 { 70 _model->set_x(x); 71 for (auto &cp: _model->_cons_name) { 72 auto nb_ins = cp.second->get_nb_inst(); 73 cp.second->_dual.resize(nb_ins); 74 auto idx = 0; 75 for (size_t inst = 0; inst < nb_ins; inst++) { 76 if (!*cp.second->_all_lazy || !cp.second->_lazy[inst]) { 77 cp.second->_dual[inst] = lambda[cp.second->_id + idx++]; 78 } 79 } 80 } 81 for (auto &vp: _model->_vars) { 82 auto nb_inst = vp.second->get_dim(); 83 vp.second->_u_dual.resize(nb_inst); 84 vp.second->_l_dual.resize(nb_inst); 85 auto vid = vp.second->get_id(); 86 for (size_t inst = 0; inst < nb_inst; inst++) { 87 vp.second->_u_dual[inst] = z_U[vid + inst]; 88 vp.second->_l_dual[inst] = z_L[vid + inst]; 89 } 90 } 91 } 92 93 /** Method to return the bounds for my problem */ get_bounds_info(Index n,Number * x_l,Number * x_u,Index m,Number * g_l,Number * g_u)94 bool get_bounds_info(Index n, Number* x_l, Number* x_u, 95 Index m, Number* g_l, Number* g_u){ 96 assert(n==_model->get_nb_vars()); 97 assert(m==_model->get_nb_cons()); 98 _model->fill_in_var_bounds(x_l , x_u); 99 _model->fill_in_cstr_bounds(g_l , g_u); 100 101 return true; 102 } 103 get_starting_point(Index n,bool init_x,Number * x,bool init_z,Number * z_L,Number * z_U,Index m,bool init_lambda,Number * lambda)104 bool get_starting_point(Index n, bool init_x, Number* x, 105 bool init_z, Number* z_L, Number* z_U, 106 Index m, bool init_lambda, 107 Number* lambda){ 108 return get_starting_point_(n, init_x, x,init_z, z_L, z_U,m, init_lambda,lambda); 109 } 110 111 /** Method to return the starting point for the algorithm */ 112 template<typename T=type,typename std::enable_if<is_arithmetic<T>::value>::type* = nullptr> get_starting_point_(Index n,bool init_x,Number * x,bool init_z,Number * z_L,Number * z_U,Index m,bool init_lambda,Number * lambda)113 bool get_starting_point_(Index n, bool init_x, Number* x, 114 bool init_z, Number* z_L, Number* z_U, 115 Index m, bool init_lambda, 116 Number* lambda){ 117 assert(n==_model->get_nb_vars()); 118 assert(m==_model->get_nb_cons()); 119 120 if (init_x) { 121 _model->fill_in_var_init(x); 122 } 123 if (init_lambda && init_z) { 124 _model->fill_in_duals(lambda,z_L,z_U); 125 } 126 return true; 127 } 128 eval_f(Index n,const Number * x,bool new_x,Number & obj_value)129 bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value){ 130 return eval_f_(n, x, new_x, obj_value); 131 } 132 133 /** Method to return the objective value */ 134 template<typename T=type,typename std::enable_if<is_arithmetic<T>::value>::type* = nullptr> eval_f_(Index n,const Number * x,bool new_x,Number & obj_value)135 bool eval_f_(Index n, const Number* x, bool new_x, Number& obj_value){ 136 137 assert(n==_model->get_nb_vars()); 138 _model->fill_in_obj(x, obj_value,new_x); 139 return true; 140 } 141 eval_grad_f(Index n,const Number * x,bool new_x,Number * grad_f)142 bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f){ 143 return eval_grad_f_(n, x, new_x, grad_f); 144 } 145 146 /** Method to return the gradient of the objective */ 147 template<typename T=type,typename std::enable_if<is_arithmetic<T>::value>::type* = nullptr> eval_grad_f_(Index n,const Number * x,bool new_x,Number * grad_f)148 bool eval_grad_f_(Index n, const Number* x, bool new_x, Number* grad_f){ 149 150 assert(n==_model->get_nb_vars()); 151 _model->fill_in_grad_obj(x, grad_f, new_x); 152 return true; 153 } 154 eval_g(Index n,const Number * x,bool new_x,Index m,Number * g)155 bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g){ 156 return eval_g_(n, x, new_x, m, g); 157 } 158 159 /** Method to return the constraint residuals */ 160 template<typename T=type,typename std::enable_if<is_arithmetic<T>::value>::type* = nullptr> eval_g_(Index n,const Number * x,bool new_x,Index m,Number * g)161 bool eval_g_(Index n, const Number* x, bool new_x, Index m, Number* g){ 162 163 assert(n==_model->get_nb_vars()); 164 _model->fill_in_cstr(x, g, new_x); 165 return true; 166 } 167 eval_jac_g(Index n,const Number * x,bool new_x,Index m,Index nele_jac,Index * iRow,Index * jCol,Number * values)168 bool eval_jac_g(Index n, const Number* x, bool new_x, 169 Index m, Index nele_jac, Index* iRow, Index *jCol, 170 Number* values){ 171 return eval_jac_g_(n, x, new_x, m, nele_jac, iRow, jCol, values); 172 } 173 174 /** Method to return: 175 * 1) The structure of the jacobian (if "values" is NULL) 176 * 2) The values of the jacobian (if "values" is not NULL) 177 */ 178 template<typename T=type,typename std::enable_if<is_arithmetic<T>::value>::type* = nullptr> eval_jac_g_(Index n,const Number * x,bool new_x,Index m,Index nele_jac,Index * iRow,Index * jCol,Number * values)179 bool eval_jac_g_(Index n, const Number* x, bool new_x, 180 Index m, Index nele_jac, Index* iRow, Index *jCol, 181 Number* values){ 182 183 assert(n==_model->get_nb_vars()); 184 assert(m==_model->get_nb_cons()); 185 assert(nele_jac==_model->get_nb_nnz_g()); 186 if (values == NULL){ 187 _model->fill_in_jac_nnz(iRow, jCol); 188 } else { 189 _model->fill_in_jac(x, values, new_x); 190 } 191 192 return true; 193 } 194 eval_h(Index n,const Number * x,bool new_x,Number obj_factor,Index m,const Number * lambda,bool new_lambda,Index nele_hess,Index * iRow,Index * jCol,Number * values)195 bool eval_h(Index n, const Number* x, bool new_x, 196 Number obj_factor, Index m, const Number* lambda, 197 bool new_lambda, Index nele_hess, Index* iRow, 198 Index* jCol, Number* values){ 199 return eval_h_(n, x, new_x, obj_factor, m, lambda, new_lambda, nele_hess, iRow, jCol, values); 200 } 201 /** Method to return: 202 * 1) The structure of the hessian of the lagrangian (if "values" is NULL) 203 * 2) The values of the hessian of the lagrangian (if "values" is not NULL) 204 */ 205 template<typename T=type,typename std::enable_if<is_arithmetic<T>::value>::type* = nullptr> eval_h_(Index n,const Number * x,bool new_x,Number obj_factor,Index m,const Number * lambda,bool new_lambda,Index nele_hess,Index * iRow,Index * jCol,Number * values)206 bool eval_h_(Index n, const Number* x, bool new_x, 207 Number obj_factor, Index m, const Number* lambda, 208 bool new_lambda, Index nele_hess, Index* iRow, 209 Index* jCol, Number* values){ 210 211 212 assert(n==_model->get_nb_vars()); 213 assert(m==_model->get_nb_cons()); 214 assert(nele_hess==_model->get_nb_nnz_h()); 215 if (values == NULL){ 216 _model->fill_in_hess_nnz(iRow, jCol); 217 } else { 218 _model->fill_in_hess(x, obj_factor, lambda, values, new_x); 219 } 220 221 return true; 222 } 223 get_variables_linearity(Index n,LinearityType * var_types)224 bool get_variables_linearity(Index n, LinearityType* var_types){ 225 assert(n==_model->get_nb_vars()); 226 return false; 227 } 228 get_constraints_linearity(Index m,LinearityType * const_types)229 bool get_constraints_linearity(Index m, LinearityType* const_types){ 230 assert(m==_model->get_nb_cons()); 231 return false; 232 } 233 234 }; 235 236 237 238 #endif /* defined(__Gravity____IpoptProgram__) */ 239