Prev Next abs_get_started.cpp

@(@\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}} }@)@This is cppad-20221105 documentation. Here is a link to its current documentation .
abs_normal Getting Started: Example and Test

Purpose
Creates an abs_normal representation @(@ g @)@ for the function @(@ f : \B{R}^3 \rightarrow \B{R} @)@ defined by @[@ f( x_0, x_1, x_2 ) = | x_0 + x_1 | + | x_1 + x_2 | @]@ The corresponding g @(@ : \B{R}^5 \rightarrow \B{R}^3 @)@ is given by @[@ \begin{array}{rclrcl} g_0 ( x_0, x_1, x_2, u_0, u_1 ) & = & u_0 + u_1 & = & y_0 (x, u) \\ g_1 ( x_0, x_1, x_2, u_0, u_1 ) & = & x_0 + x_1 & = & z_0 (x, u) \\ g_1 ( x_0, x_1, x_2, u_0, u_1 ) & = & x_1 + x_2 & = & z_1 (x, u) \end{array} @]@

Source 

# include <cppad/cppad.hpp>
namespace {
    CPPAD_TESTVECTOR(double) join(
        const CPPAD_TESTVECTOR(double)& x ,
        const CPPAD_TESTVECTOR(double)& u )
    {   size_t n = x.size();
        size_t s = u.size();
        CPPAD_TESTVECTOR(double) xu(n + s);
        for(size_t j = 0; j < n; j++)
            xu[j] = x[j];
        for(size_t j = 0; j < s; j++)
            xu[n + j] = u[j];
        return xu;
    }
}
bool get_started(void)
{   bool ok = true;
    //
    using CppAD::AD;
    using CppAD::ADFun;
    //
    size_t n = 3; // size of x
    size_t m = 1; // size of y
    size_t s = 2; // size of u and z
    //
    // record the function f(x)
    CPPAD_TESTVECTOR( AD<double> ) ax(n), ay(m);
    for(size_t j = 0; j < n; j++)
        ax[j] = double(j + 1);
    Independent( ax );
    // for this example, we ensure first absolute value is | x_0 + x_1 |
    AD<double> a0 = abs( ax[0] + ax[1] );
    // and second absolute value is | x_1 + x_2 |
    AD<double> a1 = abs( ax[1] + ax[2] );
    ay[0]         = a0 + a1;
    ADFun<double> f(ax, ay);

    // create its abs_normal representation in g, a
    ADFun<double> g, a;
    f.abs_normal_fun(g, a);

    // check dimension of domain and range space for g
    ok &= g.Domain() == n + s;
    ok &= g.Range() == m + s;

    // check dimension of domain and range space for a
    ok &= a.Domain() == n;
    ok &= a.Range() == s;

    // --------------------------------------------------------------------
    // a(x) has all the operations used to compute f(x), but the sum of the
    // absolute values is not needed for a(x), so optimize it out.
    size_t n_op = f.size_op();
    ok         &= a.size_op() == n_op;
    a.optimize();
    ok         &= a.size_op() < n_op;

    // --------------------------------------------------------------------
    // zero order forward mode calculation using g(x, u)
    CPPAD_TESTVECTOR(double) x(n), u(s), xu(n+s), yz(m+s);
    for(size_t j = 0; j < n; j++)
        x[j] = double(j + 2);
    for(size_t j = 0; j < s; j++)
        u[j] = double(j + n + 2);
    xu = join(x, u);
    yz = g.Forward(0, xu);

    // check y_0(x, u)
    double y0 = u[0] + u[1];
    ok       &= y0 == yz[0];

    // check z_0 (x, u)
    double z0 = x[0] + x[1];
    ok       &= z0 == yz[1];

    // check z_1 (x, u)
    double z1 = x[1] + x[2];
    ok       &= z1 == yz[2];


    // --------------------------------------------------------------------
    // check that y(x, a(x) ) == f(x)
    CPPAD_TESTVECTOR(double) y(m);
    y  = f.Forward(0, x);  // y  = f(x)
    u  = a.Forward(0, x);  // u  = a(x)
    xu = join(x, u);       // xu = ( x, a(x) )
    yz = g.Forward(0, xu); // yz = ( y(x, a(x)), z(x, a(x)) )
    ok &= yz[0] == y[0];

    return ok;
}
Input File: example/abs_normal/get_started.cpp