abs_get_started.cpp¶

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{split}\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}\end{split}\]

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