From c33ef74162cc6afbe77b6056c67d5b77097f5745 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Thu, 28 Apr 2022 18:11:17 -0400 Subject: [PATCH 01/20] Adds scalar root finding function --- stan/math/prim/functor.hpp | 2 +- stan/math/prim/functor/root_finder.hpp | 97 +++++++++++++++ stan/math/rev/fun/abs.hpp | 7 ++ stan/math/rev/functor.hpp | 1 + stan/math/rev/functor/root_finder.hpp | 112 ++++++++++++++++++ .../math/rev/functor/root_finder_test.cpp | 29 +++++ 6 files changed, 247 insertions(+), 1 deletion(-) create mode 100644 stan/math/prim/functor/root_finder.hpp create mode 100644 stan/math/rev/functor/root_finder.hpp create mode 100644 test/unit/math/rev/functor/root_finder_test.cpp diff --git a/stan/math/prim/functor.hpp b/stan/math/prim/functor.hpp index 0ec5c343ff7..47c7116d257 100644 --- a/stan/math/prim/functor.hpp +++ b/stan/math/prim/functor.hpp @@ -26,5 +26,5 @@ #include #include #include - +#include #endif diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp new file mode 100644 index 00000000000..9ed69436ece --- /dev/null +++ b/stan/math/prim/functor/root_finder.hpp @@ -0,0 +1,97 @@ +#ifndef STAN_MATH_PRIM_FUNCTOR_ROOT_FINDER_HPP +#define STAN_MATH_PRIM_FUNCTOR_ROOT_FINDER_HPP + +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + + +// TODO: Can just have one signature with another function for constructing tuple + +// When caller provides functions for value and derivative and second derivative +template * = nullptr> +auto root_finder_tol(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, GuessScalar guess, MinScalar min, MaxScalar max, int digits, + std::uintmax_t& max_iter, Types&&... args) { + check_bounded("root_finder", "initial guess", guess, min, max); + return_type_t ret = 0; + auto f_plus_div = [&f, &f_div, &f_div_div, &args...](auto&& g) { + return std::make_tuple(f(g, args...), f_div(g, args...), f_div_div(g, args...)); + }; + try { + ret = boost::math::tools::halley_iterate(f_plus_div, + guess, min, max, digits, max_iter); + } catch (const std::exception& e) { + throw e; + } + return ret; +} +/* +// When caller provides functions for value and derivative +template * = nullptr> +double root_finder_tol(F&& f, FDiv&& f_div, GuessScalar guess, MinScalar min, MaxScalar max, int digits, + std::uintmax_t& max_iter, Types&&... args) { + check_bounded("root_finder", "initial guess", guess, min, max); + double ret = 0; + auto f_plus_div = [&f, &f_div, &args...](auto&& g) { + return std::make_tuple(f(g, args...), f_div(g, args...)); + }; + try { + ret = boost::math::tools::halley_iterate(f_plus_div, + guess, min, max, digits, max_iter); + } catch (const std::exception& e) { + throw e; + } + return ret; +} + +// When supplied function returns pair +template * = nullptr> +double root_finder_tol(F&& f, GuessScalar guess, MinScalar min, MaxScalar max, int digits, + std::uintmax_t& max_iter, Types&&... args) { + check_bounded("root_finder", "initial guess", guess, min, max); + double ret = 0; + auto f_plus_args = [&f, &args...](auto&& g) { return f(g, args...);}; + try { + ret = boost::math::tools::halley_iterate(f_plus_args, + guess, min, max, digits, max_iter); + } catch (const std::exception& e) { + throw e; + } + return ret; +} +*/ + +// Non-tol versions +template +auto root_finder(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, GuessScalar guess, MinScalar min, MaxScalar max, Types&&... args) { + int digits = 16; + std::uintmax_t max_iter = 100; + return root_finder_tol(f, f_div, f_div_div, guess, min, max, digits, max_iter, args...); +} +/* +template +double root_finder(F&& f, FDiv&& f_div, GuessScalar guess, MinScalar min, MaxScalar max, Types&&... args) { + int digits = 16; + std::uintmax_t max_iter = 100; + return root_finder_tol(f, f_div, guess, min, max, digits, max_iter, args...); +} + +template +double root_finder(F&& f, GuessScalar guess, MinScalar min, MaxScalar max, Types&&... args) { + int digits = 16; + std::uintmax_t max_iter = 100; + return root_finder_tol(f, guess, min, max, digits, max_iter, args...); +} +*/ +} +} +#endif diff --git a/stan/math/rev/fun/abs.hpp b/stan/math/rev/fun/abs.hpp index 919fb320bf9..156271595de 100644 --- a/stan/math/rev/fun/abs.hpp +++ b/stan/math/rev/fun/abs.hpp @@ -10,6 +10,13 @@ namespace stan { namespace math { + inline auto frexp(stan::math::var x, int* exponent) noexcept { + return std::frexp(x.val(), exponent); + } + inline int sign(stan::math::var z) { + return (z == 0) ? 0 : z < 0 ? -1 : 1; + } + /** * Return the absolute value of the variable (std). * diff --git a/stan/math/rev/functor.hpp b/stan/math/rev/functor.hpp index d494bfa5c6b..5e55c3a3506 100644 --- a/stan/math/rev/functor.hpp +++ b/stan/math/rev/functor.hpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #endif diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp new file mode 100644 index 00000000000..2395d62169d --- /dev/null +++ b/stan/math/rev/functor/root_finder.hpp @@ -0,0 +1,112 @@ +#ifndef STAN_MATH_REV_FUNCTOR_ROOT_FINDER_HPP +#define STAN_MATH_REV_FUNCTOR_ROOT_FINDER_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + + /* +// TODO: Can just have one signature with another function for constructing +// tuple +template < + typename F, typename FDiv, typename FDivDiv, typename GuessScalar, + typename MinScalar, typename MaxScalar, typename... Types, + require_any_st_var* = nullptr> +auto root_finder_tol(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, + GuessScalar guess, MinScalar min, MaxScalar max, + int digits, std::uintmax_t& max_iter, Types&&... args) { + check_bounded("root_finder", "initial guess", guess, min, max); + return_type_t ret = 0; + auto f_plus_div = [&f, &f_div, &f_div_div, &args...](auto&& g) { + return std::make_tuple(f(g, args...), f_div(g, args...), + f_div_div(g, args...)); + }; + try { + ret = boost::math::tools::halley_iterate(f_plus_div, guess, min, max, + digits, max_iter); + } catch (const std::exception& e) { + throw e; + } + return ret; +} +*/ + +template < + typename F, typename FDiv, typename FDivDiv, typename GuessScalar, + typename MinScalar, typename MaxScalar, typename... Types, + require_any_st_var* = nullptr> +auto root_finder_tol(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, + GuessScalar guess, MinScalar min, MaxScalar max, + int digits, std::uintmax_t& max_iter, Types&&... args) { + check_bounded("root_finder", "initial guess", guess, min, max); + check_nonnegative("root_finder", "digits", digits); + check_positive("root_finder", "max_iter", max_iter); + auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(args)...)); + auto args_vals_tuple = apply( + [&](const auto&... args) { + return std::make_tuple(to_ref(value_of(args))...); + }, + *arena_args_tuple); + // Solve the system + double theta_dbl = apply( + [&](auto&&... vals) { + return root_finder(f, f_div, f_div_div, value_of(guess), value_of(min), + value_of(max), vals...); + }, + args_vals_tuple); + + auto f_wrt_x = [&](auto&& x) { + return apply([&](auto&&... args) { return f(x, args...); }, + args_vals_tuple); + }; + + double Jf_x; + double f_x; + + //jacobian(f_wrt_x, theta_dbl, f_x, Jf_x); + { + nested_rev_autodiff nested; + stan::math::var x_var(theta_dbl); + stan::math::var fx_var = f_wrt_x(x_var); + f_x = fx_var.val(); + fx_var.grad(); + Jf_x = x_var.adj(); + } + stan::math::var ret = theta_dbl; + //auto Jf_x_T_lu_ptr = make_unsafe_chainable_ptr(Jf_x.transpose().partialPivLu()); // Lu + auto Jf_x_T_lu_ptr = 1. / Jf_x; // Lu + + reverse_pass_callback( + [f, ret, arena_args_tuple, Jf_x_T_lu_ptr]() mutable { + // Eigen::VectorXd eta = -Jf_x_T_lu_ptr->solve(ret.adj().eval()); + double eta = -(ret.adj() * Jf_x_T_lu_ptr); + + // Contract with Jacobian of f with respect to y using a nested reverse + // autodiff pass. + { + nested_rev_autodiff rev; + + double ret_val = ret.val(); + auto x_nrad_ = apply( + [&ret_val, &f](const auto&... args) { + return eval(f(ret_val, args...)); + }, + *arena_args_tuple); + x_nrad_.adj() = eta; + grad(); + } + }); + + return ret; +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp new file mode 100644 index 00000000000..cec6e665aaf --- /dev/null +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -0,0 +1,29 @@ +#include +#include +#include + + +TEST(RevFunctor, root_finder) { + using stan::math::root_finder; + using stan::math::var; + // return cube root of x using 1st and 2nd derivatives and Halley. + //using namespace std; // Help ADL of std functions. + var x = 27; + int exponent; + std::frexp(x.val(), &exponent); // Get exponent of z (ignore mantissa). + var guess = ldexp(1., exponent / 3); // Rough guess is to divide the exponent by three. + var min = ldexp(0.5, exponent / 3); // Minimum possible value is half our guess. + var max = ldexp(2., exponent / 3); // Maximum possible value is twice our guess. + const int digits = std::numeric_limits::digits; // Maximum possible binary digits accuracy for type T. + // digits used to control how accurate to try to make the result. + int get_digits = static_cast(digits * 0.4); // Accuracy triples with each step, so stop when just + // over one third of the digits are correct. + std::uintmax_t maxit = 20; + auto f = [](const auto& g, const auto x){ return g * g * g - x; }; + auto f_g = [](const auto& g, const auto x){ return 3 * g * g; }; + auto f_gg = [](const auto& g, const auto x){ return 6 * g; }; + var result = root_finder(f, f_g, f_gg, guess, min, max, x); + result.grad(); + // val: 27 adj:0.037037 + std::cout << "\nval: " << x.val() << " adj:" << x.adj() << "\n"; +} From b015dae746385987071cd1c16c9a2e45e8614aa9 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Thu, 28 Apr 2022 18:24:10 -0400 Subject: [PATCH 02/20] cleanup --- stan/math/rev/functor/root_finder.hpp | 92 ++++++++++++--------------- 1 file changed, 40 insertions(+), 52 deletions(-) diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp index 2395d62169d..1245408f1b0 100644 --- a/stan/math/rev/functor/root_finder.hpp +++ b/stan/math/rev/functor/root_finder.hpp @@ -12,29 +12,29 @@ namespace stan { namespace math { - /* +/* // TODO: Can just have one signature with another function for constructing // tuple template < - typename F, typename FDiv, typename FDivDiv, typename GuessScalar, - typename MinScalar, typename MaxScalar, typename... Types, - require_any_st_var* = nullptr> + typename F, typename FDiv, typename FDivDiv, typename GuessScalar, + typename MinScalar, typename MaxScalar, typename... Types, + require_any_st_var* = nullptr> auto root_finder_tol(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, - GuessScalar guess, MinScalar min, MaxScalar max, - int digits, std::uintmax_t& max_iter, Types&&... args) { - check_bounded("root_finder", "initial guess", guess, min, max); - return_type_t ret = 0; - auto f_plus_div = [&f, &f_div, &f_div_div, &args...](auto&& g) { - return std::make_tuple(f(g, args...), f_div(g, args...), - f_div_div(g, args...)); - }; - try { - ret = boost::math::tools::halley_iterate(f_plus_div, guess, min, max, - digits, max_iter); - } catch (const std::exception& e) { - throw e; - } - return ret; + GuessScalar guess, MinScalar min, MaxScalar max, + int digits, std::uintmax_t& max_iter, Types&&... args) { +check_bounded("root_finder", "initial guess", guess, min, max); +return_type_t ret = 0; +auto f_plus_div = [&f, &f_div, &f_div_div, &args...](auto&& g) { + return std::make_tuple(f(g, args...), f_div(g, args...), + f_div_div(g, args...)); +}; +try { + ret = boost::math::tools::halley_iterate(f_plus_div, guess, min, max, + digits, max_iter); +} catch (const std::exception& e) { + throw e; +} +return ret; } */ @@ -58,51 +58,39 @@ auto root_finder_tol(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, double theta_dbl = apply( [&](auto&&... vals) { return root_finder(f, f_div, f_div_div, value_of(guess), value_of(min), - value_of(max), vals...); + value_of(max), vals...); }, args_vals_tuple); - - auto f_wrt_x = [&](auto&& x) { - return apply([&](auto&&... args) { return f(x, args...); }, - args_vals_tuple); - }; - double Jf_x; double f_x; - - //jacobian(f_wrt_x, theta_dbl, f_x, Jf_x); { nested_rev_autodiff nested; stan::math::var x_var(theta_dbl); - stan::math::var fx_var = f_wrt_x(x_var); + stan::math::var fx_var = apply( + [&x_var, &f](auto&&... args) { return f(x_var, std::move(args)...); }, + std::move(args_vals_tuple)); f_x = fx_var.val(); fx_var.grad(); Jf_x = x_var.adj(); } stan::math::var ret = theta_dbl; - //auto Jf_x_T_lu_ptr = make_unsafe_chainable_ptr(Jf_x.transpose().partialPivLu()); // Lu - auto Jf_x_T_lu_ptr = 1. / Jf_x; // Lu - - reverse_pass_callback( - [f, ret, arena_args_tuple, Jf_x_T_lu_ptr]() mutable { - // Eigen::VectorXd eta = -Jf_x_T_lu_ptr->solve(ret.adj().eval()); - double eta = -(ret.adj() * Jf_x_T_lu_ptr); - - // Contract with Jacobian of f with respect to y using a nested reverse - // autodiff pass. - { - nested_rev_autodiff rev; - - double ret_val = ret.val(); - auto x_nrad_ = apply( - [&ret_val, &f](const auto&... args) { - return eval(f(ret_val, args...)); - }, - *arena_args_tuple); - x_nrad_.adj() = eta; - grad(); - } - }); + reverse_pass_callback([f, ret, arena_args_tuple, Jf_x]() mutable { + // Eigen::VectorXd eta = -Jf_x_T_lu_ptr->solve(ret.adj().eval()); + double eta = -(ret.adj() / Jf_x); + // Contract with Jacobian of f with respect to y using a nested reverse + // autodiff pass. + { + nested_rev_autodiff rev; + double ret_val = ret.val(); + auto x_nrad_ = apply( + [&ret_val, &f](const auto&... args) { + return eval(f(ret_val, args...)); + }, + *arena_args_tuple); + x_nrad_.adj() = eta; + grad(); + } + }); return ret; } From bd9fc696743825aaa035f08c9b450fe7cf8df20d Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Fri, 29 Apr 2022 18:34:42 -0400 Subject: [PATCH 03/20] cleanup and have the functors pass as a set of tuples --- stan/math/fwd/fun.hpp | 2 + stan/math/fwd/fun/abs.hpp | 1 + stan/math/fwd/fun/frexp.hpp | 17 ++ stan/math/fwd/fun/sign.hpp | 19 +++ stan/math/prim/functor/root_finder.hpp | 145 +++++++++--------- stan/math/rev/fun.hpp | 2 + stan/math/rev/fun/abs.hpp | 7 - stan/math/rev/fun/frexp.hpp | 14 ++ stan/math/rev/fun/sign.hpp | 14 ++ stan/math/rev/functor/root_finder.hpp | 110 ++++++------- .../math/mix/functor/root_finder_test.cpp | 71 +++++++++ .../math/rev/functor/root_finder_test.cpp | 30 ++-- 12 files changed, 285 insertions(+), 147 deletions(-) create mode 100644 stan/math/fwd/fun/frexp.hpp create mode 100644 stan/math/fwd/fun/sign.hpp create mode 100644 stan/math/rev/fun/frexp.hpp create mode 100644 stan/math/rev/fun/sign.hpp create mode 100644 test/unit/math/mix/functor/root_finder_test.cpp diff --git a/stan/math/fwd/fun.hpp b/stan/math/fwd/fun.hpp index 0237d1a8af4..e529fcb6cd3 100644 --- a/stan/math/fwd/fun.hpp +++ b/stan/math/fwd/fun.hpp @@ -38,6 +38,7 @@ #include #include #include +#include #include #include #include @@ -103,6 +104,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/fwd/fun/abs.hpp b/stan/math/fwd/fun/abs.hpp index 1f50b149fd6..5817ca62795 100644 --- a/stan/math/fwd/fun/abs.hpp +++ b/stan/math/fwd/fun/abs.hpp @@ -6,6 +6,7 @@ #include #include #include +#include #include namespace stan { diff --git a/stan/math/fwd/fun/frexp.hpp b/stan/math/fwd/fun/frexp.hpp new file mode 100644 index 00000000000..3d0b44b121b --- /dev/null +++ b/stan/math/fwd/fun/frexp.hpp @@ -0,0 +1,17 @@ +#ifndef STAN_MATH_FWD_FUN_FREXP_HPP +#define STAN_MATH_FWD_FUN_FREXP_HPP + +#include +#include +#include + +namespace stan { +namespace math { + +template +inline auto frexp(const fvar& x, int* exponent) noexcept { + return std::frexp(value_of_rec(x), exponent); +} +} +} +#endif diff --git a/stan/math/fwd/fun/sign.hpp b/stan/math/fwd/fun/sign.hpp new file mode 100644 index 00000000000..b36067b367d --- /dev/null +++ b/stan/math/fwd/fun/sign.hpp @@ -0,0 +1,19 @@ +#ifndef STAN_MATH_FWD_FUN_SIGN_HPP +#define STAN_MATH_FWD_FUN_SIGN_HPP + +#include +#include +#include + +namespace stan { +namespace math { + +template +inline auto sign(const fvar& x) { + double x_val = value_of_rec(x); + return (0. < x_val) - (x_val < 0.); +} + +} +} +#endif diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index 9ed69436ece..fef56d3139a 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -10,88 +10,87 @@ namespace stan { namespace math { - - -// TODO: Can just have one signature with another function for constructing tuple - -// When caller provides functions for value and derivative and second derivative -template * = nullptr> -auto root_finder_tol(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, GuessScalar guess, MinScalar min, MaxScalar max, int digits, - std::uintmax_t& max_iter, Types&&... args) { - check_bounded("root_finder", "initial guess", guess, min, max); - return_type_t ret = 0; - auto f_plus_div = [&f, &f_div, &f_div_div, &args...](auto&& g) { - return std::make_tuple(f(g, args...), f_div(g, args...), f_div_div(g, args...)); - }; - try { - ret = boost::math::tools::halley_iterate(f_plus_div, - guess, min, max, digits, max_iter); - } catch (const std::exception& e) { - throw e; - } - return ret; -} -/* -// When caller provides functions for value and derivative -template * = nullptr> -double root_finder_tol(F&& f, FDiv&& f_div, GuessScalar guess, MinScalar min, MaxScalar max, int digits, - std::uintmax_t& max_iter, Types&&... args) { - check_bounded("root_finder", "initial guess", guess, min, max); - double ret = 0; - auto f_plus_div = [&f, &f_div, &args...](auto&& g) { - return std::make_tuple(f(g, args...), f_div(g, args...)); - }; - try { - ret = boost::math::tools::halley_iterate(f_plus_div, - guess, min, max, digits, max_iter); - } catch (const std::exception& e) { - throw e; - } - return ret; +namespace internal { +template +inline auto func_with_derivs(Tuple&& f_tuple, Args&&... args) { + return stan::math::apply( + [&args...](auto&&... funcs) { + return [&args..., &funcs...](auto&& g) { + return std::make_tuple(funcs(g, args...)...); + }; + }, + f_tuple); } +} // namespace internal -// When supplied function returns pair -template * = nullptr> -double root_finder_tol(F&& f, GuessScalar guess, MinScalar min, MaxScalar max, int digits, - std::uintmax_t& max_iter, Types&&... args) { +/** + * Solve for root using Boost's Halley method + * @tparam FTuple A tuple holding functors whose signatures all match + * `(GuessScalar g, Types&&... Args)`. + * @tparam GuessScalar Scalar type + * @tparam MinScalar Scalar type + * @tparam MaxScalar Scalar type + * @tparam Types Arg types to pass to functors in `f_tuple` + * @param f_tuple A tuple of functors to calculate the value and any derivates + * needed. + * @param guess An initial guess at the root value + * @param min The minimum possible value for the result, this is used as an + * initial lower bracket + * @param max The maximum possible value for the result, this is used as an + * initial upper bracket + * @param digits The desired number of binary digits precision + * @param max_iter An optional maximum number of iterations to perform. On exit, + * this is updated to the actual number of iterations performed + * @param args Parameter pack of arguments to pass the the functors in `f_tuple` + */ +template * = nullptr> +auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, const MinScalar min, + const MaxScalar max, const int digits, std::uintmax_t& max_iter, + Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); - double ret = 0; - auto f_plus_args = [&f, &args...](auto&& g) { return f(g, args...);}; + check_positive("root_finder", "digits", digits); + check_positive("root_finder", "max_iter", max_iter); + return_type_t ret = 0; + auto f_plus_div = internal::func_with_derivs(f_tuple, args...); try { - ret = boost::math::tools::halley_iterate(f_plus_args, - guess, min, max, digits, max_iter); + ret = boost::math::tools::halley_iterate(f_plus_div, guess, min, max, + digits, max_iter); } catch (const std::exception& e) { throw e; } return ret; } -*/ -// Non-tol versions -template -auto root_finder(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, GuessScalar guess, MinScalar min, MaxScalar max, Types&&... args) { - int digits = 16; - std::uintmax_t max_iter = 100; - return root_finder_tol(f, f_div, f_div_div, guess, min, max, digits, max_iter, args...); -} -/* -template -double root_finder(F&& f, FDiv&& f_div, GuessScalar guess, MinScalar min, MaxScalar max, Types&&... args) { - int digits = 16; - std::uintmax_t max_iter = 100; - return root_finder_tol(f, f_div, guess, min, max, digits, max_iter, args...); -} - -template -double root_finder(F&& f, GuessScalar guess, MinScalar min, MaxScalar max, Types&&... args) { - int digits = 16; - std::uintmax_t max_iter = 100; - return root_finder_tol(f, guess, min, max, digits, max_iter, args...); -} -*/ -} +/** + * Solve for root using Boost Halley method with default values for the + * tolerances + * @tparam FTuple A tuple holding functors whose signatures all match + * `(GuessScalar g, Types&&... Args)`. + * @tparam GuessScalar Scalar type + * @tparam MinScalar Scalar type + * @tparam MaxScalar Scalar type + * @tparam Types Arg types to pass to functors in `f_tuple` + * @param f_tuple A tuple of functors to calculate the value and any derivates + * needed. + * @param guess An initial guess at the root value + * @param min The minimum possible value for the result, this is used as an + * initial lower bracket + * @param max The maximum possible value for the result, this is used as an + * initial upper bracket + * @param args Parameter pack of arguments to pass the the functors in `f_tuple` + */ +template +auto root_finder(FTuple&& f_tuple, const GuessScalar guess, const MinScalar min, + const MaxScalar max, Types&&... args) { + constexpr int digits = 16; + std::uintmax_t max_iter = std::numeric_limits::max(); + return root_finder_tol(std::forward(f_tuple), guess, min, max, digits, + max_iter, std::forward(args)...); } +} // namespace math +} // namespace stan #endif diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index de47503372f..0cb81cf2903 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -67,6 +67,7 @@ #include #include #include +#include #include #include #include @@ -161,6 +162,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/abs.hpp b/stan/math/rev/fun/abs.hpp index 156271595de..919fb320bf9 100644 --- a/stan/math/rev/fun/abs.hpp +++ b/stan/math/rev/fun/abs.hpp @@ -10,13 +10,6 @@ namespace stan { namespace math { - inline auto frexp(stan::math::var x, int* exponent) noexcept { - return std::frexp(x.val(), exponent); - } - inline int sign(stan::math::var z) { - return (z == 0) ? 0 : z < 0 ? -1 : 1; - } - /** * Return the absolute value of the variable (std). * diff --git a/stan/math/rev/fun/frexp.hpp b/stan/math/rev/fun/frexp.hpp new file mode 100644 index 00000000000..ba5437d1f43 --- /dev/null +++ b/stan/math/rev/fun/frexp.hpp @@ -0,0 +1,14 @@ +#ifndef STAN_MATH_REV_FUN_FREXP_HPP +#define STAN_MATH_REV_FUN_FREXP_HPP + +#include +#include + +namespace stan { +namespace math { +inline auto frexp(stan::math::var x, int* exponent) noexcept { + return std::frexp(x.val(), exponent); +} +} +} +#endif diff --git a/stan/math/rev/fun/sign.hpp b/stan/math/rev/fun/sign.hpp new file mode 100644 index 00000000000..0c3bb816084 --- /dev/null +++ b/stan/math/rev/fun/sign.hpp @@ -0,0 +1,14 @@ +#ifndef STAN_MATH_REV_FUN_SIGN_HPP +#define STAN_MATH_REV_FUN_SIGN_HPP + +#include +#include + +namespace stan { +namespace math { +inline int sign(stan::math::var z) { + return (z == 0) ? 0 : z < 0 ? -1 : 1; +} +} +} +#endif diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp index 1245408f1b0..435c3ea13a1 100644 --- a/stan/math/rev/functor/root_finder.hpp +++ b/stan/math/rev/functor/root_finder.hpp @@ -12,41 +12,36 @@ namespace stan { namespace math { -/* -// TODO: Can just have one signature with another function for constructing -// tuple +/** + * var specialization for root solving using Boost's Halley method + * @tparam FTuple A tuple holding functors whose signatures all match + * `(GuessScalar g, Types&&... Args)`. + * @tparam GuessScalar Scalar type + * @tparam MinScalar Scalar type + * @tparam MaxScalar Scalar type + * @tparam Types Arg types to pass to functors in `f_tuple` + * @param f_tuple A tuple of functors to calculate the value and any derivates + * needed. + * @param guess An initial guess at the root value + * @param min The minimum possible value for the result, this is used as an + * initial lower bracket + * @param max The maximum possible value for the result, this is used as an + * initial upper bracket + * @param digits The desired number of binary digits precision + * @param max_iter An optional maximum number of iterations to perform. On exit, + * this is updated to the actual number of iterations performed + * @param args Parameter pack of arguments to pass the the functors in `f_tuple` + */ template < - typename F, typename FDiv, typename FDivDiv, typename GuessScalar, - typename MinScalar, typename MaxScalar, typename... Types, - require_any_st_var* = nullptr> -auto root_finder_tol(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, - GuessScalar guess, MinScalar min, MaxScalar max, - int digits, std::uintmax_t& max_iter, Types&&... args) { -check_bounded("root_finder", "initial guess", guess, min, max); -return_type_t ret = 0; -auto f_plus_div = [&f, &f_div, &f_div_div, &args...](auto&& g) { - return std::make_tuple(f(g, args...), f_div(g, args...), - f_div_div(g, args...)); -}; -try { - ret = boost::math::tools::halley_iterate(f_plus_div, guess, min, max, - digits, max_iter); -} catch (const std::exception& e) { - throw e; -} -return ret; -} -*/ - -template < - typename F, typename FDiv, typename FDivDiv, typename GuessScalar, - typename MinScalar, typename MaxScalar, typename... Types, - require_any_st_var* = nullptr> -auto root_finder_tol(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, - GuessScalar guess, MinScalar min, MaxScalar max, - int digits, std::uintmax_t& max_iter, Types&&... args) { + typename FTuple, typename GuessScalar, typename MinScalar, + typename MaxScalar, typename... Types, + require_any_st_var* = nullptr, + require_all_stan_scalar_t* = nullptr> +auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, const MinScalar min, + const MaxScalar max, const int digits, std::uintmax_t& max_iter, + Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); - check_nonnegative("root_finder", "digits", digits); + check_positive("root_finder", "digits", digits); check_positive("root_finder", "max_iter", max_iter); auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(args)...)); auto args_vals_tuple = apply( @@ -56,13 +51,14 @@ auto root_finder_tol(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, *arena_args_tuple); // Solve the system double theta_dbl = apply( - [&](auto&&... vals) { - return root_finder(f, f_div, f_div_div, value_of(guess), value_of(min), - value_of(max), vals...); + [&f_tuple, guess_val = value_of(guess), min_val = value_of(min), + max_val = value_of(max)](auto&&... vals) { + return root_finder(f_tuple, guess_val, min_val, max_val, vals...); }, args_vals_tuple); double Jf_x; double f_x; + auto&& f = std::get<0>(f_tuple); { nested_rev_autodiff nested; stan::math::var x_var(theta_dbl); @@ -73,26 +69,30 @@ auto root_finder_tol(F&& f, FDiv&& f_div, FDivDiv&& f_div_div, fx_var.grad(); Jf_x = x_var.adj(); } - stan::math::var ret = theta_dbl; - reverse_pass_callback([f, ret, arena_args_tuple, Jf_x]() mutable { - // Eigen::VectorXd eta = -Jf_x_T_lu_ptr->solve(ret.adj().eval()); - double eta = -(ret.adj() / Jf_x); - // Contract with Jacobian of f with respect to y using a nested reverse - // autodiff pass. - { - nested_rev_autodiff rev; - double ret_val = ret.val(); - auto x_nrad_ = apply( - [&ret_val, &f](const auto&... args) { - return eval(f(ret_val, args...)); - }, - *arena_args_tuple); - x_nrad_.adj() = eta; - grad(); - } - }); - return ret; + /* + * Note: Because we put this on the callback stack, if `f` is a lambda + * its captures must be in Stan's arena memory or trivially destructable. + */ + return make_callback_var(theta_dbl, + [f, arena_args_tuple, Jf_x](auto& ret) mutable { + // Eigen::VectorXd eta = + // -Jf_x_T_lu_ptr->solve(ret.adj().eval()); + double eta = -(ret.adj() / Jf_x); + // Contract with Jacobian of f with respect to y + // using a nested reverse autodiff pass. + { + nested_rev_autodiff rev; + double ret_val = ret.val(); + auto x_nrad_ = apply( + [&ret_val, &f](const auto&... args) { + return eval(f(ret_val, args...)); + }, + *arena_args_tuple); + x_nrad_.adj() = eta; + grad(); + } + }); } } // namespace math diff --git a/test/unit/math/mix/functor/root_finder_test.cpp b/test/unit/math/mix/functor/root_finder_test.cpp new file mode 100644 index 00000000000..454ec8af0a5 --- /dev/null +++ b/test/unit/math/mix/functor/root_finder_test.cpp @@ -0,0 +1,71 @@ +#include + +TEST(MixFun, root_finder_cubed) { + using stan::math::root_finder; + // return cube root of x using 1st and 2nd derivatives and Halley. + //using namespace std; // Help ADL of std functions. + double x = 27; + int exponent; + // Get exponent of z (ignore mantissa). + std::frexp(x, &exponent); + // Rough guess is to divide the exponent by three. + double guess = ldexp(1., exponent / 3); + // Minimum possible value is half our guess. + double min = ldexp(0.5, exponent / 3); + // Maximum possible value is twice our guess. + double max = ldexp(2., exponent / 3); + // Maximum possible binary digits accuracy for type T. + const int digits = std::numeric_limits::digits; + /* digits used to control how accurate to try to make the result. + * Accuracy triples with each step, so stop when just + * over one third of the digits are correct. + */ + int get_digits = static_cast(digits * 0.4); + std::uintmax_t maxit = 20; + auto f = [](const auto& g, const auto& x){ return g * g * g - x; }; + auto f_g = [](const auto& g, const auto& x){ return 3 * g * g; }; + auto f_gg = [](const auto& g, const auto& x){ return 6 * g; }; + auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& xx) { + using x_t = std::decay_t; + return root_finder(std::make_tuple(f, f_g, f_gg), + x_t(guess), + x_t(min), + x_t(max), xx); + }; + stan::test::expect_ad(full_f, x); +} + +TEST(MixFun, root_finder_fifth) { + using stan::math::root_finder; + // return cube root of x using 1st and 2nd derivatives and Halley. + //using namespace std; // Help ADL of std functions. + double x = 27; + int exponent; + // Get exponent of z (ignore mantissa). + std::frexp(x, &exponent); + // Rough guess is to divide the exponent by three. + double guess = ldexp(1., exponent / 5); + // Minimum possible value is half our guess. + double min = ldexp(0.5, exponent / 5); + // Maximum possible value is twice our guess. + double max = ldexp(2., exponent / 5); + // Maximum possible binary digits accuracy for type T. + const int digits = std::numeric_limits::digits; + /* digits used to control how accurate to try to make the result. + * Accuracy triples with each step, so stop when just + * over one third of the digits are correct. + */ + int get_digits = static_cast(digits * 0.4); + std::uintmax_t maxit = 20; + auto f = [](const auto& g, const auto& x){ return g * g * g * g * g - x; }; + auto f_g = [](const auto& g, const auto& x){ return 5 * g * g * g * g; }; + auto f_gg = [](const auto& g, const auto& x){ return 20 * g * g * g; }; + auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& xx) { + using x_t = std::decay_t; + return root_finder(std::make_tuple(f, f_g, f_gg), + x_t(guess), + x_t(min), + x_t(max), xx); + }; + stan::test::expect_ad(full_f, x); +} diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index cec6e665aaf..d0bc09a7cc8 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -1,5 +1,4 @@ #include -#include #include @@ -10,20 +9,27 @@ TEST(RevFunctor, root_finder) { //using namespace std; // Help ADL of std functions. var x = 27; int exponent; - std::frexp(x.val(), &exponent); // Get exponent of z (ignore mantissa). - var guess = ldexp(1., exponent / 3); // Rough guess is to divide the exponent by three. - var min = ldexp(0.5, exponent / 3); // Minimum possible value is half our guess. - var max = ldexp(2., exponent / 3); // Maximum possible value is twice our guess. - const int digits = std::numeric_limits::digits; // Maximum possible binary digits accuracy for type T. - // digits used to control how accurate to try to make the result. - int get_digits = static_cast(digits * 0.4); // Accuracy triples with each step, so stop when just - // over one third of the digits are correct. + // Get exponent of z (ignore mantissa). + std::frexp(x.val(), &exponent); + // Rough guess is to divide the exponent by three. + var guess = ldexp(1., exponent / 3); + // Minimum possible value is half our guess. + var min = ldexp(0.5, exponent / 3); + // Maximum possible value is twice our guess. + var max = ldexp(2., exponent / 3); + // Maximum possible binary digits accuracy for type T. + const int digits = std::numeric_limits::digits; + /* digits used to control how accurate to try to make the result. + * Accuracy triples with each step, so stop when just + * over one third of the digits are correct. + */ + int get_digits = static_cast(digits * 0.4); std::uintmax_t maxit = 20; auto f = [](const auto& g, const auto x){ return g * g * g - x; }; auto f_g = [](const auto& g, const auto x){ return 3 * g * g; }; auto f_gg = [](const auto& g, const auto x){ return 6 * g; }; - var result = root_finder(f, f_g, f_gg, guess, min, max, x); + var result = root_finder(std::make_tuple(f, f_g, f_gg), guess, min, max, x); result.grad(); - // val: 27 adj:0.037037 - std::cout << "\nval: " << x.val() << " adj:" << x.adj() << "\n"; + EXPECT_EQ(27, x.val()); + EXPECT_NEAR(0.037037, x.adj(), 1e-6); } From 9b2b47713f3166c89afb899171e31d911e32fbac Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Fri, 29 Apr 2022 18:45:13 -0400 Subject: [PATCH 04/20] adds promotion logic to root_solver_tol to avoid boost errors --- stan/math/prim/functor/root_finder.hpp | 5 +++-- test/unit/math/mix/functor/root_finder_test.cpp | 12 ++---------- 2 files changed, 5 insertions(+), 12 deletions(-) diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index fef56d3139a..e405c35cb4d 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -53,10 +53,11 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, const MinScalar check_bounded("root_finder", "initial guess", guess, min, max); check_positive("root_finder", "digits", digits); check_positive("root_finder", "max_iter", max_iter); - return_type_t ret = 0; + using ret_t = return_type_t; + ret_t ret = 0; auto f_plus_div = internal::func_with_derivs(f_tuple, args...); try { - ret = boost::math::tools::halley_iterate(f_plus_div, guess, min, max, + ret = boost::math::tools::halley_iterate(f_plus_div, ret_t(guess), ret_t(min), ret_t(max), digits, max_iter); } catch (const std::exception& e) { throw e; diff --git a/test/unit/math/mix/functor/root_finder_test.cpp b/test/unit/math/mix/functor/root_finder_test.cpp index 454ec8af0a5..52fc5ea7e7f 100644 --- a/test/unit/math/mix/functor/root_finder_test.cpp +++ b/test/unit/math/mix/functor/root_finder_test.cpp @@ -26,11 +26,7 @@ TEST(MixFun, root_finder_cubed) { auto f_g = [](const auto& g, const auto& x){ return 3 * g * g; }; auto f_gg = [](const auto& g, const auto& x){ return 6 * g; }; auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& xx) { - using x_t = std::decay_t; - return root_finder(std::make_tuple(f, f_g, f_gg), - x_t(guess), - x_t(min), - x_t(max), xx); + return root_finder(std::make_tuple(f, f_g, f_gg), guess, min, max, xx); }; stan::test::expect_ad(full_f, x); } @@ -61,11 +57,7 @@ TEST(MixFun, root_finder_fifth) { auto f_g = [](const auto& g, const auto& x){ return 5 * g * g * g * g; }; auto f_gg = [](const auto& g, const auto& x){ return 20 * g * g * g; }; auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& xx) { - using x_t = std::decay_t; - return root_finder(std::make_tuple(f, f_g, f_gg), - x_t(guess), - x_t(min), - x_t(max), xx); + return root_finder(std::make_tuple(f, f_g, f_gg), guess, min, max, xx); }; stan::test::expect_ad(full_f, x); } From d72900171622017048c9416c28cb5eb7e7bac281 Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Fri, 29 Apr 2022 19:01:27 -0400 Subject: [PATCH 05/20] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/fwd/fun/frexp.hpp | 4 +- stan/math/fwd/fun/sign.hpp | 4 +- stan/math/prim/functor/root_finder.hpp | 10 ++-- stan/math/rev/fun/frexp.hpp | 4 +- stan/math/rev/fun/sign.hpp | 8 ++-- stan/math/rev/functor/root_finder.hpp | 48 +++++++++---------- .../math/mix/functor/root_finder_test.cpp | 16 +++---- .../math/rev/functor/root_finder_test.cpp | 9 ++-- 8 files changed, 50 insertions(+), 53 deletions(-) diff --git a/stan/math/fwd/fun/frexp.hpp b/stan/math/fwd/fun/frexp.hpp index 3d0b44b121b..aab3f108e28 100644 --- a/stan/math/fwd/fun/frexp.hpp +++ b/stan/math/fwd/fun/frexp.hpp @@ -12,6 +12,6 @@ template inline auto frexp(const fvar& x, int* exponent) noexcept { return std::frexp(value_of_rec(x), exponent); } -} -} +} // namespace math +} // namespace stan #endif diff --git a/stan/math/fwd/fun/sign.hpp b/stan/math/fwd/fun/sign.hpp index b36067b367d..a46ae60ba3c 100644 --- a/stan/math/fwd/fun/sign.hpp +++ b/stan/math/fwd/fun/sign.hpp @@ -14,6 +14,6 @@ inline auto sign(const fvar& x) { return (0. < x_val) - (x_val < 0.); } -} -} +} // namespace math +} // namespace stan #endif diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index e405c35cb4d..f59a18298a7 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -47,9 +47,9 @@ template * = nullptr> -auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, const MinScalar min, - const MaxScalar max, const int digits, std::uintmax_t& max_iter, - Types&&... args) { +auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, + const MinScalar min, const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); check_positive("root_finder", "digits", digits); check_positive("root_finder", "max_iter", max_iter); @@ -57,8 +57,8 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, const MinScalar ret_t ret = 0; auto f_plus_div = internal::func_with_derivs(f_tuple, args...); try { - ret = boost::math::tools::halley_iterate(f_plus_div, ret_t(guess), ret_t(min), ret_t(max), - digits, max_iter); + ret = boost::math::tools::halley_iterate( + f_plus_div, ret_t(guess), ret_t(min), ret_t(max), digits, max_iter); } catch (const std::exception& e) { throw e; } diff --git a/stan/math/rev/fun/frexp.hpp b/stan/math/rev/fun/frexp.hpp index ba5437d1f43..300e4a67914 100644 --- a/stan/math/rev/fun/frexp.hpp +++ b/stan/math/rev/fun/frexp.hpp @@ -9,6 +9,6 @@ namespace math { inline auto frexp(stan::math::var x, int* exponent) noexcept { return std::frexp(x.val(), exponent); } -} -} +} // namespace math +} // namespace stan #endif diff --git a/stan/math/rev/fun/sign.hpp b/stan/math/rev/fun/sign.hpp index 0c3bb816084..ad7e8ff3737 100644 --- a/stan/math/rev/fun/sign.hpp +++ b/stan/math/rev/fun/sign.hpp @@ -6,9 +6,7 @@ namespace stan { namespace math { -inline int sign(stan::math::var z) { - return (z == 0) ? 0 : z < 0 ? -1 : 1; -} -} -} +inline int sign(stan::math::var z) { return (z == 0) ? 0 : z < 0 ? -1 : 1; } +} // namespace math +} // namespace stan #endif diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp index 435c3ea13a1..936a4157155 100644 --- a/stan/math/rev/functor/root_finder.hpp +++ b/stan/math/rev/functor/root_finder.hpp @@ -37,9 +37,9 @@ template < typename MaxScalar, typename... Types, require_any_st_var* = nullptr, require_all_stan_scalar_t* = nullptr> -auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, const MinScalar min, - const MaxScalar max, const int digits, std::uintmax_t& max_iter, - Types&&... args) { +auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, + const MinScalar min, const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); check_positive("root_finder", "digits", digits); check_positive("root_finder", "max_iter", max_iter); @@ -71,28 +71,28 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, const MinScalar } /* - * Note: Because we put this on the callback stack, if `f` is a lambda - * its captures must be in Stan's arena memory or trivially destructable. - */ + * Note: Because we put this on the callback stack, if `f` is a lambda + * its captures must be in Stan's arena memory or trivially destructable. + */ return make_callback_var(theta_dbl, - [f, arena_args_tuple, Jf_x](auto& ret) mutable { - // Eigen::VectorXd eta = - // -Jf_x_T_lu_ptr->solve(ret.adj().eval()); - double eta = -(ret.adj() / Jf_x); - // Contract with Jacobian of f with respect to y - // using a nested reverse autodiff pass. - { - nested_rev_autodiff rev; - double ret_val = ret.val(); - auto x_nrad_ = apply( - [&ret_val, &f](const auto&... args) { - return eval(f(ret_val, args...)); - }, - *arena_args_tuple); - x_nrad_.adj() = eta; - grad(); - } - }); + [f, arena_args_tuple, Jf_x](auto& ret) mutable { + // Eigen::VectorXd eta = + // -Jf_x_T_lu_ptr->solve(ret.adj().eval()); + double eta = -(ret.adj() / Jf_x); + // Contract with Jacobian of f with respect to y + // using a nested reverse autodiff pass. + { + nested_rev_autodiff rev; + double ret_val = ret.val(); + auto x_nrad_ = apply( + [&ret_val, &f](const auto&... args) { + return eval(f(ret_val, args...)); + }, + *arena_args_tuple); + x_nrad_.adj() = eta; + grad(); + } + }); } } // namespace math diff --git a/test/unit/math/mix/functor/root_finder_test.cpp b/test/unit/math/mix/functor/root_finder_test.cpp index 52fc5ea7e7f..1236fd8aee8 100644 --- a/test/unit/math/mix/functor/root_finder_test.cpp +++ b/test/unit/math/mix/functor/root_finder_test.cpp @@ -3,7 +3,7 @@ TEST(MixFun, root_finder_cubed) { using stan::math::root_finder; // return cube root of x using 1st and 2nd derivatives and Halley. - //using namespace std; // Help ADL of std functions. + // using namespace std; // Help ADL of std functions. double x = 27; int exponent; // Get exponent of z (ignore mantissa). @@ -22,9 +22,9 @@ TEST(MixFun, root_finder_cubed) { */ int get_digits = static_cast(digits * 0.4); std::uintmax_t maxit = 20; - auto f = [](const auto& g, const auto& x){ return g * g * g - x; }; - auto f_g = [](const auto& g, const auto& x){ return 3 * g * g; }; - auto f_gg = [](const auto& g, const auto& x){ return 6 * g; }; + auto f = [](const auto& g, const auto& x) { return g * g * g - x; }; + auto f_g = [](const auto& g, const auto& x) { return 3 * g * g; }; + auto f_gg = [](const auto& g, const auto& x) { return 6 * g; }; auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& xx) { return root_finder(std::make_tuple(f, f_g, f_gg), guess, min, max, xx); }; @@ -34,7 +34,7 @@ TEST(MixFun, root_finder_cubed) { TEST(MixFun, root_finder_fifth) { using stan::math::root_finder; // return cube root of x using 1st and 2nd derivatives and Halley. - //using namespace std; // Help ADL of std functions. + // using namespace std; // Help ADL of std functions. double x = 27; int exponent; // Get exponent of z (ignore mantissa). @@ -53,9 +53,9 @@ TEST(MixFun, root_finder_fifth) { */ int get_digits = static_cast(digits * 0.4); std::uintmax_t maxit = 20; - auto f = [](const auto& g, const auto& x){ return g * g * g * g * g - x; }; - auto f_g = [](const auto& g, const auto& x){ return 5 * g * g * g * g; }; - auto f_gg = [](const auto& g, const auto& x){ return 20 * g * g * g; }; + auto f = [](const auto& g, const auto& x) { return g * g * g * g * g - x; }; + auto f_g = [](const auto& g, const auto& x) { return 5 * g * g * g * g; }; + auto f_gg = [](const auto& g, const auto& x) { return 20 * g * g * g; }; auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& xx) { return root_finder(std::make_tuple(f, f_g, f_gg), guess, min, max, xx); }; diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index d0bc09a7cc8..9e8def06604 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -1,12 +1,11 @@ #include #include - TEST(RevFunctor, root_finder) { using stan::math::root_finder; using stan::math::var; // return cube root of x using 1st and 2nd derivatives and Halley. - //using namespace std; // Help ADL of std functions. + // using namespace std; // Help ADL of std functions. var x = 27; int exponent; // Get exponent of z (ignore mantissa). @@ -25,9 +24,9 @@ TEST(RevFunctor, root_finder) { */ int get_digits = static_cast(digits * 0.4); std::uintmax_t maxit = 20; - auto f = [](const auto& g, const auto x){ return g * g * g - x; }; - auto f_g = [](const auto& g, const auto x){ return 3 * g * g; }; - auto f_gg = [](const auto& g, const auto x){ return 6 * g; }; + auto f = [](const auto& g, const auto x) { return g * g * g - x; }; + auto f_g = [](const auto& g, const auto x) { return 3 * g * g; }; + auto f_gg = [](const auto& g, const auto x) { return 6 * g; }; var result = root_finder(std::make_tuple(f, f_g, f_gg), guess, min, max, x); result.grad(); EXPECT_EQ(27, x.val()); From f2f8c51aaca5b982799048448f6db15e226892ac Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Fri, 29 Apr 2022 19:02:28 -0400 Subject: [PATCH 06/20] fix header includes --- stan/math/prim/functor/root_finder.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index f59a18298a7..9cb79a9466e 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include From 04dbda898386e8201241da86e1450aaaae6fe11a Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Fri, 29 Apr 2022 19:19:36 -0400 Subject: [PATCH 07/20] fix header includes --- stan/math/prim/functor/root_finder.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index 9cb79a9466e..2024e5568d1 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include From c0d4a819526403d8ce7c0d8bbe47d9572ea574e4 Mon Sep 17 00:00:00 2001 From: stevebronder Date: Tue, 3 May 2022 16:47:41 -0400 Subject: [PATCH 08/20] start working on beta test --- stan/math/prim/functor/root_finder.hpp | 1 + stan/math/rev/functor/root_finder.hpp | 7 +- .../math/mix/functor/root_finder_test.cpp | 63 ++++++++++++++-- .../math/rev/functor/root_finder_test.cpp | 75 +++++++++++++++++-- 4 files changed, 130 insertions(+), 16 deletions(-) diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index 2024e5568d1..23f5fcdf3fe 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -62,6 +62,7 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, ret = boost::math::tools::halley_iterate( f_plus_div, ret_t(guess), ret_t(min), ret_t(max), digits, max_iter); } catch (const std::exception& e) { + std::cout << "err: \n" << e.what() << "\n"; throw e; } return ret; diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp index 936a4157155..a0c67cbdc0a 100644 --- a/stan/math/rev/functor/root_finder.hpp +++ b/stan/math/rev/functor/root_finder.hpp @@ -2,9 +2,10 @@ #define STAN_MATH_REV_FUNCTOR_ROOT_FINDER_HPP #include -#include #include #include +#include +#include #include #include #include @@ -51,9 +52,9 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, *arena_args_tuple); // Solve the system double theta_dbl = apply( - [&f_tuple, guess_val = value_of(guess), min_val = value_of(min), + [&f_tuple, &max_iter, digits, guess_val = value_of(guess), min_val = value_of(min), max_val = value_of(max)](auto&&... vals) { - return root_finder(f_tuple, guess_val, min_val, max_val, vals...); + return root_finder_tol(f_tuple, guess_val, min_val, max_val, digits, max_iter, vals...); }, args_vals_tuple); double Jf_x; diff --git a/test/unit/math/mix/functor/root_finder_test.cpp b/test/unit/math/mix/functor/root_finder_test.cpp index 1236fd8aee8..a97865f7a61 100644 --- a/test/unit/math/mix/functor/root_finder_test.cpp +++ b/test/unit/math/mix/functor/root_finder_test.cpp @@ -1,5 +1,6 @@ #include +/* TEST(MixFun, root_finder_cubed) { using stan::math::root_finder; // return cube root of x using 1st and 2nd derivatives and Halley. @@ -16,10 +17,6 @@ TEST(MixFun, root_finder_cubed) { double max = ldexp(2., exponent / 3); // Maximum possible binary digits accuracy for type T. const int digits = std::numeric_limits::digits; - /* digits used to control how accurate to try to make the result. - * Accuracy triples with each step, so stop when just - * over one third of the digits are correct. - */ int get_digits = static_cast(digits * 0.4); std::uintmax_t maxit = 20; auto f = [](const auto& g, const auto& x) { return g * g * g - x; }; @@ -47,10 +44,6 @@ TEST(MixFun, root_finder_fifth) { double max = ldexp(2., exponent / 5); // Maximum possible binary digits accuracy for type T. const int digits = std::numeric_limits::digits; - /* digits used to control how accurate to try to make the result. - * Accuracy triples with each step, so stop when just - * over one third of the digits are correct. - */ int get_digits = static_cast(digits * 0.4); std::uintmax_t maxit = 20; auto f = [](const auto& g, const auto& x) { return g * g * g * g * g - x; }; @@ -61,3 +54,57 @@ TEST(MixFun, root_finder_fifth) { }; stan::test::expect_ad(full_f, x); } +*/ +TEST(MixFun, root_finder_beta_log) { + auto beta_lpdf_deriv = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + using stan::math::log; + using stan::math::lgamma; + using stan::math::log1m; + auto val = -((beta * log1m(x) + (-2 + alpha) * log(x) + log(1 - 2*x - alpha + x * alpha + x * beta) + lgamma(alpha + beta)) + - (2 * log1m(x) + lgamma(alpha) + lgamma(beta))); + return val; + }; + auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) {return stan::math::beta_lcdf(x, alpha, beta) - p;}; + auto f_g = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) {return stan::math::beta_lpdf(x, alpha, beta);}; + auto f_gg = [&beta_lpdf_deriv](auto&& x, auto&& alpha, auto&& beta, auto&& p) {return beta_lpdf_deriv(x, alpha, beta, p);}; + double guess = .3; + double min = 0.0000000001; + double max = 1; + auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& alpha, auto&& beta, auto&& p) { + std::uintmax_t max_its=1000; + return stan::math::root_finder_tol(std::make_tuple(f, f_g, f_gg), guess, min, max, 16, max_its, alpha, beta, p); + }; + double p = 0.45; + double alpha = .5; + double beta = .4; + stan::test::expect_ad(full_f, alpha, beta, p); +} + +TEST(MixFun, root_finder_beta) { + auto beta_pdf_deriv = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + using stan::math::lgamma; + using stan::math::pow; + auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) * (1. - 2. * x - alpha + x * alpha + x * beta) * lgamma(alpha + beta)) / + (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); + return val; + }; + auto beta_pdf = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + using stan::math::lgamma; + using stan::math::pow; + auto val = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta))/ ((lgamma(alpha) * lgamma(beta))); + return val; + }; + + auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) {return stan::math::beta_cdf(x, alpha, beta) - p;}; + double guess = .7; + double min = 0.0000000001; + double max = 1; + auto full_f = [&f, &beta_pdf, &beta_pdf_deriv, guess, min, max](auto&& alpha, auto&& beta, auto&& p) { + std::uintmax_t max_its=1000; + return stan::math::root_finder_tol(std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, 16, max_its, alpha, beta, p); + }; + double p = 0.8; + double alpha = .4; + double beta = .5; + stan::test::expect_ad(full_f, alpha, beta, p); +} diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index 9e8def06604..fd27f22f141 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -1,6 +1,7 @@ #include #include - +#include +/* TEST(RevFunctor, root_finder) { using stan::math::root_finder; using stan::math::var; @@ -18,10 +19,6 @@ TEST(RevFunctor, root_finder) { var max = ldexp(2., exponent / 3); // Maximum possible binary digits accuracy for type T. const int digits = std::numeric_limits::digits; - /* digits used to control how accurate to try to make the result. - * Accuracy triples with each step, so stop when just - * over one third of the digits are correct. - */ int get_digits = static_cast(digits * 0.4); std::uintmax_t maxit = 20; auto f = [](const auto& g, const auto x) { return g * g * g - x; }; @@ -32,3 +29,71 @@ TEST(RevFunctor, root_finder) { EXPECT_EQ(27, x.val()); EXPECT_NEAR(0.037037, x.adj(), 1e-6); } +*/ +TEST(RevFunctor, root_finder_beta_cdf) { + using stan::math::var; + auto func = [](auto&& vals) { + auto p = vals(0); + auto alpha = vals(1); + auto beta = vals(2); + boost::math::beta_distribution my_beta(alpha, beta); + return boost::math::quantile(my_beta, p); + }; + Eigen::VectorXd vals(3); + // p, alpha, beta + vals << 0.4, 1.2, 1.9; + double fx = 0; + Eigen::VectorXd finit_grad_fx(3); + stan::math::finite_diff_gradient(func, vals, fx, finit_grad_fx, 1e-14); + std::cout << "--- Finit Diff----\n"; + std::cout << "fx: " << fx; + std::cout << "\ngrads: \n" << + "p: " << finit_grad_fx(0) << "\n" + "alpha: " << finit_grad_fx(1) << "\n" + "beta: " << finit_grad_fx(2) << "\n"; + + auto beta_pdf_deriv = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + using stan::math::lgamma; + using stan::math::pow; + auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) * (1. - 2. * x - alpha + x * alpha + x * beta) * lgamma(alpha + beta)) / + (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); + return val; + }; + auto beta_pdf = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + using stan::math::lgamma; + using stan::math::pow; + auto val = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta))/ ((lgamma(alpha) * lgamma(beta))); + return val; + }; + + auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + return stan::math::beta_cdf(x, alpha, beta) - p; + }; + double guess = .3; + double min = 0; + double max = 1; + auto full_f = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max](auto&& alpha, auto&& beta, auto&& p) { + std::uintmax_t max_its=1000; + return stan::math::root_finder_tol(std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, 16, max_its, alpha, beta, p); + }; + auto func2 = [&full_f](auto&& vals) { + auto p = vals(0); + auto alpha = vals(1); + auto beta = vals(2); + return full_f(alpha, beta, p); + }; + Eigen::VectorXd grad_fx(3); + stan::math::gradient(func2, vals, fx, grad_fx); + std::cout << "--- Auto Diff----\n"; + std::cout << "fx: " << fx; + std::cout << "\ngrads: \n" << + "p: " << grad_fx(0) << "\n" + "alpha: " << grad_fx(1) << "\n" + "beta: " << grad_fx(2) << "\n"; + + Eigen::VectorXd diff_grad_fx = finit_grad_fx - grad_fx; + std::cout << "\ngrad diffs: \n" << + "p: " << diff_grad_fx(0) << "\n" + "alpha: " << diff_grad_fx(1) << "\n" + "beta: " << diff_grad_fx(2) << "\n"; +} From b08f3adfd8ab2606a5a8bff5f6bbf7819b72bcaa Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Tue, 3 May 2022 16:48:37 -0400 Subject: [PATCH 09/20] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/rev/functor/root_finder.hpp | 7 ++- .../math/mix/functor/root_finder_test.cpp | 52 +++++++++++----- .../math/rev/functor/root_finder_test.cpp | 59 ++++++++++++------- 3 files changed, 79 insertions(+), 39 deletions(-) diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp index a0c67cbdc0a..04e94b2ebe6 100644 --- a/stan/math/rev/functor/root_finder.hpp +++ b/stan/math/rev/functor/root_finder.hpp @@ -52,9 +52,10 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, *arena_args_tuple); // Solve the system double theta_dbl = apply( - [&f_tuple, &max_iter, digits, guess_val = value_of(guess), min_val = value_of(min), - max_val = value_of(max)](auto&&... vals) { - return root_finder_tol(f_tuple, guess_val, min_val, max_val, digits, max_iter, vals...); + [&f_tuple, &max_iter, digits, guess_val = value_of(guess), + min_val = value_of(min), max_val = value_of(max)](auto&&... vals) { + return root_finder_tol(f_tuple, guess_val, min_val, max_val, digits, + max_iter, vals...); }, args_vals_tuple); double Jf_x; diff --git a/test/unit/math/mix/functor/root_finder_test.cpp b/test/unit/math/mix/functor/root_finder_test.cpp index a97865f7a61..ff6de352e68 100644 --- a/test/unit/math/mix/functor/root_finder_test.cpp +++ b/test/unit/math/mix/functor/root_finder_test.cpp @@ -60,19 +60,30 @@ TEST(MixFun, root_finder_beta_log) { using stan::math::log; using stan::math::lgamma; using stan::math::log1m; - auto val = -((beta * log1m(x) + (-2 + alpha) * log(x) + log(1 - 2*x - alpha + x * alpha + x * beta) + lgamma(alpha + beta)) - - (2 * log1m(x) + lgamma(alpha) + lgamma(beta))); + auto val = -((beta * log1m(x) + (-2 + alpha) * log(x) + + log(1 - 2 * x - alpha + x * alpha + x * beta) + + lgamma(alpha + beta)) + - (2 * log1m(x) + lgamma(alpha) + lgamma(beta))); return val; }; - auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) {return stan::math::beta_lcdf(x, alpha, beta) - p;}; - auto f_g = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) {return stan::math::beta_lpdf(x, alpha, beta);}; - auto f_gg = [&beta_lpdf_deriv](auto&& x, auto&& alpha, auto&& beta, auto&& p) {return beta_lpdf_deriv(x, alpha, beta, p);}; + auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + return stan::math::beta_lcdf(x, alpha, beta) - p; + }; + auto f_g = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + return stan::math::beta_lpdf(x, alpha, beta); + }; + auto f_gg + = [&beta_lpdf_deriv](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + return beta_lpdf_deriv(x, alpha, beta, p); + }; double guess = .3; double min = 0.0000000001; double max = 1; - auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& alpha, auto&& beta, auto&& p) { - std::uintmax_t max_its=1000; - return stan::math::root_finder_tol(std::make_tuple(f, f_g, f_gg), guess, min, max, 16, max_its, alpha, beta, p); + auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& alpha, auto&& beta, + auto&& p) { + std::uintmax_t max_its = 1000; + return stan::math::root_finder_tol(std::make_tuple(f, f_g, f_gg), guess, + min, max, 16, max_its, alpha, beta, p); }; double p = 0.45; double alpha = .5; @@ -84,24 +95,33 @@ TEST(MixFun, root_finder_beta) { auto beta_pdf_deriv = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { using stan::math::lgamma; using stan::math::pow; - auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) * (1. - 2. * x - alpha + x * alpha + x * beta) * lgamma(alpha + beta)) / - (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); - return val; + auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) + * (1. - 2. * x - alpha + x * alpha + x * beta) + * lgamma(alpha + beta)) + / (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); + return val; }; auto beta_pdf = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { using stan::math::lgamma; using stan::math::pow; - auto val = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta))/ ((lgamma(alpha) * lgamma(beta))); + auto val + = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta)) + / ((lgamma(alpha) * lgamma(beta))); return val; }; - auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) {return stan::math::beta_cdf(x, alpha, beta) - p;}; + auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + return stan::math::beta_cdf(x, alpha, beta) - p; + }; double guess = .7; double min = 0.0000000001; double max = 1; - auto full_f = [&f, &beta_pdf, &beta_pdf_deriv, guess, min, max](auto&& alpha, auto&& beta, auto&& p) { - std::uintmax_t max_its=1000; - return stan::math::root_finder_tol(std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, 16, max_its, alpha, beta, p); + auto full_f = [&f, &beta_pdf, &beta_pdf_deriv, guess, min, max]( + auto&& alpha, auto&& beta, auto&& p) { + std::uintmax_t max_its = 1000; + return stan::math::root_finder_tol( + std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, 16, + max_its, alpha, beta, p); }; double p = 0.8; double alpha = .4; diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index fd27f22f141..7f1190da882 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -47,22 +47,30 @@ TEST(RevFunctor, root_finder_beta_cdf) { stan::math::finite_diff_gradient(func, vals, fx, finit_grad_fx, 1e-14); std::cout << "--- Finit Diff----\n"; std::cout << "fx: " << fx; - std::cout << "\ngrads: \n" << - "p: " << finit_grad_fx(0) << "\n" - "alpha: " << finit_grad_fx(1) << "\n" - "beta: " << finit_grad_fx(2) << "\n"; + std::cout << "\ngrads: \n" + << "p: " << finit_grad_fx(0) + << "\n" + "alpha: " + << finit_grad_fx(1) + << "\n" + "beta: " + << finit_grad_fx(2) << "\n"; auto beta_pdf_deriv = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { using stan::math::lgamma; using stan::math::pow; - auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) * (1. - 2. * x - alpha + x * alpha + x * beta) * lgamma(alpha + beta)) / - (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); - return val; + auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) + * (1. - 2. * x - alpha + x * alpha + x * beta) + * lgamma(alpha + beta)) + / (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); + return val; }; auto beta_pdf = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { using stan::math::lgamma; using stan::math::pow; - auto val = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta))/ ((lgamma(alpha) * lgamma(beta))); + auto val + = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta)) + / ((lgamma(alpha) * lgamma(beta))); return val; }; @@ -72,9 +80,12 @@ TEST(RevFunctor, root_finder_beta_cdf) { double guess = .3; double min = 0; double max = 1; - auto full_f = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max](auto&& alpha, auto&& beta, auto&& p) { - std::uintmax_t max_its=1000; - return stan::math::root_finder_tol(std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, 16, max_its, alpha, beta, p); + auto full_f = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( + auto&& alpha, auto&& beta, auto&& p) { + std::uintmax_t max_its = 1000; + return stan::math::root_finder_tol( + std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, 16, + max_its, alpha, beta, p); }; auto func2 = [&full_f](auto&& vals) { auto p = vals(0); @@ -86,14 +97,22 @@ TEST(RevFunctor, root_finder_beta_cdf) { stan::math::gradient(func2, vals, fx, grad_fx); std::cout << "--- Auto Diff----\n"; std::cout << "fx: " << fx; - std::cout << "\ngrads: \n" << - "p: " << grad_fx(0) << "\n" - "alpha: " << grad_fx(1) << "\n" - "beta: " << grad_fx(2) << "\n"; + std::cout << "\ngrads: \n" + << "p: " << grad_fx(0) + << "\n" + "alpha: " + << grad_fx(1) + << "\n" + "beta: " + << grad_fx(2) << "\n"; - Eigen::VectorXd diff_grad_fx = finit_grad_fx - grad_fx; - std::cout << "\ngrad diffs: \n" << - "p: " << diff_grad_fx(0) << "\n" - "alpha: " << diff_grad_fx(1) << "\n" - "beta: " << diff_grad_fx(2) << "\n"; + Eigen::VectorXd diff_grad_fx = finit_grad_fx - grad_fx; + std::cout << "\ngrad diffs: \n" + << "p: " << diff_grad_fx(0) + << "\n" + "alpha: " + << diff_grad_fx(1) + << "\n" + "beta: " + << diff_grad_fx(2) << "\n"; } From 33e0b505a379e7e66f58c1f6dc8068a35fb50274 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Tue, 10 May 2022 18:03:45 -0400 Subject: [PATCH 10/20] update rev test --- stan/math/prim/functor/root_finder.hpp | 96 +++++++++++++++++-- stan/math/rev/functor/root_finder.hpp | 15 ++- .../math/mix/functor/root_finder_test.cpp | 9 +- .../math/rev/functor/root_finder_test.cpp | 29 ++++-- 4 files changed, 117 insertions(+), 32 deletions(-) diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index 23f5fcdf3fe..41ad22490e4 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -45,12 +45,13 @@ inline auto func_with_derivs(Tuple&& f_tuple, Args&&... args) { * this is updated to the actual number of iterations performed * @param args Parameter pack of arguments to pass the the functors in `f_tuple` */ -template * = nullptr> -auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, - const MinScalar min, const MaxScalar max, const int digits, +auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, + const GuessScalar guess, const MinScalar min, + const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); check_positive("root_finder", "digits", digits); @@ -59,8 +60,8 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, ret_t ret = 0; auto f_plus_div = internal::func_with_derivs(f_tuple, args...); try { - ret = boost::math::tools::halley_iterate( - f_plus_div, ret_t(guess), ret_t(min), ret_t(max), digits, max_iter); + ret = f_solver(f_plus_div, ret_t(guess), ret_t(min), ret_t(max), digits, + max_iter); } catch (const std::exception& e) { std::cout << "err: \n" << e.what() << "\n"; throw e; @@ -68,6 +69,45 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, return ret; } +template +auto root_finder_halley_tol(FTuple&& f_tuple, + const GuessScalar guess, const MinScalar min, + const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types&&... args) { + return root_finder_tol( + [](auto&&... args) { + return boost::math::tools::halley_iterate(args...); + }, std::forward(f_tuple), guess, + min, max, digits, max_iter, std::forward(args)...); +} + +template +auto root_finder_newton_raphson_tol(FTuple&& f_tuple, + const GuessScalar guess, const MinScalar min, + const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types&&... args) { + return root_finder_tol( + [](auto&&... args) { + return boost::math::tools::newton_raphson_iterate(args...); + }, std::forward(f_tuple), guess, + min, max, digits, max_iter, std::forward(args)...); +} + +template +auto root_finder_schroder_tol(FTuple&& f_tuple, + const GuessScalar guess, const MinScalar min, + const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types&&... args) { + return root_finder_tol( + [](auto&&... args) { + return boost::math::tools::schroder_iterate(args...); + }, std::forward(f_tuple), guess, + min, max, digits, max_iter, std::forward(args)...); +} + /** * Solve for root using Boost Halley method with default values for the * tolerances @@ -86,15 +126,51 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, * initial upper bracket * @param args Parameter pack of arguments to pass the the functors in `f_tuple` */ -template -auto root_finder(FTuple&& f_tuple, const GuessScalar guess, const MinScalar min, +template +auto root_finder(SolverFun&& f_solver, FTuple&& f_tuple, + const GuessScalar guess, const MinScalar min, const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); - return root_finder_tol(std::forward(f_tuple), guess, min, max, digits, + return root_finder_tol(std::forward(f_solver), + std::forward(f_tuple), guess, min, max, digits, max_iter, std::forward(args)...); } + +template +auto root_finder_hailey(FTuple&& f_tuple, const GuessScalar guess, + const MinScalar min, const MaxScalar max, + Types&&... args) { + constexpr int digits = 16; + std::uintmax_t max_iter = std::numeric_limits::max(); + return root_finder_halley_tol(std::forward(f_tuple), guess, + min, max, digits, max_iter, std::forward(args)...); +} + +template +auto root_finder_newton_raphson(FTuple&& f_tuple, const GuessScalar guess, + const MinScalar min, const MaxScalar max, + Types&&... args) { + constexpr int digits = 16; + std::uintmax_t max_iter = std::numeric_limits::max(); + return root_finder_newton_raphson_tol(std::forward(f_tuple), guess, + min, max, digits, max_iter, std::forward(args)...); +} + +template +auto root_finder_schroder(FTuple&& f_tuple, const GuessScalar guess, + const MinScalar min, const MaxScalar max, + Types&&... args) { + constexpr int digits = 16; + std::uintmax_t max_iter = std::numeric_limits::max(); + return root_finder_schroder_tol(std::forward(f_tuple), guess, + min, max, digits, max_iter, std::forward(args)...); +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp index 04e94b2ebe6..d5ac07c09ad 100644 --- a/stan/math/rev/functor/root_finder.hpp +++ b/stan/math/rev/functor/root_finder.hpp @@ -33,18 +33,18 @@ namespace math { * this is updated to the actual number of iterations performed * @param args Parameter pack of arguments to pass the the functors in `f_tuple` */ -template < +template * = nullptr, require_all_stan_scalar_t* = nullptr> -auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, +auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, const GuessScalar guess, const MinScalar min, const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); check_positive("root_finder", "digits", digits); check_positive("root_finder", "max_iter", max_iter); - auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(args)...)); + auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(std::forward(args))...)); auto args_vals_tuple = apply( [&](const auto&... args) { return std::make_tuple(to_ref(value_of(args))...); @@ -52,14 +52,13 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, *arena_args_tuple); // Solve the system double theta_dbl = apply( - [&f_tuple, &max_iter, digits, guess_val = value_of(guess), + [&f_solver, &f_tuple, &max_iter, digits, guess_val = value_of(guess), min_val = value_of(min), max_val = value_of(max)](auto&&... vals) { - return root_finder_tol(f_tuple, guess_val, min_val, max_val, digits, + return root_finder_tol(f_solver, f_tuple, guess_val, min_val, max_val, digits, max_iter, vals...); }, args_vals_tuple); double Jf_x; - double f_x; auto&& f = std::get<0>(f_tuple); { nested_rev_autodiff nested; @@ -67,7 +66,6 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, stan::math::var fx_var = apply( [&x_var, &f](auto&&... args) { return f(x_var, std::move(args)...); }, std::move(args_vals_tuple)); - f_x = fx_var.val(); fx_var.grad(); Jf_x = x_var.adj(); } @@ -80,7 +78,6 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, [f, arena_args_tuple, Jf_x](auto& ret) mutable { // Eigen::VectorXd eta = // -Jf_x_T_lu_ptr->solve(ret.adj().eval()); - double eta = -(ret.adj() / Jf_x); // Contract with Jacobian of f with respect to y // using a nested reverse autodiff pass. { @@ -91,7 +88,7 @@ auto root_finder_tol(FTuple&& f_tuple, const GuessScalar guess, return eval(f(ret_val, args...)); }, *arena_args_tuple); - x_nrad_.adj() = eta; + x_nrad_.adj() = -(ret.adj() / Jf_x); grad(); } }); diff --git a/test/unit/math/mix/functor/root_finder_test.cpp b/test/unit/math/mix/functor/root_finder_test.cpp index ff6de352e68..f37e1e5a722 100644 --- a/test/unit/math/mix/functor/root_finder_test.cpp +++ b/test/unit/math/mix/functor/root_finder_test.cpp @@ -82,8 +82,8 @@ TEST(MixFun, root_finder_beta_log) { auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& alpha, auto&& beta, auto&& p) { std::uintmax_t max_its = 1000; - return stan::math::root_finder_tol(std::make_tuple(f, f_g, f_gg), guess, - min, max, 16, max_its, alpha, beta, p); + return stan::math::root_finder_hailey(std::make_tuple(f, f_g, f_gg), guess, + min, max, alpha, beta, p); }; double p = 0.45; double alpha = .5; @@ -119,9 +119,8 @@ TEST(MixFun, root_finder_beta) { auto full_f = [&f, &beta_pdf, &beta_pdf_deriv, guess, min, max]( auto&& alpha, auto&& beta, auto&& p) { std::uintmax_t max_its = 1000; - return stan::math::root_finder_tol( - std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, 16, - max_its, alpha, beta, p); + return stan::math::root_finder_hailey( + std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, beta, p); }; double p = 0.8; double alpha = .4; diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index 7f1190da882..45af7104c4e 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -41,7 +41,7 @@ TEST(RevFunctor, root_finder_beta_cdf) { }; Eigen::VectorXd vals(3); // p, alpha, beta - vals << 0.4, 1.2, 1.9; + vals << 0.4, .5, .5; double fx = 0; Eigen::VectorXd finit_grad_fx(3); stan::math::finite_diff_gradient(func, vals, fx, finit_grad_fx, 1e-14); @@ -83,9 +83,9 @@ TEST(RevFunctor, root_finder_beta_cdf) { auto full_f = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( auto&& alpha, auto&& beta, auto&& p) { std::uintmax_t max_its = 1000; - return stan::math::root_finder_tol( - std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, 16, - max_its, alpha, beta, p); + return stan::math::root_finder_hailey( + std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, + alpha, beta, p); }; auto func2 = [&full_f](auto&& vals) { auto p = vals(0); @@ -94,7 +94,13 @@ TEST(RevFunctor, root_finder_beta_cdf) { return full_f(alpha, beta, p); }; Eigen::VectorXd grad_fx(3); - stan::math::gradient(func2, vals, fx, grad_fx); + Eigen::Matrix var_vec(vals); + stan::math::var fxvar = func2(var_vec); + fxvar.grad(); + grad_fx = var_vec.adj(); + fx = fxvar.val(); + std::cout << "fxvar adj:" << fxvar.adj() << "\n"; + //stan::math::gradient(func2, vals, fx, grad_fx); std::cout << "--- Auto Diff----\n"; std::cout << "fx: " << fx; std::cout << "\ngrads: \n" @@ -105,14 +111,21 @@ TEST(RevFunctor, root_finder_beta_cdf) { << "\n" "beta: " << grad_fx(2) << "\n"; - Eigen::VectorXd diff_grad_fx = finit_grad_fx - grad_fx; - std::cout << "\ngrad diffs: \n" - << "p: " << diff_grad_fx(0) + std::cout << "--- grad diffs----\n"; + std::cout << "p: " << diff_grad_fx(0) << "\n" "alpha: " << diff_grad_fx(1) << "\n" "beta: " << diff_grad_fx(2) << "\n"; + + double known_p_grad = 1.493914820513846; + double known_alpha_grad = 0.948540678312004; + double known_beta_grad = -0.7464041618483364; + std::cout << "--- Auto Diff----\n"; + std::cout << "p: " << grad_fx(0) - known_p_grad << "\n" << + "alpha: " << grad_fx(1) - known_alpha_grad << "\n" << + "beta: " << grad_fx(2) - known_beta_grad << "\n"; } From 78f84c5b55cfdefaec8bb536205ea5f8e5064335 Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Tue, 10 May 2022 18:04:56 -0400 Subject: [PATCH 11/20] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/functor/root_finder.hpp | 89 ++++++++++--------- stan/math/rev/functor/root_finder.hpp | 18 ++-- .../math/mix/functor/root_finder_test.cpp | 5 +- .../math/rev/functor/root_finder_test.cpp | 12 +-- 4 files changed, 66 insertions(+), 58 deletions(-) diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index 41ad22490e4..10a95a5b96c 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -69,43 +69,46 @@ auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, return ret; } -template -auto root_finder_halley_tol(FTuple&& f_tuple, - const GuessScalar guess, const MinScalar min, - const MaxScalar max, const int digits, - std::uintmax_t& max_iter, Types&&... args) { - return root_finder_tol( - [](auto&&... args) { - return boost::math::tools::halley_iterate(args...); - }, std::forward(f_tuple), guess, - min, max, digits, max_iter, std::forward(args)...); +template +auto root_finder_halley_tol(FTuple&& f_tuple, const GuessScalar guess, + const MinScalar min, const MaxScalar max, + const int digits, std::uintmax_t& max_iter, + Types&&... args) { + return root_finder_tol( + [](auto&&... args) { + return boost::math::tools::halley_iterate(args...); + }, + std::forward(f_tuple), guess, min, max, digits, max_iter, + std::forward(args)...); } -template -auto root_finder_newton_raphson_tol(FTuple&& f_tuple, - const GuessScalar guess, const MinScalar min, - const MaxScalar max, const int digits, - std::uintmax_t& max_iter, Types&&... args) { - return root_finder_tol( - [](auto&&... args) { - return boost::math::tools::newton_raphson_iterate(args...); - }, std::forward(f_tuple), guess, - min, max, digits, max_iter, std::forward(args)...); +template +auto root_finder_newton_raphson_tol(FTuple&& f_tuple, const GuessScalar guess, + const MinScalar min, const MaxScalar max, + const int digits, std::uintmax_t& max_iter, + Types&&... args) { + return root_finder_tol( + [](auto&&... args) { + return boost::math::tools::newton_raphson_iterate(args...); + }, + std::forward(f_tuple), guess, min, max, digits, max_iter, + std::forward(args)...); } -template -auto root_finder_schroder_tol(FTuple&& f_tuple, - const GuessScalar guess, const MinScalar min, - const MaxScalar max, const int digits, - std::uintmax_t& max_iter, Types&&... args) { - return root_finder_tol( - [](auto&&... args) { - return boost::math::tools::schroder_iterate(args...); - }, std::forward(f_tuple), guess, - min, max, digits, max_iter, std::forward(args)...); +template +auto root_finder_schroder_tol(FTuple&& f_tuple, const GuessScalar guess, + const MinScalar min, const MaxScalar max, + const int digits, std::uintmax_t& max_iter, + Types&&... args) { + return root_finder_tol( + [](auto&&... args) { + return boost::math::tools::schroder_iterate(args...); + }, + std::forward(f_tuple), guess, min, max, digits, max_iter, + std::forward(args)...); } /** @@ -145,30 +148,32 @@ auto root_finder_hailey(FTuple&& f_tuple, const GuessScalar guess, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); - return root_finder_halley_tol(std::forward(f_tuple), guess, - min, max, digits, max_iter, std::forward(args)...); + return root_finder_halley_tol(std::forward(f_tuple), guess, min, max, + digits, max_iter, std::forward(args)...); } template auto root_finder_newton_raphson(FTuple&& f_tuple, const GuessScalar guess, - const MinScalar min, const MaxScalar max, - Types&&... args) { + const MinScalar min, const MaxScalar max, + Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); return root_finder_newton_raphson_tol(std::forward(f_tuple), guess, - min, max, digits, max_iter, std::forward(args)...); + min, max, digits, max_iter, + std::forward(args)...); } template auto root_finder_schroder(FTuple&& f_tuple, const GuessScalar guess, - const MinScalar min, const MaxScalar max, - Types&&... args) { + const MinScalar min, const MaxScalar max, + Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); - return root_finder_schroder_tol(std::forward(f_tuple), guess, - min, max, digits, max_iter, std::forward(args)...); + return root_finder_schroder_tol(std::forward(f_tuple), guess, min, + max, digits, max_iter, + std::forward(args)...); } } // namespace math diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp index d5ac07c09ad..9bec4789ccd 100644 --- a/stan/math/rev/functor/root_finder.hpp +++ b/stan/math/rev/functor/root_finder.hpp @@ -33,18 +33,20 @@ namespace math { * this is updated to the actual number of iterations performed * @param args Parameter pack of arguments to pass the the functors in `f_tuple` */ -template * = nullptr, require_all_stan_scalar_t* = nullptr> -auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, const GuessScalar guess, - const MinScalar min, const MaxScalar max, const int digits, +auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, + const GuessScalar guess, const MinScalar min, + const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); check_positive("root_finder", "digits", digits); check_positive("root_finder", "max_iter", max_iter); - auto arena_args_tuple = make_chainable_ptr(std::make_tuple(eval(std::forward(args))...)); + auto arena_args_tuple + = make_chainable_ptr(std::make_tuple(eval(std::forward(args))...)); auto args_vals_tuple = apply( [&](const auto&... args) { return std::make_tuple(to_ref(value_of(args))...); @@ -54,8 +56,8 @@ auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, const GuessScalar g double theta_dbl = apply( [&f_solver, &f_tuple, &max_iter, digits, guess_val = value_of(guess), min_val = value_of(min), max_val = value_of(max)](auto&&... vals) { - return root_finder_tol(f_solver, f_tuple, guess_val, min_val, max_val, digits, - max_iter, vals...); + return root_finder_tol(f_solver, f_tuple, guess_val, min_val, max_val, + digits, max_iter, vals...); }, args_vals_tuple); double Jf_x; diff --git a/test/unit/math/mix/functor/root_finder_test.cpp b/test/unit/math/mix/functor/root_finder_test.cpp index f37e1e5a722..6fbd3732d58 100644 --- a/test/unit/math/mix/functor/root_finder_test.cpp +++ b/test/unit/math/mix/functor/root_finder_test.cpp @@ -83,7 +83,7 @@ TEST(MixFun, root_finder_beta_log) { auto&& p) { std::uintmax_t max_its = 1000; return stan::math::root_finder_hailey(std::make_tuple(f, f_g, f_gg), guess, - min, max, alpha, beta, p); + min, max, alpha, beta, p); }; double p = 0.45; double alpha = .5; @@ -120,7 +120,8 @@ TEST(MixFun, root_finder_beta) { auto&& alpha, auto&& beta, auto&& p) { std::uintmax_t max_its = 1000; return stan::math::root_finder_hailey( - std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, beta, p); + std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, + beta, p); }; double p = 0.8; double alpha = .4; diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index 45af7104c4e..bcb1b2df935 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -84,8 +84,8 @@ TEST(RevFunctor, root_finder_beta_cdf) { auto&& alpha, auto&& beta, auto&& p) { std::uintmax_t max_its = 1000; return stan::math::root_finder_hailey( - std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, - alpha, beta, p); + std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, + beta, p); }; auto func2 = [&full_f](auto&& vals) { auto p = vals(0); @@ -100,7 +100,7 @@ TEST(RevFunctor, root_finder_beta_cdf) { grad_fx = var_vec.adj(); fx = fxvar.val(); std::cout << "fxvar adj:" << fxvar.adj() << "\n"; - //stan::math::gradient(func2, vals, fx, grad_fx); + // stan::math::gradient(func2, vals, fx, grad_fx); std::cout << "--- Auto Diff----\n"; std::cout << "fx: " << fx; std::cout << "\ngrads: \n" @@ -125,7 +125,7 @@ TEST(RevFunctor, root_finder_beta_cdf) { double known_alpha_grad = 0.948540678312004; double known_beta_grad = -0.7464041618483364; std::cout << "--- Auto Diff----\n"; - std::cout << "p: " << grad_fx(0) - known_p_grad << "\n" << - "alpha: " << grad_fx(1) - known_alpha_grad << "\n" << - "beta: " << grad_fx(2) - known_beta_grad << "\n"; + std::cout << "p: " << grad_fx(0) - known_p_grad << "\n" + << "alpha: " << grad_fx(1) - known_alpha_grad << "\n" + << "beta: " << grad_fx(2) - known_beta_grad << "\n"; } From 5bff38b57f9f28fd8b874daf8b85ab62c13249bc Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Wed, 11 May 2022 13:00:53 -0400 Subject: [PATCH 12/20] update test --- test/unit/math/rev/functor/root_finder_test.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index 45af7104c4e..87a1f984bca 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -44,7 +44,7 @@ TEST(RevFunctor, root_finder_beta_cdf) { vals << 0.4, .5, .5; double fx = 0; Eigen::VectorXd finit_grad_fx(3); - stan::math::finite_diff_gradient(func, vals, fx, finit_grad_fx, 1e-14); + stan::math::finite_diff_gradient(func, vals, fx, finit_grad_fx, 1e-3); std::cout << "--- Finit Diff----\n"; std::cout << "fx: " << fx; std::cout << "\ngrads: \n" @@ -121,9 +121,9 @@ TEST(RevFunctor, root_finder_beta_cdf) { "beta: " << diff_grad_fx(2) << "\n"; - double known_p_grad = 1.493914820513846; + double known_p_grad = 1.493916082370778; double known_alpha_grad = 0.948540678312004; - double known_beta_grad = -0.7464041618483364; + double known_beta_grad = -0.7464041618483361; std::cout << "--- Auto Diff----\n"; std::cout << "p: " << grad_fx(0) - known_p_grad << "\n" << "alpha: " << grad_fx(1) - known_alpha_grad << "\n" << From 26ca3fc08623cd35be7e133f79d40d420350badd Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Wed, 11 May 2022 21:00:25 -0400 Subject: [PATCH 13/20] add more rev test with explicit derivative --- .../math/rev/functor/root_finder_test.cpp | 191 +++++++++++++++++- 1 file changed, 186 insertions(+), 5 deletions(-) diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index 5e30bd06f68..c5f244d36fb 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -1,6 +1,10 @@ #include #include #include +#include +#include +#include +#include /* TEST(RevFunctor, root_finder) { using stan::math::root_finder; @@ -39,9 +43,12 @@ TEST(RevFunctor, root_finder_beta_cdf) { boost::math::beta_distribution my_beta(alpha, beta); return boost::math::quantile(my_beta, p); }; + double p = 0.4; + double a = 0.5; + double b = 0.5; Eigen::VectorXd vals(3); // p, alpha, beta - vals << 0.4, .5, .5; + vals << p, a, b; double fx = 0; Eigen::VectorXd finit_grad_fx(3); stan::math::finite_diff_gradient(func, vals, fx, finit_grad_fx, 1e-3); @@ -120,12 +127,186 @@ TEST(RevFunctor, root_finder_beta_cdf) { << "\n" "beta: " << diff_grad_fx(2) << "\n"; + auto deriv_p = [](auto& p, auto& a, auto& b) { + using std::pow; + using boost::math::ibeta_inv; + using boost::math::beta; + return beta(a, b) * pow(1 - ibeta_inv(a, b, p), (1 - b)) + * pow(ibeta_inv(a, b, p), (1 - a)); + }; + auto deriv_a = [](auto& p, auto& a, auto& b) { + using std::pow; + using boost::math::ibeta_inv; + using boost::math::ibeta; + using boost::math::beta; + using boost::math::tgamma; + using boost::math::hypergeometric_pFq; + using boost::math::beta; + using boost::math::polygamma; + using std::log; + double w = ibeta_inv(a, b, p); + return pow(1 - w, (1 - b)) * pow(w, (1 - a)) * (pow(w, a) * + pow(tgamma(a), 2) * + (hypergeometric_pFq({a, a, 1 - b}, {1 + a, 1 + a}, w) / (tgamma(1 + a) * tgamma(1 + a))) - beta(a, b) * ibeta(a, b, w) * (log(w) - polygamma(0, a) + polygamma(0, a + b))); + }; - double known_p_grad = 1.493916082370778; - double known_alpha_grad = 0.948540678312004; - double known_beta_grad = -0.7464041618483361; - std::cout << "--- Auto Diff----\n"; + auto deriv_b = [](auto& p, auto& a, auto& b) { + using std::pow; + using boost::math::ibeta_inv; + using boost::math::ibeta; + using boost::math::beta; + using boost::math::tgamma; + using boost::math::hypergeometric_pFq; + using boost::math::beta; + using boost::math::polygamma; + using std::log; + return pow(1 - ibeta_inv(a, b, p), -b) * (-1 + ibeta_inv(a, b, p)) + * pow(ibeta_inv(a, b, p), (1 - a)) + * (pow(tgamma(b), 2) + * (hypergeometric_pFq({b, b, 1 - a}, {1 + b, 1 + b}, + 1 - ibeta_inv(a, b, p))/ (tgamma(1 + b) * tgamma(1 + b))) + * pow(1 - ibeta_inv(a, b, p), b) + - beta(b, a, 1 - ibeta_inv(a, b, p)) + * (log(1 - ibeta_inv(a, b, p)) - polygamma(0, b) + + polygamma(0, a + b))); + }; + double known_p_grad = deriv_p(p, a, b); + double known_alpha_grad = deriv_a(p, a, b); + double known_beta_grad = deriv_b(p, a, b); + std::cout << "--- Mathematica Calculate Grad----\n"; + std::cout << "p: " << known_p_grad << "\n" + << "alpha: " << known_alpha_grad << "\n" + << "beta: " << known_beta_grad << "\n"; + std::cout << "--- Mathematica Calculate Grad Diff----\n"; std::cout << "p: " << grad_fx(0) - known_p_grad << "\n" << "alpha: " << grad_fx(1) - known_alpha_grad << "\n" << "beta: " << grad_fx(2) - known_beta_grad << "\n"; } + +template +inline constexpr auto make_index_tuple(const std::index_sequence&) { + return std::make_tuple(I...); +} +template +void check_vs_known_grads(FTuple&& grad_tuple, FGrad&& f, double tolerance, Args&&... args) { + try{ + const stan::math::nested_rev_autodiff nested; + auto var_tuple = std::make_tuple(stan::math::var(args)...); + auto arg_num_tuple = stan::math::apply([&args...](auto&&... nums) { + return std::make_tuple(std::make_tuple(args, nums)...); + }, make_index_tuple(std::make_index_sequence())); + auto ret = stan::math::apply([&f](auto&&... var_args) { + return f(var_args...); + }, var_tuple); + ret.grad(); + auto grad_val_tuple = stan::math::apply([&args...](auto&&... grad_funcs) { + return std::make_tuple(grad_funcs(args...)...); + }, grad_tuple); + auto adj_tuple = stan::math::apply([](auto&&... var_arg) { + return std::make_tuple(var_arg.adj()...); + }, var_tuple); + stan::math::for_each([tolerance](auto&& grad, auto&& adj, auto&& num_helper) { + EXPECT_NEAR(grad, adj, tolerance) << + "Diff: (" << grad - adj << ")\nWith arg #" << std::get<1>(num_helper) << " (" << std::get<0>(num_helper) << ")"; + }, grad_val_tuple, adj_tuple, arg_num_tuple); + } catch (const std::exception& e) { + stan::math::recover_memory(); + } + stan::math::recover_memory(); +} +TEST(RevFunctor, root_finder_beta_cdf2) { + auto beta_pdf_deriv = [](auto&& x, auto&& a, auto&& b, auto&& p) { + using stan::math::lgamma; + using stan::math::pow; + double beta_ab = boost::math::beta(a, b); + auto val = ((-1 + a) * pow(1 - x, -1 + b) * pow(x, -2 + a)) / beta_ab + - ((-1 + b) * pow(1 - x, -2 + b) * pow(x, -1 + a)) / beta_ab; + return val; + }; + auto beta_pdf = [](auto&& x, auto&& a, auto&& b, auto&& p) { + using stan::math::lgamma; + using stan::math::pow; + auto val = (pow(1 - x, -1 + b) * pow(x, -1 + a)) / boost::math::beta(a, b); + return val; + }; + + auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + return stan::math::beta_cdf(x, alpha, beta) - p; + }; + constexpr double guess = .5; + constexpr double min = 0; + constexpr double max = 1; + auto f_hailey = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( + auto&& alpha, auto&& beta, auto&& p) { + return stan::math::root_finder_hailey( + std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, + beta, p); + }; + auto f_newton = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( + auto&& alpha, auto&& beta, auto&& p) { + return stan::math::root_finder_newton_raphson( + std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, + beta, p); + }; + auto f_schroder = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( + auto&& alpha, auto&& beta, auto&& p) { + return stan::math::root_finder_schroder( + std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, + beta, p); + }; + + auto deriv_a = [](auto&& a, auto&& b, auto&& p) { + using std::pow; + using boost::math::ibeta_inv; + using boost::math::ibeta; + using boost::math::beta; + using boost::math::tgamma; + using boost::math::hypergeometric_pFq; + using boost::math::beta; + using boost::math::polygamma; + using std::log; + double w = ibeta_inv(a, b, p); + return pow(1 - w, (1 - b)) * pow(w, (1 - a)) * (pow(w, a) * + pow(tgamma(a), 2) * + (hypergeometric_pFq({a, a, 1 - b}, {1 + a, 1 + a}, w) / (tgamma(1 + a) * tgamma(1 + a))) - beta(a, b) * ibeta(a, b, w) * (log(w) - polygamma(0, a) + polygamma(0, a + b))); + }; + auto deriv_b = [](auto&& a, auto&& b, auto&& p) { + using std::pow; + using boost::math::ibeta_inv; + using boost::math::ibeta; + using boost::math::beta; + using boost::math::tgamma; + using boost::math::hypergeometric_pFq; + using boost::math::beta; + using boost::math::polygamma; + using std::log; + return pow(1 - ibeta_inv(a, b, p), -b) * (-1 + ibeta_inv(a, b, p)) + * pow(ibeta_inv(a, b, p), (1 - a)) + * (pow(tgamma(b), 2) + * (hypergeometric_pFq({b, b, 1 - a}, {1 + b, 1 + b}, + 1 - ibeta_inv(a, b, p))/ (tgamma(1 + b) * tgamma(1 + b))) + * pow(1 - ibeta_inv(a, b, p), b) + - beta(b, a, 1 - ibeta_inv(a, b, p)) + * (log(1 - ibeta_inv(a, b, p)) - polygamma(0, b) + + polygamma(0, a + b))); + }; + auto deriv_p = [](auto&& a, auto&& b, auto&& p) { + using std::pow; + using boost::math::ibeta_inv; + using boost::math::beta; + return beta(a, b) * pow(1 - ibeta_inv(a, b, p), (1 - b)) + * pow(ibeta_inv(a, b, p), (1 - a)); + }; + for (double a = .1; a < .9; a += .1) { + for (double b = .1; b < .9; b += .1) { + for (double p = .1; p < .9; p += .1) { + check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), f_schroder, 5e-4, a, b, p); + if (::testing::Test::HasFailure()) { + std::cout << "--\na: " << a << "\nb: " << b << "\np: " << p << "\n"; + } + } + } + } + + +} From 91a748c4055b8dd7d7f0f8c6d324653fd0bf3f01 Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Wed, 11 May 2022 21:01:21 -0400 Subject: [PATCH 14/20] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- .../math/rev/functor/root_finder_test.cpp | 81 +++++++++++-------- 1 file changed, 48 insertions(+), 33 deletions(-) diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index c5f244d36fb..16957e2d067 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -145,9 +145,12 @@ TEST(RevFunctor, root_finder_beta_cdf) { using boost::math::polygamma; using std::log; double w = ibeta_inv(a, b, p); - return pow(1 - w, (1 - b)) * pow(w, (1 - a)) * (pow(w, a) * - pow(tgamma(a), 2) * - (hypergeometric_pFq({a, a, 1 - b}, {1 + a, 1 + a}, w) / (tgamma(1 + a) * tgamma(1 + a))) - beta(a, b) * ibeta(a, b, w) * (log(w) - polygamma(0, a) + polygamma(0, a + b))); + return pow(1 - w, (1 - b)) * pow(w, (1 - a)) + * (pow(w, a) * pow(tgamma(a), 2) + * (hypergeometric_pFq({a, a, 1 - b}, {1 + a, 1 + a}, w) + / (tgamma(1 + a) * tgamma(1 + a))) + - beta(a, b) * ibeta(a, b, w) + * (log(w) - polygamma(0, a) + polygamma(0, a + b))); }; auto deriv_b = [](auto& p, auto& a, auto& b) { @@ -164,7 +167,8 @@ TEST(RevFunctor, root_finder_beta_cdf) { * pow(ibeta_inv(a, b, p), (1 - a)) * (pow(tgamma(b), 2) * (hypergeometric_pFq({b, b, 1 - a}, {1 + b, 1 + b}, - 1 - ibeta_inv(a, b, p))/ (tgamma(1 + b) * tgamma(1 + b))) + 1 - ibeta_inv(a, b, p)) + / (tgamma(1 + b) * tgamma(1 + b))) * pow(1 - ibeta_inv(a, b, p), b) - beta(b, a, 1 - ibeta_inv(a, b, p)) * (log(1 - ibeta_inv(a, b, p)) - polygamma(0, b) @@ -188,27 +192,35 @@ inline constexpr auto make_index_tuple(const std::index_sequence&) { return std::make_tuple(I...); } template -void check_vs_known_grads(FTuple&& grad_tuple, FGrad&& f, double tolerance, Args&&... args) { - try{ +void check_vs_known_grads(FTuple&& grad_tuple, FGrad&& f, double tolerance, + Args&&... args) { + try { const stan::math::nested_rev_autodiff nested; auto var_tuple = std::make_tuple(stan::math::var(args)...); - auto arg_num_tuple = stan::math::apply([&args...](auto&&... nums) { - return std::make_tuple(std::make_tuple(args, nums)...); - }, make_index_tuple(std::make_index_sequence())); - auto ret = stan::math::apply([&f](auto&&... var_args) { - return f(var_args...); - }, var_tuple); + auto arg_num_tuple = stan::math::apply( + [&args...](auto&&... nums) { + return std::make_tuple(std::make_tuple(args, nums)...); + }, + make_index_tuple(std::make_index_sequence())); + auto ret = stan::math::apply( + [&f](auto&&... var_args) { return f(var_args...); }, var_tuple); ret.grad(); - auto grad_val_tuple = stan::math::apply([&args...](auto&&... grad_funcs) { - return std::make_tuple(grad_funcs(args...)...); - }, grad_tuple); - auto adj_tuple = stan::math::apply([](auto&&... var_arg) { - return std::make_tuple(var_arg.adj()...); - }, var_tuple); - stan::math::for_each([tolerance](auto&& grad, auto&& adj, auto&& num_helper) { - EXPECT_NEAR(grad, adj, tolerance) << - "Diff: (" << grad - adj << ")\nWith arg #" << std::get<1>(num_helper) << " (" << std::get<0>(num_helper) << ")"; - }, grad_val_tuple, adj_tuple, arg_num_tuple); + auto grad_val_tuple = stan::math::apply( + [&args...](auto&&... grad_funcs) { + return std::make_tuple(grad_funcs(args...)...); + }, + grad_tuple); + auto adj_tuple = stan::math::apply( + [](auto&&... var_arg) { return std::make_tuple(var_arg.adj()...); }, + var_tuple); + stan::math::for_each( + [tolerance](auto&& grad, auto&& adj, auto&& num_helper) { + EXPECT_NEAR(grad, adj, tolerance) + << "Diff: (" << grad - adj << ")\nWith arg #" + << std::get<1>(num_helper) << " (" << std::get<0>(num_helper) + << ")"; + }, + grad_val_tuple, adj_tuple, arg_num_tuple); } catch (const std::exception& e) { stan::math::recover_memory(); } @@ -220,7 +232,7 @@ TEST(RevFunctor, root_finder_beta_cdf2) { using stan::math::pow; double beta_ab = boost::math::beta(a, b); auto val = ((-1 + a) * pow(1 - x, -1 + b) * pow(x, -2 + a)) / beta_ab - - ((-1 + b) * pow(1 - x, -2 + b) * pow(x, -1 + a)) / beta_ab; + - ((-1 + b) * pow(1 - x, -2 + b) * pow(x, -1 + a)) / beta_ab; return val; }; auto beta_pdf = [](auto&& x, auto&& a, auto&& b, auto&& p) { @@ -237,19 +249,19 @@ TEST(RevFunctor, root_finder_beta_cdf2) { constexpr double min = 0; constexpr double max = 1; auto f_hailey = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( - auto&& alpha, auto&& beta, auto&& p) { + auto&& alpha, auto&& beta, auto&& p) { return stan::math::root_finder_hailey( std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, beta, p); }; auto f_newton = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( - auto&& alpha, auto&& beta, auto&& p) { + auto&& alpha, auto&& beta, auto&& p) { return stan::math::root_finder_newton_raphson( std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, beta, p); }; auto f_schroder = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( - auto&& alpha, auto&& beta, auto&& p) { + auto&& alpha, auto&& beta, auto&& p) { return stan::math::root_finder_schroder( std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, beta, p); @@ -266,9 +278,12 @@ TEST(RevFunctor, root_finder_beta_cdf2) { using boost::math::polygamma; using std::log; double w = ibeta_inv(a, b, p); - return pow(1 - w, (1 - b)) * pow(w, (1 - a)) * (pow(w, a) * - pow(tgamma(a), 2) * - (hypergeometric_pFq({a, a, 1 - b}, {1 + a, 1 + a}, w) / (tgamma(1 + a) * tgamma(1 + a))) - beta(a, b) * ibeta(a, b, w) * (log(w) - polygamma(0, a) + polygamma(0, a + b))); + return pow(1 - w, (1 - b)) * pow(w, (1 - a)) + * (pow(w, a) * pow(tgamma(a), 2) + * (hypergeometric_pFq({a, a, 1 - b}, {1 + a, 1 + a}, w) + / (tgamma(1 + a) * tgamma(1 + a))) + - beta(a, b) * ibeta(a, b, w) + * (log(w) - polygamma(0, a) + polygamma(0, a + b))); }; auto deriv_b = [](auto&& a, auto&& b, auto&& p) { using std::pow; @@ -284,7 +299,8 @@ TEST(RevFunctor, root_finder_beta_cdf2) { * pow(ibeta_inv(a, b, p), (1 - a)) * (pow(tgamma(b), 2) * (hypergeometric_pFq({b, b, 1 - a}, {1 + b, 1 + b}, - 1 - ibeta_inv(a, b, p))/ (tgamma(1 + b) * tgamma(1 + b))) + 1 - ibeta_inv(a, b, p)) + / (tgamma(1 + b) * tgamma(1 + b))) * pow(1 - ibeta_inv(a, b, p), b) - beta(b, a, 1 - ibeta_inv(a, b, p)) * (log(1 - ibeta_inv(a, b, p)) - polygamma(0, b) @@ -300,13 +316,12 @@ TEST(RevFunctor, root_finder_beta_cdf2) { for (double a = .1; a < .9; a += .1) { for (double b = .1; b < .9; b += .1) { for (double p = .1; p < .9; p += .1) { - check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), f_schroder, 5e-4, a, b, p); + check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), + f_schroder, 5e-4, a, b, p); if (::testing::Test::HasFailure()) { std::cout << "--\na: " << a << "\nb: " << b << "\np: " << p << "\n"; } } } } - - } From ee38928119358e92e83fb47962d3fe88ee0f19bb Mon Sep 17 00:00:00 2001 From: stevebronder Date: Mon, 23 May 2022 17:01:55 -0400 Subject: [PATCH 15/20] update root finder to pass the function F at compile time --- stan/math/prim/functor/root_finder.hpp | 176 ++++++++++-------- stan/math/rev/functor/root_finder.hpp | 65 ++++--- .../math/rev/functor/root_finder_test.cpp | 135 +++++++------- 3 files changed, 193 insertions(+), 183 deletions(-) diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index 10a95a5b96c..7be54423261 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -13,28 +13,60 @@ namespace stan { namespace math { namespace internal { -template -inline auto func_with_derivs(Tuple&& f_tuple, Args&&... args) { - return stan::math::apply( - [&args...](auto&&... funcs) { - return [&args..., &funcs...](auto&& g) { - return std::make_tuple(funcs(g, args...)...); - }; - }, - f_tuple); +template * = nullptr> +inline auto make_root_func(Args&&... args) { + return [&args...](auto&& x) { + return std::decay_t::template run(x, args...); + }; } + +template * = nullptr> +inline auto make_root_func() { + return [](auto&&... args) { + return std::decay_t::template run(args...); + }; +} + +struct NewtonRootSolver { + template + static inline auto run(Types&&... args) { + return boost::math::tools::newton_raphson_iterate( + std::forward(args)...); + } +}; + +struct HalleyRootSolver { + template + static inline auto run(Types&&... args) { + return boost::math::tools::halley_iterate(std::forward(args)...); + } +}; + +struct SchroderRootSolver { + template + static inline auto run(Types&&... args) { + return boost::math::tools::schroder_iterate(std::forward(args)...); + } +}; + } // namespace internal /** * Solve for root using Boost's Halley method - * @tparam FTuple A tuple holding functors whose signatures all match - * `(GuessScalar g, Types&&... Args)`. + * @tparam FRootFunc A struct or class with a static function called `run`. + * The structs `run` function must have a boolean template parameter that + * when `true` returns a tuple containing the function result and the + * derivatives needed for the root finder. When the boolean template parameter + * is `false` the function should return a single value containing the function + * result. + * @tparam SolverFun One of the three struct types used to call the root solver. + * (`NewtonRootSolver`, `HalleyRootSolver`, `SchroderRootSolver`). * @tparam GuessScalar Scalar type * @tparam MinScalar Scalar type * @tparam MaxScalar Scalar type * @tparam Types Arg types to pass to functors in `f_tuple` - * @param f_tuple A tuple of functors to calculate the value and any derivates - * needed. * @param guess An initial guess at the root value * @param min The minimum possible value for the result, this is used as an * initial lower bracket @@ -45,12 +77,11 @@ inline auto func_with_derivs(Tuple&& f_tuple, Args&&... args) { * this is updated to the actual number of iterations performed * @param args Parameter pack of arguments to pass the the functors in `f_tuple` */ -template * = nullptr> -auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, - const GuessScalar guess, const MinScalar min, +auto root_finder_tol(const GuessScalar guess, const MinScalar min, const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); @@ -58,10 +89,11 @@ auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, check_positive("root_finder", "max_iter", max_iter); using ret_t = return_type_t; ret_t ret = 0; - auto f_plus_div = internal::func_with_derivs(f_tuple, args...); + auto f_plus_div + = internal::make_root_func(std::forward(args)...); try { - ret = f_solver(f_plus_div, ret_t(guess), ret_t(min), ret_t(max), digits, - max_iter); + ret = std::decay_t::run(f_plus_div, ret_t(guess), ret_t(min), + ret_t(max), digits, max_iter); } catch (const std::exception& e) { std::cout << "err: \n" << e.what() << "\n"; throw e; @@ -69,59 +101,48 @@ auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, return ret; } -template -auto root_finder_halley_tol(FTuple&& f_tuple, const GuessScalar guess, - const MinScalar min, const MaxScalar max, - const int digits, std::uintmax_t& max_iter, - Types&&... args) { - return root_finder_tol( - [](auto&&... args) { - return boost::math::tools::halley_iterate(args...); - }, - std::forward(f_tuple), guess, min, max, digits, max_iter, - std::forward(args)...); +auto root_finder_halley_tol(const GuessScalar guess, const MinScalar min, + const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types&&... args) { + return root_finder_tol( + guess, min, max, digits, max_iter, std::forward(args)...); } -template -auto root_finder_newton_raphson_tol(FTuple&& f_tuple, const GuessScalar guess, +auto root_finder_newton_raphson_tol(const GuessScalar guess, const MinScalar min, const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { - return root_finder_tol( - [](auto&&... args) { - return boost::math::tools::newton_raphson_iterate(args...); - }, - std::forward(f_tuple), guess, min, max, digits, max_iter, - std::forward(args)...); + return root_finder_tol( + guess, min, max, digits, max_iter, std::forward(args)...); } -template -auto root_finder_schroder_tol(FTuple&& f_tuple, const GuessScalar guess, - const MinScalar min, const MaxScalar max, - const int digits, std::uintmax_t& max_iter, - Types&&... args) { - return root_finder_tol( - [](auto&&... args) { - return boost::math::tools::schroder_iterate(args...); - }, - std::forward(f_tuple), guess, min, max, digits, max_iter, - std::forward(args)...); +auto root_finder_schroder_tol(const GuessScalar guess, const MinScalar min, + const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types&&... args) { + return root_finder_tol( + guess, min, max, digits, max_iter, std::forward(args)...); } /** - * Solve for root using Boost Halley method with default values for the - * tolerances - * @tparam FTuple A tuple holding functors whose signatures all match - * `(GuessScalar g, Types&&... Args)`. + * Solve for root with default values for the tolerances + * @tparam FRootFunc A struct or class with a static function called `run`. + * The structs `run` function must have a boolean template parameter that + * when `true` returns a tuple containing the function result and the + * derivatives needed for the root finder. When the boolean template parameter + * is `false` the function should return a single value containing the function + * result. + * @tparam SolverFun One of the three struct types used to call the root solver. + * (`NewtonRootSolver`, `HalleyRootSolver`, `SchroderRootSolver`). * @tparam GuessScalar Scalar type * @tparam MinScalar Scalar type * @tparam MaxScalar Scalar type * @tparam Types Arg types to pass to functors in `f_tuple` - * @param f_tuple A tuple of functors to calculate the value and any derivates - * needed. * @param guess An initial guess at the root value * @param min The minimum possible value for the result, this is used as an * initial lower bracket @@ -129,51 +150,44 @@ auto root_finder_schroder_tol(FTuple&& f_tuple, const GuessScalar guess, * initial upper bracket * @param args Parameter pack of arguments to pass the the functors in `f_tuple` */ -template -auto root_finder(SolverFun&& f_solver, FTuple&& f_tuple, - const GuessScalar guess, const MinScalar min, +auto root_finder(const GuessScalar guess, const MinScalar min, const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); - return root_finder_tol(std::forward(f_solver), - std::forward(f_tuple), guess, min, max, digits, - max_iter, std::forward(args)...); + return root_finder_tol( + guess, min, max, digits, max_iter, std::forward(args)...); } -template -auto root_finder_hailey(FTuple&& f_tuple, const GuessScalar guess, - const MinScalar min, const MaxScalar max, - Types&&... args) { +auto root_finder_hailey(const GuessScalar guess, const MinScalar min, + const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); - return root_finder_halley_tol(std::forward(f_tuple), guess, min, max, - digits, max_iter, std::forward(args)...); + return root_finder_halley_tol(guess, min, max, digits, max_iter, + std::forward(args)...); } -template -auto root_finder_newton_raphson(FTuple&& f_tuple, const GuessScalar guess, - const MinScalar min, const MaxScalar max, - Types&&... args) { +auto root_finder_newton_raphson(const GuessScalar guess, const MinScalar min, + const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); - return root_finder_newton_raphson_tol(std::forward(f_tuple), guess, - min, max, digits, max_iter, - std::forward(args)...); + return root_finder_newton_raphson_tol( + guess, min, max, digits, max_iter, std::forward(args)...); } -template -auto root_finder_schroder(FTuple&& f_tuple, const GuessScalar guess, - const MinScalar min, const MaxScalar max, - Types&&... args) { +auto root_finder_schroder(const GuessScalar guess, const MinScalar min, + const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); - return root_finder_schroder_tol(std::forward(f_tuple), guess, min, - max, digits, max_iter, - std::forward(args)...); + return root_finder_schroder_tol(guess, min, max, digits, max_iter, + std::forward(args)...); } } // namespace math diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp index 9bec4789ccd..313928f2f6e 100644 --- a/stan/math/rev/functor/root_finder.hpp +++ b/stan/math/rev/functor/root_finder.hpp @@ -15,14 +15,18 @@ namespace math { /** * var specialization for root solving using Boost's Halley method - * @tparam FTuple A tuple holding functors whose signatures all match - * `(GuessScalar g, Types&&... Args)`. + * @tparam FRootFunc A struct or class with a static function called `run`. + * The structs `run` function must have a boolean template parameter that + * when `true` returns a tuple containing the function result and the + * derivatives needed for the root finder. When the boolean template parameter + * is `false` the function should return a single value containing the function + * result. + * @tparam SolverFun One of the three struct types used to call the root solver. + * (`NewtonRootSolver`, `HalleyRootSolver`, `SchroderRootSolver`). * @tparam GuessScalar Scalar type * @tparam MinScalar Scalar type * @tparam MaxScalar Scalar type * @tparam Types Arg types to pass to functors in `f_tuple` - * @param f_tuple A tuple of functors to calculate the value and any derivates - * needed. * @param guess An initial guess at the root value * @param min The minimum possible value for the result, this is used as an * initial lower bracket @@ -34,12 +38,11 @@ namespace math { * @param args Parameter pack of arguments to pass the the functors in `f_tuple` */ template < - typename SolverFun, typename FTuple, typename GuessScalar, + typename FRootFunc, typename SolverFun, typename GuessScalar, typename MinScalar, typename MaxScalar, typename... Types, require_any_st_var* = nullptr, require_all_stan_scalar_t* = nullptr> -auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, - const GuessScalar guess, const MinScalar min, +auto root_finder_tol(const GuessScalar guess, const MinScalar min, const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); @@ -54,19 +57,21 @@ auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, *arena_args_tuple); // Solve the system double theta_dbl = apply( - [&f_solver, &f_tuple, &max_iter, digits, guess_val = value_of(guess), - min_val = value_of(min), max_val = value_of(max)](auto&&... vals) { - return root_finder_tol(f_solver, f_tuple, guess_val, min_val, max_val, - digits, max_iter, vals...); + [&max_iter, digits, guess_val = value_of(guess), min_val = value_of(min), + max_val = value_of(max)](auto&&... vals) { + return root_finder_tol( + guess_val, min_val, max_val, digits, max_iter, vals...); }, args_vals_tuple); double Jf_x; - auto&& f = std::get<0>(f_tuple); { nested_rev_autodiff nested; stan::math::var x_var(theta_dbl); stan::math::var fx_var = apply( - [&x_var, &f](auto&&... args) { return f(x_var, std::move(args)...); }, + [&x_var](auto&&... args) { + return std::decay_t::template run( + x_var, std::move(args)...); + }, std::move(args_vals_tuple)); fx_var.grad(); Jf_x = x_var.adj(); @@ -76,24 +81,22 @@ auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, * Note: Because we put this on the callback stack, if `f` is a lambda * its captures must be in Stan's arena memory or trivially destructable. */ - return make_callback_var(theta_dbl, - [f, arena_args_tuple, Jf_x](auto& ret) mutable { - // Eigen::VectorXd eta = - // -Jf_x_T_lu_ptr->solve(ret.adj().eval()); - // Contract with Jacobian of f with respect to y - // using a nested reverse autodiff pass. - { - nested_rev_autodiff rev; - double ret_val = ret.val(); - auto x_nrad_ = apply( - [&ret_val, &f](const auto&... args) { - return eval(f(ret_val, args...)); - }, - *arena_args_tuple); - x_nrad_.adj() = -(ret.adj() / Jf_x); - grad(); - } - }); + return make_callback_var( + theta_dbl, [arena_args_tuple, Jf_x](auto& ret) mutable { + { + nested_rev_autodiff rev; + double ret_val = ret.val(); + auto x_nrad_ = apply( + [&ret_val](const auto&... args) { + auto f = internal::make_root_func(); + return eval(std::decay_t::template run( + ret_val, args...)); + }, + *arena_args_tuple); + x_nrad_.adj() = -(ret.adj() / Jf_x); + grad(); + } + }); } } // namespace math diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index 16957e2d067..7aafe239cbd 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -34,6 +34,39 @@ TEST(RevFunctor, root_finder) { EXPECT_NEAR(0.037037, x.adj(), 1e-6); } */ + +struct BetaCdfRoot { + template * = nullptr> + static auto run(T1&& x, T2&& alpha, T3&& beta, T4&& p) { + auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + return stan::math::beta_cdf(x, alpha, beta) - p; + }(x, alpha, beta, p); + auto beta_pdf_deriv = [](auto&& x, auto&& alpha, auto&& beta) { + using stan::math::lgamma; + using stan::math::pow; + auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) + * (1. - 2. * x - alpha + x * alpha + x * beta) + * lgamma(alpha + beta)) + / (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); + return val; + }(x, alpha, beta); + auto beta_pdf = [](auto&& x, auto&& alpha, auto&& beta) { + using stan::math::lgamma; + using stan::math::pow; + auto val + = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta)) + / ((lgamma(alpha) * lgamma(beta))); + return val; + }(x, alpha, beta); + return std::make_tuple(f, beta_pdf, beta_pdf_deriv); + } + template * = nullptr> + static auto run(T1&& x, T2&& alpha, T3&& beta, T4&& p) { + return stan::math::beta_cdf(x, alpha, beta) - p; + } +}; TEST(RevFunctor, root_finder_beta_cdf) { using stan::math::var; auto func = [](auto&& vals) { @@ -62,37 +95,13 @@ TEST(RevFunctor, root_finder_beta_cdf) { << "\n" "beta: " << finit_grad_fx(2) << "\n"; - - auto beta_pdf_deriv = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - using stan::math::lgamma; - using stan::math::pow; - auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) - * (1. - 2. * x - alpha + x * alpha + x * beta) - * lgamma(alpha + beta)) - / (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); - return val; - }; - auto beta_pdf = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - using stan::math::lgamma; - using stan::math::pow; - auto val - = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta)) - / ((lgamma(alpha) * lgamma(beta))); - return val; - }; - - auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - return stan::math::beta_cdf(x, alpha, beta) - p; - }; double guess = .3; double min = 0; double max = 1; - auto full_f = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( - auto&& alpha, auto&& beta, auto&& p) { + auto full_f = [guess, min, max](auto&& alpha, auto&& beta, auto&& p) { std::uintmax_t max_its = 1000; - return stan::math::root_finder_hailey( - std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, - beta, p); + return stan::math::root_finder_hailey(guess, min, max, alpha, + beta, p); }; auto func2 = [&full_f](auto&& vals) { auto p = vals(0); @@ -197,11 +206,13 @@ void check_vs_known_grads(FTuple&& grad_tuple, FGrad&& f, double tolerance, try { const stan::math::nested_rev_autodiff nested; auto var_tuple = std::make_tuple(stan::math::var(args)...); + // For pretty printer index of incorrect values auto arg_num_tuple = stan::math::apply( [&args...](auto&&... nums) { return std::make_tuple(std::make_tuple(args, nums)...); }, make_index_tuple(std::make_index_sequence())); + auto ret = stan::math::apply( [&f](auto&&... var_args) { return f(var_args...); }, var_tuple); ret.grad(); @@ -227,45 +238,9 @@ void check_vs_known_grads(FTuple&& grad_tuple, FGrad&& f, double tolerance, stan::math::recover_memory(); } TEST(RevFunctor, root_finder_beta_cdf2) { - auto beta_pdf_deriv = [](auto&& x, auto&& a, auto&& b, auto&& p) { - using stan::math::lgamma; - using stan::math::pow; - double beta_ab = boost::math::beta(a, b); - auto val = ((-1 + a) * pow(1 - x, -1 + b) * pow(x, -2 + a)) / beta_ab - - ((-1 + b) * pow(1 - x, -2 + b) * pow(x, -1 + a)) / beta_ab; - return val; - }; - auto beta_pdf = [](auto&& x, auto&& a, auto&& b, auto&& p) { - using stan::math::lgamma; - using stan::math::pow; - auto val = (pow(1 - x, -1 + b) * pow(x, -1 + a)) / boost::math::beta(a, b); - return val; - }; - - auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - return stan::math::beta_cdf(x, alpha, beta) - p; - }; constexpr double guess = .5; constexpr double min = 0; constexpr double max = 1; - auto f_hailey = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( - auto&& alpha, auto&& beta, auto&& p) { - return stan::math::root_finder_hailey( - std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, - beta, p); - }; - auto f_newton = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( - auto&& alpha, auto&& beta, auto&& p) { - return stan::math::root_finder_newton_raphson( - std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, - beta, p); - }; - auto f_schroder = [&f, &beta_pdf_deriv, &beta_pdf, guess, min, max]( - auto&& alpha, auto&& beta, auto&& p) { - return stan::math::root_finder_schroder( - std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, - beta, p); - }; auto deriv_a = [](auto&& a, auto&& b, auto&& p) { using std::pow; @@ -313,15 +288,33 @@ TEST(RevFunctor, root_finder_beta_cdf2) { return beta(a, b) * pow(1 - ibeta_inv(a, b, p), (1 - b)) * pow(ibeta_inv(a, b, p), (1 - a)); }; - for (double a = .1; a < .9; a += .1) { - for (double b = .1; b < .9; b += .1) { - for (double p = .1; p < .9; p += .1) { - check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), - f_schroder, 5e-4, a, b, p); - if (::testing::Test::HasFailure()) { - std::cout << "--\na: " << a << "\nb: " << b << "\np: " << p << "\n"; - } + auto f_newton = [guess, min, max](auto&& alpha, auto&& beta, auto&& p) { + return stan::math::root_finder_newton_raphson(guess, min, max, + alpha, beta, p); + }; + auto f_schroder = [guess, min, max](auto&& alpha, auto&& beta, auto&& p) { + return stan::math::root_finder_schroder(guess, min, max, alpha, + beta, p); + }; + auto f_hailey = [](auto&& alpha, auto&& beta, auto&& p) { + return stan::math::root_finder_hailey(p, min, max, alpha, beta, + p); + }; + check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), f_schroder, + 5e-2, .5, .5, .3); + + /* For some reason this fails after a while?? + for (double p = .3; p < .7; p += .1) { + for (double a = .3; a < .7; a += .1) { + for (double b = .3; b < .7; b += .1) { + + check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), + f_schroder, 5e-2, a, b, p); + if (::testing::Test::HasFailure()) { + std::cout << "--\na: " << a << "\nb: " << b << "\np: " << p << "\n"; + } } } } + */ } From 344a450ca62d537b29833c668d946ef572344f9e Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Tue, 24 May 2022 13:10:09 -0400 Subject: [PATCH 16/20] fix sign --- stan/math/fwd/fun/sign.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/fwd/fun/sign.hpp b/stan/math/fwd/fun/sign.hpp index a46ae60ba3c..fcdfb1805a8 100644 --- a/stan/math/fwd/fun/sign.hpp +++ b/stan/math/fwd/fun/sign.hpp @@ -10,8 +10,8 @@ namespace math { template inline auto sign(const fvar& x) { - double x_val = value_of_rec(x); - return (0. < x_val) - (x_val < 0.); + double z = value_of_rec(x); + return (z == 0) ? 0 : z < 0 ? -1 : 1; } } // namespace math From ef8d4e0e55dab0e17c642d81bd83cd9407def714 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Tue, 24 May 2022 13:10:59 -0400 Subject: [PATCH 17/20] fix sign --- stan/math/prim/functor/root_finder.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index 10a95a5b96c..eceb63c83b3 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -63,7 +63,6 @@ auto root_finder_tol(SolverFun&& f_solver, FTuple&& f_tuple, ret = f_solver(f_plus_div, ret_t(guess), ret_t(min), ret_t(max), digits, max_iter); } catch (const std::exception& e) { - std::cout << "err: \n" << e.what() << "\n"; throw e; } return ret; From 19588e44c5a0004be225d8f5785f05f3489746f0 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Wed, 25 May 2022 16:22:10 -0400 Subject: [PATCH 18/20] turn on fifth root mix check --- .../math/mix/functor/root_finder_test.cpp | 144 +++++++++--------- 1 file changed, 68 insertions(+), 76 deletions(-) diff --git a/test/unit/math/mix/functor/root_finder_test.cpp b/test/unit/math/mix/functor/root_finder_test.cpp index 6fbd3732d58..8e51d316747 100644 --- a/test/unit/math/mix/functor/root_finder_test.cpp +++ b/test/unit/math/mix/functor/root_finder_test.cpp @@ -1,5 +1,4 @@ #include - /* TEST(MixFun, root_finder_cubed) { using stan::math::root_finder; @@ -27,6 +26,27 @@ TEST(MixFun, root_finder_cubed) { }; stan::test::expect_ad(full_f, x); } +*/ + + +struct FifthRootFinder { + template * = nullptr> + static auto run(T1&& g, T2&& x) { + auto g_pow = g * g * g; + auto second_deriv = 20 * g_pow; + g_pow *= g; + auto first_deriv = 5 * g_pow; + g_pow *= g; + auto func_val = g_pow - x; + return std::make_tuple(func_val, first_deriv, second_deriv); + } + template * = nullptr> + static auto run(T1&& g, T2&& x) { + return g * g * g * g * g - x; + } +}; TEST(MixFun, root_finder_fifth) { using stan::math::root_finder; @@ -46,85 +66,57 @@ TEST(MixFun, root_finder_fifth) { const int digits = std::numeric_limits::digits; int get_digits = static_cast(digits * 0.4); std::uintmax_t maxit = 20; - auto f = [](const auto& g, const auto& x) { return g * g * g * g * g - x; }; - auto f_g = [](const auto& g, const auto& x) { return 5 * g * g * g * g; }; - auto f_gg = [](const auto& g, const auto& x) { return 20 * g * g * g; }; - auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& xx) { - return root_finder(std::make_tuple(f, f_g, f_gg), guess, min, max, xx); - }; - stan::test::expect_ad(full_f, x); -} -*/ -TEST(MixFun, root_finder_beta_log) { - auto beta_lpdf_deriv = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - using stan::math::log; - using stan::math::lgamma; - using stan::math::log1m; - auto val = -((beta * log1m(x) + (-2 + alpha) * log(x) - + log(1 - 2 * x - alpha + x * alpha + x * beta) - + lgamma(alpha + beta)) - - (2 * log1m(x) + lgamma(alpha) + lgamma(beta))); - return val; - }; - auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - return stan::math::beta_lcdf(x, alpha, beta) - p; - }; - auto f_g = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - return stan::math::beta_lpdf(x, alpha, beta); + auto f_hailey = [guess, min, max](auto&& x) { + return stan::math::root_finder_hailey(guess, min, max, x); }; - auto f_gg - = [&beta_lpdf_deriv](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - return beta_lpdf_deriv(x, alpha, beta, p); - }; - double guess = .3; - double min = 0.0000000001; - double max = 1; - auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& alpha, auto&& beta, - auto&& p) { - std::uintmax_t max_its = 1000; - return stan::math::root_finder_hailey(std::make_tuple(f, f_g, f_gg), guess, - min, max, alpha, beta, p); - }; - double p = 0.45; - double alpha = .5; - double beta = .4; - stan::test::expect_ad(full_f, alpha, beta, p); + stan::test::expect_ad(f_hailey, x); } -TEST(MixFun, root_finder_beta) { - auto beta_pdf_deriv = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - using stan::math::lgamma; - using stan::math::pow; - auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) - * (1. - 2. * x - alpha + x * alpha + x * beta) - * lgamma(alpha + beta)) - / (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); - return val; - }; - auto beta_pdf = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - using stan::math::lgamma; - using stan::math::pow; - auto val - = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta)) - / ((lgamma(alpha) * lgamma(beta))); - return val; - }; - - auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { +/* +struct BetaCdfRoot { + template * = nullptr> + static auto run(T1&& x, T2&& alpha, T3&& beta, T4&& p) { + auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { + return stan::math::beta_cdf(x, alpha, beta) - p; + }(x, alpha, beta, p); + auto beta_pdf_deriv = [](auto&& x, auto&& alpha, auto&& beta) { + using stan::math::lgamma; + using stan::math::pow; + auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) + * (1. - 2. * x - alpha + x * alpha + x * beta) + * lgamma(alpha + beta)) + / (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); + return val; + }(x, alpha, beta); + auto beta_pdf = [](auto&& x, auto&& alpha, auto&& beta) { + using stan::math::lgamma; + using stan::math::pow; + auto val + = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta)) + / ((lgamma(alpha) * lgamma(beta))); + return val; + }(x, alpha, beta); + return std::make_tuple(f, beta_pdf, beta_pdf_deriv); + } + template * = nullptr> + static auto run(T1&& x, T2&& alpha, T3&& beta, T4&& p) { return stan::math::beta_cdf(x, alpha, beta) - p; + } +}; + +TEST(MixFun, root_finder_beta) { + constexpr double guess = .5; + constexpr double min = 0; + constexpr double max = 1; + auto f_hailey = [](auto&& alpha, auto&& beta, auto&& p) { + return stan::math::root_finder_hailey(p, min, max, alpha, beta, + p); }; - double guess = .7; - double min = 0.0000000001; - double max = 1; - auto full_f = [&f, &beta_pdf, &beta_pdf_deriv, guess, min, max]( - auto&& alpha, auto&& beta, auto&& p) { - std::uintmax_t max_its = 1000; - return stan::math::root_finder_hailey( - std::make_tuple(f, beta_pdf, beta_pdf_deriv), guess, min, max, alpha, - beta, p); - }; - double p = 0.8; - double alpha = .4; + double alpha = .5; double beta = .5; - stan::test::expect_ad(full_f, alpha, beta, p); + double p = .3; + stan::test::expect_ad(f_hailey, alpha, beta, p); } +*/ From 13a8886f0ca164cd19d861ab5a64c155c4224055 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Fri, 3 Jun 2022 17:46:18 -0400 Subject: [PATCH 19/20] fixup tests to pass --- stan/math/prim/functor/root_finder.hpp | 16 ++-- stan/math/rev/functor/root_finder.hpp | 78 ++++++++++++++++++- .../math/mix/functor/root_finder_test.cpp | 62 +++++++-------- .../math/rev/functor/root_finder_test.cpp | 76 +++++------------- 4 files changed, 132 insertions(+), 100 deletions(-) diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index 2ec801d9291..1102f1f9bb4 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -81,7 +81,7 @@ template * = nullptr> -auto root_finder_tol(const GuessScalar guess, const MinScalar min, +inline auto root_finder_tol(const GuessScalar guess, const MinScalar min, const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); @@ -102,7 +102,7 @@ auto root_finder_tol(const GuessScalar guess, const MinScalar min, template -auto root_finder_halley_tol(const GuessScalar guess, const MinScalar min, +inline auto root_finder_halley_tol(const GuessScalar guess, const MinScalar min, const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { return root_finder_tol( @@ -111,7 +111,7 @@ auto root_finder_halley_tol(const GuessScalar guess, const MinScalar min, template -auto root_finder_newton_raphson_tol(const GuessScalar guess, +inline auto root_finder_newton_raphson_tol(const GuessScalar guess, const MinScalar min, const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { @@ -121,7 +121,7 @@ auto root_finder_newton_raphson_tol(const GuessScalar guess, template -auto root_finder_schroder_tol(const GuessScalar guess, const MinScalar min, +inline auto root_finder_schroder_tol(const GuessScalar guess, const MinScalar min, const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { return root_finder_tol( @@ -151,7 +151,7 @@ auto root_finder_schroder_tol(const GuessScalar guess, const MinScalar min, */ template -auto root_finder(const GuessScalar guess, const MinScalar min, +inline auto root_finder(const GuessScalar guess, const MinScalar min, const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); @@ -161,7 +161,7 @@ auto root_finder(const GuessScalar guess, const MinScalar min, template -auto root_finder_hailey(const GuessScalar guess, const MinScalar min, +inline auto root_finder_hailey(const GuessScalar guess, const MinScalar min, const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); @@ -171,7 +171,7 @@ auto root_finder_hailey(const GuessScalar guess, const MinScalar min, template -auto root_finder_newton_raphson(const GuessScalar guess, const MinScalar min, +inline auto root_finder_newton_raphson(const GuessScalar guess, const MinScalar min, const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); @@ -181,7 +181,7 @@ auto root_finder_newton_raphson(const GuessScalar guess, const MinScalar min, template -auto root_finder_schroder(const GuessScalar guess, const MinScalar min, +inline auto root_finder_schroder(const GuessScalar guess, const MinScalar min, const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp index 313928f2f6e..3966bb0ecb3 100644 --- a/stan/math/rev/functor/root_finder.hpp +++ b/stan/math/rev/functor/root_finder.hpp @@ -41,8 +41,9 @@ template < typename FRootFunc, typename SolverFun, typename GuessScalar, typename MinScalar, typename MaxScalar, typename... Types, require_any_st_var* = nullptr, - require_all_stan_scalar_t* = nullptr> -auto root_finder_tol(const GuessScalar guess, const MinScalar min, + require_all_stan_scalar_t* = nullptr, + require_any_not_stan_scalar_t* = nullptr> +inline auto root_finder_tol(const GuessScalar guess, const MinScalar min, const MaxScalar max, const int digits, std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); @@ -60,7 +61,7 @@ auto root_finder_tol(const GuessScalar guess, const MinScalar min, [&max_iter, digits, guess_val = value_of(guess), min_val = value_of(min), max_val = value_of(max)](auto&&... vals) { return root_finder_tol( - guess_val, min_val, max_val, digits, max_iter, vals...); + guess_val, min_val, max_val, digits > 20 ? digits : 21, max_iter, vals...); }, args_vals_tuple); double Jf_x; @@ -85,6 +86,7 @@ auto root_finder_tol(const GuessScalar guess, const MinScalar min, theta_dbl, [arena_args_tuple, Jf_x](auto& ret) mutable { { nested_rev_autodiff rev; + double eta = -(ret.adj() / Jf_x); double ret_val = ret.val(); auto x_nrad_ = apply( [&ret_val](const auto&... args) { @@ -93,7 +95,75 @@ auto root_finder_tol(const GuessScalar guess, const MinScalar min, ret_val, args...)); }, *arena_args_tuple); - x_nrad_.adj() = -(ret.adj() / Jf_x); + x_nrad_.adj() = eta; + grad(); + } + }); +} + +/** + * var specialization for root solving using Boost's Halley method + * @tparam FRootFunc A struct or class with a static function called `run`. + * The structs `run` function must have a boolean template parameter that + * when `true` returns a tuple containing the function result and the + * derivatives needed for the root finder. When the boolean template parameter + * is `false` the function should return a single value containing the function + * result. + * @tparam SolverFun One of the three struct types used to call the root solver. + * (`NewtonRootSolver`, `HalleyRootSolver`, `SchroderRootSolver`). + * @tparam GuessScalar Scalar type + * @tparam MinScalar Scalar type + * @tparam MaxScalar Scalar type + * @tparam Types Arg types to pass to functors in `f_tuple` + * @param guess An initial guess at the root value + * @param min The minimum possible value for the result, this is used as an + * initial lower bracket + * @param max The maximum possible value for the result, this is used as an + * initial upper bracket + * @param digits The desired number of binary digits precision + * @param max_iter An optional maximum number of iterations to perform. On exit, + * this is updated to the actual number of iterations performed + * @param args Parameter pack of arguments to pass the the functors in `f_tuple` + */ +template < + typename FRootFunc, typename SolverFun, typename GuessScalar, + typename MinScalar, typename MaxScalar, typename... Types, + require_any_st_var* = nullptr, + require_all_stan_scalar_t* = nullptr, + require_all_stan_scalar_t* = nullptr> +inline auto root_finder_tol(const GuessScalar guess, const MinScalar min, + const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types... args) { + check_bounded("root_finder", "initial guess", guess, min, max); + check_positive("root_finder", "digits", digits); + check_positive("root_finder", "max_iter", max_iter); + // Solve the system + double theta_dbl = root_finder_tol( + value_of(guess), value_of(min), value_of(max), digits > 20 ? digits : 21, max_iter, value_of(args)...); + double Jf_x; + { + nested_rev_autodiff nested; + stan::math::var x_var(theta_dbl); + stan::math::var fx_var = std::decay_t::template run( + x_var, value_of(args)...); + fx_var.grad(); + Jf_x = x_var.adj(); + } + + /* + * Note: Because we put this on the callback stack, if `f` is a lambda + * its captures must be in Stan's arena memory or trivially destructable. + */ + return make_callback_var( + theta_dbl, [Jf_x, args...](auto& ret) mutable { + { + nested_rev_autodiff rev; + double eta = -(ret.adj() / Jf_x); + double ret_val = ret.val(); + auto f = internal::make_root_func(); + auto x_nrad = std::decay_t::template run( + ret_val, args...); + x_nrad.adj() = eta; grad(); } }); diff --git a/test/unit/math/mix/functor/root_finder_test.cpp b/test/unit/math/mix/functor/root_finder_test.cpp index 8e51d316747..ce3242bede5 100644 --- a/test/unit/math/mix/functor/root_finder_test.cpp +++ b/test/unit/math/mix/functor/root_finder_test.cpp @@ -1,5 +1,24 @@ #include -/* + +struct CubedRootFinder { + template * = nullptr> + static auto run(T1&& g, T2&& x) { + auto g_pow = g; + auto second_deriv = 20 * g_pow; + g_pow *= g; + auto first_deriv = 3 * g_pow; + g_pow *= g; + auto func_val = g_pow - x; + return std::make_tuple(func_val, first_deriv, second_deriv); + } + template * = nullptr> + static auto run(T1&& g, T2&& x) { + return g * g * g - x; + } +}; + TEST(MixFun, root_finder_cubed) { using stan::math::root_finder; // return cube root of x using 1st and 2nd derivatives and Halley. @@ -18,15 +37,12 @@ TEST(MixFun, root_finder_cubed) { const int digits = std::numeric_limits::digits; int get_digits = static_cast(digits * 0.4); std::uintmax_t maxit = 20; - auto f = [](const auto& g, const auto& x) { return g * g * g - x; }; - auto f_g = [](const auto& g, const auto& x) { return 3 * g * g; }; - auto f_gg = [](const auto& g, const auto& x) { return 6 * g; }; - auto full_f = [&f, &f_g, &f_gg, guess, min, max](auto&& xx) { - return root_finder(std::make_tuple(f, f_g, f_gg), guess, min, max, xx); + auto full_f = [guess, min, max](auto&& xx) { + return stan::math::root_finder_hailey(guess, min, max, xx); }; stan::test::expect_ad(full_f, x); } -*/ + struct FifthRootFinder { @@ -72,37 +88,22 @@ TEST(MixFun, root_finder_fifth) { stan::test::expect_ad(f_hailey, x); } -/* struct BetaCdfRoot { template * = nullptr> static auto run(T1&& x, T2&& alpha, T3&& beta, T4&& p) { - auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - return stan::math::beta_cdf(x, alpha, beta) - p; - }(x, alpha, beta, p); - auto beta_pdf_deriv = [](auto&& x, auto&& alpha, auto&& beta) { - using stan::math::lgamma; - using stan::math::pow; - auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) - * (1. - 2. * x - alpha + x * alpha + x * beta) - * lgamma(alpha + beta)) - / (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); - return val; - }(x, alpha, beta); - auto beta_pdf = [](auto&& x, auto&& alpha, auto&& beta) { - using stan::math::lgamma; - using stan::math::pow; - auto val - = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta)) - / ((lgamma(alpha) * lgamma(beta))); - return val; - }(x, alpha, beta); - return std::make_tuple(f, beta_pdf, beta_pdf_deriv); + auto f_val = stan::math::inc_beta(alpha, beta, x) - p; + auto beta_ab = stan::math::beta(alpha, beta); + auto f_val_d = (stan::math::pow(1.0 - x,-1.0 + beta)*stan::math::pow(x,-1.0 + alpha)) / beta_ab; + auto f_val_dd = (stan::math::pow(1 - x,-1 + beta)*stan::math::pow(x,-2 + alpha)*(-1 + alpha))/beta_ab - + (stan::math::pow(1 - x,-2 + beta)*stan::math::pow(x,-1 + alpha)*(-1 + beta))/beta_ab; + + return std::make_tuple(f_val, f_val_d, f_val_dd); } template * = nullptr> static auto run(T1&& x, T2&& alpha, T3&& beta, T4&& p) { - return stan::math::beta_cdf(x, alpha, beta) - p; + return stan::math::inc_beta(alpha, beta, x) - p; } }; @@ -119,4 +120,3 @@ TEST(MixFun, root_finder_beta) { double p = .3; stan::test::expect_ad(f_hailey, alpha, beta, p); } -*/ diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index 7aafe239cbd..0bfa096c47a 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -5,66 +5,24 @@ #include #include #include -/* -TEST(RevFunctor, root_finder) { - using stan::math::root_finder; - using stan::math::var; - // return cube root of x using 1st and 2nd derivatives and Halley. - // using namespace std; // Help ADL of std functions. - var x = 27; - int exponent; - // Get exponent of z (ignore mantissa). - std::frexp(x.val(), &exponent); - // Rough guess is to divide the exponent by three. - var guess = ldexp(1., exponent / 3); - // Minimum possible value is half our guess. - var min = ldexp(0.5, exponent / 3); - // Maximum possible value is twice our guess. - var max = ldexp(2., exponent / 3); - // Maximum possible binary digits accuracy for type T. - const int digits = std::numeric_limits::digits; - int get_digits = static_cast(digits * 0.4); - std::uintmax_t maxit = 20; - auto f = [](const auto& g, const auto x) { return g * g * g - x; }; - auto f_g = [](const auto& g, const auto x) { return 3 * g * g; }; - auto f_gg = [](const auto& g, const auto x) { return 6 * g; }; - var result = root_finder(std::make_tuple(f, f_g, f_gg), guess, min, max, x); - result.grad(); - EXPECT_EQ(27, x.val()); - EXPECT_NEAR(0.037037, x.adj(), 1e-6); -} -*/ + struct BetaCdfRoot { template * = nullptr> static auto run(T1&& x, T2&& alpha, T3&& beta, T4&& p) { - auto f = [](auto&& x, auto&& alpha, auto&& beta, auto&& p) { - return stan::math::beta_cdf(x, alpha, beta) - p; - }(x, alpha, beta, p); - auto beta_pdf_deriv = [](auto&& x, auto&& alpha, auto&& beta) { - using stan::math::lgamma; - using stan::math::pow; - auto val = -(pow((1 - x), beta) * pow(x, (-2 + alpha)) - * (1. - 2. * x - alpha + x * alpha + x * beta) - * lgamma(alpha + beta)) - / (pow((-1 + x), 2) * lgamma(alpha) * lgamma(beta)); - return val; - }(x, alpha, beta); - auto beta_pdf = [](auto&& x, auto&& alpha, auto&& beta) { - using stan::math::lgamma; - using stan::math::pow; - auto val - = (pow(x, alpha - 1) * pow((1 - x), beta - 1) * lgamma(alpha + beta)) - / ((lgamma(alpha) * lgamma(beta))); - return val; - }(x, alpha, beta); - return std::make_tuple(f, beta_pdf, beta_pdf_deriv); + auto f_val = boost::math::ibeta(alpha, beta, x) - p; + auto beta_ab = boost::math::beta(alpha, beta); + double f_val_d = (std::pow(1.0 - x,-1.0 + beta)*std::pow(x,-1.0 + alpha)) / beta_ab; + double f_val_dd = (std::pow(1 - x,-1 + beta)*std::pow(x,-2 + alpha)*(-1 + alpha))/beta_ab - + (std::pow(1 - x,-2 + beta)*std::pow(x,-1 + alpha)*(-1 + beta))/beta_ab; + + return std::make_tuple(f_val, f_val_d, f_val_dd); } template * = nullptr> static auto run(T1&& x, T2&& alpha, T3&& beta, T4&& p) { - return stan::math::beta_cdf(x, alpha, beta) - p; + return stan::math::inc_beta(alpha, beta, x) - p; } }; TEST(RevFunctor, root_finder_beta_cdf) { @@ -293,6 +251,8 @@ TEST(RevFunctor, root_finder_beta_cdf2) { alpha, beta, p); }; auto f_schroder = [guess, min, max](auto&& alpha, auto&& beta, auto&& p) { + constexpr int digits = 16; + std::uintmax_t max_iter = std::numeric_limits::max(); return stan::math::root_finder_schroder(guess, min, max, alpha, beta, p); }; @@ -303,18 +263,20 @@ TEST(RevFunctor, root_finder_beta_cdf2) { check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), f_schroder, 5e-2, .5, .5, .3); - /* For some reason this fails after a while?? - for (double p = .3; p < .7; p += .1) { - for (double a = .3; a < .7; a += .1) { - for (double b = .3; b < .7; b += .1) { + check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), f_schroder, + 5e-2, .5, .6, .5); + // For some reason this fails after a while?? + for (double p = .1; p < .9; p += .1) { + for (double a = .1; a < .9; a += .1) { + for (double b = .1; b < .9; b += .1) { check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), - f_schroder, 5e-2, a, b, p); + f_schroder, 1e-9, a, b, p); if (::testing::Test::HasFailure()) { std::cout << "--\na: " << a << "\nb: " << b << "\np: " << p << "\n"; } } } } - */ + } From 8cd506b68ed30b2f3434cde14224ea86be802ab4 Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Fri, 3 Jun 2022 17:47:17 -0400 Subject: [PATCH 20/20] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/functor/root_finder.hpp | 33 ++++++++-------- stan/math/rev/functor/root_finder.hpp | 39 ++++++++++--------- .../math/mix/functor/root_finder_test.cpp | 22 ++++++----- .../math/rev/functor/root_finder_test.cpp | 27 +++++++------ 4 files changed, 65 insertions(+), 56 deletions(-) diff --git a/stan/math/prim/functor/root_finder.hpp b/stan/math/prim/functor/root_finder.hpp index 1102f1f9bb4..81e14456305 100644 --- a/stan/math/prim/functor/root_finder.hpp +++ b/stan/math/prim/functor/root_finder.hpp @@ -82,8 +82,8 @@ template * = nullptr> inline auto root_finder_tol(const GuessScalar guess, const MinScalar min, - const MaxScalar max, const int digits, - std::uintmax_t& max_iter, Types&&... args) { + const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); check_positive("root_finder", "digits", digits); check_positive("root_finder", "max_iter", max_iter); @@ -103,27 +103,27 @@ inline auto root_finder_tol(const GuessScalar guess, const MinScalar min, template inline auto root_finder_halley_tol(const GuessScalar guess, const MinScalar min, - const MaxScalar max, const int digits, - std::uintmax_t& max_iter, Types&&... args) { + const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types&&... args) { return root_finder_tol( guess, min, max, digits, max_iter, std::forward(args)...); } template -inline auto root_finder_newton_raphson_tol(const GuessScalar guess, - const MinScalar min, const MaxScalar max, - const int digits, std::uintmax_t& max_iter, - Types&&... args) { +inline auto root_finder_newton_raphson_tol( + const GuessScalar guess, const MinScalar min, const MaxScalar max, + const int digits, std::uintmax_t& max_iter, Types&&... args) { return root_finder_tol( guess, min, max, digits, max_iter, std::forward(args)...); } template -inline auto root_finder_schroder_tol(const GuessScalar guess, const MinScalar min, - const MaxScalar max, const int digits, - std::uintmax_t& max_iter, Types&&... args) { +inline auto root_finder_schroder_tol(const GuessScalar guess, + const MinScalar min, const MaxScalar max, + const int digits, std::uintmax_t& max_iter, + Types&&... args) { return root_finder_tol( guess, min, max, digits, max_iter, std::forward(args)...); } @@ -152,7 +152,7 @@ inline auto root_finder_schroder_tol(const GuessScalar guess, const MinScalar mi template inline auto root_finder(const GuessScalar guess, const MinScalar min, - const MaxScalar max, Types&&... args) { + const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); return root_finder_tol( @@ -162,7 +162,7 @@ inline auto root_finder(const GuessScalar guess, const MinScalar min, template inline auto root_finder_hailey(const GuessScalar guess, const MinScalar min, - const MaxScalar max, Types&&... args) { + const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); return root_finder_halley_tol(guess, min, max, digits, max_iter, @@ -171,8 +171,9 @@ inline auto root_finder_hailey(const GuessScalar guess, const MinScalar min, template -inline auto root_finder_newton_raphson(const GuessScalar guess, const MinScalar min, - const MaxScalar max, Types&&... args) { +inline auto root_finder_newton_raphson(const GuessScalar guess, + const MinScalar min, const MaxScalar max, + Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); return root_finder_newton_raphson_tol( @@ -182,7 +183,7 @@ inline auto root_finder_newton_raphson(const GuessScalar guess, const MinScalar template inline auto root_finder_schroder(const GuessScalar guess, const MinScalar min, - const MaxScalar max, Types&&... args) { + const MaxScalar max, Types&&... args) { constexpr int digits = 16; std::uintmax_t max_iter = std::numeric_limits::max(); return root_finder_schroder_tol(guess, min, max, digits, max_iter, diff --git a/stan/math/rev/functor/root_finder.hpp b/stan/math/rev/functor/root_finder.hpp index 3966bb0ecb3..7de5124df05 100644 --- a/stan/math/rev/functor/root_finder.hpp +++ b/stan/math/rev/functor/root_finder.hpp @@ -44,8 +44,8 @@ template < require_all_stan_scalar_t* = nullptr, require_any_not_stan_scalar_t* = nullptr> inline auto root_finder_tol(const GuessScalar guess, const MinScalar min, - const MaxScalar max, const int digits, - std::uintmax_t& max_iter, Types&&... args) { + const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types&&... args) { check_bounded("root_finder", "initial guess", guess, min, max); check_positive("root_finder", "digits", digits); check_positive("root_finder", "max_iter", max_iter); @@ -61,7 +61,8 @@ inline auto root_finder_tol(const GuessScalar guess, const MinScalar min, [&max_iter, digits, guess_val = value_of(guess), min_val = value_of(min), max_val = value_of(max)](auto&&... vals) { return root_finder_tol( - guess_val, min_val, max_val, digits > 20 ? digits : 21, max_iter, vals...); + guess_val, min_val, max_val, digits > 20 ? digits : 21, max_iter, + vals...); }, args_vals_tuple); double Jf_x; @@ -132,14 +133,15 @@ template < require_all_stan_scalar_t* = nullptr, require_all_stan_scalar_t* = nullptr> inline auto root_finder_tol(const GuessScalar guess, const MinScalar min, - const MaxScalar max, const int digits, - std::uintmax_t& max_iter, Types... args) { + const MaxScalar max, const int digits, + std::uintmax_t& max_iter, Types... args) { check_bounded("root_finder", "initial guess", guess, min, max); check_positive("root_finder", "digits", digits); check_positive("root_finder", "max_iter", max_iter); // Solve the system double theta_dbl = root_finder_tol( - value_of(guess), value_of(min), value_of(max), digits > 20 ? digits : 21, max_iter, value_of(args)...); + value_of(guess), value_of(min), value_of(max), digits > 20 ? digits : 21, + max_iter, value_of(args)...); double Jf_x; { nested_rev_autodiff nested; @@ -154,19 +156,18 @@ inline auto root_finder_tol(const GuessScalar guess, const MinScalar min, * Note: Because we put this on the callback stack, if `f` is a lambda * its captures must be in Stan's arena memory or trivially destructable. */ - return make_callback_var( - theta_dbl, [Jf_x, args...](auto& ret) mutable { - { - nested_rev_autodiff rev; - double eta = -(ret.adj() / Jf_x); - double ret_val = ret.val(); - auto f = internal::make_root_func(); - auto x_nrad = std::decay_t::template run( - ret_val, args...); - x_nrad.adj() = eta; - grad(); - } - }); + return make_callback_var(theta_dbl, [Jf_x, args...](auto& ret) mutable { + { + nested_rev_autodiff rev; + double eta = -(ret.adj() / Jf_x); + double ret_val = ret.val(); + auto f = internal::make_root_func(); + auto x_nrad + = std::decay_t::template run(ret_val, args...); + x_nrad.adj() = eta; + grad(); + } + }); } } // namespace math diff --git a/test/unit/math/mix/functor/root_finder_test.cpp b/test/unit/math/mix/functor/root_finder_test.cpp index ce3242bede5..a87f491072e 100644 --- a/test/unit/math/mix/functor/root_finder_test.cpp +++ b/test/unit/math/mix/functor/root_finder_test.cpp @@ -2,7 +2,7 @@ struct CubedRootFinder { template * = nullptr> + std::enable_if_t* = nullptr> static auto run(T1&& g, T2&& x) { auto g_pow = g; auto second_deriv = 20 * g_pow; @@ -13,7 +13,7 @@ struct CubedRootFinder { return std::make_tuple(func_val, first_deriv, second_deriv); } template * = nullptr> + std::enable_if_t* = nullptr> static auto run(T1&& g, T2&& x) { return g * g * g - x; } @@ -43,11 +43,9 @@ TEST(MixFun, root_finder_cubed) { stan::test::expect_ad(full_f, x); } - - struct FifthRootFinder { template * = nullptr> + std::enable_if_t* = nullptr> static auto run(T1&& g, T2&& x) { auto g_pow = g * g * g; auto second_deriv = 20 * g_pow; @@ -58,7 +56,7 @@ struct FifthRootFinder { return std::make_tuple(func_val, first_deriv, second_deriv); } template * = nullptr> + std::enable_if_t* = nullptr> static auto run(T1&& g, T2&& x) { return g * g * g * g * g - x; } @@ -94,9 +92,15 @@ struct BetaCdfRoot { static auto run(T1&& x, T2&& alpha, T3&& beta, T4&& p) { auto f_val = stan::math::inc_beta(alpha, beta, x) - p; auto beta_ab = stan::math::beta(alpha, beta); - auto f_val_d = (stan::math::pow(1.0 - x,-1.0 + beta)*stan::math::pow(x,-1.0 + alpha)) / beta_ab; - auto f_val_dd = (stan::math::pow(1 - x,-1 + beta)*stan::math::pow(x,-2 + alpha)*(-1 + alpha))/beta_ab - - (stan::math::pow(1 - x,-2 + beta)*stan::math::pow(x,-1 + alpha)*(-1 + beta))/beta_ab; + auto f_val_d = (stan::math::pow(1.0 - x, -1.0 + beta) + * stan::math::pow(x, -1.0 + alpha)) + / beta_ab; + auto f_val_dd = (stan::math::pow(1 - x, -1 + beta) + * stan::math::pow(x, -2 + alpha) * (-1 + alpha)) + / beta_ab + - (stan::math::pow(1 - x, -2 + beta) + * stan::math::pow(x, -1 + alpha) * (-1 + beta)) + / beta_ab; return std::make_tuple(f_val, f_val_d, f_val_dd); } diff --git a/test/unit/math/rev/functor/root_finder_test.cpp b/test/unit/math/rev/functor/root_finder_test.cpp index 0bfa096c47a..971fc9eefd1 100644 --- a/test/unit/math/rev/functor/root_finder_test.cpp +++ b/test/unit/math/rev/functor/root_finder_test.cpp @@ -6,16 +6,20 @@ #include #include - struct BetaCdfRoot { template * = nullptr> static auto run(T1&& x, T2&& alpha, T3&& beta, T4&& p) { auto f_val = boost::math::ibeta(alpha, beta, x) - p; auto beta_ab = boost::math::beta(alpha, beta); - double f_val_d = (std::pow(1.0 - x,-1.0 + beta)*std::pow(x,-1.0 + alpha)) / beta_ab; - double f_val_dd = (std::pow(1 - x,-1 + beta)*std::pow(x,-2 + alpha)*(-1 + alpha))/beta_ab - - (std::pow(1 - x,-2 + beta)*std::pow(x,-1 + alpha)*(-1 + beta))/beta_ab; + double f_val_d + = (std::pow(1.0 - x, -1.0 + beta) * std::pow(x, -1.0 + alpha)) + / beta_ab; + double f_val_dd + = (std::pow(1 - x, -1 + beta) * std::pow(x, -2 + alpha) * (-1 + alpha)) + / beta_ab + - (std::pow(1 - x, -2 + beta) * std::pow(x, -1 + alpha) * (-1 + beta)) + / beta_ab; return std::make_tuple(f_val, f_val_d, f_val_dd); } @@ -263,20 +267,19 @@ TEST(RevFunctor, root_finder_beta_cdf2) { check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), f_schroder, 5e-2, .5, .5, .3); - check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), f_schroder, - 5e-2, .5, .6, .5); + check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), f_schroder, + 5e-2, .5, .6, .5); // For some reason this fails after a while?? for (double p = .1; p < .9; p += .1) { for (double a = .1; a < .9; a += .1) { for (double b = .1; b < .9; b += .1) { - check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), - f_schroder, 1e-9, a, b, p); - if (::testing::Test::HasFailure()) { - std::cout << "--\na: " << a << "\nb: " << b << "\np: " << p << "\n"; - } + check_vs_known_grads(std::make_tuple(deriv_a, deriv_b, deriv_p), + f_schroder, 1e-9, a, b, p); + if (::testing::Test::HasFailure()) { + std::cout << "--\na: " << a << "\nb: " << b << "\np: " << p << "\n"; + } } } } - }