Changeset 1447 for trunk/test_more


Ignore:
Timestamp:
Jul 4, 2009 3:15:09 PM (11 years ago)
Author:
bradbell
Message:

trunk: Merge in branches/sweep from revision 1404 to 1446.

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk

  • trunk/test_more/cond_exp.cpp

    r1370 r1447  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    1717
    1818# include <cppad/cppad.hpp>
     19
     20namespace { // Begin empty namespace
    1921
    2022bool CondExp_pvvv(void)
     
    330332        return ok;
    331333}
     334
     335# include <limits>
     336bool SecondOrderReverse(void)
     337{       // Bradley M. Bell 2009-07-04
     338        // Reverse mode for CExpOp was only modifying the highest order partial
     339        // This test demonstrated the bug
     340        bool ok = true;
     341        using CppAD::AD;
     342        using CppAD::NearEqual;
     343        double eps = 10. * std::numeric_limits<double>::epsilon();
     344
     345        size_t n = 1;
     346        CPPAD_TEST_VECTOR< AD<double> > X(n);
     347        X[0] = 2.;
     348        CppAD::Independent(X);
     349
     350        size_t m = 2;
     351        CPPAD_TEST_VECTOR< AD<double> > Y(m);
     352
     353        AD<double> left = X[0];
     354        AD<double> right = X[0] * X[0];
     355        AD<double> exp_if_true  = left;
     356        AD<double> exp_if_false = right;
     357
     358        // Test of reverse mode using exp_if_true case
     359        // For this value of X, should be the same as Z = X[0]
     360        AD<double> Z = CondExpLt(left, right, exp_if_true, exp_if_false);
     361        Y[0] = Z * Z;
     362
     363        // Test of reverse mode  using exp_if_false case
     364        exp_if_false = left;
     365        exp_if_true  = right;
     366        Z            = CondExpGt(left, right, exp_if_true, exp_if_false);
     367        Y[1]         = Z * Z;
     368       
     369        CppAD::ADFun<double> f(X, Y);
     370
     371        // first order forward
     372        CPPAD_TEST_VECTOR<double> dx(n);
     373        size_t p = 1;
     374        dx[0]    = 1.;
     375        f.Forward(p, dx);
     376
     377        // second order reverse (test exp_if_true case)
     378        CPPAD_TEST_VECTOR<double> w(m), dw(2 * n);
     379        w[0] = 1.;
     380        w[1] = 0.;
     381        p    = 2;
     382        dw = f.Reverse(p, w);
     383
     384        // check first derivative in dw
     385        double check = 2. * Value( X[0] );
     386        ok &= NearEqual(dw[0], check, eps, eps);
     387
     388        // check second derivative in dw
     389        check = 2.;
     390        ok &= NearEqual(dw[1], check, eps, eps);
     391
     392        // test exp_if_false case
     393        w[0] = 0.;
     394        w[1] = 1.;
     395        p    = 2;
     396        dw = f.Reverse(p, w);
     397
     398        // check first derivative in dw
     399        check = 2. * Value( X[0] );
     400        ok &= NearEqual(dw[0], check, eps, eps);
     401
     402        // check second derivative in dw
     403        check = 2.;
     404        ok &= NearEqual(dw[1], check, eps, eps);
     405
     406        return ok;
     407}
     408
     409double Infinity(double zero)
     410{       return 1. / zero; }
     411
     412bool OldExample(void)
     413{       bool ok = true;
     414
     415        using CppAD::AD;
     416        using CppAD::NearEqual;
     417        using CppAD::log;
     418        using CppAD::abs;
     419        double eps = 100. * std::numeric_limits<double>::epsilon();
     420
     421        // domain space vector
     422        size_t n = 5;
     423        CPPAD_TEST_VECTOR< AD<double> > X(n);
     424        size_t j;
     425        for(j = 0; j < n; j++)
     426                X[j] = 1.;
     427
     428        // declare independent variables and start tape recording
     429        CppAD::Independent(X);
     430
     431        // sum with respect to j of log of absolute value of X[j]
     432        // sould be - infinity if any of the X[j] are zero
     433        AD<double> MinusInfinity = - Infinity(0.);
     434        AD<double> Sum           = 0.;
     435        AD<double> Zero(0);
     436        for(j = 0; j < n; j++)
     437        {       // if X[j] > 0
     438                Sum += CppAD::CondExpGt(X[j], Zero, log(X[j]),     Zero);
     439
     440                // if X[j] < 0
     441                Sum += CppAD::CondExpLt(X[j], Zero, log(-X[j]),    Zero);
     442
     443                // if X[j] == 0
     444                Sum += CppAD::CondExpEq(X[j], Zero, MinusInfinity, Zero);
     445        }
     446
     447        // range space vector
     448        size_t m = 1;
     449        CPPAD_TEST_VECTOR< AD<double> > Y(m);
     450        Y[0] = Sum;
     451
     452        // create f: X -> Y and stop tape recording
     453        CppAD::ADFun<double> f(X, Y);
     454
     455        // vectors for arguments to the function object f
     456        CPPAD_TEST_VECTOR<double> x(n);   // argument values
     457        CPPAD_TEST_VECTOR<double> y(m);   // function values
     458        CPPAD_TEST_VECTOR<double> w(m);   // function weights
     459        CPPAD_TEST_VECTOR<double> dw(n);  // derivative of weighted function
     460
     461        // a case where abs( x[j] ) > 0 for all j
     462        double check  = 0.;
     463        double sign   = 1.;
     464        for(j = 0; j < n; j++)
     465        {       sign *= -1.;
     466                x[j] = sign * double(j + 1);
     467                check += log( abs( x[j] ) );
     468        }
     469
     470        // function value
     471        y  = f.Forward(0, x);
     472        ok &= ( y[0] == check );
     473
     474        // compute derivative of y[0]
     475        w[0] = 1.;
     476        dw   = f.Reverse(1, w);
     477        for(j = 0; j < n; j++)
     478        {       if( x[j] > 0. )
     479                        ok &= NearEqual(dw[j], 1./abs( x[j] ), eps, eps);
     480                else    ok &= NearEqual(dw[j], -1./abs( x[j] ), eps, eps);
     481        }
     482
     483        // a case where x[0] is equal to zero
     484        sign = 1.;
     485        for(j = 0; j < n; j++)
     486        {       sign *= -1.;
     487                x[j] = sign * double(j);
     488        }
     489
     490        // function value
     491        y   = f.Forward(0, x);
     492        ok &= ( y[0] == -Infinity(0.) );
     493
     494        // compute derivative of y[0]
     495        w[0] = 1.;
     496        dw   = f.Reverse(1, w);
     497        for(j = 0; j < n; j++)
     498        {       if( x[j] > 0. )
     499                        ok &= NearEqual(dw[j], 1./abs( x[j] ), eps, eps);
     500                else if( x[j] < 0. )
     501                        ok &= NearEqual(dw[j], -1./abs( x[j] ), eps, eps);
     502                else
     503                {       // in this case computing dw[j] ends up multiplying
     504                        // -infinity * zero and hence results in Nan
     505                }
     506        }
     507       
     508        return ok;
     509}
     510} // end empty namespace
     511
    332512bool CondExp(void)
    333513{       bool ok  = true;
     
    336516        ok      &= CondExp_vvpv();
    337517        ok      &= CondExp_vvvp();
     518        ok      &= SecondOrderReverse();
     519        ok      &= OldExample();
    338520        return ok;
    339521}
  • trunk/test_more/div.cpp

    r1370 r1447  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    162162}
    163163
     164bool DivTestThree(void)
     165{       bool ok = true;
     166        using namespace CppAD;
     167
     168        // more testing of variable / variable case
     169        double x0 = 2.;
     170        double x1 = 3.;
     171        size_t n  = 2;
     172        CPPAD_TEST_VECTOR< AD<double> > X(n);
     173        X[0]      = x0;
     174        X[1]      = x1;
     175        Independent(X);
     176        size_t m  = 1;
     177        CPPAD_TEST_VECTOR< AD<double> > Y(m);
     178        Y[0]      = X[0] / X[1];
     179        ADFun<double> f(X, Y);
     180
     181        CPPAD_TEST_VECTOR<double> dx(n), dy(m);
     182        double check;
     183        dx[0] = 1.;
     184        dx[1] = 1.;
     185        dy    = f.Forward(1, dx);
     186        check = 1. / x1 - x0 / (x1 * x1);
     187        ok   &= NearEqual(dy[0], check, 1e-10 , 1e-10);
     188
     189        CPPAD_TEST_VECTOR<double> w(m), dw(n);
     190        w[0]  = 1.;
     191        dw    = f.Reverse(1, w);
     192        check = 1. / x1;
     193        ok   &= NearEqual(dw[0], check, 1e-10 , 1e-10);
     194        check = - x0 / (x1 * x1);
     195        ok   &= NearEqual(dw[1], check, 1e-10 , 1e-10);
     196
     197        return ok;
     198}
     199
    164200} // END empty namespace
    165201
  • trunk/test_more/sub.cpp

    r1370 r1447  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
     
    132132}
    133133
     134
     135bool SubTestThree(void)
     136{       bool ok = true;
     137        using namespace CppAD;
     138
     139        // special cases where tests above check OK and SubpvOp
     140        // implementation is known to be worng.
     141        // Probably two minuses make a plus.
     142        size_t n = 1;
     143        CPPAD_TEST_VECTOR< AD<double> > X(n);
     144        X[0] = 1.;
     145        Independent(X);
     146        size_t m = 1;
     147        CPPAD_TEST_VECTOR< AD<double> > Y(m);
     148        Y[0] = 1. - X[0];
     149        ADFun<double> f(X, Y);
     150       
     151        CPPAD_TEST_VECTOR<double> w(m), dw(n);
     152        w[0] = 1.;
     153        dw = f.Reverse(1, w);
     154        ok &= (dw[0] == -1.);
     155
     156        return ok;
     157}
     158
    134159} // END empty namespace
    135160
     
    138163        ok &= SubTestOne();
    139164        ok &= SubTestTwo();
     165        ok &= SubTestThree();
    140166        return ok;
    141167}
Note: See TracChangeset for help on using the changeset viewer.