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

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

    void testDynamicFull(std::vector<ADCG>& u,
                         const std::vector<double>& x,
                         const std::vector<double>& xNorm,
                         const std::vector<double>& eqNorm,
                         size_t maxAssignPerFunc = 100,
                         double epsilonR = 1e-14,
                         double epsilonA = 1e-14) {
        ASSERT_EQ(u.size(), x.size());
        ASSERT_EQ(x.size(), xNorm.size());

        using namespace std;

        // use a special object for source code generation
        CppAD::Independent(u);

        for (size_t i = 0; i < u.size(); i++)
            u[i] *= xNorm[i];

        // dependent variable vector 
        std::vector<ADCG> Z = model(u);

        if (eqNorm.size() > 0) {
            ASSERT_EQ(Z.size(), eqNorm.size());
            for (size_t i = 0; i < Z.size(); i++)
                Z[i] /= eqNorm[i];
        }

        /**
         * create the CppAD tape as usual
         */
        // create f: U -> Z and vectors used for derivative calculations
        ADFun<CGD> fun;
        fun.Dependent(Z);

        /**
         * Create the dynamic library
         * (generate and compile source code)
         */
        ModelCSourceGen<double> compHelp(fun, _name + "dynamic");

        compHelp.setCreateForwardZero(true);
        compHelp.setCreateJacobian(true);
        compHelp.setCreateHessian(true);
        compHelp.setCreateSparseJacobian(true);
        compHelp.setCreateSparseHessian(true);
        compHelp.setCreateForwardOne(true);
        compHelp.setCreateReverseOne(true);
        compHelp.setCreateReverseTwo(true);
        compHelp.setMaxAssignmentsPerFunc(maxAssignPerFunc);

        ModelLibraryCSourceGen<double> compDynHelp(compHelp);

        SaveFilesModelLibraryProcessor<double>::saveLibrarySourcesTo(compDynHelp, "sources_" + _name + "_1");

        DynamicModelLibraryProcessor<double> p(compDynHelp);
        GccCompiler<double> compiler;
        DynamicLib<double>* dynamicLib = p.createDynamicLibrary(compiler);

        /**
         * test the library
         */
        GenericModel<double>* model = dynamicLib->model(_name + "dynamic");
        ASSERT_TRUE(model != nullptr);

        testModelResults(*model, fun, x, epsilonR, epsilonA);

        delete model;
        delete dynamicLib;
    }

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);
	// --------------------------------------------------------------------
//.........这里部分代码省略.........

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]
//.........这里部分代码省略.........