/home/coin/SVN-release/OS-2.0.0/OS/CoinAllExamples/CppAD/ipopt_cppad_nlp.cpp

Go to the documentation of this file.
00001 /* --------------------------------------------------------------------------
00002 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-08 Bradley M. Bell
00003 
00004 CppAD is distributed under multiple licenses. This distribution is under
00005 the terms of the 
00006                     Common Public License Version 1.0.
00007 
00008 A copy of this license is included in the COPYING file of this distribution.
00009 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00010 -------------------------------------------------------------------------- */
00011 # include "ipopt_cppad_nlp.hpp"
00012 
00013 // define as 0 for false or 1 for true
00014 # define  CPPAD_NLP_TRACE 0
00015 
00016 # if CPPAD_NLP_TRACE
00017 # include <cstdio>
00018 # endif
00019 
00020 /* Constructor. */
00021 ipopt_cppad_nlp::ipopt_cppad_nlp(
00022         size_t                     n           , 
00023         size_t                     m           , 
00024         const NumberVector&        x_i         ,
00025         const NumberVector&        x_l         ,
00026         const NumberVector&        x_u         ,
00027         const NumberVector&        g_l         ,
00028         const NumberVector&        g_u         ,
00029         ipopt_cppad_fg_info*       fg_info     ,
00030         ipopt_cppad_solution*      solution    )
00031         : n_ ( n ),
00032           m_ ( m ),
00033           x_i_ ( x_i ),
00034           x_l_ ( x_l ),
00035           x_u_ ( x_u ),
00036           g_l_ ( g_l ),
00037           g_u_ ( g_u ),
00038           fg_info_ ( fg_info ) ,
00039           solution_ (solution)
00040 {       size_t i, j, k;
00041 
00042         // set information needed in ipopt_cppad_fg_info
00043         fg_info->set_n(n);
00044         fg_info->set_m(m);
00045 
00046         // get information from derived class version of fg_info
00047         K_ = fg_info->number_functions();
00048         L_.resize(K_);
00049         p_.resize(K_);
00050         q_.resize(K_);
00051         r_fun_.resize(K_);
00052         retape_.resize(K_);
00053         pattern_jac_r_.resize(K_);
00054         pattern_r_lag_.resize(K_);
00055         size_t max_p      = 0;
00056         size_t max_q      = 0;
00057         bool   retape_any = false;
00058         for(k = 0; k < K_; k++)
00059         {       L_[k]       = fg_info->number_terms(k);
00060                 p_[k]       = fg_info->range_size(k);
00061                 q_[k]       = fg_info->domain_size(k);
00062                 retape_[k]  = fg_info->retape(k);
00063                 max_p       = std::max(max_p, p_[k]);
00064                 max_q       = std::max(max_q, q_[k]);
00065                 retape_any |= retape_[k];
00066                 pattern_jac_r_[k].resize( p_[k] * q_[k] );
00067                 pattern_r_lag_[k].resize( q_[k] * q_[k] );
00068         }
00069         I_.resize(max_p);
00070         J_.resize(max_q);
00071 # ifndef NDEBUG
00072         size_t ell;
00073         // check for valid range and domain indices
00074         for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00075         {
00076                 for( i = 0; i < p_[k]; i++)
00077                         I_[i] = m+1; // an invalid range index
00078                 for( j = 0; j < q_[k]; j++)
00079                         J_[j] = n; // an invalid domain index
00080                 fg_info->index(k, ell, I_, J_); 
00081                 for( i = 0; i < p_[k]; i++) if( I_[i] > m )
00082                 {       std::cerr << "k=" << k << ", ell=" << ell 
00083                         << ", I[" << i << "]=" << I_[i] << std::endl;
00084                         CPPAD_ASSERT_KNOWN( I_[i] <= m,
00085                         "ipopt_cppad_nlp: invalid value in index vector I"
00086                         );
00087                 }
00088                 for( j = 0; j < q_[k]; j++) if( J_[j] >= n )
00089                 {       std::cerr << "k=" << k << ", ell=" << ell 
00090                         << ", J[" << j << "]=" << J_[j] << std::endl;
00091                         CPPAD_ASSERT_KNOWN( J_[j] < n,
00092                         "ipopt_cppad_nlp: invalid value in index vector J"
00093                         );
00094                 }
00095         }
00096 # endif
00097         for(k = 0; k < K_; k++) if( ! retape_[k] )
00098         {       // Record r_k (u): operation sequence does not depend on value 
00099                 // of u but record at initial value to make debugging easier.
00100                 fg_info->index(k, 0, I_, J_);
00101                 ADVector u_ad(q_[k]);
00102                 for(j = 0; j < q_[k]; j++)
00103                         u_ad[j] = x_i[ J_[j] ];
00104                 record_r_fun(
00105                         fg_info_, k, p_, q_, u_ad, // inputs
00106                         r_fun_                     // outputs
00107                 );
00108         }
00109         if ( retape_any )
00110         {       // true sparsity pattern valid for all x is unknown
00111                 BoolVector pattern_jac_fg((m+1) * n);
00112                 for(i = 0; i <= m; i++)
00113                 {       for(j = 0; j < n; j++)
00114                                 pattern_jac_fg[i * n + j] = true;
00115                 }
00116                 compute_index_jac_fg(m, n, pattern_jac_fg, index_jac_fg_);
00117                 //
00118                 BoolVector pattern_h_lag(n * n);
00119                 for(j = 0; j < n; j++)
00120                 {       for(k = 0; k < n; k++)
00121                                 pattern_h_lag[j * n + k] = true;
00122                 }
00123                 compute_index_h_lag(m, n, pattern_h_lag, index_h_lag_);
00124         }
00125         else
00126         {       // compute index map for Jacobian of fg
00127                 compute_index_jac_fg(
00128                 fg_info_, I_, J_, K_, L_, m_, n_, p_, q_, r_fun_,   // inputs
00129                 pattern_jac_r_, index_jac_fg_                       // outputs
00130                 );
00131 
00132                 // compute index map for Hessian of Lagragian
00133                 compute_index_h_lag(
00134                 fg_info_, I_, J_, K_, L_, m_, n_, p_, q_, r_fun_,   // inputs
00135                 pattern_r_lag_, index_h_lag_                        // outputs
00136                 );
00137         }
00138 
00139         // Compute Ipopt sparsity structure for Jacobian of g 
00140         compute_structure_jac_g(
00141                 index_jac_fg_, m, n,                  // inputs
00142                 nnz_jac_g_, iRow_jac_g_, jCol_jac_g_  // outputs
00143         );
00144 
00145         // Compute Ipopt sparsity structure for Hessian of Lagragian
00146         compute_structure_h_lag(
00147                 index_h_lag_, m, n,                   // inputs
00148                 nnz_h_lag_, iRow_h_lag_, jCol_h_lag_  // outputs
00149         );
00150 
00151 }
00152 
00153 // static member function that records operation sequence
00154 void ipopt_cppad_nlp::record_r_fun( 
00155         ipopt_cppad_fg_info*   fg_info  , 
00156         size_t                 k        ,
00157         SizeVector&            p        ,
00158         SizeVector&            q        ,
00159         ADVector&              u_ad     ,
00160         ADFunVector&           r_fun    )
00161 /*
00162 fg_info: input
00163 the ipopt_cppad_fg_info object that is used to information
00164 about the representation of fg(x).
00165 
00166 k: input
00167 index of the function r_k (u)
00168 
00169 p: input
00170 p[k] is number of components in the range space for r_k (u).
00171 
00172 q: input
00173 q[k] number of components in domain space for r_k (u).
00174 
00175 u_ad: input
00176 vector of independent variable values at which to record r_k (u).
00177 This is an input except for 
00178 the fact that its CppAD private data changes.
00179 
00180 r_fun: output 
00181 The CppAD operation sequence corresponding to the value of u_ad,
00182 and the algorithm used by fg_info, is stored in r_fun[k]. (Any operation
00183 seqeunce that was previously in r_fun[k] is deleted.)
00184 */
00185 {       CPPAD_ASSERT_UNKNOWN( u_ad.size() == size_t(q[k]) );
00186         // start the recording
00187         CppAD::Independent(u_ad);
00188         // vector of dependent variables during function recording
00189         ADVector r_ad = fg_info->eval_r(k, u_ad);
00190         CPPAD_ASSERT_KNOWN( r_ad.size() == p[k] ,
00191                 "ipopt_cppad_nlp: eval_r return value size not equal to p[k]."
00192         );
00193         // stop the recording and store operation sequence in r_fun
00194         r_fun[k].Dependent(u_ad, r_ad);
00195 }
00196 
00197 // static member function that computes CppAD sparsity pattern for 
00198 // Jacobian of fg
00199 void ipopt_cppad_nlp::compute_index_jac_fg(
00200         ipopt_cppad_fg_info*  fg_info        , 
00201         SizeVector&           I              ,
00202         SizeVector&           J              ,
00203         size_t                K              ,
00204         SizeVector&           L              ,
00205         size_t                m              ,
00206         size_t                n              ,
00207         SizeVector&           p              ,
00208         SizeVector&           q              ,
00209         ADFunVector&          r_fun          ,
00210         BoolVectorVector&     pattern_jac_r  ,
00211         IndexMap&             index_jac_fg   )
00212 /*
00213 fg_info: input
00214 the ipopt_cppad_fg_info object that is used to compute 
00215 information about the representation of fg(x).
00216 
00217 I: scratch space
00218 Vector of length greater than or equal p[k] for all k.
00219 
00220 J: scratch space
00221 Vector of length greatern than or equal q[k] for all k.
00222 
00223 K: input
00224 Number of functions in the representation for fg(x) in terms of r_k (u).
00225 
00226 L: input
00227 L[k] is the number of terms, in the representaiton for fg(x), that use r_k (u).
00228 
00229 m: input
00230 The number of components in the constraint function g.
00231 
00232 n: input
00233 Number of indpendent variables.
00234 
00235 p: input
00236 for k = 0 , ... , K-1,
00237 p[k] is dimension of the range space for r_k (u).
00238 
00239 q: input
00240 for k = 0 , ... , K-1,
00241 p[k] is dimension of the domain space for r_k (u).
00242 
00243 r_fun: input
00244 for k = 0 , ... , K-1,
00245 r_fun[k] is the CppAD function object that is used to compute 
00246 the sparsity patterns for r_k (u).
00247 The state of these objects actually changes because some forward mode
00248 routines are used.
00249 The operation sequence correspopnding to these object do not change.
00250 
00251 pattern_jac_r: output
00252 For k = 0 , ... , K-1,
00253 on input pattern_jac_r[k], this must be a vector of length p[k] * q[k].
00254 On output it is the CppAD sparsity pattern for the Jacobian of r_k (u).
00255 
00256 index_jac_fg:
00257 On input, this is empty; i.e., index_jac_fg.size() == 0.
00258 On output, it is the index mapping from (i, j) in the Jacobian of fg
00259 to the corresponding values array index in Ipopt. 
00260 Furthermore, if index_jac_fg[i].find(j) == index_jac_fg[i].end(),
00261 then either i = 0 or the (i, j) entry in the Jacobian of fg is always zero.
00262 */
00263 {
00264         size_t i, j, k, ell, ir, ifg;
00265 
00266         for(k = 0; k < K; k++)
00267         {       // compute sparsity pattern for r_k (u)
00268                 CPPAD_ASSERT_UNKNOWN( pattern_jac_r[k].size() == p[k] * q[k] );
00269                 CPPAD_ASSERT_UNKNOWN( r_fun[k].Range() == p[k] );
00270                 CPPAD_ASSERT_UNKNOWN( r_fun[k].Domain() == q[k] );
00271                 if( q[k] < p[k] )
00272                 {       // use forward mode
00273                         BoolVector pattern_domain(q[k] * q[k]);
00274                         for(i = 0; i < q[k]; i++)
00275                         {       for(j = 0; j < q[k]; j++) 
00276                                         pattern_domain[i * q[k] + j] = false;
00277                                 pattern_domain[i * q[k] + i] = true;
00278                         }
00279                         pattern_jac_r[k] = 
00280                                 r_fun[k].ForSparseJac(q[k], pattern_domain);
00281                 }
00282                 else
00283                 {       // use reverse mode
00284                         BoolVector pattern_range(p[k] * p[k]);
00285                         for(i = 0; i < p[k]; i++)
00286                         {       for(j = 0; j < p[k]; j++) 
00287                                         pattern_range[i * p[k] + j] = false;
00288                                 pattern_range[i * p[k] + i] = true;
00289                         }
00290                         pattern_jac_r[k] = 
00291                                 r_fun[k].RevSparseJac(p[k], pattern_range);
00292                 }
00293         }
00294 
00295         // now compute pattern for fg
00296         BoolVector pattern_jac_fg((m+1) * n);
00297         j = (m+1) * n;
00298         while(j--)
00299                 pattern_jac_fg[j] = false;
00300         for(k = 0; k < K; k++) for(ell = 0; ell < L[k]; ell++)
00301         {       fg_info->index(k, ell, I, J);   
00302                 for(i = 0; i < p[k]; i++)
00303                 {       for(j = 0; j < q[k]; j++)
00304                         {       ir  = i * q[k] + j;
00305                                 ifg = I[i] * n + J[j];
00306                                 pattern_jac_fg[ifg] = ( pattern_jac_fg[ifg] 
00307                                                     | pattern_jac_r[k][ir]    );
00308                         }
00309                 }
00310         }
00311         compute_index_jac_fg(m, n, pattern_jac_fg, index_jac_fg);
00312 }
00313 
00314 // static member function that computes index map from array indices
00315 // for Jacobian of fg
00316 void ipopt_cppad_nlp::compute_index_jac_fg(
00317         size_t                m              ,
00318         size_t                n              ,
00319         const BoolVector&     pattern_jac_fg ,
00320         IndexMap&             index_jac_fg
00321 )
00322 /*
00323 m: input
00324 The number of components in the constraint function g.
00325 
00326 n: input
00327 Number of indpendent variables.
00328 
00329 pattern_jac_fg:
00330 The CppAD sparsity pattern for the Jacobian of fg(x).
00331 
00332 index_jac_fg:
00333 On input, this is empty; i.e., index_jac_fg.size() == 0.
00334 On output, it is the index mapping from (i, j) in the Jacobian of fg
00335 to the corresponding values array index in Ipopt. 
00336 Furthermore, if index_jac_fg[i].find(j) == index_jac_fg[i].end(),
00337 then either i = 0 or the (i, j) entry in the Jacobian of fg is always zero.
00338 */
00339 {       CPPAD_ASSERT_UNKNOWN( index_jac_fg.size() == 0 );
00340         index_jac_fg.resize(m+1);
00341         size_t i, j, l = 0;
00342         for(i = 1; i <= m; i++)
00343         {       for(j = 0; j < n; j++)
00344                 {       if( pattern_jac_fg[ i * n + j ] )
00345                                 index_jac_fg[i][j] = l++;
00346                 }
00347         }
00348 }
00349 
00350 
00351 // static member function that computes CppAD sparsity pattern for 
00352 // Hessian of Lagragian
00353 void ipopt_cppad_nlp::compute_index_h_lag(
00354         ipopt_cppad_fg_info  *fg_info        , 
00355         SizeVector&           I              ,
00356         SizeVector&           J              ,
00357         size_t                K              ,
00358         SizeVector&           L              ,
00359         size_t                m              ,
00360         size_t                n              ,
00361         SizeVector&           p              ,
00362         SizeVector&           q              ,
00363         ADFunVector&          r_fun          ,
00364         BoolVectorVector&     pattern_r_lag  ,
00365         IndexMap&             index_h_lag    )
00366 /*
00367 fg_info: input
00368 the ipopt_cppad_fg_info object that is used to compute 
00369 information about the representation of fg(x).
00370 
00371 I: scratch space
00372 Vector of length greater than or equal p[k] for all k.
00373 
00374 J: scratch space
00375 Vector of length greatern than or equal q[k] for all k.
00376 
00377 K: input
00378 Number of functions in the representation for fg(x) in terms of r_k (u).
00379 
00380 L: input
00381 L[k] is the number of terms, in the representaiton for fg(x), that use r_k (u).
00382 
00383 m: input
00384 The number of components in the constraint function g.
00385 
00386 n: input
00387 Number of indpendent variables.
00388 
00389 p: input
00390 for k = 0 , ... , K-1,
00391 p[k] is dimension of the range space for r_k (u).
00392 
00393 q: input
00394 for k = 0 , ... , K-1,
00395 p[k] is dimension of the domain space for r_k (u).
00396 
00397 r_fun: input
00398 for k = 0 , ... , K-1,
00399 r_fun[k] is the CppAD function object that is used to compute 
00400 the sparsity patterns for r_k (u).
00401 The state of these objects actually changes because some forward mode
00402 routines are used.
00403 The operation sequence correspopnding to these object do not change.
00404 
00405 pattern_r_lag: output
00406 For k = 0 , ... , K-1,
00407 on input pattern_r_lag[k], this must be a vector of length q[k] * q[k].
00408 On output it is the CppAD sparsity pattern for the Hessian of 
00409 a Lagragian that sums components of r_k (u).
00410 
00411 index_h_lag:
00412 On input, this is empty; i.e., index_h_lag.size() == 0.
00413 On output, it is the index mapping from (i, j) in the Hessian of the Lagragian
00414 to the corresponding values array index in Ipopt. 
00415 Furthermore, if index_h_lag[i].find(j) == index_h_lag[i].end(),
00416 then either i < j or the (i, j) entry in the Hessian of the Lagragian is
00417 always zero.
00418 */
00419 {
00420         size_t i, j, k, ell;
00421 
00422         for(k = 0; k < K; k++)
00423         {       CPPAD_ASSERT_UNKNOWN( pattern_r_lag[k].size() == q[k] * q[k] );
00424                 CPPAD_ASSERT_UNKNOWN( r_fun[k].Range() == p[k] );
00425                 CPPAD_ASSERT_UNKNOWN( r_fun[k].Domain() == q[k] );
00426 
00427                 BoolVector pattern_domain(q[k] * q[k]);
00428                 BoolVector pattern_ones(p[k]);
00429 
00430                 for(i = 0; i < q[k]; i++)
00431                 {       for(j = 0; j < q[k]; j++) 
00432                                 pattern_domain[i * q[k] + j] = false;
00433                         pattern_domain[i * q[k] + i] = true;
00434                 }
00435                 r_fun[k].ForSparseJac(q[k], pattern_domain);
00436                 for(i = 0; i < p[k]; i++)
00437                         pattern_ones[i] = true;
00438                 pattern_r_lag[k] = r_fun[k].RevSparseHes(q[k], pattern_ones);
00439         }
00440 
00441         // now compute pattern for fg
00442         BoolVector pattern_h_lag(n * n);
00443         j = (n * n);
00444         while(j--)
00445                 pattern_h_lag[j] = false;
00446         for(k = 0; k < K; k++) for(ell = 0; ell < L[k]; ell++)
00447         {       fg_info->index(k, ell, I, J);
00448                 for(i = 0; i < q[k]; i++)
00449                 {       for(j = 0; j < q[k]; j++)
00450                         {       size_t ir, ifg;
00451                                 ir  = i * q[k] + j;
00452                                 ifg = J[i] * n + J[j];
00453                                 pattern_h_lag[ifg] = ( pattern_h_lag[ifg] 
00454                                                     | pattern_r_lag[k][ir]    );
00455                         }
00456                 }
00457         }
00458         compute_index_h_lag(m, n, pattern_h_lag, index_h_lag);
00459 }
00460 
00461 // static member function that computes index map from array indices
00462 // in Hessian of Lagragian
00463 void ipopt_cppad_nlp::compute_index_h_lag(
00464         size_t                m              ,
00465         size_t                n              ,
00466         const BoolVector&     pattern_h_lag  ,
00467         IndexMap&             index_h_lag
00468 )
00469 /*
00470 m: input
00471 The number of components in the constraint function g.
00472 
00473 n: input
00474 Number of indpendent variables.
00475 
00476 pattern_h_lag:
00477 The CppAD sparsity pattern for the Hessian of the Lagragian.
00478 
00479 index_h_lag:
00480 On input, this is empty; i.e., index_h_lag.size() == 0.
00481 On output, it is the index mapping from (i, j) in the Hessian of the Lagragian
00482 to the corresponding values array index in Ipopt. 
00483 Furthermore, if index_h_lag[i].find(j) == index_h_lag[i].end(),
00484 then either i < j or the (i, j) entry in the Hessian of the Lagragian is
00485 always zero.
00486 */
00487 {       CPPAD_ASSERT_UNKNOWN( index_h_lag.size() == 0 );
00488         index_h_lag.resize(n);
00489         size_t i, j, l = 0;
00490         for(i = 0; i < n; i++)
00491         {       for(j = 0; j <= i; j++)
00492                 {       if( pattern_h_lag[ i * n + j ] )
00493                                 index_h_lag[i][j] = l++;
00494                 }
00495         }
00496 }
00497 
00498 // static member function that computes the Ipopt sparsity structure for 
00499 // Jacobian of g
00500 void ipopt_cppad_nlp::compute_structure_jac_g(
00501         IndexMap&         index_jac_fg   , // const does not work
00502         size_t            m              ,
00503         size_t            n              ,
00504         size_t&           nnz_jac_g      ,
00505         SizeVector&       iRow_jac_g     ,
00506         SizeVector&       jCol_jac_g     )
00507 /*
00508 index_jac_fg:
00509 is the index mapping from (i, j) in the Jacobian of fg
00510 to the corresponding values array index in Ipopt. 
00511 If index_jac_fg[i].find(j) == index_jac_fg[i].end(),
00512 then either i = 0 or the (i, j) entry in the Jacobian of fg is always zero.
00513 
00514 m: input
00515 The number of components in the constraint function g.
00516 
00517 n: input
00518 Number of indpendent variables.
00519 
00520 nnz_jac_g: output
00521 The number of possibly non-zero entries in the Jacobian of g.
00522 
00523 iRow_jac_g: output
00524 The input size of this vector does not matter.
00525 On output it has size nnz_jac_g.
00526 It specifies the C row index (i.e. base one) 
00527 corresponding to each non-zero entry in the Jacobian of g.
00528 
00529 jCol_jac_g: output
00530 The input size of this vector does not matter.
00531 On output it has size nnz_jac_g.
00532 It specifies the C column index (i.e. base one) 
00533 corresponding to each non-zero entry in the Jacobian of g.
00534 */
00535 {       size_t i, j, l;
00536         std::map<size_t,size_t>::iterator index_ij;
00537 
00538         nnz_jac_g = 0;
00539         for(i = 1; i <= m; i++)
00540         {       for(j = 0; j < n; j++)
00541                 {       index_ij = index_jac_fg[i].find(j);
00542                         if( index_ij != index_jac_fg[i].end() )
00543                                 ++nnz_jac_g;
00544                 }
00545         }
00546         iRow_jac_g.resize( nnz_jac_g );
00547         jCol_jac_g.resize( nnz_jac_g );
00548         l = 0;
00549         for(i = 1; i <= m; i++)
00550         {       for(j = 0; j < n; j++)
00551                 {       index_ij = index_jac_fg[i].find(j);
00552                         if( index_ij != index_jac_fg[i].end() )
00553                         {       iRow_jac_g[l] = i - 1;
00554                                 jCol_jac_g[l] = j;
00555                                 l++;
00556                         }
00557                 }
00558         }
00559 }
00560 
00561 // static member function that computes the Ipopt sparsity structure for 
00562 // Hessian of Lagragian
00563 void ipopt_cppad_nlp::compute_structure_h_lag(
00564         IndexMap&         index_h_lag    , // const does not work
00565         size_t             m             ,
00566         size_t             n             ,
00567         size_t&           nnz_h_lag      ,
00568         SizeVector&       iRow_h_lag     ,
00569         SizeVector&       jCol_h_lag     )
00570 /*
00571 index_h_lag:
00572 is the index mapping from (i, j) in the Hessian of the Lagragian
00573 to the corresponding values array index in Ipopt. 
00574 If index_h_lag[i].find(j) == index_h_lag[i].end(),
00575 then either i < j or the (i, j) entry in the Hessian of the Lagragian is
00576 always zero.
00577 
00578 m: input
00579 The number of components in the constraint function g.
00580 
00581 n: input
00582 Number of indpendent variables.
00583 
00584 nnz_h_lag: output
00585 The number of possibly non-zero entries in the Hessian of the Lagragian.
00586 
00587 iRow_h_lag: output
00588 The input size of this vector does not matter.
00589 On output it has size nnz_h_lag.
00590 It specifies the C row index (i.e. base one) 
00591 corresponding to each non-zero entry in the Hessian of the Lagragian.
00592 
00593 jCol_h_lag: output
00594 The input size of this vector does not matter.
00595 On output it has size nnz_h_lag.
00596 It specifies the C column index (i.e. base one) 
00597 corresponding to each non-zero entry in the Hessian of the Lagragian.
00598 */
00599 {       size_t i, j, l;
00600         std::map<size_t,size_t>::iterator index_ij;
00601 
00602         nnz_h_lag = 0;
00603         for(i = 0; i < n; i++)
00604         {       for(j = 0; j <= i; j++)
00605                 {       index_ij = index_h_lag[i].find(j);
00606                         if( index_ij != index_h_lag[i].end() )
00607                                 ++nnz_h_lag;
00608                 }
00609         }
00610         iRow_h_lag.resize( nnz_h_lag );
00611         jCol_h_lag.resize( nnz_h_lag );
00612         l = 0;
00613         for(i = 0; i < n; i++)
00614         {       for(j = 0; j <= i; j++)
00615                 {       index_ij = index_h_lag[i].find(j);
00616                         if( index_ij != index_h_lag[i].end() )
00617                         {       iRow_h_lag[l] = i;
00618                                 jCol_h_lag[l] = j;
00619                                 l++;
00620                         }
00621                 }
00622         }
00623 }
00624 
00625 ipopt_cppad_nlp::~ipopt_cppad_nlp()
00626 {}
00627 
00628 bool ipopt_cppad_nlp::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00629                          Index& nnz_h_lag, IndexStyleEnum& index_style)
00630 {
00631         n = n_;
00632         m = m_;
00633         nnz_jac_g = nnz_jac_g_;
00634         nnz_h_lag = nnz_h_lag_;
00635 
00636         // use the fortran index style for row/col entries
00637         index_style = C_STYLE;
00638 
00639         return true;
00640 }
00641 
00642 bool ipopt_cppad_nlp::get_bounds_info(Index n, Number* x_l, Number* x_u,
00643                             Index m, Number* g_l, Number* g_u)
00644 {       size_t i, j;
00645         // here, the n and m we gave IPOPT in get_nlp_info are passed back 
00646         CPPAD_ASSERT_UNKNOWN(size_t(n) == n_);
00647         CPPAD_ASSERT_UNKNOWN(size_t(m) == m_);
00648 
00649         // pass back bounds
00650         for(j = 0; j < n_; j++)
00651         {       x_l[j] = x_l_[j];
00652                 x_u[j] = x_u_[j];
00653         }
00654         for(i = 0; i < m_; i++)
00655         {       g_l[i] = g_l_[i];
00656                 g_u[i] = g_u_[i];
00657         }
00658         
00659         return true;
00660 }
00661 
00662 bool ipopt_cppad_nlp::get_starting_point(Index n, bool init_x, Number* x,
00663                                bool init_z, Number* z_L, Number* z_U,
00664                                Index m, bool init_lambda,
00665                                Number* lambda)
00666 {       size_t j;
00667 
00668         CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00669         CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ );
00670         CPPAD_ASSERT_UNKNOWN(init_x == true);
00671         CPPAD_ASSERT_UNKNOWN(init_z == false);
00672         CPPAD_ASSERT_UNKNOWN(init_lambda == false);
00673 
00674         for(j = 0; j < n_; j++)
00675                 x[j] = x_i_[j];
00676 
00677         return true;
00678 }
00679 
00680 bool ipopt_cppad_nlp::eval_f(
00681         Index n, const Number* x, bool new_x, Number& obj_value
00682 )
00683 {       // This routine is not efficient. It would be better if users version
00684         // of the function was templated so that we could use double to 
00685         // calculate function values instead of retaping (when necessary).
00686         CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00687 
00688         size_t iobj, j, k, ell;
00689 
00690         // initialize summation
00691         obj_value = 0.;
00692 
00693         for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00694         {       fg_info_->index(k, ell, I_, J_);
00695                 for(iobj = 0; iobj < p_[k]; iobj++) if( I_[iobj] == 0 )
00696                 {       if( (new_x || K_ > 1)  && retape_[k] )
00697                         {       // Record r_k for value of u corresponding to x
00698                                 ADVector u_ad(q_[k]);
00699                                 for(j = 0; j < q_[k]; j++)
00700                                 {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00701                                         u_ad[j] = x[ J_[j] ];
00702                                 }
00703                                 record_r_fun(
00704                                         fg_info_, k, p_, q_, u_ad,  // inputs
00705                                         r_fun_                      // outputs
00706                                 );
00707                         }
00708                         NumberVector u(q_[k]);
00709                         NumberVector r(p_[k]);
00710                         for(j = 0; j < q_[k]; j++)
00711                         {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00712                                 u[j]   = x[ J_[j] ];
00713                         }
00714                         r          = r_fun_[k].Forward(0, u);
00715                         obj_value += r[iobj];
00716                 }
00717         }
00718 # if CPPAD_NLP_TRACE
00719         using std::printf;
00720         for(j = 0; j < n_; j++)
00721                 printf("ipopt_cppad_nlp::eval_f::x[%d] = %20.14g\n", j, x[j]);
00722         printf("ipopt_cppad_nlp::eval_f::obj_value = %20.14g\n", obj_value);
00723 # endif
00724         return true;
00725 }
00726 
00727 bool ipopt_cppad_nlp::eval_grad_f(
00728         Index n, const Number* x, bool new_x, Number* grad_f
00729 )
00730 {       CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00731 
00732         size_t iobj, i, j, k, ell;
00733 
00734         // initialize summation
00735         for(j = 0; j < n_; j++)
00736                 grad_f[j] = 0.;
00737 
00738         for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00739         {       fg_info_->index(k, ell, I_, J_);
00740                 for(iobj = 0; iobj < p_[k]; iobj++) if( I_[iobj] == 0 )
00741                 {       if( (new_x || K_ > 1)  && retape_[k] )
00742                         {       // Record r_k for value of u corresponding to x
00743                                 // If we has stored all recordings in f_eval
00744                                 // we might not need to do this over again.
00745                                 ADVector u_ad(q_[k]);
00746                                 for(j = 0; j < q_[k]; j++)
00747                                 {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00748                                         u_ad[j] = x[ J_[j] ];
00749                                 }
00750                                 record_r_fun(
00751                                         fg_info_, k, p_, q_, u_ad,  // inputs
00752                                         r_fun_                      // outputs
00753                                 );
00754                         }
00755                         NumberVector u(q_[k]);
00756                         NumberVector w(p_[k]);
00757                         NumberVector r_grad(q_[k]);
00758                         for(j = 0; j < q_[k]; j++)
00759                         {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00760                                 u[j]   = x[ J_[j] ];
00761                         }
00762                         r_fun_[k].Forward(0, u);
00763                         for(i = 0; i < p_[k]; i++)
00764                                 w[i] = 0.;
00765                         w[iobj]    = 1.;
00766                         r_grad     = r_fun_[k].Reverse(1, w);
00767                         for(j = 0; j < q_[k]; j++)
00768                         {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00769                                 grad_f[ J_[j] ]  += r_grad[j];
00770                         }
00771                 }
00772         }
00773 # if CPPAD_NLP_TRACE
00774         using std::printf;
00775         for(j = 0; j < n_; j++) printf(
00776         "ipopt_cppad_nlp::eval_grad_f::x[%d] = %20.14g\n", j, x[j]
00777         );
00778         for(j = 0; j < n_; j++) printf(
00779         "ipopt_cppad_nlp::eval_grad_f::grad_f[%d] = %20.14g\n", j, grad_f[j]
00780         );
00781 # endif
00782         return true;
00783 }
00784 
00785 bool ipopt_cppad_nlp::eval_g(
00786         Index n, const Number* x, bool new_x, Index m, Number* g
00787 )
00788 {       CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00789 
00790         size_t i, j, k, ell;
00791 
00792         // initialize summation
00793         for(i = 0; i < m_; i++)
00794                 g[i] = 0.;
00795 
00796         for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00797         {       fg_info_->index(k, ell, I_, J_);
00798                 if( (new_x || K_ > 1)  && retape_[k] )
00799                 {       // Record r_k for value of u corresponding to x
00800                         ADVector     u_ad(q_[k]);
00801                         for(j = 0; j < q_[k]; j++)
00802                         {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00803                                 u_ad[j] = x[ J_[j] ];
00804                         }
00805                         record_r_fun(
00806                                 fg_info_, k, p_, q_, u_ad,  // inputs
00807                                 r_fun_                      // outputs
00808                         );
00809                 }
00810                 NumberVector u(q_[k]);
00811                 NumberVector r(p_[k]);
00812                 for(j = 0; j < q_[k]; j++)
00813                 {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00814                         u[j]   = x[ J_[j] ];
00815                 }
00816                 r   = r_fun_[k].Forward(0, u);
00817                 for(i = 0; i < p_[k]; i++)
00818                 {       CPPAD_ASSERT_UNKNOWN( I_[i] <= m_ );
00819                         if( I_[i] >= 1 )
00820                                 g[ I_[i] - 1 ] += r[i];
00821                 }
00822         }
00823 # if CPPAD_NLP_TRACE
00824         using std::printf;
00825         for(j = 0; j < n_; j++)
00826                 printf("ipopt_cppad_nlp::eval_g::x[%d] = %20.14g\n", j, x[j]);
00827         for(i = 0; i < m_; i++)
00828                 printf("ipopt_cppad_nlp::eval_g::g[%d] = %20.14g\n", i, g[i]);
00829 # endif
00830         return true;
00831 }
00832 
00833 bool ipopt_cppad_nlp::eval_jac_g(Index n, const Number* x, bool new_x,
00834                        Index m, Index nele_jac, Index* iRow, Index *jCol,
00835                        Number* values)
00836 {       CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ );
00837         CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00838 
00839         size_t i, j, k, ell, l;
00840         std::map<size_t,size_t>::iterator index_ij;
00841 
00842 
00843         if (values == NULL) 
00844         {       for(k = 0; k < nnz_jac_g_; k++)
00845                 {       iRow[k] = iRow_jac_g_[k];
00846                         jCol[k] = jCol_jac_g_[k];
00847                 }
00848                 return true;
00849         }
00850 
00851         // initialize summation
00852         l = nnz_jac_g_;
00853         while(l--)
00854                 values[l] = 0.;
00855 
00856         for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00857         {       fg_info_->index(k, ell, I_, J_);
00858                 if( (new_x || K_ > 1)  && retape_[k] )
00859                 {       // Record r_k for value of u corresponding to x
00860                         ADVector     u_ad(q_[k]);
00861                         for(j = 0; j < q_[k]; j++)
00862                         {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00863                                 u_ad[j] = x[ J_[j] ];
00864                         }
00865                         record_r_fun(
00866                                 fg_info_, k, p_, q_, u_ad,  // inputs
00867                                 r_fun_                      // outputs
00868                         );
00869                 }
00870                 NumberVector u(q_[k]);
00871                 NumberVector jac_r(p_[k] * q_[k]);
00872                 for(j = 0; j < q_[k]; j++)
00873                 {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00874                         u[j]   = x[ J_[j] ];
00875                 }
00876                 if( retape_[k] )
00877                         jac_r = r_fun_[k].Jacobian(u);
00878                 else    jac_r = r_fun_[k].SparseJacobian(u, pattern_jac_r_[k]);
00879                 for(i = 0; i < p_[k]; i++) if( I_[i] != 0 )
00880                 {       CPPAD_ASSERT_UNKNOWN( I_[i] <= m_ );
00881                         for(j = 0; j < q_[k]; j++)
00882                         {       index_ij = index_jac_fg_[I_[i]].find(J_[j]);
00883                                 if( index_ij != index_jac_fg_[I_[i]].end() )
00884                                 {       l          = index_ij->second;
00885                                         values[l] += jac_r[i * q_[k] + j];
00886                                 }
00887                                 else    CPPAD_ASSERT_UNKNOWN(
00888                                         jac_r[i * q_[k] + j] == 0.
00889                                 );
00890                         }
00891                 }
00892         }
00893         return true;
00894 }
00895 
00896 bool ipopt_cppad_nlp::eval_h(Index n, const Number* x, bool new_x,
00897                    Number obj_factor, Index m, const Number* lambda,
00898                    bool new_lambda, Index nele_hess, Index* iRow,
00899                    Index* jCol, Number* values)
00900 {       CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ );
00901         CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00902 
00903         size_t i, j, k, ell, l;
00904         std::map<size_t,size_t>::iterator index_ij;
00905 
00906         if (values == NULL) 
00907         {       for(k = 0; k < nnz_h_lag_; k++)
00908                 {       iRow[k] = iRow_h_lag_[k];
00909                         jCol[k] = jCol_h_lag_[k];
00910                 }
00911                 return true;
00912         }
00913 
00914         // initialize summation
00915         l = nnz_h_lag_;
00916         while(l--)
00917                 values[l] = 0.;
00918 
00919         for(k = 0; k < K_; k++) for(ell = 0; ell < L_[k]; ell++)
00920         {       fg_info_->index(k, ell, I_, J_);
00921                 if( (new_x || K_ > 1)  && retape_[k] )
00922                 {       // Record r_k for value of u corresponding to x
00923                         ADVector     u_ad(q_[k]);
00924                         for(j = 0; j < q_[k]; j++)
00925                         {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00926                                 u_ad[j] = x[ J_[j] ];
00927                         }
00928                         record_r_fun(
00929                                 fg_info_, k, p_, q_, u_ad,  // inputs
00930                                 r_fun_                      // outputs
00931                         );
00932                 }
00933                 NumberVector w(p_[k]);
00934                 NumberVector r_hes(q_[k] * q_[k]);
00935                 NumberVector u(q_[k]);
00936                 for(j = 0; j < q_[k]; j++)
00937                 {       CPPAD_ASSERT_UNKNOWN( J_[j] < n_ );
00938                         u[j]   = x[ J_[j] ];
00939                 }
00940                 for(i = 0; i < p_[k]; i++)
00941                 {       CPPAD_ASSERT_UNKNOWN( I_[i] <= m_ );
00942                         if( I_[i] == 0 )
00943                                 w[i] = obj_factor;
00944                         else    w[i] = lambda[ I_[i] - 1 ];
00945                 }
00946                 if( retape_[k] )
00947                         r_hes = r_fun_[k].Hessian(u, w);
00948                 else    r_hes = 
00949                         r_fun_[k].SparseHessian(u, w, pattern_r_lag_[k]);
00950                 for(i = 0; i < q_[k]; i++)
00951                 {       for(j = 0; j < q_[k]; j++) if( J_[j] <= J_[i] ) 
00952                         {       index_ij = index_h_lag_[J_[i]].find(J_[j]);
00953                                 if( index_ij != index_h_lag_[J_[i]].end() )
00954                                 {       l          = index_ij->second;
00955                                         values[l] += r_hes[i * q_[k] + j];
00956                                 }
00957                                 else    CPPAD_ASSERT_UNKNOWN(
00958                                         r_hes[i * q_[k] + j] == 0.
00959                                 );
00960                         }
00961                 }
00962         }
00963         return true;
00964 }
00965 
00966 void ipopt_cppad_nlp::finalize_solution(
00967         Ipopt::SolverReturn               status    ,
00968         Index                             n         , 
00969         const Number*                     x         , 
00970         const Number*                     z_L       , 
00971         const Number*                     z_U       ,
00972         Index                             m         , 
00973         const Number*                     g         , 
00974         const Number*                     lambda    ,
00975         Number                            obj_value ,
00976         const Ipopt::IpoptData*           ip_data   ,
00977         Ipopt::IpoptCalculatedQuantities* ip_cq
00978 )
00979 {       size_t i, j;
00980 
00981         CPPAD_ASSERT_UNKNOWN(size_t(n) == n_ );
00982         CPPAD_ASSERT_UNKNOWN(size_t(m) == m_ );
00983 
00984         switch(status)
00985         {       // convert status from Ipopt enum to ipopt_cppad_solution enum
00986                 case Ipopt::SUCCESS:
00987                 solution_->status = 
00988                         ipopt_cppad_solution::success;
00989                 break;
00990 
00991                 case Ipopt::MAXITER_EXCEEDED:
00992                 solution_->status = 
00993                         ipopt_cppad_solution::maxiter_exceeded;
00994                 break;
00995 
00996                 case Ipopt::STOP_AT_TINY_STEP:
00997                 solution_->status = 
00998                         ipopt_cppad_solution::stop_at_tiny_step;
00999                 break;
01000 
01001                 case Ipopt::STOP_AT_ACCEPTABLE_POINT:
01002                 solution_->status = 
01003                         ipopt_cppad_solution::stop_at_acceptable_point;
01004                 break;
01005 
01006                 case Ipopt::LOCAL_INFEASIBILITY:
01007                 solution_->status = 
01008                         ipopt_cppad_solution::local_infeasibility;
01009                 break;
01010 
01011                 case Ipopt::USER_REQUESTED_STOP:
01012                 solution_->status = 
01013                         ipopt_cppad_solution::user_requested_stop;
01014                 break;
01015 
01016                 case Ipopt::DIVERGING_ITERATES:
01017                 solution_->status = 
01018                         ipopt_cppad_solution::diverging_iterates;
01019                 break;
01020 
01021                 case Ipopt::RESTORATION_FAILURE:
01022                 solution_->status = 
01023                         ipopt_cppad_solution::restoration_failure;
01024                 break;
01025 
01026                 case Ipopt::ERROR_IN_STEP_COMPUTATION:
01027                 solution_->status = 
01028                         ipopt_cppad_solution::error_in_step_computation;
01029                 break;
01030 
01031                 case Ipopt::INVALID_NUMBER_DETECTED:
01032                 solution_->status = 
01033                         ipopt_cppad_solution::invalid_number_detected;
01034                 break;
01035 
01036                 case Ipopt::INTERNAL_ERROR:
01037                 solution_->status = 
01038                         ipopt_cppad_solution::internal_error;
01039                 break;
01040 
01041                 default:
01042                 solution_->status = 
01043                         ipopt_cppad_solution::unknown;
01044         }
01045 
01046         solution_->x.resize(n_);
01047         solution_->z_l.resize(n_);
01048         solution_->z_u.resize(n_);
01049         for(j = 0; j < n_; j++)
01050         {       solution_->x[j]    = x[j];
01051                 solution_->z_l[j]  = z_L[j];
01052                 solution_->z_u[j]  = z_U[j];
01053         }
01054         solution_->g.resize(m_);
01055         solution_->lambda.resize(m_);
01056         for(i = 0; i < m_; i++)
01057         {       solution_->g[i]      = g[i];
01058                 solution_->lambda[i] = lambda[i];
01059         }
01060         solution_->obj_value = obj_value;
01061         return;
01062 }

Generated on Mon Aug 3 03:02:21 2009 by  doxygen 1.4.7