Prev

Next

Index-> contents reference index search external

Up-> cppad_py library py_lib py_fun py_fun_reverse fun_reverse_xam.py

library-> py_lib cpp_lib

py_lib-> py_fun py_sparse py_utility more_py

py_fun-> py_independent py_abort_recording py_fun_ctor py_fun_property py_fun_new_dynamic py_fun_jacobian py_fun_hessian py_fun_forward py_fun_reverse py_fun_optimize

py_fun_reverse-> fun_reverse_xam.py

fun_reverse_xam.py

Headings

def fun_reverse_xam() :
     #
     import numpy
     import cppad_py
     #
     # initialize return variable
     ok = True
     # ---------------------------------------------------------------------
     # number of dependent and independent variables
     n_dep = 1
     n_ind = 3
     #
     # create the independent variables ax
     xp = numpy.empty(n_ind, dtype=float)
     for i in range( n_ind  ) :
          xp[i] = i
     #
     ax = cppad_py.independent(xp)
     #
     # create dependent variables ay with ay0 = ax_0 * ax_1 * ax_2
     ax_0  = ax[0]
     ax_1  = ax[1]
     ax_2  = ax[2]
     ay    = numpy.empty(n_dep, dtype=cppad_py.a_double)
     ay[0] = ax_0 * ax_1 * ax_2
     #
     # define af corresponding to f(x) = x_0 * x_1 * x_2
     f  = cppad_py.d_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
     p     = 0
     xp[0] = 2.0
     xp[1] = 3.0
     xp[2] = 4.0
     yp = f.forward(p, xp)
     ok = ok and yp[0] == 24.0
     # -----------------------------------------------------------------------
     # first order reverse (derivative of zero order forward)
     # define G( Y ) = y_0 = x_0 * x_1 * x_2
     m         = f.size_range()
     q         = 1
     yq1       = numpy.empty( (m, q), dtype=float)
     yq1[0, 0] = 1.0
     xq1       = f.reverse(q, yq1)
     # partial G w.r.t x_0
     ok = ok and xq1[0,0] == 3.0 * 4.0
     # partial G w.r.t x_1
     ok = ok and xq1[1,0] == 2.0 * 4.0
     # partial G w.r.t x_2
     ok = ok and xq1[2,0] == 2.0 * 3.0
     # -----------------------------------------------------------------------
     # first order forward mode
     p     = 1
     xp[0] = 1.0
     xp[1] = 1.0
     xp[2] = 1.0
     yp    = f.forward(p, xp)
     ok    = ok and 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
     yq2       = numpy.empty( (m, q), dtype=float)
     yq2[0, 0] = 0.0 # partial of G w.r.t y_0^0
     yq2[0, 1] = 1.0 # partial of G w.r.t y_0^1
     xq2       = f.reverse(q, yq2)
     # partial G w.r.t x_0^0
     ok = ok and xq2[0, 0] == 3.0 + 4.0
     # partial G w.r.t x_1^0
     ok = ok and xq2[1, 0] == 2.0 + 4.0
     # partial G w.r.t x_2^0
     ok = ok and xq2[2, 0] == 2.0 + 3.0
     # -----------------------------------------------------------------------
     af = cppad_py.a_fun(f)
     #
     # zero order forward
     axp   = numpy.empty(n_ind, dtype=cppad_py.a_double)
     p     = 0
     axp[0] = 2.0
     axp[1] = 3.0
     axp[2] = 4.0
     ayp = af.forward(p, axp)
     ok = ok and ayp[0] == cppad_py.a_double(24.0)
     #
     # first order reverse
     q          = 1
     ayq1       = numpy.empty( (m, q), dtype=cppad_py.a_double)
     ayq1[0, 0] = 1.0
     axq1       = af.reverse(q, ayq1)
     # partial G w.r.t x_0
     ok = ok and axq1[0,0] == cppad_py.a_double(3.0 * 4.0)
     # partial G w.r.t x_1
     ok = ok and axq1[1,0] == cppad_py.a_double(2.0 * 4.0)
     # partial G w.r.t x_2
     ok = ok and axq1[2,0] == cppad_py.a_double(2.0 * 3.0)
     #
     return( ok )
#