CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_FORWARD_INCLUDED 00003 # define CPPAD_FORWARD_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell 00007 00008 CppAD is distributed under multiple licenses. This distribution is under 00009 the terms of the 00010 Eclipse Public License Version 1.0. 00011 00012 A copy of this license is included in the COPYING file of this distribution. 00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00014 -------------------------------------------------------------------------- */ 00015 00016 // documened after Forward but included here so easy to see 00017 # include <cppad/local/capacity_order.hpp> 00018 # include <cppad/local/num_skip.hpp> 00019 00020 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 00021 /*! 00022 \file forward.hpp 00023 User interface to forward mode computations 00024 */ 00025 00026 /*! 00027 Arbitrary order, one direction, forward mode Taylor coefficieints. 00028 00029 \tparam Base 00030 The type used during the forward mode computations; i.e., the corresponding 00031 recording of operations used the type AD<Base>. 00032 00033 \tparam Vector 00034 is a Simple Vector class with eleements of type Base. 00035 00036 \param q 00037 is the hightest order for this forward mode computation; i.e., 00038 after this calculation there will be <code>q+1</code> 00039 Taylor coefficients per variable. 00040 00041 \param xq 00042 contains Taylor coefficients for the independent variables. 00043 The size of xq must either be n or <code>(q+1)*n</code>, 00044 We define <code>p = q + 1 - xq.size()/n</code>. 00045 For <code>j = 0 , ... , n-1</code>, 00046 <code>k = p, ... , q</code>, are 00047 <code>xq[ (q+1-p)*j + k - p ]</code> 00048 is the k-th order coefficient for the j-th independent variable. 00049 00050 \param s 00051 Is the stream where output corresponding to PriOp operations will written. 00052 00053 \return 00054 contains Taylor coefficients for the independent variables. 00055 The size of the return value y is <code>m*(q+1-p)</code>. 00056 For <code>i = 0, ... , m-1</code>, 00057 <code>k = p, ..., q</code>, 00058 <code>y[(q+1-p)*i + (k-p)]</code> 00059 is the k-th order coefficient for the i-th dependent variable. 00060 00061 \par taylor_ 00062 The Taylor coefficients up to order p-1 are inputs 00063 and the coefficents from order p through q are outputs. 00064 Let <code>N = num_var_tape_</code>, and 00065 <code>C = cap_order_taylor_</code>. 00066 Note that for 00067 <code>i = 1 , ..., N-1</code>, 00068 <code>k = 0 , ..., q</code>, 00069 <code>taylor_[ C*i + k ]</code> 00070 is the k-th order cofficent, 00071 for the i-th varaible on the tape. 00072 (The first independent variable has index one on the tape 00073 and there is no variable with index zero.) 00074 */ 00075 00076 template <typename Base> 00077 template <typename Vector> 00078 Vector ADFun<Base>::Forward( 00079 size_t q , 00080 const Vector& xq , 00081 std::ostream& s ) 00082 { // temporary indices 00083 size_t i, j, k; 00084 00085 // number of independent variables 00086 size_t n = ind_taddr_.size(); 00087 00088 // number of dependent variables 00089 size_t m = dep_taddr_.size(); 00090 00091 // check Vector is Simple Vector class with Base type elements 00092 CheckSimpleVector<Base, Vector>(); 00093 00094 00095 CPPAD_ASSERT_KNOWN( 00096 size_t(xq.size()) == n || size_t(xq.size()) == n*(q+1), 00097 "Forward(q, xq): xq.size() is not equal n or n*(q+1)" 00098 ); 00099 00100 // lowest order we are computing 00101 size_t p = q + 1 - size_t(xq.size()) / n; 00102 CPPAD_ASSERT_UNKNOWN( p == 0 || p == q ); 00103 CPPAD_ASSERT_KNOWN( 00104 q <= num_order_taylor_ || p == 0, 00105 "Forward(q, xq): Number of Taylor coefficient orders stored in this" 00106 " ADFun\nis less than q and xq.size() != n*(q+1)." 00107 ); 00108 // does taylor_ need more orders or different number of directions 00109 if( cap_order_taylor_ <= q ) 00110 { if( p == 0 ) 00111 { // no need to copy old values during capacity_order 00112 num_order_taylor_ = 0; 00113 } 00114 size_t c = std::max(q + 1, cap_order_taylor_); 00115 capacity_order(c); 00116 } 00117 CPPAD_ASSERT_UNKNOWN( cap_order_taylor_ > q ); 00118 00119 // short hand notation for order capacity 00120 size_t C = cap_order_taylor_; 00121 00122 // set Taylor coefficients for independent variables 00123 for(j = 0; j < n; j++) 00124 { CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < num_var_tape_ ); 00125 00126 // ind_taddr_[j] is operator taddr for j-th independent variable 00127 CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == InvOp ); 00128 00129 if( p == q ) 00130 taylor_[ C * ind_taddr_[j] + q] = xq[j]; 00131 else 00132 { for(k = 0; k <= q; k++) 00133 taylor_[ C * ind_taddr_[j] + k] = xq[ (q+1)*j + k]; 00134 } 00135 } 00136 00137 // evaluate the derivatives 00138 CPPAD_ASSERT_UNKNOWN( cskip_op_.size() == play_.num_op_rec() ); 00139 CPPAD_ASSERT_UNKNOWN( load_op_.size() == play_.num_load_op_rec() ); 00140 if( q == 0 ) 00141 { 00142 compare_change_ = forward0sweep(s, true, 00143 n, num_var_tape_, &play_, C, 00144 taylor_.data(), cskip_op_.data(), load_op_ 00145 ); 00146 } 00147 else if( p == 0 ) 00148 { compare_change_ = forward1sweep(s, true, p, q, 00149 n, num_var_tape_, &play_, C, 00150 taylor_.data(), cskip_op_.data(), load_op_ 00151 ); 00152 } 00153 else 00154 { forward1sweep(s, true, p, q, 00155 n, num_var_tape_, &play_, C, 00156 taylor_.data(), cskip_op_.data(), load_op_ 00157 ); 00158 } 00159 00160 // return Taylor coefficients for dependent variables 00161 Vector yq; 00162 if( p == q ) 00163 { yq.resize(m); 00164 for(i = 0; i < m; i++) 00165 { CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_ ); 00166 yq[i] = taylor_[ C * dep_taddr_[i] + q]; 00167 } 00168 } 00169 else 00170 { yq.resize(m * (q+1) ); 00171 for(i = 0; i < m; i++) 00172 { for(k = 0; k <= q; k++) 00173 yq[ (q+1) * i + k] = 00174 taylor_[ C * dep_taddr_[i] + k ]; 00175 } 00176 } 00177 # ifndef NDEBUG 00178 if( check_for_nan_ ) 00179 { bool ok = true; 00180 if( p == 0 ) 00181 { for(i = 0; i < m; i++) 00182 { // Visual Studio 2012, CppAD required in front of isnan ? 00183 ok &= ! CppAD::isnan( yq[ (q+1) * i + 0 ] ); 00184 } 00185 } 00186 CPPAD_ASSERT_KNOWN(ok, 00187 "yq = f.Forward(q, xq): has a zero order Taylor coefficient " 00188 "with the value nan." 00189 ); 00190 if( 0 < q ) 00191 { for(i = 0; i < m; i++) 00192 { for(k = p; k <= q; k++) 00193 { // Studio 2012, CppAD required in front of isnan ? 00194 ok &= ! CppAD::isnan( yq[ (q+1-p)*i + k-p ] ); 00195 } 00196 } 00197 } 00198 CPPAD_ASSERT_KNOWN(ok, 00199 "yq = f.Forward(q, xq): has a non-zero order Taylor coefficient\n" 00200 "with the value nan (but zero order coefficients are not nan)." 00201 ); 00202 } 00203 # endif 00204 00205 // now we have q + 1 taylor_ coefficient orders per variable 00206 num_order_taylor_ = q + 1; 00207 00208 return yq; 00209 } 00210 00211 00212 } // END_CPPAD_NAMESPACE 00213 # endif