Prev Next mul_level.cpp

Multiple Tapes: Example and Test

Purpose
This is an example and test of using the AD<double> type, together with the AD< AD<double> > type, for multiple levels of taping. The example computes the value  \[
     \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_adolc.cpp computes the same values using Adolc's type adouble and CppAD's type AD<adouble>.

# 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(void) 
{    bool ok = true;                          // initialize test result

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

     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
     for(j = 0; j < n; j++)
          a_x[j] = x[j] = double(j); // x[j] = j
     Independent(a_x);                  // a_x is indedendent 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 ADdouble calculation)
     CppAD::ADFun<double> df(a_x, a_df);

     // compute the d/dx of f'(x) * v = f''(x) * v
     CPPAD_TEST_VECTOR<double> v(n);
     CPPAD_TEST_VECTOR<double> ddf_v(n);
     for(j = 0; j < n; j++)
          v[j] = double(n - j);
     ddf_v = df.Reverse(p, 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);

     return ok;
}

Input File: example/mul_level.cpp