abs_normal Getting Started: Example and Test

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}} }@)@ 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