Prev Next Index-> contents reference index search external Up-> CppAD AD ADValued atomic atomic_base atomic_reciprocal.cpp ADValued-> Arithmetic unary_standard_math binary_math CondExp Discrete numeric_limits atomic atomic-> checkpoint atomic_base atomic_base-> atomic_ctor atomic_option atomic_afun atomic_forward atomic_reverse atomic_for_sparse_jac atomic_rev_sparse_jac atomic_for_sparse_hes atomic_rev_sparse_hes atomic_base_clear atomic_get_started.cpp atomic_norm_sq.cpp atomic_reciprocal.cpp atomic_set_sparsity.cpp atomic_tangent.cpp atomic_eigen_mat_mul.cpp atomic_eigen_mat_inv.cpp atomic_eigen_cholesky.cpp atomic_mat_mul.cpp atomic_reciprocal.cpp Headings-> Theory sparsity Start Class Definition Constructor forward reverse for_sparse_jac rev_sparse_jac rev_sparse_hes End Class Definition Use Atomic Function ---..Constructor ---..Recording ---..forward ---..reverse ---..for_sparse_jac ---..rev_sparse_jac ---..rev_sparse_hes

$\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}} }$
Reciprocal as an Atomic Operation: Example and Test

Theory
This example demonstrates using atomic_base to define the operation $f : \B{R}^n \rightarrow \B{R}^m$ where $n = 1$, $m = 1$, and $f(x) = 1 / x$.

sparsity
This example only uses set sparsity patterns.

Start Class Definition
# include <cppad/cppad.hpp>
namespace {           // isolate items below to this file
using CppAD::vector;  // abbreviate as vector
//
// a utility to compute the union of two sets.
//
class atomic_reciprocal : public CppAD::atomic_base<double> {

Constructor
public:
// constructor (could use const char* for name)
atomic_reciprocal(const std::string& name) :
// this exmaple only uses set sparsity patterns
{ }
private:

forward
     // forward mode routine called by CppAD
virtual bool forward(
size_t                    p ,
size_t                    q ,
const vector<bool>&      vx ,
vector<bool>&      vy ,
const vector<double>&    tx ,
vector<double>&    ty
)
{
# ifndef NDEBUG
size_t n = tx.size() / (q + 1);
size_t m = ty.size() / (q + 1);
# endif
assert( n == 1 );
assert( m == 1 );
assert( p <= q );

// return flag
bool ok = q <= 2;

// check for defining variable information
// This case must always be implemented
if( vx.size() > 0 )
vy[0] = vx[0];

// Order zero forward mode.
// This case must always be implemented
// y^0 = f( x^0 ) = 1 / x^0
double f = 1. / tx[0];
if( p <= 0 )
ty[0] = f;
if( q <= 0 )
return ok;
assert( vx.size() == 0 );

// Order one forward mode.
// This case needed if first order forward mode is used.
// y^1 = f'( x^0 ) x^1
double fp = - f / tx[0];
if( p <= 1 )
ty[1] = fp * tx[1];
if( q <= 1 )
return ok;

// Order two forward mode.
// This case needed if second order forward mode is used.
// 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
double fpp  = - 2.0 * fp / tx[0];
ty[2] = tx[1] * fpp * tx[1] / 2.0 + fp * tx[2];
if( q <= 2 )
return ok;

// Assume we are not using forward mode with order > 2
assert( ! ok );
return ok;
}

reverse
     // reverse mode routine called by CppAD
virtual bool reverse(
size_t                    q ,
const vector<double>&    tx ,
const vector<double>&    ty ,
vector<double>&    px ,
const vector<double>&    py
)
{
# ifndef NDEBUG
size_t n = tx.size() / (q + 1);
size_t m = ty.size() / (q + 1);
# endif
assert( px.size() == n * (q + 1) );
assert( py.size() == m * (q + 1) );
assert( n == 1 );
assert( m == 1 );
bool ok = q <= 2;

double f, fp, fpp, fppp;
switch(q)
{     case 0:
// This case 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;;
assert(ok);
break;

case 1:
// This case 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;
assert(ok);
break;

case 2:
// This 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;
assert(ok);
break;

default:
assert(!ok);
}
return ok;
}

for_sparse_jac
     // forward Jacobian set sparsity routine called by CppAD
virtual bool for_sparse_jac(
size_t                                p ,
const vector< std::set<size_t> >&     r ,
vector< std::set<size_t> >&     s ,
const vector<double>&                 x )
{     // This function needed if using f.ForSparseJac
# ifndef NDEBUG
size_t n = r.size();
size_t m = s.size();
# endif
assert( n == x.size() );
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;
}

rev_sparse_jac
     // reverse Jacobian set sparsity routine called by CppAD
virtual bool rev_sparse_jac(
size_t                                p  ,
const vector< std::set<size_t> >&     rt ,
vector< std::set<size_t> >&     st ,
const vector<double>&                 x  )
{     // This function needed if using RevSparseJac or optimize
# ifndef NDEBUG
size_t n = st.size();
size_t m = rt.size();
# endif
assert( n == x.size() );
assert( n == 1 );
assert( m == 1 );

// sparsity for S(x)^T = f'(x)^T * R^T is same as sparsity for R^T
st[0] = rt[0];

return true;
}

rev_sparse_hes
     // reverse Hessian set sparsity routine called by CppAD
virtual bool rev_sparse_hes(
const vector<bool>&                   vx,
const vector<bool>&                   s ,
vector<bool>&                   t ,
size_t                                p ,
const vector< std::set<size_t> >&     r ,
const vector< std::set<size_t> >&     u ,
vector< std::set<size_t> >&     v ,
const vector<double>&                 x )
{     // This function needed if using RevSparseHes
# ifndef NDEBUG
size_t n = vx.size();
size_t m = s.size();
# endif
assert( x.size() == n );
assert( t.size() == n );
assert( r.size() == n );
assert( u.size() == m );
assert( v.size() == n );
assert( n == 1 );
assert( m == 1 );

// There are no cross term second derivatives for this case,
// so it is not necessary to vx.

// 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) * R  +  g'(y) * f''(x) * R
// U(x) = g''(y) * f'(x) * R
// S(x) = g'(y)

// back propagate the sparsity for U, note f'(x) may be non-zero;
v[0] = u[0];

// include forward Jacobian sparsity in Hessian sparsity
// (note sparsty for f''(x) * R same as for R)
if( s[0] )
v[0] = set_union(v[0], r[0] );

return true;
}

End Class Definition

}; // End of atomic_reciprocal class
}  // End empty namespace


Use Atomic Function
bool reciprocal(void)
{     bool ok = true;
double eps = 10. * CppAD::numeric_limits<double>::epsilon();

Constructor

// --------------------------------------------------------------------
// Create the atomic reciprocal object
atomic_reciprocal afun("atomic_reciprocal");


Recording
     // Create the function f(x)
//
// domain space vector
size_t n  = 1;
double  x0 = 0.5;
ax[0]     = x0;

// declare independent variables and start tape recording

// range space vector
size_t m = 1;

// call user function and store reciprocal(x) in au[0]
afun(ax, au);        // u = 1 / x

// now use AD division to invert to invert the operation
ay[0] = 1.0 / au[0]; // y = 1 / u = x

// create f: x -> y and stop tape recording
f.Dependent (ax, ay);  // f(x) = x

forward
     // 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);

reverse
     // 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);

for_sparse_jac
     // 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]

rev_sparse_jac
     // 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]

rev_sparse_hes
     // Hessian sparsity (using previous ForSparseJac call)
}