\(\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_three_rev_depend.cpp¶
View page sourceAtomic Functions Reverse Dependency Analysis: Example and Test¶
Purpose¶
This example demonstrates using atomic_three function in the definition of a function that is optimized.
Function¶
For this example, the atomic function \(g : \B{R}^3 \rightarrow \B{R}^3\) is defined by \(g_0 (x) = x_0 * x_0\), \(g_1 (x) = x_0 * x_1\), \(g_2 (x) = x_1 * x_2\).
Start Class Definition¶
# include <cppad/cppad.hpp> // CppAD include file
namespace { // start empty namespace
using CppAD::vector; // abbreviate CppAD::vector using vector
// start definition of atomic derived class using atomic_three interface
class atomic_optimize : public CppAD::atomic_three<double> {
Constructor¶
public:
// can use const char* name when calling this constructor
atomic_optimize(const std::string& name) : // can have more arguments
CppAD::atomic_three<double>(name) // inform base class of name
{ }
private:
for_type¶
// calculate type_y
bool for_type(
const vector<double>& parameter_x ,
const vector<CppAD::ad_type_enum>& type_x ,
vector<CppAD::ad_type_enum>& type_y ) override
{ assert( parameter_x.size() == type_x.size() );
bool ok = type_x.size() == 3; // n
ok &= type_y.size() == 3; // m
if( ! ok )
return false;
type_y[0] = type_x[0];
type_y[1] = std::max( type_x[0], type_x[1] );
type_y[2] = std::max( type_x[1], type_x[2] );
return true;
}
rev_depend¶
// calculate depend_x
bool rev_depend(
const vector<double>& parameter_x ,
const vector<CppAD::ad_type_enum>& type_x ,
vector<bool>& depend_x ,
const vector<bool>& depend_y ) override
{ assert( parameter_x.size() == depend_x.size() );
bool ok = depend_x.size() == 3; // n
ok &= depend_y.size() == 3; // m
if( ! ok )
return false;
depend_x[0] = depend_y[0] | depend_y[1];
depend_x[1] = depend_y[1] | depend_y[2];
depend_x[2] = depend_y[2];
return true;
}
forward¶
// forward mode routine called by CppAD
bool forward(
const vector<double>& parameter_x ,
const vector<CppAD::ad_type_enum>& type_x ,
size_t need_y ,
size_t order_low ,
size_t order_up ,
const vector<double>& taylor_x ,
vector<double>& taylor_y
) override
{
# ifndef NDEBUG
size_t n = taylor_x.size() / (order_up + 1);
size_t m = taylor_y.size() / (order_up + 1);
# endif
assert( n == 3 );
assert( m == 3 );
assert( order_low <= order_up );
// return flag
bool ok = order_up == 0;
if( ! ok )
return ok;
// Order zero forward mode.
// This case must always be implemented
if( need_y > size_t(CppAD::variable_enum) )
{ // g_0 = x_0 * x_0
taylor_y[0] = taylor_x[0] * taylor_x[0];
// g_1 = x_0 * x_1
taylor_y[1] = taylor_x[0] * taylor_x[1];
// g_2 = x_1 * x_2
taylor_y[2] = taylor_x[1] * taylor_x[2];
}
else
{ // This uses need_y to reduce amount of computation.
// It is probably faster, for this case, to ignore need_y.
vector<CppAD::ad_type_enum> type_y( taylor_y.size() );
for_type(taylor_x, type_x, type_y);
// g_0 = x_0 * x_0
if( size_t(type_y[0]) == need_y )
taylor_y[0] = taylor_x[0] * taylor_x[0];
// g_1 = x_0 * x_1
if( size_t(type_y[1]) == need_y )
taylor_y[1] = taylor_x[0] * taylor_x[1];
// g_2 = x_1 * x_2
if( size_t(type_y[2]) == need_y )
taylor_y[2] = taylor_x[1] * taylor_x[2];
}
return ok;
}
End Class Definition¶
}; // End of atomic_optimize class
} // End empty namespace
Use Atomic Function¶
bool rev_depend(void)
{ bool ok = true;
using CppAD::AD;
using CppAD::NearEqual;
double eps = 10. * CppAD::numeric_limits<double>::epsilon();
Constructor¶
// Create the atomic dynamic object corresponding to g(x)
atomic_optimize afun("atomic_optimize");
Recording¶
// Create the function f(u) = g(c, p, u) for this example.
//
// constant parameter
double c_0 = 2.0;
//
// indepndent dynamic parameter vector
size_t np = 1;
CPPAD_TESTVECTOR(double) p(np);
CPPAD_TESTVECTOR( AD<double> ) ap(np);
ap[0] = p[0] = 3.0;
//
// independent variable vector
size_t nu = 1;
double u_0 = 0.5;
CPPAD_TESTVECTOR( AD<double> ) au(nu);
au[0] = u_0;
// declare independent variables and start tape recording
CppAD::Independent(au, ap);
// range space vector
size_t ny = 3;
CPPAD_TESTVECTOR( AD<double> ) ay(ny);
// call atomic function and store result in ay
// y = ( c * c, c * p, p * u )
CPPAD_TESTVECTOR( AD<double> ) ax(3);
ax[0] = c_0; // x_0 = c
ax[1] = ap[0]; // x_1 = p
ax[2] = au[0]; // x_2 = u
afun(ax, ay);
// check type of result
ok &= Constant( ay[0] ); // c * c
ok &= Dynamic( ay[1] ); // c * p
ok &= Variable( ay[2] ); // p * u
// create f: u -> y and stop tape recording
CppAD::ADFun<double> f;
f.Dependent (au, ay); // f(u) = (c * c, c * p, p * u)
optimize¶
This operation does a callback to rev_depend defined above
f.optimize();
forward¶
// check function value
double check = c_0 * c_0;
ok &= NearEqual( Value(ay[0]) , check, eps, eps);
check = c_0 * p[0];
ok &= NearEqual( Value(ay[1]) , check, eps, eps);
check = p[0] * u_0;
ok &= NearEqual( Value(ay[2]) , check, eps, eps);
// check zero order forward mode
size_t q;
CPPAD_TESTVECTOR( double ) u_q(nu), y_q(ny);
q = 0;
u_q[0] = u_0;
y_q = f.Forward(q, u_q);
check = c_0 * c_0;
ok &= NearEqual(y_q[0] , check, eps, eps);
check = c_0 * p[0];
ok &= NearEqual(y_q[1] , check, eps, eps);
check = p[0] * u_0;
ok &= NearEqual(y_q[2] , check, eps, eps);
// set new value for dynamic parameters
p[0] = 2.0 * p[0];
f.new_dynamic(p);
y_q = f.Forward(q, u_q);
check = c_0 * c_0;
ok &= NearEqual(y_q[0] , check, eps, eps);
check = c_0 * p[0];
ok &= NearEqual(y_q[1] , check, eps, eps);
check = p[0] * u_0;
ok &= NearEqual(y_q[2] , check, eps, eps);
Return Test Result¶
return ok;
}