00001 # ifndef CPPAD_VEC_AD_INCLUDED 00002 # define CPPAD_VEC_AD_INCLUDED 00003 00004 /* -------------------------------------------------------------------------- 00005 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell 00006 00007 CppAD is distributed under multiple licenses. This distribution is under 00008 the terms of the 00009 Common Public License Version 1.0. 00010 00011 A copy of this license is included in the COPYING file of this distribution. 00012 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00013 -------------------------------------------------------------------------- */ 00014 /* 00015 $begin VecAD$$ 00016 $spell 00017 cppad.hpp 00018 CondExpGt 00019 grep 00020 Ld 00021 vp 00022 Lu 00023 wc 00024 op 00025 Ldp 00026 Ldv 00027 Taylor 00028 VecAD 00029 const 00030 Cpp 00031 $$ 00032 00033 $index VecAD$$ 00034 $index vector, AD index$$ 00035 $index record, AD index$$ 00036 $index tape, AD index$$ 00037 $index index, AD record$$ 00038 00039 $section AD Vectors that Record Index Operations$$ 00040 00041 00042 $head Syntax$$ 00043 $syntax%VecAD<%Base%> %v%(%n%)%$$ 00044 $pre 00045 $$ 00046 $syntax%%v%.size()%$$ 00047 $pre 00048 $$ 00049 $syntax%%b% = %v%[%i%]%$$ 00050 $pre 00051 $$ 00052 $syntax%%r% = %v%[%x%]%$$ 00053 00054 00055 $head Purpose$$ 00056 If either $italic v$$ or $italic x$$ is a 00057 $cref/variable/glossary/Variable/$$, 00058 the indexing operation 00059 $syntax% 00060 %y% = %v%[%x%] 00061 %$$ 00062 is recorded in the corresponding 00063 AD of $italic Base$$ 00064 $xref/glossary/Operation/Sequence/operation sequence/1/$$ and 00065 transferred to the corresponding $xref/ADFun/$$ object $italic f$$. 00066 Such an index can change each time 00067 zero order $xref/Forward//f.Forward/$$ is used; i.e., 00068 $italic f$$ is evaluated with new value for the 00069 $cref/independent variables/glossary/Tape/Independent Variable/$$. 00070 Note that the value of $italic y$$ depends on the value of $italic x$$ 00071 in a discrete fashion and CppAD computes its partial derivative with 00072 respect to $italic x$$ as zero. 00073 00074 $head Alternatives$$ 00075 If only the values in the vector, 00076 and not the indices, 00077 depend on the independent variables, 00078 the class $syntax%%Vector%< AD<%Base%> >%$$ is much more efficient for 00079 storing AD values where $italic Vector$$ is any 00080 $xref/SimpleVector/$$ template class, 00081 If only the indices, 00082 and not the values in the vector, 00083 depend on the independent variables, 00084 The $xref/Discrete/$$ functions are a much more efficient 00085 way to represent these vectors. 00086 00087 $head VecAD<Base>::reference$$ 00088 $index VecAD<Base>::reference$$ 00089 $index reference, VecAD<Base>$$ 00090 The result $italic y$$ has type 00091 $syntax% 00092 VecAD<%Base%>::reference 00093 %$$ 00094 which is very much like the $syntax%AD<%Base%>%$$ type 00095 with some notable exceptions: 00096 00097 $subhead Exceptions$$ 00098 $list number$$ 00099 The object $italic y$$ cannot be used with the 00100 $xref/Value/$$ function to compute the corresponding $italic Base$$ value. 00101 If $italic v$$ is not a $cref/variable/glossary/Variable/$$ 00102 $syntax% 00103 v[%i%] 00104 %$$ 00105 can be used to compute the corresponding $italic Base$$ value. 00106 00107 $lnext 00108 The object $italic y$$ cannot be used 00109 with the $xref/Arithmetic//computed assignments operators/$$ 00110 $code +=$$, 00111 $code -=$$, 00112 $code *=$$, or 00113 $code /=$$. 00114 For example, the following syntax is not valid: 00115 $syntax% 00116 %v%[%x%] += %z%; 00117 %$$ 00118 no matter what the types of $italic z$$. 00119 00120 $lnext 00121 Assignment to $italic y$$ returns a $code void$$. 00122 For example, the following syntax is not valid: 00123 $syntax% 00124 %z% = %v%[%x%] = %u%; 00125 %$$ 00126 no matter what the types of $italic z$$, and $italic u$$. 00127 00128 $lnext 00129 The $xref/CondExp/$$ functions do not accept 00130 $syntax%VecAD<%Base%>::reference%$$ arguments. 00131 For example, the following syntax is not valid: 00132 $syntax% 00133 CondExpGt(%y%, %z%, %u%, %v%) 00134 %$$ 00135 no matter what the types of $italic z$$, $italic u$$, and $italic v$$. 00136 00137 $lnext 00138 The $xref/ParVar//Parameter and Variable/$$ functions cannot be used with 00139 $syntax%VecAD<%Base%>::reference%$$ arguments 00140 (use the entire $syntax%VecAD<%Base%>%$$ vector instead). 00141 00142 $lnext 00143 The vectors passed to $xref/Independent/$$ must have elements 00144 of type $syntax%AD<%Base%>%$$; i.e., $xref/VecAD/$$ vectors 00145 cannot be passed to $code Independent$$. 00146 00147 $lnext 00148 If one uses this type in a 00149 AD of $italic Base$$ 00150 $xref/glossary/Operation/Sequence/operation sequence/1/$$, 00151 $xref/glossary/Sparsity Pattern/sparsity pattern/$$ calculations 00152 ($xref/Sparse/$$) 00153 are only valid for the current independent variable values, 00154 instead of valid for all independent variables. 00155 00156 $lend 00157 00158 $head Constructor$$ 00159 00160 $subhead v$$ 00161 The syntax 00162 $syntax% 00163 VecAD<%Base%> %v%(%n%) 00164 %$$ 00165 creates an $code VecAD$$ object $italic v$$ with 00166 $italic n$$ elements. 00167 The initial value of the elements of $italic v$$ is unspecified. 00168 00169 $head n$$ 00170 The argument $italic n$$ has prototype 00171 $syntax% 00172 size_t %n% 00173 %$$ 00174 00175 $head size$$ 00176 The syntax 00177 $syntax% 00178 %v%.size() 00179 %$$ 00180 returns the number of elements in the vector $italic v$$; 00181 i.e., the value of $italic n$$ when it was constructed. 00182 00183 $head size_t Indexing$$ 00184 We refer to the syntax 00185 $syntax% 00186 %b% = %v%[%i%] 00187 %$$ 00188 as $code size_t$$ indexing of a $code VecAD$$ object. 00189 This indexing is only valid if the vector $italic v$$ is a 00190 $cref/parameter/ParVar/$$; i.e., 00191 it does not depend on the independent variables. 00192 00193 $subhead i$$ 00194 The operand $italic i$$ has prototype 00195 $syntax% 00196 size_t %i% 00197 %$$ 00198 It must be greater than or equal zero 00199 and less than $italic n$$; i.e., less than 00200 the number of elements in $italic v$$. 00201 00202 $subhead b$$ 00203 The result $italic b$$ has prototype 00204 $syntax% 00205 %Base% %b% 00206 %$$ 00207 and is a reference to the $th i$$ element in the vector $italic v$$. 00208 It can be used to change the element value; 00209 for example, 00210 $syntax% 00211 %v%[%i%] = %c% 00212 %$$ 00213 is valid where $italic c$$ is a $italic Base$$ object. 00214 The reference $italic b$$ is no longer valid once the 00215 destructor for $italic v$$ is called; for example, 00216 when $italic v$$ falls out of scope. 00217 00218 $head AD Indexing$$ 00219 We refer to the syntax 00220 $syntax% 00221 %r% = %v%[%x%] 00222 %$$ 00223 as AD indexing of a $code VecAD$$ object. 00224 00225 $subhead x$$ 00226 The argument $italic x$$ has prototype 00227 $syntax% 00228 const AD<%Base%> &%x% 00229 %$$ 00230 The value of $italic x$$ must be greater than or equal zero 00231 and less than $italic n$$; i.e., less than 00232 the number of elements in $italic v$$. 00233 00234 $subhead r$$ 00235 The result $italic y$$ has prototype 00236 $syntax% 00237 VecAD<%Base%>::reference %r% 00238 %$$ 00239 The object $italic r$$ has an AD type and its 00240 operations are recorded as part of the same 00241 AD of $italic Base$$ 00242 $xref/glossary/Operation/Sequence/operation sequence/1/$$ as 00243 for $syntax%AD<%Base%>%$$ objects. 00244 It acts as a reference to the 00245 element with index $latex {\rm floor} (x)$$ in the vector $italic v$$ 00246 ($latex {\rm floor} (x)$$ is 00247 the greatest integer less than or equal $italic x$$). 00248 Because it is a reference, it can be used to change the element value; 00249 for example, 00250 $syntax% 00251 %v%[%x%] = %z% 00252 %$$ 00253 is valid where $italic z$$ is an 00254 $syntax%VecAD<%Base%>::reference%$$ object. 00255 As a reference, $italic r$$ is no longer valid once the 00256 destructor for $italic v$$ is called; for example, 00257 when $italic v$$ falls out of scope. 00258 00259 $head Example$$ 00260 $children% 00261 example/vec_ad.cpp 00262 %$$ 00263 The file 00264 $xref/VecAD.cpp/$$ 00265 contains an example and test using $code VecAD$$ vectors. 00266 It returns true if it succeeds and false otherwise. 00267 00268 00269 $head Speed and Memory$$ 00270 The $xref/VecAD/$$ vector type is inefficient because every 00271 time an element of a vector is accessed, a new CppAD 00272 $xref/glossary/Variable/variable/$$ is created on the tape 00273 using either the $code Ldp$$ or $code Ldv$$ operation 00274 (unless all of the elements of the vector are 00275 $xref/glossary/Parameter/parameters/$$). 00276 The effect of this can be seen by executing the following steps: 00277 00278 $list number$$ 00279 In the file $code cppad/local/ForwardSweep.h$$, 00280 change the definition of $code CppADForwardSweepTrace$$ to 00281 $codep 00282 # define CppADForwardSweepTrace 1 00283 $$ 00284 $lnext 00285 In the $code Example$$ directory, execute the command 00286 $codep 00287 ./OneTest LuVecADOk "lu_vec_ad.cpp -DNDEBUG" > LuVecADOk.log 00288 $$ 00289 This will write a trace of all the forward tape operations, 00290 for the test case $xref/LuVecADOk.cpp/$$, 00291 to the file $code LuVecADOk.log$$. 00292 $lnext 00293 In the $code Example$$ directory execute the commands 00294 $codep 00295 grep "op=" LuVecADOk.log | wc -l 00296 grep "op=Ld[vp]" LuVecADOk.log | wc -l 00297 grep "op=St[vp][vp]" LuVecADOk.log | wc -l 00298 $$ 00299 The first command counts the number of operators in the tracing, 00300 the second counts the number of VecAD load operations, 00301 and the third counts the number of VecAD store operations. 00302 (For CppAD version 05-11-20 these counts were 956, 348, and 118 00303 respectively.) 00304 $lend 00305 00306 $end 00307 ------------------------------------------------------------------------ 00308 */ 00309 00310 # define CPPAD_VEC_AD_COMPUTED_ASSIGNMENT(op, name) \ 00311 VecAD_reference& operator op (const VecAD_reference<Base> &right) \ 00312 { CPPAD_ASSERT_KNOWN( \ 00313 0, \ 00314 "Cannot use a ADVec element on left side of" name \ 00315 ); \ 00316 return *this; \ 00317 } \ 00318 VecAD_reference& operator op (const AD<Base> &right) \ 00319 { CPPAD_ASSERT_KNOWN( \ 00320 0, \ 00321 "Cannot use a ADVec element on left side of" name \ 00322 ); \ 00323 return *this; \ 00324 } \ 00325 VecAD_reference& operator op (const Base &right) \ 00326 { CPPAD_ASSERT_KNOWN( \ 00327 0, \ 00328 "Cannot use a ADVec element on left side of" name \ 00329 ); \ 00330 return *this; \ 00331 } 00332 00333 00334 // BEGIN CppAD namespace 00335 namespace CppAD { 00336 00337 // Element of VecAD 00338 template <class Base> 00339 class VecAD_reference { 00340 friend bool Parameter<Base> (const VecAD<Base> &vec); 00341 friend bool Variable<Base> (const VecAD<Base> &vec); 00342 friend class VecAD<Base>; 00343 friend class ADTape<Base>; 00344 00345 public: 00346 VecAD_reference(VecAD<Base> *v, const AD<Base> &x) 00347 : vec_( v ) , x_(x) 00348 { } 00349 00350 // assignment operators 00351 inline void operator = (const VecAD_reference<Base> &right); 00352 void operator = (const AD<Base> &right); 00353 void operator = (const Base &right); 00354 void operator = (int right); 00355 00356 // computed assignments 00357 CPPAD_VEC_AD_COMPUTED_ASSIGNMENT( += , " += " ) 00358 CPPAD_VEC_AD_COMPUTED_ASSIGNMENT( -= , " -= " ) 00359 CPPAD_VEC_AD_COMPUTED_ASSIGNMENT( *= , " *= " ) 00360 CPPAD_VEC_AD_COMPUTED_ASSIGNMENT( /= , " /= " ) 00361 00362 00363 // AD<Base> constructor 00364 AD<Base> ADBase(void) const 00365 { AD<Base> result; 00366 CPPAD_ASSERT_UNKNOWN( Parameter(result) ); 00367 00368 size_t i = static_cast<size_t>( Integer(x_) ); 00369 CPPAD_ASSERT_UNKNOWN( i < vec_->length_ ); 00370 00371 // value_ corresponding to this element 00372 result.value_ = *(vec_->data_ + i); 00373 00374 // index corresponding to this element 00375 if( Variable(*vec_) ) 00376 { 00377 ADTape<Base> *tape = AD<Base>::tape_ptr(vec_->id_); 00378 CPPAD_ASSERT_UNKNOWN( tape != CPPAD_NULL ); 00379 00380 CPPAD_ASSERT_UNKNOWN( vec_->offset_ > 0 ); 00381 00382 if( IdenticalPar(x_) ) 00383 { // use parameter indexing 00384 tape->RecordLoadOp( 00385 LdpOp, 00386 result, 00387 vec_->offset_, 00388 static_cast<size_t>(i) 00389 ); 00390 CPPAD_ASSERT_UNKNOWN( Parameter(x_) ); 00391 } 00392 else 00393 { // check if we need to convert x to a variable 00394 // note that x is mutable 00395 if( Parameter(x_) ) 00396 { x_.id_ = vec_->id_; 00397 x_.taddr_ = 00398 tape->RecordParOp(x_.value_); 00399 } 00400 00401 // use variable indexing 00402 tape->RecordLoadOp( 00403 LdvOp, 00404 result, 00405 vec_->offset_, 00406 x_.taddr_ 00407 ); 00408 CPPAD_ASSERT_UNKNOWN( 00409 x_.taddr_ > 0 && Variable(x_) 00410 ); 00411 } 00412 } 00413 return result; 00414 } 00415 00416 private: 00417 VecAD<Base> *vec_; // pointer to entire vector 00418 mutable AD<Base> x_; // index for this element 00419 }; 00420 00421 // VecAD 00422 template <class Base> 00423 class VecAD { 00424 friend bool Parameter<Base> (const VecAD<Base> &vec); 00425 friend bool Variable<Base> (const VecAD<Base> &vec); 00426 friend class ADTape<Base>; 00427 friend class VecAD_reference<Base>; 00428 00429 friend std::ostream& operator << <Base> 00430 (std::ostream &os, const VecAD<Base> &vec_); 00431 public: 00432 // declare the user's view of this type here 00433 typedef VecAD_reference<Base> reference; 00434 00435 // default constructor 00436 // initialize id_ to one so not a variable 00437 VecAD(void) : length_(0) , data_(CPPAD_NULL), offset_(0), id_(1) 00438 { CPPAD_ASSERT_UNKNOWN( Parameter(*this) ); } 00439 00440 // constructor 00441 // initialize id_ to one so not a variable 00442 VecAD(size_t n) : length_(n), offset_(0), id_(1) 00443 { data_ = CPPAD_NULL; 00444 if( length_ > 0 ) 00445 { size_t i; 00446 Base zero(0); 00447 data_ = CPPAD_TRACK_NEW_VEC(length_, data_); 00448 00449 // Initialize data to zero so all have same value. 00450 // This uses less memory and avoids a valgrind error 00451 // during TapeRec<Base>::PutPar 00452 for(i = 0; i < length_; i++) 00453 data_[i] = zero; 00454 } 00455 CPPAD_ASSERT_UNKNOWN( Parameter(*this) ); 00456 } 00457 00458 // destructor 00459 ~VecAD(void) 00460 { if( data_ != CPPAD_NULL ) 00461 CPPAD_TRACK_DEL_VEC(data_); 00462 } 00463 00464 // size 00465 size_t size(void) 00466 { return length_; } 00467 00468 // not taped element access 00469 Base &operator[](size_t i) 00470 { 00471 CPPAD_ASSERT_KNOWN( 00472 Parameter(*this), 00473 "VecAD: cannot use size_t indexing because this" 00474 " VecAD vector is a variable." 00475 ); 00476 CPPAD_ASSERT_KNOWN( 00477 i < length_, 00478 "VecAD: element index is >= vector length" 00479 ); 00480 00481 return data_[i]; 00482 } 00483 00484 // taped elemement access 00485 VecAD_reference<Base> operator[](const AD<Base> &x) 00486 { 00487 CPPAD_ASSERT_KNOWN( 00488 0 <= Integer(x), 00489 "VecAD: element index is < zero" 00490 ); 00491 CPPAD_ASSERT_KNOWN( 00492 static_cast<size_t>( Integer(x) ) < length_, 00493 "VecAD: element index is >= vector length" 00494 ); 00495 00496 // if no need to track indexing operation, return now 00497 if( Parameter(*this) & Parameter(x) ) 00498 return VecAD_reference<Base>(this, x); 00499 00500 CPPAD_ASSERT_KNOWN( 00501 Parameter(*this) | Parameter(x) | (id_ == x.id_), 00502 "VecAD: vector and index are variables for" 00503 " different tapes." 00504 ); 00505 00506 if( Parameter(*this) ) 00507 { // must place a copy of vector in tape 00508 offset_ = 00509 AD<Base>::tape_ptr(x.id_)->AddVec(length_, data_); 00510 00511 // advance pointer by one so is always > 0 00512 offset_++; 00513 00514 // tape id corresponding to this offest 00515 id_ = x.id_; 00516 } 00517 00518 return VecAD_reference<Base>(this, x); 00519 } 00520 00521 private: 00522 const size_t length_; // size of this VecAD vector 00523 Base *data_; // value_ of elements of this vector 00524 00525 // offset in cumulate vector corresponding to this object 00526 size_t offset_; 00527 00528 // tape id corresponding to the offset 00529 size_t id_; 00530 }; 00531 00532 00533 template <class Base> 00534 void VecAD_reference<Base>::operator=(const AD<Base> &y) 00535 { 00536 if( Parameter(y) ) 00537 { // fold into the Base type assignment 00538 *this = y.value_; 00539 return; 00540 } 00541 00542 CPPAD_ASSERT_KNOWN( 00543 Parameter(*vec_) | (vec_->id_ == y.id_), 00544 "VecAD: vector and new element value are variables" 00545 "\nfor different tapes." 00546 ); 00547 00548 ADTape<Base> *tape = AD<Base>::tape_ptr(y.id_); 00549 CPPAD_ASSERT_UNKNOWN( tape != CPPAD_NULL ); 00550 if( Parameter(*vec_) ) 00551 { // must place a copy of vector in tape 00552 vec_->offset_ = tape->AddVec(vec_->length_, vec_->data_); 00553 00554 // advance pointer by one so is always > 0 00555 (vec_->offset_)++; 00556 00557 // tape id corresponding to this offest 00558 vec_->id_ = y.id_; 00559 } 00560 CPPAD_ASSERT_UNKNOWN( Variable(*vec_) ); 00561 00562 00563 size_t i = static_cast<size_t>( Integer(x_) ); 00564 CPPAD_ASSERT_UNKNOWN( i < vec_->length_ ); 00565 00566 // assign value both in the element and the original array 00567 *(vec_->data_ + i) = y.value_; 00568 00569 // record the setting of this array element 00570 CPPAD_ASSERT_UNKNOWN( vec_->offset_ > 0 ); 00571 if( Parameter(x_) ) tape->RecordStoreOp( 00572 StpvOp, vec_->offset_, i, y.taddr_ ); 00573 else tape->RecordStoreOp( 00574 StvvOp, vec_->offset_, x_.taddr_, y.taddr_ ); 00575 } 00576 00577 template <class Base> 00578 void VecAD_reference<Base>::operator=(const Base &y) 00579 { 00580 size_t y_taddr; 00581 00582 size_t i = static_cast<size_t>( Integer(x_) ); 00583 CPPAD_ASSERT_UNKNOWN( i < vec_->length_ ); 00584 00585 // assign value both in the element and the original array 00586 *(vec_->data_ + i) = y; 00587 00588 // check if this ADVec object is a parameter 00589 if( Parameter(*vec_) ) 00590 return; 00591 00592 ADTape<Base> *tape = AD<Base>::tape_ptr(vec_->id_); 00593 CPPAD_ASSERT_UNKNOWN( tape != CPPAD_NULL ); 00594 00595 // place a copy of y in the tape 00596 y_taddr = tape->Rec.PutPar(y); 00597 00598 // record the setting of this array element 00599 CPPAD_ASSERT_UNKNOWN( vec_->offset_ > 0 ); 00600 if( Parameter(x_) ) tape->RecordStoreOp( 00601 StppOp, vec_->offset_, i, y_taddr ); 00602 else tape->RecordStoreOp( 00603 StvpOp, vec_->offset_, x_.taddr_, y_taddr ); 00604 } 00605 00606 // fold this case into AD<Base> case above 00607 template <class Base> 00608 inline void VecAD_reference<Base>::operator=(const VecAD_reference<Base> &y) 00609 { *this = y.ADBase(); } 00610 00611 // fold this case into Base case above 00612 template <class Base> 00613 inline void VecAD_reference<Base>::operator=(int y) 00614 { *this = Base(y); } 00615 00616 template <class Base> 00617 inline std::ostream& operator << (std::ostream &os, const VecAD<Base> &v) 00618 { 00619 os << "vecAD( length = " << v.length_ 00620 << ", offset = " << v.offset_ << ")"; 00621 return os; 00622 } 00623 00624 00625 } // END CppAD namespace 00626 00627 # undef CPPAD_VEC_AD_COMPUTED_ASSIGNMENT 00628 00629 # endif