libs/math/test/test_autodiff.hpp

   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;
  77   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>
 111     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>
 118     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>
 127   RandomSample(U start, V finish)
 128       : rng_(std::random_device{}()),
 129         dist_(static_cast<distribution_param_t>(start),
 130               get_endpoint_t{}(finish)) {}
 131
 132   T next() noexcept { return static_cast<T>(dist_(rng_)); }
 133   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>
 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>
 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>
 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>
 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>
 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>
 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>
 221 T uncast_return(const T& x) {
 222   return x == 0 ? 0 : 1;
 223 }
 224
 225 #endif  // BOOST_MATH_TEST_AUTODIFF_HPP