Prev Next mul_level.cpp

Multiple Tapes: Example and Test

Purpose
In this example, we use AD< AD<double> > (level two taping), the compute values of the function  f : \R^n \rightarrow \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 \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 HesTimesDir.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>. 
 

# include <cppad/cppad.hpp>

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

		sum  = 0.;
		size_t i = x.size();
		for(i = 0; i < 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>    A1_double;  // for one level of taping
	typedef CppAD::AD<A1_double> A2_double;  // 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<A1_double>  a1_x(n), a1_v(n), a1_dy(1) ;
	CPPAD_TEST_VECTOR<A2_double>  a2_x(n), a2_y(1);

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

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

	// Declare a1_F as the corresponding ADFun< AD<double> > value a2_y
	// (make sure we do not run zero order forward during constructor)
	CppAD::ADFun<A1_double> a1_F;
	a1_F.Dependent(a2_x, a2_y);

	// Values for the independent variables while taping f'(x) * v
	// Declare the independent variable for taping f'(x) * v
	CppAD::Independent(a1_x); 
	// set the argument value x for computing f'(x) * v
	a1_F.Forward(0, a1_x);
	// compute f'(x) * v
	for(j = 0; j < n; j++)
		a1_v[j] = double(n - j);
	a1_dy = a1_F.Forward(1, a1_v); 

	// declare dF as ADFun<double> function corresponding to f'(x) * v
	CppAD::ADFun<double> dF;
	dF.Dependent(a1_x, a1_dy);

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

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

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

	for(j = 0; j < n; j++)
		ok &= CppAD::NearEqual(dw[j], a1_v[j], 1e-10, 1e-10);

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