C++ CppAD类代码示例

bool interp_retape(void)
{	bool ok = true;

	using CppAD::AD;
	using CppAD::NearEqual;

	// domain space vector
	size_t n = 1;
	CPPAD_TESTVECTOR(AD<double>) X(n);

	// loop over argument values
	size_t k;
	for(k = 0; k < TableLength - 1; k++)
	{
		X[0] = .4 * ArgumentValue[k] + .6 * ArgumentValue[k+1];

		// declare independent variables and start tape recording
		// (use a different tape for each argument value)
		CppAD::Independent(X);

		// evaluate piecewise linear interpolant at X[0]
		AD<double> A = Argument(X[0]);
		AD<double> F = Function(X[0]);
		AD<double> S = Slope(X[0]);
		AD<double> I = F + (X[0] - A) * S;

		// range space vector
		size_t m = 1;
		CPPAD_TESTVECTOR(AD<double>) Y(m);
		Y[0] = I;

		// create f: X -> Y and stop tape recording
		CppAD::ADFun<double> f(X, Y);

		// vectors for arguments to the function object f
		CPPAD_TESTVECTOR(double) x(n);   // argument values
		CPPAD_TESTVECTOR(double) y(m);   // function values
		CPPAD_TESTVECTOR(double) dx(n);  // differentials in x space
		CPPAD_TESTVECTOR(double) dy(m);  // differentials in y space

		// to check function value we use the fact that X[0] is between
		// ArgumentValue[k] and ArgumentValue[k+1]
		double delta, check;
		x[0]   = Value(X[0]);
		delta  = ArgumentValue[k+1] - ArgumentValue[k];
		check  = FunctionValue[k+1] * (x[0]-ArgumentValue[k]) / delta
	               + FunctionValue[k] * (ArgumentValue[k+1]-x[0]) / delta;
		ok    &= NearEqual(Y[0], check, 1e-10, 1e-10);

		// evaluate partials w.r.t. x[0]
		dx[0] = 1.;
		dy    = f.Forward(1, dx);

		// check that the derivative is the slope
		check = (FunctionValue[k+1] - FunctionValue[k])
		      / (ArgumentValue[k+1] - ArgumentValue[k]);
		ok   &= NearEqual(dy[0], check, 1e-10, 1e-10);
	}
	return ok;
}

