diff --git a/stan/math/rev/fun/log_sum_exp.hpp b/stan/math/rev/fun/log_sum_exp.hpp index 9ca331aab74..4ef0ece4af3 100644 --- a/stan/math/rev/fun/log_sum_exp.hpp +++ b/stan/math/rev/fun/log_sum_exp.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -22,8 +23,10 @@ class log_sum_exp_vv_vari : public op_vv_vari { log_sum_exp_vv_vari(vari* avi, vari* bvi) : op_vv_vari(log_sum_exp(avi->val_, bvi->val_), avi, bvi) {} void chain() { - avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_); - bvi_->adj_ += adj_ * calculate_chain(bvi_->val_, val_); + // avi_->adj_ += adj_ * calculate_chain(avi_->val_, val_); + // bvi_->adj_ += adj_ * calculate_chain(bvi_->val_, val_); + avi_->adj_ += adj_ * inv_logit(avi_->val_ - bvi_->val_); + bvi_->adj_ += adj_ * inv_logit(bvi_->val_ - avi_->val_); } }; class log_sum_exp_vd_vari : public op_vd_vari { diff --git a/stan/math/rev/functor.hpp b/stan/math/rev/functor.hpp index 0af65fe88e5..348fb087c02 100644 --- a/stan/math/rev/functor.hpp +++ b/stan/math/rev/functor.hpp @@ -20,5 +20,6 @@ #include #include #include +#include #endif diff --git a/stan/math/rev/functor/pack_params.hpp b/stan/math/rev/functor/pack_params.hpp new file mode 100644 index 00000000000..dbae8d788ed --- /dev/null +++ b/stan/math/rev/functor/pack_params.hpp @@ -0,0 +1,47 @@ +#ifndef STAN_MATH_REV_FUNCTOR_PACK_PARAMS_HPP +#define STAN_MATH_REV_FUNCTOR_PACK_PARAMS_HPP + +#include +#include +#include +#include + +namespace stan { +namespace math { + +template +class invoke_packed { + public: + static var invoke(const F& f, + const Eigen::Matrix& x_arr, + Ts... xs) { + return invoke_packed::invoke( + f, x_arr, x_arr(x_arr.size() - N), xs...); + } +}; + +template +class invoke_packed<0, F, Ts...> { + public: + static var invoke(const F& f, const Eigen::Matrix&, + Ts... xs) { + return f(xs...); + } +}; + +template +class pack_params { + public: + explicit pack_params(const F& f) : f_(f) {} + + var operator()(const Eigen::Matrix& x) const { + return invoke_packed::invoke(f_, x); + } + + private: + const F& f_; +}; + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/ad_tolerances.hpp b/test/unit/math/ad_tolerances.hpp index 5c7eca72bc1..c34a3dd412c 100644 --- a/test/unit/math/ad_tolerances.hpp +++ b/test/unit/math/ad_tolerances.hpp @@ -1,6 +1,8 @@ #ifndef TEST_UNIT_MATH_AD_TOLERANCES_HPP #define TEST_UNIT_MATH_AD_TOLERANCES_HPP +#include + namespace stan { namespace test { @@ -13,34 +15,35 @@ namespace test { * if the function is implemented using only forward mode, and end * with the quantity being calculated; for example, * `hessian_fvar_grad_` is the gradient calculated by the `hessian` - * function using forward-mode autodiff. + * function using forward-mode autodiff. Those get interpreted as + * relative tolerances, see `relative_tolerance` class more details. * * `gradient_val_`: 1e-8; `gradient_grad_`: 1e-4 * * `gradient_fvar_val_`: 1e-8; `gradient_fvar_grad_`: 1e-4 * - * `hessian_val_` : 1e-8; `hessian_grad_`: 1e-4; `hessian_hessian_`: 1e-3 + * `hessian_val_` : 1e-8; `hessian_grad_`: 1e-4; `hessian_hessian_`: (1e-4,1e-3) * * `hessian_fvar_val_` : 1e-8; `hessian_fvar_grad_`: 1e-4; - * `hessian_fvar_hessian_`: 1e-3 + * `hessian_fvar_hessian_`: (1e-4,1e-3) * * `grad_hessian_val_` : 1e-8; `grad_hessian_hessian_`: 1e-3; * `grad_hessian_grad_hessian_`: 1e-2 */ struct ad_tolerances { - double gradient_val_; - double gradient_grad_; - double gradient_fvar_val_; - double gradient_fvar_grad_; - double hessian_val_; - double hessian_grad_; - double hessian_hessian_; - double hessian_fvar_val_; - double hessian_fvar_grad_; - double hessian_fvar_hessian_; - double grad_hessian_val_; - double grad_hessian_hessian_; - double grad_hessian_grad_hessian_; + relative_tolerance gradient_val_; + relative_tolerance gradient_grad_; + relative_tolerance gradient_fvar_val_; + relative_tolerance gradient_fvar_grad_; + relative_tolerance hessian_val_; + relative_tolerance hessian_grad_; + relative_tolerance hessian_hessian_; + relative_tolerance hessian_fvar_val_; + relative_tolerance hessian_fvar_grad_; + relative_tolerance hessian_fvar_hessian_; + relative_tolerance grad_hessian_val_; + relative_tolerance grad_hessian_hessian_; + relative_tolerance grad_hessian_grad_hessian_; ad_tolerances() : gradient_val_(1e-8), gradient_grad_(1e-4), @@ -50,11 +53,11 @@ struct ad_tolerances { hessian_val_(1e-8), hessian_grad_(1e-4), - hessian_hessian_(1e-3), + hessian_hessian_(1e-4, 1e-3), hessian_fvar_val_(1e-8), hessian_fvar_grad_(1e-4), - hessian_fvar_hessian_(1e-3), + hessian_fvar_hessian_(1e-4, 1e-3), grad_hessian_val_(1e-8), grad_hessian_hessian_(1e-3), diff --git a/test/unit/math/expect_near_rel.hpp b/test/unit/math/expect_near_rel.hpp index 37d12c266cb..9f7f7ae0696 100644 --- a/test/unit/math/expect_near_rel.hpp +++ b/test/unit/math/expect_near_rel.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -15,13 +16,7 @@ namespace internal { * Test that the specified values are within the specified tolerance * on relative error, and if not, fail the embedded google test. * - *

Relative error is defined to be the error `u - v` rescaled by the - * average absolute value, - * `rel_err(u, v) = (u - v) / (0.5 * (abs(u) * + abs(v))).` - * - *

If at least one of `u` or `v` is zero, the absolute error is - * tested at the specified tolerance, because the relative error - * reduces to a constant. + * Uses relative_tolerance::inexact to compute the tolerance. * * @tparam T1 type of first argument * @tparam T2 type of second argument @@ -32,23 +27,11 @@ namespace internal { */ template ...> void expect_near_rel_finite(const std::string& msg, const T1& x1, const T2& x2, - double tol = 1e-8) { - using stan::math::fabs; - // if either arg near zero, must use absolute tolerance as rel tol -> 2 - if (fabs(x1) < tol || fabs(x2) < tol) { - EXPECT_NEAR(x1, x2, tol) << "expect_near_rel_finite(" << x1 << ", " << x2 - << ", absolute tolerance = " << tol << ")" - << "; absolute diff = " << fabs(x1 - x2) - << " in: " << msg << std::endl; - return; - } - auto avg = 0.5 * (fabs(x1) + fabs(x2)); - auto relative_diff = (x1 - x2) / avg; - EXPECT_NEAR(0, relative_diff, tol) - << "expect_near_rel_finite(" << x1 << ", " << x2 - << ", relative tolerance = " << tol << ")" - << "; relative diff = " << relative_diff << std::endl - << " in: " << msg << std::endl; + const relative_tolerance tol + = relative_tolerance()) { + double tol_val = tol.inexact(x1, x2); + EXPECT_NEAR(x1, x2, tol_val) + << "expect_near_rel_finite in: " << msg << std::endl; } template & x1, */ template ...> void expect_near_rel(const std::string& msg, const T1& x1, const T2& x2, - double tol = 1e-8) { + relative_tolerance tol = relative_tolerance()) { if (stan::math::is_nan(x1) || stan::math::is_nan(x2)) EXPECT_TRUE(stan::math::is_nan(x1) && stan::math::is_nan(x2)) << "expect_near_rel(" << x1 << ", " << x2 << ")" << std::endl @@ -120,7 +103,7 @@ void expect_near_rel(const std::string& msg, const T1& x1, const T2& x2, template ...> void expect_near_rel(const std::string& msg, EigMat1&& x1, EigMat2&& x2, - double tol = 1e-8) { + relative_tolerance tol = relative_tolerance()) { EXPECT_EQ(x1.rows(), x2.rows()) << "expect_near_rel (Eigen::Matrix)" << " rows must be same size." << " x1.rows() = " << x1.rows() @@ -141,7 +124,8 @@ void expect_near_rel(const std::string& msg, EigMat1&& x1, EigMat2&& x2, template void expect_near_rel(const std::string& msg, const std::vector& x1, - const std::vector& x2, double tol = 1e-8) { + const std::vector& x2, + relative_tolerance tol = relative_tolerance()) { EXPECT_EQ(x1.size(), x2.size()) << "expect_near_rel (std::vector):" << " vectors must be same size." << " x1.size() = " << x1.size() diff --git a/test/unit/math/expect_near_rel_test.cpp b/test/unit/math/expect_near_rel_test.cpp index a07f5195f68..58e778541d9 100644 --- a/test/unit/math/expect_near_rel_test.cpp +++ b/test/unit/math/expect_near_rel_test.cpp @@ -9,10 +9,10 @@ TEST(testUnitMath, ExpectNearRelScalar) { expect_near_rel("test A1", 0, 0, 1e-16); expect_near_rel("test A2", 0, 0.0, 1e-16); expect_near_rel("test A3", 0.0, 0, 1e-16); - expect_near_rel("test B1", 0, 1e-8, 1e-5); - expect_near_rel("test B2", 0.0, 1e-8, 1e-5); - expect_near_rel("test C1", 1e-8, 0, 1e-5); - expect_near_rel("test C2", 1e-8, 0.0, 1e-5); + expect_near_rel("test B1", 0, 1e-8, 1e-4); + expect_near_rel("test B2", 0.0, 1e-8, 1e-4); + expect_near_rel("test C1", 1e-8, 0, 1e-4); + expect_near_rel("test C2", 1e-8, 0.0, 1e-4); // non-zero examples expect_near_rel("test D1", 1, 1, 1e-16); @@ -33,6 +33,7 @@ TEST(testUnitMath, ExpectNearRelScalar) { TEST(testUnitMath, ExpectNearRelMatrix) { using Eigen::Matrix; using stan::test::expect_near_rel; + using stan::test::relative_tolerance; typedef Matrix v_t; typedef Matrix rv_t; typedef Matrix m_t; @@ -57,7 +58,7 @@ TEST(testUnitMath, ExpectNearRelMatrix) { m_t d2(2, 3); d1 << 1, 2, 3, 0, 0, 0 - 1e-8; d2 << 1 + 1e-8, 2 - 1e-8, 3, 0, 0 + 1e-8, 0; - expect_near_rel("test D", d1, d2, 1e-6); + expect_near_rel("test D", d1, d2, relative_tolerance(1e-6, 1e-8)); // these will fail // v_t e1(1); @@ -82,7 +83,7 @@ TEST(testUnitMath, ExpectNearRelVector) { expect_near_rel("test C", v_t{1, 1, 1}, v_t{1 + 1e-8, 1 - 1e-9, 1}, 1e-6); - expect_near_rel("test D", v_t{0, 0, 0}, v_t{0, 0 + 1e-6, 0 - 1e-6}, 1e-4); + expect_near_rel("test D", v_t{0, 0, 0}, v_t{0, 0 + 9e-9, 0 - 9e-9}, 1e-4); // ones after here fail // expect_near_rel("test E", v_t{1}, v_t{1, 2}, 1e-6); @@ -91,6 +92,7 @@ TEST(testUnitMath, ExpectNearRelVector) { TEST(testUnitMath, ExpectNearRelVectorNesting) { using stan::test::expect_near_rel; + using stan::test::relative_tolerance; using std::vector; typedef vector v_t; typedef vector vv_t; @@ -102,19 +104,19 @@ TEST(testUnitMath, ExpectNearRelVectorNesting) { expect_near_rel("test A", vv_t{}, vv_t{}, 1e-10); expect_near_rel("test B", vv_t{v_t{1, 2, 3}, v_t{0, 0, 0}}, - vv_t{v_t{1, 2, 3}, v_t{0, 0 + 1e-6, 0 - 1e-6}}, 1e-5); + vv_t{v_t{1, 2, 3}, v_t{0, 0 + 1e-6, 0 - 1e-6}}, 1e-3); expect_near_rel( "test C", vvv_t{vv_t{v_t{1, 2, 3}, v_t{0, 0, 0}}, vv_t{v_t{1, 2, 3}, v_t{0, 0, 0}}}, vvv_t{vv_t{v_t{1, 2, 3}, v_t{0, 0 + 1e-6, 0 - 1e-6}}, vv_t{v_t{1, 2, 3}, v_t{0, 0 + 1e-6, 0 - 1e-6}}}, - 1e-5); + relative_tolerance(1e-5, 1e-6)); ev_t d1(3); ev_t d2(3); d1 << 1, 0, 0; - d2 << 1 + 1e-8, 0 + 1e-8, 0; + d2 << 1 + 1e-8, 0 + 1e-12, 0; vev_t e1{d1, d2}; vev_t e2{d2, d1}; expect_near_rel("test E", e1, e2, 1e-6); diff --git a/test/unit/math/mix/fun/add_diag_test.cpp b/test/unit/math/mix/fun/add_diag_test.cpp index 728e9da61b5..8602a4ea48b 100644 --- a/test/unit/math/mix/fun/add_diag_test.cpp +++ b/test/unit/math/mix/fun/add_diag_test.cpp @@ -1,6 +1,8 @@ #include +#include TEST(MathMixMatFun, addDiag) { + using stan::test::relative_tolerance; auto f = [](const auto& x, const auto& y) { return stan::math::add_diag(x, y); }; @@ -22,18 +24,21 @@ TEST(MathMixMatFun, addDiag) { Eigen::MatrixXd m23(2, 3); m23 << 1, 1, 1, 1, 1, 1; + stan::test::ad_tolerances tol; + tol.hessian_hessian_ = relative_tolerance(1e-4, 1e-3); + tol.hessian_fvar_hessian_ = relative_tolerance(1e-4, 1e-3); // these are OK calls - stan::test::expect_ad(f, m00, d); - stan::test::expect_ad(f, m00, v0); - stan::test::expect_ad(f, m11, d); - stan::test::expect_ad(f, m11, v1); - stan::test::expect_ad(f, m22, d); - stan::test::expect_ad(f, m22, v2); - stan::test::expect_ad(f, m23, d); - stan::test::expect_ad(f, m23, v2); + stan::test::expect_ad(tol, f, m00, d); + stan::test::expect_ad(tol, f, m00, v0); + stan::test::expect_ad(tol, f, m11, d); + stan::test::expect_ad(tol, f, m11, v1); + stan::test::expect_ad(tol, f, m22, d); + stan::test::expect_ad(tol, f, m22, v2); + stan::test::expect_ad(tol, f, m23, d); + stan::test::expect_ad(tol, f, m23, v2); // these throw - stan::test::expect_ad(f, m11, v2); - stan::test::expect_ad(f, m22, v1); - stan::test::expect_ad(f, m23, v1); + stan::test::expect_ad(tol, f, m11, v2); + stan::test::expect_ad(tol, f, m22, v1); + stan::test::expect_ad(tol, f, m23, v1); } diff --git a/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp b/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp index 0d995e0e898..1520fcb340d 100644 --- a/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp +++ b/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp @@ -33,9 +33,10 @@ typename stan::scalar_type::type g3(const T& x) { template void expect_cov_matrix_transform(const T& x) { + using stan::test::relative_tolerance; stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; + tols.hessian_hessian_ = relative_tolerance(1e-3, 1e-3); + tols.hessian_fvar_hessian_ = relative_tolerance(1e-3, 1e-3); auto f1 = [](const auto& x) { return g1(x); }; auto f2 = [](const auto& x) { return g2(x); }; diff --git a/test/unit/math/mix/fun/determinant_test.cpp b/test/unit/math/mix/fun/determinant_test.cpp index a384f1999a1..925c8c90c38 100644 --- a/test/unit/math/mix/fun/determinant_test.cpp +++ b/test/unit/math/mix/fun/determinant_test.cpp @@ -3,11 +3,12 @@ #include TEST(MathMixMatFun, determinant) { + using stan::test::relative_tolerance; auto f = [](const auto& y) { return stan::math::determinant(y); }; stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; // default 1e-3 - tols.hessian_fvar_hessian_ = 1e-2; // default 1e-3 + tols.hessian_hessian_ = relative_tolerance(1e-4, 1e-2); + tols.hessian_fvar_hessian_ = relative_tolerance(1e-4, 1e-2); Eigen::MatrixXd z(0, 0); diff --git a/test/unit/math/mix/fun/diag_post_multiply_test.cpp b/test/unit/math/mix/fun/diag_post_multiply_test.cpp index 04d6d4fd737..5350ab72332 100644 --- a/test/unit/math/mix/fun/diag_post_multiply_test.cpp +++ b/test/unit/math/mix/fun/diag_post_multiply_test.cpp @@ -17,6 +17,7 @@ void expect_diag_post_multiply(const Eigen::MatrixXd& a, } TEST(MathMixMatFun, diagPostMultiply) { + using stan::test::relative_tolerance; // 0 x 0 Eigen::MatrixXd a00(0, 0); Eigen::VectorXd u0(0); @@ -50,8 +51,8 @@ TEST(MathMixMatFun, diagPostMultiply) { expect_diag_post_multiply(a33c, u3c); stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 2e-2; // default 1e-3 - tols.hessian_fvar_hessian_ = 2e-2; // default 1e-3 + tols.hessian_hessian_ = relative_tolerance(1e-4, 2e-2); + tols.hessian_fvar_hessian_ = relative_tolerance(1e-4, 2e-2); Eigen::MatrixXd a33(3, 3); a33 << 1, 10, 100, 1000, 2, -4, 8, -16, 32; diff --git a/test/unit/math/mix/fun/diag_pre_multiply_test.cpp b/test/unit/math/mix/fun/diag_pre_multiply_test.cpp index 5bf716e0825..5cbe90c99b6 100644 --- a/test/unit/math/mix/fun/diag_pre_multiply_test.cpp +++ b/test/unit/math/mix/fun/diag_pre_multiply_test.cpp @@ -16,6 +16,7 @@ void expect_diag_pre_multiply(const Eigen::VectorXd& v, expect_diag_pre_multiply(v, a, tols); } TEST(MathMixMatFun, diagPreMultiply) { + using stan::test::relative_tolerance; // 0 x 0 Eigen::MatrixXd a00(0, 0); Eigen::VectorXd u0(0); @@ -49,8 +50,8 @@ TEST(MathMixMatFun, diagPreMultiply) { expect_diag_pre_multiply(u3c, a33c); stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 2e-2; // default 1e-3 - tols.hessian_fvar_hessian_ = 2e-2; // default 1e-3 + tols.hessian_hessian_ = relative_tolerance(1e-4, 2e-2); + tols.hessian_fvar_hessian_ = relative_tolerance(1e-4, 2e-2); Eigen::MatrixXd a33(3, 3); a33 << 1, 10, 100, 1000, 2, -4, 8, -16, 32; diff --git a/test/unit/math/mix/fun/dot_product_test.cpp b/test/unit/math/mix/fun/dot_product_test.cpp index 3b1f40be7bd..301a691a6e1 100644 --- a/test/unit/math/mix/fun/dot_product_test.cpp +++ b/test/unit/math/mix/fun/dot_product_test.cpp @@ -4,9 +4,10 @@ template void test_dot_product(const F& f) { + using stan::test::relative_tolerance; stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; // default is 1e-3 - tols.hessian_fvar_hessian_ = 1e-2; // default is 1e-3 + tols.hessian_hessian_ = relative_tolerance(1e-3, 1e-3); + tols.hessian_fvar_hessian_ = relative_tolerance(1e-3, 1e-3); // size 0 Eigen::VectorXd v0(0); diff --git a/test/unit/math/mix/fun/eigenvalues_sym_test.cpp b/test/unit/math/mix/fun/eigenvalues_sym_test.cpp index a8eb80cd1c6..71913b70f0e 100644 --- a/test/unit/math/mix/fun/eigenvalues_sym_test.cpp +++ b/test/unit/math/mix/fun/eigenvalues_sym_test.cpp @@ -1,6 +1,7 @@ #include TEST(MathMixMatFun, eigenvaluesSym) { + using stan::test::relative_tolerance; auto f = [](const auto& y) { // need to maintain symmetry for finite diffs auto a = ((y + y.transpose()) * 0.5).eval(); @@ -35,6 +36,6 @@ TEST(MathMixMatFun, eigenvaluesSym) { a22 << 1, 2, 2, 1; tols.hessian_hessian_ = 5e-1; tols.hessian_fvar_hessian_ = 5e-1; - tols.grad_hessian_grad_hessian_ = 5e-1; + tols.grad_hessian_grad_hessian_ = relative_tolerance(1e-1, 5e-1); stan::test::expect_ad(tols, f, a22); } diff --git a/test/unit/math/mix/fun/matrix_exp_pade_test.cpp b/test/unit/math/mix/fun/matrix_exp_pade_test.cpp index e55a4617894..5854823a638 100644 --- a/test/unit/math/mix/fun/matrix_exp_pade_test.cpp +++ b/test/unit/math/mix/fun/matrix_exp_pade_test.cpp @@ -1,6 +1,7 @@ #include TEST(MathMixMatFun, matrixExpPade) { + using stan::test::relative_tolerance; auto f = [](const auto& x) { return stan::math::matrix_exp_pade(x); }; Eigen::MatrixXd m00(0, 0); @@ -19,9 +20,9 @@ TEST(MathMixMatFun, matrixExpPade) { d << -2 * a + 3 * b, 1.5 * a - 1.5 * b, -4 * a + 4 * b, 3 * a - 2 * b; stan::test::expect_ad(f, d); - stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; + stan::test::ad_tolerances tols2; + tols2.hessian_hessian_ = relative_tolerance(5e-4, 1e-3); + tols2.hessian_fvar_hessian_ = relative_tolerance(5e-4, 1e-3); a = -1; b = 2; @@ -30,13 +31,13 @@ TEST(MathMixMatFun, matrixExpPade) { e << -24 * a + 40 * b - 15 * c, 18 * a - 30 * b + 12 * c, 5 * a - 8 * b + 3 * c, 20 * b - 20 * c, -15 * b + 16 * c, -4 * b + 4 * c, -120 * a + 120 * b, 90 * a - 90 * b, 25 * a - 24 * b; - stan::test::expect_ad(f, e); + stan::test::expect_ad(tols2, f, e); // replace original random tests for (const auto& x : stan::test::ar_test_cov_matrices(1, 3, 0.0)) { - stan::test::expect_ad(tols, f, x); + stan::test::expect_ad(tols2, f, x); } for (const auto& x : stan::test::ar_test_cov_matrices(1, 3, 0.9)) { - stan::test::expect_ad(tols, f, x); + stan::test::expect_ad(tols2, f, x); } } diff --git a/test/unit/math/mix/fun/matrix_exp_test.cpp b/test/unit/math/mix/fun/matrix_exp_test.cpp index f86d390bf8b..07eb2d8c723 100644 --- a/test/unit/math/mix/fun/matrix_exp_test.cpp +++ b/test/unit/math/mix/fun/matrix_exp_test.cpp @@ -1,6 +1,7 @@ #include TEST(MathMixMatFun, matrixExp) { + using stan::test::relative_tolerance; auto f = [](const auto& x) { return stan::math::matrix_exp(x); }; Eigen::MatrixXd m00(0, 0); @@ -17,9 +18,8 @@ TEST(MathMixMatFun, matrixExp) { stan::test::expect_ad(f, d); stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; - tols.grad_hessian_grad_hessian_ = 1e1; // bad 3rd order ad around 0 + tols.hessian_hessian_ = relative_tolerance(5e-4, 1e-3); + tols.hessian_fvar_hessian_ = relative_tolerance(5e-4, 1e-3); a = -1; b = 2; @@ -28,18 +28,23 @@ TEST(MathMixMatFun, matrixExp) { e << -24 * a + 40 * b - 15 * c, 18 * a - 30 * b + 12 * c, 5 * a - 8 * b + 3 * c, 20 * b - 20 * c, -15 * b + 16 * c, -4 * b + 4 * c, -120 * a + 120 * b, 90 * a - 90 * b, 25 * a - 24 * b; - stan::test::expect_ad(f, e); + stan::test::expect_ad(tols, f, e); + + stan::test::ad_tolerances tols2; + tols2.hessian_hessian_ = 1e-2; + tols2.hessian_fvar_hessian_ = 1e-2; + tols2.grad_hessian_grad_hessian_ = 1e1; // bad 3rd order ad around 0 // replace original random tests for (const auto& x : stan::test::ar_test_cov_matrices(1, 3, 0.0)) { - stan::test::expect_ad(tols, f, x); + stan::test::expect_ad(tols2, f, x); } for (const auto& x : stan::test::ar_test_cov_matrices(1, 3, 0.9)) { - stan::test::expect_ad(tols, f, x); + stan::test::expect_ad(tols2, f, x); } // special tests vs. matrix_exp_2x2 Eigen::MatrixXd mm(2, 2); mm << -0.999984, 0.511211, -0.736924, -0.0826997; - stan::test::expect_ad(tols, f, mm); + stan::test::expect_ad(tols2, f, mm); } diff --git a/test/unit/math/mix/fun/mdivide_right_ldlt_test.cpp b/test/unit/math/mix/fun/mdivide_right_ldlt_test.cpp index cd2a36c327d..78d7778abb8 100644 --- a/test/unit/math/mix/fun/mdivide_right_ldlt_test.cpp +++ b/test/unit/math/mix/fun/mdivide_right_ldlt_test.cpp @@ -2,6 +2,7 @@ #include TEST(MathMixMatFun, mdivideRightLdlt) { + using stan::test::relative_tolerance; auto f = [](const auto& x, const auto& y) { return stan::math::mdivide_right_ldlt(x, stan::test::ldlt_factor(y)); }; @@ -50,8 +51,8 @@ TEST(MathMixMatFun, mdivideRightLdlt) { rv3 << 1, 2, 3; stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; + tols.hessian_hessian_ = relative_tolerance(2e-3, 2e-4); + tols.hessian_fvar_hessian_ = relative_tolerance(2e-3, 2e-4); Eigen::RowVectorXd u(5); u << 62, 84, 84, 76, 108; Eigen::MatrixXd v(5, 5); diff --git a/test/unit/math/mix/fun/mdivide_right_test.cpp b/test/unit/math/mix/fun/mdivide_right_test.cpp index 7232fc2a38d..0a87a830076 100644 --- a/test/unit/math/mix/fun/mdivide_right_test.cpp +++ b/test/unit/math/mix/fun/mdivide_right_test.cpp @@ -5,6 +5,7 @@ #include TEST(MathMixMatFun, mdivideRight) { + using stan::test::relative_tolerance; auto f = [](const auto& x, const auto& y) { return stan::math::mdivide_right(x, y); }; @@ -61,8 +62,8 @@ TEST(MathMixMatFun, mdivideRight) { rv3 << 1, 2, 3; stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; + tols.hessian_hessian_ = relative_tolerance(1e-2, 1e-3); + tols.hessian_fvar_hessian_ = relative_tolerance(1e-2, 1e-3); Eigen::RowVectorXd u(5); u << 62, 84, 84, 76, 108; Eigen::MatrixXd v(5, 5); diff --git a/test/unit/math/mix/fun/quad_form_diag_test.cpp b/test/unit/math/mix/fun/quad_form_diag_test.cpp index bdcf2dc7898..5426d96b387 100644 --- a/test/unit/math/mix/fun/quad_form_diag_test.cpp +++ b/test/unit/math/mix/fun/quad_form_diag_test.cpp @@ -1,6 +1,8 @@ #include TEST(MathMixMatFun, quadFormDiag) { + using stan::test::relative_tolerance; + auto f = [](const auto& x, const auto& y) { return stan::math::quad_form_diag(x, y); }; @@ -22,13 +24,17 @@ TEST(MathMixMatFun, quadFormDiag) { Eigen::VectorXd v3(3); v3 << 1, 2, 3; + stan::test::ad_tolerances tols; + tols.hessian_hessian_ = relative_tolerance(5e-4, 1e-3); + tols.hessian_fvar_hessian_ = relative_tolerance(5e-4, 1e-3); + // matched sizes - stan::test::expect_ad(f, m00, v0); - stan::test::expect_ad(f, m11, v1); - stan::test::expect_ad(f, m22, v2); - stan::test::expect_ad(f, m33, v3); + stan::test::expect_ad(tols, f, m00, v0); + stan::test::expect_ad(tols, f, m11, v1); + stan::test::expect_ad(tols, f, m22, v2); + stan::test::expect_ad(tols, f, m33, v3); // exceptions from mismached sizes - stan::test::expect_ad(f, m33, v2); - stan::test::expect_ad(f, m22, v3); + stan::test::expect_ad(tols, f, m33, v2); + stan::test::expect_ad(tols, f, m22, v3); } diff --git a/test/unit/math/mix/fun/scale_matrix_exp_multiply_test.cpp b/test/unit/math/mix/fun/scale_matrix_exp_multiply_test.cpp index 96cd953564b..f8b03a30bac 100644 --- a/test/unit/math/mix/fun/scale_matrix_exp_multiply_test.cpp +++ b/test/unit/math/mix/fun/scale_matrix_exp_multiply_test.cpp @@ -1,6 +1,7 @@ #include TEST(mathMixMatFun, scaleMatrixExpMultiply) { + using stan::test::relative_tolerance; auto f = [](const auto& t, const auto& a, const auto& b) { return stan::math::scale_matrix_exp_multiply(t, a, b); }; @@ -27,8 +28,8 @@ TEST(mathMixMatFun, scaleMatrixExpMultiply) { stan::test::expect_ad(f, t, a11, b14); stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; + tols.hessian_hessian_ = relative_tolerance(1e-4, 1e-2); + tols.hessian_fvar_hessian_ = relative_tolerance(1e-4, 1e-2); // 3 x 1 Eigen::MatrixXd a33(3, 3); diff --git a/test/unit/math/mix/fun/tcrossprod_test.cpp b/test/unit/math/mix/fun/tcrossprod_test.cpp index c060bbcc959..aad9b4ab926 100644 --- a/test/unit/math/mix/fun/tcrossprod_test.cpp +++ b/test/unit/math/mix/fun/tcrossprod_test.cpp @@ -1,6 +1,7 @@ #include TEST(MathMixMatFun, tcrossprod) { + using stan::test::relative_tolerance; auto f = [](const auto& y) { return stan::math::tcrossprod(y); }; Eigen::MatrixXd a00(0, 0); @@ -33,8 +34,8 @@ TEST(MathMixMatFun, tcrossprod) { stan::test::expect_ad(f, a32); stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; + tols.hessian_hessian_ = relative_tolerance(2e-3, 2e-4); + tols.hessian_fvar_hessian_ = relative_tolerance(2e-3, 2e-4); Eigen::MatrixXd a33(3, 3); a33 << 1, 0, 0, 2, 3, 0, 4, 5, 6; diff --git a/test/unit/math/mix/fun/trace_gen_quad_form_test.cpp b/test/unit/math/mix/fun/trace_gen_quad_form_test.cpp index e233c5d5eff..e18fe59f360 100644 --- a/test/unit/math/mix/fun/trace_gen_quad_form_test.cpp +++ b/test/unit/math/mix/fun/trace_gen_quad_form_test.cpp @@ -1,6 +1,7 @@ #include TEST(mathMixMatFun, traceGenQuadForm) { + using stan::test::relative_tolerance; auto f = [](const auto& x1, const auto& x2, const auto& x3) { return stan::math::trace_gen_quad_form(x1, x2, x3); }; @@ -19,8 +20,8 @@ TEST(mathMixMatFun, traceGenQuadForm) { stan::test::expect_ad(f, c11, a11, b11); stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-1; - tols.hessian_fvar_hessian_ = 1e-1; + tols.hessian_hessian_ = relative_tolerance(1e-1, 1e-1); + tols.hessian_fvar_hessian_ = relative_tolerance(1e-1, 1e-1); Eigen::MatrixXd a(4, 4); a << 2, 3, 4, 5, 6, 10, 2, 2, 7, 2, 7, 1, 8, 2, 1, 112; diff --git a/test/unit/math/mix/fun/trigamma_test.cpp b/test/unit/math/mix/fun/trigamma_test.cpp index 012eb1e13ca..bc2da319f0c 100644 --- a/test/unit/math/mix/fun/trigamma_test.cpp +++ b/test/unit/math/mix/fun/trigamma_test.cpp @@ -4,29 +4,30 @@ TEST(mathMixMatFun, trigamma) { using stan::math::trigamma; using stan::test::ad_tolerances; using stan::test::expect_unary_vectorized; + using stan::test::relative_tolerance; auto f = [](const auto& x1) { return trigamma(x1); }; - // reduce second and third order tests one order of magnitude + // reduce the tol_min for second and third order tests one order of magnitude stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; - tols.grad_hessian_hessian_ = 1e-2; - tols.grad_hessian_grad_hessian_ = 1e-1; + tols.hessian_hessian_ = relative_tolerance(1e-4, 1e-2); + tols.hessian_fvar_hessian_ = relative_tolerance(1e-4, 1e-2); + tols.grad_hessian_hessian_ = relative_tolerance(1e-3, 1e-2); + tols.grad_hessian_grad_hessian_ = relative_tolerance(1e-2, 1e-1); expect_unary_vectorized(tols, f, -103.52, -0.9, -0.5, 0, 0.5, 1.3, 5.1, 19.2); - // reduce first deriv tests one order, second derivs four orders, + // reduce tol_min for first deriv tests one order, second derivs four orders, // and third derivs three orders stan::test::ad_tolerances tols2; - tols2.gradient_grad_ = 1e-3; - tols2.gradient_fvar_grad_ = 1e-3; - tols2.hessian_grad_ = 1e-3; - tols2.hessian_fvar_grad_ = 1e-3; - tols2.hessian_hessian_ = 1e2; - tols2.hessian_fvar_hessian_ = 1e2; - tols2.grad_hessian_hessian_ = 1e1; - tols2.grad_hessian_grad_hessian_ = 1e1; + tols2.gradient_grad_ = relative_tolerance(1e-4, 1e-3); + tols2.gradient_fvar_grad_ = relative_tolerance(1e-4, 1e-3); + tols2.hessian_grad_ = relative_tolerance(1e-4, 1e-3); + tols2.hessian_fvar_grad_ = relative_tolerance(1e-4, 1e-3); + tols2.hessian_hessian_ = relative_tolerance(1e2, 1); + tols2.hessian_fvar_hessian_ = relative_tolerance(1e2, 1); + tols2.grad_hessian_hessian_ = relative_tolerance(1, 1); + tols2.grad_hessian_grad_hessian_ = relative_tolerance(1, 1); expect_unary_vectorized(tols2, f, -20, 1, 5); } diff --git a/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp b/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp index 866ac14766b..e288331aae0 100644 --- a/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp +++ b/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp @@ -1,6 +1,6 @@ #include +#include #include -#include #include template @@ -15,8 +15,8 @@ void expect_match_autodiff(const F& f, Eigen::VectorXd x) { EXPECT_FLOAT_EQ(fx, fx_fd); EXPECT_EQ(grad_fx.size(), grad_fx_fd.size()); - for (size_t i = 0; i < grad_fx.size(); ++i) - expect_near_relative(grad_fx(i), grad_fx_fd(i)); + stan::test::expect_near_rel("expect_match_autodiff", grad_fx, grad_fx_fd, + stan::test::relative_tolerance(1e-7, 1e-9)); } TEST(MathMixMatFunctor, FiniteDiffGradientAuto) { diff --git a/test/unit/math/relative_tolerance.hpp b/test/unit/math/relative_tolerance.hpp new file mode 100644 index 00000000000..d6ea28b2293 --- /dev/null +++ b/test/unit/math/relative_tolerance.hpp @@ -0,0 +1,95 @@ +#ifndef TEST_UNIT_MATH_RELATIVE_TOLERANCE_HPP +#define TEST_UNIT_MATH_RELATIVE_TOLERANCE_HPP + +#include +#include + +namespace stan { +namespace test { + +/** + * Class holding information about relative tolerance and the minimal absolute + * tolerance that should be tested against. The final tolerance is computed as + * max(tol * fabs(x), tol_min) + * Where x is either an exact value to be tested against or average of + * two inexact values to be compared. + */ +class relative_tolerance { + public: + /** + * Construct with default tolerances + */ + relative_tolerance() : tol_(1e-8), tol_min_(1e-14) {} + + /** + * Construct with default tol_min (max(tol_ * tol_, 1e-14)) + * @param tol_ the relative tolerance + */ + relative_tolerance(const double tol) // NOLINT + : tol_(tol), tol_min_(std::max(tol * tol, 1e-14)) {} + + /** + * Construct fully specified + * @param[in] tol_ the relative tolerance + * @param[in] tol_min_ the minimum absolute tolerance + */ + relative_tolerance(const double tol, const double tol_min) + : tol_(tol), tol_min_(tol_min) {} + + double tol() const { return tol_; } + double tol_min() const { return tol_min_; } + + relative_tolerance change_tol(double tol) const { + return relative_tolerance(tol, tol_min_); + } + + relative_tolerance change_tol_min(double tol_min) const { + return relative_tolerance(tol_, tol_min); + } + + relative_tolerance operator*(double a) const { + return relative_tolerance(a * tol_, a * tol_min_); + } + + /** + * Computes tolerance around an exact target value. + * + * @tparam T1 the type of argument, must be scalar + * @param[in] x the target value + * @return tolerance + */ + template ...> + double exact(const T1& x) const { + using stan::math::fabs; + return std::max(tol_ * fabs(x), tol_min_); + } + + /** + * Computes tolerance for comparing two inexact values. + * + * @tparam T1 the type of first argument, must be scalar + * @tparam T2 the type of second argument, must be scalar + * @param[in] x first value that will be compared + * @param[in] y second value that will be compared + * @return tolerance + */ + template ...> + double inexact(const T1& x, const T2& y) const { + using stan::math::fabs; + return std::max(tol_ * 0.5 * (fabs(x) + fabs(y)), tol_min_); + } + + private: + /** + * The relative tolerance + */ + double tol_; + /** + * The minimal absolute tolerance + */ + double tol_min_; +}; + +} // namespace test +} // namespace stan +#endif diff --git a/test/unit/math/relative_tolerance_test.cpp b/test/unit/math/relative_tolerance_test.cpp new file mode 100644 index 00000000000..b6f20375172 --- /dev/null +++ b/test/unit/math/relative_tolerance_test.cpp @@ -0,0 +1,44 @@ +#include +#include +#include +#include + +TEST(testUnitMath, RelativeTolerance) { + using stan::test::relative_tolerance; + + // Note that default tol_min is 1e-14 in this case + relative_tolerance t_8(1e-8); + EXPECT_FLOAT_EQ(t_8.exact(1e3), 1e-5); + EXPECT_FLOAT_EQ(t_8.exact(1), 1e-8); + EXPECT_FLOAT_EQ(t_8.exact(1e-6), 1e-14); + EXPECT_FLOAT_EQ(t_8.exact(1e-14), 1e-14); + EXPECT_FLOAT_EQ(t_8.exact(0), 1e-14); + + EXPECT_FLOAT_EQ(t_8.inexact(-1e3, 3e3), 2e-5); + EXPECT_FLOAT_EQ(t_8.inexact(1, 1000), 5.005e-6); + EXPECT_FLOAT_EQ(t_8.inexact(1e-5, -8e-6), 9e-14); + EXPECT_FLOAT_EQ(t_8.inexact(0, 2e3), 1e-5); + EXPECT_FLOAT_EQ(t_8.inexact(0, 1e-10), 1e-14); + + // Default tol_min is 1e-8 + relative_tolerance t_4(1e-4); + EXPECT_FLOAT_EQ(t_4.exact(1e3), 1e-1); + EXPECT_FLOAT_EQ(t_4.exact(1), 1e-4); + EXPECT_FLOAT_EQ(t_4.exact(-1e-3), 1e-7); + EXPECT_FLOAT_EQ(t_4.exact(1e-5), 1e-8); + EXPECT_FLOAT_EQ(t_4.exact(0), 1e-8); + + EXPECT_FLOAT_EQ(t_4.inexact(-1e3, 3e3), 2e-1); + EXPECT_FLOAT_EQ(t_4.inexact(2e-5, 0), 1e-8); + EXPECT_FLOAT_EQ(t_4.inexact(0, 1), 5e-5); + EXPECT_FLOAT_EQ(t_4.inexact(0, -1e-10), 1e-8); + + relative_tolerance t_16_10(1e-16, 1e-10); + EXPECT_FLOAT_EQ(t_16_10.exact(-30), 1e-10); + EXPECT_FLOAT_EQ(t_16_10.exact(-500), 1e-10); + EXPECT_FLOAT_EQ(t_16_10.exact(1e-3), 1e-10); + EXPECT_FLOAT_EQ(t_16_10.exact(-1e-5), 1e-10); + + EXPECT_FLOAT_EQ(t_16_10.inexact(-5e6, 1e-14), 2.5e-10); + EXPECT_FLOAT_EQ(t_16_10.inexact(1e-315, 1e-2), 1e-10); +} diff --git a/test/unit/math/rev/expect_identity.hpp b/test/unit/math/rev/expect_identity.hpp new file mode 100644 index 00000000000..b223f37e3aa --- /dev/null +++ b/test/unit/math/rev/expect_identity.hpp @@ -0,0 +1,54 @@ +#ifndef TEST_UNIT_MATH_REV_EXPECT_IDENTITY_HPP +#define TEST_UNIT_MATH_REV_EXPECT_IDENTITY_HPP + +#include +#include +#include +#include + +namespace stan { +namespace test { + +template +void expect_identity(const std::string& msg, const F& f, const G& g, Ts... xs) { + using stan::math::gradient; + using stan::math::pack_params; + + Eigen::Matrix x_vec(xs...); + + auto packed_f = pack_params(f); + auto packed_g = pack_params(g); + + double value_f; + Eigen::Matrix grad_f; + + double value_g; + Eigen::Matrix grad_g; + + gradient(packed_f, x_vec, value_f, grad_f); + gradient(packed_g, x_vec, value_g, grad_g); + + std::stringstream x_strstr; + x_strstr << std::setprecision(24) << " for params: "; + for (size_t i = 0; i < x_vec.size(); ++i) { + if (i > 0) { + x_strstr << ", "; + } + x_strstr << x_vec(i); + } + std::string x_str = x_strstr.str(); + + std::stringstream msg_val; + msg_val << msg << ": value" << x_str; + expect_near_rel(msg_val.str(), value_f, value_g); + for (size_t i = 0; i < grad_f.size(); ++i) { + std::stringstream msg_grad; + msg_grad << msg << ": grad(" << i << ")" << x_str; + expect_near_rel(msg_grad.str(), grad_f(i), grad_g(i)); + } +} + +} // namespace test +} // namespace stan + +#endif diff --git a/test/unit/math/rev/fun/log_sum_exp_test.cpp b/test/unit/math/rev/fun/log_sum_exp_test.cpp new file mode 100644 index 00000000000..65469bcaeb6 --- /dev/null +++ b/test/unit/math/rev/fun/log_sum_exp_test.cpp @@ -0,0 +1,44 @@ +#include +#include +#include +#include +#include + +TEST(MathFunctions, log_sum_exp_identities_rev) { + using stan::math::log_sum_exp; + using stan::math::var; + using stan::test::expect_identity; + + std::vector values_to_check + = {1e-8, 3e-5, 8e-1, 1, 8.6, 4.3e3, 9.6e6, 9e9, 1.8e15, 4.3e30}; + + auto lh_duplicate = [](const var& x) { return log_sum_exp(x, x); }; + auto rh_duplicate = [](const var& x) { return stan::math::LOG_TWO + x; }; + + auto lh_plus_one = [](const var& x) { return log_sum_exp(x, x + 1); }; + auto rh_plus_one = [](const var& x) { return x + log(1 + stan::math::e()); }; + + for (double x : values_to_check) { + expect_identity("Duplicate argument", lh_duplicate, rh_duplicate, x); + expect_identity("Plus one", lh_plus_one, rh_plus_one, x); + } + + auto lh_multiply = [](const var& x, const var& y) { + return log_sum_exp(1e5 + x, 1e5 + y); + }; + + auto rh_multiply + = [](const var& x, const var& y) { return 1e5 + log_sum_exp(x, y); }; + auto lh_orig = [](const var& x, const var& y) { return log_sum_exp(x, y); }; + + auto rh_distribute + = [](const var& x, const var& y) { return x + log_sum_exp(0, y - x); }; + + for (double x : values_to_check) { + for (double y : values_to_check) { + expect_identity("Multiply", lh_multiply, rh_multiply, x, y); + expect_identity("Distribute", lh_orig, rh_distribute, x, y); + } + } + // log(exp(1e5)) +} diff --git a/test/unit/math/rev/fun/util.hpp b/test/unit/math/rev/fun/util.hpp index 62cc02e4569..2cdc459b6ba 100644 --- a/test/unit/math/rev/fun/util.hpp +++ b/test/unit/math/rev/fun/util.hpp @@ -88,19 +88,6 @@ VEC cgradvec(AVAR f, AVEC x) { return g; } -double relative_diff(double u, double v) { - return 2 * (u - v) / (fabs(u) + fabs(v)); -} - -void expect_near_relative(double u, double v) { - if (u == v) - SUCCEED(); - else if (u == 0 || v == 0) - EXPECT_NEAR(0, (u - v), 1e-7); - else - EXPECT_NEAR(0, relative_diff(u, v), 1e-7); -} - using size_type = stan::math::index_type_t>; // Returns a matrix with the contents of a