|
|
Prev | Next | a_fun_reverse_xam.m | Headings |
function ok = a_fun_reverse_xam()
%
% load the Cppad Swig library
m_cppad
%
% initialize return variable
ok = true;
% -----------------------------------------------------------------------
% number of dependent and independent variables
n_dep = 1;
n_ind = 3;
%
% create the independent variables ax
xp = m_cppad.vec_double(n_ind);
for i = [ 0 :(n_ind -1) ]
xp(i) = i;
end
ax = m_cppad.independent(xp);
%
% create dependent variables ay with ay0 = ax_0 * ax_1 * ax_2
ax_0 = ax(0);
ax_1 = ax(1);
ax_2 = ax(2);
ay = m_cppad.vec_a_double(n_dep);
ay(0) = ax_0 * ax_1 * ax_2;
%
% define af corresponding to f(x) = x_0 * x_1 * x_2
af = m_cppad.a_fun(ax, ay);
% -----------------------------------------------------------------------
% define X(t) = (x_0 + t, x_1 + t, x_2 + t)
% it follows that Y(t) = f(X(t)) = (x_0 + t) * (x_1 + t) * (x_2 + t)
% and that Y'(0) = x_1 * x_2 + x_0 * x_2 + x_0 * x_1
% -----------------------------------------------------------------------
% zero order forward mode
p = 0;
xp(0) = 2.0;
xp(1) = 3.0;
xp(2) = 4.0;
yp = af.forward(p, xp);
ok = ok && yp(0) == 24.0;
% -----------------------------------------------------------------------
% first order reverse (derivative of zero order forward)
% define G( Y ) = y_0 = x_0 * x_1 * x_2
q = 1;
yq1 = m_cppad.vec_double(n_dep);
yq1(0) = 1.0;
xq1 = af.reverse(q, yq1);
% partial G w.r.t x_0
ok = ok && xq1(0) == 3.0 * 4.0 ;
% partial G w.r.t x_1
ok = ok && xq1(1) == 2.0 * 4.0 ;
% partial G w.r.t x_2
ok = ok && xq1(2) == 2.0 * 3.0 ;
% -----------------------------------------------------------------------
% first order forward mode
p = 1;
xp(0) = 1.0;
xp(1) = 1.0;
xp(2) = 1.0;
yp = af.forward(p, xp);
ok = ok && yp(0) == 3.0*4.0 + 2.0*4.0 + 2.0*3.0;
% -----------------------------------------------------------------------
% second order reverse (derivative of first order forward)
% define G( y_0^0 , y_0^1 ) = y_0^1
% = x_1^0 * x_2^0 + x_0^0 * x_2^0 + x_0^0 * x_1^0
q = 2;
yq2 = m_cppad.vec_double(n_dep * q);
yq2(0 * q + 0) = 0.0; % partial of G w.r.t y_0^0
yq2(0 * q + 1) = 1.0; % partial of G w.r.t y_0^1
xq2 = af.reverse(q, yq2);
% partial G w.r.t x_0^0
ok = ok && xq2(0 * q + 0) == 3.0 + 4.0;
% partial G w.r.t x_1^0
ok = ok && xq2(1 * q + 0) == 2.0 + 4.0;
% partial G w.r.t x_2^0
ok = ok && xq2(2 * q + 0) == 2.0 + 3.0;
% -----------------------------------------------------------------------
%
return;
end