From 47cdde2f9b8628bb3f76ac15c4eb1c088f303e75 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Oct 2023 09:31:34 +0300 Subject: [PATCH 01/36] Initial natural scale impl --- stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/grad_pFq_2.hpp | 88 +++++ test/unit/math/rev/fun/grad_pFq_2_test.cpp | 430 +++++++++++++++++++++ 3 files changed, 519 insertions(+) create mode 100644 stan/math/prim/fun/grad_pFq_2.hpp create mode 100644 test/unit/math/rev/fun/grad_pFq_2_test.cpp diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 11ba14123bd..364a56c6e85 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -130,6 +130,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/grad_pFq_2.hpp b/stan/math/prim/fun/grad_pFq_2.hpp new file mode 100644 index 00000000000..442c5f34170 --- /dev/null +++ b/stan/math/prim/fun/grad_pFq_2.hpp @@ -0,0 +1,88 @@ +#ifndef STAN_MATH_PRIM_FUN_GRAD_PFQ_2_HPP +#define STAN_MATH_PRIM_FUN_GRAD_PFQ_2_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +template +auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, + double precision = 1e-16, + int max_steps = 1e8) { + using T_partials = partials_return_t; + using Ta_partials = promote_scalar_t>; + using Tb_partials = promote_scalar_t>; + using Tz_partials = promote_scalar_t>; + using T_ArrayPartials = promote_scalar_t>; + using T_return = std::tuple; + + T_ArrayPartials a_val = as_value_column_array_or_scalar(a); + T_ArrayPartials b_val = as_value_column_array_or_scalar(b); + const auto& z_val = value_of(z); + T_ArrayPartials a_k = a_val; + T_ArrayPartials b_k = b_val; + T_ArrayPartials phammer_a_k = T_ArrayPartials::Ones(a.size()); + T_ArrayPartials phammer_b_k = T_ArrayPartials::Ones(b.size()); + T_ArrayPartials digamma_a = inv(a_k); + T_ArrayPartials digamma_b = inv(b_k); + T_ArrayPartials digamma_a_k = digamma_a; + T_ArrayPartials digamma_b_k = digamma_b; + T_ArrayPartials a_infsum = T_ArrayPartials::Zero(a.size()); + T_ArrayPartials b_infsum = T_ArrayPartials::Zero(b.size()); + + int k = 0; + double curr_log_prec = 1.0; + double curr_prec = 1.0; + double factorial_k = 1; + while (k <= max_steps && curr_prec > precision) { + T_partials num = prod(phammer_a_k) * pow(z_val, k); + T_partials den = factorial_k * prod(phammer_b_k); + T_ArrayPartials a_grad = (digamma_a_k * num) / den; + T_ArrayPartials b_grad = (digamma_b_k * num) / den; + +/* std::cout << "k: " << k << "\n" + << "num: " << num << "\n" + << "den: " << den << "\n" + << "phammer_a_k: " << phammer_a_k.matrix().transpose() << "\n" + << "phammer_b_k: " << phammer_b_k.matrix().transpose() << "\n" + << "digamma_a_k: " << digamma_a_k.matrix().transpose() << "\n" + << "digamma_b_k: " << digamma_b_k.matrix().transpose() << "\n" + << "a: " << a_grad.matrix().transpose() << "\n" + << "b: " << b_grad.matrix().transpose() << "\n" + << std::endl; */ + + curr_prec = max(a_grad.abs().maxCoeff(), b_grad.abs().maxCoeff()); + a_infsum += a_grad; + b_infsum += b_grad; + + + phammer_a_k *= a_k; + phammer_b_k *= b_k; + digamma_a_k += select(a_k == 0.0, 0.0, inv(a_k)); + digamma_b_k += select(b_k == 0.0, 0.0, inv(b_k)); + a_k += 1.0; + b_k += 1.0; + k += 1; + factorial_k *= k; + } + + T_partials pfq_val = hypergeometric_pFq(a_val, b_val, z_val); + T_partials pfq_p1_val = hypergeometric_pFq(a_val + 1, b_val + 1, z_val); + + T_return ret_tuple; + std::get<0>(ret_tuple) = a_infsum - digamma_a * pfq_val; + std::get<1>(ret_tuple) = digamma_b * pfq_val - b_infsum; + std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; + return ret_tuple; +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/rev/fun/grad_pFq_2_test.cpp b/test/unit/math/rev/fun/grad_pFq_2_test.cpp new file mode 100644 index 00000000000..3d1951840f8 --- /dev/null +++ b/test/unit/math/rev/fun/grad_pFq_2_test.cpp @@ -0,0 +1,430 @@ +#include +#include +#include + +TEST(RevMath, grad_2F2) { + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_d; + using stan::math::vector_v; + + vector_v a_v(2); + vector_d a_d(2); + a_v << 4, 2; + a_d << 4, 2; + + vector_v b_v(2); + vector_d b_d(2); + b_v << 6, 3; + b_d << 6, 3; + + var z_v = 4; + double z_d = 4; + + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); + + EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); + + EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); + + grad_tuple = grad_pFq_2(a_v, b_d, z_d); + + EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); + + grad_tuple = grad_pFq_2(a_d, b_v, z_d); + + EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); + + grad_tuple = grad_pFq_2(a_d, b_d, z_v); + + EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); + + grad_tuple = grad_pFq_2(a_v, b_v, z_d); + + EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); + + grad_tuple = grad_pFq_2(a_v, b_d, z_v); + + EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); + + grad_tuple = grad_pFq_2(a_d, b_v, z_v); + + EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); +} +TEST(RevMath, grad_2F3) { + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << 2, 3; + vector_v b_v(3); + b_v << 2, 4, 5; + var z_v = 1; + + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(0.08377717301140296, std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(0.05615450733193106, std::get<0>(grad_tuple)[1]); + + EXPECT_FLOAT_EQ(-0.08377717301140296, std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(-0.04225296806591615, std::get<1>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(-0.03387575989873739, std::get<1>(grad_tuple)[2]); + + EXPECT_FLOAT_EQ(0.1712340452215524, std::get<2>(grad_tuple)); +} + +TEST(RevMath, grad_4F3) { + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(4); + a_v << 1, 2, 3, 4; + vector_v b_v(3); + b_v << 5, 6, 7; + var z_v = 1; + + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(0.1587098625610631, std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(0.08249338029396848, std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(0.05611368752226367, std::get<0>(grad_tuple)[2]); + EXPECT_FLOAT_EQ(0.04261209968272329, std::get<0>(grad_tuple)[3]); + + EXPECT_FLOAT_EQ(-0.03438035893346993, std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(-0.02882791253333995, std::get<1>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(-0.02482622713079612, std::get<1>(grad_tuple)[2]); + + EXPECT_FLOAT_EQ(0.1800529055890911, std::get<2>(grad_tuple)); +} + +TEST(RevMath, grad_2F1_derivs_match) { + using stan::math::grad_2F1; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << 1, 1; + vector_v b_v(1); + b_v << 1; + var z_v = 0.6; + + auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); +} + +TEST(RevMath, grad2F1_2) { + using stan::math::grad_2F1; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << 1, 31; + vector_v b_v(1); + b_v << 41; + var z_v = 1; + + auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); +} + + +TEST(RevMath, grad2F1_3) { + using stan::math::grad_2F1; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << 1.0, -2.1; + vector_v b_v(1); + b_v << 41.0; + var z_v = 1.0; + + auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); +} + +TEST(RevMath, grad2F1_6) { + using stan::math::grad_2F1; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << 1, -0.5; + vector_v b_v(1); + b_v << 10.6; + var z_v = 0.3; + + auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); +} + +TEST(RevMath, grad2F1_7) { + using stan::math::grad_2F1; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << 1, -0.5; + vector_v b_v(1); + b_v << 10; + var z_v = 0.3; + + auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); +} + +TEST(RevMath, grad2F1_8) { + using stan::math::grad_2F1; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << -0.5, -4.5; + vector_v b_v(1); + b_v << 11; + var z_v = 0.3; + + auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); +} + +TEST(RevMath, grad2F1_9) { + using stan::math::grad_2F1; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << -0.5, -4.5; + vector_v b_v(1); + b_v << -3.2; + var z_v = 0.9; + + auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); +} + +TEST(RevMath, grad2F1_10) { + using stan::math::grad_2F1; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << 2, 1; + vector_v b_v(1); + b_v << 2; + var z_v = 0.4; + + auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); +} + +/* +TEST(RevMath, grad2F1_11) { + using stan::math::grad_2F1; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << 3.70975, 1; + vector_v b_v(1); + b_v << 2.70975; + var z_v = 0.999696; + + double std::get<0>(grad_2F1_tuple); + double std::get<2>(grad_2F1_tuple); + + grad_2F1(std::get<0>(grad_2F1_tuple), std::get<2>(grad_2F1_tuple), a_v[0], +a_v[1], b_v[0], z_v); auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); +}*/ +/* +TEST(RevMath, F32_converges_by_z) { + using stan::math::grad_F32; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(3); + a_v << 1.0, 1.0, 1.0; + vector_v b_v(2); + b_v << 1.0, 1.0; + var z_v = 0.6; + + double g_calc[6]; + + grad_F32(g_calc, 1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 1e-10); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(g_calc[2], std::get<0>(grad_tuple)[2]); + EXPECT_FLOAT_EQ(g_calc[3], std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(g_calc[4], std::get<1>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); +} +*/ + +TEST(RevMath, grad_F32_double_sign_flip_1) { + using stan::math::grad_F32; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(3); + a_v << 1.0, -0.5, -2.5; + vector_v b_v(2); + b_v << 10.0, 1.0; + var z_v = 0.3; + + double g_calc[6]; + + grad_F32(g_calc, 1.0, -.5, -2.5, 10.0, 1.0, 0.3, 1e-10); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(g_calc[2], std::get<0>(grad_tuple)[2]); + EXPECT_FLOAT_EQ(g_calc[3], std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(g_calc[4], std::get<1>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); +} + +TEST(RevMath, grad_F32_double_sign_flip_2) { + using stan::math::grad_F32; + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(3); + a_v << 1.0, -0.5, -4.5; + vector_v b_v(2); + b_v << 10.0, 1.0; + var z_v = 0.3; + + double g_calc[6]; + + grad_F32(g_calc, 1.0, -.5, -4.5, 10.0, 1.0, 0.3, 1e-10); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(g_calc[2], std::get<0>(grad_tuple)[2]); + EXPECT_FLOAT_EQ(g_calc[3], std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(g_calc[4], std::get<1>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); +} + +TEST(RevMath, grad_2F1_negative_z) { + using stan::math::grad_pFq_2; + using stan::math::var; + using stan::math::vector_v; + + vector_v a_v(2); + a_v << 3.70975, 1.0; + vector_v b_v(1); + b_v << 2.70975; + var z_v = -0.2; + + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(-0.0488658806159776, std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(-0.193844936204681, std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(0.0677809985598383, std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(0.865295247272367, std::get<2>(grad_tuple)); +} + + +TEST(RevMath, grad_3F2_cross_zero) { + using stan::math::grad_F32; + using stan::math::grad_pFq_2; + using stan::math::hypergeometric_3F2; + using stan::math::var; + using stan::math::vector_d; + using stan::math::vector_v; + + vector_v a_v(3); + a_v << 1, 1, -1; + + vector_v b_v(2); + b_v << 2, 2; + + var z_v = 0.292893; + + double g_calc[6]; + grad_F32(g_calc, 1.0, 1.0, -1.0, 2.0, 2.0, 0.292893); + + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + + EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(g_calc[2], std::get<0>(grad_tuple)[2]); + EXPECT_FLOAT_EQ(g_calc[3], std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(g_calc[4], std::get<1>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); +} From f091e0c61af0f3a90b926f9340429f952d8eb205 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Oct 2023 10:20:14 +0300 Subject: [PATCH 02/36] Absolute impl on natural scale --- stan/math/prim/fun/grad_pFq_2.hpp | 79 +++++++++++++++++----- test/unit/math/rev/fun/grad_pFq_2_test.cpp | 4 +- 2 files changed, 63 insertions(+), 20 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq_2.hpp b/stan/math/prim/fun/grad_pFq_2.hpp index 442c5f34170..9e7c8073bd2 100644 --- a/stan/math/prim/fun/grad_pFq_2.hpp +++ b/stan/math/prim/fun/grad_pFq_2.hpp @@ -11,6 +11,13 @@ namespace stan { namespace math { +namespace internal { + template + auto bin_sign(const T& x) { + decltype(auto) sign1 = sign(value_of_rec(x)); + return select(sign1 == 0, 1, sign1); + } +} template auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, @@ -22,18 +29,38 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, using Tz_partials = promote_scalar_t>; using T_ArrayPartials = promote_scalar_t>; using T_return = std::tuple; + using IntArray = Eigen::Array; T_ArrayPartials a_val = as_value_column_array_or_scalar(a); T_ArrayPartials b_val = as_value_column_array_or_scalar(b); - const auto& z_val = value_of(z); + Tz_partials z_val = value_of(z); + Tz_partials z_k = abs(z_val); + + IntArray a_signs = internal::bin_sign(a_val); + IntArray b_signs = internal::bin_sign(b_val); + int z_sign = internal::bin_sign(z_val); + + IntArray a_k_signs = a_signs; + IntArray b_k_signs = b_signs; + int z_k_sign = 1; + T_ArrayPartials a_k = a_val; T_ArrayPartials b_k = b_val; + Tz_partials z_pow_k = 1; + T_ArrayPartials phammer_a_k = T_ArrayPartials::Ones(a.size()); T_ArrayPartials phammer_b_k = T_ArrayPartials::Ones(b.size()); + IntArray phammer_a_k_signs = IntArray::Ones(a.size()); + IntArray phammer_b_k_signs = IntArray::Ones(b.size()); + T_ArrayPartials digamma_a = inv(a_k); T_ArrayPartials digamma_b = inv(b_k); T_ArrayPartials digamma_a_k = digamma_a; T_ArrayPartials digamma_b_k = digamma_b; + T_ArrayPartials abs_digamma_a_k = abs(digamma_a); + T_ArrayPartials abs_digamma_b_k = abs(digamma_b); + IntArray digamma_a_k_signs = a_k_signs; + IntArray digamma_b_k_signs = b_k_signs; T_ArrayPartials a_infsum = T_ArrayPartials::Zero(a.size()); T_ArrayPartials b_infsum = T_ArrayPartials::Zero(b.size()); @@ -42,33 +69,36 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, double curr_prec = 1.0; double factorial_k = 1; while (k <= max_steps && curr_prec > precision) { - T_partials num = prod(phammer_a_k) * pow(z_val, k); - T_partials den = factorial_k * prod(phammer_b_k); - T_ArrayPartials a_grad = (digamma_a_k * num) / den; - T_ArrayPartials b_grad = (digamma_b_k * num) / den; - -/* std::cout << "k: " << k << "\n" - << "num: " << num << "\n" - << "den: " << den << "\n" - << "phammer_a_k: " << phammer_a_k.matrix().transpose() << "\n" - << "phammer_b_k: " << phammer_b_k.matrix().transpose() << "\n" - << "digamma_a_k: " << digamma_a_k.matrix().transpose() << "\n" - << "digamma_b_k: " << digamma_b_k.matrix().transpose() << "\n" - << "a: " << a_grad.matrix().transpose() << "\n" - << "b: " << b_grad.matrix().transpose() << "\n" - << std::endl; */ + int base_sign = z_k_sign * phammer_a_k_signs.prod() * phammer_b_k_signs.prod(); + T_partials base = base_sign * (prod(phammer_a_k) * z_pow_k) + / (factorial_k * prod(phammer_b_k)); + T_ArrayPartials a_grad = abs_digamma_a_k * base * digamma_a_k_signs; + T_ArrayPartials b_grad = abs_digamma_b_k * base * digamma_b_k_signs; curr_prec = max(a_grad.abs().maxCoeff(), b_grad.abs().maxCoeff()); a_infsum += a_grad; b_infsum += b_grad; + phammer_a_k *= abs(a_k); + phammer_b_k *= abs(b_k); + phammer_a_k_signs *= a_k_signs; + phammer_b_k_signs *= b_k_signs; + + z_pow_k *= z_k; + z_k_sign *= z_sign; - phammer_a_k *= a_k; - phammer_b_k *= b_k; digamma_a_k += select(a_k == 0.0, 0.0, inv(a_k)); digamma_b_k += select(b_k == 0.0, 0.0, inv(b_k)); + abs_digamma_a_k = abs(digamma_a_k); + abs_digamma_b_k = abs(digamma_b_k); + digamma_a_k_signs = internal::bin_sign(digamma_a_k); + digamma_b_k_signs = internal::bin_sign(digamma_b_k); + a_k += 1.0; b_k += 1.0; + a_k_signs = internal::bin_sign(a_k); + b_k_signs = internal::bin_sign(b_k); + k += 1; factorial_k *= k; } @@ -86,3 +116,16 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, } // namespace math } // namespace stan #endif + + + +/* std::cout << "k: " << k << "\n" + << "num: " << num << "\n" + << "den: " << den << "\n" + << "phammer_a_k: " << phammer_a_k.matrix().transpose() << "\n" + << "phammer_b_k: " << phammer_b_k.matrix().transpose() << "\n" + << "digamma_a_k: " << digamma_a_k.matrix().transpose() << "\n" + << "digamma_b_k: " << digamma_b_k.matrix().transpose() << "\n" + << "a: " << a_grad.matrix().transpose() << "\n" + << "b: " << b_grad.matrix().transpose() << "\n" + << std::endl; */ diff --git a/test/unit/math/rev/fun/grad_pFq_2_test.cpp b/test/unit/math/rev/fun/grad_pFq_2_test.cpp index 3d1951840f8..423c952c5aa 100644 --- a/test/unit/math/rev/fun/grad_pFq_2_test.cpp +++ b/test/unit/math/rev/fun/grad_pFq_2_test.cpp @@ -132,7 +132,7 @@ TEST(RevMath, grad_2F1_derivs_match) { EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } - +/* TEST(RevMath, grad2F1_2) { using stan::math::grad_2F1; using stan::math::grad_pFq_2; @@ -153,7 +153,7 @@ TEST(RevMath, grad2F1_2) { EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } - + */ TEST(RevMath, grad2F1_3) { using stan::math::grad_2F1; From 7d58254f7caf47e4b95a1fe916c20eda7e9f7de6 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Oct 2023 11:26:20 +0300 Subject: [PATCH 03/36] Log scale initial impl --- stan/math/prim/fun/grad_pFq_2.hpp | 42 +++++++++++----------- test/unit/math/rev/fun/grad_pFq_2_test.cpp | 17 ++++----- 2 files changed, 28 insertions(+), 31 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq_2.hpp b/stan/math/prim/fun/grad_pFq_2.hpp index 9e7c8073bd2..8fc7337f955 100644 --- a/stan/math/prim/fun/grad_pFq_2.hpp +++ b/stan/math/prim/fun/grad_pFq_2.hpp @@ -34,7 +34,6 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, T_ArrayPartials a_val = as_value_column_array_or_scalar(a); T_ArrayPartials b_val = as_value_column_array_or_scalar(b); Tz_partials z_val = value_of(z); - Tz_partials z_k = abs(z_val); IntArray a_signs = internal::bin_sign(a_val); IntArray b_signs = internal::bin_sign(b_val); @@ -46,10 +45,11 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, T_ArrayPartials a_k = a_val; T_ArrayPartials b_k = b_val; - Tz_partials z_pow_k = 1; + Tz_partials log_z_pow_k = 0.0; + Tz_partials log_z_k = log(abs(z_val)); - T_ArrayPartials phammer_a_k = T_ArrayPartials::Ones(a.size()); - T_ArrayPartials phammer_b_k = T_ArrayPartials::Ones(b.size()); + T_ArrayPartials log_phammer_a_k = T_ArrayPartials::Zero(a.size()); + T_ArrayPartials log_phammer_b_k = T_ArrayPartials::Zero(b.size()); IntArray phammer_a_k_signs = IntArray::Ones(a.size()); IntArray phammer_b_k_signs = IntArray::Ones(b.size()); @@ -57,8 +57,8 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, T_ArrayPartials digamma_b = inv(b_k); T_ArrayPartials digamma_a_k = digamma_a; T_ArrayPartials digamma_b_k = digamma_b; - T_ArrayPartials abs_digamma_a_k = abs(digamma_a); - T_ArrayPartials abs_digamma_b_k = abs(digamma_b); + T_ArrayPartials log_digamma_a_k = log(abs(digamma_a)); + T_ArrayPartials log_digamma_b_k = log(abs(digamma_b)); IntArray digamma_a_k_signs = a_k_signs; IntArray digamma_b_k_signs = b_k_signs; T_ArrayPartials a_infsum = T_ArrayPartials::Zero(a.size()); @@ -67,30 +67,30 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, int k = 0; double curr_log_prec = 1.0; double curr_prec = 1.0; - double factorial_k = 1; - while (k <= max_steps && curr_prec > precision) { + double log_factorial_k = 0; + while (k <= max_steps && curr_log_prec > log(precision)) { int base_sign = z_k_sign * phammer_a_k_signs.prod() * phammer_b_k_signs.prod(); - T_partials base = base_sign * (prod(phammer_a_k) * z_pow_k) - / (factorial_k * prod(phammer_b_k)); - T_ArrayPartials a_grad = abs_digamma_a_k * base * digamma_a_k_signs; - T_ArrayPartials b_grad = abs_digamma_b_k * base * digamma_b_k_signs; + T_partials log_base = (sum(log_phammer_a_k) + log_z_pow_k) + - (log_factorial_k + sum(log_phammer_b_k)); + T_ArrayPartials a_grad = log_digamma_a_k + log_base; + T_ArrayPartials b_grad = log_digamma_b_k + log_base; - curr_prec = max(a_grad.abs().maxCoeff(), b_grad.abs().maxCoeff()); - a_infsum += a_grad; - b_infsum += b_grad; + curr_log_prec = max(a_grad.maxCoeff(), b_grad.maxCoeff()); + a_infsum += exp(a_grad) * base_sign * digamma_a_k_signs; + b_infsum += exp(b_grad) * base_sign * digamma_b_k_signs; - phammer_a_k *= abs(a_k); - phammer_b_k *= abs(b_k); + log_phammer_a_k += log(abs(a_k)); + log_phammer_b_k += log(abs(b_k)); phammer_a_k_signs *= a_k_signs; phammer_b_k_signs *= b_k_signs; - z_pow_k *= z_k; + log_z_pow_k += log_z_k; z_k_sign *= z_sign; digamma_a_k += select(a_k == 0.0, 0.0, inv(a_k)); digamma_b_k += select(b_k == 0.0, 0.0, inv(b_k)); - abs_digamma_a_k = abs(digamma_a_k); - abs_digamma_b_k = abs(digamma_b_k); + log_digamma_a_k = log(abs(digamma_a_k)); + log_digamma_b_k = log(abs(digamma_b_k)); digamma_a_k_signs = internal::bin_sign(digamma_a_k); digamma_b_k_signs = internal::bin_sign(digamma_b_k); @@ -100,7 +100,7 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, b_k_signs = internal::bin_sign(b_k); k += 1; - factorial_k *= k; + log_factorial_k += log(k); } T_partials pfq_val = hypergeometric_pFq(a_val, b_val, z_val); diff --git a/test/unit/math/rev/fun/grad_pFq_2_test.cpp b/test/unit/math/rev/fun/grad_pFq_2_test.cpp index 423c952c5aa..2db6bb8ee47 100644 --- a/test/unit/math/rev/fun/grad_pFq_2_test.cpp +++ b/test/unit/math/rev/fun/grad_pFq_2_test.cpp @@ -132,7 +132,7 @@ TEST(RevMath, grad_2F1_derivs_match) { EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -/* + TEST(RevMath, grad2F1_2) { using stan::math::grad_2F1; using stan::math::grad_pFq_2; @@ -153,7 +153,6 @@ TEST(RevMath, grad2F1_2) { EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } - */ TEST(RevMath, grad2F1_3) { using stan::math::grad_2F1; @@ -281,7 +280,6 @@ TEST(RevMath, grad2F1_10) { EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -/* TEST(RevMath, grad2F1_11) { using stan::math::grad_2F1; using stan::math::grad_pFq_2; @@ -294,16 +292,16 @@ TEST(RevMath, grad2F1_11) { b_v << 2.70975; var z_v = 0.999696; - double std::get<0>(grad_2F1_tuple); - double std::get<2>(grad_2F1_tuple); - grad_2F1(std::get<0>(grad_2F1_tuple), std::get<2>(grad_2F1_tuple), a_v[0], -a_v[1], b_v[0], z_v); auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); + auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); -}*/ -/* + EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); +} + TEST(RevMath, F32_converges_by_z) { using stan::math::grad_F32; using stan::math::grad_pFq_2; @@ -328,7 +326,6 @@ TEST(RevMath, F32_converges_by_z) { EXPECT_FLOAT_EQ(g_calc[4], std::get<1>(grad_tuple)[1]); EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); } -*/ TEST(RevMath, grad_F32_double_sign_flip_1) { using stan::math::grad_F32; From af83ccde73245f9a0a990e42c5631111e652ea5f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Oct 2023 12:46:04 +0300 Subject: [PATCH 04/36] Log scale increment --- stan/math/prim/fun/grad_pFq_2.hpp | 44 +++++++++++++++---------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq_2.hpp b/stan/math/prim/fun/grad_pFq_2.hpp index 8fc7337f955..a73d754289d 100644 --- a/stan/math/prim/fun/grad_pFq_2.hpp +++ b/stan/math/prim/fun/grad_pFq_2.hpp @@ -8,6 +8,7 @@ #include #include #include +#include namespace stan { namespace math { @@ -61,24 +62,36 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, T_ArrayPartials log_digamma_b_k = log(abs(digamma_b)); IntArray digamma_a_k_signs = a_k_signs; IntArray digamma_b_k_signs = b_k_signs; - T_ArrayPartials a_infsum = T_ArrayPartials::Zero(a.size()); - T_ArrayPartials b_infsum = T_ArrayPartials::Zero(b.size()); + T_ArrayPartials a_log_infsum = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); + T_ArrayPartials b_log_infsum = T_ArrayPartials::Constant(b.size(), NEGATIVE_INFTY); + + IntArray a_infsum_sign = IntArray::Ones(a.size()); + IntArray b_infsum_sign = IntArray::Ones(b.size()); int k = 0; double curr_log_prec = 1.0; - double curr_prec = 1.0; double log_factorial_k = 0; while (k <= max_steps && curr_log_prec > log(precision)) { int base_sign = z_k_sign * phammer_a_k_signs.prod() * phammer_b_k_signs.prod(); T_partials log_base = (sum(log_phammer_a_k) + log_z_pow_k) - (log_factorial_k + sum(log_phammer_b_k)); + T_ArrayPartials a_grad = log_digamma_a_k + log_base; T_ArrayPartials b_grad = log_digamma_b_k + log_base; - + IntArray a_grad_signs = base_sign * digamma_a_k_signs; + IntArray b_grad_signs = base_sign * digamma_b_k_signs; + + for (int i = 0; i < a.size(); i++) { + std::forward_as_tuple(a_log_infsum(i), a_infsum_sign(i)) + = log_sum_exp_signed(a_log_infsum(i), a_infsum_sign(i), + a_grad(i), a_grad_signs(i)); + } + for (int i = 0; i < b.size(); i++) { + std::forward_as_tuple(b_log_infsum(i), b_infsum_sign(i)) + = log_sum_exp_signed(b_log_infsum(i), b_infsum_sign(i), + b_grad(i), b_grad_signs(i)); + } curr_log_prec = max(a_grad.maxCoeff(), b_grad.maxCoeff()); - a_infsum += exp(a_grad) * base_sign * digamma_a_k_signs; - b_infsum += exp(b_grad) * base_sign * digamma_b_k_signs; - log_phammer_a_k += log(abs(a_k)); log_phammer_b_k += log(abs(b_k)); phammer_a_k_signs *= a_k_signs; @@ -107,8 +120,8 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, T_partials pfq_p1_val = hypergeometric_pFq(a_val + 1, b_val + 1, z_val); T_return ret_tuple; - std::get<0>(ret_tuple) = a_infsum - digamma_a * pfq_val; - std::get<1>(ret_tuple) = digamma_b * pfq_val - b_infsum; + std::get<0>(ret_tuple) = a_infsum_sign * exp(a_log_infsum) - digamma_a * pfq_val; + std::get<1>(ret_tuple) = digamma_b * pfq_val - b_infsum_sign * exp(b_log_infsum); std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; return ret_tuple; } @@ -116,16 +129,3 @@ auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, } // namespace math } // namespace stan #endif - - - -/* std::cout << "k: " << k << "\n" - << "num: " << num << "\n" - << "den: " << den << "\n" - << "phammer_a_k: " << phammer_a_k.matrix().transpose() << "\n" - << "phammer_b_k: " << phammer_b_k.matrix().transpose() << "\n" - << "digamma_a_k: " << digamma_a_k.matrix().transpose() << "\n" - << "digamma_b_k: " << digamma_b_k.matrix().transpose() << "\n" - << "a: " << a_grad.matrix().transpose() << "\n" - << "b: " << b_grad.matrix().transpose() << "\n" - << std::endl; */ From c9cd18dfb8e98bcf79a9380f300ba45276ff0888 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Oct 2023 12:51:58 +0300 Subject: [PATCH 05/36] Replace primary definition --- stan/math/prim/fun.hpp | 1 - stan/math/prim/fun/grad_pFq.hpp | 422 ++++++-------------- stan/math/prim/fun/grad_pFq_2.hpp | 131 ------- test/unit/math/rev/fun/grad_pFq_2_test.cpp | 427 --------------------- test/unit/math/rev/fun/grad_pFq_test.cpp | 12 +- 5 files changed, 114 insertions(+), 879 deletions(-) delete mode 100644 stan/math/prim/fun/grad_pFq_2.hpp delete mode 100644 test/unit/math/rev/fun/grad_pFq_2_test.cpp diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 364a56c6e85..11ba14123bd 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -130,7 +130,6 @@ #include #include #include -#include #include #include #include diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 0228f1c7481..12de42df362 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -2,314 +2,17 @@ #define STAN_MATH_PRIM_FUN_GRAD_PFQ_HPP #include -#include -#include -#include -#include -#include -#include -#include +#include #include -#include -#include +#include #include -#include -#include -#include -#include +#include +#include +#include +#include namespace stan { namespace math { -namespace internal { -/** - * Returns the gradient of generalized hypergeometric function wrt to the - * input arguments: - * \f$ _pF_q(a_1,...,a_p;b_1,...,b_q;z) \f$ - * - * The derivatives wrt a and b are defined using a Kampé de Fériet function, - * (https://en.wikipedia.org/wiki/Kamp%C3%A9_de_F%C3%A9riet_function). - * This is implemented below as an infinite sum (on the log scale) until - * convergence. - * - * \f$ \frac{\partial}{\partial a_1} = - * \frac{z\prod_{j=2}^pa_j}{\prod_{j=1}^qb_j} - * F_{q+1\:0\;1}^{p\:1\:2}\left(\begin{array}&a_1+1,...,a_p+1;1;1,a_1\\ - * 2, b_1+1,...,b_1+1;;a_1+1\end{array};z,z\right) \f$ - * - * \f$ \frac{\partial}{\partial b_1}= - * -\frac{z\prod_{j=1}^pa_j}{b_1\prod_{j=1}^qb_j} - * F_{q+1\:0\;1}^{p\:1\:2}\left(\begin{array}&a_1+1,...,a_p+1;1;1,b_1\\ - * 2, b_1+1,...,b_1+1;;b_1+1\end{array};z,z\right) \f$ - * - * \f$ \frac{\partial}{\partial z}= \frac{\prod_{j=1}^pa_j}{\prod_{j=1}^qb_j} - * {}_pF_q(a_1 + 1,...,a_p + 1; b_1 +1,...,b_q+1;z) \f$ - * - * @tparam calc_a Boolean for whether to calculate derivatives wrt to 'a' - * @tparam calc_b Boolean for whether to calculate derivatives wrt to 'b' - * @tparam calc_z Boolean for whether to calculate derivatives wrt to 'z' - * @tparam TupleT Type of tuple containing objects to evaluate gradients into - * @tparam Ta Eigen type with either one row or column at compile time - * @tparam Tb Eigen type with either one row or column at compile time - * @tparam Tz Scalar type - * @param[in] grad_tuple Tuple of references to evaluate gradients into - * @param[in] a Vector of 'a' arguments to function - * @param[in] b Vector of 'b' arguments to function - * @param[in] z Scalar z argument - * @param[in] precision Convergence criteria for infinite sum - * @param[in] outer_steps Maximum number of iterations for infinite sum - * @param[in] inner_steps Maximum number of iterations for infinite sum - * @return Generalized hypergeometric function - */ -template * = nullptr, - require_stan_scalar_t* = nullptr> -void grad_pFq_impl(TupleT&& grad_tuple, const Ta& a, const Tb& b, const Tz& z, - double precision, int outer_steps, int inner_steps) { - using std::max; - using scalar_t = return_type_t; - using Ta_plain = plain_type_t; - using Tb_plain = plain_type_t; - using T_vec = Eigen::Matrix; - ref_type_t a_ref_in = a; - ref_type_t b_ref_in = b; - - // Replace any zero inputs with the smallest number representable, so that - // taking log and aggregating does not return -Inf - Ta_plain a_ref - = (a_ref_in.array() == 0).select(EPSILON, a_ref_in.array()).matrix(); - Tb_plain b_ref - = (b_ref_in.array() == 0).select(EPSILON, b_ref_in.array()).matrix(); - int a_size = a.size(); - int b_size = b.size(); - - // Convergence criteria for the gradients are the same as for the - // Hypergeometric pFq itself - bool condition_1 = (a_size > (b_size + 1)) && (z != 0); - bool condition_2 = (a_size == (b_size + 1)) && (fabs(z) > 1); - - if (condition_1 || condition_2) { - std::stringstream msg; - msg << "hypergeometric pFq gradient does not meet convergence " - << "conditions with given arguments. " - << "a: " << a_ref << ", b: " << b_ref << ", " - << "z: " << z; - throw std::domain_error(msg.str()); - } - - // As the gradients will be aggregating on the log scale, we will track the - // the values and the signs separately - to avoid taking the log of a - // negative input - Eigen::VectorXi a_signs = sign(value_of_rec(a_ref)); - Eigen::VectorXi b_signs = sign(value_of_rec(b_ref)); - int z_sign = sign(value_of_rec(z)); - - Ta_plain ap1 = (a.array() + 1).matrix(); - Tb_plain bp1 = (b.array() + 1).matrix(); - scalar_type_t log_z = log(fabs(z)); - scalar_type_t a_prod = prod(a_ref); - scalar_type_t b_prod = prod(b_ref); - - // Only need the infinite sum for partials wrt a & b - if (calc_a || calc_b) { - double log_precision = log(precision); - - T_vec da_mn = T_vec::Constant(a_size, NEGATIVE_INFTY); - T_vec db_mn = T_vec::Constant(b_size, NEGATIVE_INFTY); - T_vec da = T_vec::Constant(a_size, 0.0); - T_vec db = T_vec::Constant(b_size, 0.0); - - int m = 0; - int n_iter = 2; - - double lgamma_mp1 = 0; - double log_phammer_1m = 0; - double log_phammer_2m = 0; - T_vec log_phammer_ap1_m = T_vec::Zero(ap1.size()); - T_vec log_phammer_bp1_m = T_vec::Zero(bp1.size()); - Tz log_z_m = 0; - Ta_plain ap1m = ap1; - Tb_plain bp1m = bp1; - - Ta_plain ap1n = ap1; - Tb_plain bp1n = bp1; - Ta_plain an = a_ref; - Tb_plain bn = b_ref; - Ta_plain ap1mn = ap1m; - Tb_plain bp1mn = bp1m; - - double log_phammer_1n; - double log_phammer_2_mpn; - double lgamma_np1; - T_vec log_phammer_an(a_size); - T_vec log_phammer_ap1_n(a_size); - T_vec log_phammer_bp1_n(b_size); - T_vec log_phammer_bn(b_size); - T_vec log_phammer_ap1_mpn(a_size); - T_vec log_phammer_bp1_mpn(b_size); - - int z_pow_m_sign = 1; - Eigen::VectorXi curr_signs_da(a_size); - Eigen::VectorXi curr_signs_db(b_size); - Eigen::VectorXi log_phammer_an_sign(a_size); - Eigen::VectorXi log_phammer_ap1n_sign(a_size); - Eigen::VectorXi log_phammer_bp1n_sign(b_size); - Eigen::VectorXi log_phammer_bn_sign(b_size); - Eigen::VectorXi log_phammer_ap1mpn_sign(a_size); - Eigen::VectorXi log_phammer_bp1mpn_sign(b_size); - Eigen::VectorXi log_phammer_ap1m_sign = Eigen::VectorXi::Ones(a_size); - Eigen::VectorXi log_phammer_bp1m_sign = Eigen::VectorXi::Ones(b_size); - - // If the inner loop converges in 1 iteration, then the sum has coverged - // and another iteration of the outer loop is not needed - while ((n_iter > 1) && (m < outer_steps)) { - ap1n = ap1; - bp1n = bp1; - an = a_ref; - bn = b_ref; - ap1mn = ap1m; - bp1mn = bp1m; - - int n = 0; - Tz log_z_mn = log_z_m; - int z_pow_mn_sign = z_pow_m_sign; - scalar_t inner_diff = 0; - lgamma_np1 = 0; - - log_phammer_1n = 0; - log_phammer_an.setZero(); - log_phammer_ap1_n.setZero(); - log_phammer_bp1_n.setZero(); - log_phammer_bn.setZero(); - log_phammer_ap1_mpn = log_phammer_ap1_m; - log_phammer_bp1_mpn = log_phammer_bp1_m; - log_phammer_2_mpn = log_phammer_2m; - log_phammer_an_sign.setOnes(); - log_phammer_ap1n_sign.setOnes(); - log_phammer_bp1n_sign.setOnes(); - log_phammer_bn_sign.setOnes(); - log_phammer_ap1mpn_sign = log_phammer_ap1m_sign; - log_phammer_bp1mpn_sign = log_phammer_bp1m_sign; - - while ((inner_diff > log_precision) && (n < inner_steps)) { - // Numerator term - scalar_t term1_mn = log_z_mn + sum(log_phammer_ap1_mpn) + log_phammer_1m - + log_phammer_1n; - // Denominator term - scalar_t term2_mn = lgamma_mp1 + lgamma_np1 + sum(log_phammer_bp1_mpn) - + log_phammer_2_mpn; - int base_sign = z_pow_mn_sign * log_phammer_ap1mpn_sign.prod() - * log_phammer_bp1mpn_sign.prod(); - - if (calc_a) { - // Division (on log scale) for the a & b partials - // Determine signs of each element - curr_signs_da = (base_sign * log_phammer_an_sign.array() - * log_phammer_ap1n_sign.array()) - .matrix(); - da_mn = (term1_mn + log_phammer_an.array()) - - (term2_mn + log_phammer_ap1_n.array()); - - // Aggregate the sums on the natural scale, so that the sign can be - // applied before aggregation - da += exp(da_mn).cwiseProduct(curr_signs_da); - } - - if (calc_b) { - curr_signs_db = (base_sign * log_phammer_bn_sign.array() - * log_phammer_bp1n_sign.array()) - .matrix(); - db_mn = (term1_mn + log_phammer_bn.array()) - - (term2_mn + log_phammer_bp1_n.array()); - db += exp(db_mn).cwiseProduct(curr_signs_db); - } - - // Series convergence assessed by whether the maximum term is - // smaller than the specified criteria (precision) - inner_diff = max(da_mn.maxCoeff(), db_mn.maxCoeff()); - - // Increment the input arguments and rising factorials - log_z_mn += log_z; - log_phammer_1n += log1p(n); - log_phammer_2_mpn += log(2 + m + n); - - log_phammer_ap1_n.array() - += log(math::fabs((ap1n.array() == 0).select(1.0, ap1n.array()))); - log_phammer_bp1_n.array() - += log(math::fabs((bp1n.array() == 0).select(1.0, bp1n.array()))); - log_phammer_an.array() - += log(math::fabs((an.array() == 0).select(1.0, an.array()))); - log_phammer_bn.array() - += log(math::fabs((bn.array() == 0).select(1.0, bn.array()))); - log_phammer_ap1_mpn.array() - += log(math::fabs((ap1mn.array() == 0).select(1.0, ap1mn.array()))); - log_phammer_bp1_mpn.array() - += log(math::fabs((bp1mn.array() == 0).select(1.0, bp1mn.array()))); - - z_pow_mn_sign *= z_sign; - log_phammer_ap1n_sign.array() *= sign(value_of_rec(ap1n)).array(); - log_phammer_bp1n_sign.array() *= sign(value_of_rec(bp1n)).array(); - log_phammer_an_sign.array() *= sign(value_of_rec(an)).array(); - log_phammer_bn_sign.array() *= sign(value_of_rec(bn)).array(); - log_phammer_ap1mpn_sign.array() *= sign(value_of_rec(ap1mn)).array(); - log_phammer_bp1mpn_sign.array() *= sign(value_of_rec(bp1mn)).array(); - - n += 1; - lgamma_np1 += log(n); - ap1n.array() += 1; - bp1n.array() += 1; - an.array() += 1; - bn.array() += 1; - ap1mn.array() += 1; - bp1mn.array() += 1; - } - - z_pow_m_sign *= z_sign; - - n_iter = n; - - log_z_m += log_z; - log_phammer_1m += log1p(m); - log_phammer_2m += log(2 + m); - log_phammer_ap1_m += log(math::fabs(ap1m)); - log_phammer_ap1m_sign.array() *= sign(value_of_rec(ap1m)).array(); - log_phammer_bp1_m += log(math::fabs(bp1m)); - log_phammer_bp1m_sign.array() *= sign(value_of_rec(bp1m)).array(); - - m += 1; - - lgamma_mp1 += log(m); - ap1m.array() += 1; - bp1m.array() += 1; - ap1mn.array() += 1; - bp1mn.array() += 1; - } - - if (m == outer_steps) { - throw_domain_error("grad_pFq", "k (internal counter)", outer_steps, - "exceeded ", - " iterations, hypergeometric function gradient " - "did not converge."); - } - - if (calc_a) { - auto pre_mult_a = (z * a_prod / a_ref.array() / b_prod).matrix(); - std::get<0>(grad_tuple) = std::move(pre_mult_a.cwiseProduct(da)); - } - - if (calc_b) { - auto pre_mult_b = ((z * a_prod) / (b.array() * b_prod)).matrix(); - std::get<1>(grad_tuple) = std::move(-pre_mult_b.cwiseProduct(db)); - } - } - - if (calc_z) { - std::get<2>(grad_tuple) - = std::move((a_prod / b_prod) * hypergeometric_pFq(ap1, bp1, z)); - } -} -} // namespace internal - /** * Wrapper function for calculating gradients for the generalized * hypergeometric function. The function always returns a tuple with @@ -329,17 +32,108 @@ void grad_pFq_impl(TupleT&& grad_tuple, const Ta& a, const Tb& b, const Tz& z, * @return Tuple of gradients */ template -auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, double precision = 1e-10, - int outer_steps = 1e6, int inner_steps = 1e6) { - using partials_t = partials_return_t; - std::tuple>, - promote_scalar_t>, - promote_scalar_t>> - ret_tuple; - internal::grad_pFq_impl::value, !is_constant::value, - !is_constant::value>( - ret_tuple, value_of(a), value_of(b), value_of(z), precision, outer_steps, - inner_steps); +auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, + double precision = 1e-16, + int max_steps = 1e8) { + using T_partials = partials_return_t; + using Ta_partials = promote_scalar_t>; + using Tb_partials = promote_scalar_t>; + using Tz_partials = promote_scalar_t>; + using T_ArrayPartials = promote_scalar_t>; + using T_return = std::tuple; + using IntArray = Eigen::Array; + + T_ArrayPartials a_val = as_value_column_array_or_scalar(a); + T_ArrayPartials b_val = as_value_column_array_or_scalar(b); + Tz_partials z_val = value_of(z); + + IntArray a_signs = sign(value_of_rec(a_val)); + IntArray b_signs = sign(value_of_rec(b_val)); + int z_sign = sign(value_of_rec(z_val)); + + IntArray a_k_signs = a_signs; + IntArray b_k_signs = b_signs; + int z_k_sign = 1; + + T_ArrayPartials a_k = a_val; + T_ArrayPartials b_k = b_val; + Tz_partials log_z_pow_k = 0.0; + Tz_partials log_z_k = log(abs(z_val)); + + T_ArrayPartials log_phammer_a_k = T_ArrayPartials::Zero(a.size()); + T_ArrayPartials log_phammer_b_k = T_ArrayPartials::Zero(b.size()); + IntArray phammer_a_k_signs = IntArray::Ones(a.size()); + IntArray phammer_b_k_signs = IntArray::Ones(b.size()); + + T_ArrayPartials digamma_a = inv(a_k); + T_ArrayPartials digamma_b = inv(b_k); + T_ArrayPartials digamma_a_k = digamma_a; + T_ArrayPartials digamma_b_k = digamma_b; + T_ArrayPartials log_digamma_a_k = log(abs(digamma_a)); + T_ArrayPartials log_digamma_b_k = log(abs(digamma_b)); + IntArray digamma_a_k_signs = a_k_signs; + IntArray digamma_b_k_signs = b_k_signs; + T_ArrayPartials a_log_infsum = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); + T_ArrayPartials b_log_infsum = T_ArrayPartials::Constant(b.size(), NEGATIVE_INFTY); + + IntArray a_infsum_sign = IntArray::Ones(a.size()); + IntArray b_infsum_sign = IntArray::Ones(b.size()); + + int k = 0; + double curr_log_prec = 1.0; + double log_factorial_k = 0; + while (k <= max_steps && curr_log_prec > log(precision)) { + int base_sign = z_k_sign * phammer_a_k_signs.prod() * phammer_b_k_signs.prod(); + T_partials log_base = (sum(log_phammer_a_k) + log_z_pow_k) + - (log_factorial_k + sum(log_phammer_b_k)); + + T_ArrayPartials a_grad = log_digamma_a_k + log_base; + T_ArrayPartials b_grad = log_digamma_b_k + log_base; + IntArray a_grad_signs = base_sign * digamma_a_k_signs; + IntArray b_grad_signs = base_sign * digamma_b_k_signs; + + for (int i = 0; i < a.size(); i++) { + std::forward_as_tuple(a_log_infsum(i), a_infsum_sign(i)) + = log_sum_exp_signed(a_log_infsum(i), a_infsum_sign(i), + a_grad(i), a_grad_signs(i)); + } + for (int i = 0; i < b.size(); i++) { + std::forward_as_tuple(b_log_infsum(i), b_infsum_sign(i)) + = log_sum_exp_signed(b_log_infsum(i), b_infsum_sign(i), + b_grad(i), b_grad_signs(i)); + } + curr_log_prec = max(a_grad.maxCoeff(), b_grad.maxCoeff()); + log_phammer_a_k += log(abs(a_k)); + log_phammer_b_k += log(abs(b_k)); + phammer_a_k_signs *= a_k_signs; + phammer_b_k_signs *= b_k_signs; + + log_z_pow_k += log_z_k; + z_k_sign *= z_sign; + + digamma_a_k += select(a_k == 0.0, 0.0, inv(a_k)); + digamma_b_k += select(b_k == 0.0, 0.0, inv(b_k)); + log_digamma_a_k = log(abs(digamma_a_k)); + log_digamma_b_k = log(abs(digamma_b_k)); + digamma_a_k_signs = sign(value_of_rec(digamma_a_k)); + digamma_b_k_signs = sign(value_of_rec(digamma_b_k)); + + a_k += 1.0; + b_k += 1.0; + a_k_signs = sign(value_of_rec(a_k)); + b_k_signs = sign(value_of_rec(b_k)); + + k += 1; + log_factorial_k += log(k); + } + + T_partials pfq_val = hypergeometric_pFq(a_val, b_val, z_val); + T_partials pfq_p1_val = hypergeometric_pFq(a_val + 1, b_val + 1, z_val); + + T_return ret_tuple; + std::get<0>(ret_tuple) = a_infsum_sign * exp(a_log_infsum) - digamma_a * pfq_val; + std::get<1>(ret_tuple) = digamma_b * pfq_val - b_infsum_sign * exp(b_log_infsum); + std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; return ret_tuple; } diff --git a/stan/math/prim/fun/grad_pFq_2.hpp b/stan/math/prim/fun/grad_pFq_2.hpp deleted file mode 100644 index a73d754289d..00000000000 --- a/stan/math/prim/fun/grad_pFq_2.hpp +++ /dev/null @@ -1,131 +0,0 @@ -#ifndef STAN_MATH_PRIM_FUN_GRAD_PFQ_2_HPP -#define STAN_MATH_PRIM_FUN_GRAD_PFQ_2_HPP - -#include -#include -#include -#include -#include -#include -#include -#include - -namespace stan { -namespace math { -namespace internal { - template - auto bin_sign(const T& x) { - decltype(auto) sign1 = sign(value_of_rec(x)); - return select(sign1 == 0, 1, sign1); - } -} - -template -auto grad_pFq_2(const Ta& a, const Tb& b, const Tz& z, - double precision = 1e-16, - int max_steps = 1e8) { - using T_partials = partials_return_t; - using Ta_partials = promote_scalar_t>; - using Tb_partials = promote_scalar_t>; - using Tz_partials = promote_scalar_t>; - using T_ArrayPartials = promote_scalar_t>; - using T_return = std::tuple; - using IntArray = Eigen::Array; - - T_ArrayPartials a_val = as_value_column_array_or_scalar(a); - T_ArrayPartials b_val = as_value_column_array_or_scalar(b); - Tz_partials z_val = value_of(z); - - IntArray a_signs = internal::bin_sign(a_val); - IntArray b_signs = internal::bin_sign(b_val); - int z_sign = internal::bin_sign(z_val); - - IntArray a_k_signs = a_signs; - IntArray b_k_signs = b_signs; - int z_k_sign = 1; - - T_ArrayPartials a_k = a_val; - T_ArrayPartials b_k = b_val; - Tz_partials log_z_pow_k = 0.0; - Tz_partials log_z_k = log(abs(z_val)); - - T_ArrayPartials log_phammer_a_k = T_ArrayPartials::Zero(a.size()); - T_ArrayPartials log_phammer_b_k = T_ArrayPartials::Zero(b.size()); - IntArray phammer_a_k_signs = IntArray::Ones(a.size()); - IntArray phammer_b_k_signs = IntArray::Ones(b.size()); - - T_ArrayPartials digamma_a = inv(a_k); - T_ArrayPartials digamma_b = inv(b_k); - T_ArrayPartials digamma_a_k = digamma_a; - T_ArrayPartials digamma_b_k = digamma_b; - T_ArrayPartials log_digamma_a_k = log(abs(digamma_a)); - T_ArrayPartials log_digamma_b_k = log(abs(digamma_b)); - IntArray digamma_a_k_signs = a_k_signs; - IntArray digamma_b_k_signs = b_k_signs; - T_ArrayPartials a_log_infsum = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); - T_ArrayPartials b_log_infsum = T_ArrayPartials::Constant(b.size(), NEGATIVE_INFTY); - - IntArray a_infsum_sign = IntArray::Ones(a.size()); - IntArray b_infsum_sign = IntArray::Ones(b.size()); - - int k = 0; - double curr_log_prec = 1.0; - double log_factorial_k = 0; - while (k <= max_steps && curr_log_prec > log(precision)) { - int base_sign = z_k_sign * phammer_a_k_signs.prod() * phammer_b_k_signs.prod(); - T_partials log_base = (sum(log_phammer_a_k) + log_z_pow_k) - - (log_factorial_k + sum(log_phammer_b_k)); - - T_ArrayPartials a_grad = log_digamma_a_k + log_base; - T_ArrayPartials b_grad = log_digamma_b_k + log_base; - IntArray a_grad_signs = base_sign * digamma_a_k_signs; - IntArray b_grad_signs = base_sign * digamma_b_k_signs; - - for (int i = 0; i < a.size(); i++) { - std::forward_as_tuple(a_log_infsum(i), a_infsum_sign(i)) - = log_sum_exp_signed(a_log_infsum(i), a_infsum_sign(i), - a_grad(i), a_grad_signs(i)); - } - for (int i = 0; i < b.size(); i++) { - std::forward_as_tuple(b_log_infsum(i), b_infsum_sign(i)) - = log_sum_exp_signed(b_log_infsum(i), b_infsum_sign(i), - b_grad(i), b_grad_signs(i)); - } - curr_log_prec = max(a_grad.maxCoeff(), b_grad.maxCoeff()); - log_phammer_a_k += log(abs(a_k)); - log_phammer_b_k += log(abs(b_k)); - phammer_a_k_signs *= a_k_signs; - phammer_b_k_signs *= b_k_signs; - - log_z_pow_k += log_z_k; - z_k_sign *= z_sign; - - digamma_a_k += select(a_k == 0.0, 0.0, inv(a_k)); - digamma_b_k += select(b_k == 0.0, 0.0, inv(b_k)); - log_digamma_a_k = log(abs(digamma_a_k)); - log_digamma_b_k = log(abs(digamma_b_k)); - digamma_a_k_signs = internal::bin_sign(digamma_a_k); - digamma_b_k_signs = internal::bin_sign(digamma_b_k); - - a_k += 1.0; - b_k += 1.0; - a_k_signs = internal::bin_sign(a_k); - b_k_signs = internal::bin_sign(b_k); - - k += 1; - log_factorial_k += log(k); - } - - T_partials pfq_val = hypergeometric_pFq(a_val, b_val, z_val); - T_partials pfq_p1_val = hypergeometric_pFq(a_val + 1, b_val + 1, z_val); - - T_return ret_tuple; - std::get<0>(ret_tuple) = a_infsum_sign * exp(a_log_infsum) - digamma_a * pfq_val; - std::get<1>(ret_tuple) = digamma_b * pfq_val - b_infsum_sign * exp(b_log_infsum); - std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; - return ret_tuple; -} - -} // namespace math -} // namespace stan -#endif diff --git a/test/unit/math/rev/fun/grad_pFq_2_test.cpp b/test/unit/math/rev/fun/grad_pFq_2_test.cpp deleted file mode 100644 index 2db6bb8ee47..00000000000 --- a/test/unit/math/rev/fun/grad_pFq_2_test.cpp +++ /dev/null @@ -1,427 +0,0 @@ -#include -#include -#include - -TEST(RevMath, grad_2F2) { - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_d; - using stan::math::vector_v; - - vector_v a_v(2); - vector_d a_d(2); - a_v << 4, 2; - a_d << 4, 2; - - vector_v b_v(2); - vector_d b_d(2); - b_v << 6, 3; - b_d << 6, 3; - - var z_v = 4; - double z_d = 4; - - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); - - EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); - - EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); - - grad_tuple = grad_pFq_2(a_v, b_d, z_d); - - EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); - - grad_tuple = grad_pFq_2(a_d, b_v, z_d); - - EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); - - grad_tuple = grad_pFq_2(a_d, b_d, z_v); - - EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); - - grad_tuple = grad_pFq_2(a_v, b_v, z_d); - - EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); - - grad_tuple = grad_pFq_2(a_v, b_d, z_v); - - EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); - - grad_tuple = grad_pFq_2(a_d, b_v, z_v); - - EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); -} -TEST(RevMath, grad_2F3) { - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << 2, 3; - vector_v b_v(3); - b_v << 2, 4, 5; - var z_v = 1; - - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(0.08377717301140296, std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(0.05615450733193106, std::get<0>(grad_tuple)[1]); - - EXPECT_FLOAT_EQ(-0.08377717301140296, std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(-0.04225296806591615, std::get<1>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(-0.03387575989873739, std::get<1>(grad_tuple)[2]); - - EXPECT_FLOAT_EQ(0.1712340452215524, std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad_4F3) { - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(4); - a_v << 1, 2, 3, 4; - vector_v b_v(3); - b_v << 5, 6, 7; - var z_v = 1; - - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(0.1587098625610631, std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(0.08249338029396848, std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(0.05611368752226367, std::get<0>(grad_tuple)[2]); - EXPECT_FLOAT_EQ(0.04261209968272329, std::get<0>(grad_tuple)[3]); - - EXPECT_FLOAT_EQ(-0.03438035893346993, std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(-0.02882791253333995, std::get<1>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(-0.02482622713079612, std::get<1>(grad_tuple)[2]); - - EXPECT_FLOAT_EQ(0.1800529055890911, std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad_2F1_derivs_match) { - using stan::math::grad_2F1; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << 1, 1; - vector_v b_v(1); - b_v << 1; - var z_v = 0.6; - - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad2F1_2) { - using stan::math::grad_2F1; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << 1, 31; - vector_v b_v(1); - b_v << 41; - var z_v = 1; - - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad2F1_3) { - using stan::math::grad_2F1; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << 1.0, -2.1; - vector_v b_v(1); - b_v << 41.0; - var z_v = 1.0; - - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad2F1_6) { - using stan::math::grad_2F1; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << 1, -0.5; - vector_v b_v(1); - b_v << 10.6; - var z_v = 0.3; - - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad2F1_7) { - using stan::math::grad_2F1; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << 1, -0.5; - vector_v b_v(1); - b_v << 10; - var z_v = 0.3; - - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad2F1_8) { - using stan::math::grad_2F1; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << -0.5, -4.5; - vector_v b_v(1); - b_v << 11; - var z_v = 0.3; - - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad2F1_9) { - using stan::math::grad_2F1; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << -0.5, -4.5; - vector_v b_v(1); - b_v << -3.2; - var z_v = 0.9; - - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad2F1_10) { - using stan::math::grad_2F1; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << 2, 1; - vector_v b_v(1); - b_v << 2; - var z_v = 0.4; - - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad2F1_11) { - using stan::math::grad_2F1; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << 3.70975, 1; - vector_v b_v(1); - b_v << 2.70975; - var z_v = 0.999696; - - - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); -} - -TEST(RevMath, F32_converges_by_z) { - using stan::math::grad_F32; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(3); - a_v << 1.0, 1.0, 1.0; - vector_v b_v(2); - b_v << 1.0, 1.0; - var z_v = 0.6; - - double g_calc[6]; - - grad_F32(g_calc, 1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 1e-10); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(g_calc[2], std::get<0>(grad_tuple)[2]); - EXPECT_FLOAT_EQ(g_calc[3], std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(g_calc[4], std::get<1>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad_F32_double_sign_flip_1) { - using stan::math::grad_F32; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(3); - a_v << 1.0, -0.5, -2.5; - vector_v b_v(2); - b_v << 10.0, 1.0; - var z_v = 0.3; - - double g_calc[6]; - - grad_F32(g_calc, 1.0, -.5, -2.5, 10.0, 1.0, 0.3, 1e-10); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(g_calc[2], std::get<0>(grad_tuple)[2]); - EXPECT_FLOAT_EQ(g_calc[3], std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(g_calc[4], std::get<1>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad_F32_double_sign_flip_2) { - using stan::math::grad_F32; - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(3); - a_v << 1.0, -0.5, -4.5; - vector_v b_v(2); - b_v << 10.0, 1.0; - var z_v = 0.3; - - double g_calc[6]; - - grad_F32(g_calc, 1.0, -.5, -4.5, 10.0, 1.0, 0.3, 1e-10); - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(g_calc[2], std::get<0>(grad_tuple)[2]); - EXPECT_FLOAT_EQ(g_calc[3], std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(g_calc[4], std::get<1>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); -} - -TEST(RevMath, grad_2F1_negative_z) { - using stan::math::grad_pFq_2; - using stan::math::var; - using stan::math::vector_v; - - vector_v a_v(2); - a_v << 3.70975, 1.0; - vector_v b_v(1); - b_v << 2.70975; - var z_v = -0.2; - - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(-0.0488658806159776, std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(-0.193844936204681, std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(0.0677809985598383, std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(0.865295247272367, std::get<2>(grad_tuple)); -} - - -TEST(RevMath, grad_3F2_cross_zero) { - using stan::math::grad_F32; - using stan::math::grad_pFq_2; - using stan::math::hypergeometric_3F2; - using stan::math::var; - using stan::math::vector_d; - using stan::math::vector_v; - - vector_v a_v(3); - a_v << 1, 1, -1; - - vector_v b_v(2); - b_v << 2, 2; - - var z_v = 0.292893; - - double g_calc[6]; - grad_F32(g_calc, 1.0, 1.0, -1.0, 2.0, 2.0, 0.292893); - - auto grad_tuple = grad_pFq_2(a_v, b_v, z_v); - - EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(g_calc[2], std::get<0>(grad_tuple)[2]); - EXPECT_FLOAT_EQ(g_calc[3], std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(g_calc[4], std::get<1>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); -} diff --git a/test/unit/math/rev/fun/grad_pFq_test.cpp b/test/unit/math/rev/fun/grad_pFq_test.cpp index 4e66e0273be..b6fe6e656b6 100644 --- a/test/unit/math/rev/fun/grad_pFq_test.cpp +++ b/test/unit/math/rev/fun/grad_pFq_test.cpp @@ -279,7 +279,7 @@ TEST(RevMath, grad2F1_10) { EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -/* + TEST(RevMath, grad2F1_11) { using stan::math::grad_2F1; using stan::math::grad_pFq; @@ -292,15 +292,15 @@ TEST(RevMath, grad2F1_11) { b_v << 2.70975; var z_v = 0.999696; - double std::get<0>(grad_2F1_tuple); - double std::get<2>(grad_2F1_tuple); - grad_2F1(std::get<0>(grad_2F1_tuple), std::get<2>(grad_2F1_tuple), a_v[0], -a_v[1], b_v[0], z_v); auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); + auto grad_tuple = grad_pFq(a_v, b_v, z_v); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); EXPECT_FLOAT_EQ(std::get<2>(grad_2F1_tuple), std::get<1>(grad_tuple)[0]); -}*/ + EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); +} TEST(RevMath, F32_converges_by_z) { using stan::math::grad_F32; From d5ce3db8ee7320e3868afc72d9ff8822c6af0b3f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Oct 2023 13:07:13 +0300 Subject: [PATCH 06/36] Don't calc for non-autodiff --- stan/math/prim/fun/grad_pFq.hpp | 93 ++++++++++++++++++++------------- 1 file changed, 58 insertions(+), 35 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 12de42df362..b06d12d0d3a 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -36,15 +36,16 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, double precision = 1e-16, int max_steps = 1e8) { using T_partials = partials_return_t; - using Ta_partials = promote_scalar_t>; - using Tb_partials = promote_scalar_t>; - using Tz_partials = promote_scalar_t>; + using Ta_partials = partials_return_t; + using Tb_partials = partials_return_t; + using Tz_partials = partials_return_t; using T_ArrayPartials = promote_scalar_t>; - using T_return = std::tuple; using IntArray = Eigen::Array; + using Ta_PartialsArray = promote_scalar_t>; + using Tb_PartialsArray = promote_scalar_t>; - T_ArrayPartials a_val = as_value_column_array_or_scalar(a); - T_ArrayPartials b_val = as_value_column_array_or_scalar(b); + Ta_PartialsArray a_val = as_value_column_array_or_scalar(a); + Tb_PartialsArray b_val = as_value_column_array_or_scalar(b); Tz_partials z_val = value_of(z); IntArray a_signs = sign(value_of_rec(a_val)); @@ -55,30 +56,35 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, IntArray b_k_signs = b_signs; int z_k_sign = 1; - T_ArrayPartials a_k = a_val; - T_ArrayPartials b_k = b_val; + Ta_PartialsArray a_k = a_val; + Tb_PartialsArray b_k = b_val; Tz_partials log_z_pow_k = 0.0; Tz_partials log_z_k = log(abs(z_val)); - T_ArrayPartials log_phammer_a_k = T_ArrayPartials::Zero(a.size()); - T_ArrayPartials log_phammer_b_k = T_ArrayPartials::Zero(b.size()); + Ta_PartialsArray log_phammer_a_k = T_ArrayPartials::Zero(a.size()); + Tb_PartialsArray log_phammer_b_k = T_ArrayPartials::Zero(b.size()); IntArray phammer_a_k_signs = IntArray::Ones(a.size()); IntArray phammer_b_k_signs = IntArray::Ones(b.size()); - T_ArrayPartials digamma_a = inv(a_k); - T_ArrayPartials digamma_b = inv(b_k); - T_ArrayPartials digamma_a_k = digamma_a; - T_ArrayPartials digamma_b_k = digamma_b; - T_ArrayPartials log_digamma_a_k = log(abs(digamma_a)); - T_ArrayPartials log_digamma_b_k = log(abs(digamma_b)); + Ta_PartialsArray digamma_a = inv(a_k); + Tb_PartialsArray digamma_b = inv(b_k); + Ta_PartialsArray digamma_a_k = digamma_a; + Tb_PartialsArray digamma_b_k = digamma_b; + Ta_PartialsArray log_digamma_a_k = log(abs(digamma_a)); + Tb_PartialsArray log_digamma_b_k = log(abs(digamma_b)); IntArray digamma_a_k_signs = a_k_signs; IntArray digamma_b_k_signs = b_k_signs; + T_ArrayPartials a_log_infsum = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); T_ArrayPartials b_log_infsum = T_ArrayPartials::Constant(b.size(), NEGATIVE_INFTY); - IntArray a_infsum_sign = IntArray::Ones(a.size()); IntArray b_infsum_sign = IntArray::Ones(b.size()); + T_ArrayPartials a_grad = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); + IntArray a_grad_signs = IntArray::Ones(a.size()); + T_ArrayPartials b_grad = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); + IntArray b_grad_signs = IntArray::Ones(a.size()); + int k = 0; double curr_log_prec = 1.0; double log_factorial_k = 0; @@ -87,21 +93,28 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, T_partials log_base = (sum(log_phammer_a_k) + log_z_pow_k) - (log_factorial_k + sum(log_phammer_b_k)); - T_ArrayPartials a_grad = log_digamma_a_k + log_base; - T_ArrayPartials b_grad = log_digamma_b_k + log_base; - IntArray a_grad_signs = base_sign * digamma_a_k_signs; - IntArray b_grad_signs = base_sign * digamma_b_k_signs; + if (!is_constant_all::value) { + a_grad = log_digamma_a_k + log_base; + a_grad_signs = base_sign * digamma_a_k_signs; - for (int i = 0; i < a.size(); i++) { - std::forward_as_tuple(a_log_infsum(i), a_infsum_sign(i)) - = log_sum_exp_signed(a_log_infsum(i), a_infsum_sign(i), - a_grad(i), a_grad_signs(i)); + for (int i = 0; i < a.size(); i++) { + std::forward_as_tuple(a_log_infsum(i), a_infsum_sign(i)) + = log_sum_exp_signed(a_log_infsum(i), a_infsum_sign(i), + a_grad(i), a_grad_signs(i)); + } } - for (int i = 0; i < b.size(); i++) { - std::forward_as_tuple(b_log_infsum(i), b_infsum_sign(i)) - = log_sum_exp_signed(b_log_infsum(i), b_infsum_sign(i), - b_grad(i), b_grad_signs(i)); + + if (!is_constant_all::value) { + b_grad = log_digamma_b_k + log_base; + b_grad_signs = base_sign * digamma_b_k_signs; + + for (int i = 0; i < b.size(); i++) { + std::forward_as_tuple(b_log_infsum(i), b_infsum_sign(i)) + = log_sum_exp_signed(b_log_infsum(i), b_infsum_sign(i), + b_grad(i), b_grad_signs(i)); + } } + curr_log_prec = max(a_grad.maxCoeff(), b_grad.maxCoeff()); log_phammer_a_k += log(abs(a_k)); log_phammer_b_k += log(abs(b_k)); @@ -127,13 +140,23 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, log_factorial_k += log(k); } - T_partials pfq_val = hypergeometric_pFq(a_val, b_val, z_val); - T_partials pfq_p1_val = hypergeometric_pFq(a_val + 1, b_val + 1, z_val); + std::tuple>, + promote_scalar_t>, + T_partials> ret_tuple; - T_return ret_tuple; - std::get<0>(ret_tuple) = a_infsum_sign * exp(a_log_infsum) - digamma_a * pfq_val; - std::get<1>(ret_tuple) = digamma_b * pfq_val - b_infsum_sign * exp(b_log_infsum); - std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; + if (!is_constant_all::value) { + T_partials pfq_val = hypergeometric_pFq(a_val, b_val, z_val); + if (!is_constant_all::value) { + std::get<0>(ret_tuple) = a_infsum_sign * exp(a_log_infsum) - digamma_a * pfq_val; + } + if (!is_constant_all::value) { + std::get<1>(ret_tuple) = digamma_b * pfq_val - b_infsum_sign * exp(b_log_infsum); + } + } + if (!is_constant_all::value) { + T_partials pfq_p1_val = hypergeometric_pFq(a_val + 1, b_val + 1, z_val); + std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; + } return ret_tuple; } From 668f7fdd14833ca6ce704fb8adfb81940e061324 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Oct 2023 18:29:42 +0300 Subject: [PATCH 07/36] Initial mix working --- stan/math/prim/fun/grad_pFq.hpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index b06d12d0d3a..7f915e8816b 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -35,6 +35,7 @@ template auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, double precision = 1e-16, int max_steps = 1e8) { + using std::max; using T_partials = partials_return_t; using Ta_partials = partials_return_t; using Tb_partials = partials_return_t; @@ -61,8 +62,8 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, Tz_partials log_z_pow_k = 0.0; Tz_partials log_z_k = log(abs(z_val)); - Ta_PartialsArray log_phammer_a_k = T_ArrayPartials::Zero(a.size()); - Tb_PartialsArray log_phammer_b_k = T_ArrayPartials::Zero(b.size()); + Ta_PartialsArray log_phammer_a_k = Ta_PartialsArray::Zero(a.size()); + Tb_PartialsArray log_phammer_b_k = Tb_PartialsArray::Zero(b.size()); IntArray phammer_a_k_signs = IntArray::Ones(a.size()); IntArray phammer_b_k_signs = IntArray::Ones(b.size()); @@ -81,17 +82,17 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, IntArray b_infsum_sign = IntArray::Ones(b.size()); T_ArrayPartials a_grad = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); - IntArray a_grad_signs = IntArray::Ones(a.size()); T_ArrayPartials b_grad = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); + IntArray a_grad_signs = IntArray::Ones(a.size()); IntArray b_grad_signs = IntArray::Ones(a.size()); int k = 0; - double curr_log_prec = 1.0; + T_partials curr_log_prec = 1.0; double log_factorial_k = 0; while (k <= max_steps && curr_log_prec > log(precision)) { int base_sign = z_k_sign * phammer_a_k_signs.prod() * phammer_b_k_signs.prod(); T_partials log_base = (sum(log_phammer_a_k) + log_z_pow_k) - - (log_factorial_k + sum(log_phammer_b_k)); + - (log_factorial_k + sum(log_phammer_b_k)); if (!is_constant_all::value) { a_grad = log_digamma_a_k + log_base; @@ -145,7 +146,7 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, T_partials> ret_tuple; if (!is_constant_all::value) { - T_partials pfq_val = hypergeometric_pFq(a_val, b_val, z_val); + T_partials pfq_val = hypergeometric_pFq(a_val.matrix(), b_val.matrix(), z_val); if (!is_constant_all::value) { std::get<0>(ret_tuple) = a_infsum_sign * exp(a_log_infsum) - digamma_a * pfq_val; } @@ -154,7 +155,7 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, } } if (!is_constant_all::value) { - T_partials pfq_p1_val = hypergeometric_pFq(a_val + 1, b_val + 1, z_val); + T_partials pfq_p1_val = hypergeometric_pFq((a_val + 1).matrix(), (b_val + 1).matrix(), z_val); std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; } return ret_tuple; From d98a289f82636b396417fa77cbc00bd9cd125c63 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Oct 2023 19:12:08 +0300 Subject: [PATCH 08/36] Begin debugging autodiff --- stan/math/mix/functor/grad_hessian.hpp | 10 ++++--- stan/math/rev/fun/hypergeometric_pFq.hpp | 12 ++++---- .../math/mix/fun/hypergeometric_pFq_test.cpp | 28 +++++++++++++++++++ test/unit/math/test_ad.hpp | 8 ++++-- 4 files changed, 45 insertions(+), 13 deletions(-) diff --git a/stan/math/mix/functor/grad_hessian.hpp b/stan/math/mix/functor/grad_hessian.hpp index 9fcf47ba1f6..f9ff5fa27f0 100644 --- a/stan/math/mix/functor/grad_hessian.hpp +++ b/stan/math/mix/functor/grad_hessian.hpp @@ -50,8 +50,10 @@ void grad_hessian( int d = x.size(); H.resize(d, d); grad_H.resize(d, Matrix(d, d)); - for (int i = 0; i < d; ++i) { - for (int j = i; j < d; ++j) { + //for (int i = 0; i < d; ++i) { + int i = 0; + //for (int j = i; j < d; ++j) { + int j = 0; // Run nested autodiff in this scope nested_rev_autodiff nested; @@ -68,8 +70,8 @@ void grad_hessian( grad_H[i](j, k) = x_ffvar(k).val_.val_.adj(); grad_H[j](i, k) = grad_H[i](j, k); } - } - } + //} + //} } } // namespace math diff --git a/stan/math/rev/fun/hypergeometric_pFq.hpp b/stan/math/rev/fun/hypergeometric_pFq.hpp index 281a2f79e62..360deabab89 100644 --- a/stan/math/rev/fun/hypergeometric_pFq.hpp +++ b/stan/math/rev/fun/hypergeometric_pFq.hpp @@ -32,16 +32,16 @@ inline var hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { [arena_a, arena_b, z](auto& vi) mutable { auto grad_tuple = grad_pFq(arena_a, arena_b, z); if (!is_constant::value) { - forward_as>(arena_a).adj() - += vi.adj() * std::get<0>(grad_tuple); + //forward_as>(arena_a).adj() + // += vi.adj() * std::get<0>(grad_tuple); } if (!is_constant::value) { - forward_as>(arena_b).adj() - += vi.adj() * std::get<1>(grad_tuple); + //forward_as>(arena_b).adj() + // += vi.adj() * std::get<1>(grad_tuple); } if (!is_constant::value) { - forward_as>(z).adj() - += vi.adj() * std::get<2>(grad_tuple); + //forward_as>(z).adj() + // += vi.adj() * std::get<2>(grad_tuple); } }); } diff --git a/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp b/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp index f11971c8fbc..afec57accd6 100644 --- a/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp +++ b/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp @@ -1,6 +1,33 @@ #include #include +TEST(mathMixScalFun, hyper_pfq) { + using stan::math::hypergeometric_pFq; + auto f = [](const auto& a, const auto& b, const auto& z) { + using stan::math::hypergeometric_pFq; + return hypergeometric_pFq(a, b, z); + }; + + Eigen::VectorXd in1(2); + in1 << 4, 2; + Eigen::VectorXd in2(2); + in2 << 6, 3; + double z = 4; + + Eigen::VectorXd x = stan::math::to_vector(stan::test::serialize(in1, in2, z)); + auto serial_functor = [&](const auto& v) { + auto v_deserializer = stan::test::to_deserializer(v); + return f(v_deserializer.read(in1), v_deserializer.read(in2), v_deserializer.read(z)); + }; + double fx_ad; + Eigen::MatrixXd H_ad; + std::vector grad_H_ad; + stan::math::grad_hessian(serial_functor, x, fx_ad, H_ad, grad_H_ad); + //std::cout << H_ad << "\n" << std::endl; + std::cout << grad_H_ad[0] << std::endl; + //stan::test::expect_ad(f, in1, in2, z); +} +/* TEST(mixScalFun, grad_2F2_ffv) { using stan::math::fvar; using stan::math::hypergeometric_pFq; @@ -48,3 +75,4 @@ TEST(mixScalFun, grad_2F2_ffv) { // double, double, fvar EXPECT_FLOAT_EQ(hypergeometric_pFq(d_a, d_b, ffv_z).val_.d_.val(), z_adj); } + */ diff --git a/test/unit/math/test_ad.hpp b/test/unit/math/test_ad.hpp index 6144336ad84..f9a4dabc96c 100644 --- a/test/unit/math/test_ad.hpp +++ b/test/unit/math/test_ad.hpp @@ -307,9 +307,11 @@ void test_grad_hessian(const ad_tolerances& tols, const F& f, expect_near_rel("grad_hessian() Hessian", H_fd, H_ad, tols.grad_hessian_hessian_); EXPECT_EQ(x.size(), grad_H_fd.size()); - for (size_t i = 0; i < grad_H_fd.size(); ++i) - expect_near_rel("grad_hessian() grad Hessian", grad_H_fd[i], grad_H_ad[i], - tols.grad_hessian_grad_hessian_); + for (size_t i = 0; i < grad_H_fd.size(); ++i) { + std::cout << i << "\n\n" << grad_H_ad[i] << "\n" << std::endl; + //expect_near_rel("grad_hessian() grad Hessian", grad_H_fd[i], grad_H_ad[i], + // tols.grad_hessian_grad_hessian_); + } } #endif From 0795e5d25286f49a17b97ab52456446d12f5c1ce Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 6 Oct 2023 15:37:39 +0300 Subject: [PATCH 09/36] Initial ad testing working --- stan/math/fwd/fun/hypergeometric_pFq.hpp | 6 +- stan/math/mix/functor/grad_hessian.hpp | 10 +-- stan/math/prim/fun/grad_pFq.hpp | 88 +++++++++++-------- stan/math/rev/fun/hypergeometric_pFq.hpp | 20 ++--- .../math/mix/fun/hypergeometric_pFq_test.cpp | 13 +-- test/unit/math/test_ad.hpp | 8 +- 6 files changed, 71 insertions(+), 74 deletions(-) diff --git a/stan/math/fwd/fun/hypergeometric_pFq.hpp b/stan/math/fwd/fun/hypergeometric_pFq.hpp index 9738557d9e5..75e0e229768 100644 --- a/stan/math/fwd/fun/hypergeometric_pFq.hpp +++ b/stan/math/fwd/fun/hypergeometric_pFq.hpp @@ -29,7 +29,8 @@ inline return_type_t hypergeometric_pFq(const Ta& a, const Tb& b, using fvar_t = return_type_t; ref_type_t a_ref = a; ref_type_t b_ref = b; - auto grad_tuple = grad_pFq(a_ref, b_ref, z); + auto pfq_val = hypergeometric_pFq(value_of(a_ref), value_of(b_ref), value_of(z)); + auto grad_tuple = grad_pFq(pfq_val, a_ref, b_ref, z); typename fvar_t::Scalar grad = 0; @@ -46,8 +47,7 @@ inline return_type_t hypergeometric_pFq(const Ta& a, const Tb& b, * std::get<2>(grad_tuple); } - return fvar_t( - hypergeometric_pFq(value_of(a_ref), value_of(b_ref), value_of(z)), grad); + return fvar_t(pfq_val, grad); } } // namespace math diff --git a/stan/math/mix/functor/grad_hessian.hpp b/stan/math/mix/functor/grad_hessian.hpp index f9ff5fa27f0..9fcf47ba1f6 100644 --- a/stan/math/mix/functor/grad_hessian.hpp +++ b/stan/math/mix/functor/grad_hessian.hpp @@ -50,10 +50,8 @@ void grad_hessian( int d = x.size(); H.resize(d, d); grad_H.resize(d, Matrix(d, d)); - //for (int i = 0; i < d; ++i) { - int i = 0; - //for (int j = i; j < d; ++j) { - int j = 0; + for (int i = 0; i < d; ++i) { + for (int j = i; j < d; ++j) { // Run nested autodiff in this scope nested_rev_autodiff nested; @@ -70,8 +68,8 @@ void grad_hessian( grad_H[i](j, k) = x_ffvar(k).val_.val_.adj(); grad_H[j](i, k) = grad_H[i](j, k); } - //} - //} + } + } } } // namespace math diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 7f915e8816b..d9bf4dc873a 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -7,12 +7,14 @@ #include #include #include +#include #include #include #include namespace stan { namespace math { +namespace internal { /** * Wrapper function for calculating gradients for the generalized * hypergeometric function. The function always returns a tuple with @@ -31,24 +33,21 @@ namespace math { * @param[in] inner_steps Maximum number of iterations for inner_steps * @return Tuple of gradients */ -template -auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, +template +auto grad_pFq_impl(const TpFq& pfq_val, const Ta& a_val, const Tb& b_val, const Tz& z_val, double precision = 1e-16, int max_steps = 1e8) { using std::max; - using T_partials = partials_return_t; - using Ta_partials = partials_return_t; - using Tb_partials = partials_return_t; - using Tz_partials = partials_return_t; + using T_partials = return_type_t; + using Ta_partials = return_type_t; + using Tb_partials = return_type_t; + using Tz_partials = Tz; using T_ArrayPartials = promote_scalar_t>; using IntArray = Eigen::Array; using Ta_PartialsArray = promote_scalar_t>; using Tb_PartialsArray = promote_scalar_t>; - Ta_PartialsArray a_val = as_value_column_array_or_scalar(a); - Tb_PartialsArray b_val = as_value_column_array_or_scalar(b); - Tz_partials z_val = value_of(z); - IntArray a_signs = sign(value_of_rec(a_val)); IntArray b_signs = sign(value_of_rec(b_val)); int z_sign = sign(value_of_rec(z_val)); @@ -62,10 +61,10 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, Tz_partials log_z_pow_k = 0.0; Tz_partials log_z_k = log(abs(z_val)); - Ta_PartialsArray log_phammer_a_k = Ta_PartialsArray::Zero(a.size()); - Tb_PartialsArray log_phammer_b_k = Tb_PartialsArray::Zero(b.size()); - IntArray phammer_a_k_signs = IntArray::Ones(a.size()); - IntArray phammer_b_k_signs = IntArray::Ones(b.size()); + Ta_PartialsArray log_phammer_a_k = Ta_PartialsArray::Zero(a_val.size()); + Tb_PartialsArray log_phammer_b_k = Tb_PartialsArray::Zero(b_val.size()); + IntArray phammer_a_k_signs = IntArray::Ones(a_val.size()); + IntArray phammer_b_k_signs = IntArray::Ones(b_val.size()); Ta_PartialsArray digamma_a = inv(a_k); Tb_PartialsArray digamma_b = inv(b_k); @@ -76,15 +75,15 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, IntArray digamma_a_k_signs = a_k_signs; IntArray digamma_b_k_signs = b_k_signs; - T_ArrayPartials a_log_infsum = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); - T_ArrayPartials b_log_infsum = T_ArrayPartials::Constant(b.size(), NEGATIVE_INFTY); - IntArray a_infsum_sign = IntArray::Ones(a.size()); - IntArray b_infsum_sign = IntArray::Ones(b.size()); + T_ArrayPartials a_log_infsum = T_ArrayPartials::Constant(a_val.size(), NEGATIVE_INFTY); + T_ArrayPartials b_log_infsum = T_ArrayPartials::Constant(b_val.size(), NEGATIVE_INFTY); + IntArray a_infsum_sign = IntArray::Ones(a_val.size()); + IntArray b_infsum_sign = IntArray::Ones(b_val.size()); - T_ArrayPartials a_grad = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); - T_ArrayPartials b_grad = T_ArrayPartials::Constant(a.size(), NEGATIVE_INFTY); - IntArray a_grad_signs = IntArray::Ones(a.size()); - IntArray b_grad_signs = IntArray::Ones(a.size()); + T_ArrayPartials a_grad = T_ArrayPartials::Constant(a_val.size(), NEGATIVE_INFTY); + T_ArrayPartials b_grad = T_ArrayPartials::Constant(a_val.size(), NEGATIVE_INFTY); + IntArray a_grad_signs = IntArray::Ones(a_val.size()); + IntArray b_grad_signs = IntArray::Ones(a_val.size()); int k = 0; T_partials curr_log_prec = 1.0; @@ -94,22 +93,22 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, T_partials log_base = (sum(log_phammer_a_k) + log_z_pow_k) - (log_factorial_k + sum(log_phammer_b_k)); - if (!is_constant_all::value) { + if (CalcA) { a_grad = log_digamma_a_k + log_base; a_grad_signs = base_sign * digamma_a_k_signs; - for (int i = 0; i < a.size(); i++) { + for (int i = 0; i < a_val.size(); i++) { std::forward_as_tuple(a_log_infsum(i), a_infsum_sign(i)) = log_sum_exp_signed(a_log_infsum(i), a_infsum_sign(i), a_grad(i), a_grad_signs(i)); } } - if (!is_constant_all::value) { + if (CalcB) { b_grad = log_digamma_b_k + log_base; b_grad_signs = base_sign * digamma_b_k_signs; - for (int i = 0; i < b.size(); i++) { + for (int i = 0; i < b_val.size(); i++) { std::forward_as_tuple(b_log_infsum(i), b_infsum_sign(i)) = log_sum_exp_signed(b_log_infsum(i), b_infsum_sign(i), b_grad(i), b_grad_signs(i)); @@ -140,26 +139,39 @@ auto grad_pFq(const Ta& a, const Tb& b, const Tz& z, k += 1; log_factorial_k += log(k); } - - std::tuple>, - promote_scalar_t>, + std::tuple, + Eigen::Matrix, T_partials> ret_tuple; - if (!is_constant_all::value) { - T_partials pfq_val = hypergeometric_pFq(a_val.matrix(), b_val.matrix(), z_val); - if (!is_constant_all::value) { - std::get<0>(ret_tuple) = a_infsum_sign * exp(a_log_infsum) - digamma_a * pfq_val; - } - if (!is_constant_all::value) { - std::get<1>(ret_tuple) = digamma_b * pfq_val - b_infsum_sign * exp(b_log_infsum); - } + if (CalcA) { + std::get<0>(ret_tuple) = a_infsum_sign * exp(a_log_infsum) - digamma_a * pfq_val; + } + if (CalcB) { + std::get<1>(ret_tuple) = digamma_b * pfq_val - b_infsum_sign * exp(b_log_infsum); } - if (!is_constant_all::value) { + if (CalcZ) { T_partials pfq_p1_val = hypergeometric_pFq((a_val + 1).matrix(), (b_val + 1).matrix(), z_val); std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; } return ret_tuple; } +} + + +template +auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, + double precision = 1e-16, + int max_steps = 1e8) { + return internal::grad_pFq_impl< + !is_constant_all::value, + !is_constant_all::value, + !is_constant_all::value>( + pfq_val, + to_ref(as_value_column_array_or_scalar(a)), + to_ref(as_value_column_array_or_scalar(b)), + value_of(z) + ); +} } // namespace math } // namespace stan diff --git a/stan/math/rev/fun/hypergeometric_pFq.hpp b/stan/math/rev/fun/hypergeometric_pFq.hpp index 360deabab89..48ea68a6abd 100644 --- a/stan/math/rev/fun/hypergeometric_pFq.hpp +++ b/stan/math/rev/fun/hypergeometric_pFq.hpp @@ -27,21 +27,21 @@ template arena_a = a; arena_t arena_b = b; - return make_callback_var( - hypergeometric_pFq(value_of(arena_a), value_of(arena_b), value_of(z)), - [arena_a, arena_b, z](auto& vi) mutable { - auto grad_tuple = grad_pFq(arena_a, arena_b, z); + auto pfq_val = hypergeometric_pFq(value_of(a), value_of(b), value_of(z)); + return make_callback_var(pfq_val, + [arena_a, arena_b, z, pfq_val](auto& vi) mutable { + auto grad_tuple = grad_pFq(pfq_val, arena_a, arena_b, z); if (!is_constant::value) { - //forward_as>(arena_a).adj() - // += vi.adj() * std::get<0>(grad_tuple); + forward_as>(arena_a).adj() + += vi.adj() * std::get<0>(grad_tuple); } if (!is_constant::value) { - //forward_as>(arena_b).adj() - // += vi.adj() * std::get<1>(grad_tuple); + forward_as>(arena_b).adj() + += vi.adj() * std::get<1>(grad_tuple); } if (!is_constant::value) { - //forward_as>(z).adj() - // += vi.adj() * std::get<2>(grad_tuple); + forward_as>(z).adj() + += vi.adj() * std::get<2>(grad_tuple); } }); } diff --git a/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp b/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp index afec57accd6..52d1953d40f 100644 --- a/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp +++ b/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp @@ -14,18 +14,7 @@ TEST(mathMixScalFun, hyper_pfq) { in2 << 6, 3; double z = 4; - Eigen::VectorXd x = stan::math::to_vector(stan::test::serialize(in1, in2, z)); - auto serial_functor = [&](const auto& v) { - auto v_deserializer = stan::test::to_deserializer(v); - return f(v_deserializer.read(in1), v_deserializer.read(in2), v_deserializer.read(z)); - }; - double fx_ad; - Eigen::MatrixXd H_ad; - std::vector grad_H_ad; - stan::math::grad_hessian(serial_functor, x, fx_ad, H_ad, grad_H_ad); - //std::cout << H_ad << "\n" << std::endl; - std::cout << grad_H_ad[0] << std::endl; - //stan::test::expect_ad(f, in1, in2, z); + stan::test::expect_ad(f, in1, in2, z); } /* TEST(mixScalFun, grad_2F2_ffv) { diff --git a/test/unit/math/test_ad.hpp b/test/unit/math/test_ad.hpp index ff7504d6a41..78eae739334 100644 --- a/test/unit/math/test_ad.hpp +++ b/test/unit/math/test_ad.hpp @@ -306,11 +306,9 @@ void test_grad_hessian(const ad_tolerances& tols, const F& f, expect_near_rel("grad_hessian() Hessian", H_fd, H_ad, tols.grad_hessian_hessian_); EXPECT_EQ(x.size(), grad_H_fd.size()); - for (size_t i = 0; i < grad_H_fd.size(); ++i) { - std::cout << i << "\n\n" << grad_H_ad[i] << "\n" << std::endl; - //expect_near_rel("grad_hessian() grad Hessian", grad_H_fd[i], grad_H_ad[i], - // tols.grad_hessian_grad_hessian_); - } + for (size_t i = 0; i < grad_H_fd.size(); ++i) + expect_near_rel("grad_hessian() grad Hessian", grad_H_fd[i], grad_H_ad[i], + tols.grad_hessian_grad_hessian_); } #endif From 67fdd57f04d5a5a4adb015db2a196441d09efc36 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 6 Oct 2023 19:54:03 +0300 Subject: [PATCH 10/36] Updating test --- stan/math/prim/fun/grad_pFq.hpp | 7 +- .../math/mix/fun/hypergeometric_pFq_test.cpp | 120 ++++++++++-------- test/unit/math/rev/fun/grad_pFq_test.cpp | 19 +-- 3 files changed, 84 insertions(+), 62 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index d9bf4dc873a..b0570450787 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -98,9 +98,9 @@ auto grad_pFq_impl(const TpFq& pfq_val, const Ta& a_val, const Tb& b_val, const a_grad_signs = base_sign * digamma_a_k_signs; for (int i = 0; i < a_val.size(); i++) { - std::forward_as_tuple(a_log_infsum(i), a_infsum_sign(i)) - = log_sum_exp_signed(a_log_infsum(i), a_infsum_sign(i), - a_grad(i), a_grad_signs(i)); + std::forward_as_tuple(a_log_infsum.coeffRef(i), a_infsum_sign.coeffRef(i)) + = log_sum_exp_signed(a_log_infsum.coeffRef(i), a_infsum_sign.coeffRef(i), + a_grad.coeffRef(i), a_grad_signs.coeffRef(i)); } } @@ -172,7 +172,6 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, value_of(z) ); } - } // namespace math } // namespace stan #endif diff --git a/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp b/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp index 52d1953d40f..b7fb374a3dc 100644 --- a/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp +++ b/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp @@ -1,8 +1,7 @@ #include #include -TEST(mathMixScalFun, hyper_pfq) { - using stan::math::hypergeometric_pFq; +TEST(mathMixScalFun, hyper_2f2) { auto f = [](const auto& a, const auto& b, const auto& z) { using stan::math::hypergeometric_pFq; return hypergeometric_pFq(a, b, z); @@ -16,52 +15,73 @@ TEST(mathMixScalFun, hyper_pfq) { stan::test::expect_ad(f, in1, in2, z); } -/* -TEST(mixScalFun, grad_2F2_ffv) { - using stan::math::fvar; - using stan::math::hypergeometric_pFq; - using stan::math::var; - using stan::math::vector_d; - using stan::math::vector_ffv; - - vector_ffv ffv_a(2); - vector_d d_a(2); - ffv_a.val().val() << 4, 2; - d_a << 4, 2; - ffv_a.val().d() << 1, 1; - - vector_ffv ffv_b(2); - vector_d d_b(2); - ffv_b.val().val() << 6, 3; - d_b << 6, 3; - ffv_b.val().d() << 1, 1; - - fvar> ffv_z; - ffv_z.val_.val_ = 4; - ffv_z.val_.d_ = 1; - double d_z = 4; - - double a_adj = 3.924636646666071 + 6.897245961898751; - double b_adj = -2.775051002566842 - 4.980095849781222; - double z_adj = 4.916522138006060; - - // fvar, fvar, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(ffv_a, ffv_b, ffv_z).val_.d_.val(), - a_adj + b_adj + z_adj); - // fvar, fvar, double - EXPECT_FLOAT_EQ(hypergeometric_pFq(ffv_a, ffv_b, d_z).val_.d_.val(), - a_adj + b_adj); - // fvar, double, double - EXPECT_FLOAT_EQ(hypergeometric_pFq(ffv_a, d_b, d_z).val_.d_.val(), a_adj); - // fvar, double, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(ffv_a, d_b, ffv_z).val_.d_.val(), - a_adj + z_adj); - // double, fvar, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(d_a, ffv_b, ffv_z).val_.d_.val(), - b_adj + z_adj); - // double, fvar, double - EXPECT_FLOAT_EQ(hypergeometric_pFq(d_a, ffv_b, d_z).val_.d_.val(), b_adj); - // double, double, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(d_a, d_b, ffv_z).val_.d_.val(), z_adj); + +TEST(mathMixScalFun, hyper_2f3) { + auto f = [](const auto& a, const auto& b, const auto& z) { + using stan::math::hypergeometric_pFq; + return hypergeometric_pFq(a, b, z); + }; + + Eigen::VectorXd in1(2); + in1 << 2, 3; + Eigen::VectorXd in2(3); + in2 << 2, 4, 5; + double z = 1; + + stan::test::expect_ad(f, in1, in2, z); +} + +TEST(mathMixScalFun, hyper_4f3) { + auto f = [](const auto& a, const auto& b, const auto& z) { + using stan::math::hypergeometric_pFq; + return hypergeometric_pFq(a, b, z); + }; + + Eigen::VectorXd in1(4); + in1 << 1, 2, 3, 4; + Eigen::VectorXd in2(3); + in2 << 5, 6, 7; + double z = 0.8; + + stan::test::expect_ad(f, in1, in2, z); +} + +TEST(mathMixScalFun, hyper_2f1) { + auto f = [](const auto& a, const auto& b, const auto& z) { + using stan::math::hypergeometric_pFq; + return hypergeometric_pFq(a, b, z); + }; + + Eigen::VectorXd in1(2); + in1 << 1, 1; + Eigen::VectorXd in2(1); + in2 << 1; + double z = 0.6; + + stan::test::expect_ad(f, in1, in2, z); + + in1 << 1, 31; + in2 << 41; + z = 0.6; + stan::test::expect_ad(f, in1, in2, z); + + in1 << 1, -2.1; + in2 << 41; + z = 0.6; + stan::test::expect_ad(f, in1, in2, z); + + in1 << 1, -0.5; + in2 << 10; + z = 0.3; + //stan::test::expect_ad(f, in1, in2, z); + + in1 << -0.5, -4.5; + in2 << 1; + z = 0.3; + //stan::test::expect_ad(f, in1, in2, z); + + in1 << -0.5, -4.5; + in2 << -3.2; + z = 0.9; + //stan::test::expect_ad(f, in1, in2, z); } - */ diff --git a/test/unit/math/rev/fun/grad_pFq_test.cpp b/test/unit/math/rev/fun/grad_pFq_test.cpp index b6fe6e656b6..43921c4a6fd 100644 --- a/test/unit/math/rev/fun/grad_pFq_test.cpp +++ b/test/unit/math/rev/fun/grad_pFq_test.cpp @@ -4,6 +4,7 @@ TEST(RevMath, grad_2F2) { using stan::math::grad_pFq; + using stan::math::hypergeometric_pFq; using stan::math::var; using stan::math::vector_d; using stan::math::vector_v; @@ -21,7 +22,8 @@ TEST(RevMath, grad_2F2) { var z_v = 4; double z_d = 4; - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_v, b_v, z_v); EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); @@ -31,34 +33,34 @@ TEST(RevMath, grad_2F2) { EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); - grad_tuple = grad_pFq(a_v, b_d, z_d); + grad_tuple = grad_pFq(pfq_val, a_v, b_d, z_d); EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); - grad_tuple = grad_pFq(a_d, b_v, z_d); + grad_tuple = grad_pFq(pfq_val, a_d, b_v, z_d); EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); - grad_tuple = grad_pFq(a_d, b_d, z_v); + grad_tuple = grad_pFq(pfq_val, pfq_val, a_d, b_d, z_v); EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); - grad_tuple = grad_pFq(a_v, b_v, z_d); + grad_tuple = grad_pFq(pfq_val, a_v, b_v, z_d); EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); - grad_tuple = grad_pFq(a_v, b_d, z_v); + grad_tuple = grad_pFq(pfq_val, a_v, b_d, z_v); EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); - grad_tuple = grad_pFq(a_d, b_v, z_v); + grad_tuple = grad_pFq(pfq_val, a_d, b_v, z_v); EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); @@ -75,7 +77,8 @@ TEST(RevMath, grad_2F3) { b_v << 2, 4, 5; var z_v = 1; - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto pfq_val = hypergeometric_pFq(a_v.val(), b_v.val(), z_v.val()); + auto grad_tuple = grad_pFq(pfq_val, a_v, b_v, z_v); EXPECT_FLOAT_EQ(0.08377717301140296, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(0.05615450733193106, std::get<0>(grad_tuple)[1]); From d52db964711da2071a918b52b12ac0908f21715e Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 7 Oct 2023 17:13:48 +0300 Subject: [PATCH 11/36] Workind ad tests --- stan/math/prim/fun/grad_pFq.hpp | 127 +++++++----------- .../math/mix/fun/hypergeometric_pFq_test.cpp | 45 ++++++- 2 files changed, 87 insertions(+), 85 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index b0570450787..4b2f2203c15 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -15,6 +15,11 @@ namespace stan { namespace math { namespace internal { + template + inline auto binarysign(const T& x) { + return select(x == 0.0, 1, sign(value_of_rec(x))); + } + /** * Wrapper function for calculating gradients for the generalized * hypergeometric function. The function always returns a tuple with @@ -36,121 +41,79 @@ namespace internal { template auto grad_pFq_impl(const TpFq& pfq_val, const Ta& a_val, const Tb& b_val, const Tz& z_val, - double precision = 1e-16, - int max_steps = 1e8) { + double precision = 1e-16, + int max_steps = 1e8) { using std::max; using T_partials = return_type_t; using Ta_partials = return_type_t; using Tb_partials = return_type_t; - using Tz_partials = Tz; + using Tz_partials = return_type_t; using T_ArrayPartials = promote_scalar_t>; using IntArray = Eigen::Array; using Ta_PartialsArray = promote_scalar_t>; using Tb_PartialsArray = promote_scalar_t>; - IntArray a_signs = sign(value_of_rec(a_val)); - IntArray b_signs = sign(value_of_rec(b_val)); - int z_sign = sign(value_of_rec(z_val)); - - IntArray a_k_signs = a_signs; - IntArray b_k_signs = b_signs; - int z_k_sign = 1; + int z_sign = binarysign(z_val); Ta_PartialsArray a_k = a_val; Tb_PartialsArray b_k = b_val; - Tz_partials log_z_pow_k = 0.0; - Tz_partials log_z_k = log(abs(z_val)); - - Ta_PartialsArray log_phammer_a_k = Ta_PartialsArray::Zero(a_val.size()); - Tb_PartialsArray log_phammer_b_k = Tb_PartialsArray::Zero(b_val.size()); - IntArray phammer_a_k_signs = IntArray::Ones(a_val.size()); - IntArray phammer_b_k_signs = IntArray::Ones(b_val.size()); + Tz_partials log_z = log(abs(z_val)); Ta_PartialsArray digamma_a = inv(a_k); Tb_PartialsArray digamma_b = inv(b_k); - Ta_PartialsArray digamma_a_k = digamma_a; - Tb_PartialsArray digamma_b_k = digamma_b; - Ta_PartialsArray log_digamma_a_k = log(abs(digamma_a)); - Tb_PartialsArray log_digamma_b_k = log(abs(digamma_b)); - IntArray digamma_a_k_signs = a_k_signs; - IntArray digamma_b_k_signs = b_k_signs; - - T_ArrayPartials a_log_infsum = T_ArrayPartials::Constant(a_val.size(), NEGATIVE_INFTY); - T_ArrayPartials b_log_infsum = T_ArrayPartials::Constant(b_val.size(), NEGATIVE_INFTY); - IntArray a_infsum_sign = IntArray::Ones(a_val.size()); - IntArray b_infsum_sign = IntArray::Ones(b_val.size()); - - T_ArrayPartials a_grad = T_ArrayPartials::Constant(a_val.size(), NEGATIVE_INFTY); - T_ArrayPartials b_grad = T_ArrayPartials::Constant(a_val.size(), NEGATIVE_INFTY); - IntArray a_grad_signs = IntArray::Ones(a_val.size()); - IntArray b_grad_signs = IntArray::Ones(a_val.size()); - int k = 0; - T_partials curr_log_prec = 1.0; - double log_factorial_k = 0; - while (k <= max_steps && curr_log_prec > log(precision)) { - int base_sign = z_k_sign * phammer_a_k_signs.prod() * phammer_b_k_signs.prod(); - T_partials log_base = (sum(log_phammer_a_k) + log_z_pow_k) - - (log_factorial_k + sum(log_phammer_b_k)); + std::tuple, Eigen::Matrix, + T_partials> ret_tuple; + std::get<0>(ret_tuple).setConstant(a_val.size(), 0.0); + std::get<1>(ret_tuple).setConstant(b_val.size(), 0.0); + std::get<2>(ret_tuple) = 0.0; + T_ArrayPartials a_grad(a_val.size()); + T_ArrayPartials b_grad(b_val.size()); + + int k = 0; + int base_sign = 1; + T_partials curr_log_prec = NEGATIVE_INFTY; + T_partials log_base = 0; + while ((k < 10 || curr_log_prec > log(precision)) && (k <= max_steps)) { + curr_log_prec = NEGATIVE_INFTY; if (CalcA) { - a_grad = log_digamma_a_k + log_base; - a_grad_signs = base_sign * digamma_a_k_signs; - - for (int i = 0; i < a_val.size(); i++) { - std::forward_as_tuple(a_log_infsum.coeffRef(i), a_infsum_sign.coeffRef(i)) - = log_sum_exp_signed(a_log_infsum.coeffRef(i), a_infsum_sign.coeffRef(i), - a_grad.coeffRef(i), a_grad_signs.coeffRef(i)); - } - } + a_grad = log(select(digamma_a == 0.0, 1.0, abs(digamma_a))) + log_base; + std::get<0>(ret_tuple).array() + += exp(a_grad) * base_sign * binarysign(digamma_a); - if (CalcB) { - b_grad = log_digamma_b_k + log_base; - b_grad_signs = base_sign * digamma_b_k_signs; - - for (int i = 0; i < b_val.size(); i++) { - std::forward_as_tuple(b_log_infsum(i), b_infsum_sign(i)) - = log_sum_exp_signed(b_log_infsum(i), b_infsum_sign(i), - b_grad(i), b_grad_signs(i)); - } + curr_log_prec = max(curr_log_prec, a_grad.maxCoeff()); + digamma_a += select(a_k == 0.0, 0.0, inv(a_k)); } - curr_log_prec = max(a_grad.maxCoeff(), b_grad.maxCoeff()); - log_phammer_a_k += log(abs(a_k)); - log_phammer_b_k += log(abs(b_k)); - phammer_a_k_signs *= a_k_signs; - phammer_b_k_signs *= b_k_signs; + if (CalcB) { + b_grad = log(select(digamma_b == 0.0, 1.0, abs(digamma_b))) + log_base; + std::get<1>(ret_tuple).array() + -= exp(b_grad) * base_sign * binarysign(digamma_b); - log_z_pow_k += log_z_k; - z_k_sign *= z_sign; + curr_log_prec = max(curr_log_prec, b_grad.maxCoeff()); + digamma_b += select(b_k == 0.0, 0.0, inv(b_k)); + } - digamma_a_k += select(a_k == 0.0, 0.0, inv(a_k)); - digamma_b_k += select(b_k == 0.0, 0.0, inv(b_k)); - log_digamma_a_k = log(abs(digamma_a_k)); - log_digamma_b_k = log(abs(digamma_b_k)); - digamma_a_k_signs = sign(value_of_rec(digamma_a_k)); - digamma_b_k_signs = sign(value_of_rec(digamma_b_k)); + log_base += sum(log(select(a_k == 0.0, 0.0, abs(a_k)))) + log_z; + log_base -= sum(log(select(b_k == 0.0, 0.0, abs(b_k)))) + log1p(k); + base_sign *= z_sign * binarysign(a_k).prod() * binarysign(b_k).prod(); a_k += 1.0; b_k += 1.0; - a_k_signs = sign(value_of_rec(a_k)); - b_k_signs = sign(value_of_rec(b_k)); k += 1; - log_factorial_k += log(k); } - std::tuple, - Eigen::Matrix, - T_partials> ret_tuple; if (CalcA) { - std::get<0>(ret_tuple) = a_infsum_sign * exp(a_log_infsum) - digamma_a * pfq_val; + std::get<0>(ret_tuple).array() -= pfq_val / a_val; } if (CalcB) { - std::get<1>(ret_tuple) = digamma_b * pfq_val - b_infsum_sign * exp(b_log_infsum); + std::get<1>(ret_tuple).array() += pfq_val / b_val; } if (CalcZ) { - T_partials pfq_p1_val = hypergeometric_pFq((a_val + 1).matrix(), (b_val + 1).matrix(), z_val); + T_partials pfq_p1_val = hypergeometric_pFq((a_val + 1).matrix(), + (b_val + 1).matrix(), z_val); std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; } return ret_tuple; @@ -169,7 +132,9 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, pfq_val, to_ref(as_value_column_array_or_scalar(a)), to_ref(as_value_column_array_or_scalar(b)), - value_of(z) + value_of(z), + precision, + max_steps ); } } // namespace math diff --git a/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp b/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp index b7fb374a3dc..ba498378459 100644 --- a/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp +++ b/test/unit/math/mix/fun/hypergeometric_pFq_test.cpp @@ -73,15 +73,52 @@ TEST(mathMixScalFun, hyper_2f1) { in1 << 1, -0.5; in2 << 10; z = 0.3; - //stan::test::expect_ad(f, in1, in2, z); + stan::test::expect_ad(f, in1, in2, z); + + in1 << 1, -0.5; + in2 << 10.6; + z = 0.3; + stan::test::expect_ad(f, in1, in2, z); in1 << -0.5, -4.5; - in2 << 1; + in2 << 11; z = 0.3; - //stan::test::expect_ad(f, in1, in2, z); + stan::test::expect_ad(f, in1, in2, z); in1 << -0.5, -4.5; in2 << -3.2; z = 0.9; - //stan::test::expect_ad(f, in1, in2, z); + stan::test::expect_ad(f, in1, in2, z); + + in1 << 2, 1; + in2 << 2; + z = 0.4; + stan::test::expect_ad(f, in1, in2, z); + + in1 << 3.70975, 1.0; + in2 << 2.70975; + z = -0.2; + stan::test::expect_ad(f, in1, in2, z); +} + +TEST(mathMixScalFun, hyper_3f2) { + using stan::math::var; + + auto f = [](const auto& a, const auto& b, const auto& z) { + using stan::math::hypergeometric_pFq; + return hypergeometric_pFq(a, b, z); + }; + + Eigen::VectorXd in1(3); + in1 << 1.0, 1.0, 1.0; + Eigen::VectorXd in2(2); + in2 << 1.0, 1.0; + double z = 0.6; + + stan::test::expect_ad(f, in1, in2, z); + + in1 << 1.0, -0.5, -2.5; + in2 << 10.0, 1.0; + z = 0.3; + stan::test::expect_ad(f, in1, in2, z); } From 01988f1116753d39069366bfee3c5061babae722 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 7 Oct 2023 18:03:47 +0300 Subject: [PATCH 12/36] Simplified and tidied --- stan/math/fwd/fun/hypergeometric_pFq.hpp | 16 +- stan/math/prim/fun/grad_pFq.hpp | 95 ++--- stan/math/rev/fun/hypergeometric_pFq.hpp | 16 +- .../math/{rev => prim}/fun/grad_pFq_test.cpp | 346 +++++++++--------- 4 files changed, 235 insertions(+), 238 deletions(-) rename test/unit/math/{rev => prim}/fun/grad_pFq_test.cpp (58%) diff --git a/stan/math/fwd/fun/hypergeometric_pFq.hpp b/stan/math/fwd/fun/hypergeometric_pFq.hpp index 75e0e229768..f63164a9d13 100644 --- a/stan/math/fwd/fun/hypergeometric_pFq.hpp +++ b/stan/math/fwd/fun/hypergeometric_pFq.hpp @@ -22,6 +22,9 @@ namespace math { * @return Generalized hypergeometric function */ template ::value, + bool GradB = !is_constant::value, + bool GradZ = !is_constant::value, require_all_matrix_t* = nullptr, require_return_type_t* = nullptr> inline return_type_t hypergeometric_pFq(const Ta& a, const Tb& b, @@ -29,20 +32,23 @@ inline return_type_t hypergeometric_pFq(const Ta& a, const Tb& b, using fvar_t = return_type_t; ref_type_t a_ref = a; ref_type_t b_ref = b; - auto pfq_val = hypergeometric_pFq(value_of(a_ref), value_of(b_ref), value_of(z)); - auto grad_tuple = grad_pFq(pfq_val, a_ref, b_ref, z); + auto&& a_val = value_of(a_ref); + auto&& b_val = value_of(b_ref); + auto&& z_val = value_of(z); + auto pfq_val = hypergeometric_pFq(a_val, b_val, z_val); + auto grad_tuple = grad_pFq(pfq_val, a_val, b_val, z_val); typename fvar_t::Scalar grad = 0; - if (!is_constant::value) { + if (GradA) { grad += dot_product(forward_as>(a_ref).d(), std::get<0>(grad_tuple)); } - if (!is_constant::value) { + if (GradB) { grad += dot_product(forward_as>(b_ref).d(), std::get<1>(grad_tuple)); } - if (!is_constant::value) { + if (GradZ) { grad += forward_as>(z).d_ * std::get<2>(grad_tuple); } diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 4b2f2203c15..0b31a970e16 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -19,6 +19,7 @@ namespace internal { inline auto binarysign(const T& x) { return select(x == 0.0, 1, sign(value_of_rec(x))); } +} /** * Wrapper function for calculating gradients for the generalized @@ -34,53 +35,49 @@ namespace internal { * @param[in] b Vector of 'b' arguments to function * @param[in] z Scalar z argument * @param[in] precision Convergence criteria for infinite sum - * @param[in] outer_steps Maximum number of iterations for outer_steps - * @param[in] inner_steps Maximum number of iterations for inner_steps + * @param[in] outer_steps Maximum number of iterations * @return Tuple of gradients */ -template -auto grad_pFq_impl(const TpFq& pfq_val, const Ta& a_val, const Tb& b_val, const Tz& z_val, - double precision = 1e-16, - int max_steps = 1e8) { +auto grad_pFq(const TpFq& pfq_val, const Ta& a_val, const Tb& b_val, + const Tz& z_val, double precision = 1e-14, int max_steps = 1e6) { using std::max; - using T_partials = return_type_t; - using Ta_partials = return_type_t; - using Tb_partials = return_type_t; - using Tz_partials = return_type_t; - using T_ArrayPartials = promote_scalar_t>; - using IntArray = Eigen::Array; - using Ta_PartialsArray = promote_scalar_t>; - using Tb_PartialsArray = promote_scalar_t>; - - int z_sign = binarysign(z_val); - - Ta_PartialsArray a_k = a_val; - Tb_PartialsArray b_k = b_val; - Tz_partials log_z = log(abs(z_val)); - - Ta_PartialsArray digamma_a = inv(a_k); - Tb_PartialsArray digamma_b = inv(b_k); - - std::tuple, Eigen::Matrix, - T_partials> ret_tuple; + using T_Rtn = return_type_t; + using Ta_Array = Eigen::Array, -1, 1>; + using Tb_Array = Eigen::Array, -1, 1>; + using T_RtnVector = Eigen::Matrix; + + Ta_Array a_array = as_column_vector_or_scalar(a_val).array(); + Ta_Array a_k = a_array; + Tb_Array b_array = as_column_vector_or_scalar(b_val).array(); + Tb_Array b_k = b_array; + Tz log_z = log(abs(z_val)); + int z_sign = internal::binarysign(z_val); + + Ta_Array digamma_a = inv(a_k); + Tb_Array digamma_b = inv(b_k); + + std::tuple< + promote_scalar_t>, + promote_scalar_t>, T_Rtn> ret_tuple; std::get<0>(ret_tuple).setConstant(a_val.size(), 0.0); std::get<1>(ret_tuple).setConstant(b_val.size(), 0.0); std::get<2>(ret_tuple) = 0.0; - T_ArrayPartials a_grad(a_val.size()); - T_ArrayPartials b_grad(b_val.size()); + Eigen::Array a_grad(a_val.size()); + Eigen::Array b_grad(b_val.size()); int k = 0; int base_sign = 1; - T_partials curr_log_prec = NEGATIVE_INFTY; - T_partials log_base = 0; + T_Rtn curr_log_prec = NEGATIVE_INFTY; + T_Rtn log_base = 0; while ((k < 10 || curr_log_prec > log(precision)) && (k <= max_steps)) { curr_log_prec = NEGATIVE_INFTY; if (CalcA) { a_grad = log(select(digamma_a == 0.0, 1.0, abs(digamma_a))) + log_base; std::get<0>(ret_tuple).array() - += exp(a_grad) * base_sign * binarysign(digamma_a); + += exp(a_grad) * base_sign * internal::binarysign(digamma_a); curr_log_prec = max(curr_log_prec, a_grad.maxCoeff()); digamma_a += select(a_k == 0.0, 0.0, inv(a_k)); @@ -89,54 +86,36 @@ auto grad_pFq_impl(const TpFq& pfq_val, const Ta& a_val, const Tb& b_val, const if (CalcB) { b_grad = log(select(digamma_b == 0.0, 1.0, abs(digamma_b))) + log_base; std::get<1>(ret_tuple).array() - -= exp(b_grad) * base_sign * binarysign(digamma_b); + -= exp(b_grad) * base_sign * internal::binarysign(digamma_b); curr_log_prec = max(curr_log_prec, b_grad.maxCoeff()); digamma_b += select(b_k == 0.0, 0.0, inv(b_k)); } - log_base += sum(log(select(a_k == 0.0, 0.0, abs(a_k)))) + log_z; - log_base -= sum(log(select(b_k == 0.0, 0.0, abs(b_k)))) + log1p(k); - base_sign *= z_sign * binarysign(a_k).prod() * binarysign(b_k).prod(); + log_base += (sum(log(select(a_k == 0.0, 0.0, abs(a_k)))) + log_z) + - (sum(log(select(b_k == 0.0, 0.0, abs(b_k)))) + log1p(k)); + base_sign *= z_sign * internal::binarysign(a_k).prod() + * internal::binarysign(b_k).prod(); a_k += 1.0; b_k += 1.0; - k += 1; } if (CalcA) { - std::get<0>(ret_tuple).array() -= pfq_val / a_val; + std::get<0>(ret_tuple).array() -= pfq_val / a_array; } if (CalcB) { - std::get<1>(ret_tuple).array() += pfq_val / b_val; + std::get<1>(ret_tuple).array() += pfq_val / b_array; } if (CalcZ) { - T_partials pfq_p1_val = hypergeometric_pFq((a_val + 1).matrix(), - (b_val + 1).matrix(), z_val); + T_Rtn pfq_p1_val = hypergeometric_pFq((a_array + 1).matrix(), + (b_array + 1).matrix(), z_val); std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; } return ret_tuple; } -} - -template -auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, - double precision = 1e-16, - int max_steps = 1e8) { - return internal::grad_pFq_impl< - !is_constant_all::value, - !is_constant_all::value, - !is_constant_all::value>( - pfq_val, - to_ref(as_value_column_array_or_scalar(a)), - to_ref(as_value_column_array_or_scalar(b)), - value_of(z), - precision, - max_steps - ); -} } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/fun/hypergeometric_pFq.hpp b/stan/math/rev/fun/hypergeometric_pFq.hpp index 48ea68a6abd..92dc7baff71 100644 --- a/stan/math/rev/fun/hypergeometric_pFq.hpp +++ b/stan/math/rev/fun/hypergeometric_pFq.hpp @@ -22,24 +22,28 @@ namespace math { * @return Generalized hypergeometric function */ template ::value, + bool GradB = !is_constant::value, + bool GradZ = !is_constant::value, require_all_matrix_t* = nullptr, require_return_type_t* = nullptr> inline var hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { arena_t arena_a = a; arena_t arena_b = b; - auto pfq_val = hypergeometric_pFq(value_of(a), value_of(b), value_of(z)); + auto pfq_val = hypergeometric_pFq(a.val(), b.val(), value_of(z)); + auto grad_tuple + = grad_pFq(pfq_val, a.val(), b.val(), value_of(z)); return make_callback_var(pfq_val, - [arena_a, arena_b, z, pfq_val](auto& vi) mutable { - auto grad_tuple = grad_pFq(pfq_val, arena_a, arena_b, z); - if (!is_constant::value) { + [arena_a, arena_b, z, grad_tuple](auto& vi) mutable { + if (GradA) { forward_as>(arena_a).adj() += vi.adj() * std::get<0>(grad_tuple); } - if (!is_constant::value) { + if (GradB) { forward_as>(arena_b).adj() += vi.adj() * std::get<1>(grad_tuple); } - if (!is_constant::value) { + if (GradZ) { forward_as>(z).adj() += vi.adj() * std::get<2>(grad_tuple); } diff --git a/test/unit/math/rev/fun/grad_pFq_test.cpp b/test/unit/math/prim/fun/grad_pFq_test.cpp similarity index 58% rename from test/unit/math/rev/fun/grad_pFq_test.cpp rename to test/unit/math/prim/fun/grad_pFq_test.cpp index 43921c4a6fd..2f8c0cd5ef1 100644 --- a/test/unit/math/rev/fun/grad_pFq_test.cpp +++ b/test/unit/math/prim/fun/grad_pFq_test.cpp @@ -2,28 +2,21 @@ #include #include -TEST(RevMath, grad_2F2) { +TEST(PrimMath, grad_2F2) { using stan::math::grad_pFq; using stan::math::hypergeometric_pFq; - using stan::math::var; using stan::math::vector_d; - using stan::math::vector_v; - vector_v a_v(2); vector_d a_d(2); - a_v << 4, 2; a_d << 4, 2; - vector_v b_v(2); vector_d b_d(2); - b_v << 6, 3; b_d << 6, 3; - var z_v = 4; double z_d = 4; auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); - auto grad_tuple = grad_pFq(pfq_val, a_v, b_v, z_v); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); @@ -33,52 +26,53 @@ TEST(RevMath, grad_2F2) { EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); - grad_tuple = grad_pFq(pfq_val, a_v, b_d, z_d); + grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); - grad_tuple = grad_pFq(pfq_val, a_d, b_v, z_d); + grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); - grad_tuple = grad_pFq(pfq_val, pfq_val, a_d, b_d, z_v); + grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); - grad_tuple = grad_pFq(pfq_val, a_v, b_v, z_d); + grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); - grad_tuple = grad_pFq(pfq_val, a_v, b_d, z_v); + grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(3.924636646666071, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(6.897245961898751, std::get<0>(grad_tuple)[1]); EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); - grad_tuple = grad_pFq(pfq_val, a_d, b_v, z_v); + grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(-2.775051002566842, std::get<1>(grad_tuple)[0]); EXPECT_FLOAT_EQ(-4.980095849781222, std::get<1>(grad_tuple)[1]); EXPECT_FLOAT_EQ(4.916522138006060, std::get<2>(grad_tuple)); } -TEST(RevMath, grad_2F3) { + +TEST(PrimMath, grad_2F3) { using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << 2, 3; - vector_v b_v(3); - b_v << 2, 4, 5; - var z_v = 1; + vector_d a_d(2); + a_d << 2, 3; + vector_d b_d(3); + b_d << 2, 4, 5; + double z_d = 1; - auto pfq_val = hypergeometric_pFq(a_v.val(), b_v.val(), z_v.val()); - auto grad_tuple = grad_pFq(pfq_val, a_v, b_v, z_v); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(0.08377717301140296, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(0.05615450733193106, std::get<0>(grad_tuple)[1]); @@ -90,18 +84,19 @@ TEST(RevMath, grad_2F3) { EXPECT_FLOAT_EQ(0.1712340452215524, std::get<2>(grad_tuple)); } -TEST(RevMath, grad_4F3) { +TEST(PrimMath, grad_4F3) { using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(4); - a_v << 1, 2, 3, 4; - vector_v b_v(3); - b_v << 5, 6, 7; - var z_v = 1; + vector_d a_d(4); + a_d << 1, 2, 3, 4; + vector_d b_d(3); + b_d << 5, 6, 7; + double z_d = 1; - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(0.1587098625610631, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(0.08249338029396848, std::get<0>(grad_tuple)[1]); @@ -115,20 +110,21 @@ TEST(RevMath, grad_4F3) { EXPECT_FLOAT_EQ(0.1800529055890911, std::get<2>(grad_tuple)); } -TEST(RevMath, grad_2F1_derivs_match) { +TEST(PrimMath, grad_2F1_derivs_match) { using stan::math::grad_2F1; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << 1, 1; - vector_v b_v(1); - b_v << 1; - var z_v = 0.6; + vector_d a_d(2); + a_d << 1, 1; + vector_d b_d(1); + b_d << 1; + double z_d = 0.6; - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_d[0], a_d[1], b_d[0], z_d); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); @@ -136,20 +132,21 @@ TEST(RevMath, grad_2F1_derivs_match) { EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -TEST(RevMath, grad2F1_2) { +TEST(PrimMath, grad2F1_2) { using stan::math::grad_2F1; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << 1, 31; - vector_v b_v(1); - b_v << 41; - var z_v = 1; + vector_d a_d(2); + a_d << 1, 31; + vector_d b_d(1); + b_d << 41; + double z_d = 1; - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_d[0], a_d[1], b_d[0], z_d); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); @@ -157,20 +154,21 @@ TEST(RevMath, grad2F1_2) { EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -TEST(RevMath, grad2F1_3) { +TEST(PrimMath, grad2F1_3) { using stan::math::grad_2F1; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << 1.0, -2.1; - vector_v b_v(1); - b_v << 41.0; - var z_v = 1.0; + vector_d a_d(2); + a_d << 1.0, -2.1; + vector_d b_d(1); + b_d << 41.0; + double z_d = 1.0; - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_d[0], a_d[1], b_d[0], z_d); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); @@ -178,20 +176,21 @@ TEST(RevMath, grad2F1_3) { EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -TEST(RevMath, grad2F1_6) { +TEST(PrimMath, grad2F1_6) { using stan::math::grad_2F1; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << 1, -0.5; - vector_v b_v(1); - b_v << 10.6; - var z_v = 0.3; + vector_d a_d(2); + a_d << 1, -0.5; + vector_d b_d(1); + b_d << 10.6; + double z_d = 0.3; - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_d[0], a_d[1], b_d[0], z_d); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); @@ -199,20 +198,21 @@ TEST(RevMath, grad2F1_6) { EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -TEST(RevMath, grad2F1_7) { +TEST(PrimMath, grad2F1_7) { using stan::math::grad_2F1; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << 1, -0.5; - vector_v b_v(1); - b_v << 10; - var z_v = 0.3; + vector_d a_d(2); + a_d << 1, -0.5; + vector_d b_d(1); + b_d << 10; + double z_d = 0.3; - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_d[0], a_d[1], b_d[0], z_d); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); @@ -220,20 +220,21 @@ TEST(RevMath, grad2F1_7) { EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -TEST(RevMath, grad2F1_8) { +TEST(PrimMath, grad2F1_8) { using stan::math::grad_2F1; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << -0.5, -4.5; - vector_v b_v(1); - b_v << 11; - var z_v = 0.3; + vector_d a_d(2); + a_d << -0.5, -4.5; + vector_d b_d(1); + b_d << 11; + double z_d = 0.3; - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_d[0], a_d[1], b_d[0], z_d); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); @@ -241,20 +242,21 @@ TEST(RevMath, grad2F1_8) { EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -TEST(RevMath, grad2F1_9) { +TEST(PrimMath, grad2F1_9) { using stan::math::grad_2F1; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << -0.5, -4.5; - vector_v b_v(1); - b_v << -3.2; - var z_v = 0.9; + vector_d a_d(2); + a_d << -0.5, -4.5; + vector_d b_d(1); + b_d << -3.2; + double z_d = 0.9; - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_d[0], a_d[1], b_d[0], z_d); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); @@ -262,20 +264,21 @@ TEST(RevMath, grad2F1_9) { EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -TEST(RevMath, grad2F1_10) { +TEST(PrimMath, grad2F1_10) { using stan::math::grad_2F1; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << 2, 1; - vector_v b_v(1); - b_v << 2; - var z_v = 0.4; + vector_d a_d(2); + a_d << 2, 1; + vector_d b_d(1); + b_d << 2; + double z_d = 0.4; - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_d[0], a_d[1], b_d[0], z_d); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); @@ -283,21 +286,22 @@ TEST(RevMath, grad2F1_10) { EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -TEST(RevMath, grad2F1_11) { +TEST(PrimMath, grad2F1_11) { using stan::math::grad_2F1; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << 3.70975, 1; - vector_v b_v(1); - b_v << 2.70975; - var z_v = 0.999696; + vector_d a_d(2); + a_d << 3.70975, 1; + vector_d b_d(1); + b_d << 2.70975; + double z_d = 0.999696; - auto grad_2F1_tuple = grad_2F1(a_v[0], a_v[1], b_v[0], z_v); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto grad_2F1_tuple = grad_2F1(a_d[0], a_d[1], b_d[0], z_d); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(std::get<0>(grad_2F1_tuple), std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(std::get<1>(grad_2F1_tuple), std::get<0>(grad_tuple)[1]); @@ -305,22 +309,23 @@ TEST(RevMath, grad2F1_11) { EXPECT_FLOAT_EQ(std::get<3>(grad_2F1_tuple), std::get<2>(grad_tuple)); } -TEST(RevMath, F32_converges_by_z) { +TEST(PrimMath, F32_converges_by_z) { using stan::math::grad_F32; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(3); - a_v << 1.0, 1.0, 1.0; - vector_v b_v(2); - b_v << 1.0, 1.0; - var z_v = 0.6; + vector_d a_d(3); + a_d << 1.0, 1.0, 1.0; + vector_d b_d(2); + b_d << 1.0, 1.0; + double z_d = 0.6; double g_calc[6]; grad_F32(g_calc, 1.0, 1.0, 1.0, 1.0, 1.0, 0.6, 1e-10); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); @@ -330,22 +335,23 @@ TEST(RevMath, F32_converges_by_z) { EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); } -TEST(RevMath, grad_F32_double_sign_flip_1) { +TEST(PrimMath, grad_F32_double_sign_flip_1) { using stan::math::grad_F32; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(3); - a_v << 1.0, -0.5, -2.5; - vector_v b_v(2); - b_v << 10.0, 1.0; - var z_v = 0.3; + vector_d a_d(3); + a_d << 1.0, -0.5, -2.5; + vector_d b_d(2); + b_d << 10.0, 1.0; + double z_d = 0.3; double g_calc[6]; grad_F32(g_calc, 1.0, -.5, -2.5, 10.0, 1.0, 0.3, 1e-10); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); @@ -355,22 +361,23 @@ TEST(RevMath, grad_F32_double_sign_flip_1) { EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); } -TEST(RevMath, grad_F32_double_sign_flip_2) { +TEST(PrimMath, grad_F32_double_sign_flip_2) { using stan::math::grad_F32; using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(3); - a_v << 1.0, -0.5, -4.5; - vector_v b_v(2); - b_v << 10.0, 1.0; - var z_v = 0.3; + vector_d a_d(3); + a_d << 1.0, -0.5, -4.5; + vector_d b_d(2); + b_d << 10.0, 1.0; + double z_d = 0.3; double g_calc[6]; grad_F32(g_calc, 1.0, -.5, -4.5, 10.0, 1.0, 0.3, 1e-10); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); @@ -380,18 +387,19 @@ TEST(RevMath, grad_F32_double_sign_flip_2) { EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); } -TEST(RevMath, grad_2F1_negative_z) { +TEST(PrimMath, grad_2F1_negative_z) { using stan::math::grad_pFq; - using stan::math::var; - using stan::math::vector_v; + using stan::math::hypergeometric_pFq; + using stan::math::vector_d; - vector_v a_v(2); - a_v << 3.70975, 1.0; - vector_v b_v(1); - b_v << 2.70975; - var z_v = -0.2; + vector_d a_d(2); + a_d << 3.70975, 1.0; + vector_d b_d(1); + b_d << 2.70975; + double z_d = -0.2; - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(-0.0488658806159776, std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(-0.193844936204681, std::get<0>(grad_tuple)[1]); @@ -399,26 +407,26 @@ TEST(RevMath, grad_2F1_negative_z) { EXPECT_FLOAT_EQ(0.865295247272367, std::get<2>(grad_tuple)); } -TEST(RevMath, grad_3F2_cross_zero) { +TEST(PrimMath, grad_3F2_cross_zero) { using stan::math::grad_F32; using stan::math::grad_pFq; + using stan::math::hypergeometric_pFq; using stan::math::hypergeometric_3F2; - using stan::math::var; using stan::math::vector_d; - using stan::math::vector_v; - vector_v a_v(3); - a_v << 1, 1, -1; + vector_d a_d(3); + a_d << 1, 1, -1; - vector_v b_v(2); - b_v << 2, 2; + vector_d b_d(2); + b_d << 2, 2; - var z_v = 0.292893; + double z_d = 0.292893; double g_calc[6]; grad_F32(g_calc, 1.0, 1.0, -1.0, 2.0, 2.0, 0.292893); - auto grad_tuple = grad_pFq(a_v, b_v, z_v); + auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); + auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); From 644f9cb84f610c1c53391f4e4ae5337b6d818556 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 7 Oct 2023 18:15:45 +0300 Subject: [PATCH 13/36] Update doc --- stan/math/prim/fun/grad_pFq.hpp | 52 +++++++++++++++++---------------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 0b31a970e16..d0fad3c3a64 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -2,58 +2,60 @@ #define STAN_MATH_PRIM_FUN_GRAD_PFQ_HPP #include -#include #include -#include #include +#include +#include +#include #include -#include #include -#include -#include +#include +#include namespace stan { namespace math { namespace internal { template inline auto binarysign(const T& x) { - return select(x == 0.0, 1, sign(value_of_rec(x))); + return select(x == 0.0, 1.0, sign(value_of_rec(x))); } } /** - * Wrapper function for calculating gradients for the generalized - * hypergeometric function. The function always returns a tuple with - * three elements (gradients wrt a, b, and z, respectively), but the - * elements will only be defined/calculated when the respective parameter - * is not a primitive type. + * Returns the gradient of generalized hypergeometric function wrt to the + * input arguments: + * \f$ _pF_q(a_1,...,a_p;b_1,...,b_q;z) \f$ * + * @tparam CalcA Boolean for whether to calculate derivatives wrt to 'a' + * @tparam CalcB Boolean for whether to calculate derivatives wrt to 'b' + * @tparam CalcZ Boolean for whether to calculate derivatives wrt to 'z' + * @tparam TpFq Scalar type of hypergeometric_pFq return value * @tparam Ta Eigen type with either one row or column at compile time * @tparam Tb Eigen type with either one row or column at compile time * @tparam Tz Scalar type + * @param[in] pfq_val Return value from the hypergeometric_pFq function * @param[in] a Vector of 'a' arguments to function * @param[in] b Vector of 'b' arguments to function * @param[in] z Scalar z argument * @param[in] precision Convergence criteria for infinite sum - * @param[in] outer_steps Maximum number of iterations + * @param[in] max_steps Maximum number of iterations * @return Tuple of gradients */ template -auto grad_pFq(const TpFq& pfq_val, const Ta& a_val, const Tb& b_val, - const Tz& z_val, double precision = 1e-14, int max_steps = 1e6) { +auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, + const Tz& z, double precision = 1e-14, int max_steps = 1e6) { using std::max; using T_Rtn = return_type_t; using Ta_Array = Eigen::Array, -1, 1>; using Tb_Array = Eigen::Array, -1, 1>; - using T_RtnVector = Eigen::Matrix; - Ta_Array a_array = as_column_vector_or_scalar(a_val).array(); + Ta_Array a_array = as_column_vector_or_scalar(a).array(); Ta_Array a_k = a_array; - Tb_Array b_array = as_column_vector_or_scalar(b_val).array(); + Tb_Array b_array = as_column_vector_or_scalar(b).array(); Tb_Array b_k = b_array; - Tz log_z = log(abs(z_val)); - int z_sign = internal::binarysign(z_val); + Tz log_z = log(abs(z)); + int z_sign = internal::binarysign(z); Ta_Array digamma_a = inv(a_k); Tb_Array digamma_b = inv(b_k); @@ -61,12 +63,12 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a_val, const Tb& b_val, std::tuple< promote_scalar_t>, promote_scalar_t>, T_Rtn> ret_tuple; - std::get<0>(ret_tuple).setConstant(a_val.size(), 0.0); - std::get<1>(ret_tuple).setConstant(b_val.size(), 0.0); + std::get<0>(ret_tuple).setConstant(a.size(), 0.0); + std::get<1>(ret_tuple).setConstant(b.size(), 0.0); std::get<2>(ret_tuple) = 0.0; - Eigen::Array a_grad(a_val.size()); - Eigen::Array b_grad(b_val.size()); + Eigen::Array a_grad(a.size()); + Eigen::Array b_grad(b.size()); int k = 0; int base_sign = 1; @@ -110,8 +112,8 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a_val, const Tb& b_val, } if (CalcZ) { T_Rtn pfq_p1_val = hypergeometric_pFq((a_array + 1).matrix(), - (b_array + 1).matrix(), z_val); - std::get<2>(ret_tuple) = prod(a_val) / prod(b_val) * pfq_p1_val; + (b_array + 1).matrix(), z); + std::get<2>(ret_tuple) = prod(a) / prod(b) * pfq_p1_val; } return ret_tuple; } From f7131531072b5f0098283ee30d024178bf7a12d0 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 7 Oct 2023 11:22:57 -0400 Subject: [PATCH 14/36] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/fwd/fun/hypergeometric_pFq.hpp | 11 ++++----- stan/math/prim/fun/grad_pFq.hpp | 28 +++++++++++------------ stan/math/rev/fun/hypergeometric_pFq.hpp | 17 +++++++------- test/unit/math/prim/fun/grad_pFq_test.cpp | 3 +-- 4 files changed, 28 insertions(+), 31 deletions(-) diff --git a/stan/math/fwd/fun/hypergeometric_pFq.hpp b/stan/math/fwd/fun/hypergeometric_pFq.hpp index f63164a9d13..bb8a87afeb1 100644 --- a/stan/math/fwd/fun/hypergeometric_pFq.hpp +++ b/stan/math/fwd/fun/hypergeometric_pFq.hpp @@ -21,12 +21,11 @@ namespace math { * @param[in] z Scalar z argument * @return Generalized hypergeometric function */ -template ::value, - bool GradB = !is_constant::value, - bool GradZ = !is_constant::value, - require_all_matrix_t* = nullptr, - require_return_type_t* = nullptr> +template < + typename Ta, typename Tb, typename Tz, bool GradA = !is_constant::value, + bool GradB = !is_constant::value, bool GradZ = !is_constant::value, + require_all_matrix_t* = nullptr, + require_return_type_t* = nullptr> inline return_type_t hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { using fvar_t = return_type_t; diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index d0fad3c3a64..36845739f95 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -15,11 +15,11 @@ namespace stan { namespace math { namespace internal { - template - inline auto binarysign(const T& x) { - return select(x == 0.0, 1.0, sign(value_of_rec(x))); - } +template +inline auto binarysign(const T& x) { + return select(x == 0.0, 1.0, sign(value_of_rec(x))); } +} // namespace internal /** * Returns the gradient of generalized hypergeometric function wrt to the @@ -43,8 +43,8 @@ namespace internal { */ template -auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, - const Tz& z, double precision = 1e-14, int max_steps = 1e6) { +auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, + double precision = 1e-14, int max_steps = 1e6) { using std::max; using T_Rtn = return_type_t; using Ta_Array = Eigen::Array, -1, 1>; @@ -60,9 +60,9 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, Ta_Array digamma_a = inv(a_k); Tb_Array digamma_b = inv(b_k); - std::tuple< - promote_scalar_t>, - promote_scalar_t>, T_Rtn> ret_tuple; + std::tuple>, + promote_scalar_t>, T_Rtn> + ret_tuple; std::get<0>(ret_tuple).setConstant(a.size(), 0.0); std::get<1>(ret_tuple).setConstant(b.size(), 0.0); std::get<2>(ret_tuple) = 0.0; @@ -79,7 +79,7 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, if (CalcA) { a_grad = log(select(digamma_a == 0.0, 1.0, abs(digamma_a))) + log_base; std::get<0>(ret_tuple).array() - += exp(a_grad) * base_sign * internal::binarysign(digamma_a); + += exp(a_grad) * base_sign * internal::binarysign(digamma_a); curr_log_prec = max(curr_log_prec, a_grad.maxCoeff()); digamma_a += select(a_k == 0.0, 0.0, inv(a_k)); @@ -88,7 +88,7 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, if (CalcB) { b_grad = log(select(digamma_b == 0.0, 1.0, abs(digamma_b))) + log_base; std::get<1>(ret_tuple).array() - -= exp(b_grad) * base_sign * internal::binarysign(digamma_b); + -= exp(b_grad) * base_sign * internal::binarysign(digamma_b); curr_log_prec = max(curr_log_prec, b_grad.maxCoeff()); digamma_b += select(b_k == 0.0, 0.0, inv(b_k)); @@ -97,7 +97,7 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, log_base += (sum(log(select(a_k == 0.0, 0.0, abs(a_k)))) + log_z) - (sum(log(select(b_k == 0.0, 0.0, abs(b_k)))) + log1p(k)); base_sign *= z_sign * internal::binarysign(a_k).prod() - * internal::binarysign(b_k).prod(); + * internal::binarysign(b_k).prod(); a_k += 1.0; b_k += 1.0; @@ -111,8 +111,8 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, std::get<1>(ret_tuple).array() += pfq_val / b_array; } if (CalcZ) { - T_Rtn pfq_p1_val = hypergeometric_pFq((a_array + 1).matrix(), - (b_array + 1).matrix(), z); + T_Rtn pfq_p1_val + = hypergeometric_pFq((a_array + 1).matrix(), (b_array + 1).matrix(), z); std::get<2>(ret_tuple) = prod(a) / prod(b) * pfq_p1_val; } return ret_tuple; diff --git a/stan/math/rev/fun/hypergeometric_pFq.hpp b/stan/math/rev/fun/hypergeometric_pFq.hpp index 92dc7baff71..fab219b4f04 100644 --- a/stan/math/rev/fun/hypergeometric_pFq.hpp +++ b/stan/math/rev/fun/hypergeometric_pFq.hpp @@ -21,20 +21,19 @@ namespace math { * @param[in] z Scalar z argument * @return Generalized hypergeometric function */ -template ::value, - bool GradB = !is_constant::value, - bool GradZ = !is_constant::value, - require_all_matrix_t* = nullptr, - require_return_type_t* = nullptr> +template < + typename Ta, typename Tb, typename Tz, bool GradA = !is_constant::value, + bool GradB = !is_constant::value, bool GradZ = !is_constant::value, + require_all_matrix_t* = nullptr, + require_return_type_t* = nullptr> inline var hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { arena_t arena_a = a; arena_t arena_b = b; auto pfq_val = hypergeometric_pFq(a.val(), b.val(), value_of(z)); auto grad_tuple - = grad_pFq(pfq_val, a.val(), b.val(), value_of(z)); - return make_callback_var(pfq_val, - [arena_a, arena_b, z, grad_tuple](auto& vi) mutable { + = grad_pFq(pfq_val, a.val(), b.val(), value_of(z)); + return make_callback_var( + pfq_val, [arena_a, arena_b, z, grad_tuple](auto& vi) mutable { if (GradA) { forward_as>(arena_a).adj() += vi.adj() * std::get<0>(grad_tuple); diff --git a/test/unit/math/prim/fun/grad_pFq_test.cpp b/test/unit/math/prim/fun/grad_pFq_test.cpp index 2f8c0cd5ef1..bb206a4cf8c 100644 --- a/test/unit/math/prim/fun/grad_pFq_test.cpp +++ b/test/unit/math/prim/fun/grad_pFq_test.cpp @@ -298,7 +298,6 @@ TEST(PrimMath, grad2F1_11) { b_d << 2.70975; double z_d = 0.999696; - auto grad_2F1_tuple = grad_2F1(a_d[0], a_d[1], b_d[0], z_d); auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); @@ -410,8 +409,8 @@ TEST(PrimMath, grad_2F1_negative_z) { TEST(PrimMath, grad_3F2_cross_zero) { using stan::math::grad_F32; using stan::math::grad_pFq; - using stan::math::hypergeometric_pFq; using stan::math::hypergeometric_3F2; + using stan::math::hypergeometric_pFq; using stan::math::vector_d; vector_d a_d(3); From 3c5b15d8d6bd5037fcfbdd04af19bf67856e63dd Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 7 Oct 2023 18:25:33 +0300 Subject: [PATCH 15/36] Skip infsum if not needed --- stan/math/prim/fun/grad_pFq.hpp | 72 +++++++++++++++++---------------- 1 file changed, 37 insertions(+), 35 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 36845739f95..9f2c30eec3a 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -67,48 +67,50 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, std::get<1>(ret_tuple).setConstant(b.size(), 0.0); std::get<2>(ret_tuple) = 0.0; - Eigen::Array a_grad(a.size()); - Eigen::Array b_grad(b.size()); + if (CalcA || CalcB) { + Eigen::Array a_grad(a.size()); + Eigen::Array b_grad(b.size()); - int k = 0; - int base_sign = 1; - T_Rtn curr_log_prec = NEGATIVE_INFTY; - T_Rtn log_base = 0; - while ((k < 10 || curr_log_prec > log(precision)) && (k <= max_steps)) { - curr_log_prec = NEGATIVE_INFTY; - if (CalcA) { - a_grad = log(select(digamma_a == 0.0, 1.0, abs(digamma_a))) + log_base; - std::get<0>(ret_tuple).array() - += exp(a_grad) * base_sign * internal::binarysign(digamma_a); + int k = 0; + int base_sign = 1; + T_Rtn curr_log_prec = NEGATIVE_INFTY; + T_Rtn log_base = 0; + while ((k < 10 || curr_log_prec > log(precision)) && (k <= max_steps)) { + curr_log_prec = NEGATIVE_INFTY; + if (CalcA) { + a_grad = log(select(digamma_a == 0.0, 1.0, abs(digamma_a))) + log_base; + std::get<0>(ret_tuple).array() + += exp(a_grad) * base_sign * internal::binarysign(digamma_a); - curr_log_prec = max(curr_log_prec, a_grad.maxCoeff()); - digamma_a += select(a_k == 0.0, 0.0, inv(a_k)); - } + curr_log_prec = max(curr_log_prec, a_grad.maxCoeff()); + digamma_a += select(a_k == 0.0, 0.0, inv(a_k)); + } - if (CalcB) { - b_grad = log(select(digamma_b == 0.0, 1.0, abs(digamma_b))) + log_base; - std::get<1>(ret_tuple).array() - -= exp(b_grad) * base_sign * internal::binarysign(digamma_b); + if (CalcB) { + b_grad = log(select(digamma_b == 0.0, 1.0, abs(digamma_b))) + log_base; + std::get<1>(ret_tuple).array() + -= exp(b_grad) * base_sign * internal::binarysign(digamma_b); - curr_log_prec = max(curr_log_prec, b_grad.maxCoeff()); - digamma_b += select(b_k == 0.0, 0.0, inv(b_k)); - } + curr_log_prec = max(curr_log_prec, b_grad.maxCoeff()); + digamma_b += select(b_k == 0.0, 0.0, inv(b_k)); + } - log_base += (sum(log(select(a_k == 0.0, 0.0, abs(a_k)))) + log_z) - - (sum(log(select(b_k == 0.0, 0.0, abs(b_k)))) + log1p(k)); - base_sign *= z_sign * internal::binarysign(a_k).prod() - * internal::binarysign(b_k).prod(); + log_base += (sum(log(select(a_k == 0.0, 0.0, abs(a_k)))) + log_z) + - (sum(log(select(b_k == 0.0, 0.0, abs(b_k)))) + log1p(k)); + base_sign *= z_sign * internal::binarysign(a_k).prod() + * internal::binarysign(b_k).prod(); - a_k += 1.0; - b_k += 1.0; - k += 1; - } + a_k += 1.0; + b_k += 1.0; + k += 1; + } - if (CalcA) { - std::get<0>(ret_tuple).array() -= pfq_val / a_array; - } - if (CalcB) { - std::get<1>(ret_tuple).array() += pfq_val / b_array; + if (CalcA) { + std::get<0>(ret_tuple).array() -= pfq_val / a_array; + } + if (CalcB) { + std::get<1>(ret_tuple).array() += pfq_val / b_array; + } } if (CalcZ) { T_Rtn pfq_p1_val From 50acf603d545f4e77a4c7e1e2baf065d7e344bef Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 7 Oct 2023 11:26:36 -0400 Subject: [PATCH 16/36] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/grad_pFq.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 9f2c30eec3a..f9d5b219c7a 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -98,7 +98,7 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, log_base += (sum(log(select(a_k == 0.0, 0.0, abs(a_k)))) + log_z) - (sum(log(select(b_k == 0.0, 0.0, abs(b_k)))) + log1p(k)); base_sign *= z_sign * internal::binarysign(a_k).prod() - * internal::binarysign(b_k).prod(); + * internal::binarysign(b_k).prod(); a_k += 1.0; b_k += 1.0; From c1cbe26b86dfc7a1d061b348794a9d9400f6f373 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 7 Oct 2023 18:42:12 +0300 Subject: [PATCH 17/36] Test fix asan error --- stan/math/rev/fun/hypergeometric_pFq.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/stan/math/rev/fun/hypergeometric_pFq.hpp b/stan/math/rev/fun/hypergeometric_pFq.hpp index fab219b4f04..d720d619fc4 100644 --- a/stan/math/rev/fun/hypergeometric_pFq.hpp +++ b/stan/math/rev/fun/hypergeometric_pFq.hpp @@ -30,10 +30,11 @@ inline var hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { arena_t arena_a = a; arena_t arena_b = b; auto pfq_val = hypergeometric_pFq(a.val(), b.val(), value_of(z)); - auto grad_tuple - = grad_pFq(pfq_val, a.val(), b.val(), value_of(z)); return make_callback_var( - pfq_val, [arena_a, arena_b, z, grad_tuple](auto& vi) mutable { + pfq_val, [arena_a, arena_b, z, pfq_val](auto& vi) mutable { + auto grad_tuple = grad_pFq(pfq_val, arena_a.val(), + arena_b.val(), + value_of(z)); if (GradA) { forward_as>(arena_a).adj() += vi.adj() * std::get<0>(grad_tuple); From 07b5f4c8f772844db3b7a5c1074da10af922f083 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 7 Oct 2023 11:43:20 -0400 Subject: [PATCH 18/36] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/rev/fun/hypergeometric_pFq.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/stan/math/rev/fun/hypergeometric_pFq.hpp b/stan/math/rev/fun/hypergeometric_pFq.hpp index d720d619fc4..b0a55e2e11d 100644 --- a/stan/math/rev/fun/hypergeometric_pFq.hpp +++ b/stan/math/rev/fun/hypergeometric_pFq.hpp @@ -32,9 +32,8 @@ inline var hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { auto pfq_val = hypergeometric_pFq(a.val(), b.val(), value_of(z)); return make_callback_var( pfq_val, [arena_a, arena_b, z, pfq_val](auto& vi) mutable { - auto grad_tuple = grad_pFq(pfq_val, arena_a.val(), - arena_b.val(), - value_of(z)); + auto grad_tuple = grad_pFq( + pfq_val, arena_a.val(), arena_b.val(), value_of(z)); if (GradA) { forward_as>(arena_a).adj() += vi.adj() * std::get<0>(grad_tuple); From f810c558969cd3a501316a6e70a0144ec3ae9d31 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 8 Oct 2023 11:06:42 +0300 Subject: [PATCH 19/36] Fix handling of zero-inputs --- stan/math/prim/fun/grad_pFq.hpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index f9d5b219c7a..2847a76cf1e 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -57,8 +57,8 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, Tz log_z = log(abs(z)); int z_sign = internal::binarysign(z); - Ta_Array digamma_a = inv(a_k); - Tb_Array digamma_b = inv(b_k); + Ta_Array digamma_a = select(a_k == 0.0, 0.0, inv(a_k)); + Tb_Array digamma_b = select(b_k == 0.0, 0.0, inv(b_k)); std::tuple>, promote_scalar_t>, T_Rtn> @@ -106,10 +106,12 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, } if (CalcA) { - std::get<0>(ret_tuple).array() -= pfq_val / a_array; + std::get<0>(ret_tuple).array() + -= select(a_array == 0.0, 0.0, pfq_val / a_array); } if (CalcB) { - std::get<1>(ret_tuple).array() += pfq_val / b_array; + std::get<1>(ret_tuple).array() + += select(b_array == 0.0, 0.0, pfq_val / b_array); } } if (CalcZ) { From 535dadb665381c6e7a71557a930e1e4a6a368457 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 8 Oct 2023 04:07:50 -0400 Subject: [PATCH 20/36] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/grad_pFq.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 2847a76cf1e..19e04573210 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -107,11 +107,11 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, if (CalcA) { std::get<0>(ret_tuple).array() - -= select(a_array == 0.0, 0.0, pfq_val / a_array); + -= select(a_array == 0.0, 0.0, pfq_val / a_array); } if (CalcB) { std::get<1>(ret_tuple).array() - += select(b_array == 0.0, 0.0, pfq_val / b_array); + += select(b_array == 0.0, 0.0, pfq_val / b_array); } } if (CalcZ) { From e76963032a7c087fab4b31c28c245847eb6525e3 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 9 Oct 2023 12:27:50 +0300 Subject: [PATCH 21/36] Update doc, additional simplification --- stan/math/prim/fun/grad_pFq.hpp | 58 +++++++++++++++++++++++++++++---- 1 file changed, 52 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 19e04573210..11feaa9b0ac 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -26,6 +26,54 @@ inline auto binarysign(const T& x) { * input arguments: * \f$ _pF_q(a_1,...,a_p;b_1,...,b_q;z) \f$ * + * Where: + * \f$ + * \frac{\partial }{\partial a_1} = + * \sum_{k=0}^{\infty}{ + * \frac + * {\psi\left(k+a_1\right)\left(\prod_{j=1}^p\left(a_j\right)_k\right)z^k} + * {k!\prod_{j=1}^q\left(b_j\right)_k}} + * - \psi\left(a_1\right)_pF_q(a_1,...,a_p;b_1,...,b_q;z) + * \f$ + * \f$ + * \frac{\partial }{\partial b_1} = + * \psi\left(b_1\right)_pF_q(a_1,...,a_p;b_1,...,b_q;z) - + * \sum_{k=0}^{\infty}{ + * \frac + * {\psi\left(k+b_1\right)\left(\prod_{j=1}^p\left(a_j\right)_k\right)z^k} + * {k!\prod_{j=1}^q\left(b_j\right)_k}} + * \f$ + * + * \f$ + * \frac{\partial }{\partial z} = + * \frac{\prod_{j=1}^{p}(a_j)}{\prod_{j=1}^{q} (b_j)}\ + * * _pF_q(a_1+1,...,a_p+1;b_1+1,...,b_q+1;z) + * \f$ + * + * Noting the the recurrence relation for the digamma function: + * \f$ \psi(x + 1) = \psi(x) + \frac{1}{x} \f$, and as such the presence of the + * digamma function in both operands of the subtraction, this then becomes a + * scaling factor and can be removed. The gradients for the function w.r.t a & b + * then simplify to: + * \f$ + * \frac{\partial }{\partial a_1} = + * \sum_{k=1}^{\infty}{ + * \frac + * {\left(1 + \sum_{m=0}^{k-1}\frac{1}{m+a_1}\right) + * * (\prod_{j=1}^p\left(a_j\right)_k\right)z^k} + * {k!\prod_{j=1}^q\left(b_j\right)_k}} + * - _pF_q(a_1,...,a_p;b_1,...,b_q;z) + * \f$ + * \f$ + * \frac{\partial }{\partial b_1} = + * _pF_q(a_1,...,a_p;b_1,...,b_q;z) - + * \sum_{k=1}^{\infty}{ + * \frac + * {\left(1 + \sum_{m=0}^{k-1}\frac{1}{m+b_1}\right) + * * \left(\prod_{j=1}^p\left(a_j\right)_k\right)z^k} + * {k!\prod_{j=1}^q\left(b_j\right)_k}} + * \f$ + * * @tparam CalcA Boolean for whether to calculate derivatives wrt to 'a' * @tparam CalcB Boolean for whether to calculate derivatives wrt to 'b' * @tparam CalcZ Boolean for whether to calculate derivatives wrt to 'z' @@ -57,8 +105,8 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, Tz log_z = log(abs(z)); int z_sign = internal::binarysign(z); - Ta_Array digamma_a = select(a_k == 0.0, 0.0, inv(a_k)); - Tb_Array digamma_b = select(b_k == 0.0, 0.0, inv(b_k)); + Ta_Array digamma_a = Ta_Array::Ones(a.size()); + Tb_Array digamma_b = Tb_Array::Ones(b.size()); std::tuple>, promote_scalar_t>, T_Rtn> @@ -106,12 +154,10 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, } if (CalcA) { - std::get<0>(ret_tuple).array() - -= select(a_array == 0.0, 0.0, pfq_val / a_array); + std::get<0>(ret_tuple).array() -= pfq_val; } if (CalcB) { - std::get<1>(ret_tuple).array() - += select(b_array == 0.0, 0.0, pfq_val / b_array); + std::get<1>(ret_tuple).array() += pfq_val; } } if (CalcZ) { From a0546651b7bb1be7f13e31bbc798657252ecdd7c Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 9 Oct 2023 12:36:21 +0300 Subject: [PATCH 22/36] Additional simplification --- stan/math/prim/fun/grad_pFq.hpp | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 11feaa9b0ac..77040ea849e 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -111,8 +111,8 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, std::tuple>, promote_scalar_t>, T_Rtn> ret_tuple; - std::get<0>(ret_tuple).setConstant(a.size(), 0.0); - std::get<1>(ret_tuple).setConstant(b.size(), 0.0); + std::get<0>(ret_tuple).setConstant(a.size(), -pfq_val); + std::get<1>(ret_tuple).setConstant(b.size(), pfq_val); std::get<2>(ret_tuple) = 0.0; if (CalcA || CalcB) { @@ -152,13 +152,6 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, b_k += 1.0; k += 1; } - - if (CalcA) { - std::get<0>(ret_tuple).array() -= pfq_val; - } - if (CalcB) { - std::get<1>(ret_tuple).array() += pfq_val; - } } if (CalcZ) { T_Rtn pfq_p1_val From b9388b1b8f4775be3278b7aa1646999603f839d8 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 9 Oct 2023 13:07:41 +0300 Subject: [PATCH 23/36] Fix sign-flip case --- stan/math/prim/fun/grad_pFq.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 77040ea849e..448d8915607 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -126,7 +126,7 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, while ((k < 10 || curr_log_prec > log(precision)) && (k <= max_steps)) { curr_log_prec = NEGATIVE_INFTY; if (CalcA) { - a_grad = log(select(digamma_a == 0.0, 1.0, abs(digamma_a))) + log_base; + a_grad = log(abs(digamma_a)) + log_base; std::get<0>(ret_tuple).array() += exp(a_grad) * base_sign * internal::binarysign(digamma_a); @@ -135,7 +135,7 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, } if (CalcB) { - b_grad = log(select(digamma_b == 0.0, 1.0, abs(digamma_b))) + log_base; + b_grad = log(abs(digamma_b)) + log_base; std::get<1>(ret_tuple).array() -= exp(b_grad) * base_sign * internal::binarysign(digamma_b); @@ -143,8 +143,8 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, digamma_b += select(b_k == 0.0, 0.0, inv(b_k)); } - log_base += (sum(log(select(a_k == 0.0, 0.0, abs(a_k)))) + log_z) - - (sum(log(select(b_k == 0.0, 0.0, abs(b_k)))) + log1p(k)); + log_base += (sum(log(abs(a_k))) + log_z) + - (sum(log(abs(b_k))) + log1p(k)); base_sign *= z_sign * internal::binarysign(a_k).prod() * internal::binarysign(b_k).prod(); From fbcec85476ecd50aec96644affda454bf8cd180a Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 9 Oct 2023 06:08:45 -0400 Subject: [PATCH 24/36] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/grad_pFq.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 448d8915607..017cd29bb96 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -143,8 +143,8 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, digamma_b += select(b_k == 0.0, 0.0, inv(b_k)); } - log_base += (sum(log(abs(a_k))) + log_z) - - (sum(log(abs(b_k))) + log1p(k)); + log_base + += (sum(log(abs(a_k))) + log_z) - (sum(log(abs(b_k))) + log1p(k)); base_sign *= z_sign * internal::binarysign(a_k).prod() * internal::binarysign(b_k).prod(); From c149fdacf1111acb23eb6a9907773b757a0bcb1d Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 21 Feb 2024 14:19:12 +0200 Subject: [PATCH 25/36] Avoid unneeded computation --- stan/math/prim/fun/grad_pFq.hpp | 60 ++++++++++++++++----------------- 1 file changed, 29 insertions(+), 31 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 017cd29bb96..b71bffabc77 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -14,12 +14,6 @@ namespace stan { namespace math { -namespace internal { -template -inline auto binarysign(const T& x) { - return select(x == 0.0, 1.0, sign(value_of_rec(x))); -} -} // namespace internal /** * Returns the gradient of generalized hypergeometric function wrt to the @@ -51,10 +45,8 @@ inline auto binarysign(const T& x) { * \f$ * * Noting the the recurrence relation for the digamma function: - * \f$ \psi(x + 1) = \psi(x) + \frac{1}{x} \f$, and as such the presence of the - * digamma function in both operands of the subtraction, this then becomes a - * scaling factor and can be removed. The gradients for the function w.r.t a & b - * then simplify to: + * \f$ \psi(x + 1) = \psi(x) + \frac{1}{x} \f$, the gradients for the + * function w.r.t a & b then simplify to: * \f$ * \frac{\partial }{\partial a_1} = * \sum_{k=1}^{\infty}{ @@ -90,37 +82,43 @@ inline auto binarysign(const T& x) { * @return Tuple of gradients */ template -auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, - double precision = 1e-14, int max_steps = 1e6) { + typename TpFq, typename Ta, typename Tb, typename Tz, + typename T_Rtn = return_type_t, + typename Ta_Rtn = promote_scalar_t>, + typename Tb_Rtn = promote_scalar_t>> +std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, + const Tb& b, const Tz& z, + double precision = 1e-14, + int max_steps = 1e6) { using std::max; - using T_Rtn = return_type_t; using Ta_Array = Eigen::Array, -1, 1>; using Tb_Array = Eigen::Array, -1, 1>; Ta_Array a_array = as_column_vector_or_scalar(a).array(); - Ta_Array a_k = a_array; Tb_Array b_array = as_column_vector_or_scalar(b).array(); - Tb_Array b_k = b_array; - Tz log_z = log(abs(z)); - int z_sign = internal::binarysign(z); - - Ta_Array digamma_a = Ta_Array::Ones(a.size()); - Tb_Array digamma_b = Tb_Array::Ones(b.size()); - - std::tuple>, - promote_scalar_t>, T_Rtn> - ret_tuple; - std::get<0>(ret_tuple).setConstant(a.size(), -pfq_val); - std::get<1>(ret_tuple).setConstant(b.size(), pfq_val); + + std::tuple ret_tuple; + std::get<0>(ret_tuple).setConstant(a.size(), 0); + std::get<1>(ret_tuple).setConstant(b.size(), 0); std::get<2>(ret_tuple) = 0.0; if (CalcA || CalcB) { + std::get<0>(ret_tuple).setConstant(-pfq_val); + std::get<1>(ret_tuple).setConstant(pfq_val); Eigen::Array a_grad(a.size()); Eigen::Array b_grad(b.size()); int k = 0; int base_sign = 1; + + Ta_Array a_k = a_array; + Tb_Array b_k = b_array; + Tz log_z = log(abs(z)); + int z_sign = sign(value_of_rec(z)); + + Ta_Array digamma_a = Ta_Array::Ones(a.size()); + Tb_Array digamma_b = Tb_Array::Ones(b.size()); + T_Rtn curr_log_prec = NEGATIVE_INFTY; T_Rtn log_base = 0; while ((k < 10 || curr_log_prec > log(precision)) && (k <= max_steps)) { @@ -128,7 +126,7 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, if (CalcA) { a_grad = log(abs(digamma_a)) + log_base; std::get<0>(ret_tuple).array() - += exp(a_grad) * base_sign * internal::binarysign(digamma_a); + += exp(a_grad) * base_sign * sign(value_of_rec(digamma_a)); curr_log_prec = max(curr_log_prec, a_grad.maxCoeff()); digamma_a += select(a_k == 0.0, 0.0, inv(a_k)); @@ -137,7 +135,7 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, if (CalcB) { b_grad = log(abs(digamma_b)) + log_base; std::get<1>(ret_tuple).array() - -= exp(b_grad) * base_sign * internal::binarysign(digamma_b); + -= exp(b_grad) * base_sign * sign(value_of_rec(digamma_b)); curr_log_prec = max(curr_log_prec, b_grad.maxCoeff()); digamma_b += select(b_k == 0.0, 0.0, inv(b_k)); @@ -145,8 +143,8 @@ auto grad_pFq(const TpFq& pfq_val, const Ta& a, const Tb& b, const Tz& z, log_base += (sum(log(abs(a_k))) + log_z) - (sum(log(abs(b_k))) + log1p(k)); - base_sign *= z_sign * internal::binarysign(a_k).prod() - * internal::binarysign(b_k).prod(); + base_sign *= z_sign * sign(value_of_rec(a_k)).prod() + * sign(value_of_rec(b_k)).prod(); a_k += 1.0; b_k += 1.0; From 2869827e1e72368b3c3cec0260f0fa5fa1b599be Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 21 Feb 2024 07:19:55 -0500 Subject: [PATCH 26/36] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/grad_pFq.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index b71bffabc77..619490f7308 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -87,9 +87,9 @@ template >, typename Tb_Rtn = promote_scalar_t>> std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, - const Tb& b, const Tz& z, - double precision = 1e-14, - int max_steps = 1e6) { + const Tb& b, const Tz& z, + double precision = 1e-14, + int max_steps = 1e6) { using std::max; using Ta_Array = Eigen::Array, -1, 1>; using Tb_Array = Eigen::Array, -1, 1>; From b97f02a43a942da1fbedba868c691e677ad35f31 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 24 Feb 2024 23:33:35 +0200 Subject: [PATCH 27/36] Tidying --- stan/math/fwd/fun/hypergeometric_pFq.hpp | 41 +++++++++++++----------- stan/math/prim/fun/grad_pFq.hpp | 7 ++-- 2 files changed, 25 insertions(+), 23 deletions(-) diff --git a/stan/math/fwd/fun/hypergeometric_pFq.hpp b/stan/math/fwd/fun/hypergeometric_pFq.hpp index bb8a87afeb1..8e5636478f6 100644 --- a/stan/math/fwd/fun/hypergeometric_pFq.hpp +++ b/stan/math/fwd/fun/hypergeometric_pFq.hpp @@ -22,37 +22,42 @@ namespace math { * @return Generalized hypergeometric function */ template < - typename Ta, typename Tb, typename Tz, bool GradA = !is_constant::value, - bool GradB = !is_constant::value, bool GradZ = !is_constant::value, - require_all_matrix_t* = nullptr, - require_return_type_t* = nullptr> -inline return_type_t hypergeometric_pFq(const Ta& a, const Tb& b, - const Tz& z) { - using fvar_t = return_type_t; - ref_type_t a_ref = a; - ref_type_t b_ref = b; + typename Ta, typename Tb, typename Tz, + typename FvarT = return_type_t, + bool GradA = !is_constant::value, + bool GradB = !is_constant::value, + bool GradZ = !is_constant::value, + require_all_vector_t* = nullptr, + require_return_type_t* = nullptr> +inline FvarT hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { + using PartialsT = partials_type_t; + using ARefT = ref_type_t; + using BRefT = ref_type_t; + + ARefT a_ref = a; + BRefT b_ref = b; auto&& a_val = value_of(a_ref); auto&& b_val = value_of(b_ref); auto&& z_val = value_of(z); - auto pfq_val = hypergeometric_pFq(a_val, b_val, z_val); + PartialsT pfq_val = hypergeometric_pFq(a_val, b_val, z_val); auto grad_tuple = grad_pFq(pfq_val, a_val, b_val, z_val); - typename fvar_t::Scalar grad = 0; + FvarT rtn = FvarT(pfq_val, 0.0); if (GradA) { - grad += dot_product(forward_as>(a_ref).d(), - std::get<0>(grad_tuple)); + rtn.d_ += dot_product(forward_as>(a_ref).d(), + std::get<0>(grad_tuple)); } if (GradB) { - grad += dot_product(forward_as>(b_ref).d(), - std::get<1>(grad_tuple)); + rtn.d_ += dot_product(forward_as>(b_ref).d(), + std::get<1>(grad_tuple)); } if (GradZ) { - grad += forward_as>(z).d_ - * std::get<2>(grad_tuple); + rtn.d_ += forward_as>(z).d_ + * std::get<2>(grad_tuple); } - return fvar_t(pfq_val, grad); + return rtn; } } // namespace math diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 619490f7308..2313b7972f5 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -98,13 +98,10 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, Tb_Array b_array = as_column_vector_or_scalar(b).array(); std::tuple ret_tuple; - std::get<0>(ret_tuple).setConstant(a.size(), 0); - std::get<1>(ret_tuple).setConstant(b.size(), 0); - std::get<2>(ret_tuple) = 0.0; if (CalcA || CalcB) { - std::get<0>(ret_tuple).setConstant(-pfq_val); - std::get<1>(ret_tuple).setConstant(pfq_val); + std::get<0>(ret_tuple).setConstant(a.size(), -pfq_val); + std::get<1>(ret_tuple).setConstant(b.size(), pfq_val); Eigen::Array a_grad(a.size()); Eigen::Array b_grad(b.size()); From af23ee62f60d417b94a5c111eef35f03b46983b4 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 24 Feb 2024 16:34:07 -0500 Subject: [PATCH 28/36] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/fwd/fun/hypergeometric_pFq.hpp | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/stan/math/fwd/fun/hypergeometric_pFq.hpp b/stan/math/fwd/fun/hypergeometric_pFq.hpp index 8e5636478f6..907d4ea32d5 100644 --- a/stan/math/fwd/fun/hypergeometric_pFq.hpp +++ b/stan/math/fwd/fun/hypergeometric_pFq.hpp @@ -21,14 +21,13 @@ namespace math { * @param[in] z Scalar z argument * @return Generalized hypergeometric function */ -template < - typename Ta, typename Tb, typename Tz, - typename FvarT = return_type_t, - bool GradA = !is_constant::value, - bool GradB = !is_constant::value, - bool GradZ = !is_constant::value, - require_all_vector_t* = nullptr, - require_return_type_t* = nullptr> +template , + bool GradA = !is_constant::value, + bool GradB = !is_constant::value, + bool GradZ = !is_constant::value, + require_all_vector_t* = nullptr, + require_return_type_t* = nullptr> inline FvarT hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { using PartialsT = partials_type_t; using ARefT = ref_type_t; @@ -46,15 +45,15 @@ inline FvarT hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { if (GradA) { rtn.d_ += dot_product(forward_as>(a_ref).d(), - std::get<0>(grad_tuple)); + std::get<0>(grad_tuple)); } if (GradB) { rtn.d_ += dot_product(forward_as>(b_ref).d(), - std::get<1>(grad_tuple)); + std::get<1>(grad_tuple)); } if (GradZ) { rtn.d_ += forward_as>(z).d_ - * std::get<2>(grad_tuple); + * std::get<2>(grad_tuple); } return rtn; From 30273c819150b5d23349bbab4f2d63c03f169ccc Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 26 Feb 2024 02:35:52 +0200 Subject: [PATCH 29/36] Update handling integers --- stan/math/prim/fun/grad_pFq.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 2313b7972f5..e2ecf20c1f1 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -126,7 +126,7 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, += exp(a_grad) * base_sign * sign(value_of_rec(digamma_a)); curr_log_prec = max(curr_log_prec, a_grad.maxCoeff()); - digamma_a += select(a_k == 0.0, 0.0, inv(a_k)); + digamma_a += inv(a_k); } if (CalcB) { @@ -135,7 +135,7 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, -= exp(b_grad) * base_sign * sign(value_of_rec(digamma_b)); curr_log_prec = max(curr_log_prec, b_grad.maxCoeff()); - digamma_b += select(b_k == 0.0, 0.0, inv(b_k)); + digamma_b += inv(b_k); } log_base @@ -143,8 +143,8 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, base_sign *= z_sign * sign(value_of_rec(a_k)).prod() * sign(value_of_rec(b_k)).prod(); - a_k += 1.0; - b_k += 1.0; + a_k = (a_k == 0.0).select(std::numeric_limits::min(), a_k + 1.0); + b_k = (b_k == 0.0).select(std::numeric_limits::min(), b_k + 1.0); k += 1; } } From b631517f121ec7bc5a1dcf9ef15529966dc72377 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 26 Feb 2024 09:29:14 +0200 Subject: [PATCH 30/36] Optimise handling of negative args --- stan/math/prim/fun/grad_pFq.hpp | 70 +++++++++--- .../math/fwd/fun/hypergeometric_pFq_test.cpp | 101 ------------------ .../math/rev/fun/hypergeometric_pFq_test.cpp | 28 ----- 3 files changed, 57 insertions(+), 142 deletions(-) delete mode 100644 test/unit/math/fwd/fun/hypergeometric_pFq_test.cpp delete mode 100644 test/unit/math/rev/fun/hypergeometric_pFq_test.cpp diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index e2ecf20c1f1..01c9a787bd9 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -3,12 +3,15 @@ #include #include +#include #include #include #include #include +#include #include #include +#include #include #include @@ -108,9 +111,24 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, int k = 0; int base_sign = 1; - Ta_Array a_k = a_array; - Tb_Array b_k = b_array; + static constexpr double dbl_min = std::numeric_limits::min(); + Ta_Array a_k = (a_array == 0.0).select(dbl_min, abs(a_array)); + Tb_Array b_k = (b_array == 0.0).select(dbl_min, abs(b_array)); Tz log_z = log(abs(z)); + + // Identify the number of iterations to needed for each element to sign + // flip from negative to positive - rather than checking at each iteration + Eigen::ArrayXi a_pos_k + = (a_array < 0.0).select(-floor(value_of_rec(a_array)), 0) + .template cast(); + int all_a_pos_k = a_pos_k.maxCoeff(); + Eigen::ArrayXi b_pos_k + = (b_array < 0.0).select(-floor(value_of_rec(b_array)), 0) + .template cast(); + int all_b_pos_k = b_pos_k.maxCoeff(); + Eigen::ArrayXi a_sign = select(a_pos_k == 0, 1, -1); + Eigen::ArrayXi b_sign = select(b_pos_k == 0, 1, -1); + int z_sign = sign(value_of_rec(z)); Ta_Array digamma_a = Ta_Array::Ones(a.size()); @@ -126,7 +144,7 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, += exp(a_grad) * base_sign * sign(value_of_rec(digamma_a)); curr_log_prec = max(curr_log_prec, a_grad.maxCoeff()); - digamma_a += inv(a_k); + digamma_a += inv(a_k) * a_sign; } if (CalcB) { @@ -135,23 +153,49 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, -= exp(b_grad) * base_sign * sign(value_of_rec(digamma_b)); curr_log_prec = max(curr_log_prec, b_grad.maxCoeff()); - digamma_b += inv(b_k); + digamma_b += inv(b_k) * b_sign; } - log_base - += (sum(log(abs(a_k))) + log_z) - (sum(log(abs(b_k))) + log1p(k)); - base_sign *= z_sign * sign(value_of_rec(a_k)).prod() - * sign(value_of_rec(b_k)).prod(); + log_base += (sum(log(a_k)) + log_z) - (sum(log(b_k)) + log1p(k)); + base_sign *= z_sign * a_sign.prod() * b_sign.prod(); + + // Wrap negative value handling in a conditional on iteration number so + // branch prediction likely to ignore once positive + if (k < all_a_pos_k) { + // Avoid log(0) and 1/0 in next iteration by using smallest double + // - This is smaller than EPSILON, so the following iteration will + // still be 1.0 + a_k = (a_k == 1.0 && a_sign == -1) + .select(dbl_min, (a_k < 1.0 && a_sign == -1) + .select(1.0 - a_k, a_k + 1.0 * a_sign)); + a_sign = select(k == a_pos_k - 1, 1, a_sign); + } else { + a_k += 1.0; + + if (k == all_a_pos_k) { + a_sign.setOnes(); + } + } + + if (k < all_a_pos_k) { + b_k = (b_k == 1.0 && b_sign == -1) + .select(dbl_min, (b_k < 1.0 && b_sign == -1) + .select(1.0 - b_k, b_k + 1.0 * b_sign)); + b_sign = select(k == b_pos_k - 1, 1, b_sign); + } else { + b_k += 1.0; + + if (k == all_b_pos_k) { + b_sign.setOnes(); + } + } - a_k = (a_k == 0.0).select(std::numeric_limits::min(), a_k + 1.0); - b_k = (b_k == 0.0).select(std::numeric_limits::min(), b_k + 1.0); k += 1; } } if (CalcZ) { - T_Rtn pfq_p1_val - = hypergeometric_pFq((a_array + 1).matrix(), (b_array + 1).matrix(), z); - std::get<2>(ret_tuple) = prod(a) / prod(b) * pfq_p1_val; + std::get<2>(ret_tuple) + = hypergeometric_pFq(add(a, 1.0), add(b, 1.0), z) * prod(a) / prod(b); } return ret_tuple; } diff --git a/test/unit/math/fwd/fun/hypergeometric_pFq_test.cpp b/test/unit/math/fwd/fun/hypergeometric_pFq_test.cpp deleted file mode 100644 index f8b6edb16b8..00000000000 --- a/test/unit/math/fwd/fun/hypergeometric_pFq_test.cpp +++ /dev/null @@ -1,101 +0,0 @@ -#include -#include -#include - -TEST(primScalFun, grad_2F2_fd) { - using stan::math::fvar; - using stan::math::hypergeometric_pFq; - using stan::math::vector_d; - using stan::math::vector_fd; - - vector_fd fd_a(2); - vector_d d_a(2); - fd_a.val() << 4, 2; - d_a << 4, 2; - fd_a.d() << 1, 1; - - vector_fd fd_b(2); - vector_d d_b(2); - fd_b << 6, 3; - d_b << 6, 3; - fd_b.d() << 1, 1; - - fvar fd_z = fvar(4, 1); - double d_z = 4; - - double a_adj = 3.924636646666071 + 6.897245961898751; - double b_adj = -2.775051002566842 - 4.980095849781222; - double z_adj = 4.916522138006060; - - // fvar, fvar, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(fd_a, fd_b, fd_z).d_, - a_adj + b_adj + z_adj); - - // fvar, fvar, double - EXPECT_FLOAT_EQ(hypergeometric_pFq(fd_a, fd_b, d_z).d_, a_adj + b_adj); - - // fvar, double, double - EXPECT_FLOAT_EQ(hypergeometric_pFq(fd_a, d_b, d_z).d_, a_adj); - - // fvar, double, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(fd_a, d_b, fd_z).d_, a_adj + z_adj); - - // double, fvar, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(d_a, fd_b, fd_z).d_, b_adj + z_adj); - - // double, fvar, double - EXPECT_FLOAT_EQ(hypergeometric_pFq(d_a, fd_b, d_z).d_, b_adj); - - // double, double, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(d_a, d_b, fd_z).d_, z_adj); -} - -TEST(primScalFun, grad_2F2_ffd) { - using stan::math::fvar; - using stan::math::hypergeometric_pFq; - using stan::math::vector_d; - using stan::math::vector_ffd; - - vector_ffd ffd_a(2); - vector_d d_a(2); - ffd_a.val() << 4, 2; - d_a << 4, 2; - ffd_a.val().d() << 1, 1; - - vector_ffd ffd_b(2); - vector_d d_b(2); - ffd_b.val() << 6, 3; - d_b << 6, 3; - ffd_b.val().d() << 1, 1; - - fvar> ffd_z; - ffd_z.val_.val_ = 4; - ffd_z.val_.d_ = 1; - double d_z = 4; - - double a_adj = 3.924636646666071 + 6.897245961898751; - double b_adj = -2.775051002566842 - 4.980095849781222; - double z_adj = 4.916522138006060; - - // fvar, fvar, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(ffd_a, ffd_b, ffd_z).val_.d_, - a_adj + b_adj + z_adj); - - // fvar, fvar, double - EXPECT_FLOAT_EQ(hypergeometric_pFq(ffd_a, ffd_b, d_z).val_.d_, a_adj + b_adj); - - // fvar, double, double - EXPECT_FLOAT_EQ(hypergeometric_pFq(ffd_a, d_b, d_z).val_.d_, a_adj); - - // fvar, double, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(ffd_a, d_b, ffd_z).val_.d_, a_adj + z_adj); - - // double, fvar, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(d_a, ffd_b, ffd_z).val_.d_, b_adj + z_adj); - - // double, fvar, double - EXPECT_FLOAT_EQ(hypergeometric_pFq(d_a, ffd_b, d_z).val_.d_, b_adj); - - // double, double, fvar - EXPECT_FLOAT_EQ(hypergeometric_pFq(d_a, d_b, ffd_z).val_.d_, z_adj); -} diff --git a/test/unit/math/rev/fun/hypergeometric_pFq_test.cpp b/test/unit/math/rev/fun/hypergeometric_pFq_test.cpp deleted file mode 100644 index c323401f233..00000000000 --- a/test/unit/math/rev/fun/hypergeometric_pFq_test.cpp +++ /dev/null @@ -1,28 +0,0 @@ -#include -#include -#include -#include - -TEST(RevMath, hypergeometric_pFq) { - using stan::math::hypergeometric_pFq; - using stan::math::var; - using stan::math::vector_d; - using stan::math::vector_v; - - vector_v a(2); - a << 4, 2; - vector_v b(2); - b << 6, 3; - var z = 4; - - var result = hypergeometric_pFq(a, b, z); - result.grad(); - - EXPECT_FLOAT_EQ(3.924636646666071, a[0].adj()); - EXPECT_FLOAT_EQ(6.897245961898751, a[1].adj()); - - EXPECT_FLOAT_EQ(-2.775051002566842, b[0].adj()); - EXPECT_FLOAT_EQ(-4.980095849781222, b[1].adj()); - - EXPECT_FLOAT_EQ(4.916522138006060, z.adj()); -} From a35dc3d5498f0851fe928fc9eac4b5454e3f2d84 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 26 Feb 2024 02:30:13 -0500 Subject: [PATCH 31/36] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/grad_pFq.hpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 01c9a787bd9..f2ac364bcc3 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -118,13 +118,13 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, // Identify the number of iterations to needed for each element to sign // flip from negative to positive - rather than checking at each iteration - Eigen::ArrayXi a_pos_k - = (a_array < 0.0).select(-floor(value_of_rec(a_array)), 0) - .template cast(); + Eigen::ArrayXi a_pos_k = (a_array < 0.0) + .select(-floor(value_of_rec(a_array)), 0) + .template cast(); int all_a_pos_k = a_pos_k.maxCoeff(); - Eigen::ArrayXi b_pos_k - = (b_array < 0.0).select(-floor(value_of_rec(b_array)), 0) - .template cast(); + Eigen::ArrayXi b_pos_k = (b_array < 0.0) + .select(-floor(value_of_rec(b_array)), 0) + .template cast(); int all_b_pos_k = b_pos_k.maxCoeff(); Eigen::ArrayXi a_sign = select(a_pos_k == 0, 1, -1); Eigen::ArrayXi b_sign = select(b_pos_k == 0, 1, -1); @@ -166,8 +166,8 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, // - This is smaller than EPSILON, so the following iteration will // still be 1.0 a_k = (a_k == 1.0 && a_sign == -1) - .select(dbl_min, (a_k < 1.0 && a_sign == -1) - .select(1.0 - a_k, a_k + 1.0 * a_sign)); + .select(dbl_min, (a_k < 1.0 && a_sign == -1) + .select(1.0 - a_k, a_k + 1.0 * a_sign)); a_sign = select(k == a_pos_k - 1, 1, a_sign); } else { a_k += 1.0; @@ -179,8 +179,8 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, if (k < all_a_pos_k) { b_k = (b_k == 1.0 && b_sign == -1) - .select(dbl_min, (b_k < 1.0 && b_sign == -1) - .select(1.0 - b_k, b_k + 1.0 * b_sign)); + .select(dbl_min, (b_k < 1.0 && b_sign == -1) + .select(1.0 - b_k, b_k + 1.0 * b_sign)); b_sign = select(k == b_pos_k - 1, 1, b_sign); } else { b_k += 1.0; @@ -195,7 +195,7 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, } if (CalcZ) { std::get<2>(ret_tuple) - = hypergeometric_pFq(add(a, 1.0), add(b, 1.0), z) * prod(a) / prod(b); + = hypergeometric_pFq(add(a, 1.0), add(b, 1.0), z) * prod(a) / prod(b); } return ret_tuple; } From 83edc2976179906ad32e27ea2d5ed392de698fe5 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 26 Feb 2024 09:51:00 +0200 Subject: [PATCH 32/36] New impl more accurate than 3F2 --- test/unit/math/prim/fun/grad_pFq_test.cpp | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/test/unit/math/prim/fun/grad_pFq_test.cpp b/test/unit/math/prim/fun/grad_pFq_test.cpp index bb206a4cf8c..da124aa1a31 100644 --- a/test/unit/math/prim/fun/grad_pFq_test.cpp +++ b/test/unit/math/prim/fun/grad_pFq_test.cpp @@ -407,9 +407,7 @@ TEST(PrimMath, grad_2F1_negative_z) { } TEST(PrimMath, grad_3F2_cross_zero) { - using stan::math::grad_F32; using stan::math::grad_pFq; - using stan::math::hypergeometric_3F2; using stan::math::hypergeometric_pFq; using stan::math::vector_d; @@ -421,16 +419,14 @@ TEST(PrimMath, grad_3F2_cross_zero) { double z_d = 0.292893; - double g_calc[6]; - grad_F32(g_calc, 1.0, 1.0, -1.0, 2.0, 2.0, 0.292893); - auto pfq_val = hypergeometric_pFq(a_d, b_d, z_d); auto grad_tuple = grad_pFq(pfq_val, a_d, b_d, z_d); - EXPECT_FLOAT_EQ(g_calc[0], std::get<0>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(g_calc[1], std::get<0>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(g_calc[2], std::get<0>(grad_tuple)[2]); - EXPECT_FLOAT_EQ(g_calc[3], std::get<1>(grad_tuple)[0]); - EXPECT_FLOAT_EQ(g_calc[4], std::get<1>(grad_tuple)[1]); - EXPECT_FLOAT_EQ(g_calc[5], std::get<2>(grad_tuple)); + // Values from Mathematica + EXPECT_FLOAT_EQ(-0.0732232500000000, std::get<0>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(-0.0732232500000000, std::get<0>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(0.0681675749282880, std::get<0>(grad_tuple)[2]); + EXPECT_FLOAT_EQ(0.0366116250000000, std::get<1>(grad_tuple)[0]); + EXPECT_FLOAT_EQ(0.0366116250000000, std::get<1>(grad_tuple)[1]); + EXPECT_FLOAT_EQ(-0.250000000000000, std::get<2>(grad_tuple)); } From 8da99c512ab28e8fde79411b5f649c9f19d067da Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 27 Feb 2024 06:32:07 +0200 Subject: [PATCH 33/36] Revert bool camel case --- stan/math/fwd/fun/hypergeometric_pFq.hpp | 15 ++++++++------- stan/math/prim/fun/grad_pFq.hpp | 16 ++++++++-------- stan/math/rev/fun/hypergeometric_pFq.hpp | 14 ++++++++------ 3 files changed, 24 insertions(+), 21 deletions(-) diff --git a/stan/math/fwd/fun/hypergeometric_pFq.hpp b/stan/math/fwd/fun/hypergeometric_pFq.hpp index 907d4ea32d5..50f5ac904fc 100644 --- a/stan/math/fwd/fun/hypergeometric_pFq.hpp +++ b/stan/math/fwd/fun/hypergeometric_pFq.hpp @@ -23,9 +23,9 @@ namespace math { */ template , - bool GradA = !is_constant::value, - bool GradB = !is_constant::value, - bool GradZ = !is_constant::value, + bool grad_a = !is_constant::value, + bool grad_b = !is_constant::value, + bool grad_z = !is_constant::value, require_all_vector_t* = nullptr, require_return_type_t* = nullptr> inline FvarT hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { @@ -39,19 +39,20 @@ inline FvarT hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { auto&& b_val = value_of(b_ref); auto&& z_val = value_of(z); PartialsT pfq_val = hypergeometric_pFq(a_val, b_val, z_val); - auto grad_tuple = grad_pFq(pfq_val, a_val, b_val, z_val); + auto grad_tuple + = grad_pFq(pfq_val, a_val, b_val, z_val); FvarT rtn = FvarT(pfq_val, 0.0); - if (GradA) { + if (grad_a) { rtn.d_ += dot_product(forward_as>(a_ref).d(), std::get<0>(grad_tuple)); } - if (GradB) { + if (grad_b) { rtn.d_ += dot_product(forward_as>(b_ref).d(), std::get<1>(grad_tuple)); } - if (GradZ) { + if (grad_z) { rtn.d_ += forward_as>(z).d_ * std::get<2>(grad_tuple); } diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index f2ac364bcc3..7ecd20e5819 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -69,9 +69,9 @@ namespace math { * {k!\prod_{j=1}^q\left(b_j\right)_k}} * \f$ * - * @tparam CalcA Boolean for whether to calculate derivatives wrt to 'a' - * @tparam CalcB Boolean for whether to calculate derivatives wrt to 'b' - * @tparam CalcZ Boolean for whether to calculate derivatives wrt to 'z' + * @tparam calc_a Boolean for whether to calculate derivatives wrt to 'a' + * @tparam calc_b Boolean for whether to calculate derivatives wrt to 'b' + * @tparam calc_z Boolean for whether to calculate derivatives wrt to 'z' * @tparam TpFq Scalar type of hypergeometric_pFq return value * @tparam Ta Eigen type with either one row or column at compile time * @tparam Tb Eigen type with either one row or column at compile time @@ -84,7 +84,7 @@ namespace math { * @param[in] max_steps Maximum number of iterations * @return Tuple of gradients */ -template , typename Ta_Rtn = promote_scalar_t>, @@ -102,7 +102,7 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, std::tuple ret_tuple; - if (CalcA || CalcB) { + if (calc_a || calc_b) { std::get<0>(ret_tuple).setConstant(a.size(), -pfq_val); std::get<1>(ret_tuple).setConstant(b.size(), pfq_val); Eigen::Array a_grad(a.size()); @@ -138,7 +138,7 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, T_Rtn log_base = 0; while ((k < 10 || curr_log_prec > log(precision)) && (k <= max_steps)) { curr_log_prec = NEGATIVE_INFTY; - if (CalcA) { + if (calc_a) { a_grad = log(abs(digamma_a)) + log_base; std::get<0>(ret_tuple).array() += exp(a_grad) * base_sign * sign(value_of_rec(digamma_a)); @@ -147,7 +147,7 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, digamma_a += inv(a_k) * a_sign; } - if (CalcB) { + if (calc_b) { b_grad = log(abs(digamma_b)) + log_base; std::get<1>(ret_tuple).array() -= exp(b_grad) * base_sign * sign(value_of_rec(digamma_b)); @@ -193,7 +193,7 @@ std::tuple grad_pFq(const TpFq& pfq_val, const Ta& a, k += 1; } } - if (CalcZ) { + if (calc_z) { std::get<2>(ret_tuple) = hypergeometric_pFq(add(a, 1.0), add(b, 1.0), z) * prod(a) / prod(b); } diff --git a/stan/math/rev/fun/hypergeometric_pFq.hpp b/stan/math/rev/fun/hypergeometric_pFq.hpp index b0a55e2e11d..81be176534a 100644 --- a/stan/math/rev/fun/hypergeometric_pFq.hpp +++ b/stan/math/rev/fun/hypergeometric_pFq.hpp @@ -22,8 +22,10 @@ namespace math { * @return Generalized hypergeometric function */ template < - typename Ta, typename Tb, typename Tz, bool GradA = !is_constant::value, - bool GradB = !is_constant::value, bool GradZ = !is_constant::value, + typename Ta, typename Tb, typename Tz, + bool grad_a = !is_constant::value, + bool grad_b = !is_constant::value, + bool grad_z = !is_constant::value, require_all_matrix_t* = nullptr, require_return_type_t* = nullptr> inline var hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { @@ -32,17 +34,17 @@ inline var hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { auto pfq_val = hypergeometric_pFq(a.val(), b.val(), value_of(z)); return make_callback_var( pfq_val, [arena_a, arena_b, z, pfq_val](auto& vi) mutable { - auto grad_tuple = grad_pFq( + auto grad_tuple = grad_pFq( pfq_val, arena_a.val(), arena_b.val(), value_of(z)); - if (GradA) { + if (grad_a) { forward_as>(arena_a).adj() += vi.adj() * std::get<0>(grad_tuple); } - if (GradB) { + if (grad_b) { forward_as>(arena_b).adj() += vi.adj() * std::get<1>(grad_tuple); } - if (GradZ) { + if (grad_z) { forward_as>(z).adj() += vi.adj() * std::get<2>(grad_tuple); } From 1f1ff16335496790fd1720a0e6df8f80202985b1 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 26 Feb 2024 23:39:27 -0500 Subject: [PATCH 34/36] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/fwd/fun/hypergeometric_pFq.hpp | 2 +- stan/math/rev/fun/hypergeometric_pFq.hpp | 13 ++++++------- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/stan/math/fwd/fun/hypergeometric_pFq.hpp b/stan/math/fwd/fun/hypergeometric_pFq.hpp index 50f5ac904fc..7418ac8186e 100644 --- a/stan/math/fwd/fun/hypergeometric_pFq.hpp +++ b/stan/math/fwd/fun/hypergeometric_pFq.hpp @@ -40,7 +40,7 @@ inline FvarT hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { auto&& z_val = value_of(z); PartialsT pfq_val = hypergeometric_pFq(a_val, b_val, z_val); auto grad_tuple - = grad_pFq(pfq_val, a_val, b_val, z_val); + = grad_pFq(pfq_val, a_val, b_val, z_val); FvarT rtn = FvarT(pfq_val, 0.0); diff --git a/stan/math/rev/fun/hypergeometric_pFq.hpp b/stan/math/rev/fun/hypergeometric_pFq.hpp index 81be176534a..15c5616e17e 100644 --- a/stan/math/rev/fun/hypergeometric_pFq.hpp +++ b/stan/math/rev/fun/hypergeometric_pFq.hpp @@ -21,13 +21,12 @@ namespace math { * @param[in] z Scalar z argument * @return Generalized hypergeometric function */ -template < - typename Ta, typename Tb, typename Tz, - bool grad_a = !is_constant::value, - bool grad_b = !is_constant::value, - bool grad_z = !is_constant::value, - require_all_matrix_t* = nullptr, - require_return_type_t* = nullptr> +template ::value, + bool grad_b = !is_constant::value, + bool grad_z = !is_constant::value, + require_all_matrix_t* = nullptr, + require_return_type_t* = nullptr> inline var hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { arena_t arena_a = a; arena_t arena_b = b; From 8919824d8fb873209e2734400eee1f3bb7168f28 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 28 Feb 2024 22:29:46 +0200 Subject: [PATCH 35/36] Trigger CI From ecc713bd44ad1fbcee94aa7396c21466d59489a4 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 12 Apr 2024 11:24:13 +0300 Subject: [PATCH 36/36] Apply review suggestions --- stan/math/fwd/fun/hypergeometric_pFq.hpp | 2 +- stan/math/prim/fun/grad_pFq.hpp | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/stan/math/fwd/fun/hypergeometric_pFq.hpp b/stan/math/fwd/fun/hypergeometric_pFq.hpp index 7418ac8186e..cec91aaa90b 100644 --- a/stan/math/fwd/fun/hypergeometric_pFq.hpp +++ b/stan/math/fwd/fun/hypergeometric_pFq.hpp @@ -27,7 +27,7 @@ template ::value, bool grad_z = !is_constant::value, require_all_vector_t* = nullptr, - require_return_type_t* = nullptr> + require_fvar_t* = nullptr> inline FvarT hypergeometric_pFq(const Ta& a, const Tb& b, const Tz& z) { using PartialsT = partials_type_t; using ARefT = ref_type_t; diff --git a/stan/math/prim/fun/grad_pFq.hpp b/stan/math/prim/fun/grad_pFq.hpp index 7ecd20e5819..2275662a443 100644 --- a/stan/math/prim/fun/grad_pFq.hpp +++ b/stan/math/prim/fun/grad_pFq.hpp @@ -30,11 +30,11 @@ namespace math { * \frac * {\psi\left(k+a_1\right)\left(\prod_{j=1}^p\left(a_j\right)_k\right)z^k} * {k!\prod_{j=1}^q\left(b_j\right)_k}} - * - \psi\left(a_1\right)_pF_q(a_1,...,a_p;b_1,...,b_q;z) + * - \psi\left(a_1\right){}_pF_q(a_1,...,a_p;b_1,...,b_q;z) * \f$ * \f$ * \frac{\partial }{\partial b_1} = - * \psi\left(b_1\right)_pF_q(a_1,...,a_p;b_1,...,b_q;z) - + * \psi\left(b_1\right){}_pF_q(a_1,...,a_p;b_1,...,b_q;z) - * \sum_{k=0}^{\infty}{ * \frac * {\psi\left(k+b_1\right)\left(\prod_{j=1}^p\left(a_j\right)_k\right)z^k} @@ -44,24 +44,24 @@ namespace math { * \f$ * \frac{\partial }{\partial z} = * \frac{\prod_{j=1}^{p}(a_j)}{\prod_{j=1}^{q} (b_j)}\ - * * _pF_q(a_1+1,...,a_p+1;b_1+1,...,b_q+1;z) + * {}_pF_q(a_1+1,...,a_p+1;b_1+1,...,b_q+1;z) * \f$ * * Noting the the recurrence relation for the digamma function: * \f$ \psi(x + 1) = \psi(x) + \frac{1}{x} \f$, the gradients for the - * function w.r.t a & b then simplify to: + * function with respect to a and b then simplify to: * \f$ * \frac{\partial }{\partial a_1} = * \sum_{k=1}^{\infty}{ * \frac * {\left(1 + \sum_{m=0}^{k-1}\frac{1}{m+a_1}\right) - * * (\prod_{j=1}^p\left(a_j\right)_k\right)z^k} + * * \left(\prod_{j=1}^p\left(a_j\right)_k\right)z^k} * {k!\prod_{j=1}^q\left(b_j\right)_k}} - * - _pF_q(a_1,...,a_p;b_1,...,b_q;z) + * - {}_pF_q(a_1,...,a_p;b_1,...,b_q;z) * \f$ * \f$ * \frac{\partial }{\partial b_1} = - * _pF_q(a_1,...,a_p;b_1,...,b_q;z) - + * {}_pF_q(a_1,...,a_p;b_1,...,b_q;z) - * \sum_{k=1}^{\infty}{ * \frac * {\left(1 + \sum_{m=0}^{k-1}\frac{1}{m+b_1}\right)