atomic_four_lin_ode_rev_depend.cpp

View page source

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