bool opt_val_hes(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;

	// temporary indices
	size_t j, k;

	// x space vector
	size_t n = 1;
	BaseVector x(n);
	x[0] = 2. * 3.141592653;

	// y space vector
	size_t m = 1;
	BaseVector y(m);
	y[0] = 1.;

	// t and z vectors
	size_t ell = 10;
	BaseVector t(ell);
	BaseVector z(ell);
	for(k = 0; k < ell; k++)
	{	t[k] = double(k) / double(ell);       // time of measurement
		z[k] = y[0] * sin( x[0] * t[k] );     // data without noise
	}

	// construct the function object
	Fun fun(t, z);

	// evaluate the Jacobian and Hessian
	BaseVector jac(n), hes(n * n);
	int signdet = CppAD::opt_val_hes(x, y, fun, jac, hes);

	// we know that F_yy is positive definate for this case
	assert( signdet == 1 );

	// create ADFun object g corresponding to V(x)
	ADVector a_x(n), a_v(1);
	for(j = 0; j < n; j++)
		a_x[j] = x[j];
	Independent(a_x);
	a_v[0] = V(a_x, t, z);
	CppAD::ADFun<double> g(a_x, a_v);

	// accuracy for checks
	double eps = 10. * CppAD::numeric_limits<double>::epsilon();

	// check Jacobian
	BaseVector check_jac = g.Jacobian(x);
	for(j = 0; j < n; j++)
		ok &= NearEqual(jac[j], check_jac[j], eps, eps);

	// check Hessian
	BaseVector check_hes = g.Hessian(x, 0);
	for(j = 0; j < n*n; j++)
		ok &= NearEqual(hes[j], check_hes[j], eps, eps);

	return ok;
}

	/*!
	Fetch the next operator during a forward sweep.

	Use forward_start to initialize to the first operator; i.e.,
	the BeginOp at the beginning of the recording. 
	We use the notation forward_routine to denote the set
	forward_start, forward_next, forward_csum, forward_cskip.

	\param op [in,out]
	The input value of \c op must be its output value from the 
	previous call to a forward_routine.
	Its output value is the next operator in the recording.
	For speed, \c forward_next does not check for the special cases
	where  <tt>op == CSumOp</tt> or <tt>op == CSkipOp</tt>. In these cases, 
	the other return values from \c forward_next must be corrected by a call 
	to \c forward_csum or \c forward_cskip respectively.

	\param op_arg [in,out]
	The input value of \c op_arg must be its output value form the
	previous call to a forward routine. 
	Its output value is the
	beginning of the vector of argument indices for this operation.

	\param op_index [in,out]
	The input value of \c op_index must be its output value form the
	previous call to a forward routine. 
	Its output value is the index of the next operator in the recording. 
	Thus the ouput value following the previous call to forward_start is one.
	In addition,
	the output value increases by one with each call to forward_next. 

	\param var_index [in,out]
	The input value of \c var_index must be its output value form the
	previous call to a forward routine. 
	Its output value is the
	index of the primary (last) result corresponding to the operator op.
	*/
	void forward_next(
	OpCode& op, const addr_t*& op_arg, size_t& op_index, size_t& var_index)
	{	using CppAD::NumRes;
		using CppAD::NumArg;
		CPPAD_ASSERT_UNKNOWN( op_       == op );
		CPPAD_ASSERT_UNKNOWN( op_arg    == op_arg_ );
		CPPAD_ASSERT_UNKNOWN( op_index  == op_index_ );
		CPPAD_ASSERT_UNKNOWN( var_index == var_index_ );

		// index for the next operator 
		op_index    = ++op_index_;

		// first argument for next operator 
		op_arg      = op_arg_    += NumArg(op_);

		// next operator
		op          = op_         = OpCode( op_rec_[ op_index_ ] );

		// index for last result for next operator
		var_index   = var_index_ += NumRes(op);
		

		CPPAD_ASSERT_UNKNOWN( op_arg_rec_.data() <= op_arg_ );
		CPPAD_ASSERT_UNKNOWN( 
			op_arg_ + NumArg(op) <= op_arg_rec_.data() + op_arg_rec_.size() 
		);
		CPPAD_ASSERT_UNKNOWN( var_index_  < num_var_rec_ );
	}

	/*!
	Correct \c forward_next return values when <tt>op == CSkipOp</tt>.

	\param op [in]
	The input value of op must be the return value from the previous
	call to \c forward_next and must be \c CSkipOp. It is not modified.

	\param op_arg [in,out]
	The input value of \c op_arg must be the return value from the 
	previous call to \c forward_next. Its output value is the
	beginning of the vector of argument indices for the next operation.

	\param op_index [in]
	The input value of \c op_index does must be the return value from the
	previous call to \c forward_next. Its is not modified.

	\param var_index [in,out]
	The input value of \c var_index must be the return value from the
	previous call to \c forward_next. It is not modified.
	*/
	void forward_cskip(
	OpCode& op, const addr_t*& op_arg, size_t& op_index, size_t& var_index)
	{	using CppAD::NumRes;
		using CppAD::NumArg;
		CPPAD_ASSERT_UNKNOWN( op_       == op );
		CPPAD_ASSERT_UNKNOWN( op_arg    == op_arg_ );
		CPPAD_ASSERT_UNKNOWN( op_index  == op_index_ );
		CPPAD_ASSERT_UNKNOWN( var_index == var_index_ );

		CPPAD_ASSERT_UNKNOWN( op == CSkipOp );
		CPPAD_ASSERT_UNKNOWN( NumArg(CSkipOp) == 0 );
		CPPAD_ASSERT_UNKNOWN(
		op_arg[4] + op_arg[5] == op_arg[ 6 + op_arg[4] + op_arg[5] ]
		);
		/*
		The only thing that really needs fixing is op_arg_.
		Actual number of arugments for this operator is
			7 + op_arg[4] + op_arg[5]
 		We must change op_arg_ so that when you add NumArg(CSkipOp)
		you get first argument for next operator in sequence.
		*/
		op_arg  = op_arg_  += 7 + op_arg[4] + op_arg[5];

		CPPAD_ASSERT_UNKNOWN( op_arg_rec_.data() <= op_arg_ );
		CPPAD_ASSERT_UNKNOWN( 
			op_arg_ + NumArg(op) <= op_arg_rec_.data() + op_arg_rec_.size()
		);
		CPPAD_ASSERT_UNKNOWN( var_index_  < num_var_rec_ );
	}

	/*!
	Fetch the next operator during a reverse sweep.

	Use reverse_start to initialize to reverse play back.
	The first call to reverse_next (after reverse_start) will give the 
	last operator in the recording.
	We use the notation reverse_routine to denote the set
	reverse_start, reverse_next, reverse_csum, reverse_cskip.

	\param op [in,out]
	The input value of \c op must be its output value from the 
	previous call to a reverse_routine.
	Its output value is the next operator in the recording (in reverse order).
	The last operator sets op equal to EndOp.

	\param op_arg [in,out]
	The input value of \c op_arg must be its output value from the 
	previous call to a reverse_routine.
	Its output value is the
	beginning of the vector of argument indices for this operation.
	The last operator sets op_arg equal to the beginning of the 
	argument indices for the entire recording.
	For speed, \c reverse_next does not check for the special cases
	<tt>op == CSumOp</tt> or <tt>op == CSkipOp</tt>. In these cases, the other
	return values from \c reverse_next must be corrected by a call to 
	\c reverse_csum or \c reverse_cskip respectively.


	\param op_index [in,out]
	The input value of \c op_index must be its output value from the 
	previous call to a reverse_routine.
	Its output value
	is the index of this operator in the recording. Thus the output
	value following the previous call to reverse_start is equal to 
	the number of variables in the recording minus one.
	In addition, the output value decreases by one with each call to
	reverse_next.
	The last operator sets op_index equal to 0.

	\param var_index [in,out]
	The input value of \c var_index must be its output value from the 
	previous call to a reverse_routine.
	Its output value is the
	index of the primary (last) result corresponding to the operator op.
	The last operator sets var_index equal to 0 (corresponding to BeginOp
	at beginning of operation sequence).
	*/
	void reverse_next(
	OpCode& op, const addr_t*& op_arg, size_t& op_index, size_t& var_index)
	{	using CppAD::NumRes;
		using CppAD::NumArg;
		CPPAD_ASSERT_UNKNOWN( op_       == op );
		CPPAD_ASSERT_UNKNOWN( op_arg    == op_arg_ );
		CPPAD_ASSERT_UNKNOWN( op_index  == op_index_ );
		CPPAD_ASSERT_UNKNOWN( var_index == var_index_ );

		// index of the last result for the next operator
		CPPAD_ASSERT_UNKNOWN( var_index_ >= NumRes(op_) );
		var_index   = var_index_ -= NumRes(op_);

		// next operator
		CPPAD_ASSERT_UNKNOWN( op_index_  > 0 );
		op_index    = --op_index_;                                  // index
		op          = op_         = OpCode( op_rec_[ op_index_ ] ); // value

		// first argument for next operator
		op_arg      = op_arg_    -= NumArg(op);
		CPPAD_ASSERT_UNKNOWN( op_arg_rec_.data() <= op_arg_ );
		CPPAD_ASSERT_UNKNOWN( 
			op_arg_ + NumArg(op) <= op_arg_rec_.data() + op_arg_rec_.size() 
		);
	}

