Changeset 2859 for trunk/test_more


Ignore:
Timestamp:
May 28, 2013 2:03:21 AM (7 years ago)
Author:
bradbell
Message:

merge in changes from branches/atomic; see bin/svn_merge.sh

Location:
trunk
Files:
11 edited
2 copied

Legend:

Unmodified
Added
Removed
  • trunk

  • trunk/test_more/CMakeLists.txt

    r2756 r2859  
    8888        base_alloc.cpp
    8989        check_simple_vector.cpp
     90        checkpoint.cpp
    9091        compare.cpp
    9192        compare_change.cpp
     
    105106        for_sparse_jac.cpp
    106107        forward.cpp
     108        forward_mul.cpp
    107109        from_base.cpp
    108110        fun_check.cpp
  • trunk/test_more/for_sparse_jac.cpp

    r2506 r2859  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    394394}
    395395
     396bool case_four()
     397{       bool ok = true;
     398        using namespace CppAD;
     399
     400        // dimension of the domain space
     401        size_t n = 2;
     402
     403        // dimension of the range space
     404        size_t m = 3;
     405
     406        // independent variable vector
     407        CPPAD_TESTVECTOR(AD<double>) X(n);
     408        X[0] = 2.;
     409        X[1] = 3.;
     410        Independent(X);
     411
     412        // dependent variable vector
     413        CPPAD_TESTVECTOR(AD<double>) Y(m);
     414
     415        // check results vector
     416        CPPAD_TESTVECTOR( bool )       Check(m * n);
     417
     418        // initialize index into Y
     419        size_t index = 0;
     420
     421        // Y[0] only depends on X[0];
     422        Y[index]             = pow(X[0], 2.);
     423        Check[index * n + 0] = true;
     424        Check[index * n + 1] = false;
     425        index++;
     426
     427        // Y[1] depends on X[1]
     428        Y[index]             = pow(2., X[1]);
     429        Check[index * n + 0] = false;
     430        Check[index * n + 1] = true;
     431        index++;
     432
     433        // Y[2] depends on X[0] and X[1]
     434        Y[index]             = pow(X[0], X[1]);
     435        Check[index * n + 0] = true;
     436        Check[index * n + 1] = true;
     437        index++;
     438
     439        // check final index
     440        assert( index == m );
     441
     442        // create function object F : X -> Y
     443        ADFun<double> F(X, Y);
     444
     445        // -----------------------------------------------------------------
     446        // dependency matrix for the identity function
     447        CPPAD_TESTVECTOR( bool ) Px(n * n);
     448        size_t i, j;
     449        for(i = 0; i < n; i++)
     450        {       for(j = 0; j < n; j++)
     451                        Px[ i * n + j ] = false;
     452                Px[ i * n + i ] = true;
     453        }
     454
     455        // evaluate the dependency matrix for F(X(x))
     456        bool transpose = true;
     457        CPPAD_TESTVECTOR( bool ) Py(n * m);
     458        Py = F.ForSparseJac(n, Px, transpose);
     459
     460        // check values
     461        for(i = 0; i < m; i++)
     462        {       for(j = 0; j < n; j++)
     463                        ok &= (Py[j * m + i] == Check[i * n + j]);
     464        }       
     465
     466        // ---------------------------------------------------------
     467        // dependency matrix for the identity function
     468        CPPAD_TESTVECTOR(std::set<size_t>) Sx(n);
     469        for(i = 0; i < n; i++)
     470        {       assert( Sx[i].empty() );
     471                Sx[i].insert(i);
     472        }
     473
     474        // evaluate the dependency matrix for F(X(x))
     475        CPPAD_TESTVECTOR(std::set<size_t>) Sy(n);
     476        Sy = F.ForSparseJac(n, Sx, transpose);
     477
     478        // check values
     479        bool found;
     480        for(i = 0; i < m; i++)
     481        {       for(j = 0; j < n; j++)
     482                {       found = Sy[j].find(i) != Sy[j].end();
     483                        ok &= (found == Check[i * n + j]);
     484                }
     485        }       
     486
     487        return ok;
     488}
     489
    396490} // End empty namespace
    397491
     
    402496        ok &= case_two();
    403497        ok &= case_three();
     498        ok &= case_four();
    404499
    405500        return ok;
  • trunk/test_more/makefile.am

    r2756 r2859  
    101101        base_alloc.cpp \
    102102        check_simple_vector.cpp \
     103        checkpoint.cpp \
    103104        compare.cpp \
    104105        compare_change.cpp \
     
    119120        for_sparse_jac.cpp \
    120121        forward.cpp \
     122        forward_mul.cpp \
    121123        from_base.cpp \
    122124        fun_check.cpp \
  • trunk/test_more/makefile.in

    r2778 r2859  
    7979        alloc_openmp.cpp test_more.cpp abs.cpp acos.cpp asin.cpp \
    8080        assign.cpp add.cpp add_eq.cpp add_zero.cpp atan.cpp atan2.cpp \
    81         base_alloc.cpp check_simple_vector.cpp compare.cpp \
    82         compare_change.cpp cond_exp.cpp cond_exp_ad.cpp copy.cpp \
    83         cos.cpp cosh.cpp dbl_epsilon.cpp div.cpp div_eq.cpp \
     81        base_alloc.cpp check_simple_vector.cpp checkpoint.cpp \
     82        compare.cpp compare_change.cpp cond_exp.cpp cond_exp_ad.cpp \
     83        copy.cpp cos.cpp cosh.cpp dbl_epsilon.cpp div.cpp div_eq.cpp \
    8484        div_zero_one.cpp erf.cpp exp.cpp extern_value.cpp \
    8585        extern_value.hpp for_hess.cpp for_sparse_jac.cpp forward.cpp \
    86         from_base.cpp fun_check.cpp jacobian.cpp limits.cpp log.cpp \
    87         log10.cpp mul.cpp mul_eq.cpp mul_level.cpp mul_zero_one.cpp \
    88         ndebug.cpp near_equal_ext.cpp neg.cpp ode_err_control.cpp \
    89         optimize.cpp parameter.cpp poly.cpp pow.cpp pow_int.cpp \
    90         print_for.cpp romberg_one.cpp rosen_34.cpp runge_45.cpp \
    91         reverse.cpp rev_sparse_hes.cpp rev_sparse_jac.cpp rev_two.cpp \
    92         simple_vector.cpp sin.cpp sin_cos.cpp sinh.cpp \
     86        forward_mul.cpp from_base.cpp fun_check.cpp jacobian.cpp \
     87        limits.cpp log.cpp log10.cpp mul.cpp mul_eq.cpp mul_level.cpp \
     88        mul_zero_one.cpp ndebug.cpp near_equal_ext.cpp neg.cpp \
     89        ode_err_control.cpp optimize.cpp parameter.cpp poly.cpp \
     90        pow.cpp pow_int.cpp print_for.cpp romberg_one.cpp rosen_34.cpp \
     91        runge_45.cpp reverse.cpp rev_sparse_hes.cpp rev_sparse_jac.cpp \
     92        rev_two.cpp simple_vector.cpp sin.cpp sin_cos.cpp sinh.cpp \
    9393        sparse_hessian.cpp sparse_jacobian.cpp sparse_vec_ad.cpp \
    9494        sqrt.cpp std_math.cpp sub.cpp sub_eq.cpp sub_zero.cpp tan.cpp \
     
    103103        add_eq.$(OBJEXT) add_zero.$(OBJEXT) atan.$(OBJEXT) \
    104104        atan2.$(OBJEXT) base_alloc.$(OBJEXT) \
    105         check_simple_vector.$(OBJEXT) compare.$(OBJEXT) \
    106         compare_change.$(OBJEXT) cond_exp.$(OBJEXT) \
     105        check_simple_vector.$(OBJEXT) checkpoint.$(OBJEXT) \
     106        compare.$(OBJEXT) compare_change.$(OBJEXT) cond_exp.$(OBJEXT) \
    107107        cond_exp_ad.$(OBJEXT) copy.$(OBJEXT) cos.$(OBJEXT) \
    108108        cosh.$(OBJEXT) dbl_epsilon.$(OBJEXT) div.$(OBJEXT) \
    109109        div_eq.$(OBJEXT) div_zero_one.$(OBJEXT) erf.$(OBJEXT) \
    110110        exp.$(OBJEXT) extern_value.$(OBJEXT) for_hess.$(OBJEXT) \
    111         for_sparse_jac.$(OBJEXT) forward.$(OBJEXT) from_base.$(OBJEXT) \
    112         fun_check.$(OBJEXT) jacobian.$(OBJEXT) limits.$(OBJEXT) \
    113         log.$(OBJEXT) log10.$(OBJEXT) mul.$(OBJEXT) mul_eq.$(OBJEXT) \
     111        for_sparse_jac.$(OBJEXT) forward.$(OBJEXT) \
     112        forward_mul.$(OBJEXT) from_base.$(OBJEXT) fun_check.$(OBJEXT) \
     113        jacobian.$(OBJEXT) limits.$(OBJEXT) log.$(OBJEXT) \
     114        log10.$(OBJEXT) mul.$(OBJEXT) mul_eq.$(OBJEXT) \
    114115        mul_level.$(OBJEXT) mul_zero_one.$(OBJEXT) ndebug.$(OBJEXT) \
    115116        near_equal_ext.$(OBJEXT) neg.$(OBJEXT) \
     
    411412        base_alloc.cpp \
    412413        check_simple_vector.cpp \
     414        checkpoint.cpp \
    413415        compare.cpp \
    414416        compare_change.cpp \
     
    429431        for_sparse_jac.cpp \
    430432        forward.cpp \
     433        forward_mul.cpp \
    431434        from_base.cpp \
    432435        fun_check.cpp \
     
    545548@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/base_alloc.Po@am__quote@
    546549@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/check_simple_vector.Po@am__quote@
     550@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/checkpoint.Po@am__quote@
    547551@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/compare.Po@am__quote@
    548552@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/compare_change.Po@am__quote@
     
    562566@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/for_sparse_jac.Po@am__quote@
    563567@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/forward.Po@am__quote@
     568@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/forward_mul.Po@am__quote@
    564569@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/from_base.Po@am__quote@
    565570@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/fun_check.Po@am__quote@
  • trunk/test_more/optimize.cpp

    r2506 r2859  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    1010Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
    1111-------------------------------------------------------------------------- */
     12// 2DO: Test that optimize.hpp use of base_atomic<Base>::rev_sparse_jac works.
    1213
    1314# include <limits>
    1415# include <cppad/cppad.hpp>
     16
    1517
    1618namespace {
     
    955957        }
    956958
    957         bool user_atomic_forward(
     959        bool old_atomic_forward(
    958960                size_t                         id ,
    959961                size_t                          k ,
     
    981983        }
    982984
    983         bool user_atomic_reverse(
     985        bool old_atomic_reverse(
    984986                size_t                         id ,
    985987                size_t                          k ,
     
    992994        {       return false; }
    993995
    994         bool user_atomic_for_jac_sparse(
     996        bool old_atomic_for_jac_sparse(
    995997                size_t                                  id ,
    996998                size_t                                   n ,
     
    10011003        {       return false; }
    10021004
    1003         bool user_atomic_rev_jac_sparse(
     1005        bool old_atomic_rev_jac_sparse(
    10041006                size_t                                  id ,
    10051007                size_t                                   n ,
     
    10221024        }
    10231025
    1024         bool user_atomic_rev_hes_sparse(
     1026        bool old_atomic_rev_hes_sparse(
    10251027                size_t                                  id ,
    10261028                size_t                                   n ,
     
    10351037
    10361038        CPPAD_USER_ATOMIC(
    1037                 my_user_atomic             ,
     1039                my_old_atomic             ,
    10381040                CppAD::vector              ,
    10391041                double                     ,
    1040                 user_atomic_forward        ,
    1041                 user_atomic_reverse        ,
    1042                 user_atomic_for_jac_sparse ,
    1043                 user_atomic_rev_jac_sparse ,
    1044                 user_atomic_rev_hes_sparse
     1042                old_atomic_forward        ,
     1043                old_atomic_reverse        ,
     1044                old_atomic_for_jac_sparse ,
     1045                old_atomic_rev_jac_sparse ,
     1046                old_atomic_rev_hes_sparse
    10451047        )
    10461048
    1047         bool user_atomic_test(void)
     1049        bool old_atomic_test(void)
    10481050        {       bool ok = true;
    10491051
     
    10591061                size_t id = 0;
    10601062                // first call should stay in the tape
    1061                 my_user_atomic(id++, ax, ay);
     1063                my_old_atomic(id++, ax, ay);
    10621064                // second call will not get used
    1063                 my_user_atomic(id++, ax, az);
     1065                my_old_atomic(id++, ax, az);
    10641066                // create function
    10651067                CppAD::ADFun<double> g(ax, ay);
     
    11331135        // check that CondExp properly detects dependencies
    11341136        ok     &= cond_exp_depend();
    1135         // check user_atomic functions
    1136         ok     &= user_atomic_test();
     1137        // check old_atomic functions
     1138        ok     &= old_atomic_test();
    11371139        // case where results are not identically equal
    11381140        ok     &= not_identically_equal();
  • trunk/test_more/rev_sparse_hes.cpp

    r2506 r2859  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    441441}
    442442
     443bool case_six()
     444{       bool ok = true;
     445        using namespace CppAD;
     446
     447        // dimension of the domain space
     448        size_t n = 3;
     449
     450        // dimension of the range space
     451        size_t m = 1;
     452
     453        // independent variable vector
     454        CPPAD_TESTVECTOR(AD<double>) X(n);
     455        X[0] = 0.;
     456        X[1] = 1.;
     457        X[2] = 2.;
     458        Independent(X);
     459        // y = z[j] where j might be zero or one.
     460        CPPAD_TESTVECTOR(AD<double>) Y(m);
     461        Y[0]  =  X[1] * X[2];
     462        // create function object F : X -> Y
     463        ADFun<double> F(X, Y);
     464
     465        // sparsity pattern for hessian of F^2
     466        CPPAD_TESTVECTOR(bool) F2(n * n);
     467        F2[0 * n + 0] = false; // partial w.r.t x[0], x[0]
     468        F2[0 * n + 1] = false; // partial w.r.t x[0], x[1]
     469        F2[0 * n + 2] = false; // partial w.r.t x[0], x[2]
     470
     471        F2[1 * n + 0] = false; // partial w.r.t x[1], x[0]
     472        F2[1 * n + 1] = false; // partial w.r.t x[1], x[1]
     473        F2[1 * n + 2] = true;  // partial w.r.t x[1], x[2]
     474
     475        F2[2 * n + 0] = false; // partial w.r.t x[2], x[0]
     476        F2[2 * n + 1] = true;  // partial w.r.t x[2], x[1]
     477        F2[2 * n + 2] = false; // partial w.r.t x[2], x[2]
     478
     479
     480        // choose a non-symmetric sparsity patter for R
     481        CPPAD_TESTVECTOR( bool ) r(n * n);
     482        size_t i, j, k;
     483        for(i = 0; i < n; i++)
     484        {       for(j = 0; j < n; j++)
     485                        r[ i * n + j ] = false;
     486                j = n - i - 1;
     487                r[ j * n + j ] = true;
     488        }
     489
     490        // sparsity pattern for H^T
     491        CPPAD_TESTVECTOR(bool) Check(n * n);
     492        for(i = 0; i < n; i++)
     493        {       for(j = 0; j < n; j++)
     494                {       Check[ i * n + j] = false;
     495                        for(k = 0; k < n; k++)
     496                                Check[i * n + j] |= F2[i * n + k] & r[ k * n + j];
     497                }
     498        }
     499
     500        // compute the reverse Hessian sparsity pattern for F^2 * R
     501        F.ForSparseJac(n, r);
     502        CPPAD_TESTVECTOR( bool ) s(m), h(n * n);
     503        s[0] = 1.;
     504        bool transpose = true;
     505        h = F.RevSparseHes(n, s, transpose);
     506
     507        // check values
     508        for(i = 0; i < n; i++)
     509        {       for(j = 0; j < n; j++)
     510                        ok &= (h[i * n + j] == Check[i * n + j]);
     511        }       
     512
     513        // compute the reverse Hessian sparsity pattern for R^T * F^2
     514        transpose = false;
     515        h = F.RevSparseHes(n, s, transpose);
     516
     517        // check values
     518        for(i = 0; i < n; i++)
     519        {       for(j = 0; j < n; j++)
     520                        ok &= (h[j * n + i] == Check[i * n + j]);
     521        }       
     522
     523        return ok;
     524}
     525
    443526} // End of empty namespace
    444527
     
    451534        ok &= case_four();
    452535        ok &= case_five();
    453 
    454         return ok;
    455 }
     536        ok &= case_six();
     537
     538        return ok;
     539}
  • trunk/test_more/rev_sparse_jac.cpp

    r2811 r2859  
    423423        r = F.RevSparseJac(q, s);
    424424
    425         ok &= r.size() == q * n;
     425        ok &= size_t( r.size() ) == q * n;
    426426        ok &= r[0] == false;
    427427        ok &= r[1] == true;
     
    429429        return ok;
    430430}
     431
     432bool case_five()
     433{       bool ok = true;
     434        using namespace CppAD;
     435
     436        // dimension of the domain space
     437        size_t n = 2;
     438
     439        // dimension of the range space
     440        size_t m = 3;
     441
     442        // independent variable vector
     443        CPPAD_TESTVECTOR(AD<double>) X(n);
     444        X[0] = 2.;
     445        X[1] = 3.;
     446        Independent(X);
     447
     448        // dependent variable vector
     449        CPPAD_TESTVECTOR(AD<double>) Y(m);
     450
     451        // check results vector
     452        CPPAD_TESTVECTOR( bool )       Check(m * n);
     453
     454        // initialize index into Y
     455        size_t index = 0;
     456
     457        // Y[0] only depends on X[0];
     458        Y[index]             = pow(X[0], 2.);
     459        Check[index * n + 0] = true;
     460        Check[index * n + 1] = false;
     461        index++;
     462
     463        // Y[1] depends on X[1]
     464        Y[index]             = pow(2., X[1]);
     465        Check[index * n + 0] = false;
     466        Check[index * n + 1] = true;
     467        index++;
     468
     469        // Y[2] depends on X[0] and X[1]
     470        Y[index]             = pow(X[0], X[1]);
     471        Check[index * n + 0] = true;
     472        Check[index * n + 1] = true;
     473        index++;
     474
     475        // check final index
     476        assert( index == m );
     477
     478        // create function object F : X -> Y
     479        ADFun<double> F(X, Y);
     480
     481        // -----------------------------------------------------------------
     482        // dependency matrix for the identity function
     483        CPPAD_TESTVECTOR( bool ) Py(m * m);
     484        size_t i, j;
     485        for(i = 0; i < m; i++)
     486        {       for(j = 0; j < m; j++)
     487                        Py[ i * m + j ] = (i == j);
     488        }
     489
     490        // evaluate the dependency matrix for F(x)
     491        bool transpose = true;
     492        CPPAD_TESTVECTOR( bool ) Px(n * m);
     493        Px = F.RevSparseJac(m, Py, transpose);
     494
     495        // check values
     496        for(i = 0; i < m; i++)
     497        {       for(j = 0; j < n; j++)
     498                        ok &= (Px[j * m + i] == Check[i * n + j]);
     499        }       
     500
     501        // ---------------------------------------------------------
     502        // dependency matrix for the identity function
     503        CPPAD_TESTVECTOR(std::set<size_t>) Sy(m);
     504        for(i = 0; i < m; i++)
     505        {       assert( Sy[i].empty() );
     506                Sy[i].insert(i);
     507        }
     508
     509        // evaluate the dependency matrix for F(x)
     510        CPPAD_TESTVECTOR(std::set<size_t>) Sx(n);
     511        Sx = F.RevSparseJac(m, Sy, transpose);
     512
     513        // check values
     514        bool found;
     515        for(i = 0; i < m; i++)
     516        {       for(j = 0; j < n; j++)
     517                {       found = Sx[j].find(i) != Sx[j].end();
     518                        ok &= (found == Check[i * n + j]);
     519                }
     520        }       
     521
     522        return ok;
     523}
     524
     525
    431526
    432527} // END empty namespace
  • trunk/test_more/reverse.cpp

    r2506 r2859  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    1818namespace { // ----------------------------------------------------------
    1919
    20 bool Reverse(void)
     20bool reverse_one(void)
    2121{       bool ok = true;
    2222
     
    204204        return ok;
    205205}
     206/*
     207$comment reverse_any.cpp$$
     208$spell
     209        Taylor
     210$$
     211
     212$section Reverse Mode General Case: Example and Test$$
     213
     214$index general, reverse example$$
     215$index reverse, general example$$
     216$index example, general reverse$$
     217$index test, general reverse$$
     218
     219$index composition, example$$
     220$index example, composition$$
     221$index test, composition$$
     222
     223$head Purpose$$
     224Break a derivative computation into pieces and only store values at the
     225interface of the pieces.
     226In actual applications, there may be many functions, but
     227for this example there are only two.
     228The functions
     229$latex F : \B{R}^2 \rightarrow \B{R}^2$$
     230and
     231$latex G : \B{R}^2 \rightarrow \B{R}^2$$
     232defined by
     233$latex \[
     234        F(x) = \left( \begin{array}{c} x_0 x_1   \\ x_1 - x_0 \end{array} \right)
     235        \; , \;
     236        G(y) = \left( \begin{array}{c} y_0 - y_1 \\ y_1  y_0   \end{array} \right)
     237\] $$
     238Another difference is that in actual applications,
     239the memory corresponding to function objects not currently being used
     240is sometimes returned to the system (see $cref checkpoint.cpp$$).
     241
     242$head Processing Steps$$
     243We apply reverse mode to compute the derivative of
     244$latex H : \B{R}^2 \rightarrow \B{R}$$
     245is defined by
     246$latex \[
     247\begin{array}{rcl}
     248        H(x)
     249        & = & G_0 [ F(x) ] + G_1 [ F(x)  ]
     250        \\
     251        & = & x_0 x_1 - ( x_1 - x_0 ) + x_0 x_1 ( x_1 - x_0 )
     252        \\
     253        & = & x_0 x_1 ( 1 - x_0 + x_1 ) - x_1 + x_0
     254\end{array}
     255\] $$
     256Given the zero and first order Taylor coefficients
     257$latex x^{(0)} $$ and $latex x^{(1)}$$,
     258we use $latex X(t)$$, $latex Y(t)$$ and $latex Z(t)$$
     259for the corresponding functions; i.e.,
     260$latex \[
     261\begin{array}{rcl}
     262        X(t) & = & x^{(0)} + x^{(1)} t
     263        \\
     264        Y(t) & = & F[X(t)] = y^{(0)} + y^{(1)} t  + O(t^2)
     265        \\
     266        Z(t) & = & G \{ F [ X(t) ] \} = z^{(0)} + z^{(1)} t  + O(t^2)
     267        \\
     268        h^{(0)} & = & z^{(0)}_0 + z^{(0)}_1
     269        \\
     270        h^{(1)} & = & z^{(1)}_0 + z^{(1)}_1
     271\end{array}
     272\] $$
     273Here are the processing steps:
     274$list number$$
     275Use forward mode on $latex F(x)$$ to compute
     276$latex y^{(0)}$$ and $latex y^{(1)}$$
     277$lnext
     278Use forward mode on $latex G(y)$$ to compute
     279$latex z^{(0)}$$ and $latex z^{(1)}$$
     280$lnext
     281Use reverse mode on $latex G(y)$$ to compute the derivative of
     282$latex h^{(k)}$$ with respect to
     283$latex y^{(0)}$$ and $latex y^{(1)}$$.
     284$lnext
     285Use reverse mode on $latex F(x)$$ to compute the derivative of
     286$latex h^{(k)}$$ with respect to
     287$latex x^{(0)}$$ and $latex x^{(1)}$$.
     288$lend
     289This uses the following relations for $latex k = 0 , 1$$:
     290$latex \[
     291\begin{array}{rcl}
     292        \partial_{x(0)} h^{(k)} [ x^{(0)} , x^{(1)} ]
     293        & = &
     294        \partial_{y(0)} h^{(k)} [ y^{(0)} , y^{(1)} ]
     295        \partial_{x(0)} y^{(0)} [ x^{(0)} , x^{(1)} ]
     296        \\
     297        & + &
     298        \partial_{y(1)} h^{(k)} [ y^{(0)} , y^{(1)} ]
     299        \partial_{x(0)} y^{(1)} [ x^{(0)} , x^{(1)} ]
     300        \\
     301        \partial_{x(1)} h^{(k)} [ x^{(0)} , x^{(1)} ]
     302        & = &
     303        \partial_{y(0)} h^{(k)} [ y^{(0)} , y^{(1)} ]
     304        \partial_{x(1)} y^{(0)} [ x^{(0)} , x^{(1)} ]
     305        \\
     306        & + &
     307        \partial_{y(1)} h^{(k)} [ y^{(0)} , y^{(1)} ]
     308        \partial_{x(1)} y^{(1)} [ x^{(0)} , x^{(1)} ]
     309\end{array}
     310\] $$
     311where $latex \partial_{x(0)}$$ denotes the partial with respect
     312to $latex x^{(0)}$$.
     313
     314$code
     315$comment%example/reverse_any.cpp%0%// BEGIN C++%// END C++%1%$$
     316$$
     317$end
     318*/
     319template <class Vector>
     320Vector F_reverse_mul(const Vector& x)
     321{       Vector y(2);
     322        y[0] = x[0] * x[1];
     323        y[1] = x[1] - x[0];
     324        return y;
     325}
     326template <class Vector>
     327Vector G_reverse_mul(const Vector& y)
     328{       Vector z(2);
     329        z[0] = y[0] - y[1];
     330        z[1] = y[1] * y[0];
     331        return z;
     332}
     333bool reverse_mul(void)
     334{
     335        bool ok = true;
     336     double eps = 10. * CppAD::numeric_limits<double>::epsilon();
     337
     338        using CppAD::AD;
     339        using CppAD::NearEqual;
     340        CppAD::ADFun<double> f, g;
     341
     342        // Record the function F(x)
     343        size_t n    = 2;
     344        CPPAD_TESTVECTOR(AD<double>) X(n), Y(n);
     345        X[0] = X[1] = 0.;
     346        CppAD::Independent(X);
     347        Y = F_reverse_mul(X);
     348        f.Dependent(X, Y);
     349
     350        // Record the function G(x)
     351        CPPAD_TESTVECTOR(AD<double>) Z(n);
     352        Y[0] = Y[1] = 0.;
     353        CppAD::Independent(Y);
     354        Z = G_reverse_mul(Y);
     355        g.Dependent(Y, Z);
     356
     357        // argument and function values
     358        CPPAD_TESTVECTOR(double) x0(n), y0(n), z0(n);
     359        x0[0] = 1.;
     360        x0[1] = 2.;
     361        y0    = f.Forward(0, x0);
     362        z0    = g.Forward(0, y0);
     363
     364        // check function value
     365        double check = x0[0] * x0[1] * (1. - x0[0] + x0[1]) - x0[1] + x0[0];
     366        double h0    = z0[0] + z0[1];
     367        ok          &= NearEqual(h0, check, eps, eps);
     368
     369        // first order Taylor coefficients
     370        CPPAD_TESTVECTOR(double) x1(n), y1(n), z1(n);
     371        x1[0] = 3.;
     372        x1[1] = 4.;
     373        y1    = f.Forward(1, x1);
     374        z1    = g.Forward(1, y1);
     375
     376        // check first order Taylor coefficients
     377        check     = x0[0] * x0[1] * (- x1[0] + x1[1]) - x1[1] + x1[0];
     378        check    += x1[0] * x0[1] * (1. - x0[0] + x0[1]);
     379        check    += x0[0] * x1[1] * (1. - x0[0] + x0[1]);
     380        double h1 = z1[0] + z1[1];
     381        ok       &= NearEqual(h1, check, eps, eps);
     382
     383        // ----------------------------------------------------------------
     384        // dw^0 (y) = \partial_y^0 h^0 (y)
     385        // dw^1 (y) = \partial_y^1 h^0 (y)
     386        size_t p = 2;
     387        CPPAD_TESTVECTOR(double) w(n*p), dw(n*p);
     388        w[0*p+0] = 1.; // coefficient for z^0_0
     389        w[1*p+0] = 1.; // coefficient for z^0_1
     390        w[0*p+1] = 0.; // coefficient for z^1_0
     391        w[1*p+1] = 0.; // coefficient for z^1_1
     392        dw       = g.Reverse(p, w);
     393
     394        // dv^0 = dw^0 * \partial_x^0 y^0 (x) + dw^1 * \partial_x^0 y^1 (x) 
     395        // dv^1 = dw^0 * \partial_x^1 y^0 (x) + dw^1 * \partial_x^1 y^1 (x) 
     396        CPPAD_TESTVECTOR(double) dv(n*p);
     397        dv   = f.Reverse(p, dw);
     398
     399        // check partial of h^0 w.r.t x^0_0
     400        check  = x0[1] * (1. - x0[0] + x0[1]) + 1.;
     401        check -= x0[0] * x0[1];
     402        ok    &= NearEqual(dv[0*p+0], check, eps, eps);
     403
     404        // check partial of h^0 w.r.t x^0_1
     405        check  = x0[0] * (1. - x0[0] + x0[1]) - 1.;
     406        check += x0[0] * x0[1];
     407        ok    &= NearEqual(dv[1*p+0], check, eps, eps);
     408
     409        // check partial of h^0 w.r.t x^1_0 and x^1_1
     410        check  = 0.;
     411        ok    &= NearEqual(dv[0*p+1], check, eps, eps);
     412        ok    &= NearEqual(dv[1*p+1], check, eps, eps);
     413
     414        // ----------------------------------------------------------------
     415        // dw^0 (y) = \partial_y^0 h^1 (y)
     416        // dw^1 (y) = \partial_y^1 h^1 (y)
     417        w[0*p+0] = 0.; // coefficient for z^0_0
     418        w[1*p+0] = 0.; // coefficient for z^0_1
     419        w[0*p+1] = 1.; // coefficient for z^1_0
     420        w[1*p+1] = 1.; // coefficient for z^1_1
     421        dw       = g.Reverse(p, w);
     422
     423        // dv^0 = dw^0 * \partial_x^0 y^0 (x) + dw^1 * \partial_x^0 y^1 (x) 
     424        // dv^1 = dw^0 * \partial_x^1 y^0 (x) + dw^1 * \partial_x^1 y^1 (x) 
     425        dv   = f.Reverse(p, dw);
     426
     427        // check partial of h^1 w.r.t x^0_0
     428        check  = x0[1] * (- x1[0] + x1[1]);
     429        check -= x1[0] * x0[1];
     430        check += x1[1] * (1. - x0[0] + x0[1]) - x0[0] * x1[1];
     431        ok    &= NearEqual(dv[0*p+0], check, eps, eps);
     432
     433        // check partial of h^1 w.r.t x^0_1
     434        check  = x0[0] * (- x1[0] + x1[1]);
     435        check += x1[0] * (1. - x0[0] + x0[1]) + x1[0] * x0[1];
     436        check += x0[0] * x1[1];
     437        ok    &= NearEqual(dv[1*p+0], check, eps, eps);
     438
     439        // check partial of h^1 w.r.t x^1_0
     440        // (by reverse mode identity is equal to partial h^0 w.r.t. x^0_0)
     441        check  = 1. - x0[0] * x0[1];
     442        check += x0[1] * (1. - x0[0] + x0[1]);
     443        ok    &= NearEqual(dv[0*p+1], check, eps, eps);
     444
     445        // check partial of h^1 w.r.t x^1_1
     446        // (by reverse mode identity is equal to partial h^0 w.r.t. x^0_1)
     447        check  = x0[0] * x0[1] - 1.;
     448        check += x0[0] * (1. - x0[0] + x0[1]);
     449        ok    &= NearEqual(dv[1*p+1], check, eps, eps);
     450
     451        return ok;
     452}
     453// ----------------------------------------------------------------------------
    206454} // End empty namespace
    207455
     
    210458bool reverse(void)
    211459{       bool ok = true;
    212         ok &= Reverse();
     460        ok &= reverse_one();
     461        ok &= reverse_mul();
    213462
    214463        ok &= reverse_any_cases< CppAD::vector  <double> >();
  • trunk/test_more/test_more.cpp

    r2722 r2859  
    3131extern bool base_alloc_test(void);
    3232extern bool check_simple_vector(void);
     33extern bool checkpoint(void);
    3334extern bool Compare(void);
    3435extern bool CompareChange(void);
     
    4849extern bool for_sparse_jac(void);
    4950extern bool Forward(void);
     51extern bool forward_mul(void);
    5052extern bool FromBase(void);
    5153extern bool FunCheck(void);
     
    139141        ok &= Run( Atan2,           "Atan2"          );
    140142        ok &= Run( check_simple_vector, "check_simple_vector" );
     143        ok &= Run( checkpoint,      "checkpoint"     );
    141144        ok &= Run( Compare,         "Compare"        );
    142145        ok &= Run( CompareChange,   "CompareChange"  );
     
    155158        ok &= Run( for_sparse_jac,  "for_sparse_jac" );
    156159        ok &= Run( Forward,         "Forward"        );
     160        ok &= Run( forward_mul,     "forward_mul"    );
    157161        ok &= Run( FromBase,        "FromBase"       );
    158162        ok &= Run( FunCheck,        "FunCheck"       );
  • trunk/test_more/test_one.sh.in

    r2683 r2859  
    22# $Id$
    33# -----------------------------------------------------------------------------
    4 # CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     4# CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-13 Bradley M. Bell
    55#
    66# CppAD is distributed under multiple licenses. This distribution is under
     
    4545        -o test_one.exe
    4646        $cxxflags
    47         -std=c++98
     47        -std=c++11
    4848        -DCPPAD_ADOLC_TEST
    4949        -DCPPAD_OPENMP_TEST
Note: See TracChangeset for help on using the changeset viewer.