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