bool BenderQuad(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;

	// temporary indices
	size_t i, j;

	// x space vector
	size_t n = 1;
	BAvector x(n);
	x[0] = 2. * 3.141592653;

	// y space vector
	size_t m = 1;
	BAvector y(m);
	y[0] = 1.;

	// t and z vectors
	size_t N = 10;
	BAvector t(N);
	BAvector z(N);
	for(i = 0; i < N; i++)
	{	t[i] = double(i) / double(N);       // time of measurement
		z[i] = y[0] * sin( x[0] * t[i] );   // data without noise
	}

	// construct the function object
	Fun fun(t, z);

	// evaluate the G(x), G'(x) and G''(x)
	BAvector g(1), gx(n), gxx(n * n);
	CppAD::BenderQuad(x, y, fun, g, gx, gxx);


	// create ADFun object Gfun corresponding to G(x)
	ADvector a_x(n), a_g(1);
	for(j = 0; j < n; j++)
		a_x[j] = x[j];
	Independent(a_x);
	a_g[0] = G(a_x, t, z);
	CppAD::ADFun<double> Gfun(a_x, a_g);

	// accuracy for checks
	double eps = 10. * CppAD::numeric_limits<double>::epsilon();

	// check Jacobian
	BAvector check_gx = Gfun.Jacobian(x);
	for(j = 0; j < n; j++)
		ok &= NearEqual(gx[j], check_gx[j], eps, eps);

	// check Hessian
	BAvector check_gxx = Gfun.Hessian(x, 0);
	for(j = 0; j < n*n; j++)
		ok &= NearEqual(gxx[j], check_gxx[j], eps, eps);

	return ok;
}

