atomic_four_forward¶

Atomic Function Forward Mode¶

Syntax¶

Base¶

ok = afun . forward (
      call_id , select_y ,
      order_low , order_up , type_x , taylor_x , taylor_y
)

AD<Base>¶

ok = afun . forward (
      call_id , select_y ,
      order_low , order_up , type_x , ataylor_x , ataylor_y
)

Prototype¶

Base¶

template <class Base>
bool atomic_four<Base>::forward(
   size_t                       call_id     ,
   const vector<bool>&          select_y    ,
   size_t                       order_low   ,
   size_t                       order_up    ,
   const vector<Base>&          taylor_x    ,
   vector<Base>&                taylor_y    )

AD<Base>¶

template <class Base>
bool atomic_four<Base>::forward(
   size_t                       call_id      ,
   const vector<bool>&          select_y    ,
   size_t                       order_low    ,
   size_t                       order_up     ,
   const vector< AD<Base> >&    ataylor_x    ,
   vector< AD<Base> >&          ataylor_y    )

Base¶

see Base .

vector¶

is the CppAD_vector template class.

Usage¶

Base¶

The Base syntax and prototype are used by a call to the atomic function afun . They are also used by f . Forward and f . new_dynamic where f has prototype

ADFun < Base > f

and afun is used during the recording of f .

AD<Base>¶

The AD < Base > syntax and prototype are used by af . Forward and af . new_dynamic where af has prototype

ADFun< AD< Base > , Base > af

and afun is used in a function af , created from f using base2ad .

Implementation¶

The taylor_x , taylor_y version of this function must be defined by the atomic_user class. It can return ok == false (and not compute anything) for values of order_up that are greater than those used by your Forward mode calculations. Order zero must be implemented.

call_id¶

See call_id .

select_y¶

This argument has size equal to the number of results to this atomic function; i.e. the size of ay . It specifies which components of y the corresponding Taylor coefficients must be computed.

order_low¶

This argument specifies the lowest order Taylor coefficient that we are computing.

p¶

We sometimes use the notation p = order_low below.

order_up¶

This argument is the highest order Taylor coefficient that we are computing ( order_low <= order_up ).

q¶

We use the notation q = order_up + 1 below. This is the number of Taylor coefficients for each component of x and y .

taylor_x¶

The size of taylor_x is q * n . For \(j = 0 , \ldots , n-1\) and \(k = 0 , \ldots , q-1\), we use the Taylor coefficient notation

\begin{eqnarray} x_j^k & = & \R{taylor\_x} [ j * q + k ] \\ X_j (t) & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^{q-1} t^{q-1} \end{eqnarray}

Note that superscripts represent an index for \(x_j^k\) and an exponent for \(t^k\). Also note that the Taylor coefficients for \(X(t)\) correspond to the derivatives of \(X(t)\) at \(t = 0\) in the following way:

\[x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)\]

parameters¶

If the j-th component of x is a parameter,

type_x [ j ] < CppAD::variable_enum

In this case, for k > 0 ,

taylor_x [ j * q + k ] == 0

ataylor_x¶

The specifications for ataylor_x is the same as for taylor_x (only the type of ataylor_x is different).

taylor_y¶

The size of taylor_y is q * m . Upon return, For \(i = 0 , \ldots , m-1\) and \(k = 0 , \ldots , q-1\), if select_y [ i ] is true,

\begin{eqnarray} Y_i (t) & = & g_i [ X(t) ] \\ Y_i (t) & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^{q-1} t^{q-1} + o( t^{q-1} ) \\ \R{taylor\_y} [ i * q + k ] & = & y_i^k \end{eqnarray}

where \(o( t^{q-1} ) / t^{q-1} \rightarrow 0\) as \(t \rightarrow 0\). Note that superscripts represent an index for \(y_j^k\) and an exponent for \(t^k\). Also note that the Taylor coefficients for \(Y(t)\) correspond to the derivatives of \(Y(t)\) at \(t = 0\) in the following way:

\[y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)\]

If \(p > 0\), for \(i = 0 , \ldots , m-1\) and \(k = 0 , \ldots , p-1\), the input of taylor_y satisfies

\[\R{taylor\_y} [ i * q + k ] = y_i^k\]

These values do not need to be recalculated and can be used during the computation of the higher order coefficients.

ataylor_y¶

The specifications for ataylor_y is the same as for taylor_y (only the type of ataylor_y is different).

ok¶

If this calculation succeeded, ok is true. Otherwise, it is false.

Discussion¶

For example, suppose that order_up == 2 , and you know how to compute the function \(g(x)\), its first derivative \(g^{(1)} (x)\), and it component wise Hessian \(g_i^{(2)} (x)\). Then you can compute taylor_x using the following formulas:

\begin{eqnarray} y_i^0 & = & Y(0) = g_i ( x^0 ) \\ y_i^1 & = & Y^{(1)} ( 0 ) = g_i^{(1)} ( x^0 ) X^{(1)} ( 0 ) = g_i^{(1)} ( x^0 ) x^1 \\ y_i^2 & = & \frac{1}{2 !} Y^{(2)} (0) \\ & = & \frac{1}{2} X^{(1)} (0)^\R{T} g_i^{(2)} ( x^0 ) X^{(1)} ( 0 ) + \frac{1}{2} g_i^{(1)} ( x^0 ) X^{(2)} ( 0 ) \\ & = & \frac{1}{2} (x^1)^\R{T} g_i^{(2)} ( x^0 ) x^1 + g_i^{(1)} ( x^0 ) x^2 \end{eqnarray}

For \(i = 0 , \ldots , m-1\), and \(k = 0 , 1 , 2\),

\[\R{taylor\_y} [ i * q + k ] = y_i^k\]

Example¶

The following is an example forward definition taken from atomic_four_norm_sq.cpp :

      bool forward(
         size_t                             call_id     ,
         const CppAD::vector<bool>&         select_y    ,
         size_t                             order_low   ,
         size_t                             order_up    ,
         const CppAD::vector<double>&       tx          ,
         CppAD::vector<double>&             ty          ) override
      {
         size_t q = order_up + 1;
         size_t n = tx.size() / q;
   # ifndef NDEBUG
         size_t m = ty.size() / q;
         assert( call_id == 0 );
         assert( m == 1 );
         assert( m == select_y.size() );
   # endif
         // ok
         bool ok = order_up <= 1 && order_low <= order_up;
         if ( ! ok )
            return ok;
         //
         // sum = x_0^0 * x_0^0 + x_1^0 * x_1^0 + ...
         double sum = 0.0;
         for(size_t j = 0; j < n; ++j)
         {  double xj0 = tx[ j * q + 0];
            sum       += xj0 * xj0;
         }
         //
         // ty[0] = sum
         if( order_low <= 0 )
            ty[0] = sum;
         if( order_up < 1 )
            return ok;

         // sum = x_0^0 * x_0^1 + x_1^0 ^ x_1^1 + ...
         sum   = 0.0;
         for(size_t j = 0; j < n; ++j)
         {  double xj0 = tx[ j * q + 0];
            double xj1 = tx[ j * q + 1];
            sum       += xj0 * xj1;
         }
         // ty[1] = 2.0 * sum
         assert( order_up == 1 );
         ty[1] = 2.0 * sum;
         return ok;
      }