Prev

Next

Index-> contents reference index search external

Up-> cppad_swig library a_fun a_fun_reverse a_fun_reverse_xam.pm

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.pm

Headings

package a_fun_reverse_xam;
sub a_fun_reverse_xam() {
     # check for standard perl programming conventions
     use strict;
     use warnings;
     #
     # load the Cppad Swig library
     use pm_cppad;
     #
     # initilaize return variable
     my $ok = 1;
     # ---------------------------------------------------------------------
     # number of dependent and independent variables
     my $n_dep = 1;
     my $n_ind = 3;
     #
     # create the independent variables ax
     my $xp = new pm_cppad::vec_double($n_ind);
     for(my $i = 0; $i < $n_ind ; $i++) {
          $xp->set($i, $i);
     }
     my $ax = pm_cppad::independent($xp);
     #
     # create dependent variables ay with ay0 = ax_0 * ax_1 * ax_2
     my $ax_0 = $ax->get(0);
     my $ax_1 = $ax->get(1);
     my $ax_2 = $ax->get(2);
     my $ay = new pm_cppad::vec_a_double($n_dep);
     $ay->set(0, $ax_0 * $ax_1 * $ax_2);
     #
     # define af corresponding to f(x) = x_0 * x_1 * x_2
     my $af = new pm_cppad::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
     my $p = 0;
     $xp->set(0, 2.0);
     $xp->set(1, 3.0);
     $xp->set(2, 4.0);
     my $yp = $af->forward($p, $xp);
     $ok = $ok && $yp->get(0) == 24.0;
     # -----------------------------------------------------------------------
     # first order reverse (derivative of zero order forward)
     # define G( Y ) = y_0 = x_0 * x_1 * x_2
     my $q = 1;
     my $yq1 = new pm_cppad::vec_double($n_dep);
     $yq1->set(0, 1.0);
     my $xq1 = $af->reverse($q, $yq1);
     # partial G w.r.t x_0
     $ok = $ok && $xq1->get(0) == 3.0 * 4.0 ;
     # partial G w.r.t x_1
     $ok = $ok && $xq1->get(1) == 2.0 * 4.0 ;
     # partial G w.r.t x_2
     $ok = $ok && $xq1->get(2) == 2.0 * 3.0 ;
     # -----------------------------------------------------------------------
     # first order forward mode
     $p = 1;
     $xp->set(0, 1.0);
     $xp->set(1, 1.0);
     $xp->set(2, 1.0);
     $yp = $af->forward($p, $xp);
     $ok = $ok && $yp->get(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;
     my $yq2 = new pm_cppad::vec_double($n_dep * $q);
     $yq2->set(0 * $q + 0, 0.0); # partial of G w.r.t y_0^0
     $yq2->set(0 * $q + 1, 1.0); # partial of G w.r.t y_0^1
     my $xq2 = $af->reverse($q, $yq2);
     # partial G w.r.t x_0^0
     $ok = $ok && $xq2->get(0 * $q + 0) == 3.0 + 4.0;
     # partial G w.r.t x_1^0
     $ok = $ok && $xq2->get(1 * $q + 0) == 2.0 + 4.0;
     # partial G w.r.t x_2^0
     $ok = $ok && $xq2->get(2 * $q + 0) == 2.0 + 3.0;
     # -----------------------------------------------------------------------
     #
     return( $ok );
}