Prev Next Index-> contents reference index search external Up-> CppAD AD ADValued atomic atomic_base atomic_for_sparse_jac atomic_for_sparse_jac.cpp 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_for_sparse_jac-> atomic_for_sparse_jac.cpp atomic_for_sparse_jac.cpp Headings-> Purpose function Start Class Definition Constructor forward for_sparse_jac Use Atomic Function Test with x_1 Both a Variable and a Parameter

$\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}} }$
Atomic Forward Jacobian Sparsity: Example and Test

Purpose
This example demonstrates calculation of the forward Jacobian sparsity pattern for an atomic operation.

function
For this example, the atomic function $f : \B{R}^3 \rightarrow \B{R}^2$ is defined by $$f(x) = \left( \begin{array}{c} x_2 * x_2 \\ x_0 * x_1 \end{array} \right)$$ The corresponding Jacobian is $$f^{(1)} (x) = \left( \begin{array}{ccc} 0 & 0 & 2 x_2 \\ x_1 & x_0 & 0 \end{array} \right)$$

Start Class Definition
# include <cppad/cppad.hpp>
namespace {          // isolate items below to this file
using CppAD::vector; // abbreviate as vector
//
class atomic_for_sparse_jac : public CppAD::atomic_base<double> {

Constructor
public:
// constructor (could use const char* for name)
atomic_for_sparse_jac(const std::string& name) :
// this example only uses pack sparsty 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 == 3 );
assert( m == 2 );

// return flag
bool ok = q == 0;
if( ! ok )
return ok;

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

// Order zero forward mode.
// This case must always be implemented
// f(x) = [ x_2 * x_2 ]
//        [ x_0 * x_1 ]
assert( p <= 0 );
if( p <= 0 )
{     ty[0] = tx[2] * tx[2];
ty[1] = tx[0] * tx[1];
}
return ok;
}

for_sparse_jac
     // forward Jacobian sparsity routine called by CppAD
virtual bool for_sparse_jac(
size_t                     q ,
const vector<double>&      x )
{     // This function needed because we are using ForSparseJac
# ifndef NDEBUG
size_t n = r.size() / q;
size_t m = s.size() / q;
# endif
assert( x.size() == n );
assert( n == 3 );
assert( m == 2 );

// f'(x) = [   0,   0, 2 x_2 ]
//         [ x_1, x_0,     0 ]

// sparsity for first row of S(x) = f'(x) * R
size_t i = 0;
for(size_t j = 0; j < q; j++)
s[ i * q + j ] = r[ 2 * q + j ];

// sparsity for second row of S(x) = f'(x) * R
i = 1;
for(size_t j = 0; j < q; j++)
s[ i * q + j ] = r[ 0 * q + j ] | r[ 1 * q + j];

return true;
}
}; // End of atomic_for_sparse_jac class

Use Atomic Function
bool use_atomic_for_sparse_jac(bool x_1_variable)
{     bool ok = true;
double eps = 10. * CppAD::numeric_limits<double>::epsilon();
//
// Create the atomic for_sparse_jac object
atomic_for_sparse_jac afun("atomic_for_sparse_jac");
//
// Create the function f(u)
//
// domain space vector
size_t n  = 3;
double x_0 = 1.00;
double x_1 = 2.00;
double x_2 = 3.00;
au[0] = x_0;
au[1] = x_1;
au[2] = x_2;

// declare independent variables and start tape recording

// range space vector
size_t m = 2;

// call user function
ax[0] = au[0];
ax[2] = au[2];
if( x_1_variable )
ax[1] = au[1];
else
ax[1] = x_1;
afun(ax, ay);          // y = [ x_2 * x_2 ,  x_0 * x_1 ]^T

// create f: u -> y and stop tape recording
f.Dependent (au, ay);  // f(u) = y
//
// check function value
double check = x_2 * x_2;
ok &= NearEqual( Value(ay[0]) , check,  eps, eps);
check = x_0 * x_1;
ok &= NearEqual( Value(ay[1]) , check,  eps, eps);

// check zero order forward mode
size_t q;
vector<double> xq(n), yq(m);
q     = 0;
xq[0] = x_0;
xq[1] = x_1;
xq[2] = x_2;
yq    = f.Forward(q, xq);
check = x_2 * x_2;
ok &= NearEqual(yq[0] , check,  eps, eps);
check = x_0 * x_1;
ok &= NearEqual(yq[1] , check,  eps, eps);

// forward sparse Jacobian
CppAD::vectorBool r(n * n), s(m * n);
// r = identity matrix
for(size_t i = 0; i < n; i++)
for(size_t j = 0; j < n; j++)
r[ i * n + j] = i == j;
s = f.ForSparseJac(n, r);

// check result
check_s[ 0 * n + 0 ] = false;
check_s[ 0 * n + 1 ] = false;
check_s[ 0 * n + 2 ] = true;
check_s[ 1 * n + 0 ] = true;
check_s[ 1 * n + 1 ] = x_1_variable;
check_s[ 1 * n + 2 ] = false;
//
for(size_t i = 0; i < m * n; i++)
ok &= s[ i ] == check_s[ i ];
//
return ok;
}
}  // End empty namespace

Test with x_1 Both a Variable and a Parameter
bool for_sparse_jac(void)
{     bool ok = true;
// test with x_1 a variable
ok     &= use_atomic_for_sparse_jac(true);
// test with x_1 a parameter
ok     &= use_atomic_for_sparse_jac(false);
return ok;
}

Input File: example/atomic/for_sparse_jac.cpp