bool forward_order(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;
	double eps = 10. * std::numeric_limits<double>::epsilon();

	// domain space vector
	size_t n = 2;
	CPPAD_TESTVECTOR(AD<double>) ax(n);
	ax[0] = 0.;
	ax[1] = 1.;

	// declare independent variables and starting recording
	CppAD::Independent(ax);

	// range space vector
	size_t m = 1;
	CPPAD_TESTVECTOR(AD<double>) ay(m);
	ay[0] = ax[0] * ax[0] * ax[1];

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f(ax, ay);

	// initially, the variable values during taping are stored in f
	ok &= f.size_order() == 1;

	// Compute three forward orders at one
	size_t q = 2, q1 = q+1;
	CPPAD_TESTVECTOR(double) xq(n*q1), yq;
	xq[q1*0 + 0] = 3.;    xq[q1*1 + 0] = 4.; // x^0 (order zero)
	xq[q1*0 + 1] = 1.;    xq[q1*1 + 1] = 0.; // x^1 (order one)
	xq[q1*0 + 2] = 0.;    xq[q1*1 + 2] = 0.; // x^2 (order two)
	// X(t) =   x^0 + x^1 * t + x^2 * t^2
	//      = [ 3 + t, 4 ]
	yq  = f.Forward(q, xq);
	ok &= size_t( yq.size() ) == m*q1;
	// Y(t) = F[X(t)]
	//      = (3 + t) * (3 + t) * 4
	//      = y^0 + y^1 * t + y^2 * t^2 + o(t^3)
	//
	// check y^0 (order zero)
	CPPAD_TESTVECTOR(double) x0(n);
	x0[0] = xq[q1*0 + 0];
	x0[1] = xq[q1*1 + 0];
	ok  &= NearEqual(yq[q1*0 + 0] , x0[0]*x0[0]*x0[1], eps, eps);
	//
	// check y^1 (order one)
	ok  &= NearEqual(yq[q1*0 + 1] , 2.*x0[0]*x0[1], eps, eps);
	//
	// check y^2 (order two)
	double F_00 = 2. * yq[q1*0 + 2]; // second partial F w.r.t. x_0, x_0
	ok   &= NearEqual(F_00, 2.*x0[1], eps, eps);

	// check number of orders per variable
	ok   &= f.size_order() == 3;

	return ok;
}

bool mul_level(void)
{	bool ok = true;
	using CppAD::checkpoint;
	using CppAD::ADFun;
	using CppAD::Independent;

	// domain dimension for this problem
	size_t n = 10;
	size_t m = 1;

	// checkpoint version of the function F(x)
	a2vector a2x(n), a2y(m);
	for(size_t j = 0; j < n; j++)
		a2x[j] = a2double(j + 1);
	//
	// could also use bool_sparsity_enum or set_sparsity_enum
	checkpoint<a1double> atom_f("atom_f", f_algo, a2x, a2y);
	//
	// Record a version of y = f(x) without checkpointing
	Independent(a2x);
	f_algo(a2x, a2y);
	ADFun<a1double> check_not(a2x, a2y);
	//
	// number of variables in a tape of f_algo that does not use checkpointing
	size_t size_not = check_not.size_var();
	//
	// Record a version of y = f(x) with checkpointing
	Independent(a2x);
	atom_f(a2x, a2y);
	ADFun<a1double> check_yes(a2x, a2y);
	//
	// f_algo is represented by one atomic operation in this tape
	ok &= check_yes.size_var() < size_not;
	//
	// now record operations at a1double level
	a1vector a1x(n), a1y(m);
	for(size_t j = 0; j < n; j++)
		a1x[j] = a1double(j + 1);
	//
	// without checkpointing
	Independent(a1x);
	a1y = check_not.Forward(0, a1x);
	ADFun<double> with_not(a1x, a1y);
	//
	// should have the same size
	ok &= with_not.size_var() == size_not;
	//
	// with checkpointing
	Independent(a1x);
	a1y = check_yes.Forward(0, a1x);
	ADFun<double> with_yes(a1x, a1y);
	//
	// f_algo is nolonger represented by one atomic operation in this tape
	ok &= with_yes.size_var() == size_not;
	//
	return ok;
}

