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