C++ ADFun::Forward方法代码示例

	/*!
 	Link from user_atomic to reverse mode 

	\copydetails atomic_base::reverse
 	*/
	virtual bool reverse(
		size_t                    q  ,
		const vector<Base>&       tx ,
		const vector<Base>&       ty ,
		      vector<Base>&       px ,
		const vector<Base>&       py )
	{
		CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 );
		CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
		CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
		bool ok  = true;	

		// put proper forward mode coefficients in f_
# ifdef NDEBUG
		f_.Forward(q, tx);
# else
		size_t n = tx.size() / (q+1);
		size_t m = ty.size() / (q+1);
		CPPAD_ASSERT_UNKNOWN( px.size() == n * (q+1) );
		CPPAD_ASSERT_UNKNOWN( py.size() == m * (q+1) );
		size_t i, j, k;
		//
		vector<Base> check_ty = f_.Forward(q, tx);
		for(i = 0; i < m; i++)
		{	for(k = 0; k <= q; k++)
			{	j = i * (q+1) + k;
				CPPAD_ASSERT_UNKNOWN( check_ty[j] == ty[j] );
			}
		}
# endif
		// now can run reverse mode
		px = f_.Reverse(q+1, py);

		// no longer need the Taylor coefficients in f_
		// (have to reconstruct them every time)
		size_t c = 0;
		size_t r = 0;
		f_.capacity_order(c, r);
		return ok;
	}

	/*!
 	Link from user_atomic to forward mode 

	\copydetails atomic_base::forward
 	*/
	virtual bool forward(
		size_t                    p ,
		size_t                    q ,
		const vector<bool>&      vx , 
		      vector<bool>&      vy , 
		const vector<Base>&      tx ,
		      vector<Base>&      ty )
	{
		CPPAD_ASSERT_UNKNOWN( f_.size_var() > 0 );
		CPPAD_ASSERT_UNKNOWN( tx.size() % (q+1) == 0 );
		CPPAD_ASSERT_UNKNOWN( ty.size() % (q+1) == 0 );
		size_t n = tx.size() / (q+1);
		size_t m = ty.size() / (q+1);
		bool ok  = true;	
		size_t i, j;

		// 2DO: test both forward and reverse vy information
		if( vx.size() > 0 )
		{	//Compute Jacobian sparsity pattern.
			vector< std::set<size_t> > s(m);
			if( n <= m )
			{	vector< std::set<size_t> > r(n);
				for(j = 0; j < n; j++)
					r[j].insert(j);
				s = f_.ForSparseJac(n, r);
			}
			else
			{	vector< std::set<size_t> > r(m);
				for(i = 0; i < m; i++)
					r[i].insert(i);
				s = f_.RevSparseJac(m, r);
			}
			std::set<size_t>::const_iterator itr;
			for(i = 0; i < m; i++)
			{	vy[i] = false;
				for(itr = s[i].begin(); itr != s[i].end(); itr++)
				{	j = *itr;
					assert( j < n );
					// y[i] depends on the value of x[j]
					vy[i] |= vx[j];
				}
			}
		}
		ty = f_.Forward(q, tx);

		// no longer need the Taylor coefficients in f_
		// (have to reconstruct them every time)
		size_t c = 0;
		size_t r = 0;
		f_.capacity_order(c, r);
		return ok;
	}

	bool FunCheck(
		ADFun<Base>  &f , 
		Fun          &g , 
		const Vector &x , 
		const Base   &r ,
		const Base   &a )
	{	bool ok = true;
	
		size_t m   = f.Range();
		Vector yf  = f.Forward(0, x); 
		Vector yg  = g(x);
	
		size_t i;
		for(i = 0; i < m; i++)
			ok  &= NearEqual(yf[i], yg[i], r, a);
		return ok;
	}

    void testModel(ADFun<CGD>& f,
                   size_t expectedTmp,
                   size_t expectedArraySize) {
        using CppAD::vector;

        size_t n = f.Domain();
        //size_t m = f.Range();

        CodeHandler<double> handler(10 + n * n);

        vector<CGD> indVars(n);
        handler.makeVariables(indVars);

        vector<CGD> dep = f.Forward(0, indVars);

        LanguageC<double> langC("double");
        LangCDefaultVariableNameGenerator<double> nameGen;

        handler.generateCode(std::cout, langC, dep, nameGen);

        ASSERT_EQ(handler.getTemporaryVariableCount(), expectedTmp);
        ASSERT_EQ(handler.getTemporaryArraySize(), expectedArraySize);
    }

