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

	/*!
 	Link from user_atomic to forward sparse Jacobian 

	\copydetails atomic_base::rev_sparse_hes
 	*/
	virtual bool rev_sparse_hes(
		const vector<bool>&                     vx ,
		const vector<bool>&                     s  ,
		      vector<bool>&                     t  ,
		size_t                                  q  ,
		const vector< std::set<size_t> >&       r  ,
		const vector< std::set<size_t> >&       u  ,
		      vector< std::set<size_t> >&       v  )
	{	size_t n       = v.size();
		size_t m       = u.size();
		CPPAD_ASSERT_UNKNOWN( r.size() == v.size() );
		CPPAD_ASSERT_UNKNOWN( s.size() == m );
		CPPAD_ASSERT_UNKNOWN( t.size() == n );
		bool ok        = true;
		bool transpose = true;
		std::set<size_t>::const_iterator itr;

		// compute sparsity pattern for T(x) = S(x) * f'(x)
		t = f_.RevSparseJac(1, s);
# ifndef NDEBUG
		for(size_t j = 0; j < n; j++)
			CPPAD_ASSERT_UNKNOWN( vx[j] || ! t[j] )
# endif

		// V(x) = f'(x)^T * g''(y) * f'(x) * R  +  g'(y) * f''(x) * R 
		// U(x) = g''(y) * f'(x) * R
		// S(x) = g'(y)
		
		// compute sparsity pattern for A(x) = f'(x)^T * U(x)
		vector< std::set<size_t> > a(n);
		a = f_.RevSparseJac(q, u, transpose);

		// set version of s
		vector< std::set<size_t> > set_s(1);
		CPPAD_ASSERT_UNKNOWN( set_s[0].empty() );
		size_t i;
		for(i = 0; i < m; i++)
			if( s[i] )
				set_s[0].insert(i);

		// compute sparsity pattern for H(x) = (S(x) * F)''(x) * R
		// (store it in v)
		f_.ForSparseJac(q, r);
		v = f_.RevSparseHes(q, set_s, transpose);

		// compute sparsity pattern for V(x) = A(x) + H(x)
		for(i = 0; i < n; i++)
		{	for(itr = a[i].begin(); itr != a[i].end(); itr++)
			{	size_t j = *itr;
				CPPAD_ASSERT_UNKNOWN( j < q );
				v[i].insert(j);
			}
		}

		// no longer need the forward mode sparsity pattern
		// (have to reconstruct them every time)
		f_.size_forward_set(0);

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

	/*!
 	Link from user_atomic to forward sparse Jacobian 

	\copydetails atomic_base::rev_sparse_jac
 	*/
	virtual bool rev_sparse_jac(
		size_t                                  q  ,
		const vector<bool>&                     rt ,
		      vector<bool>&                     st )
	{
		bool ok  = true;

		// compute rt
		bool transpose  = true;
		bool nz_compare = true;
		// 2DO: remove need for nz_compare all the time. It is only really
		// necessary when optimizer calls this member function.
		st = f_.RevSparseJac(q, rt, transpose, nz_compare);

		return ok; 
	}

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

//.........这里部分代码省略.........
	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);
	// --------------------------------------------------------------------
	// forward mode sparsity pattern for the Jacobian
	// f_u = [   1, 0,   1 ]
	//       [ u_2, 1, u_2 ]
	size_t i, j, p = n;
	CppAD::vectorBool r(n * p), s(m * p);
	// r = identity sparsity pattern
	for(i = 0; i < n; i++)
		for(j = 0; j < p; j++)
			r[i*n +j] = (i == j);
	s   = f.ForSparseJac(p, r);
	ok &= s[ 0 * p + 0] == true;
	ok &= s[ 0 * p + 1] == false;
	ok &= s[ 0 * p + 2] == true;
	ok &= s[ 1 * p + 0] == true;
	ok &= s[ 1 * p + 1] == true;
	ok &= s[ 1 * p + 2] == true;
	// --------------------------------------------------------------------
	// reverse mode sparsity pattern for the Jacobian
	q = m;
	s.resize(q * m);
	r.resize(q * n);
	// s = identity sparsity pattern
	for(i = 0; i < q; i++)
		for(j = 0; j < m; j++)
			s[i*m +j] = (i == j);
	r   = f.RevSparseJac(q, s);
	ok &= r[ 0 * n + 0] == true;
	ok &= r[ 0 * n + 1] == false;
	ok &= r[ 0 * n + 2] == true;
	ok &= r[ 1 * n + 0] == true;
	ok &= r[ 1 * n + 1] == true;
	ok &= r[ 1 * n + 2] == true;

	// --------------------------------------------------------------------
	// Hessian sparsity for y_1 (u) = u_1 + u_0 * u_2 + u_2^2 / 2
	s.resize(m);
	s[0] = false;
	s[1] = true;
	r.resize(n * n);
	for(i = 0; i < n; i++)
		for(j = 0; j < n; j++)
			r[ i * n + j ] = (i == j);
	CppAD::vectorBool h(n * n);
	h   = f.RevSparseHes(n, s);
	ok &= h[0 * n + 0] == false;
	ok &= h[0 * n + 1] == false;
	ok &= h[0 * n + 2] == true;
	ok &= h[1 * n + 0] == false;
	ok &= h[1 * n + 1] == false;
	ok &= h[1 * n + 2] == false;
	ok &= h[2 * n + 0] == true;
	ok &= h[2 * n + 1] == false;
	ok &= h[2 * n + 2] == true;

	// --------------------------------------------------------------------
	destroy_r();

	// Free all temporary work space associated with old_atomic objects.
	// (If there are future calls to user atomic functions, they will
	// create new temporary work space.)
	CppAD::user_atomic<double>::clear();

	return ok;
}