Prev Next Index-> contents reference index search external Up-> CppAD Appendix deprecated old_atomic old_reciprocal.cpp Appendix-> Faq directory Theory glossary Bib wish_list whats_new deprecated compare_c numeric_ad addon License deprecated-> include_deprecated FunDeprecated CompareChange omp_max_thread TrackNewDel omp_alloc memory_leak epsilon test_vector cppad_ipopt_nlp old_atomic zdouble autotools old_atomic-> old_reciprocal.cpp old_usead_1.cpp old_usead_2.cpp old_tan.cpp old_mat_mul.cpp old_reciprocal.cpp Headings-> Deprecated 2013-05-27 Theory

$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
Old Atomic Operation Reciprocal: Example and Test

Deprecated 2013-05-27
This example has been deprecated; see atomic_reciprocal.cpp instead.

Theory
The example below defines the user atomic function $f : \B{R}^n \rightarrow \B{R}^m$ where $n = 1$, $m = 1$, and $f(x) = 1 / x$.  # include <cppad/cppad.hpp> namespace { // Begin empty namespace using CppAD::vector; // ---------------------------------------------------------------------- // a utility to compute the union of two sets. using CppAD::set_union; // ---------------------------------------------------------------------- // forward mode routine called by CppAD bool reciprocal_forward( size_t id , size_t k , size_t n , size_t m , const vector<bool>& vx , vector<bool>& vy , const vector<double>& tx , vector<double>& ty ) { assert( id == 0 ); assert( n == 1 ); assert( m == 1 ); assert( k == 0 || vx.size() == 0 ); bool ok = false; double f, fp, fpp; // Must always define the case k = 0. // Do not need case k if not using f.Forward(q, xp) for q >= k. switch(k) { case 0: // this case must be implemented if( vx.size() > 0 ) vy[0] = vx[0]; // y^0 = f( x^0 ) = 1 / x^0 ty[0] = 1. / tx[0]; ok = true; break; case 1: // needed if first order forward mode is used assert( vx.size() == 0 ); // y^1 = f'( x^0 ) x^1 f = ty[0]; fp = - f / tx[0]; ty[1] = fp * tx[1]; ok = true; break; case 2: // needed if second order forward mode is used assert( vx.size() == 0 ); // Y''(t) = X'(t)^\R{T} f''[X(t)] X'(t) + f'[X(t)] X''(t) // 2 y^2 = x^1 * f''( x^0 ) x^1 + 2 f'( x^0 ) x^2 f = ty[0]; fp = - f / tx[0]; fpp = - 2.0 * fp / tx[0]; ty[2] = tx[1] * fpp * tx[1] / 2.0 + fp * tx[2]; ok = true; break; } return ok; } // ---------------------------------------------------------------------- // reverse mode routine called by CppAD bool reciprocal_reverse( size_t id , size_t k , size_t n , size_t m , const vector<double>& tx , const vector<double>& ty , vector<double>& px , const vector<double>& py ) { // Do not need case k if not using f.Reverse(k+1, w). assert( id == 0 ); assert( n == 1 ); assert( m == 1 ); bool ok = false; double f, fp, fpp, fppp; switch(k) { case 0: // needed if first order reverse mode is used // reverse: F^0 ( tx ) = y^0 = f( x^0 ) f = ty[0]; fp = - f / tx[0]; px[0] = py[0] * fp;; ok = true; break; case 1: // needed if second order reverse mode is used // reverse: F^1 ( tx ) = y^1 = f'( x^0 ) x^1 f = ty[0]; fp = - f / tx[0]; fpp = - 2.0 * fp / tx[0]; px[1] = py[1] * fp; px[0] = py[1] * fpp * tx[1]; // reverse: F^0 ( tx ) = y^0 = f( x^0 ); px[0] += py[0] * fp; ok = true; break; case 2: // needed if third order reverse mode is used // reverse: F^2 ( tx ) = y^2 = // = x^1 * f''( x^0 ) x^1 / 2 + f'( x^0 ) x^2 f = ty[0]; fp = - f / tx[0]; fpp = - 2.0 * fp / tx[0]; fppp = - 3.0 * fpp / tx[0]; px[2] = py[2] * fp; px[1] = py[2] * fpp * tx[1]; px[0] = py[2] * tx[1] * fppp * tx[1] / 2.0 + fpp * tx[2]; // reverse: F^1 ( tx ) = y^1 = f'( x^0 ) x^1 px[1] += py[1] * fp; px[0] += py[1] * fpp * tx[1]; // reverse: F^0 ( tx ) = y^0 = f( x^0 ); px[0] += py[0] * fp; ok = true; break; } return ok; } // ---------------------------------------------------------------------- // forward Jacobian sparsity routine called by CppAD bool reciprocal_for_jac_sparse( size_t id , size_t n , size_t m , size_t p , const vector< std::set<size_t> >& r , vector< std::set<size_t> >& s ) { // Can just return false if not using f.ForSparseJac assert( id == 0 ); assert( n == 1 ); assert( m == 1 ); // sparsity for S(x) = f'(x) * R is same as sparsity for R s[0] = r[0]; return true; } // ---------------------------------------------------------------------- // reverse Jacobian sparsity routine called by CppAD bool reciprocal_rev_jac_sparse( size_t id , size_t n , size_t m , size_t p , vector< std::set<size_t> >& r , const vector< std::set<size_t> >& s ) { // Can just return false if not using RevSparseJac. assert( id == 0 ); assert( n == 1 ); assert( m == 1 ); // sparsity for R(x) = S * f'(x) is same as sparsity for S for(size_t q = 0; q < p; q++) r[q] = s[q]; return true; } // ---------------------------------------------------------------------- // reverse Hessian sparsity routine called by CppAD bool reciprocal_rev_hes_sparse( size_t id , size_t n , size_t m , size_t p , const vector< std::set<size_t> >& r , const vector<bool>& s , vector<bool>& t , const vector< std::set<size_t> >& u , vector< std::set<size_t> >& v ) { // Can just return false if not use RevSparseHes. assert( id == 0 ); assert( n == 1 ); assert( m == 1 ); // sparsity for T(x) = S(x) * f'(x) is same as sparsity for S t[0] = s[0]; // V(x) = [ f'(x)^T * g''(y) * f'(x) + g'(y) * f''(x) ] * R // U(x) = g''(y) * f'(x) * R // S(x) = g'(y) // back propagate the sparsity for U because derivative of // reciprocal may be non-zero v[0] = u[0]; // convert forward Jacobian sparsity to Hessian sparsity // because second derivative of reciprocal may be non-zero if( s[0] ) v[0] = set_union(v[0], r[0] ); return true; } // --------------------------------------------------------------------- // Declare the AD<double> routine reciprocal(id, ax, ay) CPPAD_USER_ATOMIC( reciprocal , CppAD::vector , double , reciprocal_forward , reciprocal_reverse , reciprocal_for_jac_sparse , reciprocal_rev_jac_sparse , reciprocal_rev_hes_sparse ) } // End empty namespace bool old_reciprocal(void) { bool ok = true; using CppAD::AD; using CppAD::NearEqual; double eps = 10. * CppAD::numeric_limits<double>::epsilon(); // -------------------------------------------------------------------- // Create the function f(x) // // domain space vector size_t n = 1; double x0 = 0.5; vector< AD<double> > ax(n); ax[0] = x0; // declare independent variables and start tape recording CppAD::Independent(ax); // range space vector size_t m = 1; vector< AD<double> > ay(m); // call user function and store reciprocal(x) in au[0] vector< AD<double> > au(m); size_t id = 0; // not used reciprocal(id, ax, au); // u = 1 / x // call user function and store reciprocal(u) in ay[0] reciprocal(id, au, ay); // y = 1 / u = x // create f: x -> y and stop tape recording CppAD::ADFun<double> f; f.Dependent (ax, ay); // f(x) = x // -------------------------------------------------------------------- // Check forward mode results // // check function value double check = x0; ok &= NearEqual( Value(ay[0]) , check, eps, eps); // check zero order forward mode size_t q; vector<double> x_q(n), y_q(m); q = 0; x_q[0] = x0; y_q = f.Forward(q, x_q); ok &= NearEqual(y_q[0] , check, eps, eps); // check first order forward mode q = 1; x_q[0] = 1; y_q = f.Forward(q, x_q); check = 1.; ok &= NearEqual(y_q[0] , check, eps, eps); // check second order forward mode q = 2; x_q[0] = 0; y_q = f.Forward(q, x_q); check = 0.; ok &= NearEqual(y_q[0] , check, eps, eps); // -------------------------------------------------------------------- // Check reverse mode results // // third order reverse mode q = 3; vector<double> w(m), dw(n * q); w[0] = 1.; dw = f.Reverse(q, w); check = 1.; ok &= NearEqual(dw[0] , check, eps, eps); check = 0.; ok &= NearEqual(dw[1] , check, eps, eps); ok &= NearEqual(dw[2] , check, eps, eps); // -------------------------------------------------------------------- // forward mode sparstiy pattern size_t p = n; CppAD::vectorBool r1(n * p), s1(m * p); r1[0] = true; // compute sparsity pattern for x[0] s1 = f.ForSparseJac(p, r1); ok &= s1[0] == true; // f[0] depends on x[0] // -------------------------------------------------------------------- // reverse mode sparstiy pattern q = m; CppAD::vectorBool s2(q * m), r2(q * n); s2[0] = true; // compute sparsity pattern for f[0] r2 = f.RevSparseJac(q, s2); ok &= r2[0] == true; // f[0] depends on x[0] // -------------------------------------------------------------------- // Hessian sparsity (using previous ForSparseJac call) CppAD::vectorBool s3(m), h(p * n); s3[0] = true; // compute sparsity pattern for f[0] h = f.RevSparseHes(p, s3); ok &= h[0] == true; // second partial of f[0] w.r.t. x[0] may be non-zero // ----------------------------------------------------------------- // Free all temporary work space associated with old_atomic objects. // (If there are future calls to user atomic functions, they will // create new temporary work space.) CppAD::user_atomic<double>::clear(); return ok; } 
Input File: example/deprecated/old_reciprocal.cpp