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.
# include <cppad/cppad.hpp>
namespace { // isolate items below to this fileusing CppAD::vector; // abbreviate as vector//// a utility to compute the union of two sets.using CppAD::set_union;
//class atomic_reciprocal : public CppAD::atomic_base<double> {
public:
// constructor (could use const char* for name)atomic_reciprocal(const std::string& name) :
// this exmaple only uses set sparsity patterns
CppAD::atomic_base<double>(name, atomic_base<double>::set_sparsity_enum)
{ }
private:
// forward Jacobian set sparsity routine called by CppADvirtual 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();
# endifassert( 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];
returntrue;
}
// reverse Jacobian set sparsity routine called by CppADvirtual 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();
# endifassert( 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];
returntrue;
}
// reverse Hessian set sparsity routine called by CppADvirtual 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();
# endifassert( 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] );
returntrue;
}
// 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);
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
CppAD::ADFun<double> f;
f.Dependent (ax, ay); // f(x) = x