![]() |
Prev | Next | mul_level_adolc.cpp |
adouble
type,
together with CppAD's AD<adouble>
type,
for multiple levels of taping.
The example computes
\[
\frac{d}{dx} \left[ f^{(1)} (x) * v \right]
\]
where
f : \R^n \rightarrow \R
and
v \in \R^n
.
The example HesTimesDir.cpp
computes the same value using only
one level of taping (more efficient) and the identity
\[
\frac{d}{dx} \left[ f^{(1)} (x) * v \right] = f^{(2)} (x) * v
\]
The example mul_level.cpp
computes the same values using
AD<double>
and AD< AD<double> >
.
new
and delete
are used to allocate this memory.
The preprocessor macros
CPPAD_TRACK_NEW_VEC
and
CPPAD_TRACK_DEL_VEC
are used to check for errors in the
use of new
and delete
when the example is compiled for
debugging (when NDEBUG
is not defined).
# include <adolc/adouble.h>
# include <adolc/interfaces.h>
// adouble definitions not in Adolc distribution and
// required in order to use CppAD::AD<adouble>
# include "base_adolc.hpp"
# include <cppad/cppad.hpp>
namespace { // put this function in the empty namespace
// f(x) = |x|^2 = .5 * ( x[0]^2 + ... + x[n-1]^2 + .5 )
template <class Type>
Type f(CPPAD_TEST_VECTOR<Type> &x)
{ Type sum;
// check assignment of AD< AD<double> > = double
sum = .5;
sum += .5;
size_t i = x.size();
while(i--)
sum += x[i] * x[i];
// check computed assignment AD< AD<double> > -= int
sum -= 1;
// check double * AD< AD<double> >
return .5 * sum;
}
}
bool mul_level_adolc(void)
{ bool ok = true; // initialize test result
typedef adouble ADdouble; // for first level of taping
typedef CppAD::AD<ADdouble> ADDdouble; // for second level of taping
size_t n = 5; // number independent variables
CPPAD_TEST_VECTOR<double> x(n);
CPPAD_TEST_VECTOR<ADdouble> a_x(n);
CPPAD_TEST_VECTOR<ADDdouble> aa_x(n);
// value of the independent variables
int tag = 0; // Adolc setup
int keep = 1;
trace_on(tag, keep);
size_t j;
for(j = 0; j < n; j++)
{ x[j] = double(j); // x[j] = j
a_x[j] <<= x[j]; // a_x is independent for ADdouble
}
for(j = 0; j < n; j++)
aa_x[j] = a_x[j]; // track how aa_x depends on a_x
CppAD::Independent(aa_x); // aa_x is independent for ADDdouble
// compute function
CPPAD_TEST_VECTOR<ADDdouble> aa_f(1); // scalar valued function
aa_f[0] = f(aa_x); // has only one component
// declare inner function (corresponding to ADDdouble calculation)
CppAD::ADFun<ADdouble> a_F(aa_x, aa_f);
// compute f'(x)
size_t p = 1; // order of derivative of a_F
CPPAD_TEST_VECTOR<ADdouble> a_w(1); // weight vector for a_F
CPPAD_TEST_VECTOR<ADdouble> a_df(n); // value of derivative
a_w[0] = 1; // weighted function same as a_F
a_df = a_F.Reverse(p, a_w); // gradient of f
// declare outter function
// (corresponding to the tape of adouble operations)
double df_j;
for(j = 0; j < n; j++)
a_df[j] >>= df_j;
trace_off();
// compute the d/dx of f'(x) * v = f''(x) * v
size_t m = n; // # dependent in f'(x)
double *v, *ddf_v;
v = CPPAD_TRACK_NEW_VEC(m, v); // track v = new double[m]
ddf_v = CPPAD_TRACK_NEW_VEC(n, ddf_v); // track ddf_v = new double[n]
for(j = 0; j < n; j++)
v[j] = double(n - j);
fos_reverse(tag, int(m), int(n), v, ddf_v);
// f(x) = .5 * ( x[0]^2 + x[1]^2 + ... + x[n-1]^2 )
// f'(x) = (x[0], x[1], ... , x[n-1])
// f''(x) * v = ( v[0], v[1], ... , x[n-1] )
for(j = 0; j < n; j++)
ok &= CppAD::NearEqual(ddf_v[j], v[j], 1e-10, 1e-10);
CPPAD_TRACK_DEL_VEC(v); // check usage of delete
CPPAD_TRACK_DEL_VEC(ddf_v);
return ok;
}