source: branches/dev/Algorithm/IpOrigIpoptNLP.cpp @ 530

Last change on this file since 530 was 530, checked in by andreasw, 16 years ago

print unscaled initial and final variables for J_VECTOR

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.5 KB
Line 
1// Copyright (C) 2004, 2005 International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpOrigIpoptNLP.cpp 530 2005-10-02 21:27:29Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpOrigIpoptNLP.hpp"
10
11#ifdef HAVE_CMATH
12# include <cmath>
13#else
14# ifdef HAVE_MATH_H
15#  include <math.h>
16# else
17#  error "don't have header file for math"
18# endif
19#endif
20
21#ifdef HAVE_CASSERT
22# include <cassert>
23#else
24# ifdef HAVE_ASSERT_H
25#  include <assert.h>
26# else
27#  error "don't have header file for assert"
28# endif
29#endif
30
31namespace Ipopt
32{
33#ifdef IP_DEBUG
34  static const Index dbg_verbosity = 0;
35#endif
36
37  OrigIpoptNLP::OrigIpoptNLP(const SmartPtr<const Journalist>& jnlst,
38                             const SmartPtr<NLP>& nlp,
39                             const SmartPtr<NLPScalingObject>& nlp_scaling)
40      :
41      IpoptNLP(nlp_scaling),
42      jnlst_(jnlst),
43      nlp_(nlp),
44      x_space_(NULL),
45      f_cache_(1),
46      grad_f_cache_(1),
47      c_cache_(1),
48      jac_c_cache_(1),
49      d_cache_(1),
50      jac_d_cache_(1),
51      h_cache_(1),
52      f_evals_(0),
53      grad_f_evals_(0),
54      c_evals_(0),
55      jac_c_evals_(0),
56      d_evals_(0),
57      jac_d_evals_(0),
58      h_evals_(0),
59      initialized_(false)
60  {}
61
62  OrigIpoptNLP::~OrigIpoptNLP()
63  {}
64
65  void OrigIpoptNLP::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
66  {
67    roptions->AddLowerBoundedNumberOption(
68      "bound_relax_factor",
69      "Factor for initial relaxation of the bounds.",
70      0, false,
71      1e-8,
72      "Before start of the optimization, the bounds given by the user are "
73      "relaxed.  This option sets the factor for this relaxation.  If it "
74      "is set to zero, then then bounds relaxation is disabled. "
75      "(See Eqn.(35) in implmentation paper.)");
76    roptions->AddStringOption2(
77      "warm_start_same_structure",
78      "Indicates whether a problem with a structure identical to the previous one is to be solved.",
79      "no",
80      "no", "Assume this is a new problem.",
81      "yes", "Assume this is problem has known structure",
82      "If \"yes\" is chosen, then the algorithm assumes that an NLP is now to "
83      "be solved, whose strcture is identical to one that already was "
84      "considered (with the same NLP object).");
85  }
86
87  bool OrigIpoptNLP::Initialize(const Journalist& jnlst,
88                                const OptionsList& options,
89                                const std::string& prefix)
90  {
91    options.GetNumericValue("bound_relax_factor", bound_relax_factor_, prefix);
92    options.GetBoolValue("warm_start_same_structure",
93                         warm_start_same_structure_, prefix);
94
95    // Reset the function evaluation counters (for warm start)
96    f_evals_=0;
97    grad_f_evals_=0;
98    c_evals_=0;
99    jac_c_evals_=0;
100    d_evals_=0;
101    jac_d_evals_=0;
102    h_evals_=0;
103
104    // Reset the cache entries belonging to a dummy dependency.  This
105    // is required for repeated solve, since the cache is not updated
106    // if a dimension is zero
107    std::vector<const TaggedObject*> deps(1);
108    deps[0] = NULL;
109    std::vector<Number> sdeps(0);
110    c_cache_.InvalidateResult(deps, sdeps);
111    d_cache_.InvalidateResult(deps, sdeps);
112    jac_c_cache_.InvalidateResult(deps, sdeps);
113    jac_d_cache_.InvalidateResult(deps, sdeps);
114
115    if (!nlp_->ProcessOptions(options, prefix)) {
116      return false;
117    }
118
119    initialized_ = true;
120    return IpoptNLP::Initialize(jnlst, options, prefix);
121  }
122
123  bool OrigIpoptNLP::InitializeStructures(SmartPtr<Vector>& x,
124                                          bool init_x,
125                                          SmartPtr<Vector>& y_c,
126                                          bool init_y_c,
127                                          SmartPtr<Vector>& y_d,
128                                          bool init_y_d,
129                                          SmartPtr<Vector>& z_L,
130                                          bool init_z_L,
131                                          SmartPtr<Vector>& z_U,
132                                          bool init_z_U,
133                                          SmartPtr<Vector>& v_L,
134                                          SmartPtr<Vector>& v_U
135                                         )
136  {
137    DBG_START_METH("OrigIpoptNLP::InitializeStructures", dbg_verbosity);
138    DBG_ASSERT(initialized_);
139    bool retValue;
140
141    if (!warm_start_same_structure_) {
142
143      retValue = nlp_->GetSpaces(x_space_, c_space_, d_space_,
144                                 x_l_space_, px_l_space_,
145                                 x_u_space_, px_u_space_,
146                                 d_l_space_, pd_l_space_,
147                                 d_u_space_, pd_u_space_,
148                                 jac_c_space_, jac_d_space_,
149                                 h_space_);
150
151      if (!retValue) {
152        return false;
153      }
154
155      NLP_scaling()->DetermineScaling(x_space_,
156                                      c_space_, d_space_,
157                                      jac_c_space_, jac_d_space_,
158                                      h_space_,
159                                      scaled_jac_c_space_, scaled_jac_d_space_,
160                                      scaled_h_space_);
161
162      ASSERT_EXCEPTION(x_space_->Dim() >= c_space_->Dim(), TOO_FEW_DOF,
163                       "Too few degrees of freedom!");
164
165      // cannot have any null pointers, want zero length vectors
166      // instead of null - this will later need to be changed for _h;
167      retValue = (IsValid(x_space_) && IsValid(c_space_) && IsValid(d_space_)
168                  && IsValid(x_l_space_) && IsValid(px_l_space_)
169                  && IsValid(x_u_space_) && IsValid(px_u_space_)
170                  && IsValid(d_u_space_) && IsValid(pd_u_space_)
171                  && IsValid(d_l_space_) && IsValid(pd_l_space_)
172                  && IsValid(jac_c_space_) && IsValid(jac_d_space_)
173                  && IsValid(h_space_)
174                  && IsValid(scaled_jac_c_space_)
175                  && IsValid(scaled_jac_d_space_)
176                  && IsValid(scaled_h_space_));
177
178      DBG_ASSERT(retValue && "Model cannot return null vector or matrix prototypes or spaces,"
179                 " please return zero length vectors instead");
180    }
181    else {
182      ASSERT_EXCEPTION(IsValid(x_space_), INVALID_WARMSTART,
183                       "OrigIpoptNLP called with warm_start_same_structure, but the problem is solved for the first time.");
184    }
185
186    // Create the bounds structures
187    SmartPtr<Vector> x_L = x_l_space_->MakeNew();
188    SmartPtr<Matrix> Px_L = px_l_space_->MakeNew();
189    SmartPtr<Vector> x_U = x_u_space_->MakeNew();
190    SmartPtr<Matrix> Px_U = px_u_space_->MakeNew();
191    SmartPtr<Vector> d_L = d_l_space_->MakeNew();
192    SmartPtr<Matrix> Pd_L = pd_l_space_->MakeNew();
193    SmartPtr<Vector> d_U = d_u_space_->MakeNew();
194    SmartPtr<Matrix> Pd_U = pd_u_space_->MakeNew();
195
196    retValue = nlp_->GetBoundsInformation(*Px_L, *x_L, *Px_U, *x_U,
197                                          *Pd_L, *d_L, *Pd_U, *d_U);
198
199    if (!retValue) {
200      return false;
201    }
202
203    relax_bounds(-bound_relax_factor_, *x_L);
204    relax_bounds( bound_relax_factor_, *x_U);
205    relax_bounds(-bound_relax_factor_, *d_L);
206    relax_bounds( bound_relax_factor_, *d_U);
207
208    x_L_ = ConstPtr(x_L);
209    Px_L_ = ConstPtr(Px_L);
210    x_U_ = ConstPtr(x_U);
211    Px_U_ = ConstPtr(Px_U);
212    d_L_ = ConstPtr(d_L);
213    Pd_L_ = ConstPtr(Pd_L);
214    d_U_ = ConstPtr(d_U);
215    Pd_U_ = ConstPtr(Pd_U);
216
217    // now create and store the scaled bounds
218    x_L_ = NLP_scaling()->apply_vector_scaling_x_LU(*Px_L_, x_L_, *x_space_);
219    x_U_ = NLP_scaling()->apply_vector_scaling_x_LU(*Px_U_, x_U_, *x_space_);
220    d_L_ = NLP_scaling()->apply_vector_scaling_d_LU(*Pd_L_, d_L_, *d_space_);
221    d_U_ = NLP_scaling()->apply_vector_scaling_d_LU(*Pd_U_, d_U_, *d_space_);
222
223    // Create the iterates structures
224    x = x_space_->MakeNew();
225    y_c = c_space_->MakeNew();
226    y_d = d_space_->MakeNew();
227    z_L = x_l_space_->MakeNew();
228    z_U = x_u_space_->MakeNew();
229    v_L = d_l_space_->MakeNew();
230    v_U = d_u_space_->MakeNew();
231
232    retValue = nlp_->GetStartingPoint(GetRawPtr(x), init_x,
233                                      GetRawPtr(y_c), init_y_c,
234                                      GetRawPtr(y_d), init_y_d,
235                                      GetRawPtr(z_L), init_z_L,
236                                      GetRawPtr(z_U), init_z_U);
237
238    if (!retValue) {
239      return false;
240    }
241
242
243    Number obj_scal = NLP_scaling()->apply_obj_scaling(1.);
244    if (init_x) {
245      x->Print(*jnlst_, J_VECTOR, J_INITIALIZATION, "initial x unscaled");
246      if (NLP_scaling()->have_x_scaling()) {
247        x = NLP_scaling()->apply_vector_scaling_x_NonConst(ConstPtr(x));
248      }
249    }
250    if (init_y_c) {
251      y_c->Print(*jnlst_, J_VECTOR, J_INITIALIZATION, "initial y_c unscaled");
252      if (NLP_scaling()->have_c_scaling()) {
253        y_c = NLP_scaling()->unapply_vector_scaling_c_NonConst(ConstPtr(y_c));
254      }
255      if (obj_scal!=1.) {
256        y_c->Scal(obj_scal);
257      }
258    }
259    if (init_y_d) {
260      y_d->Print(*jnlst_, J_VECTOR, J_INITIALIZATION, "initial y_d unscaled");
261      if (NLP_scaling()->have_d_scaling()) {
262        y_d = NLP_scaling()->unapply_vector_scaling_d_NonConst(ConstPtr(y_d));
263      }
264      if (obj_scal!=1.) {
265        y_d->Scal(obj_scal);
266      }
267    }
268    if (init_z_L) {
269      z_L->Print(*jnlst_, J_VECTOR, J_INITIALIZATION, "initial z_L unscaled");
270      if (NLP_scaling()->have_x_scaling()) {
271        z_L = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_L_, ConstPtr(z_L), *x_space_);
272      }
273      if (obj_scal!=1.) {
274        z_L->Scal(obj_scal);
275      }
276    }
277    if (init_z_U) {
278      z_U->Print(*jnlst_, J_VECTOR, J_INITIALIZATION, "initial z_U unscaled");
279      if (NLP_scaling()->have_x_scaling()) {
280        z_U = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_U_, ConstPtr(z_U), *x_space_);
281      }
282      if (obj_scal!=1.) {
283        z_U->Scal(obj_scal);
284      }
285    }
286
287    return true;
288  }
289
290  void
291  OrigIpoptNLP::relax_bounds(Number bound_relax_factor, Vector& bounds)
292  {
293    DBG_START_METH("OrigIpoptNLP::relax_bounds", dbg_verbosity);
294    if (bound_relax_factor!=0.) {
295      SmartPtr<Vector> tmp = bounds.MakeNew();
296      tmp->Copy(bounds);
297      tmp->ElementWiseAbs();
298      SmartPtr<Vector> ones = bounds.MakeNew();
299      ones->Set(1.);
300      tmp->ElementWiseMax(*ones);
301      DBG_PRINT((1, "bound_relax_factor = %e", bound_relax_factor));
302      DBG_PRINT_VECTOR(2, "tmp", *tmp);
303      DBG_PRINT_VECTOR(2, "bounds before", bounds);
304      bounds.Axpy(bound_relax_factor, *tmp);
305      DBG_PRINT_VECTOR(2, "bounds after", bounds);
306    }
307  }
308
309  Number OrigIpoptNLP::f(const Vector& x)
310  {
311    DBG_START_METH("OrigIpoptNLP::f", dbg_verbosity);
312    Number ret = 0.0;
313    DBG_PRINT((2, "x.Tag = %d\n", x.GetTag()));
314    if (!f_cache_.GetCachedResult1Dep(ret, &x)) {
315      f_evals_++;
316      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
317      f_eval_time_.Start();
318      bool success = nlp_->Eval_f(*unscaled_x, ret);
319      f_eval_time_.End();
320      DBG_PRINT((1, "success = %d ret = %e\n", success, ret));
321      ASSERT_EXCEPTION(success && IsFiniteNumber(ret), Eval_Error,
322                       "Error evaluating the objective function");
323      ret = NLP_scaling()->apply_obj_scaling(ret);
324      f_cache_.AddCachedResult1Dep(ret, &x);
325    }
326
327    return ret;
328  }
329
330  Number OrigIpoptNLP::f(const Vector& x, Number mu)
331  {
332    assert(false && "ERROR: This method is only a placeholder for f(mu) and should not be called");
333    return 0.;
334  }
335
336  SmartPtr<const Vector> OrigIpoptNLP::grad_f(const Vector& x)
337  {
338    SmartPtr<Vector> unscaled_grad_f;
339    SmartPtr<const Vector> retValue;
340    if (!grad_f_cache_.GetCachedResult1Dep(retValue, &x)) {
341      grad_f_evals_++;
342      unscaled_grad_f = x_space_->MakeNew();
343
344      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
345      grad_f_eval_time_.Start();
346      bool success = nlp_->Eval_grad_f(*unscaled_x, *unscaled_grad_f);
347      grad_f_eval_time_.End();
348      ASSERT_EXCEPTION(success && IsFiniteNumber(unscaled_grad_f->Nrm2()),
349                       Eval_Error, "Error evaluating the gradient of the objective function");
350      retValue = NLP_scaling()->apply_grad_obj_scaling(ConstPtr(unscaled_grad_f));
351      grad_f_cache_.AddCachedResult1Dep(retValue, &x);
352    }
353
354    return retValue;
355  }
356
357  SmartPtr<const Vector> OrigIpoptNLP::grad_f(const Vector& x, Number mu)
358  {
359    assert(false && "ERROR: This method is only a placeholder for grad_f(mu) and should not be called");
360    return NULL;
361  }
362
363  /** Equality constraint residual */
364  SmartPtr<const Vector> OrigIpoptNLP::c(const Vector& x)
365  {
366    SmartPtr<const Vector> retValue;
367    if (c_space_->Dim()==0) {
368      // We do this caching of an empty vector so that the returned
369      // Vector has always the same tag (this might make a difference
370      // in cases where only the constraints are supposed to change...
371      SmartPtr<const Vector> dep = NULL;
372      if (!c_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
373        retValue = c_space_->MakeNew();
374        c_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
375      }
376    }
377    else {
378      if (!c_cache_.GetCachedResult1Dep(retValue, x)) {
379        SmartPtr<Vector> unscaled_c = c_space_->MakeNew();
380        c_evals_++;
381        SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
382        c_eval_time_.Start();
383        bool success = nlp_->Eval_c(*unscaled_x, *unscaled_c);
384        c_eval_time_.End();
385        ASSERT_EXCEPTION(success && IsFiniteNumber(unscaled_c->Nrm2()),
386                         Eval_Error, "Error evaluating the equality constraints");
387        retValue = NLP_scaling()->apply_vector_scaling_c(ConstPtr(unscaled_c));
388        c_cache_.AddCachedResult1Dep(retValue, x);
389      }
390    }
391
392    return retValue;
393  }
394
395  SmartPtr<const Vector> OrigIpoptNLP::d(const Vector& x)
396  {
397    DBG_START_METH("OrigIpoptNLP::d", dbg_verbosity);
398    SmartPtr<const Vector> retValue;
399    if (d_space_->Dim()==0) {
400      // We do this caching of an empty vector so that the returned
401      // Vector has always the same tag (this might make a difference
402      // in cases where only the constraints are supposed to change...
403      SmartPtr<const Vector> dep = NULL;
404      if (!d_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
405        retValue = d_space_->MakeNew();
406        d_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
407      }
408    }
409    else {
410      if (!d_cache_.GetCachedResult1Dep(retValue, x)) {
411        d_evals_++;
412        SmartPtr<Vector> unscaled_d = d_space_->MakeNew();
413
414        DBG_PRINT_VECTOR(2, "scaled_x", x);
415        SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
416        d_eval_time_.Start();
417        bool success = nlp_->Eval_d(*unscaled_x, *unscaled_d);
418        d_eval_time_.End();
419        DBG_PRINT_VECTOR(2, "unscaled_d", *unscaled_d);
420        ASSERT_EXCEPTION(success && IsFiniteNumber(unscaled_d->Nrm2()),
421                         Eval_Error, "Error evaluating the inequality constraints");
422        retValue = NLP_scaling()->apply_vector_scaling_d(ConstPtr(unscaled_d));
423        d_cache_.AddCachedResult1Dep(retValue, x);
424      }
425    }
426
427    return retValue;
428  }
429
430  SmartPtr<const Matrix> OrigIpoptNLP::jac_c(const Vector& x)
431  {
432    SmartPtr<const Matrix> retValue;
433    if (c_space_->Dim()==0) {
434      // We do this caching of an empty vector so that the returned
435      // Matrix has always the same tag
436      SmartPtr<const Vector> dep = NULL;
437      if (!jac_c_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
438        SmartPtr<Matrix> unscaled_jac_c = jac_c_space_->MakeNew();
439        retValue = NLP_scaling()->apply_jac_c_scaling(ConstPtr(unscaled_jac_c));
440        jac_c_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
441      }
442    }
443    else {
444      if (!jac_c_cache_.GetCachedResult1Dep(retValue, x)) {
445        jac_c_evals_++;
446        SmartPtr<Matrix> unscaled_jac_c = jac_c_space_->MakeNew();
447
448        SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
449        jac_c_eval_time_.Start();
450        bool success = nlp_->Eval_jac_c(*unscaled_x, *unscaled_jac_c);
451        jac_c_eval_time_.End();
452        ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the jacobian of the equality constraints");
453        retValue = NLP_scaling()->apply_jac_c_scaling(ConstPtr(unscaled_jac_c));
454        jac_c_cache_.AddCachedResult1Dep(retValue, x);
455      }
456    }
457
458    return retValue;
459  }
460
461  SmartPtr<const Matrix> OrigIpoptNLP::jac_d(const Vector& x)
462  {
463    SmartPtr<const Matrix> retValue;
464    if (d_space_->Dim()==0) {
465      // We do this caching of an empty vector so that the returned
466      // Matrix has always the same tag
467      SmartPtr<const Vector> dep = NULL;
468      if (!jac_d_cache_.GetCachedResult1Dep(retValue, GetRawPtr(dep))) {
469        SmartPtr<Matrix> unscaled_jac_d = jac_d_space_->MakeNew();
470        retValue = NLP_scaling()->apply_jac_d_scaling(ConstPtr(unscaled_jac_d));
471        jac_d_cache_.AddCachedResult1Dep(retValue, GetRawPtr(dep));
472      }
473    }
474    else {
475      if (!jac_d_cache_.GetCachedResult1Dep(retValue, x)) {
476        jac_d_evals_++;
477        SmartPtr<Matrix> unscaled_jac_d = jac_d_space_->MakeNew();
478
479        SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
480        jac_d_eval_time_.Start();
481        bool success = nlp_->Eval_jac_d(*unscaled_x, *unscaled_jac_d);
482        jac_d_eval_time_.End();
483        ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the jacobian of the inequality constraints");
484        retValue = NLP_scaling()->apply_jac_d_scaling(ConstPtr(unscaled_jac_d));
485        jac_d_cache_.AddCachedResult1Dep(retValue, x);
486      }
487    }
488
489    return retValue;
490  }
491
492  SmartPtr<const SymMatrix> OrigIpoptNLP::h(const Vector& x,
493      Number obj_factor,
494      const Vector& yc,
495      const Vector& yd)
496  {
497    std::vector<const TaggedObject*> deps(3);
498    deps[0] = &x;
499    deps[1] = &yc;
500    deps[2] = &yd;
501    std::vector<Number> scalar_deps(1);
502    scalar_deps[0] = obj_factor;
503
504    SmartPtr<SymMatrix> unscaled_h;
505    SmartPtr<const SymMatrix> retValue;
506    if (!h_cache_.GetCachedResult(retValue, deps, scalar_deps)) {
507      h_evals_++;
508      unscaled_h = h_space_->MakeNewSymMatrix();
509
510      SmartPtr<const Vector> unscaled_x = NLP_scaling()->unapply_vector_scaling_x(&x);
511      SmartPtr<const Vector> unscaled_yc = NLP_scaling()->apply_vector_scaling_c(&yc);
512      SmartPtr<const Vector> unscaled_yd = NLP_scaling()->apply_vector_scaling_d(&yd);
513      Number scaled_obj_factor = NLP_scaling()->apply_obj_scaling(obj_factor);
514      h_eval_time_.Start();
515      bool success = nlp_->Eval_h(*unscaled_x, scaled_obj_factor, *unscaled_yc, *unscaled_yd, *unscaled_h);
516      h_eval_time_.End();
517      ASSERT_EXCEPTION(success, Eval_Error, "Error evaluating the hessian of the lagrangian");
518      retValue = NLP_scaling()->apply_hessian_scaling(ConstPtr(unscaled_h));
519      h_cache_.AddCachedResult(retValue, deps, scalar_deps);
520    }
521
522    return retValue;
523  }
524
525  SmartPtr<const SymMatrix> OrigIpoptNLP::h(const Vector& x,
526      Number obj_factor,
527      const Vector& yc,
528      const Vector& yd,
529      Number mu)
530  {
531    assert(false &&
532           "ERROR: This method is only a for h(mu) and should not be called");
533    return NULL;
534  }
535
536
537  void OrigIpoptNLP::GetSpaces(SmartPtr<const VectorSpace>& x_space,
538                               SmartPtr<const VectorSpace>& c_space,
539                               SmartPtr<const VectorSpace>& d_space,
540                               SmartPtr<const VectorSpace>& x_l_space,
541                               SmartPtr<const MatrixSpace>& px_l_space,
542                               SmartPtr<const VectorSpace>& x_u_space,
543                               SmartPtr<const MatrixSpace>& px_u_space,
544                               SmartPtr<const VectorSpace>& d_l_space,
545                               SmartPtr<const MatrixSpace>& pd_l_space,
546                               SmartPtr<const VectorSpace>& d_u_space,
547                               SmartPtr<const MatrixSpace>& pd_u_space,
548                               SmartPtr<const MatrixSpace>& Jac_c_space,
549                               SmartPtr<const MatrixSpace>& Jac_d_space,
550                               SmartPtr<const SymMatrixSpace>& Hess_lagrangian_space)
551  {
552    // Make sure that we already have all the pointers
553    DBG_ASSERT(IsValid(x_space_) &&
554               IsValid(c_space_) &&
555               IsValid(d_space_) &&
556               IsValid(x_l_space_) &&
557               IsValid(px_l_space_) &&
558               IsValid(x_u_space_) &&
559               IsValid(px_u_space_) &&
560               IsValid(d_l_space_) &&
561               IsValid(pd_l_space_) &&
562               IsValid(d_u_space_) &&
563               IsValid(pd_u_space_) &&
564               IsValid(scaled_jac_c_space_) &&
565               IsValid(scaled_jac_d_space_) &&
566               IsValid(scaled_h_space_));
567
568    DBG_ASSERT(IsValid(NLP_scaling()));
569
570    x_space = x_space_;
571    c_space = c_space_;
572    d_space = d_space_;
573    x_l_space = x_l_space_;
574    px_l_space = px_l_space_;
575    x_u_space = x_u_space_;
576    px_u_space = px_u_space_;
577    d_l_space = d_l_space_;
578    pd_l_space = pd_l_space_;
579    d_u_space = d_u_space_;
580    pd_u_space = pd_u_space_;
581    Jac_c_space = scaled_jac_c_space_;
582    Jac_d_space = scaled_jac_d_space_;
583    Hess_lagrangian_space = scaled_h_space_;
584  }
585
586  void OrigIpoptNLP::FinalizeSolution(SolverReturn status,
587                                      const Vector& x, const Vector& z_L, const Vector& z_U,
588                                      const Vector& c, const Vector& d,
589                                      const Vector& y_c, const Vector& y_d,
590                                      Number obj_value)
591  {
592    DBG_START_METH("OrigIpoptNLP::FinalizeSolution", dbg_verbosity);
593    // need to submit the unscaled solution back to the nlp
594    SmartPtr<const Vector> unscaled_x =
595      NLP_scaling()->unapply_vector_scaling_x(&x);
596    SmartPtr<const Vector> unscaled_c =
597      NLP_scaling()->unapply_vector_scaling_c(&c);
598    SmartPtr<const Vector> unscaled_d =
599      NLP_scaling()->unapply_vector_scaling_d(&d);
600    const Number unscaled_obj = NLP_scaling()->unapply_obj_scaling(obj_value);
601
602    SmartPtr<const Vector> unscaled_z_L;
603    SmartPtr<const Vector> unscaled_z_U;
604    SmartPtr<const Vector> unscaled_y_c;
605    SmartPtr<const Vector> unscaled_y_d;
606
607    // The objective function scaling factor also appears in the constraints
608    Number obj_unscale_factor = NLP_scaling()->unapply_obj_scaling(1.);
609    if (obj_unscale_factor!=1.) {
610      SmartPtr<Vector> tmp = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_L_, &z_L, *x_space_);
611      tmp->Scal(obj_unscale_factor);
612      unscaled_z_L = ConstPtr(tmp);
613
614      tmp = NLP_scaling()->apply_vector_scaling_x_LU_NonConst(*Px_U_, &z_U, *x_space_);
615      tmp->Scal(obj_unscale_factor);
616      unscaled_z_U = ConstPtr(tmp);
617
618      tmp = NLP_scaling()->apply_vector_scaling_c_NonConst(&y_c);
619      tmp->Scal(obj_unscale_factor);
620      unscaled_y_c = ConstPtr(tmp);
621
622      tmp = NLP_scaling()->apply_vector_scaling_d_NonConst(&y_d);
623      tmp->Scal(obj_unscale_factor);
624      unscaled_y_d = ConstPtr(tmp);
625    }
626    else {
627      unscaled_z_L = NLP_scaling()->apply_vector_scaling_x_LU(*Px_L_, &z_L, *x_space_);
628      unscaled_z_U = NLP_scaling()->apply_vector_scaling_x_LU(*Px_U_, &z_U, *x_space_);
629      unscaled_y_c = NLP_scaling()->apply_vector_scaling_c(&y_c);
630      unscaled_y_d = NLP_scaling()->apply_vector_scaling_d(&y_d);
631    }
632
633    unscaled_x->Print(*jnlst_, J_VECTOR, J_SOLUTION, "final x unscaled");
634    unscaled_y_c->Print(*jnlst_, J_VECTOR, J_SOLUTION, "final y_c unscaled");
635    unscaled_y_d->Print(*jnlst_, J_VECTOR, J_SOLUTION, "final y_d unscaled");
636    unscaled_z_L->Print(*jnlst_, J_VECTOR, J_SOLUTION, "final z_L unscaled");
637    unscaled_z_U->Print(*jnlst_, J_VECTOR, J_SOLUTION, "final z_U unscaled");
638
639    nlp_->FinalizeSolution(status, *unscaled_x,
640                           *unscaled_z_L, *unscaled_z_U,
641                           *unscaled_c, *unscaled_d,
642                           *unscaled_y_c, *unscaled_y_d,
643                           unscaled_obj);
644  }
645
646  void OrigIpoptNLP::AdjustVariableBounds(const Vector& new_x_L, const Vector& new_x_U,
647                                          const Vector& new_d_L, const Vector& new_d_U)
648  {
649    x_L_ = new_x_L.MakeNewCopy();
650    x_U_ = new_x_U.MakeNewCopy();
651    d_L_ = new_d_L.MakeNewCopy();
652    d_U_ = new_d_U.MakeNewCopy();
653  }
654
655  void
656  OrigIpoptNLP::PrintTimingStatistics(
657    Journalist& jnlst,
658    EJournalLevel level,
659    EJournalCategory category) const
660  {
661    if (!jnlst.ProduceOutput(level, category))
662      return;
663
664    jnlst.Printf(level, category,
665                 "Function Evaluations................: %10.3f\n",
666                 f_eval_time_.TotalTime()+
667                 c_eval_time_.TotalTime()+
668                 d_eval_time_.TotalTime()+
669                 jac_c_eval_time_.TotalTime()+
670                 jac_d_eval_time_.TotalTime()+
671                 h_eval_time_.TotalTime());
672    jnlst.Printf(level, category,
673                 " Objective function.................: %10.3f\n",
674                 f_eval_time_.TotalTime());
675    jnlst.Printf(level, category,
676                 " Equality constraints...............: %10.3f\n",
677                 c_eval_time_.TotalTime());
678    jnlst.Printf(level, category,
679                 " Inequality constraints.............: %10.3f\n",
680                 d_eval_time_.TotalTime());
681    jnlst.Printf(level, category,
682                 " Equality constraint Jacobian.......: %10.3f\n",
683                 jac_c_eval_time_.TotalTime());
684    jnlst.Printf(level, category,
685                 " Inequality constraint Jacobian.....: %10.3f\n",
686                 jac_d_eval_time_.TotalTime());
687    jnlst.Printf(level, category,
688                 " Lagrangian Hessian.................: %10.3f\n",
689                 h_eval_time_.TotalTime());
690  }
691
692} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.