Prev

Next

Index-> contents reference index search external

Up-> cppad_swig library a_fun a_fun_reverse a_fun_reverse_xam.cpp

cppad_swig-> testing swig_xam.i library whats_new_2018

library-> module a_double vector a_fun sparse error xam.m4

a_fun-> independent abort_recording a_fun_ctor a_fun_jacobian a_fun_hessian a_fun_forward a_fun_reverse a_fun_optimize a_fun_property

a_fun_reverse-> a_fun_reverse_xam.cpp a_fun_reverse_xam.m a_fun_reverse_xam.pm a_fun_reverse_xam.py

a_fun_reverse_xam.cpp

Headings

# include <cstdio>
# include <string>
# include <cppad/swig/cppad_swig.hpp>

bool a_fun_reverse_xam(void) {
     using cppad_swig::a_double;
     using cppad_swig::vec_bool;
     using cppad_swig::vec_int;
     using cppad_swig::vec_double;
     using cppad_swig::vec_a_double;
     using cppad_swig::a_fun;
     using cppad_swig::sparse_rc;
     using cppad_swig::sparse_rcv;
     using cppad_swig::sparse_jac_work;
     using cppad_swig::sparse_hes_work;
     using std::string;
     //
     // initialize return variable
     bool ok = true;
     //------------------------------------------------------------------------
     // number of dependent and independent variables
     int n_dep = 1;
     int n_ind = 3;
     //
     // create the independent variables ax
     vec_double xp = cppad_swig::vec_double(n_ind);
     for(int i = 0; i < n_ind ; i++) {
          xp[i] = i;
     }
     vec_a_double ax = cppad_swig::independent(xp);
     //
     // create dependent variables ay with ay0 = ax_0 * ax_1 * ax_2
     a_double ax_0 = ax[0];
     a_double ax_1 = ax[1];
     a_double ax_2 = ax[2];
     vec_a_double ay = cppad_swig::vec_a_double(n_dep);
     ay[0] = ax_0 * ax_1 * ax_2;
     //
     // define af corresponding to f(x) = x_0 * x_1 * x_2
     a_fun af = cppad_swig::a_fun(ax, ay);
     // -----------------------------------------------------------------------
     // define          X(t) = (x_0 + t, x_1 + t, x_2 + t)
     // it follows that Y(t) = f(X(t)) = (x_0 + t) * (x_1 + t) * (x_2 + t)
     // and that       Y'(0) = x_1 * x_2 + x_0 * x_2 + x_0 * x_1
     // -----------------------------------------------------------------------
     // zero order forward mode
     int p = 0;
     xp[0] = 2.0;
     xp[1] = 3.0;
     xp[2] = 4.0;
     vec_double yp = af.forward(p, xp);
     ok = ok && yp[0] == 24.0;
     // -----------------------------------------------------------------------
     // first order reverse (derivative of zero order forward)
     // define G( Y ) = y_0 = x_0 * x_1 * x_2
     int q = 1;
     vec_double yq1 = cppad_swig::vec_double(n_dep);
     yq1[0] = 1.0;
     vec_double xq1 = af.reverse(q, yq1);
     // partial G w.r.t x_0
     ok = ok && xq1[0] == 3.0 * 4.0 ;
     // partial G w.r.t x_1
     ok = ok && xq1[1] == 2.0 * 4.0 ;
     // partial G w.r.t x_2
     ok = ok && xq1[2] == 2.0 * 3.0 ;
     // -----------------------------------------------------------------------
     // first order forward mode
     p = 1;
     xp[0] = 1.0;
     xp[1] = 1.0;
     xp[2] = 1.0;
     yp = af.forward(p, xp);
     ok = ok && yp[0] == 3.0*4.0 + 2.0*4.0 + 2.0*3.0;
     // -----------------------------------------------------------------------
     // second order reverse (derivative of first order forward)
     // define G( y_0^0 , y_0^1 ) = y_0^1
     // = x_1^0 * x_2^0  +  x_0^0 * x_2^0  +  x_0^0  *  x_1^0
     q = 2;
     vec_double yq2 = cppad_swig::vec_double(n_dep * q);
     yq2[0 * q + 0] = 0.0; // partial of G w.r.t y_0^0
     yq2[0 * q + 1] = 1.0; // partial of G w.r.t y_0^1
     vec_double xq2 = af.reverse(q, yq2);
     // partial G w.r.t x_0^0
     ok = ok && xq2[0 * q + 0] == 3.0 + 4.0;
     // partial G w.r.t x_1^0
     ok = ok && xq2[1 * q + 0] == 2.0 + 4.0;
     // partial G w.r.t x_2^0
     ok = ok && xq2[2 * q + 0] == 2.0 + 3.0;
     // -----------------------------------------------------------------------
     //
     return( ok );
}