bool atanh(void)
{   bool ok = true;

    using CppAD::AD;
    using CppAD::NearEqual;

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

    // domain space vector
    size_t n  = 1;
    double x0 = 0.5;
    CPPAD_TESTVECTOR(AD<double>) ax(n);
    ax[0]     = x0;

    // declare independent variables and start tape recording
    CppAD::Independent(ax);

    // a temporary value
    AD<double> tanh_of_x0 = CppAD::tanh(ax[0]);

    // range space vector
    size_t m = 1;
    CPPAD_TESTVECTOR(AD<double>) ay(m);
    ay[0] = CppAD::atanh(tanh_of_x0);

    // create f: x -> y and stop tape recording
    CppAD::ADFun<double> f(ax, ay);

    // check value
    ok &= NearEqual(ay[0] , x0,  eps, eps);

    // forward computation of first partial w.r.t. x[0]
    CPPAD_TESTVECTOR(double) dx(n);
    CPPAD_TESTVECTOR(double) dy(m);
    dx[0] = 1.;
    dy    = f.Forward(1, dx);
    ok   &= NearEqual(dy[0], 1., eps, eps);

    // forward computation of higher order partials w.r.t. x[0]
    size_t n_order = 5;
    for(size_t order = 2; order < n_order; order++)
    {   dx[0] = 0.;
        dy    = f.Forward(order, dx);
        ok   &= NearEqual(dy[0], 0., eps, eps);
    }
    // reverse computation of derivatives
    CPPAD_TESTVECTOR(double)  w(m);
    CPPAD_TESTVECTOR(double) dw(n_order * n);
    w[0]  = 1.;
    dw    = f.Reverse(n_order, w);
    ok   &= NearEqual(dw[0], 1., eps, eps);
    for(size_t order = 1; order < n_order; order++)
        ok   &= NearEqual(dw[order * n + 0], 0., eps, eps);

    return ok;
}

bool atan2(void)
{	bool ok = true;

	using CppAD::AD;
	using CppAD::NearEqual;
	double eps99 = 99.0 * std::numeric_limits<double>::epsilon();

	// domain space vector
	size_t n  = 1;
	double x0 = 0.5;
	CPPAD_TESTVECTOR(AD<double>) x(n);
	x[0]      = x0;

	// declare independent variables and start tape recording
	CppAD::Independent(x);

	// a temporary value
	AD<double> sin_of_x0 = CppAD::sin(x[0]);
	AD<double> cos_of_x0 = CppAD::cos(x[0]);

	// range space vector
	size_t m = 1;
	CPPAD_TESTVECTOR(AD<double>) y(m);
	y[0] = CppAD::atan2(sin_of_x0, cos_of_x0);

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f(x, y);

	// check value
	ok &= NearEqual(y[0] , x0, eps99, eps99);

	// forward computation of first partial w.r.t. x[0]
	CPPAD_TESTVECTOR(double) dx(n);
	CPPAD_TESTVECTOR(double) dy(m);
	dx[0] = 1.;
	dy    = f.Forward(1, dx);
	ok   &= NearEqual(dy[0], 1., eps99, eps99);

	// reverse computation of derivative of y[0]
	CPPAD_TESTVECTOR(double)  w(m);
	CPPAD_TESTVECTOR(double) dw(n);
	w[0]  = 1.;
	dw    = f.Reverse(1, w);
	ok   &= NearEqual(dw[0], 1., eps99, eps99);

	// use a VecAD<Base>::reference object with atan2
	CppAD::VecAD<double> v(2);
	AD<double> zero(0);
	AD<double> one(1);
	v[zero]           = sin_of_x0;
	v[one]            = cos_of_x0;
	AD<double> result = CppAD::atan2(v[zero], v[one]);
	ok               &= NearEqual(result, x0, eps99, eps99);

	return ok;
}

