1 // Copyright Matthew Pulver 2018 - 2019.
2 // Distributed under the Boost Software License, Version 1.0.
3 // (See accompanying file LICENSE_1_0.txt or copy at
4 // https://www.boost.org/LICENSE_1_0.txt)
5
6 #ifndef BOOST_MATH_TEST_AUTODIFF_HPP
7 #define BOOST_MATH_TEST_AUTODIFF_HPP
8
9 #ifndef BOOST_TEST_MODULE
10 #define BOOST_TEST_MODULE test_autodiff
11 #endif
12
13 #ifndef BOOST_ALLOW_DEPRECATED_HEADERS
14 #define BOOST_ALLOW_DEPRECATED_HEADERS // artifact of sp_typeinfo.hpp inclusion from unit_test.hpp
15 #endif
16
17 #include <boost/math/tools/config.hpp>
18
19 #include <boost/math/differentiation/autodiff.hpp>
20 #include <boost/multiprecision/cpp_bin_float.hpp>
21 #include <boost/multiprecision/cpp_dec_float.hpp>
22 #include <boost/mp11/function.hpp>
23 #include <boost/mp11/integral.hpp>
24 #include <boost/mp11/list.hpp>
25 #include <boost/mp11/utility.hpp>
26 #include <boost/range/irange.hpp>
27 #include <boost/test/included/unit_test.hpp>
28
29 #include <algorithm>
30 #include <cfenv>
31 #include <cstdlib>
32 #include <random>
33
34 namespace mp11 = boost::mp11;
35 namespace bmp = boost::multiprecision;
36
37 #if defined(BOOST_USE_VALGRIND) || defined(BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS)
38 using bin_float_types = mp11::mp_list<float>;
39 #else
40 using bin_float_types = mp11::mp_list<float, double, long double>;
41 #endif
42
43 // cpp_dec_float_50 cannot be used with close_at_tolerance
44 /*using multiprecision_float_types =
45 mp_list<bmp::cpp_dec_float_50, bmp::cpp_bin_float_50>;*/
46
47 #if !defined(BOOST_VERSION) || BOOST_VERSION < 107000 || defined(BOOST_USE_VALGRIND) || defined(BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS) || defined(BOOST_NO_STRESS_TEST)
48 using multiprecision_float_types = mp11::mp_list<>;
49 #else
50 #define BOOST_AUTODIFF_TESTING_INCLUDE_MULTIPRECISION
51 using multiprecision_float_types = mp11::mp_list<bmp::cpp_bin_float_50>;
52 #endif
53
54 using all_float_types = mp11::mp_append<bin_float_types, multiprecision_float_types>;
55
56 using namespace boost::math::differentiation;
57
58 namespace test_detail {
59 template <typename T>
60 using is_multiprecision_t =
61 mp11::mp_or<bmp::is_number<T>, bmp::is_number_expression<T>>;
62
63 template<bool IfValue, typename ThenType, typename ElseType>
64 using if_c = mp11::mp_eval_if_c<IfValue, ThenType, mp11::mp_identity_t, ElseType>;
65
66 template<typename IfType, typename ThenType, typename ElseType>
67 using if_t = if_c<IfType::value, ThenType, ElseType>;
68
69 /**
70 * Simple struct to hold constants that are used in each test
71 * since BOOST_AUTO_TEST_CASE_TEMPLATE doesn't support fixtures.
72 */
73 template <typename T, std::size_t OrderValue>
74 struct test_constants_t {
75 static constexpr auto n_samples = if_t<mp11::mp_or<bmp::is_number<T>, bmp::is_number_expression<T>>, mp11::mp_int<10>, mp11::mp_int<25>>::value;
76 static constexpr auto order = OrderValue;
pct_epsilontest_detail::test_constants_t77 static constexpr T pct_epsilon() BOOST_NOEXCEPT {
78 return (is_multiprecision_t<T>::value ? 2 : 1) * std::numeric_limits<T>::epsilon() * 100;
79 }
80 };
81
82 /**
83 * struct to emit pseudo-random values from a given interval.
84 * Endpoints are closed or open depending on whether or not they're infinite).
85 */
86
87 template <typename T>
88 struct RandomSample {
89 using numeric_limits_t = std::numeric_limits<T>;
90 using is_integer_t = mp11::mp_bool<std::numeric_limits<T>::is_integer>;
91
92 using distribution_param_t = if_t<
93 is_multiprecision_t<T>,
94 if_t<is_integer_t,
95 if_c<numeric_limits_t::is_signed, int64_t, uint64_t>,
96 long double>,
97 T>;
98 static_assert((std::numeric_limits<T>::is_integer &&
99 std::numeric_limits<distribution_param_t>::is_integer) ||
100 (!std::numeric_limits<T>::is_integer &&
101 !std::numeric_limits<distribution_param_t>::is_integer),
102 "T and distribution_param_t must either both be integral or "
103 "both be not integral");
104
105 using dist_t = if_t<is_integer_t,
106 std::uniform_int_distribution<distribution_param_t>,
107 std::uniform_real_distribution<distribution_param_t>>;
108
109 struct get_integral_endpoint {
110 template <typename V>
operator ()test_detail::RandomSample::get_integral_endpoint111 constexpr distribution_param_t operator()(V finish) const noexcept {
112 return static_cast<distribution_param_t>(finish);
113 }
114 };
115
116 struct get_real_endpoint {
117 template <typename V>
operator ()test_detail::RandomSample::get_real_endpoint118 constexpr distribution_param_t operator()(V finish) const noexcept {
119 return std::nextafter(static_cast<distribution_param_t>(finish),
120 (std::numeric_limits<distribution_param_t>::max)());
121 }
122 };
123
124 using get_endpoint_t = if_t<is_integer_t, get_integral_endpoint, get_real_endpoint>;
125
126 template <typename U, typename V>
RandomSampletest_detail::RandomSample127 RandomSample(U start, V finish)
128 : rng_(std::random_device{}()),
129 dist_(static_cast<distribution_param_t>(start),
130 get_endpoint_t{}(finish)) {}
131
nexttest_detail::RandomSample132 T next() noexcept { return static_cast<T>(dist_(rng_)); }
normalizetest_detail::RandomSample133 T normalize(const T& x) noexcept {
134 return x / ((dist_.max)() - (dist_.min)());
135 }
136
137 std::mt19937 rng_;
138 dist_t dist_;
139 };
140 static_assert(std::is_same<RandomSample<float>::dist_t,
141 std::uniform_real_distribution<float>>::value,
142 "");
143 static_assert(std::is_same<RandomSample<int64_t>::dist_t,
144 std::uniform_int_distribution<int64_t>>::value,
145 "");
146 static_assert(std::is_same<RandomSample<bmp::uint512_t>::dist_t,
147 std::uniform_int_distribution<uint64_t>>::value,
148 "");
149 static_assert(std::is_same<RandomSample<bmp::cpp_bin_float_50>::dist_t,
150 std::uniform_real_distribution<long double>>::value,
151 "");
152
153 } // namespace test_detail
154
155 template<typename T>
isNearZero(const T & t)156 auto isNearZero(const T& t) noexcept -> typename std::enable_if<!detail::is_fvar<T>::value, bool>::type
157 {
158 using std::sqrt;
159 using bmp::sqrt;
160 using detail::sqrt;
161 using std::fabs;
162 using bmp::fabs;
163 using detail::fabs;
164 using boost::math::fpclassify;
165 using std::sqrt;
166 return fpclassify(fabs(t)) == FP_ZERO || fpclassify(fabs(t)) == FP_SUBNORMAL || boost::math::fpc::is_small(fabs(t), sqrt(std::numeric_limits<T>::epsilon()));
167 }
168
169 template<typename T>
isNearZero(const T & t)170 auto isNearZero(const T& t) noexcept -> typename std::enable_if<detail::is_fvar<T>::value, bool>::type
171 {
172 using root_type = typename T::root_type;
173 return isNearZero(static_cast<root_type>(t));
174 }
175
176 template <typename T, std::size_t Order = 5>
177 using test_constants_t = test_detail::test_constants_t<T, Order>;
178
179 template <typename W, typename X, typename Y, typename Z>
mixed_partials_f(const W & w,const X & x,const Y & y,const Z & z)180 promote<W, X, Y, Z> mixed_partials_f(const W& w, const X& x, const Y& y,
181 const Z& z) {
182
183 return exp(w * sin(x * log(y) / z) + sqrt(w * z / (x * y))) + w * w / tan(z);
184 }
185
186 // Equations and function/variable names are from
187 // https://en.wikipedia.org/wiki/Greeks_(finance)#Formulas_for_European_option_Greeks
188 //
189 // Standard normal probability density function
190 template <typename T>
phi(const T & x)191 T phi(const T& x) {
192 return boost::math::constants::one_div_root_two_pi<T>() * exp(-0.5 * x * x);
193 }
194
195 // Standard normal cumulative distribution function
196 template <typename T>
Phi(const T & x)197 T Phi(const T& x) {
198 return 0.5 * erfc(-boost::math::constants::one_div_root_two<T>() * x);
199 }
200
201 enum class CP { call, put };
202
203 // Assume zero annual dividend yield (q=0).
204 template <typename Price, typename Sigma, typename Tau, typename Rate>
black_scholes_option_price(CP cp,double K,const Price & S,const Sigma & sigma,const Tau & tau,const Rate & r)205 promote<Price, Sigma, Tau, Rate> black_scholes_option_price(CP cp, double K,
206 const Price& S,
207 const Sigma& sigma,
208 const Tau& tau,
209 const Rate& r) {
210 const auto d1 =
211 (log(S / K) + (r + sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
212 const auto d2 =
213 (log(S / K) + (r - sigma * sigma / 2) * tau) / (sigma * sqrt(tau));
214 if (cp == CP::call) {
215 return S * Phi(d1) - exp(-r * tau) * K * Phi(d2);
216 }
217 return exp(-r * tau) * K * Phi(-d2) - S * Phi(-d1);
218 }
219
220 template <typename T>
uncast_return(const T & x)221 T uncast_return(const T& x) {
222 return x == 0 ? 0 : 1;
223 }
224
225 #endif // BOOST_MATH_TEST_AUTODIFF_HPP
226