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

	/*!
 	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;
	}

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