void JacobianFor(ADFun<Base> &f, const Vector &x, Vector &jac)
{	size_t i;
	size_t j;

	size_t n = f.Domain();
	size_t m = f.Range();

	// check Vector is Simple Vector class with Base type elements
	CheckSimpleVector<Base, Vector>();

	CPPAD_ASSERT_UNKNOWN( size_t(x.size())   == f.Domain() );
	CPPAD_ASSERT_UNKNOWN( size_t(jac.size()) == f.Range() * f.Domain() );

	// argument and result for forward mode calculations
	Vector u(n);
	Vector v(m);

	// initialize all the components
	for(j = 0; j < n; j++)
		u[j] = Base(0);

	// loop through the different coordinate directions
	for(j = 0; j < n; j++)
	{	// set u to the j-th coordinate direction
		u[j] = Base(1);

		// compute the partial of f w.r.t. this coordinate direction
		v = f.Forward(1, u);

		// reset u to vector of all zeros
		u[j] = Base(0);

		// return the result
		for(i = 0; i < m; i++)
			jac[ i * n + j ] = v[i];
	}
}

TEST_F(CppADCGEvaluatorTest, Atomic) {
    using ADCG = AD<CGD>;

    std::vector<AD<double>> ax(2);
    std::vector<AD<double>> ay(3);
    for (size_t j = 0; j < ax.size(); j++) {
        ax[j] = j + 2;
    }

    checkpoint<double> atomicFun("func", testModel, ax, ay); // the normal atomic function
    CGAtomicFun<double> atomic(atomicFun, ax); // a wrapper used to tape with CG<Base>

    ModelType model = [&](const std::vector<CGD>& x) {
        // independent variables
        std::vector<ADCG> ax(x.size());
        for (size_t j = 0; j < ax.size(); j++) {
            ax[j] = x[j];
        }

        CppAD::Independent(ax);

        // dependent variable vector
        std::vector<ADCG> ay(3);

        atomic(ax, ay);

        ADFun<CGD> fun;
        fun.Dependent(ay);

        std::vector<CGD> y = fun.Forward(0, x);

        return y;
    };

    this->testCG(model, std::vector<double>{0.5, 1.5});
}

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

	// --------------------------------------------------------------------
	// Create the ADFun<doulbe> r_
	create_r();

	// --------------------------------------------------------------------
	// 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 reciprocal(x) in au[0]
	vector< AD<double> > au(m);
	size_t id = 0;           // not used
	reciprocal(id, ax, au);	// u = 1 / x

	// call user function and store reciprocal(u) in ay[0]
	reciprocal(id, au, ay);	// y = 1 / u = x

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

	// --------------------------------------------------------------------
	// Check function value results
	//
	// check function value
	double check = x0;
	ok &= NearEqual( Value(ay[0]) , check,  eps, eps);

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

	// check first order forward mode
	q      = 1;
	x_q[0] = 1;
	y_q    = f.Forward(q, x_q);
	check  = 1.;
	ok &= NearEqual(y_q[0] , check,  eps, eps);

	// check second order forward mode
	q      = 2;
	x_q[0] = 0;
	y_q    = f.Forward(q, x_q);
	check  = 0.;
	ok &= NearEqual(y_q[0] , check,  eps, eps);

	// --------------------------------------------------------------------
	// Check reverse mode results
	//
	// third order reverse mode
	q     = 3;
	vector<double> w(m), dw(n * q);
	w[0]  = 1.;
	dw    = f.Reverse(q, w);
	check = 1.;
	ok &= NearEqual(dw[0] , check,  eps, eps);
	check = 0.;
	ok &= NearEqual(dw[1] , check,  eps, eps);
	ok &= NearEqual(dw[2] , check,  eps, eps);

	// --------------------------------------------------------------------
	// forward mode sparstiy pattern
	size_t p = n;
	CppAD::vectorBool r1(n * p), s1(m * p);
	r1[0] = true;          // compute sparsity pattern for x[0]
	s1    = f.ForSparseJac(p, r1);
	ok  &= s1[0] == true;  // f[0] depends on x[0]

	// --------------------------------------------------------------------
	// reverse mode sparstiy pattern
	q = m;
	CppAD::vectorBool s2(q * m), r2(q * n);
	s2[0] = true;          // compute sparsity pattern for f[0]
	r2    = f.RevSparseJac(q, s2);
	ok  &= r2[0] == true;  // f[0] depends on x[0]

	// --------------------------------------------------------------------
	// Hessian sparsity (using previous ForSparseJac call)
	CppAD::vectorBool s3(m), h(p * n);
	s3[0] = true;        // compute sparsity pattern for f[0]
