# source:trunk/example/atomic/reciprocal.cpp@2903

Last change on this file since 2903 was 2903, checked in by bradbell, 7 years ago

Add lines from atomic_base documentation to corresponding examples.

package.sh: ignore the bug/build directory.

• Property svn:keywords set to Id
File size: 13.7 KB
Line
1// $Id: reciprocal.cpp 2903 2013-09-19 15:53:13Z bradbell$
2/* --------------------------------------------------------------------------
4
6the terms of the
7                    Eclipse Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
11-------------------------------------------------------------------------- */
12
13/*
14$begin atomic_reciprocal.cpp$$15 16section Reciprocal as an Atomic Operation: Example and Test$$ 17$index reciprocal, atomic operation$$18index simple, atomic operation$$
19$index atomic, simple operation$$20index operation, simple atomic$$ 21 22$head Theory$$23This example demonstrates using cref atomic_base$$
24to define the operation
25$latex f : \B{R}^n \rightarrow \B{R}^m$$where 26latex n = 1$$,$latex m = 1$$, and latex f(x) = 1 / x$$.
27
28$nospell 29 30$head Start Class Definition$$31codep */ 32# include <cppad/cppad.hpp> 33namespace { // isolate items below to this file 34using CppAD::vector; // abbreviate as vector 35// 36// a utility to compute the union of two sets. 37void my_union( 38 std::set<size_t>& result , 39 const std::set<size_t>& left , 40 const std::set<size_t>& right ) 41{ std::set<size_t> temp; 42 std::set_union( 43 left.begin() , 44 left.end() , 45 right.begin() , 46 right.end() , 47 std::inserter(temp, temp.begin()) 48 ); 49 result.swap(temp); 50} 51// 52class atomic_reciprocal : public CppAD::atomic_base<double> { 53/*$$
54$head Constructor $$55codep */ 56 public: 57 // constructor (could use const char* for name) 58 atomic_reciprocal(const std::string& name) : 59 CppAD::atomic_base<double>(name) 60 { } 61 private: 62/*$$ 63$head forward$$64codep */ 65 // forward mode routine called by CppAD 66 virtual bool forward( 67 size_t q , 68 size_t p , 69 const vector<bool>& vx , 70 vector<bool>& vy , 71 const vector<double>& tx , 72 vector<double>& ty 73 ) 74 { size_t n = tx.size() / (p + 1); 75 size_t m = ty.size() / (p + 1); 76 assert( n == 1 ); 77 assert( m == 1 ); 78 assert( q <= p ); 79 80 // return flag 81 bool ok = p <= 2; 82 83 // check for defining variable information 84 // This case must always be implemented 85 if( vx.size() > 0 ) 86 vy[0] = vx[0]; 87 88 // Order zero forward mode. 89 // This case must always be implemented 90 // y^0 = f( x^0 ) = 1 / x^0 91 double f = 1. / tx[0]; 92 if( q <= 0 ) 93 ty[0] = f; 94 if( p <= 0 ) 95 return ok; 96 assert( vx.size() == 0 ); 97 98 // Order one forward mode. 99 // This case needed if first order forward mode is used. 100 // y^1 = f'( x^0 ) x^1 101 double fp = - f / tx[0]; 102 if( q <= 1 ) 103 ty[1] = fp * tx[1]; 104 if( p <= 1 ) 105 return ok; 106 107 // Order two forward mode. 108 // This case needed if second order forward mode is used. 109 // Y''(t) = X'(t)^\R{T} f''[X(t)] X'(t) + f'[X(t)] X''(t) 110 // 2 y^2 = x^1 * f''( x^0 ) x^1 + 2 f'( x^0 ) x^2 111 double fpp = - 2.0 * fp / tx[0]; 112 ty[2] = tx[1] * fpp * tx[1] / 2.0 + fp * tx[2]; 113 if( p <= 2 ) 114 return ok; 115 116 // Assume we are not using forward mode with order > 2 117 assert( ! ok ); 118 return ok; 119 } 120/*$$
121$head reverse$$122codep */ 123 // reverse mode routine called by CppAD 124 virtual bool reverse( 125 size_t p , 126 const vector<double>& tx , 127 const vector<double>& ty , 128 vector<double>& px , 129 const vector<double>& py 130 ) 131 { size_t n = tx.size() / (p + 1); 132 size_t m = ty.size() / (p + 1); 133 assert( px.size() == n * (p + 1) ); 134 assert( py.size() == m * (p + 1) ); 135 assert( n == 1 ); 136 assert( m == 1 ); 137 bool ok = p <= 2; 138 139 double f, fp, fpp, fppp; 140 switch(p) 141 { case 0: 142 // This case needed if first order reverse mode is used 143 // reverse: F^0 ( tx ) = y^0 = f( x^0 ) 144 f = ty[0]; 145 fp = - f / tx[0]; 146 px[0] = py[0] * fp;; 147 assert(ok); 148 break; 149 150 case 1: 151 // This case needed if second order reverse mode is used 152 // reverse: F^1 ( tx ) = y^1 = f'( x^0 ) x^1 153 f = ty[0]; 154 fp = - f / tx[0]; 155 fpp = - 2.0 * fp / tx[0]; 156 px[1] = py[1] * fp; 157 px[0] = py[1] * fpp * tx[1]; 158 // reverse: F^0 ( tx ) = y^0 = f( x^0 ); 159 px[0] += py[0] * fp; 160 assert(ok); 161 break; 162 163 case 2: 164 // This needed if third order reverse mode is used 165 // reverse: F^2 ( tx ) = y^2 = 166 // = x^1 * f''( x^0 ) x^1 / 2 + f'( x^0 ) x^2 167 f = ty[0]; 168 fp = - f / tx[0]; 169 fpp = - 2.0 * fp / tx[0]; 170 fppp = - 3.0 * fpp / tx[0]; 171 px[2] = py[2] * fp; 172 px[1] = py[2] * fpp * tx[1]; 173 px[0] = py[2] * tx[1] * fppp * tx[1] / 2.0 + fpp * tx[2]; 174 // reverse: F^1 ( tx ) = y^1 = f'( x^0 ) x^1 175 px[1] += py[1] * fp; 176 px[0] += py[1] * fpp * tx[1]; 177 // reverse: F^0 ( tx ) = y^0 = f( x^0 ); 178 px[0] += py[0] * fp; 179 assert(ok); 180 break; 181 182 default: 183 assert(!ok); 184 } 185 return ok; 186 } 187/*$$ 188$head for_sparse_jac$$189codep */ 190 // forward Jacobian bool sparsity routine called by CppAD 191 virtual bool for_sparse_jac( 192 size_t q , 193 const vector<bool>& r , 194 vector<bool>& s ) 195 { // This function needed if using f.ForSparseJac 196 // with afun.option( CppAD::atomic_base<double>::bool_sparsity_enum ) 197 size_t n = r.size() / q; 198 size_t m = s.size() / q; 199 assert( n == 1 ); 200 assert( m == 1 ); 201 202 // sparsity for S(x) = f'(x) * R is same as sparsity for R 203 for(size_t j = 0; j < q; j++) 204 s[j] = r[j]; 205 206 return true; 207 } 208 // forward Jacobian set sparsity routine called by CppAD 209 virtual bool for_sparse_jac( 210 size_t q , 211 const vector< std::set<size_t> >& r , 212 vector< std::set<size_t> >& s ) 213 { // This function needed if using f.ForSparseJac 214 // with afun.option( CppAD::atomic_base<double>::set_sparsity_enum ) 215 size_t n = r.size(); 216 size_t m = s.size(); 217 assert( n == 1 ); 218 assert( m == 1 ); 219 220 // sparsity for S(x) = f'(x) * R is same as sparsity for R 221 s[0] = r[0]; 222 223 return true; 224 } 225/*$$
226$head rev_sparse_jac$$227codep */ 228 // reverse Jacobian bool sparsity routine called by CppAD 229 virtual bool rev_sparse_jac( 230 size_t q , 231 const vector<bool>& rt , 232 vector<bool>& st ) 233 { // This function needed if using RevSparseJac 234 // with afun.option( CppAD::atomic_base<double>::bool_sparsity_enum ) 235 size_t n = st.size() / q; 236 size_t m = rt.size() / q; 237 assert( n == 1 ); 238 assert( m == 1 ); 239 240 // sparsity for S(x)^T = f'(x)^T * R^T is same as sparsity for R^T 241 for(size_t i = 0; i < q; i++) 242 st[i] = rt[i]; 243 244 return true; 245 } 246 // reverse Jacobian set sparsity routine called by CppAD 247 virtual bool rev_sparse_jac( 248 size_t q , 249 const vector< std::set<size_t> >& rt , 250 vector< std::set<size_t> >& st ) 251 { // This function needed if using RevSparseJac 252 // with afun.option( CppAD::atomic_base<double>::set_sparsity_enum ) 253 size_t n = st.size(); 254 size_t m = rt.size(); 255 assert( n == 1 ); 256 assert( m == 1 ); 257 258 // sparsity for S(x)^T = f'(x)^T * R^T is same as sparsity for R^T 259 st[0] = rt[0]; 260 261 return true; 262 } 263/*$$ 264$head rev_sparse_hes$$265codep */ 266 // reverse Hessian bool sparsity routine called by CppAD 267 virtual bool rev_sparse_hes( 268 const vector<bool>& vx, 269 const vector<bool>& s , 270 vector<bool>& t , 271 size_t q , 272 const vector<bool>& r , 273 const vector<bool>& u , 274 vector<bool>& v ) 275 { // This function needed if using RevSparseHes 276 // with afun.option( CppAD::atomic_base<double>::bool_sparsity_enum ) 277 size_t m = s.size(); 278 size_t n = t.size(); 279 assert( r.size() == n * q ); 280 assert( u.size() == m * q ); 281 assert( v.size() == n * q ); 282 assert( n == 1 ); 283 assert( m == 1 ); 284 285 // There are no cross term second derivatives for this case, 286 // so it is not necessary to vx. 287 288 // sparsity for T(x) = S(x) * f'(x) is same as sparsity for S 289 t[0] = s[0]; 290 291 // V(x) = f'(x)^T * g''(y) * f'(x) * R + g'(y) * f''(x) * R 292 // U(x) = g''(y) * f'(x) * R 293 // S(x) = g'(y) 294 295 // back propagate the sparsity for U, note f'(x) may be non-zero; 296 size_t j; 297 for(j = 0; j < q; j++) 298 v[j] = u[j]; 299 300 // include forward Jacobian sparsity in Hessian sparsity 301 // (note sparsty for f''(x) * R same as for R) 302 if( s[0] ) 303 { for(j = 0; j < q; j++) 304 v[j] |= r[j]; 305 } 306 307 return true; 308 } 309 // reverse Hessian set sparsity routine called by CppAD 310 virtual bool rev_sparse_hes( 311 const vector<bool>& vx, 312 const vector<bool>& s , 313 vector<bool>& t , 314 size_t q , 315 const vector< std::set<size_t> >& r , 316 const vector< std::set<size_t> >& u , 317 vector< std::set<size_t> >& v ) 318 { // This function needed if using RevSparseHes 319 // with afun.option( CppAD::atomic_base<double>::set_sparsity_enum ) 320 size_t n = vx.size(); 321 size_t m = s.size(); 322 assert( t.size() == n ); 323 assert( r.size() == n ); 324 assert( u.size() == m ); 325 assert( v.size() == n ); 326 assert( n == 1 ); 327 assert( m == 1 ); 328 329 // There are no cross term second derivatives for this case, 330 // so it is not necessary to vx. 331 332 // sparsity for T(x) = S(x) * f'(x) is same as sparsity for S 333 t[0] = s[0]; 334 335 // V(x) = f'(x)^T * g''(y) * f'(x) * R + g'(y) * f''(x) * R 336 // U(x) = g''(y) * f'(x) * R 337 // S(x) = g'(y) 338 339 // back propagate the sparsity for U, note f'(x) may be non-zero; 340 v[0] = u[0]; 341 342 // include forward Jacobian sparsity in Hessian sparsity 343 // (note sparsty for f''(x) * R same as for R) 344 if( s[0] ) 345 my_union(v[0], v[0], r[0] ); 346 347 return true; 348 } 349/*$$
350$head End Class Definition$$351codep */ 352}; // End of atomic_reciprocal class 353} // End empty namespace 354 355/*$$ 356$head Use Atomic Function$$357codep */ 358bool reciprocal(void) 359{ bool ok = true; 360 using CppAD::AD; 361 using CppAD::NearEqual; 362 double eps = 10. * CppAD::numeric_limits<double>::epsilon(); 363 364 // -------------------------------------------------------------------- 365 // Create the atomic reciprocal object 366 atomic_reciprocal afun("atomic_reciprocal"); 367 // -------------------------------------------------------------------- 368 // Create the function f(x) 369 // 370 // domain space vector 371 size_t n = 1; 372 double x0 = 0.5; 373 vector< AD<double> > ax(n); 374 ax[0] = x0; 375 376 // declare independent variables and start tape recording 377 CppAD::Independent(ax); 378 379 // range space vector 380 size_t m = 1; 381 vector< AD<double> > ay(m); 382 383 // call user function and store reciprocal(x) in au[0] 384 vector< AD<double> > au(m); 385 afun(ax, au); // u = 1 / x 386 387 // now use AD division to invert to invert the operation 388 ay[0] = 1.0 / au[0]; // y = 1 / u = x 389 390 // create f: x -> y and stop tape recording 391 CppAD::ADFun<double> f; 392 f.Dependent (ax, ay); // f(x) = x 393 394 // -------------------------------------------------------------------- 395 // Check forward mode results 396 // 397 // check function value 398 double check = x0; 399 ok &= NearEqual( Value(ay[0]) , check, eps, eps); 400 401 // check zero order forward mode 402 size_t p; 403 vector<double> x_p(n), y_p(m); 404 p = 0; 405 x_p[0] = x0; 406 y_p = f.Forward(p, x_p); 407 ok &= NearEqual(y_p[0] , check, eps, eps); 408 409 // check first order forward mode 410 p = 1; 411 x_p[0] = 1; 412 y_p = f.Forward(p, x_p); 413 check = 1.; 414 ok &= NearEqual(y_p[0] , check, eps, eps); 415 416 // check second order forward mode 417 p = 2; 418 x_p[0] = 0; 419 y_p = f.Forward(p, x_p); 420 check = 0.; 421 ok &= NearEqual(y_p[0] , check, eps, eps); 422 423 // -------------------------------------------------------------------- 424 // Check reverse mode results 425 // 426 // third order reverse mode 427 p = 3; 428 vector<double> w(m), dw(n * p); 429 w[0] = 1.; 430 dw = f.Reverse(p, w); 431 check = 1.; 432 ok &= NearEqual(dw[0] , check, eps, eps); 433 check = 0.; 434 ok &= NearEqual(dw[1] , check, eps, eps); 435 ok &= NearEqual(dw[2] , check, eps, eps); 436 437 // -------------------------------------------------------------------- 438 // forward mode sparstiy pattern 439 size_t q = n; 440 CppAD::vectorBool r1(n * q), s1(m * q); 441 r1[0] = true; // compute sparsity pattern for x[0] 442 // 443 afun.option( CppAD::atomic_base<double>::bool_sparsity_enum ); 444 s1 = f.ForSparseJac(q, r1); 445 ok &= s1[0] == true; // f[0] depends on x[0] 446 // 447 afun.option( CppAD::atomic_base<double>::set_sparsity_enum ); 448 s1 = f.ForSparseJac(q, r1); 449 ok &= s1[0] == true; // f[0] depends on x[0] 450 451 // -------------------------------------------------------------------- 452 // reverse mode sparstiy pattern 453 p = m; 454 CppAD::vectorBool s2(p * m), r2(p * n); 455 s2[0] = true; // compute sparsity pattern for f[0] 456 // 457 afun.option( CppAD::atomic_base<double>::bool_sparsity_enum ); 458 r2 = f.RevSparseJac(p, s2); 459 ok &= r2[0] == true; // f[0] depends on x[0] 460 // 461 afun.option( CppAD::atomic_base<double>::set_sparsity_enum ); 462 r2 = f.RevSparseJac(p, s2); 463 ok &= r2[0] == true; // f[0] depends on x[0] 464 465 // -------------------------------------------------------------------- 466 // Hessian sparsity (using previous ForSparseJac call) 467 CppAD::vectorBool s3(m), h(q * n); 468 s3[0] = true; // compute sparsity pattern for f[0] 469 // 470 afun.option( CppAD::atomic_base<double>::bool_sparsity_enum ); 471 h = f.RevSparseHes(q, s3); 472 ok &= h[0] == true; // second partial of f[0] w.r.t. x[0] may be non-zero 473 // 474 afun.option( CppAD::atomic_base<double>::set_sparsity_enum ); 475 h = f.RevSparseHes(q, s3); 476 ok &= h[0] == true; // second partial of f[0] w.r.t. x[0] may be non-zero 477 478 // ----------------------------------------------------------------- 479 // Free all temporary work space associated with atomic_base objects. 480 // (If there are future calls to user atomic functions, they will 481 // create new temporary work space.) 482 CppAD::atomic_base<double>::clear(); 483 484 return ok; 485} 486/*$$
487$$comment end nospell$$
488\$end
489*/
Note: See TracBrowser for help on using the repository browser.