CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_FUN_CONSTRUCT_INCLUDED 00003 # define CPPAD_FUN_CONSTRUCT_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 Bradley M. Bell 00007 00008 CppAD is distributed under multiple licenses. This distribution is under 00009 the terms of the 00010 Common 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 $begin FunConstruct$$ 00017 $spell 00018 Jac 00019 bool 00020 taylor_ 00021 var 00022 ADvector 00023 const 00024 Jacobian 00025 $$ 00026 00027 $spell 00028 $$ 00029 00030 $section Construct an ADFun Object and Stop Recording$$ 00031 00032 $index ADFun, construct$$ 00033 $index construct, ADFun$$ 00034 $index tape, stop recording$$ 00035 $index recording, stop tape$$ 00036 00037 $head Syntax$$ 00038 $codei%ADFun<%Base%> %f%, %g% 00039 %$$ 00040 $codei%ADFun<%Base%> %f%(%x%, %y%) 00041 %$$ 00042 $icode%g% = %f% 00043 %$$ 00044 00045 00046 $head Purpose$$ 00047 The $codei%AD<%Base%>%$$ object $icode f$$ can 00048 store an AD of $icode Base$$ 00049 $xref/glossary/Operation/Sequence/operation sequence/1/$$. 00050 It can then be used to calculate derivatives of the corresponding 00051 $xref/glossary/AD Function/AD function/$$ 00052 $latex \[ 00053 F : B^n \rightarrow B^m 00054 \] $$ 00055 where $latex B$$ is the space corresponding to objects of type $icode Base$$. 00056 00057 $head x$$ 00058 If the argument $icode x$$ is present, it has prototype 00059 $codei% 00060 const %VectorAD% &%x% 00061 %$$ 00062 It must be the vector argument in the previous call to 00063 $cref/Independent/$$. 00064 Neither its size, or any of its values, are allowed to change 00065 between calling 00066 $codei% 00067 Independent(%x%) 00068 %$$ 00069 and 00070 $codei% 00071 ADFun<%Base%> %f%(%x%, %y%) 00072 %$$. 00073 00074 $head y$$ 00075 If the argument $icode y$$ is present, it has prototype 00076 $codei% 00077 const %VectorAD% &%y% 00078 %$$ 00079 The sequence of operations that map $icode x$$ 00080 to $icode y$$ are stored in the ADFun object $icode f$$. 00081 00082 $head VectorAD$$ 00083 The type $icode VectorAD$$ must be a $cref/SimpleVector/$$ class with 00084 $cref/elements of type/SimpleVector/Elements of Specified Type/$$ 00085 $codei%AD<%Base%>%$$. 00086 The routine $cref/CheckSimpleVector/$$ will generate an error message 00087 if this is not the case. 00088 00089 $head Default Constructor$$ 00090 $index default, ADFun constructor$$ 00091 $index ADFun, default constructor$$ 00092 $index constructor, ADFun constructor$$ 00093 The default constructor 00094 $codei% 00095 ADFun<%Base%> %f% 00096 %$$ 00097 creates an 00098 $codei%AD<%Base%>%$$ object with no corresponding operation sequence; i.e., 00099 $codei% 00100 %f%.size_var() 00101 %$$ 00102 returns the value zero (see $xref/seq_property/size_var/size_var/$$). 00103 00104 $head Sequence Constructor$$ 00105 $index sequence, ADFun constructor$$ 00106 $index ADFun, sequence constructor$$ 00107 $index constructor, ADFun sequence$$ 00108 The default constructor 00109 The sequence constructor 00110 $codei% 00111 ADFun<%Base%> %f%(%x%, %y%) 00112 %$$ 00113 creates the $codei%AD<%Base%>%$$ object $icode f$$, 00114 stops the recording of AD of $icode Base$$ operations 00115 corresponding to the call 00116 $codei% 00117 Independent(%x%) 00118 %$$ 00119 and stores the corresponding operation sequence in the object $icode f$$. 00120 It then stores the first order taylor_ coefficients 00121 (corresponding to the value of $icode x$$) in $icode f$$. 00122 This is equivalent to the following steps using the default constructor: 00123 00124 $list number$$ 00125 Create $icode f$$ with the default constructor 00126 $codei% 00127 ADFun<%Base%> %f%; 00128 %$$ 00129 $lnext 00130 Stop the tape and storing the operation sequence using 00131 $codei% 00132 %f%.Dependent(%x%, %y%); 00133 %$$ 00134 (see $xref/Dependent/$$). 00135 $lnext 00136 Calculating the first order taylor_ coefficients for all 00137 the variables in the operation sequence using 00138 $codei% 00139 %f%.Forward(%p%, %x_p%) 00140 %$$ 00141 with $icode p$$ equal to zero and the elements of $icode x_p$$ 00142 equal to the corresponding elements of $icode x$$ 00143 (see $xref/Forward/$$). 00144 $lend 00145 00146 $head Copy Constructor$$ 00147 $index copy, ADFun constructor$$ 00148 $index ADFun, copy constructor$$ 00149 $index constructor, ADFun copy$$ 00150 It is an error to attempt to use the $code%ADFun<%Base%>%$$ copy constructor; 00151 i.e., the following syntax is not allowed: 00152 $codei% 00153 ADFun<%Base%> g(f) 00154 %$$ 00155 where $icode f$$ is an $code%ADFun<%Base%>%$$ object. 00156 Use its $cref/default constructor/FunConstruct/Default Constructor/$$ instead 00157 and its assignment operator. 00158 00159 $head Assignment Operator$$ 00160 $index ADFun, assignment operator$$ 00161 $index assignment, ADFun operator$$ 00162 $index operator, ADFun assignment$$ 00163 The $codei%ADFun<%Base%>%$$ assignment operation 00164 $codei% 00165 %g% = %f% 00166 %$$. 00167 makes a copy of the operation sequence currently stored in $icode f$$ 00168 in the object $icode g$$. 00169 The object $icode f$$ is not affected by this operation and 00170 can be $code const$$. 00171 Any operation sequence or other information in $icode g$$ is lost. 00172 00173 $subhead Taylor Coefficients$$ 00174 The Taylor coefficient information currently stored in $icode f$$ 00175 (computed by $cref/f.Forward/Forward/$$) is 00176 copied to $icode g$$. 00177 Hence, directly after this operation 00178 $codei% 00179 %g%.size_taylor() == %f%.size_taylor() 00180 %$$ 00181 00182 $subhead Sparsity Patterns$$ 00183 The forward Jacobian sparsity pattern currently stored in $icode f$$ 00184 (computed by $cref/f.ForSparseJac/ForSparseJac/$$) is 00185 copied to $icode g$$. 00186 Hence, directly after this operation 00187 $codei% 00188 %g%.size_forward_bool() == %f%.size_forward_bool() 00189 %g%.size_forward_set() == %f%.size_forward_set() 00190 %$$ 00191 00192 $head OpenMP$$ 00193 $index OpenMP, Dependent$$ 00194 $index Dependent, OpenMP$$ 00195 $index OpenMP, ADFun$$ 00196 $index ADFun, OpenMP$$ 00197 In the case of multi-threading with OpenMP, 00198 the call to $code Independent$$ 00199 and the corresponding call to 00200 $codei% 00201 ADFun<%Base%> %f%( %x%, %y%) 00202 %$$ 00203 or 00204 $codei% 00205 %f%.Dependent( %x%, %y%) 00206 %$$ 00207 must be preformed by the same thread. 00208 00209 00210 $head Example$$ 00211 00212 $subhead Sequence Constructor$$ 00213 The file 00214 $xref/Independent.cpp/$$ 00215 contains an example and test of the sequence constructor. 00216 It returns true if it succeeds and false otherwise. 00217 00218 $subhead Default Constructor$$ 00219 The files 00220 $xref/FunCheck.cpp/$$ 00221 and 00222 $xref/HesLagrangian.cpp/$$ 00223 contain an examples and tests using the default constructor. 00224 They return true if they succeed and false otherwise. 00225 00226 $children% 00227 example/fun_assign.cpp 00228 %$$ 00229 $subhead Assignment Operator$$ 00230 The file 00231 $cref/fun_assign.cpp/$$ 00232 contains an example and test of the $codei%ADFun<%Base%>%$$ 00233 assignment operator. 00234 It returns true if it succeeds and false otherwise. 00235 00236 $end 00237 ---------------------------------------------------------------------------- 00238 */ 00239 00240 CPPAD_BEGIN_NAMESPACE 00241 /*! 00242 \file fun_construct.hpp 00243 ADFun function constructors and assignment operator. 00244 */ 00245 00246 /*! 00247 ADFun default constructor 00248 00249 The C++ syntax for this operation is 00250 \verbatim 00251 ADFun<Base> f 00252 \endverbatim 00253 An empty ADFun object is created. 00254 The Dependent member function, 00255 or the ADFun<Base> assingment operator, 00256 can then be used to put an operation sequence in this ADFun object. 00257 00258 \tparam Base 00259 is the base for the recording that can be stored in this ADFun object; 00260 i.e., operation sequences that were recorded using the type \c AD<Base>. 00261 */ 00262 template <typename Base> 00263 ADFun<Base>::ADFun(void) 00264 : total_num_var_(0), taylor_(CPPAD_NULL) 00265 { } 00266 00267 /*! 00268 ADFun assignment operator 00269 00270 The C++ syntax for this operation is 00271 \verbatim 00272 g = f 00273 \endverbatim 00274 where \c g and \c f are ADFun<Base> ADFun objects. 00275 A copy of the the operation sequence currently stored in \c f 00276 is placed in this ADFun object (called \c g above). 00277 Any information currently stored in this ADFun object is lost. 00278 00279 \tparam Base 00280 is the base for the recording that can be stored in this ADFun object; 00281 i.e., operation sequences that were recorded using the type \c AD<Base>. 00282 00283 \param f 00284 ADFun object containing the operation sequence to be copied. 00285 */ 00286 template <typename Base> 00287 void ADFun<Base>::operator=(const ADFun<Base>& f) 00288 { size_t m = f.Range(); 00289 size_t n = f.Domain(); 00290 00291 // go through member variables in order 00292 // (see ad_fun.hpp for meaning of each variable) 00293 compare_change_ = 0; 00294 taylor_per_var_ = 0; 00295 taylor_col_dim_ = 0; 00296 total_num_var_ = f.total_num_var_; 00297 ind_taddr_.resize(n); 00298 ind_taddr_ = f.ind_taddr_; 00299 dep_taddr_.resize(m); 00300 dep_taddr_ = f.dep_taddr_; 00301 dep_parameter_.resize(m); 00302 dep_parameter_ = f.dep_parameter_; 00303 play_ = f.play_; 00304 if( taylor_ != CPPAD_NULL ) 00305 CPPAD_TRACK_DEL_VEC(taylor_); 00306 taylor_ = CPPAD_NULL; 00307 for_jac_sparse_pack_.resize(0, 0); 00308 for_jac_sparse_set_.resize(0, 0); 00309 00310 // allocate and copy the Taylor coefficients 00311 taylor_per_var_ = f.taylor_per_var_; 00312 taylor_col_dim_ = f.taylor_col_dim_; 00313 size_t length = total_num_var_ * taylor_col_dim_; 00314 if( length > 0 ) taylor_ = CPPAD_TRACK_NEW_VEC(length, taylor_); 00315 size_t i, j; 00316 for(i = 0; i < total_num_var_; i++) 00317 { for(j = 0; j < taylor_per_var_; j++) 00318 { taylor_[ i * taylor_col_dim_ + j ] = 00319 f.taylor_[ i * taylor_col_dim_ + j ]; 00320 } 00321 } 00322 00323 // allocate and copy the forward sparsity information 00324 size_t n_set = f.for_jac_sparse_pack_.n_set(); 00325 size_t end = f.for_jac_sparse_pack_.end(); 00326 if( n_set > 0 ) 00327 { CPPAD_ASSERT_UNKNOWN( n_set == total_num_var_ ); 00328 CPPAD_ASSERT_UNKNOWN( f.for_jac_sparse_set_.n_set() == 0 ); 00329 for_jac_sparse_pack_.resize(n_set, end); 00330 for(i = 0; i < total_num_var_ ; i++) 00331 { for_jac_sparse_pack_.assignment( 00332 i , 00333 i , 00334 f.for_jac_sparse_pack_ 00335 ); 00336 } 00337 } 00338 n_set = f.for_jac_sparse_set_.n_set(); 00339 end = f.for_jac_sparse_set_.end(); 00340 if( n_set > 0 ) 00341 { CPPAD_ASSERT_UNKNOWN( n_set == total_num_var_ ); 00342 CPPAD_ASSERT_UNKNOWN( f.for_jac_sparse_pack_.n_set() == 0 ); 00343 for_jac_sparse_set_.resize(n_set, end); 00344 for(i = 0; i < total_num_var_; i++) 00345 { for_jac_sparse_set_.assignment( 00346 i , 00347 i , 00348 f.for_jac_sparse_set_ 00349 ); 00350 } 00351 } 00352 00353 } 00354 00355 /*! 00356 ADFun constructor from an operation sequence. 00357 00358 The C++ syntax for this operation is 00359 \verbatim 00360 ADFun<Base> f(x, y) 00361 \endverbatim 00362 The operation sequence that started with the previous call 00363 \c Independent(x), and that ends with this operation, is stored 00364 in this \c ADFun<Base> object \c f. 00365 00366 \tparam Base 00367 is the base for the recording that will be stored in the object \c f; 00368 i.e., the operations were recorded using the type \c AD<Base>. 00369 00370 \tparam VectorAD 00371 is a simple vector class with elements of typea \c AD<Base>. 00372 00373 \param x 00374 is the independent variable vector for this ADFun object. 00375 The domain dimension of this object will be the size of \a x. 00376 00377 \param y 00378 is the dependent variable vector for this ADFun object. 00379 The range dimension of this object will be the size of \a y. 00380 00381 \par Taylor Coefficients 00382 A zero order forward mode sweep is done, 00383 and if NDEBUG is not defined the resulting values for the 00384 depenedent variables are checked against the values in \a y. 00385 Thus, the zero order Taylor coefficients 00386 corresponding to the value of the \a x vector 00387 are stored in this ADFun object. 00388 */ 00389 template <typename Base> 00390 template <typename VectorAD> 00391 ADFun<Base>::ADFun(const VectorAD &x, const VectorAD &y) 00392 : total_num_var_(0), taylor_(CPPAD_NULL) 00393 { 00394 CPPAD_ASSERT_KNOWN( 00395 x.size() > 0, 00396 "ADFun<Base>: independent variable vector has size zero." 00397 ); 00398 CPPAD_ASSERT_KNOWN( 00399 Variable(x[0]), 00400 "ADFun<Base>: independent variable vector has been changed." 00401 ); 00402 ADTape<Base> *tape = AD<Base>::tape_ptr(x[0].id_); 00403 CPPAD_ASSERT_KNOWN( 00404 tape->size_independent_ == x.size(), 00405 "ADFun<Base>: independent variable vector has been changed." 00406 ); 00407 size_t j, n = x.size(); 00408 # ifndef NDEBUG 00409 size_t i, m = y.size(); 00410 for(j = 0; j < n; j++) 00411 { CPPAD_ASSERT_KNOWN( 00412 x[j].taddr_ == (j+1), 00413 "ADFun<Base>: independent variable vector has been changed." 00414 ); 00415 CPPAD_ASSERT_KNOWN( 00416 x[j].id_ == x[0].id_, 00417 "ADFun<Base>: independent variable vector has been changed." 00418 ); 00419 } 00420 for(i = 0; i < m; i++) 00421 { CPPAD_ASSERT_KNOWN( 00422 CppAD::Parameter( y[i] ) | (y[i].id_ == x[0].id_) , 00423 "ADFun<Base>: dependent vector contains variables for" 00424 "\na different tape than the independent variables." 00425 ); 00426 } 00427 # endif 00428 00429 // stop the tape and store the operation sequence 00430 Dependent(tape, y); 00431 00432 // allocate memory for one zero order taylor_ coefficient 00433 taylor_per_var_ = 1; 00434 taylor_col_dim_ = 1; 00435 taylor_ = CPPAD_TRACK_NEW_VEC(total_num_var_, taylor_); 00436 00437 // set zero order coefficients corresponding to indpendent variables 00438 CPPAD_ASSERT_UNKNOWN( n == ind_taddr_.size() ); 00439 for(j = 0; j < n; j++) 00440 { CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] == (j+1) ); 00441 CPPAD_ASSERT_UNKNOWN( x[j].taddr_ == (j+1) ); 00442 taylor_[ ind_taddr_[j] ] = x[j].value_; 00443 } 00444 00445 // use independent variable values to fill in values for others 00446 # if CPPAD_USE_FORWARD0SWEEP 00447 compare_change_ = forward0sweep( 00448 false, n, total_num_var_, &play_, taylor_col_dim_, taylor_ 00449 ); 00450 # else 00451 size_t p = 0; 00452 compare_change_ = forward_sweep( 00453 false, p, n, total_num_var_, &play_, taylor_col_dim_, taylor_ 00454 ); 00455 # endif 00456 CPPAD_ASSERT_UNKNOWN( compare_change_ == 0 ); 00457 00458 # ifndef NDEBUG 00459 for(i = 0; i < m; i++) if( taylor_[dep_taddr_[i]] != y[i].value_ ) 00460 { using std::endl; 00461 std::ostringstream buf; 00462 buf << "A dependent variable value is not equal to " 00463 << "its tape evaluation value," << endl 00464 << "perhaps it is nan." << endl 00465 << "Dependent variable value = " 00466 << y[i].value_ << endl 00467 << "Tape evaluation value = " 00468 << taylor_[dep_taddr_[i]] << endl 00469 << "Difference = " 00470 << y[i].value_ - taylor_[dep_taddr_[i]] << endl 00471 ; 00472 CPPAD_ASSERT_KNOWN( 00473 0, 00474 buf.str().c_str() 00475 ); 00476 } 00477 # endif 00478 } 00479 00480 CPPAD_END_NAMESPACE 00481 # endif