bool Div(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;
	double eps99 = 99.0 * std::numeric_limits<double>::epsilon();


	// domain space vector
	size_t n  = 1;
	double x0 = 0.5;
	CPPAD_TESTVECTOR(AD<double>) x(n);
	x[0]      = x0;

	// declare independent variables and start tape recording
	CppAD::Independent(x);

	// some binary division operations
	AD<double> a = x[0] / 1.; // AD<double> / double
	AD<double> b = a  / 2;    // AD<double> / int
	AD<double> c = 3. / b;    // double     / AD<double>
	AD<double> d = 4  / c;    // int        / AD<double>

	// range space vector
	size_t m = 1;
	CPPAD_TESTVECTOR(AD<double>) y(m);
	y[0] = (x[0] * x[0]) / d;   // AD<double> / AD<double>

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f(x, y);

	// check value
	ok &= NearEqual(y[0], x0*x0*3.*2.*1./(4.*x0), eps99, eps99);

	// forward computation of partials w.r.t. x[0]
	CPPAD_TESTVECTOR(double) dx(n);
	CPPAD_TESTVECTOR(double) dy(m);
	dx[0] = 1.;
	dy    = f.Forward(1, dx);
	ok   &= NearEqual(dy[0], 3.*2.*1./4., eps99, eps99);

	// reverse computation of derivative of y[0]
	CPPAD_TESTVECTOR(double)  w(m);
	CPPAD_TESTVECTOR(double) dw(n);
	w[0]  = 1.;
	dw    = f.Reverse(1, w);
	ok   &= NearEqual(dw[0], 3.*2.*1./4., eps99, eps99);

	// use a VecAD<Base>::reference object with division
	CppAD::VecAD<double> v(1);
	AD<double> zero(0);
	v[zero] = d;
	AD<double> result = (x[0] * x[0]) / v[zero];
	ok     &= (result == y[0]);

	return ok;
}

bool AddEq(void)
{   bool ok = true;
    using CppAD::AD;
    using CppAD::NearEqual;
    double eps99 = 99.0 * std::numeric_limits<double>::epsilon();

    // domain space vector
    size_t  n = 1;
    double x0 = .5;
    CPPAD_TESTVECTOR(AD<double>) x(n);
    x[0]      = x0;

    // declare independent variables and start tape recording
    CppAD::Independent(x);

    // range space vector
    size_t m = 2;
    CPPAD_TESTVECTOR(AD<double>) y(m);
    y[0] = x[0];         // initial value
    y[0] += 2;           // AD<double> += int
    y[0] += 4.;          // AD<double> += double
    y[1] = y[0] += x[0]; // use the result of a compound assignment

    // create f: x -> y and stop tape recording
    CppAD::ADFun<double> f(x, y);

    // check value
    ok &= NearEqual(y[0] , x0+2.+4.+x0, eps99, eps99);
    ok &= NearEqual(y[1] ,        y[0], eps99, eps99);

    // forward computation of partials w.r.t. x[0]
    CPPAD_TESTVECTOR(double) dx(n);
    CPPAD_TESTVECTOR(double) dy(m);
    dx[0] = 1.;
    dy    = f.Forward(1, dx);
    ok   &= NearEqual(dy[0], 2., eps99, eps99);
    ok   &= NearEqual(dy[1], 2., eps99, eps99);

    // reverse computation of derivative of y[0]
    CPPAD_TESTVECTOR(double)  w(m);
    CPPAD_TESTVECTOR(double) dw(n);
    w[0]  = 1.;
    w[1]  = 0.;
    dw    = f.Reverse(1, w);
    ok   &= NearEqual(dw[0], 2., eps99, eps99);

    // use a VecAD<Base>::reference object with computed addition
    CppAD::VecAD<double> v(1);
    AD<double> zero(0);
    AD<double> result = 1;
    v[zero] = 2;
    result += v[zero];
    ok     &= (result == 3);

    return ok;
}

bool MulEq(void)
{   bool ok = true;
    using CppAD::AD;
    using CppAD::NearEqual;

    // domain space vector
    size_t  n = 1;
    double x0 = .5;
    CPPAD_TESTVECTOR(AD<double>) x(n);
    x[0]      = x0;

    // declare independent variables and start tape recording
    CppAD::Independent(x);

    // range space vector
    size_t m = 2;
    CPPAD_TESTVECTOR(AD<double>) y(m);
    y[0] = x[0];         // initial value
    y[0] *= 2;           // AD<double> *= int
    y[0] *= 4.;          // AD<double> *= double
    y[1] = y[0] *= x[0]; // use the result of a computed assignment

    // create f: x -> y and stop tape recording
    CppAD::ADFun<double> f(x, y);

    // check value
    ok &= NearEqual(y[0] , x0*2.*4.*x0,  1e-10 , 1e-10);
    ok &= NearEqual(y[1] ,        y[0],  1e-10 , 1e-10);

    // forward computation of partials w.r.t. x[0]
    CPPAD_TESTVECTOR(double) dx(n);
    CPPAD_TESTVECTOR(double) dy(m);
    dx[0] = 1.;
    dy    = f.Forward(1, dx);
    ok   &= NearEqual(dy[0], 8.*2.*x0, 1e-10, 1e-10);
    ok   &= NearEqual(dy[1], 8.*2.*x0, 1e-10, 1e-10);

    // reverse computation of derivative of y[0]
    CPPAD_TESTVECTOR(double)  w(m);
    CPPAD_TESTVECTOR(double) dw(n);
    w[0]  = 1.;
    w[1]  = 0.;
    dw    = f.Reverse(1, w);
    ok   &= NearEqual(dw[0], 8.*2.*x0, 1e-10, 1e-10);

    // use a VecAD<Base>::reference object with computed multiplication
    CppAD::VecAD<double> v(1);
    AD<double> zero(0);
    AD<double> result = 1;
    v[zero] = 2;
    result *= v[zero];
    ok     &= (result == 2);

    return ok;
}

/* $$
$head Use Atomic Function$$
$codep */
bool get_started(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;
	double eps = 10. * CppAD::numeric_limits<double>::epsilon();
/* $$
$subhead Constructor$$
$codep */
	// Create the atomic get_started object
	atomic_get_started afun("atomic_get_started");
/* $$
$subhead Recording$$
$codep */
	// Create the function f(x)
	//
	// domain space vector
	size_t n  = 1;
	double  x0 = 0.5;
	vector< AD<double> > ax(n);
	ax[0]     = x0;

	// declare independent variables and start tape recording
	CppAD::Independent(ax);

	// range space vector 
	size_t m = 1;
	vector< AD<double> > ay(m);

	// call user function and store get_started(x) in au[0] 
	vector< AD<double> > au(m);
	afun(ax, au);        // u = 1 / x

	// now use AD division to invert to invert the operation
	ay[0] = 1.0 / au[0]; // y = 1 / u = x

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f;
	f.Dependent (ax, ay);  // f(x) = x
/* $$
$subhead forward$$
$codep */
	// check function value 
	double check = x0;
	ok &= NearEqual( Value(ay[0]) , check,  eps, eps);

	// check zero order forward mode
	size_t p;
	vector<double> x_p(n), y_p(m);
	p      = 0;
	x_p[0] = x0;
	y_p    = f.Forward(p, x_p);
	ok &= NearEqual(y_p[0] , check,  eps, eps);

	return ok;
}

bool Add(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;
	double eps99 = 99.0 * std::numeric_limits<double>::epsilon();

	// domain space vector
	size_t n  = 1;
	double x0 = 0.5;
	CPPAD_TESTVECTOR(AD<double>) x(n);
	x[0]      = x0;

	// declare independent variables and start tape recording
	CppAD::Independent(x);

	// some binary addition operations
	AD<double> a = x[0] + 1.; // AD<double> + double
	AD<double> b = a    + 2;  // AD<double> + int
	AD<double> c = 3.   + b;  // double     + AD<double>
	AD<double> d = 4    + c;  // int        + AD<double>

	// range space vector
	size_t m = 1;
	CPPAD_TESTVECTOR(AD<double>) y(m);
	y[0] = d + x[0];          // AD<double> + AD<double>

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f(x, y);

	// check value
	ok &= NearEqual(y[0] , 2. * x0 + 10, eps99, eps99);

	// forward computation of partials w.r.t. x[0]
	CPPAD_TESTVECTOR(double) dx(n);
	CPPAD_TESTVECTOR(double) dy(m);
	dx[0] = 1.;
	dy    = f.Forward(1, dx);
	ok   &= NearEqual(dy[0], 2., eps99, eps99);

	// reverse computation of derivative of y[0]
	CPPAD_TESTVECTOR(double)  w(m);
	CPPAD_TESTVECTOR(double) dw(n);
	w[0]  = 1.;
	dw    = f.Reverse(1, w);
	ok   &= NearEqual(dw[0], 2., eps99, eps99);

	// use a VecAD<Base>::reference object with addition
	CppAD::VecAD<double> v(1);
	AD<double> zero(0);
	v[zero] = a;
	AD<double> result = v[zero] + 2;
	ok     &= (result == b);

	return ok;
}