/* $Id: change_const.cpp 2506 2012-10-24 19:36:49Z bradbell $ */
/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
CppAD is distributed under multiple licenses. This distribution is under
the terms of the
GNU General Public License Version 3.
A copy of this license is included in the COPYING file of this distribution.
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
-------------------------------------------------------------------------- */
# include <cppad/cppad.hpp>
/*
$begin change_const.cpp$$
$spell
Jacobian
$$
$section Computing a Jacobian With Constants that Change$$
$index multiple, AD level$$
$index level, multiple AD$$
$index constant, that change$$
$index change, constant$$
$head Purpose$$
In this example we use two levels of taping so that a derivative
can have constant parameters that can be changed. To be specific,
we consider the function $latex f : \B{R}^2 \rightarrow \B{R}^2$$
$latex \[
f(x) = p \left( \begin{array}{c}
\sin( x_0 ) \\
\sin( x_1 )
\end{array} \right)
\]$$
were $latex p \in \B{R}$$ is a parameter.
The Jacobian of this function is
$latex \[
g(x,p) = p \left( \begin{array}{cc}
\cos( x_0 ) & 0 \\
0 & \cos( x_1 )
\end{array} \right)
\] $$
In this example we use two levels of AD to avoid computing
the partial of $latex f(x)$$ with respect to $latex p$$,
but still allow for the evaluation of $latex g(x, p)$$
at different values of $latex p$$.
$end
*/
bool change_const(void)
{ bool ok = true; // initialize test result
typedef CppAD::AD<double> A1_double; // for first level of taping
typedef CppAD::AD<A1_double> A2_double; // for second level of taping
size_t nu = 3; // number components in u
size_t nx = 2; // number components in x
size_t ny = 2; // num components in f(x)
size_t nJ = ny * nx; // number components in Jacobian of f(x)
// temporary indices
size_t j;
// declare first level of independent variables
CPPAD_TESTVECTOR(A1_double) a1_u(nu);
for(j = 0; j < nu; j++)
a1_u[j] = 0.;
CppAD::Independent(a1_u);
// parameter in computation of Jacobian
A1_double a1_p = a1_u[2];
// declare second level of independent variables
CPPAD_TESTVECTOR(A2_double) a2_x(nx);
for(j = 0; j < nx; j++)
a2_x[j] = 0.;
CppAD::Independent(a2_x);
// compute dependent variables at second level
CPPAD_TESTVECTOR(A2_double) a2_y(ny);
a2_y[0] = sin( a2_x[0] ) * a1_p;
a2_y[1] = sin( a2_x[1] ) * a1_p;
// declare function object that computes values at the first level
// (make sure we do not run zero order forward during constructor)
CppAD::ADFun<A1_double> a1_f;
a1_f.Dependent(a2_x, a2_y);
// compute the Jacobian of a1_f at a1_u[0], a1_u[1]
CPPAD_TESTVECTOR(A1_double) a1_x(nx);
a1_x[0] = a1_u[0];
a1_x[1] = a1_u[1];
CPPAD_TESTVECTOR(A1_double) a1_J(nJ);
a1_J = a1_f.Jacobian( a1_x );
// declare function object that maps u = (x, p) to Jacobian of f
// (make sure we do not run zero order forward during constructor)
CppAD::ADFun<double> g;
g.Dependent(a1_u, a1_J);
// remove extra variables used during the reconding of a1_f,
// but not needed any more.
g.optimize();
// compute the Jacobian of f using zero order forward
// sweep with double values
CPPAD_TESTVECTOR(double) J(nJ), u(nu);
for(j = 0; j < nu; j++)
u[j] = double(j+1);
J = g.Forward(0, u);
// accuracy for tests
double eps = 100. * CppAD::numeric_limits<double>::epsilon();
// y[0] = sin( x[0] ) * p
// y[1] = sin( x[1] ) * p
CPPAD_TESTVECTOR(double) x(nx);
x[0] = u[0];
x[1] = u[1];
double p = u[2];
// J[0] = partial y[0] w.r.t x[0] = cos( x[0] ) * p
double check = cos( x[0] ) * p;
ok &= fabs( check - J[0] ) <= eps;
// J[1] = partial y[0] w.r.t x[1] = 0.;
check = 0.;
ok &= fabs( check - J[1] ) <= eps;
// J[2] = partial y[1] w.r.t. x[0] = 0.
check = 0.;
ok &= fabs( check - J[2] ) <= eps;
// J[3] = partial y[1] w.r.t x[1] = cos( x[1] ) * p
check = cos( x[1] ) * p;
ok &= fabs( check - J[3] ) <= eps;
return ok;
}
// END PROGRAM