# include <cppad/cppad.hpp>
namespace { // Begin empty namespaceusing CppAD::vector;
//// a utility to compute the union of two sets.using CppAD::set_union;
//class atomic_tangent : public CppAD::atomic_base<float> {
// reverse Hessian 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<bool>& r ,
const vector<bool>& u ,
vector<bool>& v ,
const vector<float>& x )
{
# ifndef NDEBUG
size_t m = s.size();
size_t n = t.size();
# endifassert( x.size() == n );
assert( r.size() == n * p );
assert( u.size() == m * p );
assert( v.size() == n * p );
assert( n == 1 );
assert( m == 2 );
// 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)
t[0] = s[0] | s[1];
// 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 both components// of f'(x) may be non-zero;
size_t j;
for(j = 0; j < p; j++)
v[j] = u[ 0 * p + j ] | u[ 1 * p + j ];
// include forward Jacobian sparsity in Hessian sparsity// (note sparsty for f''(x) * R same as for R)if( s[0] | s[1] )
{ for(j = 0; j < p; j++)
{ // Visual Studio 2013 generates warning without bool below
v[j] |= bool( r[j] );
}
}
returntrue;
}
// reverse Hessian 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<float>& x )
{
# ifndef NDEBUG
size_t m = s.size();
size_t n = t.size();
# endifassert( x.size() == n );
assert( r.size() == n );
assert( u.size() == m );
assert( v.size() == n );
assert( n == 1 );
assert( m == 2 );
// 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)
t[0] = s[0] | s[1];
// 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 both components// of f'(x) may be non-zero;
v[0] = set_union(u[0], u[1]);
// include forward Jacobian sparsity in Hessian sparsity// (note sparsty for f''(x) * R same as for R)if( s[0] | s[1] )
v[0] = set_union(v[0], r[0]);
returntrue;
}
// --------------------------------------------------------------------// Creater a tan and tanh object
atomic_tangent my_tan("my_tan", false), my_tanh("my_tanh", true);
// domain space vector
size_t n = 1;
float x0 = 0.5;
CppAD::vector< AD<float> > ax(n);
ax[0] = x0;
// declare independent variables and start tape recording
CppAD::Independent(ax);
// range space vector
size_t m = 3;
CppAD::vector< AD<float> > af(m);
// temporary vector for computations// (my_tan and my_tanh computes tan or tanh and its square)
CppAD::vector< AD<float> > az(2);
// call atomic tan function and store tan(x) in f[0] (ignore tan(x)^2)my_tan(ax, az);
af[0] = az[0];
// call atomic tanh function and store tanh(x) in f[1] (ignore tanh(x)^2)my_tanh(ax, az);
af[1] = az[0];
// put a constant in f[2] = tanh(1.) (for sparsity pattern testing)
CppAD::vector< AD<float> > one(1);
one[0] = 1.;
my_tanh(one, az);
af[2] = az[0];
// create f: x -> f and stop tape recording
CppAD::ADFun<float> F;
F.Dependent(ax, af);
// Forward mode computation of sparsity pattern for F.
size_t p = n;
// user vectorBool because m and n are small
CppAD::vectorBool r1(p), s1(m * p);
r1[0] = true; // propagate sparsity for x[0]
s1 = F.ForSparseJac(p, r1);
ok &= (s1[0] == true); // f[0] depends on x[0]
ok &= (s1[1] == true); // f[1] depends on x[0]
ok &= (s1[2] == false); // f[2] does not depend on x[0]
// Reverse mode computation of sparsity pattern for F.
size_t q = m;
CppAD::vectorBool s2(q * m), r2(q * n);
// Sparsity pattern for identity matrix
size_t i, j;
for(i = 0; i < q; i++)
{ for(j = 0; j < m; j++)
s2[i * q + j] = (i == j);
}
r2 = F.RevSparseJac(q, s2);
ok &= (r2[0] == true); // f[0] depends on x[0]
ok &= (r2[1] == true); // f[1] depends on x[0]
ok &= (r2[2] == false); // f[2] does not depend on x[0]
// check tanh results for a large value of x
x[0] = std::numeric_limits<float>::max() / two;
f = F.Forward(0, x);
tanh = 1.;
ok &= NearEqual(f[1], tanh, eps, eps);
df = F.Forward(1, dx);
tanhp = 0.;
ok &= NearEqual(df[1], tanhp, eps, eps);
return ok;
}