\(\newcommand{\W}[1]{ \; #1 \; }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} }\) \(\newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} }\) \(\newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} }\) \(\newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }\)
atomic_four_lin_ode_rev_depend.cpp¶
View page sourceAtomic Linear ODE Reverse Dependency Analysis: Example and Test¶
Purpose¶
This example demonstrates calculating reverse dependency with the atomic_four_lin_ode class; see atomic_four_lin_ode_rev_depend.hpp .
y(t, x)¶
We are given a positive integer m and define \(y : \B{R} \times \B{R}^m \rightarrow \B{R}^m\) by
\[y(t, x) = \prod_{i=0}^m x_i t^i / i !\]
It follows that
\[\begin{split}\partial_t y_i (t, x) = \left \{ \begin{array}{ll}
0 & \R{if} \; i = 0 \\
x_i y_{i-1} (t, x) & \R{otherwise}
\end{array} \right .
\\
\partial_t y (t, x) =
\left( \begin{array}{cccc}
0 & 0 & \cdots & 0 \\
x_1 & 0 & \cdots & 0 \\
0 & x_2 & & 0 \\
\vdots & & \ddots & \vdots \\
0 & 0 & \cdots & x_m
\end{array} \right)
y (t, x)
\W{,}
y (0, x) =
\left( \begin{array}{c}
x_0 \\
0 \\
\vdots \\
0 \\
\end{array} \right)\end{split}\]
Problem Parameters¶
The following problem parameters can be changed:
// m, r, step
size_t m = 5; // number of components in x and y
double r = 2.0; // final time in the ODE
double step = 1.0; // step size used to approximation ODE solution
Source¶
# include <cppad/cppad.hpp>
# include <cppad/example/atomic_four/lin_ode/lin_ode.hpp>
namespace { // BEGIN_EMPTY_NAMESPACE
// y(t, x)
template <class Scalar, class Vector>
Vector Y(Scalar t, const Vector& x)
{ size_t m = x.size();
Vector y(m);
//
// y
y[0] = x[0];
for(size_t i = 1; i < m; ++i)
y[i] = x[i] * y[i-1] * t / Scalar(i);
return y;
}
} // END_EMPTY_NAMESPACE
bool rev_depend(void)
{ // ok, eps
bool ok = true;
//
// sparse_rc, AD, eps99
typedef CppAD::sparse_rc< CppAD::vector<size_t> > sparse_rc;
using CppAD::AD;
double eps99 = std::numeric_limits<double>::epsilon() * 99.0;
// -----------------------------------------------------------------------
// Record f
// -----------------------------------------------------------------------
//
// afun
CppAD::atomic_lin_ode<double> afun("atomic_lin_ode");
//
// BEGIN_PROBLEM_PARAMETERS
// m, r, step
size_t m = 5; // number of components in x and y
double r = 2.0; // final time in the ODE
double step = 1.0; // step size used to approximation ODE solution
// END_PROBLEM_PARAMETERS
//
// pattern, transpose
size_t nr = m;
size_t nc = m;
size_t nnz = m - 1;
sparse_rc pattern(nr, nc, nnz);
for(size_t k = 0; k < nnz; ++k)
{ size_t i = k + 1;
size_t j = k;
pattern.set(k, i, j);
}
bool transpose = false;
//
// ax
CPPAD_TESTVECTOR( AD<double> ) ax(m);
for(size_t k = 0; k < m; ++k)
ax[k] = double(k + 1);
CppAD::Independent(ax);
//
// au
// au = (x[1], ..., x[nnz-1], x[0], 0, ..., 0)
CPPAD_TESTVECTOR( AD<double> ) au(nnz + m);
for(size_t k = 0; k < nnz; ++k)
au[k] = ax[k+1];
for(size_t i = 0; i < m; ++i)
{ if( i == 0 )
au[nnz + i] = ax[0];
else
au[nnz + i] = 0.0;
}
//
// ay
CPPAD_TESTVECTOR( AD<double> ) ay(m);
size_t call_id = afun.set(r, step, pattern, transpose);
afun(call_id, au, ay);
//
// z_index
// Fourth order Rosen34 method is exact approximation of y[i] for i <= 4
size_t z_index = m - 1;
assert(z_index <= 4);
//
// az
CPPAD_TESTVECTOR( AD<double> ) az(1);
az[0] = ay[z_index];
//
// f
// optimize uses rev_depend
CppAD::ADFun<double> f(ax, az);
f.optimize("val_graph no_conditional_skip");
// -----------------------------------------------------------------------
// check_f
// -----------------------------------------------------------------------
CppAD::Independent(ax);
AD<double> ar = r;
ay = Y(ar, ax);
az[0] = ay[z_index];
CppAD::ADFun<double> check_f(ax, az);
// -----------------------------------------------------------------------
// rev_depend
// use test_rev_depend to call rev_depend directly
// -----------------------------------------------------------------------
//
// depend_u
CppAD::vector<bool> ident_zero_u(nnz + m), depend_u(nnz + m), depend_y(m);
for(size_t i = 0; i < m; ++i)
{ depend_y[i] = i == z_index;
ident_zero_u[i] = false;
}
for(size_t i = 1; i < m; ++i)
ident_zero_u[nnz + i] = true;
afun.test_rev_depend(call_id, ident_zero_u, depend_u, depend_y);
//
// depend_x
CppAD::vector<bool> depend_x(m);
depend_x[0] = depend_u[m-1];
for(size_t j = 1; j < m; ++j)
depend_x[j] = depend_u[j-1];
//
// x
CPPAD_TESTVECTOR(double) x(m);
for(size_t j = 0; j < m; ++j)
x[j] = double( j + 2 );
//
// dw
check_f.Forward(0, x);
CPPAD_TESTVECTOR(double) w(1), dw(m);
w[0] = 1.0;
dw = check_f.Reverse(1, w);
//
// ok
// note that for this x, partial w.r.t x[j] is non-zero if and only if
// y[z_index] depends on x[j]
for(size_t j = 0; j < m; ++j)
ok &= depend_x[j] == (dw[j] != 0.0);
//
// -----------------------------------------------------------------------
// forward mode on f
// Check that the optimized version of agrees with check_f.
// -----------------------------------------------------------------------
//
// z
// zero order forward mode computation of f(x)
CPPAD_TESTVECTOR(double) z = f.Forward(0, x);
//
// ok
CPPAD_TESTVECTOR(double) check_z = check_f.Forward(0, x);
ok &= CppAD::NearEqual(z[0], check_z[0], eps99, eps99);
//
// du, ok
CPPAD_TESTVECTOR(double) dx(m), dz(1), check_dz(1);
for(size_t j = 0; j < m; ++j)
dx[j] = 0.0;
//
for(size_t j = 0; j < m; ++j)
{ dx[j] = 1.0;
dz = f.Forward(1, dx);
check_dz = check_f.Forward(1, dx);
ok &= CppAD::NearEqual(dz[0], check_dz[0], eps99, eps99);
dx[j] = 0.0;
}
// -----------------------------------------------------------------------
return ok;
}