CppAD: A C++ Algorithmic Differentiation Package  20130918
cppad_ipopt_nlp.cpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 /* --------------------------------------------------------------------------
00003 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
00004 
00005 CppAD is distributed under multiple licenses. This distribution is under
00006 the terms of the 
00007                     Eclipse Public License Version 1.0.
00008 
00009 A copy of this license is included in the COPYING file of this distribution.
00010 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00011 -------------------------------------------------------------------------- */
00012 # include <limits>
00013 
00014 # include "cppad_ipopt_nlp.hpp"
00015 # include "sparse_map2vec.hpp"
00016 # include "jac_g_map.hpp"
00017 # include "hes_fg_map.hpp"
00018 # include "vec_fun_pattern.hpp"
00019 # include "fun_record.hpp"
00020 
00021 
00022 /// If 0 tracing is off, otherwise tracing is on.
00023 # define  CPPAD_IPOPT_NLP_TRACE 0
00024 
00025 # if CPPAD_IPOPT_NLP_TRACE
00026 # include <cstdio>
00027 # endif
00028 
00029 // ---------------------------------------------------------------------------
00030 namespace cppad_ipopt {
00031 // ---------------------------------------------------------------------------
00032 /*!
00033 \{
00034 \file cppad_ipopt_nlp.cpp
00035 \brief Member functions for the cppad_ipopt_nlp class.
00036 */
00037 
00038 
00039 /*! 
00040 Constructor for the \ref Nonlinear_Programming_Problem.
00041 
00042 \param n
00043 dimension of the domain space for f(x) and g(x).
00044 
00045 \param m
00046 dimension of the range space for g(x)
00047 
00048 \param x_i
00049 initial value of x during the optimization procedure (size n).
00050 
00051 \param x_l
00052 lower limit for x (size n).
00053 
00054 \param x_u
00055 upper limit for x (size n).
00056 
00057 \param g_l
00058 lower limit for g(x) (size m).
00059 
00060 \param g_u
00061 upper limit for g(x) (size m).
00062 
00063 \param fg_info 
00064 pointer to base class version of derived class object used to get 
00065 information about the user's representation for f(x) and g(x).
00066 (The object pointed to must not be deleted before this cppad_ipopt_nlp object).
00067 
00068 \param solution
00069 pointer to object where final results are stored.
00070 (The object pointed to must not be deleted before this cppad_ipopt_nlp object).
00071 
00072 \par Constants
00073 The following values are set by the constructor and are \c const
00074 or effectively \c const; i.e., they are set by the constructor and should
00075 not be changed:
00076 \verbatim
00077      n_, m_, x_i_, x_l_, x_u_, g_l_, g_u_, K_, L_, p_, q_, retape_,
00078      pattern_jac_r_, pattern_hes_r_, index_jac_g_, index_hes_fg_,
00079      nnz_jac_g_, iRow_jac_g_, jCol_jac_g_,
00080      nnz_h_lag_, iRow_h_lag_, jCol_h_lag_,
00081 \endverbatim
00082 In addition, the function calls <tt>fg_info->set_n(n)</tt>
00083 and <tt>fg_info->set_m(m)</tt> are used to set the values of \c n
00084 and \c m in \c fg_info. 
00085 
00086 \par Variables
00087 The following arrays have fixed size which is set during this constructor:
00088 
00089 \li \c tape_ok_ has size \c K_. It is initialized as true for indices
00090 \c k such that <tt>retape[k]</tt> is false.  
00091 
00092 \li \c r_fun_ has size \c K_. It is initilaize with the default
00093 \c ADFun constructor. Then, for indices \c k such that 
00094 <tt>retape[k]</tt> is false, the operation sequence corresponding
00095 to \f$ r_k (u) \f$ is stored in <tt>r_fun_[k]</tt>.
00096 
00097 \li \c I_ has size equal to the maximum of <tt>p[k]</tt> w.r.t \c k.
00098 
00099 \li \c J_ has size equal to the maximum of <tt>q[k]</tt> w.r.t \c k.
00100 
00101 \par NDEBUG
00102 If the preprocessor symbol \c NEBUG is not defined,
00103 certain of the assumptions about the function calls of the form
00104 \verbatim
00105      fg_info->index(k, ell, I, J)
00106 \endverbatim
00107 are checked to make sure they hold.
00108 */
00109 cppad_ipopt_nlp::cppad_ipopt_nlp(
00110      size_t n                         , 
00111      size_t m                         ,
00112      const NumberVector    &x_i       ,
00113      const NumberVector    &x_l       ,
00114      const NumberVector    &x_u       ,
00115      const NumberVector    &g_l       ,
00116      const NumberVector    &g_u       ,
00117      cppad_ipopt_fg_info*   fg_info   ,
00118      cppad_ipopt_solution*  solution )
00119      : n_ ( n ),
00120        m_ ( m ),
00121        x_i_ ( x_i ),
00122        x_l_ ( x_l ),
00123        x_u_ ( x_u ),
00124        g_l_ ( g_l ),
00125        g_u_ ( g_u ),
00126        fg_info_ ( fg_info ) ,
00127        solution_ (solution) ,
00128        infinity_ ( std::numeric_limits<Number>::infinity() )
00129 {    size_t k;
00130 
00131      // set information needed in cppad_ipopt_fg_info
00132      fg_info_->set_n(n);
00133      fg_info_->set_m(m);
00134 
00135      // get information from derived class version of fg_info
00136      K_ = fg_info_->number_functions();
00137      L_.resize(K_);
00138      p_.resize(K_);
00139      q_.resize(K_);
00140      r_fun_.resize(K_);
00141      retape_.resize(K_);
00142      tape_ok_.resize(K_);
00143      pattern_jac_r_.resize(K_);
00144      pattern_hes_r_.resize(K_);
00145      size_t max_p      = 0;
00146      size_t max_q      = 0;
00147      for(k = 0; k < K_; k++)
00148      {    L_[k]       = fg_info_->number_terms(k);
00149           p_[k]       = fg_info_->range_size(k);
00150           q_[k]       = fg_info_->domain_size(k);
00151           retape_[k]  = fg_info_->retape(k);
00152           max_p       = std::max(max_p, p_[k]);
00153           max_q       = std::max(max_q, q_[k]);
00154           pattern_jac_r_[k].resize( p_[k] * q_[k] );
00155           pattern_hes_r_[k].resize( q_[k] * q_[k] );
00156      }
00157      I_.resize(max_p);
00158      J_.resize(max_q);
00159 # ifndef NDEBUG
00160      size_t i, j, ell;
00161      // check for valid range and domain indices
00162      for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00163      {
00164           for( i = 0; i < p_[k]; i++)
00165                I_[i] = m+1; // an invalid range index
00166           for( j = 0; j < q_[k]; j++)
00167                J_[j] = n; // an invalid domain index
00168           fg_info_->index(k, ell, I_, J_);   
00169           for( i = 0; i < p_[k]; i++) if( I_[i] > m )
00170           {    std::cerr << "k=" << k << ", ell=" << ell 
00171                << ", I[" << i << "]=" << I_[i] << std::endl;
00172                CPPAD_ASSERT_KNOWN( I_[i] <= m,
00173                "cppad_ipopt_nlp: invalid value in index vector I"
00174                );
00175           }
00176           for( j = 0; j < q_[k]; j++) if( J_[j] >= n )
00177           {    std::cerr << "k=" << k << ", ell=" << ell 
00178                << ", J[" << j << "]=" << J_[j] << std::endl;
00179                CPPAD_ASSERT_KNOWN( J_[j] < n,
00180                "cppad_ipopt_nlp: invalid value in index vector J"
00181                );
00182           }
00183      }
00184 # endif
00185      // record r[k] for functions that do not need retaping
00186      for(k = 0; k < K_; k++)
00187      {    tape_ok_[k] = false;
00188           if( ! retape_[k] )
00189           {    // Operation sequence does not depend on value 
00190                // of u so record it once here in the constructor.
00191                fg_info_->index(k, 0, I_, J_);
00192                fun_record(
00193                     fg_info_        ,   // inputs
00194                     k               ,
00195                     p_              ,
00196                     q_              ,
00197                     n_              ,
00198                     x_i_            ,
00199                     J_              ,
00200                     r_fun_              // output
00201                );
00202                // take time to optimize because only recording once
00203                r_fun_[k].optimize();
00204                // ok and will stay that way
00205                tape_ok_[k] = true;
00206           }
00207      }
00208 
00209      // compute a sparsity patterns for each r_k (u)
00210      vec_fun_pattern(
00211           K_, p_, q_, retape_, r_fun_,      // inputs 
00212           pattern_jac_r_, pattern_hes_r_    // outputs
00213      );
00214 
00215      // mapping from (i,j) to Ipopt sparsity index for Jacobian of g
00216      jac_g_map(
00217           fg_info_, m_, n_, K_, L_, p_, q_, pattern_jac_r_,   // inputs
00218           I_, J_,                                             // work
00219           index_jac_g_                                        // outputs
00220      );
00221 
00222      // mapping from (i,j) to Ipopt sparsity index for Hessian of Lagragian
00223      hes_fg_map(
00224           fg_info_, m_, n_, K_, L_, p_, q_, pattern_hes_r_,   // inputs
00225           I_, J_,                                             // work
00226           index_hes_fg_                                       // outputs
00227      );
00228      
00229      // Compute Ipopt sparsity structure for Jacobian of g 
00230      sparse_map2vec(
00231           index_jac_g_,                         // inputs
00232           nnz_jac_g_, iRow_jac_g_, jCol_jac_g_  // outputs
00233      );
00234 
00235      // Compute Ipopt sparsity structure for Hessian of Lagragian
00236      sparse_map2vec(
00237           index_hes_fg_,                        // inputs
00238           nnz_h_lag_, iRow_h_lag_, jCol_h_lag_  // outputs
00239      );
00240 }
00241 
00242 /// The destructor takes no special action.
00243 cppad_ipopt_nlp::~cppad_ipopt_nlp()
00244 {}
00245 
00246 /*!
00247 Return dimension information about optimization problem.
00248 
00249 \param[out] n
00250 is set to the value \c n_.
00251 
00252 \param[out] m
00253 is set to the value \c m_.
00254 
00255 \param[out] nnz_jac_g
00256 is set to the value of \c nnz_jac_g_.
00257 
00258 \param[out] nnz_h_lag
00259 is set to the vlaue of \c nnz_h_lag_.
00260 
00261 \param[out] index_style
00262 is set to C_STYLE; i.e., zeoro based indexing is used in the
00263 information passed to Ipopt.
00264 */
00265 bool cppad_ipopt_nlp::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00266                          Index& nnz_h_lag, IndexStyleEnum& index_style)
00267 {
00268      n = n_;
00269      m = m_;
00270      nnz_jac_g = nnz_jac_g_;
00271      nnz_h_lag = nnz_h_lag_;
00272 
00273      // use the fortran index style for row/col entries
00274      index_style = C_STYLE;
00275 
00276      return true;
00277 }
00278 
00279 /*!
00280 Return bound information about optimization problem.
00281 
00282 \param[in] n
00283 is the dimension of the domain space for f(x) and g(x); i.e.,
00284 it must be equal to \c n_.
00285 
00286 \param[out] x_l
00287 is a vector of size \c n.
00288 The input value of its elements does not matter.
00289 On output, it is a copy of the lower bound for \f$ x \f$; i.e.,
00290 \c x_l_.
00291 
00292 \param[out] x_u
00293 is a vector of size \c n.
00294 The input value of its elements does not matter.
00295 On output, it is a copy of the upper bound for \f$ x \f$; i.e.,
00296 \c x_u_.
00297 
00298 \param[in] m
00299 is the dimension of the range space for g(x). i.e.,
00300 it must be equal to \c m_.
00301 
00302 \param[out] g_l
00303 is a vector of size \c m.
00304 The input value of its elements does not matter.
00305 On output, it is a copy of the lower bound for \f$ g(x) \f$; i.e., \c g_l_.
00306 
00307 \param[out] g_u
00308 is a vector of size \c m.
00309 The input value of its elements does not matter.
00310 On output, it is a copy of the upper bound for \f$ g(x) \f$; i.e, \c g_u_.
00311 */
00312 bool cppad_ipopt_nlp::get_bounds_info(Index n, Number* x_l, Number* x_u,
00313                             Index m, Number* g_l, Number* g_u)
00314 {    size_t i, j;
00315      // here, the n and m we gave IPOPT in get_nlp_info are passed back 
00316      CPPAD_ASSERT_UNKNOWN(size_t(n) == n_);
00317      CPPAD_ASSERT_UNKNOWN(size_t(m) == m_);
00318 
00319      // pass back bounds
00320      for(j = 0; j < n_; j++)
00321      {    x_l[j] = x_l_[j];
00322           x_u[j] = x_u_[j];
00323      }
00324      for(i = 0; i < m_; i++)
00325      {    g_l[i] = g_l_[i];
00326           g_u[i] = g_u_[i];
00327      }
00328      
00329      return true;
00330 }
00331 
00332 /*!
00333 Return initial x value where optimiation is started.
00334 
00335 \param[in] n
00336 must be equal to the domain dimension for f(x) and g(x); i.e.,
00337 it must be equal to \c n_.
00338 
00339 \param[in] init_x
00340 must be equal to true.
00341 
00342 \param[out] x
00343 is a vector of size \c n.
00344 The input value of its elements does not matter.
00345 On output, it is a copy of the initial value for \f$ x \f$; i.e. \c x_i_.
00346 
00347 \param[in] init_z
00348 must be equal to false.
00349 
00350 \param z_L
00351 is not used.
00352 
00353 \param z_U
00354 is not used.
00355 
00356 \param[in] m
00357 must be equal to the range dimension for g(x); i.e.,
00358 it must be equal to \c m_.
00359 
00360 \param init_lambda
00361 must be equal to false.
00362 
00363 \param lambda
00364 is not used.
00365 */
00366 bool cppad_ipopt_nlp::get_starting_point(Index n, bool init_x, Number* x,
00367                                bool init_z, Number* z_L, Number* z_U,
00368                                Index m, bool init_lambda,
00369                                Number* lambda)
00370 {    size_t j;
00371 
00372      CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00373      CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ );
00374      CPPAD_ASSERT_UNKNOWN(init_x == true);
00375      CPPAD_ASSERT_UNKNOWN(init_z == false);
00376      CPPAD_ASSERT_UNKNOWN(init_lambda == false);
00377 
00378      for(j = 0; j < n_; j++)
00379           x[j] = x_i_[j];
00380 
00381      return true;
00382 }
00383 
00384 /*!
00385 Evaluate the objective fucntion f(x).
00386 
00387 \param[in] n
00388 is the dimension of the argument space for f(x); i.e., must be equal \c n_.
00389 
00390 \param[in] x
00391 is a vector of size \c n containing the point at which to evaluate
00392 the function f(x).
00393 
00394 \param[in] new_x
00395 is false if the previous call to any one of the 
00396 \ref Evaluation_Methods used the same value for \c x.
00397 
00398 \param[out] obj_value
00399 is the value of the objective f(x) at this value of \c x.
00400 
00401 \return
00402 The return value is always true; see \ref Evaluation_Methods.
00403 
00404 \par Efficiency
00405 This routine could be more efficient 
00406 (for certain when when L[k] > 1 and retape[k] is true)
00407 if the users also provided a version 
00408 of the function <tt>fg_info->eval_r(k, u)</tt> where \c u was of type
00409 \c NumberVector.
00410 */
00411 bool cppad_ipopt_nlp::eval_f(
00412      Index n, const Number* x, bool new_x, Number& obj_value
00413 )
00414 {
00415      CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00416 
00417      size_t iobj, j, k, ell;
00418 
00419      // initialize summation
00420      obj_value = 0.;
00421 
00422      // update tape_ok_ flag
00423      for(k = 0; k < K_; k++) 
00424      {    if( retape_[k] && (new_x || L_[k] > 1) )
00425                tape_ok_[k] = false;
00426      }
00427 
00428      for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00429      {    fg_info_->index(k, ell, I_, J_);
00430           for(iobj = 0; iobj < p_[k]; iobj++) if( I_[iobj] == 0 )
00431           {    if( ! tape_ok_[k] )
00432                {    // Record r_k for value of u corresponding to x
00433                     fun_record(
00434                          fg_info_        ,   // inputs
00435                          k               ,
00436                          p_              ,
00437                          q_              ,
00438                          n_              ,
00439                          x               ,
00440                          J_              ,
00441                          r_fun_             // output
00442                     );
00443                     tape_ok_[k] = ! (retape_[k] || L_[k] > 1);
00444                }
00445                NumberVector u(q_[k]);
00446                NumberVector r(p_[k]);
00447                for(j = 0; j < q_[k]; j++)
00448                {    CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00449                     u[j]   = x[ J_[j] ];
00450                }
00451                r          = r_fun_[k].Forward(0, u);
00452                obj_value += r[iobj];
00453           }
00454      }
00455 # if CPPAD_IPOPT_NLP_TRACE
00456      using std::printf;
00457      for(j = 0; j < n_; j++)
00458           printf("cppad_ipopt_nlp::eval_f::x[%d] = %20.14g\n", j, x[j]);
00459      printf("cppad_ipopt_nlp::eval_f::obj_value = %20.14g\n", obj_value);
00460 # endif
00461 # ifndef NDEBUG
00462      CPPAD_ASSERT_KNOWN(
00463           (-infinity_ < obj_value) && (obj_value < infinity_),
00464           "cppad_ipopt_nlp::eval_f:: objective value is not finite"
00465      );
00466 # endif
00467      return true;
00468 }
00469 
00470 /*!
00471 Evaluate the gradient of f(x).
00472 
00473 \param[in] n
00474 is the dimension of the argument space for f(x); i.e., must be equal \c n_.
00475 
00476 \param[in] x
00477 has a vector of size \c n containing the point at which to evaluate
00478 the gradient of f(x).
00479 
00480 \param[in] new_x
00481 is false if the previous call to any one of the 
00482 \ref Evaluation_Methods used the same value for \c x.
00483 
00484 \param[out] grad_f
00485 is a vector of size \c n.
00486 The input value of its elements does not matter.
00487 The output value of its elements is the gradient of f(x) 
00488 at this value of.
00489 
00490 \return
00491 The return value is always true; see \ref Evaluation_Methods.
00492 */
00493 bool cppad_ipopt_nlp::eval_grad_f(
00494      Index n, const Number* x, bool new_x, Number* grad_f
00495 )
00496 {    CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00497 
00498      size_t iobj, i, j, k, ell;
00499 
00500      // initialize summation
00501      for(j = 0; j < n_; j++)
00502           grad_f[j] = 0.;
00503 
00504      // update tape_ok_ flag
00505      for(k = 0; k < K_; k++) 
00506      {    if( retape_[k] && (new_x || L_[k] > 1) )
00507                tape_ok_[k] = false;
00508      }
00509 
00510      for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00511      {    fg_info_->index(k, ell, I_, J_);
00512           for(iobj = 0; iobj < p_[k]; iobj++) if( I_[iobj] == 0 )
00513           {    if( ! tape_ok_[k] )
00514                {    // Record r_k for value of u corresponding to x
00515                     fun_record(
00516                          fg_info_        ,   // inputs
00517                          k               ,
00518                          p_              ,
00519                          q_              ,
00520                          n_              ,
00521                          x               ,
00522                          J_              ,
00523                          r_fun_              // output
00524                     );
00525                     tape_ok_[k] = ! (retape_[k] || L_[k] > 1);
00526                }
00527                NumberVector u(q_[k]);
00528                NumberVector w(p_[k]);
00529                NumberVector r_grad(q_[k]);
00530                for(j = 0; j < q_[k]; j++)
00531                {    CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00532                     u[j]   = x[ J_[j] ];
00533                }
00534                r_fun_[k].Forward(0, u);
00535                for(i = 0; i < p_[k]; i++)
00536                     w[i] = 0.;
00537                w[iobj]    = 1.;
00538                r_grad     = r_fun_[k].Reverse(1, w);
00539                for(j = 0; j < q_[k]; j++)
00540                {    CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00541                     grad_f[ J_[j] ]  += r_grad[j];
00542                }
00543           }
00544      }
00545 # if CPPAD_IPOPT_NLP_TRACE
00546      using std::printf;
00547      for(j = 0; j < n_; j++) printf(
00548      "cppad_ipopt_nlp::eval_grad_f::x[%d] = %20.14g\n", j, x[j]
00549      );
00550      for(j = 0; j < n_; j++) printf(
00551      "cppad_ipopt_nlp::eval_grad_f::grad_f[%d] = %20.14g\n", j, grad_f[j]
00552      );
00553 # endif
00554 # ifndef NDEBUG
00555      for(j = 0; j < n_; j++) CPPAD_ASSERT_KNOWN(
00556           (-infinity_ < grad_f[j]) && (grad_f[j] < infinity_),
00557           "cppad_ipopt_nlp::grad_f:: gradient of objective is not finite"
00558      );
00559 # endif
00560      return true;
00561 }
00562 
00563 /*!
00564 Evaluate the function g(x).
00565 
00566 \param[in] n
00567 is the dimension of the argument space for g(x); i.e., must be equal \c n_.
00568 
00569 \param[in] x
00570 has a vector of size \c n containing the point at which to evaluate
00571 the constraint function g(x).
00572 
00573 \param[in] new_x
00574 is false if the previous call to any one of the 
00575 \ref Evaluation_Methods used the same value for \c x.
00576 
00577 \param[in] m
00578 is the dimension of the range space for g(x); i.e., must be equal to \c m_.
00579 
00580 \param[out] g
00581 is a vector of size \c m.
00582 The input value of its elements does not matter.
00583 The output value of its elements is 
00584 the value of the function g(x) at this value of \c x.
00585 
00586 \return
00587 The return value is always true; see \ref Evaluation_Methods.
00588 */
00589 bool cppad_ipopt_nlp::eval_g(
00590      Index n, const Number* x, bool new_x, Index m, Number* g
00591 )
00592 {    CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00593 
00594      size_t i, j, k, ell;
00595 
00596      // initialize summation
00597      for(i = 0; i < m_; i++)
00598           g[i] = 0.;
00599 
00600      // update tape_ok_ flag
00601      for(k = 0; k < K_; k++) 
00602      {    if( retape_[k] && (new_x || L_[k] > 1) )
00603                tape_ok_[k] = false;
00604      }
00605 
00606      for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00607      {    fg_info_->index(k, ell, I_, J_);
00608           if( ! tape_ok_[k] )
00609           {    // Record r_k for value of u corresponding to x
00610                fun_record(
00611                     fg_info_        ,   // inputs
00612                     k               ,
00613                     p_              ,
00614                     q_              ,
00615                     n_              ,
00616                     x               ,
00617                     J_              ,
00618                     r_fun_              // output
00619                );
00620           }
00621           tape_ok_[k] = ! (retape_[k] || L_[k] > 1);
00622           NumberVector u(q_[k]);
00623           NumberVector r(p_[k]);
00624           for(j = 0; j < q_[k]; j++)
00625           {    CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00626                u[j]   = x[ J_[j] ];
00627           }
00628           r   = r_fun_[k].Forward(0, u);
00629           for(i = 0; i < p_[k]; i++)
00630           {    CPPAD_ASSERT_UNKNOWN( I_[i] <= m_ );
00631                if( I_[i] >= 1 )
00632                     g[ I_[i] - 1 ] += r[i];
00633           }
00634      }
00635 # if CPPAD_IPOPT_NLP_TRACE
00636      using std::printf;
00637      for(j = 0; j < n_; j++)
00638           printf("cppad_ipopt_nlp::eval_g::x[%d] = %20.14g\n", j, x[j]);
00639      for(i = 0; i < m_; i++)
00640           printf("cppad_ipopt_nlp::eval_g::g[%d] = %20.14g\n", i, g[i]);
00641 # endif
00642 # ifndef NDEBUG
00643      for(i = 0; i < m_; i++) CPPAD_ASSERT_KNOWN(
00644           (-infinity_ < g[i]) && (g[i] < infinity_),
00645           "cppad_ipopt_nlp::eval_g:: not all constraints are not finite"
00646      );
00647 # endif
00648      return true;
00649 }
00650 
00651 /*!
00652 Evaluate the Jacobian of g(x).
00653 
00654 \param[in] n
00655 is the dimension of the argument space for g(x); i.e., must be equal \c n_.
00656 
00657 \param x
00658 if \c values is not \c NULL,
00659 \c x is a vector of size \c n containing the point at which to evaluate
00660 the gradient of g(x).
00661 
00662 \param[in] new_x
00663 is false if the previous call to any one of the 
00664 \ref Evaluation_Methods used the same value for \c x.
00665 
00666 \param[in] m
00667 is the dimension of the range space for g(x); i.e., must be equal to \c m_.
00668 
00669 \param[in] nele_jac
00670 is the number of possibly non-zero elements in the Jacobian of g(x);
00671 i.e., must be equal to \c nnz_jac_g_.
00672 
00673 \param iRow
00674 if \c values is not \c NULL, \c iRow is not defined.
00675 if \c values is \c NULL, \c iRow
00676 is a vector with size \c nele_jac.
00677 The input value of its elements does not matter.
00678 On output, 
00679 For <tt>k = 0 , ... , nele_jac-1, iRow[k]</tt> is the 
00680 base zero row index for the 
00681 k-th possibly non-zero entry in the Jacobian of g(x).
00682 
00683 \param jCol
00684 if \c values is not \c NULL, \c jCol is not defined.
00685 if \c values is \c NULL, \c jCol
00686 is a vector with size \c nele_jac.
00687 The input value of its elements does not matter.
00688 On output, 
00689 For <tt>k = 0 , ... , nele_jac-1, jCol[k]</tt> is the 
00690 base zero column index for the 
00691 k-th possibly non-zero entry in the Jacobian of g(x).
00692 
00693 \param values
00694 if \c values is not \c NULL, \c values
00695 is a vector with size \c nele_jac.
00696 The input value of its elements does not matter.
00697 On output, 
00698 For <tt>k = 0 , ... , nele_jac-1, values[k]</tt> is the 
00699 value for the 
00700 k-th possibly non-zero entry in the Jacobian of g(x).
00701 
00702 \return
00703 The return value is always true; see \ref Evaluation_Methods.
00704 */
00705 bool cppad_ipopt_nlp::eval_jac_g(Index n, const Number* x, bool new_x,
00706                        Index m, Index nele_jac, Index* iRow, Index *jCol,
00707                        Number* values)
00708 {    CPPAD_ASSERT_UNKNOWN(size_t(m)          == m_ );
00709      CPPAD_ASSERT_UNKNOWN(size_t(n)          == n_ );
00710      CPPAD_ASSERT_UNKNOWN( size_t(nele_jac)  == nnz_jac_g_ );
00711 
00712      size_t i, j, k, ell, l;
00713      std::map<size_t,size_t>::iterator index_ij;
00714 
00715 
00716      if (values == NULL) 
00717      {    for(k = 0; k < nnz_jac_g_; k++)
00718           {    iRow[k] = iRow_jac_g_[k];
00719                jCol[k] = jCol_jac_g_[k];
00720           }
00721           return true;
00722      }
00723 
00724      // initialize summation
00725      l = nnz_jac_g_;
00726      while(l--)
00727           values[l] = 0.;
00728 
00729      // update tape_ok_ flag
00730      for(k = 0; k < K_; k++) 
00731      {    if( retape_[k] && (new_x || L_[k] > 1) )
00732                tape_ok_[k] = false;
00733      }
00734 
00735      for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00736      {    fg_info_->index(k, ell, I_, J_);
00737           if( ! tape_ok_[k] )
00738           {    // Record r_k for value of u corresponding to x
00739                fun_record(
00740                     fg_info_        ,   // inputs
00741                     k               ,
00742                     p_              ,
00743                     q_              ,
00744                     n_              ,
00745                     x               ,
00746                     J_              ,
00747                     r_fun_              // output
00748                );
00749           }
00750           tape_ok_[k] = ! (retape_[k] || L_[k] > 1);
00751           NumberVector u(q_[k]);
00752           NumberVector jac_r(p_[k] * q_[k]);
00753           for(j = 0; j < q_[k]; j++)
00754           {    CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00755                u[j]   = x[ J_[j] ];
00756           }
00757           if( retape_[k] )
00758                jac_r = r_fun_[k].Jacobian(u);
00759           else jac_r = r_fun_[k].SparseJacobian(u, pattern_jac_r_[k]);
00760           for(i = 0; i < p_[k]; i++) if( I_[i] != 0 )
00761           {    CPPAD_ASSERT_UNKNOWN( I_[i] <= m_ );
00762                for(j = 0; j < q_[k]; j++)
00763                {    index_ij = index_jac_g_[I_[i]-1].find(J_[j]);
00764                     if( index_ij != index_jac_g_[I_[i]-1].end() )
00765                     {    l          = index_ij->second;
00766                          values[l] += jac_r[i * q_[k] + j];
00767                     }
00768                     else CPPAD_ASSERT_UNKNOWN(
00769                          jac_r[i * q_[k] + j] == 0.
00770                     );
00771                }
00772           }
00773      }
00774 # ifndef NDEBUG
00775      for(l = 0; l < nnz_jac_g_; l++) CPPAD_ASSERT_KNOWN(
00776           (-infinity_ < values[l]) && (values[l] < infinity_),
00777           "cppad_ipopt_nlp::eval_jac_g:: a component of "
00778           "gradient of g is not finite"
00779      );
00780 # endif
00781 
00782      return true;
00783 }
00784 
00785 /*!
00786 Evaluate the Hessian of the Lagragian
00787 
00788 \section The_Hessian_of_the_Lagragian The Hessian of the Lagragian
00789 The Hessian of the Lagragian is defined as
00790 \f[
00791 H(x, \sigma, \lambda ) 
00792 =
00793 \sigma \nabla^2 f(x) + \sum_{i=0}^{m-1} \lambda_i \nabla^2 g(x)_i
00794 \f]
00795 
00796 \param[in] n
00797 is the dimension of the argument space for g(x); i.e., must be equal \c n_.
00798 
00799 \param x
00800 if \c values is not \c NULL, \c x
00801 is a vector of size \c n containing the point at which to evaluate
00802 the Hessian of the Lagrangian.
00803 
00804 \param[in] new_x
00805 is false if the previous call to any one of the 
00806 \ref Evaluation_Methods used the same value for \c x.
00807 
00808 \param[in] obj_factor
00809 the value \f$ \sigma \f$ multiplying the Hessian of
00810 f(x) in the expression for \ref The_Hessian_of_the_Lagragian.
00811 
00812 \param[in] m
00813 is the dimension of the range space for g(x); i.e., must be equal to \c m_.
00814 
00815 \param[in] lambda
00816 if \c values is not \c NULL, \c lambda
00817 is a vector of size \c m specifing the value of \f$ \lambda \f$
00818 in the expression for \ref The_Hessian_of_the_Lagragian.
00819 
00820 \param[in] new_lambda
00821 is true if the previous call to \c eval_h had the same value for
00822 \c lambda and false otherwise.
00823 (Not currently used.)
00824 
00825 \param[in] nele_hess
00826 is the number of possibly non-zero elements in the Hessian of the Lagragian;
00827 i.e., must be equal to \c nnz_h_lag_.
00828 
00829 \param iRow
00830 if \c values is not \c NULL, \c iRow is not defined.
00831 if \c values is \c NULL, \c iRow
00832 is a vector with size \c nele_jac.
00833 The input value of its elements does not matter.
00834 On output, 
00835 For <tt>k = 0 , ... , nele_jac-1, iRow[k]</tt> is the 
00836 base zero row index for the 
00837 k-th possibly non-zero entry in the Hessian of the Lagragian.
00838 
00839 \param jCol
00840 if \c values is not \c NULL, \c jCol is not defined.
00841 if \c values is \c NULL, \c jCol
00842 is a vector with size \c nele_jac.
00843 The input value of its elements does not matter.
00844 On output, 
00845 For <tt>k = 0 , ... , nele_jac-1, jCol[k]</tt> is the 
00846 base zero column index for the 
00847 k-th possibly non-zero entry in the Hessian of the Lagragian.
00848 
00849 \param values
00850 if \c values is not \c NULL, it
00851 is a vector with size \c nele_jac.
00852 The input value of its elements does not matter.
00853 On output, 
00854 For <tt>k = 0 , ... , nele_jac-1, values[k]</tt> is the 
00855 value for the 
00856 k-th possibly non-zero entry in the Hessian of the Lagragian.
00857 
00858 \return
00859 The return value is always true; see \ref Evaluation_Methods.
00860 */
00861 bool cppad_ipopt_nlp::eval_h(Index n, const Number* x, bool new_x,
00862                    Number obj_factor, Index m, const Number* lambda,
00863                    bool new_lambda, Index nele_hess, Index* iRow,
00864                    Index* jCol, Number* values)
00865 {    CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ );
00866      CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00867 
00868      size_t i, j, k, ell, l;
00869      std::map<size_t,size_t>::iterator index_ij;
00870 
00871      if (values == NULL) 
00872      {    for(k = 0; k < nnz_h_lag_; k++)
00873           {    iRow[k] = iRow_h_lag_[k];
00874                jCol[k] = jCol_h_lag_[k];
00875           }
00876           return true;
00877      }
00878 
00879      // initialize summation
00880      l = nnz_h_lag_;
00881      while(l--)
00882           values[l] = 0.;
00883 
00884      // update tape_ok_ flag
00885      for(k = 0; k < K_; k++) 
00886      {    if( retape_[k] && (new_x || L_[k] > 1) )
00887                tape_ok_[k] = false;
00888      }
00889 
00890      for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00891      {    fg_info_->index(k, ell, I_, J_);
00892           bool in_use = false;
00893           for(i = 0; i < p_[k]; i++)
00894           {    if( I_[i] == 0 )
00895                     in_use |= obj_factor > 0.;
00896                else in_use |= lambda[ I_[i] - 1 ] > 0;
00897           }
00898           if( in_use )
00899           {
00900                if( ! tape_ok_[k]  )
00901                {    // Record r_k for value of u corresponding to x
00902                     fun_record(
00903                          fg_info_        ,   // inputs
00904                          k               ,
00905                          p_              ,
00906                          q_              ,
00907                          n_              ,
00908                          x               ,
00909                          J_              ,
00910                          r_fun_              // output
00911                     );
00912                     tape_ok_[k] = ! (retape_[k] || L_[k] > 1);
00913                }
00914                NumberVector w(p_[k]);
00915                NumberVector r_hes(q_[k] * q_[k]);
00916                NumberVector u(q_[k]);
00917                for(j = 0; j < q_[k]; j++)
00918                {    CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00919                     u[j]   = x[ J_[j] ];
00920                }
00921                for(i = 0; i < p_[k]; i++)
00922                {    CPPAD_ASSERT_UNKNOWN( I_[i] <= m_ );
00923                     if( I_[i] == 0 )
00924                          w[i] = obj_factor;
00925                     else w[i] = lambda[ I_[i] - 1 ];
00926                }
00927                if( retape_[k] )
00928                     r_hes = r_fun_[k].Hessian(u, w);
00929                else r_hes = r_fun_[k].SparseHessian(
00930                          u, w, pattern_hes_r_[k]
00931                );
00932                for(i = 0; i < q_[k]; i++) for(j = 0; j < q_[k]; j++) 
00933                if( J_[j] <= J_[i] ) 
00934                {    index_ij = index_hes_fg_[J_[i]].find(J_[j]);
00935                     if( index_ij != index_hes_fg_[J_[i]].end() )
00936                     {    l          = index_ij->second;
00937                          values[l] += r_hes[i * q_[k] + j];
00938                     }
00939                     else CPPAD_ASSERT_UNKNOWN(
00940                          r_hes[i * q_[k] + j] == 0.
00941                     );
00942                }
00943           }
00944      }
00945 # ifndef NDEBUG
00946      for(l = 0; l < nnz_h_lag_; l++) CPPAD_ASSERT_KNOWN(
00947           (-infinity_ < values[l]) && (values[l] < infinity_),
00948           "cppad_ipopt_nlp::eval_h:: a component of "
00949           "Hessian of Lagragian is not finite"
00950      );
00951 # endif
00952      return true;
00953 }
00954 
00955 /*!
00956 Pass solution information from Ipopt to users solution structure.
00957 
00958 \param[in] status
00959 is value that the Ipopt solution status
00960 which gets mapped to a correponding value for 
00961 \n
00962 <tt>solution_->status</tt>.
00963 
00964 \param[in] n
00965 is the dimension of the domain space for f(x) and g(x); i.e.,
00966 it must be equal to \c n_.
00967 
00968 \param[in] x
00969 is a vector with size \c n specifing the final solution.
00970 \n
00971 <tt>solution_->x</tt> is set to be a vector with size \c n
00972 and to have the same element values.
00973 
00974 \param[in] z_L
00975 is a vector with size \c n specifing the Lagragian multipliers for the
00976 constraint \f$ x^l \leq x \f$.
00977 \n
00978 <tt>solution_->z_l</tt> is set to be a vector with size \c n
00979 and to have the same element values.
00980 
00981 \param[in] z_U
00982 is a vector with size \c n specifing the Lagragian multipliers for the
00983 constraint \f$ x \leq x^u \f$.
00984 \n
00985 <tt>solution_->z_u</tt> is set to be a vector with size \c n
00986 and to have the same element values.
00987 
00988 \param[in] m
00989 is the dimension of the range space for g(x). i.e.,
00990 it must be equal to \c m_.
00991 
00992 \param[in] g
00993 is a vector with size \c m containing the value of the constraint function
00994 g(x) at the final solution for \c x.
00995 \n
00996 <tt>solution_->g</tt> is set to be a vector with size \c m
00997 and to have the same element values.
00998 
00999 \param[in] lambda
01000 is a vector with size \c m specifing the Lagragian multipliers for the
01001 constraints \f$ g^l \leq g(x) \leq g^u \f$.
01002 \n
01003 <tt>solution_->lambda</tt> is set to be a vector with size \c m
01004 and to have the same element values.
01005 
01006 \param[in] obj_value
01007 is the value of the objective function f(x) at the final solution for \c x.
01008 \n
01009 <tt>solution_->obj_value</tt> is set to have the same value.
01010 
01011 \param[in] ip_data
01012 is unspecified (by Ipopt) and hence not used.
01013 
01014 \param[in] ip_cq
01015 is unspecified (by Ipopt) and hence not used.
01016 
01017 \par solution_[out]
01018 the pointer \c solution_ , which is equal to the pointer \c solution
01019 in the constructor for \c cppad_ipopt_nlp,
01020 is used to set output values (see documentation above).
01021 */
01022 void cppad_ipopt_nlp::finalize_solution(
01023      Ipopt::SolverReturn               status    ,
01024      Index                             n         , 
01025      const Number*                     x         , 
01026      const Number*                     z_L       , 
01027      const Number*                     z_U       ,
01028      Index                             m         , 
01029      const Number*                     g         , 
01030      const Number*                     lambda    ,
01031      Number                            obj_value ,
01032      const Ipopt::IpoptData*           ip_data   ,
01033      Ipopt::IpoptCalculatedQuantities* ip_cq
01034 )
01035 {    size_t i, j;
01036 
01037      CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
01038      CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ );
01039 
01040      switch(status)
01041      {    // convert status from Ipopt enum to cppad_ipopt_solution enum
01042           case Ipopt::SUCCESS:
01043           solution_->status = 
01044                cppad_ipopt_solution::success;
01045           break;
01046 
01047           case Ipopt::MAXITER_EXCEEDED:
01048           solution_->status = 
01049                cppad_ipopt_solution::maxiter_exceeded;
01050           break;
01051 
01052           case Ipopt::STOP_AT_TINY_STEP:
01053           solution_->status = 
01054                cppad_ipopt_solution::stop_at_tiny_step;
01055           break;
01056 
01057           case Ipopt::STOP_AT_ACCEPTABLE_POINT:
01058           solution_->status = 
01059                cppad_ipopt_solution::stop_at_acceptable_point;
01060           break;
01061 
01062           case Ipopt::LOCAL_INFEASIBILITY:
01063           solution_->status = 
01064                cppad_ipopt_solution::local_infeasibility;
01065           break;
01066 
01067           case Ipopt::USER_REQUESTED_STOP:
01068           solution_->status = 
01069                cppad_ipopt_solution::user_requested_stop;
01070           break;
01071 
01072           case Ipopt::DIVERGING_ITERATES:
01073           solution_->status = 
01074                cppad_ipopt_solution::diverging_iterates;
01075           break;
01076 
01077           case Ipopt::RESTORATION_FAILURE:
01078           solution_->status = 
01079                cppad_ipopt_solution::restoration_failure;
01080           break;
01081 
01082           case Ipopt::ERROR_IN_STEP_COMPUTATION:
01083           solution_->status = 
01084                cppad_ipopt_solution::error_in_step_computation;
01085           break;
01086 
01087           case Ipopt::INVALID_NUMBER_DETECTED:
01088           solution_->status = 
01089                cppad_ipopt_solution::invalid_number_detected;
01090           break;
01091 
01092           case Ipopt::INTERNAL_ERROR:
01093           solution_->status = 
01094                cppad_ipopt_solution::internal_error;
01095           break;
01096 
01097           default:
01098           solution_->status = 
01099                cppad_ipopt_solution::unknown;
01100      }
01101 
01102      solution_->x.resize(n_);
01103      solution_->z_l.resize(n_);
01104      solution_->z_u.resize(n_);
01105      for(j = 0; j < n_; j++)
01106      {    solution_->x[j]    = x[j];
01107           solution_->z_l[j]  = z_L[j];
01108           solution_->z_u[j]  = z_U[j];
01109      }
01110      solution_->g.resize(m_);
01111      solution_->lambda.resize(m_);
01112      for(i = 0; i < m_; i++)
01113      {    solution_->g[i]      = g[i];
01114           solution_->lambda[i] = lambda[i];
01115      }
01116      solution_->obj_value = obj_value;
01117      return;
01118 }
01119 
01120 // This routine is defined, but not yet used
01121 // (trying to figure out a problem with Ipopt-3.9.1 and dismod4).
01122 bool cppad_ipopt_nlp::intermediate_callback(
01123      Ipopt::AlgorithmMode              mode,
01124      Index                             iter, 
01125      Number                            obj_value,
01126      Number                            inf_pr, 
01127      Number                            inf_du,
01128      Number                            mu, 
01129      Number                            d_norm,
01130      Number                            regularization_size,
01131      Number                            alpha_du, 
01132      Number                            alpha_pr,
01133      Index                             ls_trials,
01134      const Ipopt::IpoptData*           ip_data,
01135      Ipopt::IpoptCalculatedQuantities* ip_cq
01136 )
01137 {
01138      // std::cout << "intermediate_callback called" << std::endl;
01139      return true;
01140 }
01141 
01142 // ---------------------------------------------------------------------------
01143 } // end namespace cppad_ipopt
01144 // ---------------------------------------------------------------------------
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines