C++ CppAD::NearEqual方法代码示例

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;
}

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 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 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 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 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 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 pow_int(void)
{	bool ok = true;

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

	// declare independent variables and start tape recording
	size_t n  = 1;
	double x0 = -0.5;
	CPPAD_TESTVECTOR(AD<double>) x(n);
	x[0]      = x0;
	CppAD::Independent(x);

	// dependent variable vector
	size_t m = 7;
	CPPAD_TESTVECTOR(AD<double>) y(m);
	int i;
	for(i = 0; i < int(m); i++)
		y[i] = CppAD::pow(x[0], i - 3);

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

	// check value
	double check;
	for(i = 0; i < int(m); i++)
	{	check = std::pow(x0, double(i - 3));
		ok &= NearEqual(y[i] , check,  1e-10 , 1e-10);
	}

	// 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);
	for(i = 0; i < int(m); i++)
	{	check = double(i-3) * std::pow(x0, double(i - 4));
		ok &= NearEqual(dy[i] , check,  1e-10 , 1e-10);
	}

	// reverse computation of derivative of y[i]
	CPPAD_TESTVECTOR(double)  w(m);
	CPPAD_TESTVECTOR(double) dw(n);
	for(i = 0; i < int(m); i++)
		w[i] = 0.;
	for(i = 0; i < int(m); i++)
	{	w[i] = 1.;
		dw    = f.Reverse(1, w);
		check = double(i-3) * std::pow(x0, double(i - 4));
		ok &= NearEqual(dw[0] , check,  1e-10 , 1e-10);
		w[i] = 0.;
	}

	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;
}

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

	// domain space vector
	size_t  n = 1;
	double x0 = .5;
	CPPAD_TEST_VECTOR< 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_TEST_VECTOR< AD<double> > y(m);
	y[0] = 3. * 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] , 3.*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_TEST_VECTOR<double> dx(n);
	CPPAD_TEST_VECTOR<double> dy(m);
	dx[0] = 1.;
	dy    = f.Forward(1, dx);
	ok   &= NearEqual(dy[0], 2., 1e-10, 1e-10);
	ok   &= NearEqual(dy[1], 2., 1e-10, 1e-10);

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

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

	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;
}

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

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

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

	// some binary multiplication 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_TEST_VECTOR< AD<double> > y(m);
	y[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*(4.*3.*2.*1.)*x0,  1e-10 , 1e-10);

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

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

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

	return ok;
}

bool Sub(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(1);
    x[0]      = x0;

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

    AD<double> a = 2. * 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] - 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-4.+3.+2.-2.*x0+1., 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], -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 subtraction
    CppAD::VecAD<double> v(1);
    AD<double> zero(0);
    v[zero] = b;
    AD<double> result = 3. - v[zero];
    ok     &= (result == c);

    return ok;
}