diff --git a/example/autodiff_reverse_black_scholes.cpp b/example/autodiff_reverse_black_scholes.cpp new file mode 100644 index 000000000..e02069a47 --- /dev/null +++ b/example/autodiff_reverse_black_scholes.cpp @@ -0,0 +1,149 @@ +#include + +using namespace boost::math::differentiation::reverse_mode; +using namespace boost::math::constants; + +template +X phi(X const& x) +{ + return one_div_root_two_pi() * exp(-0.5 * x * x); +} + +template +X Phi(X const& x) +{ + return 0.5 * erfc(-one_div_root_two() * x); +} + +enum class CP { call, put }; + +template +T black_scholes_option_price(CP cp, double K, T const& S, T const& sigma, T const& tau, T const& r) +{ + using namespace std; + auto const d1 = (log(S / K) + (r + sigma * sigma / 2) * tau) / (sigma * sqrt(tau)); + auto const d2 = (log(S / K) + (r - sigma * sigma / 2) * tau) / (sigma * sqrt(tau)); + switch (cp) { + case CP::call: + return S * Phi(d1) - exp(-r * tau) * K * Phi(d2); + case CP::put: + return exp(-r * tau) * K * Phi(-d2) - S * Phi(-d1); + default: + throw std::runtime_error("Invalid CP value."); + } +} + +int main() +{ + double const K = 100.0; + double S_val = 105.0; + double sigma_val = 5.0; + double tau_val = 30.0 / 365; + double r_val = 1.25 / 100; + rvar S = make_rvar(S_val); + rvar sigma = make_rvar(sigma_val); + rvar tau = make_rvar(tau_val); + rvar r = make_rvar(r_val); + + rvar call_price + = black_scholes_option_price>(CP::call, K, S, sigma, tau, r); + rvar put_price + = black_scholes_option_price>(CP::put, K, S, sigma, tau, r); + + double const d1 = ((log(S_val / K) + (r_val + sigma_val * sigma_val / 2) * tau_val) + / (sigma_val * sqrt(tau_val))); + double const d2 = ((log(S_val / K) + (r_val - sigma_val * sigma_val / 2) * tau_val) + / (sigma_val * sqrt(tau_val))); + double const formula_call_delta = +Phi(+d1); + double const formula_put_delta = -Phi(-d1); + double const formula_vega = (S_val * phi(d1) * sqrt(tau_val)); + double const formula_call_theta = (-S_val * phi(d1) * sigma_val / (2 * sqrt(tau_val)) + - r_val * K * exp(-r_val * tau_val) * Phi(+d2)); + + double const formula_put_theta = (-S_val * phi(d1) * sigma_val / (2 * sqrt(tau_val)) + + r_val * K * exp(-r_val * tau_val) * Phi(-d2)); + double const formula_call_rho = (+K * tau_val * exp(-r_val * tau_val) * Phi(+d2)); + double const formula_put_rho = (-K * tau_val * exp(-r_val * tau_val) * Phi(-d2)); + double const formula_gamma = (phi(d1) / (S_val * sigma_val * sqrt(tau_val))); + double const formula_vanna = (-phi(d1) * d2 / sigma_val); + double const formula_charm = (phi(d1) * (d2 * sigma_val * sqrt(tau_val) - 2 * r_val * tau_val) + / (2 * tau_val * sigma_val * sqrt(tau_val))); + double const formula_vomma = (S_val * phi(d1) * sqrt(tau_val) * d1 * d2 / sigma_val); + double const formula_veta = (-S_val * phi(d1) * sqrt(tau_val) + * (r_val * d1 / (sigma_val * sqrt(tau_val)) + - (1 + d1 * d2) / (2 * tau_val))); + double const formula_speed = (-phi(d1) * (d1 / (sigma_val * sqrt(tau_val)) + 1) + / (S_val * S_val * sigma_val * sqrt(tau_val))); + double const formula_zomma = (phi(d1) * (d1 * d2 - 1) + / (S_val * sigma_val * sigma_val * sqrt(tau_val))); + double const formula_color = (-phi(d1) / (2 * S_val * tau_val * sigma_val * sqrt(tau_val)) + * (1 + + (2 * r_val * tau_val - d2 * sigma_val * sqrt(tau_val)) * d1 + / (sigma_val * sqrt(tau_val)))); + double const formula_ultima = -formula_vega + * ((d1 * d2 * (1 - d1 * d2) + d1 * d1 + d2 * d2) + / (sigma_val * sigma_val)); + + auto call_greeks = grad(call_price, &S, &sigma, &tau, &r); + auto put_greeks = grad(put_price, &S, &sigma, &tau, &r); + + auto call_greeks_2nd_order = grad_nd<2>(call_price, &S, &sigma, &tau, &r); + auto put_greeks_2nd_order = grad_nd<2>(put_price, &S, &sigma, &tau, &r); + + auto call_greeks_3rd_order = grad_nd<3>(call_price, &S, &sigma, &tau, &r); + auto put_greeks_3rd_order = grad_nd<3>(put_price, &S, &sigma, &tau, &r); + + std::cout << std::setprecision(std::numeric_limits::digits10) + << "autodiff black-scholes call price = " << call_price.item() << "\n" + << "autodiff black-scholes put price = " << put_price.item() << "\n" + << "\n ## First-order Greeks \n" + << "autodiff call delta = " << call_greeks[0].item() << "\n" + << "formula call delta = " << formula_call_delta << "\n" + << "autodiff call vega = " << call_greeks[1].item() << '\n' + << " formula call vega = " << formula_vega << '\n' + << "autodiff call theta = " << call_greeks[2].item() << '\n' + << " formula call theta = " << formula_call_theta << '\n' + << "autodiff call rho = " << call_greeks[3].item() << 'n' + << " formula call rho = " << formula_call_rho << '\n' + << '\n' + << "autodiff put delta = " << put_greeks[0].item() << 'n' + << " formula put delta = " << formula_put_delta << '\n' + << "autodiff put vega = " << put_greeks[1].item() << '\n' + << " formula put vega = " << formula_vega << '\n' + << "autodiff put theta = " << -put_greeks[2].item() << '\n' + << " formula put theta = " << formula_put_theta << '\n' + << "autodiff put rho = " << put_greeks[3].item() << '\n' + << " formula put rho = " << formula_put_rho << '\n' + + << "\n## Second-order Greeks\n" + << "autodiff call gamma = " << call_greeks_2nd_order[0][0].item() << '\n' + << "autodiff put gamma = " << put_greeks_2nd_order[0][0].item() << '\n' + << " formula gamma = " << formula_gamma << '\n' + << "autodiff call vanna = " << call_greeks_2nd_order[0][1].item() << '\n' + << "autodiff put vanna = " << put_greeks_2nd_order[0][1].item() << '\n' + << " formula vanna = " << formula_vanna << '\n' + << "autodiff call charm = " << -call_greeks_2nd_order[0][2].item() << '\n' + << "autodiff put charm = " << -put_greeks_2nd_order[0][2].item() << '\n' + << " formula charm = " << formula_charm << '\n' + << "autodiff call vomma = " << call_greeks_2nd_order[1][1].item() << '\n' + << "autodiff put vomma = " << put_greeks_2nd_order[1][1].item() << '\n' + << " formula vomma = " << formula_vomma << '\n' + << "autodiff call veta = " << call_greeks_2nd_order[1][2].item() << '\n' + << "autodiff put veta = " << put_greeks_2nd_order[1][2].item() << '\n' + << " formula veta = " << formula_veta << '\n' + + << "\n## Third-order Greeks\n" + << "autodiff call speed = " << call_greeks_3rd_order[0][0][0] << '\n' + << "autodiff put speed = " << put_greeks_3rd_order[0][0][0] << '\n' + << " formula speed = " << formula_speed << '\n' + << "autodiff call zomma = " << call_greeks_3rd_order[0][0][1] << '\n' + << "autodiff put zomma = " << put_greeks_3rd_order[0][0][1] << '\n' + << " formula zomma = " << formula_zomma << '\n' + << "autodiff call color = " << call_greeks_3rd_order[0][0][2] << '\n' + << "autodiff put color = " << put_greeks_3rd_order[0][0][2] << '\n' + << " formula color = " << formula_color << '\n' + << "autodiff call ultima = " << call_greeks_3rd_order[1][1][1] << '\n' + << "autodiff put ultima = " << put_greeks_3rd_order[1][1][1] << '\n' + << " formula ultima = " << formula_ultima << '\n'; + return 0; +} diff --git a/include/boost/math/differentiation/CMakeLists.txt b/include/boost/math/differentiation/CMakeLists.txt index 4d704762c..5a381a59e 100644 --- a/include/boost/math/differentiation/CMakeLists.txt +++ b/include/boost/math/differentiation/CMakeLists.txt @@ -27,36 +27,38 @@ set(HEADERS ${CMAKE_SOURCE_DIR}/detail/reverse_mode_autodiff_expression_template_base.hpp ${CMAKE_SOURCE_DIR}/detail/reverse_mode_autodiff_utilities.hpp ${CMAKE_SOURCE_DIR}/detail/reverse_mode_autodiff_erf_overloads.hpp + ${CMAKE_SOURCE_DIR}/detail/reverse_mode_autodiff_comparison_operator_overloads.hpp ${BOOST_MATH_TEST_DIRECTORY}test_autodiff_reverse.hpp ) set(TEST_MEMORY_ALLOCATORS - ${BOOST_MATH_TEST_DIRECTORY}test_autodiff_reverse_flat_linear_allocator.cpp + ${BOOST_MATH_TEST_DIRECTORY}test_reverse_mode_autodiff_flat_linear_allocator.cpp ) set(TEST_CONSTRUCTORS - ${BOOST_MATH_TEST_DIRECTORY}test_rvar_constructors.cpp + ${BOOST_MATH_TEST_DIRECTORY}test_reverse_mode_autodiff_constructors.cpp ) set(TEST_COMPARISON_OPS - ${BOOST_MATH_TEST_DIRECTORY}test_comparison_operators.cpp + ${BOOST_MATH_TEST_DIRECTORY}test_reverse_mode_autodiff_comparison_operators.cpp ) set(TEST_STL_SUPPORT - ${BOOST_MATH_TEST_DIRECTORY}test_stl_support.cpp + ${BOOST_MATH_TEST_DIRECTORY}test_reverse_mode_autodiff_stl_support.cpp ) -SET(TEST_BASIC_OPS ${BOOST_MATH_TEST_DIRECTORY}test_autodiff_basic_math_ops.cpp) +SET(TEST_BASIC_OPS ${BOOST_MATH_TEST_DIRECTORY}test_reverse_mode_autodiff_basic_math_ops.cpp) SET(TEST_ADDITIONAL_BOOST_FUNCTIONS ${BOOST_MATH_TEST_DIRECTORY}test_autodiff_reverse_additional_boost_functions.cpp) SET(LINEAR_REGRESSION_EXAMPLE ${BOOST_MATH_EXAMPLE_DIRECTORY}/reverse_mode_linear_regression_example.cpp) set(TEST_ERF ${BOOST_MATH_TEST_DIRECTORY}/test_reverse_mode_autodiff_error_functions.cpp) +set(BLACK_SCHOLES_EXAMPLE ${BOOST_MATH_EXAMPLE_DIRECTORY}/autodiff_reverse_black_scholes.cpp) # add_executable(autodiff_reverse_tests ${TEST_CONSTRUCTORS} ${HEADERS}) # add_executable(autodiff_reverse_tests ${TEST_MEMORY_ALLOCATORS} ${HEADERS}) # add_executable(autodiff_reverse_tests ${TEST_BASIC_OPS} ${HEADERS}) -# add_executable(autodiff_reverse_tests ${TEST_COMPARISON_OPS} ${HEADERS} -# detail/reverse_mode_autodiff_comparison_operator_overloads.hpp) -# add_executable(autodiff_reverse_tests ${TEST_STL_SUPPORT} ${HEADERS}) +# add_executable(autodiff_reverse_tests ${TEST_COMPARISON_OPS} ${HEADERS}) +add_executable(autodiff_reverse_tests ${TEST_STL_SUPPORT} ${HEADERS}) # add_executable(autodiff_reverse_tests ${LINEAR_REGRESSION_EXAMPLE} ${HEADERS}) -add_executable(autodiff_reverse_tests ${TEST_ERF} ${HEADERS}) +# add_executable(autodiff_reverse_tests ${TEST_ERF} ${HEADERS}) +# add_executable(autodiff_reverse_tests ${BLACK_SCHOLES_EXAMPLE} ${HEADERS}) diff --git a/include/boost/math/differentiation/CMakeLists.txt.user b/include/boost/math/differentiation/CMakeLists.txt.user index e5f05af95..f7fab51fd 100644 --- a/include/boost/math/differentiation/CMakeLists.txt.user +++ b/include/boost/math/differentiation/CMakeLists.txt.user @@ -1,6 +1,6 @@ - + EnvironmentId diff --git a/include/boost/math/differentiation/autodiff_reverse.hpp b/include/boost/math/differentiation/autodiff_reverse.hpp index 22a436d7f..2f3cfb648 100644 --- a/include/boost/math/differentiation/autodiff_reverse.hpp +++ b/include/boost/math/differentiation/autodiff_reverse.hpp @@ -9,7 +9,6 @@ #include #include #include - #include #include #include @@ -567,6 +566,22 @@ namespace detail { template struct grad_op_impl { + std::vector> operator()(rvar &f, std::vector *> &x) + { + auto &tape = get_active_tape(); + tape.zero_grad(); + f.backward(); + + std::vector> gradient_vector; + gradient_vector.reserve(x.size()); + + for (auto xi : x) { + // make a new rvar holding the adjoint value + gradient_vector.emplace_back(xi->adjoint()); + } + return gradient_vector; + } + /* std::vector *> operator()(rvar &f, std::vector *> &x) { @@ -578,7 +593,7 @@ struct grad_op_impl gradient_vector.push_back(&(xi->adjoint())); } return gradient_vector; - } + }*/ }; /** @brief helper overload for grad implementation. * @return vector of gradients of the autodiff graph. @@ -607,6 +622,22 @@ struct grad_op_impl template struct grad_nd_impl { + auto operator()(rvar &f, std::vector *> &x) + { + static_assert(N > 1, "N must be greater than 1 for this template"); + + auto current_grad = grad(f, x); // vector> or vector + + std::vector()(current_grad[0], x))> + result; + result.reserve(current_grad.size()); + + for (auto &g_i : current_grad) { + result.push_back(grad_nd_impl()(g_i, x)); + } + return result; + } + /* auto operator()(rvar &f, std::vector *> &x) { static_assert(N > 1, "N must be greater than 1 for this template"); @@ -617,7 +648,7 @@ struct grad_nd_impl result.push_back(grad_nd_impl()(*g_i, x)); } return result; - } + }*/ }; /** @brief spcialization for order = 1, * @return vector> gradients */ diff --git a/include/boost/math/differentiation/detail/reverse_mode_autodiff_basic_operator_overloads.hpp b/include/boost/math/differentiation/detail/reverse_mode_autodiff_basic_operator_overloads.hpp index b534bcbf1..4245ebc03 100644 --- a/include/boost/math/differentiation/detail/reverse_mode_autodiff_basic_operator_overloads.hpp +++ b/include/boost/math/differentiation/detail/reverse_mode_autodiff_basic_operator_overloads.hpp @@ -219,7 +219,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> mult_const_expr operator*(const expression &arg, const U &v) { return mult_const_expr(arg, static_cast(v)); @@ -228,7 +228,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> mult_const_expr operator*(const U &v, const expression &arg) { return mult_const_expr(arg, static_cast(v)); @@ -245,7 +245,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> add_const_expr operator+(const expression &arg, const U &v) { return add_const_expr(arg, static_cast(v)); @@ -254,7 +254,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> add_const_expr operator+(const U &v, const expression &arg) { return add_const_expr(arg, static_cast(v)); @@ -284,7 +284,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> add_const_expr operator-(const expression &arg, const U &v) { /* rvar - float = rvar + (-float) */ @@ -299,7 +299,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> auto operator-(const U &v, const expression &arg) { auto neg = -arg; @@ -318,7 +318,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> const_div_by_expr operator/(const U &v, const expression &arg) { return const_div_by_expr(arg, static_cast(v)); @@ -328,7 +328,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> div_by_const_expr operator/(const expression &arg, const U &v) { return div_by_const_expr(arg, static_cast(v)); diff --git a/include/boost/math/differentiation/detail/reverse_mode_autodiff_comparison_operator_overloads.hpp b/include/boost/math/differentiation/detail/reverse_mode_autodiff_comparison_operator_overloads.hpp index addd8a2aa..a6059cb81 100644 --- a/include/boost/math/differentiation/detail/reverse_mode_autodiff_comparison_operator_overloads.hpp +++ b/include/boost/math/differentiation/detail/reverse_mode_autodiff_comparison_operator_overloads.hpp @@ -1,6 +1,5 @@ #ifndef REVERSE_MODE_AUTODIFF_COMPARISON_OPERATOR_OVERLOADS_HPP #define REVERSE_MODE_AUTODIFF_COMPARISON_OPERATOR_OVERLOADS_HPP - #include namespace boost { namespace math { @@ -16,7 +15,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator==(const expression &lhs, const U &rhs) { return lhs.evaluate() == static_cast(rhs); @@ -25,7 +24,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator==(const U &lhs, const expression &rhs) { return lhs == rhs.evaluate(); @@ -40,7 +39,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator!=(const expression &lhs, const U &rhs) { return lhs.evaluate() != rhs; @@ -49,7 +48,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator!=(const U &lhs, const expression &rhs) { return lhs != rhs.evaluate(); @@ -64,7 +63,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator<(const expression &lhs, const U &rhs) { return lhs.evaluate() < rhs; @@ -73,7 +72,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator<(const U &lhs, const expression &rhs) { return lhs < rhs.evaluate(); @@ -88,7 +87,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator>(const expression &lhs, const U &rhs) { return lhs.evaluate() > rhs; @@ -97,7 +96,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator>(const U &lhs, const expression &rhs) { return lhs > rhs.evaluate(); @@ -112,7 +111,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator<=(const expression &lhs, const U &rhs) { return lhs.evaluate() <= rhs; @@ -121,7 +120,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator<=(const U &lhs, const expression &rhs) { return lhs <= rhs.evaluate(); @@ -136,7 +135,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator>=(const expression &lhs, const U &rhs) { return lhs.evaluate() >= rhs; @@ -145,7 +144,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> bool operator>=(const U &lhs, const expression &rhs) { return lhs >= rhs.evaluate(); diff --git a/include/boost/math/differentiation/detail/reverse_mode_autodiff_erf_overloads.hpp b/include/boost/math/differentiation/detail/reverse_mode_autodiff_erf_overloads.hpp index b23e62ce3..e62434201 100644 --- a/include/boost/math/differentiation/detail/reverse_mode_autodiff_erf_overloads.hpp +++ b/include/boost/math/differentiation/detail/reverse_mode_autodiff_erf_overloads.hpp @@ -6,7 +6,7 @@ #include #include #include - +#include namespace boost { namespace math { namespace differentiation { @@ -72,8 +72,10 @@ struct erf_expr : public abstract_unary_expression()); + using std::sqrt; + return 2 * exp(-argv * argv) / sqrt(constants::pi()); } }; @@ -101,8 +103,10 @@ struct erfc_expr : public abstract_unary_expression()); + using std::sqrt; + return -2 * exp(-argv * argv) / sqrt(constants::pi()); } }; @@ -130,14 +134,19 @@ struct erf_inv_expr : public abstract_unary_expression 1)>( [](auto &&x) { - return 0.5 * std::sqrt(constants::pi()) + return 0.5 * sqrt(constants::pi()) * reverse_mode::exp(reverse_mode::pow(reverse_mode::erf_inv(x), 2.0)); }, [](auto &&x) { - return 0.5 * std::sqrt(constants::pi()) - * std::exp(std::pow(boost::math::erf_inv(x), 2)); + return 0.5 * sqrt(constants::pi()) * exp(pow(boost::math::erf_inv(x), 2)); }, argv); } @@ -167,14 +176,19 @@ struct erfc_inv_expr : public abstract_unary_expression 1)>( [](auto &&x) { - return -0.5 * std::sqrt(constants::pi()) + return -0.5 * sqrt(constants::pi()) * reverse_mode::exp(reverse_mode::pow(reverse_mode::erfc_inv(x), 2.0)); }, [](auto &&x) { - return -0.5 * std::sqrt(constants::pi()) - * std::exp(std::pow(boost::math::erfc_inv(x), 2)); + return -0.5 * sqrt(constants::pi()) * exp(pow(boost::math::erfc_inv(x), 2)); }, argv); } diff --git a/include/boost/math/differentiation/detail/reverse_mode_autodiff_expression_template_base.hpp b/include/boost/math/differentiation/detail/reverse_mode_autodiff_expression_template_base.hpp index 1177b5c0a..cc719547d 100644 --- a/include/boost/math/differentiation/detail/reverse_mode_autodiff_expression_template_base.hpp +++ b/include/boost/math/differentiation/detail/reverse_mode_autodiff_expression_template_base.hpp @@ -8,6 +8,9 @@ namespace differentiation { namespace reverse_mode { /* forward declarations for utitlity functions */ +struct expression_base +{}; + template struct expression; @@ -82,14 +85,7 @@ template constexpr std::size_t count_rvars = detail::count_rvar_impl::value; template -struct is_expression : std::false_type -{}; - -template -struct is_expression> : std::true_type -{}; -template -struct is_expression> : std::true_type +struct is_expression : std::is_base_of::type> {}; template @@ -109,9 +105,6 @@ struct rvar_type_impl template using rvar_t = typename detail::rvar_type_impl::type; -struct expression_base -{}; - template struct expression : expression_base { diff --git a/include/boost/math/differentiation/detail/reverse_mode_autodiff_memory_management.hpp b/include/boost/math/differentiation/detail/reverse_mode_autodiff_memory_management.hpp index 3f8bf3e0b..65c82110e 100644 --- a/include/boost/math/differentiation/detail/reverse_mode_autodiff_memory_management.hpp +++ b/include/boost/math/differentiation/detail/reverse_mode_autodiff_memory_management.hpp @@ -32,10 +32,10 @@ public: using difference_type = ptrdiff_t; private: - size_t begin_ = 0; - size_t index_ = 0; - size_t end_ = 0; const allocator_type *storage_ = nullptr; + size_t index_ = 0; + size_t begin_ = 0; + size_t end_ = 0; public: flat_linear_allocator_iterator() = default; diff --git a/include/boost/math/differentiation/detail/reverse_mode_autodiff_stl_overloads.hpp b/include/boost/math/differentiation/detail/reverse_mode_autodiff_stl_overloads.hpp index 3714e9eb3..44231adc2 100644 --- a/include/boost/math/differentiation/detail/reverse_mode_autodiff_stl_overloads.hpp +++ b/include/boost/math/differentiation/detail/reverse_mode_autodiff_stl_overloads.hpp @@ -119,6 +119,7 @@ struct trunc_expr : public abstract_unary_expressionarg.evaluate()); } @@ -718,6 +719,89 @@ struct atanh_expr : public abstract_unary_expression +struct fmod_expr + : public abstract_binary_expression> +{ + /** @brief + * */ + using lhs_type = LHS; + using rhs_type = RHS; + using value_type = T; + using inner_t = rvar_t; + // Explicitly define constructor to forward to base class + explicit fmod_expr(const expression &left_hand_expr, + const expression &right_hand_expr) + : abstract_binary_expression>( + left_hand_expr, right_hand_expr) + {} + + inner_t evaluate() const + { + using std::fmod; + return fmod(this->lhs.evaluate(), this->rhs.evaluate()); + }; + static const inner_t left_derivative(const inner_t &l, const inner_t &r, const inner_t &v) + { + return inner_t{1.0}; + }; + static const inner_t right_derivative(const inner_t &l, const inner_t &r, const inner_t &v) + { + using std::trunc; + return -1.0 * trunc(l / r); + }; +}; + +template +struct fmod_left_float_expr + : public abstract_unary_expression> +{ + /** @brief + * */ + using arg_type = ARG; + using value_type = T; + using inner_t = rvar_t; + + explicit fmod_left_float_expr(const expression &arg_expr, const T v) + : abstract_unary_expression>(arg_expr, + v){}; + + inner_t evaluate() const + { + using std::fmod; + return fmod(this->constant, this->arg.evaluate()); + } + static const inner_t derivative(const inner_t &argv, const inner_t &v, const T &constant) + { + return -1.0 * trunc(constant / argv); + } +}; + +template +struct fmod_right_float_expr + : public abstract_unary_expression> +{ + /** @brief + * */ + using arg_type = ARG; + using value_type = T; + using inner_t = rvar_t; + + explicit fmod_right_float_expr(const expression &arg_expr, const T v) + : abstract_unary_expression>(arg_expr, + v){}; + + inner_t evaluate() const + { + using std::fmod; + return fmod(this->arg.evaluate(), this->constant); + } + static const inner_t derivative(const inner_t &argv, const inner_t &v, const T &constant) + { + return inner_t{1.0}; + } +}; +/**************************************************************************************************/ template fabs_expr fabs(const expression &arg) { @@ -756,7 +840,7 @@ template::value>::type> + typename = typename std::enable_if::value>::type> expr_pow_float_expr pow(const expression &arg, const U &v) { return expr_pow_float_expr(arg, static_cast(v)); @@ -785,7 +869,7 @@ auto frexp(const expression &arg, int *i) { using std::frexp; using std::pow; - T tmp = frexp(arg.evaluate(), i); + frexp(arg.evaluate(), i); return arg / pow(2.0, *i); } @@ -861,19 +945,22 @@ trunc_expr trunc(const expression &arg) template auto fmod(const expression &lhs, const expression &rhs) { - return lhs - trunc(lhs / rhs) * rhs; + //return lhs - trunc(lhs / rhs) * rhs; + return fmod_expr(lhs, rhs); } template auto fmod(const expression &lhs, const T rhs) { - return lhs - trunc(lhs / rhs) * rhs; + //return lhs - trunc(lhs / rhs) * rhs; + return fmod_right_float_expr(lhs, rhs); } template auto fmod(const T lhs, const expression &rhs) { - return lhs - trunc(lhs / rhs) * rhs; + //return lhs - trunc(lhs / rhs) * rhs; + return fmod_left_float_expr(rhs, lhs); } template diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index a666c7f06..6e9b6cf70 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -1382,11 +1382,19 @@ test-suite long-running-tests : [ run test_pFq_precision.cpp ../tools//mpfr ../tools//gmp /boost/test//boost_unit_test_framework /boost/system//boost_system /boost/chrono//boost_chrono : : : [ check-target-builds ../config//has_mpfr : : no ] [ requires cxx11_hdr_initializer_list cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_smart_ptr ] release clang:-Wno-literal-range ] [ run test_constant_generate.cpp : : : release USE_CPP_FLOAT=1 off:no ] ; - build-project ../example ; # Expect policy_ref_snips13 to fail (message about no Cauchy Mean). +test-suite "test_reverse_mode_autodiff" + : + [ run test_reverse_mode_autodiff_flat_linear_allocator.cpp ] + [ run test_reverse_mode_autodiff_constructors.cpp ] + [ run test_reverse_mode_autodiff_comparison_operators.cpp ] + [ run test_reverse_mode_autodiff_stl_support.cpp ] + [ run test_reverse_mode_autodiff_basic_math_ops.cpp ] + [ run test_reverse_mode_autodiff_error_functions.cpp ] + ; rule get_float128_tests { local result ; diff --git a/test/test_autodiff_basic_math_ops.cpp b/test/test_reverse_mode_autodiff_basic_math_ops.cpp similarity index 98% rename from test/test_autodiff_basic_math_ops.cpp rename to test/test_reverse_mode_autodiff_basic_math_ops.cpp index 75470439a..90bcf082f 100644 --- a/test/test_autodiff_basic_math_ops.cpp +++ b/test/test_reverse_mode_autodiff_basic_math_ops.cpp @@ -150,9 +150,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(division, T, all_float_types) template T f(T x, T y) { - //auto z1 = x / (y + x); - //auto z2 = y / (x - y); - //T f = z1 * z2; return (x * y) / ((x + y) * (x - y)); } @@ -293,7 +290,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(first_derivative, T, all_float_types) rvar x_ad = x; rvar y_ad = y; - T fv = f(x, y); + //T fv = f(x, y); rvar f_ad = f(x_ad, y_ad); auto grad_f_analytical = grad_f_a(x, y); /* intended use case */ @@ -323,7 +320,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(second_derivative_and_hessian, T, all_float_types) rvar x_ad = x; rvar y_ad = y; - T fv = f(x, y); + //T fv = f(x, y); rvar f_ad = f(x_ad, y_ad); gradient_tape& tape = get_active_tape(); tape.zero_grad(); @@ -349,7 +346,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(order_3_der, T, all_float_types) rvar x_ad = x; rvar y_ad = y; - T fv = f(x, y); + //T fv = f(x, y); rvar f_ad = f(x_ad, y_ad); gradient_tape& tape = get_active_tape(); tape.zero_grad(); diff --git a/test/test_comparison_operators.cpp b/test/test_reverse_mode_autodiff_comparison_operators.cpp similarity index 99% rename from test/test_comparison_operators.cpp rename to test/test_reverse_mode_autodiff_comparison_operators.cpp index 022ca5f1e..85b869003 100644 --- a/test/test_comparison_operators.cpp +++ b/test/test_reverse_mode_autodiff_comparison_operators.cpp @@ -171,7 +171,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(greater_than_rvar_and_scalar, T, all_float_types) T a = T(3.14); T b = T(3.14); T c = T(2.71); - T d = T(4.0); rvar var1 = rvar(a); rvar var2 = rvar(a); @@ -271,7 +270,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(greater_than_or_equal_rvar_and_scalar, T, all_floa T a = T(3.14); T b = T(3.14); T c = T(2.71); - T d = T(4.0); rvar var1 = rvar(a); rvar var2 = rvar(a); diff --git a/test/test_rvar_constructors.cpp b/test/test_reverse_mode_autodiff_constructors.cpp similarity index 99% rename from test/test_rvar_constructors.cpp rename to test/test_reverse_mode_autodiff_constructors.cpp index a2e91890a..e5e2a73b7 100644 --- a/test/test_rvar_constructors.cpp +++ b/test/test_reverse_mode_autodiff_constructors.cpp @@ -16,7 +16,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(rvar_constructors_and_utils, T, all_float_types) BOOST_CHECK_EQUAL(x_value, x1.evaluate().evaluate()); BOOST_CHECK_EQUAL(x_value, x2.evaluate().evaluate().evaluate()); - auto v = x1.item(); /* get item helper */ BOOST_CHECK_EQUAL(x_value, x0.item()); BOOST_CHECK_EQUAL(x_value, x1.item()); @@ -251,8 +250,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(inplace_division, T, all_float_types) BOOST_AUTO_TEST_CASE_TEMPLATE(test_rvar_ostream_output, T, all_float_types) { using namespace rdiff; - rvar x = 2.0; - rvar y = 3.0; + rvar x = T{2.0}; + rvar y = T{3.0}; rvar z = x * y; z.backward(); diff --git a/test/test_reverse_mode_autodiff_error_functions.cpp b/test/test_reverse_mode_autodiff_error_functions.cpp index a10824bb8..fa086a887 100644 --- a/test/test_reverse_mode_autodiff_error_functions.cpp +++ b/test/test_reverse_mode_autodiff_error_functions.cpp @@ -26,9 +26,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_erf, T, all_float_types) auto g4 = grad_nd<4>(y, &x); BOOST_CHECK_CLOSE(y.item(), answers[0], eps); - BOOST_CHECK_CLOSE(g1[0]->item(), answers[1], eps); - BOOST_CHECK_CLOSE(g2[0][0]->item(), answers[2], eps); - BOOST_CHECK_CLOSE(g3[0][0][0]->item(), answers[3], eps); + BOOST_CHECK_CLOSE(g1[0], answers[1], eps); + BOOST_CHECK_CLOSE(g2[0][0], answers[2], eps); + BOOST_CHECK_CLOSE(g3[0][0][0], answers[3], eps); BOOST_CHECK_CLOSE(g4[0][0][0][0], answers[4], eps); } @@ -52,9 +52,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_erfc, T, all_float_types) auto g4 = grad_nd<4>(y, &x); BOOST_CHECK_CLOSE(y.item(), answers[0], eps); - BOOST_CHECK_CLOSE(g1[0]->item(), answers[1], eps); - BOOST_CHECK_CLOSE(g2[0][0]->item(), answers[2], eps); - BOOST_CHECK_CLOSE(g3[0][0][0]->item(), answers[3], eps); + BOOST_CHECK_CLOSE(g1[0], answers[1], eps); + BOOST_CHECK_CLOSE(g2[0][0], answers[2], eps); + BOOST_CHECK_CLOSE(g3[0][0][0], answers[3], eps); BOOST_CHECK_CLOSE(g4[0][0][0][0], answers[4], eps); } @@ -78,9 +78,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_erf_inv, T, all_float_types) auto g4 = grad_nd<4>(y, &x); BOOST_CHECK_CLOSE(y.item(), answers[0], eps); - BOOST_CHECK_CLOSE(g1[0]->item(), answers[1], eps); - BOOST_CHECK_CLOSE(g2[0][0]->item(), answers[2], eps); - BOOST_CHECK_CLOSE(g3[0][0][0]->item(), answers[3], eps); + BOOST_CHECK_CLOSE(g1[0], answers[1], eps); + BOOST_CHECK_CLOSE(g2[0][0], answers[2], eps); + BOOST_CHECK_CLOSE(g3[0][0][0], answers[3], eps); BOOST_CHECK_CLOSE(g4[0][0][0][0], answers[4], eps); } BOOST_AUTO_TEST_CASE_TEMPLATE(test_erfc_inv, T, all_float_types) @@ -102,9 +102,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_erfc_inv, T, all_float_types) auto g4 = grad_nd<4>(y, &x); BOOST_CHECK_CLOSE(y.item(), answers[0], eps); - BOOST_CHECK_CLOSE(g1[0]->item(), answers[1], eps); - BOOST_CHECK_CLOSE(g2[0][0]->item(), answers[2], eps); - BOOST_CHECK_CLOSE(g3[0][0][0]->item(), answers[3], eps); + BOOST_CHECK_CLOSE(g1[0], answers[1], eps); + BOOST_CHECK_CLOSE(g2[0][0], answers[2], eps); + BOOST_CHECK_CLOSE(g3[0][0][0], answers[3], eps); BOOST_CHECK_CLOSE(g4[0][0][0][0], answers[4], eps); } BOOST_AUTO_TEST_SUITE_END() diff --git a/test/test_autodiff_reverse_flat_linear_allocator.cpp b/test/test_reverse_mode_autodiff_flat_linear_allocator.cpp similarity index 85% rename from test/test_autodiff_reverse_flat_linear_allocator.cpp rename to test/test_reverse_mode_autodiff_flat_linear_allocator.cpp index 06105baee..7d6695ec4 100644 --- a/test/test_autodiff_reverse_flat_linear_allocator.cpp +++ b/test/test_reverse_mode_autodiff_flat_linear_allocator.cpp @@ -11,18 +11,19 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(flat_linear_allocator_constructors, T, all_float_t float_allocator.emplace_back(rng.next()); } - BOOST_CHECK_EQUAL(float_allocator.size(), 2 * buffer_size - buffer_size / 2); - BOOST_CHECK_EQUAL(float_allocator.capacity(), 2 * buffer_size); + BOOST_CHECK_EQUAL(float_allocator.size(), + static_cast(2 * buffer_size - buffer_size / 2)); + BOOST_CHECK_EQUAL(float_allocator.capacity(), static_cast(2 * buffer_size)); float_allocator.clear(); - BOOST_CHECK_EQUAL(float_allocator.size(), 0); + BOOST_CHECK_EQUAL(float_allocator.size(), static_cast(0)); BOOST_CHECK_EQUAL(float_allocator.capacity(), buffer_size); for (size_t i = 0; i < 2 * buffer_size - buffer_size / 2; i++) { float_allocator.emplace_back(rng.next()); } float_allocator.reset(); - BOOST_CHECK_EQUAL(float_allocator.size(), 0); + BOOST_CHECK_EQUAL(float_allocator.size(), static_cast(0)); BOOST_CHECK_EQUAL(float_allocator.capacity(), 2 * buffer_size); for (size_t i = 0; i < 2 * buffer_size - buffer_size / 2; i++) { @@ -42,20 +43,20 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(flat_linear_allocator_test_emplace, T, all_float_t rdiff_detail::flat_linear_allocator float_allocator{}; std::vector test_vector; - for (int i = 0; i < 2 * buffer_size - 1; i++) { + for (size_t i = 0; i < 2 * buffer_size - 1; i++) { test_vector.push_back(rng.next()); float_allocator.emplace_back(test_vector[i]); } auto it1 = float_allocator.template emplace_back_n<4>(); - for (int i = 0; i < 4; i++) { + for (size_t i = 0; i < 4; i++) { T literal = rng.next(); test_vector.push_back(literal); *(it1 + i) = literal; } auto it2 = float_allocator.emplace_back_n(buffer_size); - for (int i = 0; i < buffer_size; i++) { + for (size_t i = 0; i < buffer_size; i++) { T literal = rng.next(); test_vector.push_back(literal); *(it2 + i) = literal; @@ -69,7 +70,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(flat_linear_allocator_test_emplace, T, all_float_t *alloc_it); // should be ok to check floats like this since they are expected to be the same. } - for (int i = 0; i < test_vector.size(); i++) { + for (size_t i = 0; i < test_vector.size(); i++) { BOOST_CHECK_EQUAL(test_vector[i], float_allocator[i]); // check random access aswell; } @@ -85,7 +86,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(flat_linear_allocator_test_checkpointing, T, all_f std::vector expected_value_at_checkpoint; size_t ckp_id = 0; - for (int i = 0; i < 2 * buffer_size; i++) { + for (size_t i = 0; i < 2 * buffer_size; i++) { T literal = rng.next(); float_allocator.emplace_back(literal); if (ckp_id < checkpoint_indices.size() && i == checkpoint_indices[ckp_id]) { @@ -94,7 +95,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(flat_linear_allocator_test_checkpointing, T, all_f ++ckp_id; } } - for (int i = 0; i < checkpoint_indices.size(); i++) { + for (size_t i = 0; i < checkpoint_indices.size(); i++) { auto it = float_allocator.checkpoint_at(i); BOOST_CHECK_EQUAL(*it, expected_value_at_checkpoint[i]); } diff --git a/test/test_stl_support.cpp b/test/test_reverse_mode_autodiff_stl_support.cpp similarity index 97% rename from test/test_stl_support.cpp rename to test/test_reverse_mode_autodiff_stl_support.cpp index 8f9713a3b..85d5c30c3 100644 --- a/test/test_stl_support.cpp +++ b/test/test_reverse_mode_autodiff_stl_support.cpp @@ -1,6 +1,8 @@ #include "test_autodiff_reverse.hpp" +#include #include #include + BOOST_AUTO_TEST_SUITE(test_stl_supported_functions) using namespace rdiff; @@ -263,10 +265,12 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_frexp, T, all_float_types) tape.zero_grad(); x_func.backward(); - int i3, i4; + // int i3, i4; + int i4; + frexp(x_v, &i4); // T h = 0.001; // T f1 = frexp(x_v + h, &i3); - T f2 = frexp(x_v, &i4); + // T f2 = frexp(x_v, &i4); T expected_deriv = 1.0 / pow(2.0, i4); BOOST_REQUIRE_CLOSE_FRACTION(x_rvar.adjoint(), expected_deriv, boost_close_tol()); @@ -307,7 +311,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_sin, T, all_float_types) rvar x_func = sin(x_rvar); - BOOST_REQUIRE_CLOSE_FRACTION(x_func.item(), test_func_v, 1 * boost_close_tol()); + BOOST_REQUIRE_CLOSE_FRACTION(x_func.item(), test_func_v, 10 * boost_close_tol()); gradient_tape& tape = get_active_tape(); @@ -331,7 +335,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_tan, T, all_float_types) rvar x_func = tan(x_rvar); - BOOST_REQUIRE_CLOSE_FRACTION(x_func.item(), test_func_v, boost_close_tol()); + BOOST_REQUIRE_CLOSE_FRACTION(x_func.item(), test_func_v, 10 * boost_close_tol()); gradient_tape& tape = get_active_tape(); @@ -340,7 +344,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_tan, T, all_float_types) T expected_deriv = 2 * x_v * 1 / (cos(x_v * x_v) * cos(x_v * x_v)); test_func_2.backward(); - BOOST_REQUIRE_CLOSE_FRACTION(x_rvar.adjoint(), expected_deriv, 1 * boost_close_tol()); + BOOST_REQUIRE_CLOSE_FRACTION(x_rvar.adjoint(), expected_deriv, 10 * boost_close_tol()); tape.clear(); } @@ -496,6 +500,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_trunc, T, all_float_types) BOOST_AUTO_TEST_CASE_TEMPLATE(test_fmod, T, all_float_types) { + using boost::math::tools::fmod_workaround; RandomSample rng_x{-1, 1}; RandomSample rng_y{-1, 1}; @@ -505,13 +510,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_fmod, T, all_float_types) y_v = rng_y.next(); } - T expected_fmod = fmod(x_v, y_v); + T expected_fmod = fmod(x_v, y_v); rvar x_rvar = x_v; rvar y_rvar = y_v; rvar fmod_rvar = fmod(x_rvar, y_rvar); - BOOST_REQUIRE_CLOSE_FRACTION(fmod_rvar.item(), expected_fmod, boost_close_tol()); + BOOST_REQUIRE_CLOSE_FRACTION(fmod_rvar.item(), expected_fmod, 10 * boost_close_tol()); gradient_tape& tape = get_active_tape(); tape.zero_grad(); @@ -520,22 +525,22 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_fmod, T, all_float_types) T dx_expected = 1.0; T dy_expected = -1.0 * trunc(x_v / y_v); - BOOST_REQUIRE_CLOSE_FRACTION(x_rvar.adjoint(), dx_expected, boost_close_tol()); - BOOST_REQUIRE_CLOSE_FRACTION(y_rvar.adjoint(), dy_expected, boost_close_tol()); + BOOST_REQUIRE_CLOSE_FRACTION(x_rvar.adjoint(), dx_expected, 10 * boost_close_tol()); + BOOST_REQUIRE_CLOSE_FRACTION(y_rvar.adjoint(), dy_expected, 10 * boost_close_tol()); rvar fmod_rvar_float = fmod(x_rvar, y_v); - BOOST_REQUIRE_CLOSE_FRACTION(fmod_rvar_float.item(), expected_fmod, boost_close_tol()); + BOOST_REQUIRE_CLOSE_FRACTION(fmod_rvar_float.item(), expected_fmod, 10 * boost_close_tol()); tape.zero_grad(); fmod_rvar_float.backward(); - BOOST_REQUIRE_CLOSE_FRACTION(x_rvar.adjoint(), dx_expected, boost_close_tol()); + BOOST_REQUIRE_CLOSE_FRACTION(x_rvar.adjoint(), dx_expected, 10 * boost_close_tol()); rvar fmod_float_rvar = fmod(x_v, y_rvar); - BOOST_REQUIRE_CLOSE_FRACTION(fmod_float_rvar.item(), expected_fmod, boost_close_tol()); + BOOST_REQUIRE_CLOSE_FRACTION(fmod_float_rvar.item(), expected_fmod, 10 * boost_close_tol()); tape.zero_grad(); fmod_float_rvar.backward(); - BOOST_REQUIRE_CLOSE_FRACTION(y_rvar.adjoint(), dy_expected, boost_close_tol()); + BOOST_REQUIRE_CLOSE_FRACTION(y_rvar.adjoint(), dy_expected, 10 * boost_close_tol()); tape.clear(); } @@ -854,7 +859,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_func_1_first_derivative, T, all_float_types) rvar x_ad = x; rvar y_ad = y; - T fv = test_function_1(x, y); + //T fv = test_function_1(x, y); rvar f_ad = test_function_1(x_ad, y_ad); auto grad_f_analytical = grad_test_func_1_analytical(x, y); /* intended use case */ @@ -888,7 +893,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_func_1_second_derivative_and_hessian, T, all_ rvar x_ad = x; rvar y_ad = y; - T fv = test_function_1(x, y); + //T fv = test_function_1(x, y); rvar f_ad = test_function_1(x_ad, y_ad); gradient_tape& tape = get_active_tape(); tape.zero_grad(); @@ -918,7 +923,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_func_1_order_3_der, T, all_float_types) rvar x_ad = x; rvar y_ad = y; - T fv = test_function_1(x, y); + //T fv = test_function_1(x, y); rvar f_ad = test_function_1(x_ad, y_ad); gradient_tape& tape = get_active_tape(); tape.zero_grad();