1 /* --------------------------------------------------------------------------
2 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-20 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 $begin zdouble.cpp$$
14 $spell
15 zdouble
16 $$
17
18 $section zdouble: Example and Test$$
19
20 $srcthisfile%0%// BEGIN C++%// END C++%1%$$
21
22 $end
23 */
24 // BEGIN C++
25 # include <cppad/cppad.hpp>
26
27 namespace {
test_one(void)28 template <class Base> bool test_one(void)
29 { bool ok = true;
30 Base eps99 = 99. * std::numeric_limits<double>::epsilon();
31
32 typedef CppAD::AD<Base> a1type;
33 typedef CppAD::AD<a1type> a2type;
34
35 // value during taping
36 size_t n = 2;
37 CPPAD_TESTVECTOR(Base) x(n);
38 x[0] = 0.0;
39 x[1] = 0.0;
40
41 // declare independent variable
42 CPPAD_TESTVECTOR(a2type) a2x(n);
43 for (size_t j = 0; j < n; j++)
44 a2x[j] = a2type( a1type(x[j]) );
45 Independent(a2x);
46
47 // zero and one as a2type values
48 a2type a2zero = a2type(0.0);
49 a2type a2one = a2type(1.0);
50
51 // h(x) = x[0] / x[1] if x[1] > x[0] else 1.0
52 a2type h_x = CondExpGt(a2x[1], a2x[0], a2x[0] / a2x[1], a2one);
53
54 // f(x) = h(x) if x[0] > 0.0 else 0.0
55 // = x[0] / x[1] if x[1] > x[0] and x[0] > 0.0
56 // = 1.0 if x[0] >= x[1] and x[0] > 0.0
57 // = 0.0 if x[0] <= 0.0
58 a2type f_x = CondExpGt(a2x[0], a2zero, h_x, a2one);
59
60 // define the function f(x)
61 size_t m = 1;
62 CPPAD_TESTVECTOR(a2type) a2y(m);
63 a2y[0] = f_x;
64 CppAD::ADFun<a1type> af1;
65 af1.Dependent(a2x, a2y);
66
67 // Define function g(x) = gradient of f(x)
68 CPPAD_TESTVECTOR(a1type) a1x(n), a1z(n), a1w(m);
69 for (size_t j = 0; j < n; j++)
70 a1x[j] = a1type(x[j]);
71 a1w[0] = a1type(1.0);
72 Independent(a1x);
73 af1.Forward(0, a1x);
74 a1z = af1.Reverse(1, a1w);
75 CppAD::ADFun<Base> g;
76 g.Dependent(a1x, a1z);
77
78 // check result for a case where f(x) = 0.0;
79 CPPAD_TESTVECTOR(Base) z(2);
80 x[0] = 0.0;
81 x[1] = 0.0;
82 z = g.Forward(0, x);
83 ok &= z[0] == 0.0;
84 ok &= z[1] == 0.0;
85
86 // check result for a case where f(x) = 1.0;
87 x[0] = 1.0;
88 x[1] = 0.5;
89 z = g.Forward(0, x);
90 ok &= z[0] == 0.0;
91 ok &= z[1] == 0.0;
92
93 // check result for a case where f(x) = x[0] / x[1];
94 x[0] = 1.0;
95 x[1] = 2.0;
96 z = g.Forward(0, x);
97 ok &= CppAD::NearEqual(z[0], 1.0/x[1], eps99, eps99);
98 ok &= CppAD::NearEqual(z[1], - x[0]/(x[1]*x[1]), eps99, eps99);
99
100 return ok;
101 }
test_two(void)102 bool test_two(void)
103 { bool ok = true;
104 using CppAD::zdouble;
105 //
106 zdouble eps = CppAD::numeric_limits<zdouble>::epsilon();
107 ok &= eps == std::numeric_limits<double>::epsilon();
108 //
109 zdouble min = CppAD::numeric_limits<zdouble>::min();
110 ok &= min == std::numeric_limits<double>::min();
111 //
112 zdouble max = CppAD::numeric_limits<zdouble>::max();
113 ok &= max == std::numeric_limits<double>::max();
114 //
115 zdouble nan = CppAD::numeric_limits<zdouble>::quiet_NaN();
116 ok &= nan != nan;
117 //
118 int digits10 = CppAD::numeric_limits<zdouble>::digits10;
119 ok &= digits10 == std::numeric_limits<double>::digits10;
120 //
121 return ok;
122 }
123 }
124
zdouble(void)125 bool zdouble(void)
126 { bool ok = true;
127 using CppAD::AD;
128 using CppAD::NearEqual;
129 using CppAD::zdouble;
130 //
131 ok &= test_one<zdouble>();
132 ok &= test_one<double>();
133 //
134 ok &= test_two();
135 //
136 return ok;
137 }
138 // END C++
139