//.........这里部分代码省略.........

void BenderQuad(
	const BAvector   &x     , 
	const BAvector   &y     , 
	Fun               fun   , 
	BAvector         &g     ,
	BAvector         &gx    ,
	BAvector         &gxx   )
{	// determine the base type
	typedef typename BAvector::value_type Base;

	// check that BAvector is a SimpleVector class
	CheckSimpleVector<Base, BAvector>();

	// declare the ADvector type
	typedef CPPAD_TESTVECTOR(AD<Base>) ADvector;

	// size of the x and y spaces
	size_t n = size_t(x.size());
	size_t m = size_t(y.size());

	// check the size of gx and gxx
	CPPAD_ASSERT_KNOWN(
		g.size() == 1,
		"BenderQuad: size of the vector g is not equal to 1"
	);
	CPPAD_ASSERT_KNOWN(
		size_t(gx.size()) == n,
		"BenderQuad: size of the vector gx is not equal to n"
	);
	CPPAD_ASSERT_KNOWN(
		size_t(gxx.size()) == n * n,
		"BenderQuad: size of the vector gxx is not equal to n * n"
	);

	// some temporary indices
	size_t i, j;

	// variable versions x
	ADvector vx(n);
	for(j = 0; j < n; j++)
		vx[j] = x[j];
	
	// declare the independent variables
	Independent(vx);

	// evaluate h = H(x, y) 
	ADvector h(m);
	h = fun.h(vx, y);

	// evaluate dy (x) = Newton step as a function of x through h only
	ADvector dy(m);
	dy = fun.dy(x, y, h);

	// variable version of y
	ADvector vy(m);
	for(j = 0; j < m; j++)
		vy[j] = y[j] + dy[j];

	// evaluate G~ (x) = F [ x , y + dy(x) ] 
	ADvector gtilde(1);
	gtilde = fun.f(vx, vy);

	// AD function object that corresponds to G~ (x)
	// We will make heavy use of this tape, so optimize it
	ADFun<Base> Gtilde;
	Gtilde.Dependent(vx, gtilde); 
	Gtilde.optimize();

	// value of G(x)
	g = Gtilde.Forward(0, x);

	// initial forward direction vector as zero
	BAvector dx(n);
	for(j = 0; j < n; j++)
		dx[j] = Base(0);

	// weight, first and second order derivative values
	BAvector dg(1), w(1), ddw(2 * n);
	w[0] = 1.;


	// Jacobian and Hessian of G(x) is equal Jacobian and Hessian of Gtilde
	for(j = 0; j < n; j++)
	{	// compute partials in x[j] direction
		dx[j] = Base(1);
		dg    = Gtilde.Forward(1, dx);
		gx[j] = dg[0];

		// restore the dx vector to zero
		dx[j] = Base(0);

		// compute second partials w.r.t x[j] and x[l]  for l = 1, n
		ddw = Gtilde.Reverse(2, w);
		for(i = 0; i < n; i++)
			gxx[ i * n + j ] = ddw[ i * 2 + 1 ];
	}

	return;
}

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

	// --------------------------------------------------------------------
	// Create the ADFun<doulbe> r_
	create_r();

	// --------------------------------------------------------------------
	// domain and range space vectors
	size_t n = 3, m = 2;
	vector< AD<double> > au(n), ax(n), ay(m);
	au[0]         = 0.0;        // value of z_0 (t) = t, at t = 0
	ax[1]         = 0.0;        // value of z_1 (t) = t^2/2, at t = 0
	au[2]         = 1.0;        // final t
	CppAD::Independent(au);
	size_t M      = 2;          // number of r steps to take
	ax[0]         = au[0];      // value of z_0 (t) = t, at t = 0
	ax[1]         = au[1];      // value of z_1 (t) = t^2/2, at t = 0
	AD<double> dt = au[2] / double(M);  // size of each r step
	ax[2]         = dt;
	for(size_t i_step = 0; i_step < M; i_step++)
	{	size_t id = 0;               // not used
		solve_ode(id, ax, ay);
		ax[0] = ay[0];
		ax[1] = ay[1];
	}

	// create f: u -> y and stop tape recording
	// y_0(t) = u_0 + t                   = u_0 + u_2
	// y_1(t) = u_1 + u_0 * t + t^2 / 2   = u_1 + u_0 * u_2 + u_2^2 / 2
	// where t = u_2
	ADFun<double> f;
	f.Dependent(au, ay);

	// --------------------------------------------------------------------
	// Check forward mode results
	//
	// zero order forward
	vector<double> up(n), yp(m);
	size_t q  = 0;
	double u0 = 0.5;
	double u1 = 0.25;
	double u2 = 0.75;
	double check;
	up[0]     = u0;
	up[1]     = u1;
	up[2]     = u2;
	yp        = f.Forward(q, up);
	check     = u0 + u2;
	ok       &= NearEqual( yp[0], check,  eps, eps);
	check     = u1 + u0 * u2 + u2 * u2 / 2.0;
	ok       &= NearEqual( yp[1], check,  eps, eps);
	//
	// forward mode first derivative w.r.t t
	q         = 1;
	up[0]     = 0.0;
	up[1]     = 0.0;
	up[2]     = 1.0;
	yp        = f.Forward(q, up);
	check     = 1.0;
	ok       &= NearEqual( yp[0], check,  eps, eps);
	check     = u0 + u2;
	ok       &= NearEqual( yp[1], check,  eps, eps);
	//
	// forward mode second order Taylor coefficient w.r.t t
	q         = 2;
	up[0]     = 0.0;
	up[1]     = 0.0;
	up[2]     = 0.0;
	yp        = f.Forward(q, up);
	check     = 0.0;
	ok       &= NearEqual( yp[0], check,  eps, eps);
	check     = 1.0 / 2.0;
	ok       &= NearEqual( yp[1], check,  eps, eps);
	// --------------------------------------------------------------------
	// reverse mode derivatives of \partial_t y_1 (t)
	vector<double> w(m * q), dw(n * q);
	w[0 * q + 0]  = 0.0;
	w[1 * q + 0]  = 0.0;
	w[0 * q + 1]  = 0.0;
	w[1 * q + 1]  = 1.0;
	dw        = f.Reverse(q, w);
	// derivative of y_1(u) = u_1 + u_0 * u_2 + u_2^2 / 2,  w.r.t. u
	// is equal deritative of \partial_u2 y_1(u) w.r.t \partial_u2 u
	check     = u2;
	ok       &= NearEqual( dw[0 * q + 1], check,  eps, eps);
	check     = 1.0;
	ok       &= NearEqual( dw[1 * q + 1], check,  eps, eps);
	check     = u0 + u2;
	ok       &= NearEqual( dw[2 * q + 1], check,  eps, eps);
	// derivative of \partial_t y_1 w.r.t u = u_0 + t,  w.r.t u
	check     = 1.0;
	ok       &= NearEqual( dw[0 * q + 0], check,  eps, eps);
	check     = 0.0;
	ok       &= NearEqual( dw[1 * q + 0], check,  eps, eps);
	check     = 1.0;
	ok       &= NearEqual( dw[2 * q + 0], check,  eps, eps);
	// --------------------------------------------------------------------
//.........这里部分代码省略.........