Prev Next mul_level.cpp

@(@\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}} }@)@ This is cppad-20221105 documentation. Here is a link to its current documentation .
Multiple Level of AD: Example and Test

See Also
base2ad.cpp

Purpose
In this example, we use AD< AD<double> > (level two taping), the compute values of the function @(@ f : \B{R}^n \rightarrow \B{R} @)@ where @[@ f(x) = \frac{1}{2} \left( x_0^2 + \cdots + x_{n-1}^2 \right) @]@ We then use AD<double> (level one taping) to compute the directional derivative @[@ f^{(1)} (x) * v = x_0 v_0 + \cdots + x_{n-1} v_{n-1} @]@. where @(@ v \in \B{R}^n @)@. We then use double (no taping) to compute @[@ \frac{d}{dx} \left[ f^{(1)} (x) * v \right] = v @]@ This is only meant as an example of multiple levels of taping. The example hes_times_dir.cpp computes the same value more efficiently by using the identity: @[@ \frac{d}{dx} \left[ f^{(1)} (x) * v \right] = f^{(2)} (x) * v @]@ The example mul_level_adolc.cpp computes the same values using Adolc's type adouble and CppAD's type AD<adouble>.

Source 


# include <cppad/cppad.hpp>

namespace {
    // f(x) = |x|^2 / 2 = .5 * ( x[0]^2 + ... + x[n-1]^2 )
    template <class Type>
    Type f(const CPPAD_TESTVECTOR(Type)& x)
    {   Type sum;

        sum  = 0.;
        size_t i = size_t(x.size());
        for(i = 0; i < size_t(x.size()); i++)
            sum += x[i] * x[i];

        return .5 * sum;
    }
}

bool mul_level(void)
{   bool ok = true;                          // initialize test result

    typedef CppAD::AD<double>   a1type;    // for one level of taping
    typedef CppAD::AD<a1type>    a2type;    // for two levels of taping
    size_t n = 5;                           // dimension for example
    size_t j;                               // a temporary index variable

    // 10 times machine epsilon
    double eps = 10. * std::numeric_limits<double>::epsilon();

    CPPAD_TESTVECTOR(double) x(n);
    CPPAD_TESTVECTOR(a1type)  a1x(n), a1v(n), a1dy(1) ;
    CPPAD_TESTVECTOR(a2type)  a2x(n), a2y(1);

    // Values for the independent variables while taping the function f(x)
    for(j = 0; j < n; j++)
        a2x[j] = a1x[j] = x[j] = double(j);
    // Declare the independent variable for taping f(x)
    CppAD::Independent(a2x);

    // Use AD< AD<double> > to tape the evaluation of f(x)
    a2y[0] = f(a2x);

    // Declare a1f as the corresponding ADFun< AD<double> >
    // (make sure we do not run zero order forward during constructor)
    CppAD::ADFun<a1type> a1f;
    a1f.Dependent(a2x, a2y);

    // Values for the independent variables while taping f'(x) * v
    // Declare the independent variable for taping f'(x) * v
    // (Note we did not have to tape the creationg of a1f.)
    CppAD::Independent(a1x);

    // set the argument value x for computing f'(x) * v
    a1f.Forward(0, a1x);
    // compute f'(x) * v
    for(j = 0; j < n; j++)
        a1v[j] = double(n - j);
    a1dy = a1f.Forward(1, a1v);

    // declare g as ADFun<double> function corresponding to f'(x) * v
    CppAD::ADFun<double> g;
    g.Dependent(a1x, a1dy);

    // optimize out operations not necessary for function f'(x) * v
    g.optimize();

    // Evaluate f'(x) * v
    g.Forward(0, x);

    // compute the d/dx of f'(x) * v = f''(x) * v = v
    CPPAD_TESTVECTOR(double) w(1);
    CPPAD_TESTVECTOR(double) dw(n);
    w[0] = 1.;
    dw   = g.Reverse(1, w);

    for(j = 0; j < n; j++)
        ok &= CppAD::NearEqual(dw[j], a1v[j], eps, eps);

    return ok;
}
Input File: example/general/mul_level.cpp