CppAD: A C++ Algorithmic Differentiation Package 20110419
hes_fg_map.cpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 /* --------------------------------------------------------------------------
00003 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 Bradley M. Bell
00004 
00005 CppAD is distributed under multiple licenses. This distribution is under
00006 the terms of the 
00007                     Common 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 "cppad_ipopt_nlp.hpp"
00013 # include "hes_fg_map.hpp"
00014 
00015 CPPAD_BEGIN_NAMESPACE
00016 /*!
00017 \file hes_fg_map.cpp
00018 \brief Creates a mapping between two representations for Hessian of fg. 
00019 */
00020 
00021 /*!
00022 Create mapping from CppAD to Ipopt sparse representations of Hessian of F(x).
00023 
00024 The functions 
00025 \f$ f : {\bf R}^n \rightarrow {\bf R} \f$ and
00026 \f$ g : {\bf R}^n \rightarrow {\bf R}^m \f$ are defined by the 
00027 \ref Users_Representation.
00028 We define the function 
00029 \f$ F : {\bf R}^n \rightarrow {\bf R} \f$ by
00030 \f[
00031         F(x) = \sum_{i=0}^m fg(x)_i
00032 \f]
00033 
00034 \param fg_info
00035 For <tt>k = 0 , ... , K-1</tt>,
00036 for <tt>ell = 0 , ... , L[k]</tt>,
00037 the function call 
00038 \verbatim
00039         fg_info->index(k, ell, I, J);   
00040 \endverbatim
00041 is made by \c hes_fg_map. 
00042 The values \c k and \c ell are inputs.
00043 The input size of \c I ( \c J ) 
00044 is greater than or equal <tt>p[k] ( q[k] )</tt>
00045 and this size is not changed.
00046 The input values of the elements of \c I and \c J are not specified.
00047 The output value of the elements of \c I define
00048 \f[
00049 I_{k, \ell} = ( {\rm I[0]} , \cdots , {\rm I[p[k]-1]} )
00050 \f]
00051 The output value of the elements of \c J define
00052 \f[
00053 J_{k, \ell} = ( {\rm J[0]} , \cdots , {\rm J[q[k]-1]} )
00054 \f]
00055 
00056 \param m
00057 is the dimension of the range space for \f$ g(x) \f$; i.e.,
00058 \f$ g(x) \in {\bf R}^m \f$.
00059 
00060 \param n
00061 is the dimension of the domain space for \f$ f(x) \f$ and \f$ g(x) \f$;
00062 i.e., \f$ x \in {\bf R}^n \f$.
00063 
00064 \param K
00065 is the number of functions \f$ r_k ( u ) \f$ used for the representation of 
00066 \f$ f(x) \f$ and \f$ g(x) \f$.
00067 
00068 \param L
00069 is a vector with size \c K.
00070 For <tt>k = 0 , ... , K-1, L[k]</tt>
00071 is the number of terms that use \f$ r_k (u) \f$
00072 in the representation of \f$ f(x) \f$ and \f$ g(x) \f$.
00073 
00074 \param p
00075 is a vector with size \c K.
00076 For <tt>k = 0 , ... , K-1, p[k]</tt>
00077 is dimension of the range space for \f$ r_k (u) \f$; i.e.,
00078 \f$ r_k (u) \in {\bf R}^{p(k)} \f$.
00079 
00080 \param q
00081 is a vector with size \c K.
00082 For <tt>k = 0 , ... , K-1, q[k]</tt>
00083 is dimension of the domain space for \f$ r_k (u) \f$; i.e.,
00084 \f$ u \in {\bf R}^{q(k)} \f$.
00085 
00086 \param pattern_hes_r
00087 is a vector with size \c K.
00088 For <tt>k = 0 , ... , K-1, pattern_jac_r[k]</tt>
00089 is a CppAD sparsity pattern for the Hessian of the function 
00090 \f[
00091         R(u) = \sum_{i=0}^{p[k]-1} r_k (u)_i
00092 \f]
00093 As such, <tt>pattern_hes_r[k].size() == q[k] * q[k]</tt>.
00094 
00095 \param I
00096 is a work vector of length greater than or equal <tt>p[k]</tt> for all \c k.
00097 The input and output value of its elements are unspecified.
00098 The size of \c I is not changed.
00099 
00100 \param J
00101 is a work vector of length greater than or equal <tt>q[k]</tt> for all \c k.
00102 The input and output value of its elements are unspecified.
00103 The size of \c J is not changed.
00104 
00105 \param index_hes_fg:
00106 On input, this is empty; i.e., <tt>index_jac_g.size() == 0</tt>.
00107 On output, it is the index mapping from \f$ (i, j) \f$ in the Jacobian of 
00108 \f$ g(x) \f$ to the corresponding index value used by Ipopt to represent
00109 the Jacobian.
00110 Furthermore, if <tt>index_jac_g[i].find(j) == index_jac_g[i].end()</tt>,
00111 then the \f$ (i, j)\f$ entry in the Jacobian of \f$ g(x) \f$ is always zero.
00112 */
00113 void hes_fg_map(
00114         cppad_ipopt_fg_info*  fg_info                                  , 
00115         size_t                                          m              ,
00116         size_t                                          n              ,
00117         size_t                                          K              ,
00118         const CppAD::vector<size_t>&                    L              ,
00119         const CppAD::vector<size_t>&                    p              ,
00120         const CppAD::vector<size_t>&                    q              ,
00121         const CppAD::vector<CppAD::vectorBool>&         pattern_hes_r  ,
00122         CppAD::vector<size_t>&                          I              ,
00123         CppAD::vector<size_t>&                          J              ,
00124         CppAD::vector< std::map<size_t,size_t> >&       index_hes_fg   )
00125 {
00126         using CppAD::vectorBool;
00127         size_t i, j, ij, k, ell;
00128 
00129         CPPAD_ASSERT_UNKNOWN( K == L.size() );
00130         CPPAD_ASSERT_UNKNOWN( K == p.size() );
00131         CPPAD_ASSERT_UNKNOWN( K == q.size() );
00132         CPPAD_ASSERT_UNKNOWN( K == pattern_hes_r.size() );
00133 # ifndef NDEBUG
00134         for(k = 0; k < K; k++)
00135         {        CPPAD_ASSERT_UNKNOWN( p[k] <= I.size() );
00136                  CPPAD_ASSERT_UNKNOWN( q[k] <= J.size() );
00137                  CPPAD_ASSERT_UNKNOWN( q[k]*q[k] == pattern_hes_r[k].size() );
00138         }
00139 # endif
00140         
00141         // Now compute pattern for fg
00142         // (use standard set representation because can be huge).
00143         CppAD::vector< std::set<size_t> > pattern_hes_fg(n);
00144         for(k = 0; k < K; k++) for(ell = 0; ell < L[k]; ell++)
00145         {       fg_info->index(k, ell, I, J);   
00146                 for(i = 0; i < q[k]; i++) 
00147                 {       for(j = 0; j < q[k]; j++)
00148                         {       ij  = i * q[k] + j;
00149                                 if( pattern_hes_r[k][ij] )
00150                                         pattern_hes_fg[J[i]].insert(J[j]);
00151                         }
00152                 }
00153         }
00154 
00155         // Now compute the mapping from (i, j) in the Hessian of fg to the 
00156         // corresponding index value used by Ipopt to represent the Hessian.
00157         CPPAD_ASSERT_UNKNOWN( index_hes_fg.size() == 0 );
00158         index_hes_fg.resize(n);
00159         std::set<size_t>::const_iterator itr;
00160         ell = 0;
00161         for(i = 0; i < n; i++)
00162         {       for(    itr = pattern_hes_fg[i].begin(); 
00163                         itr != pattern_hes_fg[i].end(); itr++)
00164                 {
00165                         index_hes_fg[i][*itr] = ell++;
00166                 }
00167         }
00168         return;
00169 }
00170 
00171 CPPAD_END_NAMESPACE