source: trunk/test_more/sqrt.cpp @ 2506

Last change on this file since 2506 was 2506, checked in by bradbell, 8 years ago

Change Licenses: CPL-1.0 -> EPL-1.0, GPL-2.0->GPL-3.0

  • Property svn:keywords set to Id
File size: 5.1 KB
Line 
1/* $Id: sqrt.cpp 2506 2012-10-24 19:36:49Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Eclipse Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12
13/*
14Two old sqrt examples now used just for validation testing.
15*/
16# include <cppad/cppad.hpp>
17# include <cmath>
18
19namespace { // BEGIN empty namespace
20
21bool SqrtTestOne(void)
22{       bool ok = true;
23
24        using CppAD::sqrt;
25        using CppAD::pow;
26        using namespace CppAD;
27
28        // independent variable vector, indices, values, and declaration
29        CPPAD_TESTVECTOR(AD<double>) U(1);
30        size_t s = 0;
31        U[s]     = 4.;
32        Independent(U);
33
34        // dependent variable vector, indices, and values
35        CPPAD_TESTVECTOR(AD<double>) Z(2);
36        size_t x = 0;
37        size_t y = 1;
38        Z[x]     = sqrt(U[s]);
39        Z[y]     = sqrt(Z[x]);
40
41        // define f : U -> Z and vectors for derivative calculations
42        ADFun<double> f(U, Z);
43        CPPAD_TESTVECTOR(double) v( f.Domain() );
44        CPPAD_TESTVECTOR(double) w( f.Range() );
45
46        // check values
47        ok &= NearEqual(Z[x] , 2.,        1e-10 , 1e-10);
48        ok &= NearEqual(Z[y] , sqrt(2.),  1e-10 , 1e-10);
49
50        // forward computation of partials w.r.t. s
51        v[s] = 1.;
52        w    = f.Forward(1, v);
53        ok &= NearEqual(w[x], .5  * pow(4., -.5),   1e-10 , 1e-10); // dx/ds
54        ok &= NearEqual(w[y], .25 * pow(4., -.75),  1e-10 , 1e-10); // dy/ds
55
56        // reverse computation of partials of y
57        w[x] = 0.;
58        w[y] = 1.;
59        v    = f.Reverse(1,w);
60        ok &= NearEqual(v[s], .25 * pow(4., -.75),  1e-10 , 1e-10); // dy/ds
61
62        // forward computation of second partials w.r.t s
63        v[s] = 1.;
64        w    = f.Forward(1, v);
65        v[s] = 0.;
66        w    = f.Forward(2, v);
67        ok &= NearEqual(       // d^2 y / (ds ds)
68                2. * w[y] , 
69                -.75 * .25 * pow(4., -1.75),
70                1e-10 ,
71                1e-10 
72        ); 
73
74        // reverse computation of second partials of y
75        CPPAD_TESTVECTOR(double) r( f.Domain() * 2 );
76        w[x] = 0.;
77        w[y] = 1.;
78        r    = f.Reverse(2, w);
79        ok &= NearEqual(      // d^2 y / (ds ds)
80                r[2 * s + 1] , 
81                -.75 * .25 * pow(4., -1.75),
82                1e-10 ,
83                1e-10 
84        ); 
85
86        return ok;
87
88}
89bool SqrtTestTwo(void)
90{       bool ok = true;
91        using namespace CppAD;
92
93        // independent variable vector
94        CPPAD_TESTVECTOR(AD<double>) U(1);
95        U[0]     = 2.;
96        Independent(U);
97
98        // a temporary values
99        AD<double> x = U[0] * U[0]; 
100
101        // dependent variable vector
102        CPPAD_TESTVECTOR(AD<double>) Z(1);
103        Z[0] =  sqrt( x ); // z = sqrt( u * u )
104
105        // create f: U -> Z and vectors used for derivative calculations
106        ADFun<double> f(U, Z); 
107        CPPAD_TESTVECTOR(double) v(1);
108        CPPAD_TESTVECTOR(double) w(1);
109
110        // check value
111        ok &= NearEqual(U[0] , Z[0],  1e-10 , 1e-10);
112
113        // forward computation of partials w.r.t. u
114        size_t j;
115        size_t p     = 5;
116        double jfac  = 1.;
117        double value = 1.;
118        v[0]         = 1.;
119        for(j = 1; j < p; j++)
120        {       jfac *= j;
121                w     = f.Forward(j, v);       
122                ok &= NearEqual(jfac*w[0], value, 1e-10 , 1e-10); // d^jz/du^j
123                v[0]  = 0.;
124                value = 0.;
125        }
126
127        // reverse computation of partials of Taylor coefficients
128        CPPAD_TESTVECTOR(double) r(p); 
129        w[0]  = 1.;
130        r     = f.Reverse(p, w);
131        jfac  = 1.;
132        value = 1.;
133        for(j = 0; j < p; j++)
134        {       ok &= NearEqual(jfac*r[j], value, 1e-10 , 1e-10); // d^jz/du^j
135                jfac *= (j + 1);
136                value = 0.;
137        }
138
139        return ok;
140}
141bool SqrtTestThree(void)
142{       bool ok = true;
143
144        using CppAD::sqrt;
145        using CppAD::exp;
146        using namespace CppAD;
147
148        // independent variable vector, indices, values, and declaration
149        double x = 4.;
150        CPPAD_TESTVECTOR(AD<double>) X(1);
151        X[0]     = x;
152        Independent(X);
153
154        // dependent variable vector, indices, and values
155        CPPAD_TESTVECTOR(AD<double>) Y(1);
156        Y[0]     = sqrt( exp(X[0]) );
157
158        // define f : X -> Y and vectors for derivative calculations
159        ADFun<double> f(X, Y);
160
161        // forward computation of first Taylor coefficient
162        CPPAD_TESTVECTOR(double) x1( f.Domain() );
163        CPPAD_TESTVECTOR(double) y1( f.Range() );
164        x1[0] = 1.;
165        y1    = f.Forward(1, x1);
166        ok   &= NearEqual(y1[0], exp(x/2.)/2.,   1e-10 , 1e-10); 
167
168        // forward computation of second Taylor coefficient
169        CPPAD_TESTVECTOR(double) x2( f.Domain() );
170        CPPAD_TESTVECTOR(double) y2( f.Range() );
171        x2[0] = 0.;
172        y2    = f.Forward(2, x2);
173        ok   &= NearEqual(2.*y2[0] , exp(x/2.)/4., 1e-10 , 1e-10 ); 
174
175        // forward computation of third Taylor coefficient
176        CPPAD_TESTVECTOR(double) x3( f.Domain() );
177        CPPAD_TESTVECTOR(double) y3( f.Range() );
178        x3[0] = 0.;
179        y3    = f.Forward(3, x3);
180        ok   &= NearEqual(6.*y3[0] , exp(x/2.)/8., 1e-10 , 1e-10 ); 
181
182        // reverse computation of deritavitve of Taylor coefficients
183        CPPAD_TESTVECTOR(double) r( f.Domain() * 4 );
184        CPPAD_TESTVECTOR(double) w(1);
185        w[0] = 1.;
186        r    = f.Reverse(4, w);
187        ok   &= NearEqual(r[0], exp(x/2.)/2., 1e-10 , 1e-10); 
188        ok   &= NearEqual(r[1], exp(x/2.)/4., 1e-10 , 1e-10 ); 
189        ok   &= NearEqual(2.*r[2], exp(x/2.)/8., 1e-10 , 1e-10 ); 
190        ok   &= NearEqual(6.*r[3], exp(x/2.)/16., 1e-10 , 1e-10 ); 
191
192        return ok;
193
194}
195
196} // END empty namespace
197
198bool Sqrt(void)
199{       bool ok = true;
200        ok &= SqrtTestOne();
201        ok &= SqrtTestTwo(); 
202        ok &= SqrtTestThree(); 
203        return ok;
204}
Note: See TracBrowser for help on using the repository browser.