1 /* --------------------------------------------------------------------------
2 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
3
4 CppAD is distributed under the terms of the
5 Eclipse Public License Version 2.0.
6
7 This Source Code may also be made available under the following
8 Secondary License when the conditions for such availability set forth
9 in the Eclipse Public License, Version 2.0 are satisfied:
10 GNU General Public License, Version 2.0 or later.
11 ---------------------------------------------------------------------------- */
12
13 /*
14 Two old sqrt examples now used just for validation testing.
15 */
16 # include <cppad/cppad.hpp>
17 # include <cmath>
18
19 namespace { // BEGIN empty namespace
20
SqrtTestOne(void)21 bool SqrtTestOne(void)
22 { bool ok = true;
23 using CppAD::sqrt;
24 using CppAD::pow;
25 using namespace CppAD;
26 double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
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., eps99 , eps99);
48 ok &= NearEqual(Z[y] , sqrt(2.), eps99 , eps99);
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), eps99 , eps99); // dx/ds
54 ok &= NearEqual(w[y], .25 * pow(4., -.75), eps99 , eps99); // 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), eps99 , eps99); // 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 eps99 ,
71 eps99
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 eps99 ,
83 eps99
84 );
85
86 return ok;
87
88 }
SqrtTestTwo(void)89 bool SqrtTestTwo(void)
90 { bool ok = true;
91 using namespace CppAD;
92 double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
93
94 // independent variable vector
95 CPPAD_TESTVECTOR(AD<double>) U(1);
96 U[0] = 2.;
97 Independent(U);
98
99 // a temporary values
100 AD<double> x = U[0] * U[0];
101
102 // dependent variable vector
103 CPPAD_TESTVECTOR(AD<double>) Z(1);
104 Z[0] = sqrt( x ); // z = sqrt( u * u )
105
106 // create f: U -> Z and vectors used for derivative calculations
107 ADFun<double> f(U, Z);
108 CPPAD_TESTVECTOR(double) v(1);
109 CPPAD_TESTVECTOR(double) w(1);
110
111 // check value
112 ok &= NearEqual(U[0] , Z[0], eps99 , eps99);
113
114 // forward computation of partials w.r.t. u
115 size_t j;
116 size_t p = 5;
117 double jfac = 1.;
118 double value = 1.;
119 v[0] = 1.;
120 for(j = 1; j < p; j++)
121 { jfac *= double(j);
122 w = f.Forward(j, v);
123 ok &= NearEqual(w[0], value/jfac, eps99, eps99); // d^jz/du^j
124 v[0] = 0.;
125 value = 0.;
126 }
127
128 // reverse computation of partials of Taylor coefficients
129 CPPAD_TESTVECTOR(double) r(p);
130 w[0] = 1.;
131 r = f.Reverse(p, w);
132 jfac = 1.;
133 value = 1.;
134 for(j = 0; j < p; j++)
135 { ok &= NearEqual(r[j], value/jfac, eps99, eps99); // d^jz/du^j
136 jfac *= double(j + 1);
137 value = 0.;
138 }
139
140 return ok;
141 }
SqrtTestThree(void)142 bool SqrtTestThree(void)
143 { bool ok = true;
144 using CppAD::sqrt;
145 using CppAD::exp;
146 using namespace CppAD;
147 double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
148
149 // independent variable vector, indices, values, and declaration
150 double x = 4.;
151 CPPAD_TESTVECTOR(AD<double>) X(1);
152 X[0] = x;
153 Independent(X);
154
155 // dependent variable vector, indices, and values
156 CPPAD_TESTVECTOR(AD<double>) Y(1);
157 Y[0] = sqrt( exp(X[0]) );
158
159 // define f : X -> Y and vectors for derivative calculations
160 ADFun<double> f(X, Y);
161
162 // forward computation of first Taylor coefficient
163 CPPAD_TESTVECTOR(double) x1( f.Domain() );
164 CPPAD_TESTVECTOR(double) y1( f.Range() );
165 x1[0] = 1.;
166 y1 = f.Forward(1, x1);
167 ok &= NearEqual(y1[0], exp(x/2.)/2., eps99 , eps99);
168
169 // forward computation of second Taylor coefficient
170 CPPAD_TESTVECTOR(double) x2( f.Domain() );
171 CPPAD_TESTVECTOR(double) y2( f.Range() );
172 x2[0] = 0.;
173 y2 = f.Forward(2, x2);
174 ok &= NearEqual(2.*y2[0] , exp(x/2.)/4., eps99 , eps99 );
175
176 // forward computation of third Taylor coefficient
177 CPPAD_TESTVECTOR(double) x3( f.Domain() );
178 CPPAD_TESTVECTOR(double) y3( f.Range() );
179 x3[0] = 0.;
180 y3 = f.Forward(3, x3);
181 ok &= NearEqual(6.*y3[0] , exp(x/2.)/8., eps99 , eps99 );
182
183 // reverse computation of deritavitve of Taylor coefficients
184 CPPAD_TESTVECTOR(double) r( f.Domain() * 4 );
185 CPPAD_TESTVECTOR(double) w(1);
186 w[0] = 1.;
187 r = f.Reverse(4, w);
188 ok &= NearEqual(r[0], exp(x/2.)/2., eps99 , eps99);
189 ok &= NearEqual(r[1], exp(x/2.)/4., eps99 , eps99 );
190 ok &= NearEqual(2.*r[2], exp(x/2.)/8., eps99 , eps99 );
191 ok &= NearEqual(6.*r[3], exp(x/2.)/16., eps99 , eps99 );
192
193 return ok;
194
195 }
196
197 } // END empty namespace
198
Sqrt(void)199 bool Sqrt(void)
200 { bool ok = true;
201 ok &= SqrtTestOne();
202 ok &= SqrtTestTwo();
203 ok &= SqrtTestThree();
204 return ok;
205 }
206