1 // Copyright (C) 2008, 2011 International Business Machines and others.
2 // All Rights Reserved.
3 // This code is published under the Eclipse Public License.
4 //
5 // $Id$
6 //
7 // Authors:  Andreas Waechter             IBM    2008-08-31
8 
9 #include "IpInexactCq.hpp"
10 #include "IpIpoptNLP.hpp"
11 
12 namespace Ipopt
13 {
14 #if COIN_IPOPT_VERBOSITY > 0
15   static const Index dbg_verbosity = 0;
16 #endif
17 
18 
InexactCq(IpoptNLP * ip_nlp,IpoptData * ip_data,IpoptCalculatedQuantities * ip_cq)19   InexactCq::InexactCq(IpoptNLP* ip_nlp,
20                        IpoptData* ip_data,
21                        IpoptCalculatedQuantities* ip_cq)
22       :
23       ip_nlp_(ip_nlp),
24       ip_data_(ip_data),
25       ip_cq_(ip_cq),
26 
27       curr_jac_cdT_times_curr_cdminuss_cache_(1),
28       curr_scaling_slacks_cache_(1),
29       curr_slack_scaled_d_minus_s_cache_(1),
30       curr_scaled_Ac_norm_cache_(1),
31       slack_scaled_norm_cache_(6),
32       curr_W_times_vec_x_cache_(0), // ToDo: decide if we want this cached
33       curr_W_times_vec_s_cache_(0),
34       curr_Wu_x_cache_(1),
35       curr_Wu_s_cache_(1),
36       curr_uWu_cache_(1),
37       curr_jac_times_normal_c_cache_(1),
38       curr_jac_times_normal_d_cache_(1)
39   {
40     DBG_START_METH("InexactCq::InexactCq", dbg_verbosity);
41     DBG_ASSERT(ip_nlp_);
42     DBG_ASSERT(ip_data_);
43     DBG_ASSERT(ip_cq_);
44   }
45 
~InexactCq()46   InexactCq::~InexactCq()
47   {}
48 
49   void
RegisterOptions(const SmartPtr<RegisteredOptions> & roptions)50   InexactCq::RegisterOptions(const SmartPtr<RegisteredOptions>& roptions)
51   {
52     roptions->AddLowerBoundedNumberOption(
53       "slack_scale_max",
54       "Upper bound on slack-based scaling parameters.",
55       0.0, true, 1.,
56       "");
57   }
58 
59   bool
Initialize(const Journalist & jnlst,const OptionsList & options,const std::string & prefix)60   InexactCq::Initialize(const Journalist& jnlst,
61                         const OptionsList& options,
62                         const std::string& prefix)
63   {
64     options.GetNumericValue("slack_scale_max", slack_scale_max_, prefix);
65     return true;
66   }
67 
68   SmartPtr<const Vector>
curr_jac_cdT_times_curr_cdminuss()69   InexactCq::curr_jac_cdT_times_curr_cdminuss()
70   {
71     SmartPtr<const Vector> result;
72     SmartPtr<const Vector> x = ip_data_->curr()->x();
73     SmartPtr<const Vector> s = ip_data_->curr()->s();
74 
75     if (!curr_jac_cdT_times_curr_cdminuss_cache_.GetCachedResult2Dep(result, *x, *s)) {
76       // c part
77       SmartPtr<const Vector> curr_c = ip_cq_->curr_c();
78       SmartPtr<const Vector> curr_jac_cT_times_curr_c =
79         ip_cq_->curr_jac_cT_times_vec(*curr_c);
80 
81       // d minus s part
82       SmartPtr<const Vector> curr_d_minus_s = ip_cq_->curr_d_minus_s();
83       SmartPtr<const Vector> curr_jac_dT_times_curr_d_minus_s =
84         ip_cq_->curr_jac_dT_times_vec(*curr_d_minus_s);
85 
86       // add them
87       SmartPtr<Vector> tmp = curr_jac_cT_times_curr_c->MakeNew();
88       tmp->AddTwoVectors(1., *curr_jac_cT_times_curr_c,
89                          1., *curr_jac_dT_times_curr_d_minus_s, 0.);
90       result = ConstPtr(tmp);
91       curr_jac_cdT_times_curr_cdminuss_cache_.AddCachedResult2Dep(result, *x, *s);
92     }
93 
94     return result;
95   }
96 
97   SmartPtr<const Vector>
curr_scaling_slacks()98   InexactCq::curr_scaling_slacks()
99   {
100     DBG_START_METH("InexactCq::curr_scaling_slacks()",
101                    dbg_verbosity);
102     SmartPtr<const Vector> result;
103     SmartPtr<const Vector> s = ip_data_->curr()->s();
104 
105     if (!curr_scaling_slacks_cache_.GetCachedResult1Dep(result, *s)) {
106       SmartPtr<Vector> tmp = s->MakeNew();
107 
108       // Lower bounds
109       SmartPtr<const Matrix> Pd_L = ip_nlp_->Pd_L();
110       SmartPtr<const Vector> curr_slack_s_L = ip_cq_->curr_slack_s_L();
111       DBG_PRINT_MATRIX(1, "Pd_L", *Pd_L);
112       DBG_PRINT_VECTOR(1, "curr_slack_s_L", *curr_slack_s_L);
113       Pd_L->MultVector(1., *curr_slack_s_L, 0., *tmp);
114 
115       // Upper bounds
116       SmartPtr<const Matrix> Pd_U = ip_nlp_->Pd_U();
117       SmartPtr<const Vector> curr_slack_s_U = ip_cq_->curr_slack_s_U();
118       DBG_PRINT_MATRIX(1, "Pd_U", *Pd_U);
119       DBG_PRINT_VECTOR(1, "curr_slack_s_U", *curr_slack_s_U);
120       Pd_U->MultVector(1., *curr_slack_s_U, 1., *tmp);
121 
122       SmartPtr<Vector> tmp2 = tmp->MakeNew();
123       tmp2->Set(slack_scale_max_);
124       tmp->ElementWiseMin(*tmp2);
125 
126       result = ConstPtr(tmp);
127       curr_scaling_slacks_cache_.AddCachedResult1Dep(result, *s);
128     }
129 
130     return result;
131   }
132 
133   SmartPtr<const Vector>
curr_slack_scaled_d_minus_s()134   InexactCq::curr_slack_scaled_d_minus_s()
135   {
136     SmartPtr<const Vector> result;
137     SmartPtr<const Vector> x = ip_data_->curr()->x();
138     SmartPtr<const Vector> s = ip_data_->curr()->s();
139 
140     if (!curr_slack_scaled_d_minus_s_cache_.GetCachedResult2Dep(result, *x, *s)) {
141       SmartPtr<const Vector> d_minus_s = ip_cq_->curr_d_minus_s();
142       SmartPtr<const Vector> scaling_slacks = curr_scaling_slacks();
143 
144       SmartPtr<Vector> tmp = d_minus_s->MakeNewCopy();
145       tmp->ElementWiseMultiply(*scaling_slacks);
146       result = ConstPtr(tmp);
147       curr_slack_scaled_d_minus_s_cache_.AddCachedResult2Dep(result, *x, *s);
148     }
149 
150     return result;
151   }
152 
153   Number
curr_scaled_Ac_norm()154   InexactCq::curr_scaled_Ac_norm()
155   {
156     DBG_START_METH("InexactCq::curr_scaled_Ac_norm",
157                    dbg_verbosity);
158     Number result;
159 
160     SmartPtr<const Vector> x = ip_data_->curr()->x();
161     SmartPtr<const Vector> s = ip_data_->curr()->s();
162 
163     if (!curr_scaled_Ac_norm_cache_.GetCachedResult2Dep(result, *x, *s)) {
164       SmartPtr<const Vector> jac_cdT_times_curr_cdminuss =
165         curr_jac_cdT_times_curr_cdminuss();
166       DBG_PRINT_VECTOR(2, "jac_cdT_times_curr_cdminuss",
167                        *jac_cdT_times_curr_cdminuss);
168       SmartPtr<const Vector> slack_scaled_d_minus_s =
169         curr_slack_scaled_d_minus_s();
170       DBG_PRINT_VECTOR(2, "slack_scaled_d_minus_s",
171                        *slack_scaled_d_minus_s);
172       result = ip_cq_->CalcNormOfType(NORM_2, *jac_cdT_times_curr_cdminuss,
173                                       *slack_scaled_d_minus_s);
174 
175       curr_scaled_Ac_norm_cache_.AddCachedResult2Dep(result, *x, *s);
176     }
177 
178     return result;
179   }
180 
181   /** ||A||^2 */
182   Number
curr_scaled_A_norm2()183   InexactCq::curr_scaled_A_norm2()
184   {
185     DBG_START_METH("InexactCq::curr_scaled_A_norm",
186                    dbg_verbosity);
187     Number result;
188 
189     //if (!curr_scaled_A_norm_cache_.GetCachedResult(...)) {
190     result = 2;
191     //curr_scaled_A_norm_cache_.AddCachedResult(...);
192     //}
193 
194     return result;
195   }
196 
197   Number
slack_scaled_norm(const Vector & x,const Vector & s)198   InexactCq::slack_scaled_norm(const Vector& x, const Vector &s)
199   {
200     SmartPtr<const Vector> scaling_slacks = curr_scaling_slacks();
201 
202     Number result;
203     if (!slack_scaled_norm_cache_.GetCachedResult3Dep(result, *scaling_slacks, x, s)) {
204       SmartPtr<Vector> tmp = s.MakeNewCopy();
205       tmp->ElementWiseDivide(*scaling_slacks);
206       result = ip_cq_->CalcNormOfType(NORM_2, x, *tmp);
207       slack_scaled_norm_cache_.AddCachedResult3Dep(result, *scaling_slacks, x, s);
208     }
209 
210     return result;
211   }
212 
213   SmartPtr<const Vector>
curr_W_times_vec_x(const Vector & vec_x)214   InexactCq::curr_W_times_vec_x(const Vector& vec_x)
215   {
216     DBG_START_METH("InexactCq::curr_W_times_vec_x",
217                    dbg_verbosity);
218     SmartPtr<const Vector> result;
219 
220     SmartPtr<const SymMatrix> W = ip_data_->W();
221     Number pd_pert_x;
222     Number pd_pert_s;
223     Number pd_pert_c;
224     Number pd_pert_d;
225     ip_data_->getPDPert(pd_pert_x, pd_pert_s, pd_pert_c, pd_pert_d);
226 
227     std::vector<const TaggedObject*> tdeps(2);
228     tdeps[0] = &vec_x;
229     tdeps[1] = GetRawPtr(W);
230     std::vector<Number> sdeps(1);
231     sdeps[0] = pd_pert_x;
232 
233     if (!curr_W_times_vec_x_cache_.GetCachedResult(result, tdeps, sdeps)) {
234       DBG_PRINT_VECTOR(2, "vec_x", vec_x);
235       DBG_PRINT_MATRIX(2, "W", *W);
236       DBG_PRINT((2, "pd_pert_x = %e\n", pd_pert_x));
237       SmartPtr<Vector> tmp = vec_x.MakeNewCopy();
238       W->MultVector(1., vec_x, pd_pert_x, *tmp);
239       result = ConstPtr(tmp);
240       curr_W_times_vec_x_cache_.AddCachedResult(result, tdeps, sdeps);
241     }
242 
243     return result;
244   }
245 
246   SmartPtr<const Vector>
curr_W_times_vec_s(const Vector & vec_s)247   InexactCq::curr_W_times_vec_s(const Vector& vec_s)
248   {
249     SmartPtr<const Vector> result;
250 
251     SmartPtr<const Vector> sigma_s = ip_cq_->curr_sigma_s();
252     Number pd_pert_x;
253     Number pd_pert_s;
254     Number pd_pert_c;
255     Number pd_pert_d;
256     ip_data_->getPDPert(pd_pert_x, pd_pert_s, pd_pert_c, pd_pert_d);
257 
258     std::vector<const TaggedObject*> tdeps(2);
259     tdeps[0] = &vec_s;
260     tdeps[1] = GetRawPtr(sigma_s);
261     std::vector<Number> sdeps(1);
262     sdeps[0] = pd_pert_s;
263 
264     if (!curr_W_times_vec_s_cache_.GetCachedResult(result, tdeps, sdeps)) {
265       SmartPtr<Vector> tmp = vec_s.MakeNewCopy();
266       tmp->ElementWiseMultiply(*sigma_s);
267       if (pd_pert_s>0.) {
268         tmp->AddOneVector(pd_pert_s, vec_s, 1.);
269       }
270       result = ConstPtr(tmp);
271       curr_W_times_vec_s_cache_.AddCachedResult(result, tdeps, sdeps);
272     }
273 
274     return result;
275   }
276 
277   SmartPtr<const Vector>
curr_Wu_x()278   InexactCq::curr_Wu_x()
279   {
280     SmartPtr<const Vector> result;
281 
282     SmartPtr<const Vector> tangential_x = InexData().tangential_x();
283     SmartPtr<const SymMatrix> W = ip_data_->W();
284     Number pd_pert_x;
285     Number pd_pert_s;
286     Number pd_pert_c;
287     Number pd_pert_d;
288     ip_data_->getPDPert(pd_pert_x, pd_pert_s, pd_pert_c, pd_pert_d);
289 
290     std::vector<const TaggedObject*> tdeps(2);
291     tdeps[0] = GetRawPtr(tangential_x);
292     tdeps[1] = GetRawPtr(W);
293     std::vector<Number> sdeps(1);
294     sdeps[0] = pd_pert_x;
295 
296     if (!curr_Wu_x_cache_.GetCachedResult(result, tdeps, sdeps)) {
297       result = curr_W_times_vec_x(*tangential_x);
298       curr_Wu_x_cache_.AddCachedResult(result, tdeps, sdeps);
299     }
300 
301     return result;
302   }
303 
304   SmartPtr<const Vector>
curr_Wu_s()305   InexactCq::curr_Wu_s()
306   {
307     SmartPtr<const Vector> result;
308 
309     SmartPtr<const Vector> tangential_s = InexData().tangential_s();
310     SmartPtr<const Vector> sigma_s = ip_cq_->curr_sigma_s();
311     Number pd_pert_x;
312     Number pd_pert_s;
313     Number pd_pert_c;
314     Number pd_pert_d;
315     ip_data_->getPDPert(pd_pert_x, pd_pert_s, pd_pert_c, pd_pert_d);
316 
317     std::vector<const TaggedObject*> tdeps(2);
318     tdeps[0] = GetRawPtr(tangential_s);
319     tdeps[1] = GetRawPtr(sigma_s);
320     std::vector<Number> sdeps(1);
321     sdeps[0] = pd_pert_s;
322 
323     if (!curr_Wu_s_cache_.GetCachedResult(result, tdeps, sdeps)) {
324       result = curr_W_times_vec_s(*tangential_s);
325       curr_Wu_s_cache_.AddCachedResult(result, tdeps, sdeps);
326     }
327 
328     return result;
329   }
330 
331   Number
curr_uWu()332   InexactCq::curr_uWu()
333   {
334     Number result;
335 
336     SmartPtr<const Vector> tangential_x = InexData().tangential_x();
337     SmartPtr<const Vector> tangential_s = InexData().tangential_s();
338     SmartPtr<const Vector> Wu_x = curr_Wu_x();
339     SmartPtr<const Vector> Wu_s = curr_Wu_s();
340 
341     std::vector<const TaggedObject*> tdeps(4);
342     tdeps[0] = GetRawPtr(tangential_x);
343     tdeps[1] = GetRawPtr(tangential_s);
344     tdeps[2] = GetRawPtr(Wu_x);
345     tdeps[3] = GetRawPtr(Wu_s);
346 
347     if (!curr_uWu_cache_.GetCachedResult(result, tdeps)) {
348       result = Wu_x->Dot(*tangential_x) + Wu_s->Dot(*tangential_s);
349       curr_uWu_cache_.AddCachedResult(result, tdeps);
350     }
351 
352     return result;
353   }
354 
355   SmartPtr<const Vector>
curr_jac_times_normal_c()356   InexactCq::curr_jac_times_normal_c()
357   {
358     SmartPtr<const Vector> result;
359 
360     SmartPtr<const Vector> x = ip_data_->curr()->x();
361     SmartPtr<const Vector> normal_x = InexData().normal_x();
362 
363     if (!curr_jac_times_normal_c_cache_.GetCachedResult2Dep(result, *x, *normal_x)) {
364       result = ip_cq_->curr_jac_c_times_vec(*normal_x);
365       curr_jac_times_normal_c_cache_.AddCachedResult2Dep(result, *x, *normal_x);
366     }
367 
368     return result;
369   }
370 
371   SmartPtr<const Vector>
curr_jac_times_normal_d()372   InexactCq::curr_jac_times_normal_d()
373   {
374     SmartPtr<const Vector> result;
375 
376     SmartPtr<const Vector> x = ip_data_->curr()->x();
377     SmartPtr<const Vector> normal_x = InexData().normal_x();
378     SmartPtr<const Vector> normal_s = InexData().normal_s();
379 
380     if (!curr_jac_times_normal_d_cache_.GetCachedResult3Dep(result, *x, *normal_x, *normal_s)) {
381       SmartPtr<const Vector> Ax = ip_cq_->curr_jac_d_times_vec(*normal_x);
382       SmartPtr<Vector> tmp = Ax->MakeNew();
383       tmp->AddTwoVectors(1., *Ax, -1., *normal_s, 0.);
384 
385       result = ConstPtr(tmp);
386       curr_jac_times_normal_d_cache_.AddCachedResult3Dep(result, *x, *normal_x, *normal_s);
387     }
388 
389     return result;
390   }
391 
392 } // namespace Ipopt
393