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