From 2b8ae140373ce37a0c6086324ae569c5277f8651 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 8 Jul 2022 06:11:07 +0300 Subject: [PATCH 01/30] Numerically stable bernoulli cdf functions --- stan/math/prim/prob/bernoulli_cdf.hpp | 17 +++++++----- stan/math/prim/prob/bernoulli_lccdf.hpp | 9 ++++--- stan/math/prim/prob/bernoulli_lcdf.hpp | 10 +++---- .../unit/math/mix/prob/bernoulli_cdf_test.cpp | 27 +++++++++++++++++++ .../math/mix/prob/bernoulli_lccdf_test.cpp | 27 +++++++++++++++++++ .../math/mix/prob/bernoulli_lcdf_test.cpp | 27 +++++++++++++++++++ 6 files changed, 101 insertions(+), 16 deletions(-) create mode 100644 test/unit/math/mix/prob/bernoulli_cdf_test.cpp create mode 100644 test/unit/math/mix/prob/bernoulli_lccdf_test.cpp create mode 100644 test/unit/math/mix/prob/bernoulli_lcdf_test.cpp diff --git a/stan/math/prim/prob/bernoulli_cdf.hpp b/stan/math/prim/prob/bernoulli_cdf.hpp index 3e0154bb22d..63d25fc8a68 100644 --- a/stan/math/prim/prob/bernoulli_cdf.hpp +++ b/stan/math/prim/prob/bernoulli_cdf.hpp @@ -47,7 +47,10 @@ return_type_t bernoulli_cdf(const T_n& n, const T_prob& theta) { operands_and_partials ops_partials(theta_ref); scalar_seq_view n_vec(n); - scalar_seq_view theta_vec(theta_ref); + const auto& log1m_theta = to_ref(log1m(theta_ref)); + scalar_seq_view theta_vec(log1m_theta); + scalar_seq_view partials_vec(exp(-log1m_theta)); + size_t max_size_seq_view = max_size(n, theta); // Explicit return for extreme values @@ -65,21 +68,21 @@ return_type_t bernoulli_cdf(const T_n& n, const T_prob& theta) { continue; } - const T_partials_return Pi = 1 - theta_vec.val(i); - - P *= Pi; + P += theta_vec.val(i); if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] += -1 / Pi; + ops_partials.edge1_.partials_[i] -= partials_vec.val(i); } } + const auto& exp_P = to_ref(exp(P)); + if (!is_constant_all::value) { for (size_t i = 0; i < stan::math::size(theta); ++i) { - ops_partials.edge1_.partials_[i] *= P; + ops_partials.edge1_.partials_[i] *= exp_P; } } - return ops_partials.build(P); + return ops_partials.build(exp_P); } } // namespace math diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index e63119bbd69..5425bfcf8dc 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -51,7 +51,9 @@ return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { operands_and_partials ops_partials(theta_ref); scalar_seq_view n_vec(n); - scalar_seq_view theta_vec(theta_ref); + const auto& log1m_theta = to_ref(log1m(theta_ref)); + scalar_seq_view theta_vec(log1m_theta); + scalar_seq_view partials_vec(exp(-log1m_theta)); size_t max_size_seq_view = max_size(n, theta); // Explicit return for extreme values @@ -67,12 +69,11 @@ return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { } for (size_t i = 0; i < max_size_seq_view; i++) { - const T_partials_return Pi = theta_vec.val(i); - P += log(Pi); + P += theta_vec.val(i); if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] += inv(Pi); + ops_partials.edge1_.partials_[i] -= partials_vec.val(i); } } diff --git a/stan/math/prim/prob/bernoulli_lcdf.hpp b/stan/math/prim/prob/bernoulli_lcdf.hpp index d9a40898ff9..4fd0c3c9710 100644 --- a/stan/math/prim/prob/bernoulli_lcdf.hpp +++ b/stan/math/prim/prob/bernoulli_lcdf.hpp @@ -51,7 +51,9 @@ return_type_t bernoulli_lcdf(const T_n& n, const T_prob& theta) { operands_and_partials ops_partials(theta_ref); scalar_seq_view n_vec(n); - scalar_seq_view theta_vec(theta_ref); + const auto& log1m_theta = to_ref(log1m(theta_ref)); + scalar_seq_view theta_vec(log1m_theta); + scalar_seq_view partials_vec(exp(-log1m_theta)); size_t max_size_seq_view = max_size(n, theta); // Explicit return for extreme values @@ -69,12 +71,10 @@ return_type_t bernoulli_lcdf(const T_n& n, const T_prob& theta) { continue; } - const T_partials_return Pi = 1 - theta_vec.val(i); - - P += log(Pi); + P += theta_vec.val(i); if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] -= inv(Pi); + ops_partials.edge1_.partials_[i] -= partials_vec.val(i); } } diff --git a/test/unit/math/mix/prob/bernoulli_cdf_test.cpp b/test/unit/math/mix/prob/bernoulli_cdf_test.cpp new file mode 100644 index 00000000000..3fc7213ae85 --- /dev/null +++ b/test/unit/math/mix/prob/bernoulli_cdf_test.cpp @@ -0,0 +1,27 @@ +#include +#include + +TEST(mathMixScalFun, bernoulliCDF) { + // bind integer arg because can't autodiff through + auto f = [](const auto& x1) { + return + [=](const auto& x2) { return stan::math::bernoulli_cdf(x1, x2); }; + }; + stan::test::expect_ad(f(0), 0.1); + stan::test::expect_ad(f(0), std::numeric_limits::quiet_NaN()); + stan::test::expect_ad(f(1), 0.5); + stan::test::expect_ad(f(1), std::numeric_limits::quiet_NaN()); + stan::test::expect_ad(f(1), 0.2); + + + std::vector std_in1{0, 1}; + Eigen::VectorXd in2(2); + in2 << 0.5, 0.9; + + + stan::test::expect_ad(f(std_in1), 0.2); + stan::test::expect_ad(f(std_in1), std::numeric_limits::quiet_NaN()); + stan::test::expect_ad(f(1), in2); + stan::test::expect_ad(f(std_in1), in2); + stan::test::expect_ad(f(std_in1), std::numeric_limits::quiet_NaN()); +} diff --git a/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp b/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp new file mode 100644 index 00000000000..d98136e21ce --- /dev/null +++ b/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp @@ -0,0 +1,27 @@ +#include +#include + +TEST(mathMixScalFun, bernoulliLCCDF) { + // bind integer arg because can't autodiff through + auto f = [](const auto& x1) { + return + [=](const auto& x2) { return stan::math::bernoulli_lccdf(x1, x2); }; + }; + stan::test::expect_ad(f(0), 0.1); + stan::test::expect_ad(f(0), std::numeric_limits::quiet_NaN()); + stan::test::expect_ad(f(1), 0.5); + stan::test::expect_ad(f(1), std::numeric_limits::quiet_NaN()); + stan::test::expect_ad(f(1), 0.2); + + + std::vector std_in1{0, 1}; + Eigen::VectorXd in2(2); + in2 << 0.5, 0.9; + + + stan::test::expect_ad(f(std_in1), 0.2); + stan::test::expect_ad(f(std_in1), std::numeric_limits::quiet_NaN()); + stan::test::expect_ad(f(1), in2); + stan::test::expect_ad(f(std_in1), in2); + stan::test::expect_ad(f(std_in1), std::numeric_limits::quiet_NaN()); +} diff --git a/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp b/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp new file mode 100644 index 00000000000..652a752ff4c --- /dev/null +++ b/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp @@ -0,0 +1,27 @@ +#include +#include + +TEST(mathMixScalFun, bernoulliLCDF) { + // bind integer arg because can't autodiff through + auto f = [](const auto& x1) { + return + [=](const auto& x2) { return stan::math::bernoulli_lcdf(x1, x2); }; + }; + stan::test::expect_ad(f(0), 0.1); + stan::test::expect_ad(f(0), std::numeric_limits::quiet_NaN()); + stan::test::expect_ad(f(1), 0.5); + stan::test::expect_ad(f(1), std::numeric_limits::quiet_NaN()); + stan::test::expect_ad(f(1), 0.2); + + + std::vector std_in1{0, 1}; + Eigen::VectorXd in2(2); + in2 << 0.5, 0.9; + + + stan::test::expect_ad(f(std_in1), 0.2); + stan::test::expect_ad(f(std_in1), std::numeric_limits::quiet_NaN()); + stan::test::expect_ad(f(1), in2); + stan::test::expect_ad(f(std_in1), in2); + stan::test::expect_ad(f(std_in1), std::numeric_limits::quiet_NaN()); +} From c2022be347bdf64c0b79d6b9e83bb59af3a22b85 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 8 Jul 2022 06:33:15 +0300 Subject: [PATCH 02/30] CDF initial value --- stan/math/prim/prob/bernoulli_cdf.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/prob/bernoulli_cdf.hpp b/stan/math/prim/prob/bernoulli_cdf.hpp index 63d25fc8a68..259e7b44a5d 100644 --- a/stan/math/prim/prob/bernoulli_cdf.hpp +++ b/stan/math/prim/prob/bernoulli_cdf.hpp @@ -43,7 +43,7 @@ return_type_t bernoulli_cdf(const T_n& n, const T_prob& theta) { return 1.0; } - T_partials_return P(1.0); + T_partials_return P(0.0); operands_and_partials ops_partials(theta_ref); scalar_seq_view n_vec(n); @@ -75,7 +75,7 @@ return_type_t bernoulli_cdf(const T_n& n, const T_prob& theta) { } } - const auto& exp_P = to_ref(exp(P)); + plain_type_t exp_P = exp(P); if (!is_constant_all::value) { for (size_t i = 0; i < stan::math::size(theta); ++i) { From a9d0800dfe79342042d470c1ccdcdf9019a9af0a Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 8 Jul 2022 06:56:45 +0300 Subject: [PATCH 03/30] cpplint --- stan/math/prim/prob/bernoulli_lccdf.hpp | 1 - test/unit/math/mix/prob/bernoulli_cdf_test.cpp | 2 -- test/unit/math/mix/prob/bernoulli_lccdf_test.cpp | 2 -- test/unit/math/mix/prob/bernoulli_lcdf_test.cpp | 2 -- 4 files changed, 7 deletions(-) diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index 5425bfcf8dc..c1a36a64ae3 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -69,7 +69,6 @@ return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { } for (size_t i = 0; i < max_size_seq_view; i++) { - P += theta_vec.val(i); if (!is_constant_all::value) { diff --git a/test/unit/math/mix/prob/bernoulli_cdf_test.cpp b/test/unit/math/mix/prob/bernoulli_cdf_test.cpp index 3fc7213ae85..7212f2aff9e 100644 --- a/test/unit/math/mix/prob/bernoulli_cdf_test.cpp +++ b/test/unit/math/mix/prob/bernoulli_cdf_test.cpp @@ -13,12 +13,10 @@ TEST(mathMixScalFun, bernoulliCDF) { stan::test::expect_ad(f(1), std::numeric_limits::quiet_NaN()); stan::test::expect_ad(f(1), 0.2); - std::vector std_in1{0, 1}; Eigen::VectorXd in2(2); in2 << 0.5, 0.9; - stan::test::expect_ad(f(std_in1), 0.2); stan::test::expect_ad(f(std_in1), std::numeric_limits::quiet_NaN()); stan::test::expect_ad(f(1), in2); diff --git a/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp b/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp index d98136e21ce..0e0c56726f0 100644 --- a/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp +++ b/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp @@ -13,12 +13,10 @@ TEST(mathMixScalFun, bernoulliLCCDF) { stan::test::expect_ad(f(1), std::numeric_limits::quiet_NaN()); stan::test::expect_ad(f(1), 0.2); - std::vector std_in1{0, 1}; Eigen::VectorXd in2(2); in2 << 0.5, 0.9; - stan::test::expect_ad(f(std_in1), 0.2); stan::test::expect_ad(f(std_in1), std::numeric_limits::quiet_NaN()); stan::test::expect_ad(f(1), in2); diff --git a/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp b/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp index 652a752ff4c..7e24e9e130f 100644 --- a/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp +++ b/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp @@ -13,12 +13,10 @@ TEST(mathMixScalFun, bernoulliLCDF) { stan::test::expect_ad(f(1), std::numeric_limits::quiet_NaN()); stan::test::expect_ad(f(1), 0.2); - std::vector std_in1{0, 1}; Eigen::VectorXd in2(2); in2 << 0.5, 0.9; - stan::test::expect_ad(f(std_in1), 0.2); stan::test::expect_ad(f(std_in1), std::numeric_limits::quiet_NaN()); stan::test::expect_ad(f(1), in2); From bcfaec1de1679a9cfc216bac96f9caed12222716 Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Fri, 8 Jul 2022 00:01:04 -0400 Subject: [PATCH 04/30] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- test/unit/math/mix/prob/bernoulli_cdf_test.cpp | 3 +-- test/unit/math/mix/prob/bernoulli_lccdf_test.cpp | 3 +-- test/unit/math/mix/prob/bernoulli_lcdf_test.cpp | 3 +-- 3 files changed, 3 insertions(+), 6 deletions(-) diff --git a/test/unit/math/mix/prob/bernoulli_cdf_test.cpp b/test/unit/math/mix/prob/bernoulli_cdf_test.cpp index 7212f2aff9e..797b9a78880 100644 --- a/test/unit/math/mix/prob/bernoulli_cdf_test.cpp +++ b/test/unit/math/mix/prob/bernoulli_cdf_test.cpp @@ -4,8 +4,7 @@ TEST(mathMixScalFun, bernoulliCDF) { // bind integer arg because can't autodiff through auto f = [](const auto& x1) { - return - [=](const auto& x2) { return stan::math::bernoulli_cdf(x1, x2); }; + return [=](const auto& x2) { return stan::math::bernoulli_cdf(x1, x2); }; }; stan::test::expect_ad(f(0), 0.1); stan::test::expect_ad(f(0), std::numeric_limits::quiet_NaN()); diff --git a/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp b/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp index 0e0c56726f0..d8cb10398c4 100644 --- a/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp +++ b/test/unit/math/mix/prob/bernoulli_lccdf_test.cpp @@ -4,8 +4,7 @@ TEST(mathMixScalFun, bernoulliLCCDF) { // bind integer arg because can't autodiff through auto f = [](const auto& x1) { - return - [=](const auto& x2) { return stan::math::bernoulli_lccdf(x1, x2); }; + return [=](const auto& x2) { return stan::math::bernoulli_lccdf(x1, x2); }; }; stan::test::expect_ad(f(0), 0.1); stan::test::expect_ad(f(0), std::numeric_limits::quiet_NaN()); diff --git a/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp b/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp index 7e24e9e130f..8a8e9c778e9 100644 --- a/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp +++ b/test/unit/math/mix/prob/bernoulli_lcdf_test.cpp @@ -4,8 +4,7 @@ TEST(mathMixScalFun, bernoulliLCDF) { // bind integer arg because can't autodiff through auto f = [](const auto& x1) { - return - [=](const auto& x2) { return stan::math::bernoulli_lcdf(x1, x2); }; + return [=](const auto& x2) { return stan::math::bernoulli_lcdf(x1, x2); }; }; stan::test::expect_ad(f(0), 0.1); stan::test::expect_ad(f(0), std::numeric_limits::quiet_NaN()); From 143f707aa264c42da5c209c862a10bde82ba3ea0 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 10 Jul 2022 19:40:28 +0300 Subject: [PATCH 05/30] Add select function, vectorise cdf --- stan/math/opencl/kernel_generator/select.hpp | 16 ---- stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/select.hpp | 83 ++++++++++++++++++++ stan/math/prim/prob/bernoulli_cdf.hpp | 41 +++------- 4 files changed, 94 insertions(+), 47 deletions(-) create mode 100644 stan/math/prim/fun/select.hpp diff --git a/stan/math/opencl/kernel_generator/select.hpp b/stan/math/opencl/kernel_generator/select.hpp index 4894b1d8c0a..c2e6717ce91 100644 --- a/stan/math/opencl/kernel_generator/select.hpp +++ b/stan/math/opencl/kernel_generator/select.hpp @@ -150,22 +150,6 @@ select(T_condition&& condition, T_then&& then, T_else&& els) { // NOLINT as_operation_cl(std::forward(els))}; } -/** - * Scalar overload of the selection operation. - * @tparam T_then type of then scalar - * @tparam T_else type of else scalar - * @param condition condition - * @param then then result - * @param els else result - * @return `condition ? then : els` - */ -template * = nullptr> -inline std::common_type_t select(bool condition, T_then then, - T_else els) { - return condition ? then : els; -} - /** @}*/ } // namespace math } // namespace stan diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 4afa7914656..104496b1227 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -301,6 +301,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/select.hpp b/stan/math/prim/fun/select.hpp new file mode 100644 index 00000000000..e71fcf4f4fe --- /dev/null +++ b/stan/math/prim/fun/select.hpp @@ -0,0 +1,83 @@ +#ifndef STAN_MATH_PRIM_FUN_SELECT_HPP +#define STAN_MATH_PRIM_FUN_SELECT_HPP + +#include + +namespace stan { +namespace math { + +/** + * Return the second argument if the first argument is true + * and otherwise return the second argument. + * + *

This is just a convenience method to provide a function + * with the same behavior as the built-in ternary operator. + * In general, this function behaves as if defined by + * + *

select(c, y1, y0) = c ? y1 : y0. + * + * @tparam T_true type of the true argument + * @tparam T_false type of the false argument + * @param c Boolean condition value. + * @param y_true Value to return if condition is true. + * @param y_false Value to return if condition is false. + */ +template * = nullptr> +inline auto select(const bool c, const T_true y_true, const T_false y_false) { + return c ? y_true : y_false; +} + +template * = nullptr> +inline auto select(const bool c, const T_true y_true, const T_false y_false) { + return c ? y_true : y_false; +} + +template * = nullptr, + require_stan_scalar_t* = nullptr> +inline plain_type_t select(const bool c, const T_true& y_true, const T_false& y_false) { + if (c) { + return y_true; + } + + return y_true.unaryExpr([&](auto&& y){ return y_false; }); +} + + +template * = nullptr, + require_eigen_t* = nullptr> +inline plain_type_t select(const bool c, const T_true y_true, const T_false y_false) { + if (c) { + return y_false.unaryExpr([&](auto&& y){ return y_true; }); + } + + return y_false; +} + +template * = nullptr, + require_all_stan_scalar_t* = nullptr> +inline auto select(const T_bool c, const T_true y_true, + const T_false y_false) { + return c.unaryExpr([&](bool cond){ + return cond ? y_true : y_false; + }).eval(); +} + +template * = nullptr, + require_any_eigen_array_t* = nullptr> +inline auto select(const T_bool c, + const T_true y_true, + const T_false y_false) { + return c.select(y_true, y_false).eval(); +} + + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/prim/prob/bernoulli_cdf.hpp b/stan/math/prim/prob/bernoulli_cdf.hpp index 259e7b44a5d..c5cc9793ff9 100644 --- a/stan/math/prim/prob/bernoulli_cdf.hpp +++ b/stan/math/prim/prob/bernoulli_cdf.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -36,6 +37,7 @@ return_type_t bernoulli_cdf(const T_n& n, const T_prob& theta) { check_consistent_sizes(function, "Random variable", n, "Probability parameter", theta); T_theta_ref theta_ref = theta; + const auto& n_arr = as_array_or_scalar(n); check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0, 1.0); @@ -43,46 +45,23 @@ return_type_t bernoulli_cdf(const T_n& n, const T_prob& theta) { return 1.0; } - T_partials_return P(0.0); operands_and_partials ops_partials(theta_ref); - scalar_seq_view n_vec(n); - const auto& log1m_theta = to_ref(log1m(theta_ref)); - scalar_seq_view theta_vec(log1m_theta); - scalar_seq_view partials_vec(exp(-log1m_theta)); - - size_t max_size_seq_view = max_size(n, theta); - // Explicit return for extreme values // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(n); i++) { - if (n_vec.val(i) < 0) { - return ops_partials.build(0.0); - } - } - - for (size_t i = 0; i < max_size_seq_view; i++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (n_vec.val(i) >= 1) { - continue; - } - - P += theta_vec.val(i); - - if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] -= partials_vec.val(i); - } + if (sum(n_arr < 0)) { + return ops_partials.build(0.0); } + const auto& theta_arr = as_value_column_array_or_scalar(theta_ref); + const auto& log1m_theta = select(theta_arr == 1, 0.0, log1m(theta_arr)); + const auto& P1 = select(n_arr == 0, log1m_theta, 0.0); - plain_type_t exp_P = exp(P); + T_partials_return P = sum(P1); if (!is_constant_all::value) { - for (size_t i = 0; i < stan::math::size(theta); ++i) { - ops_partials.edge1_.partials_[i] *= exp_P; - } + ops_partials.edge1_.partials_ = select(n_arr == 0, -exp(P - P1), 0.0); } - return ops_partials.build(exp_P); + return ops_partials.build(exp(P)); } } // namespace math From 62d8640e2387a2116483b762037bed63992dd987 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 10 Jul 2022 20:21:17 +0300 Subject: [PATCH 06/30] Vectorise lcdf and lccdf --- stan/math/prim/fun/select.hpp | 10 ++++++-- stan/math/prim/prob/bernoulli_lccdf.hpp | 31 ++++++++----------------- stan/math/prim/prob/bernoulli_lcdf.hpp | 31 +++++++------------------ 3 files changed, 26 insertions(+), 46 deletions(-) diff --git a/stan/math/prim/fun/select.hpp b/stan/math/prim/fun/select.hpp index e71fcf4f4fe..f23af49bdf1 100644 --- a/stan/math/prim/fun/select.hpp +++ b/stan/math/prim/fun/select.hpp @@ -31,13 +31,19 @@ inline auto select(const bool c, const T_true y_true, const T_false y_false) { template * = nullptr> inline auto select(const bool c, const T_true y_true, const T_false y_false) { - return c ? y_true : y_false; + using return_scalar_t = return_type_t; + if (c) { + return y_true.template cast().eval(); + } + + return y_false.template cast().eval(); } template * = nullptr, require_stan_scalar_t* = nullptr> -inline plain_type_t select(const bool c, const T_true& y_true, const T_false& y_false) { +inline plain_type_t select(const bool c, const T_true& y_true, + const T_false& y_false) { if (c) { return y_true; } diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index c1a36a64ae3..b54e813f511 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -40,6 +40,7 @@ return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { check_consistent_sizes(function, "Random variable", n, "Probability parameter", theta); T_theta_ref theta_ref = theta; + const auto& n_arr = as_array_or_scalar(n); check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0, 1.0); @@ -47,36 +48,24 @@ return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { return 0.0; } - T_partials_return P(0.0); operands_and_partials ops_partials(theta_ref); - scalar_seq_view n_vec(n); - const auto& log1m_theta = to_ref(log1m(theta_ref)); - scalar_seq_view theta_vec(log1m_theta); - scalar_seq_view partials_vec(exp(-log1m_theta)); - size_t max_size_seq_view = max_size(n, theta); - // Explicit return for extreme values // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(n); i++) { - const double n_dbl = n_vec.val(i); - if (n_dbl < 0) { - return ops_partials.build(0.0); - } - if (n_dbl >= 1) { - return ops_partials.build(NEGATIVE_INFTY); - } + if (sum(n_arr < 0)) { + return ops_partials.build(0.0); + } + if (sum(n_arr >= 1)) { + return ops_partials.build(NEGATIVE_INFTY); } - for (size_t i = 0; i < max_size_seq_view; i++) { - P += theta_vec.val(i); + const auto& theta_arr = as_value_column_array_or_scalar(theta_ref); - if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] -= partials_vec.val(i); - } + if (!is_constant_all::value) { + ops_partials.edge1_.partials_ = select(true, inv(theta_arr), n_arr); } - return ops_partials.build(P); + return ops_partials.build(sum(select(true, log(theta_arr), n_arr))); } } // namespace math diff --git a/stan/math/prim/prob/bernoulli_lcdf.hpp b/stan/math/prim/prob/bernoulli_lcdf.hpp index 4fd0c3c9710..800a94b4a7d 100644 --- a/stan/math/prim/prob/bernoulli_lcdf.hpp +++ b/stan/math/prim/prob/bernoulli_lcdf.hpp @@ -40,6 +40,7 @@ return_type_t bernoulli_lcdf(const T_n& n, const T_prob& theta) { check_consistent_sizes(function, "Random variable", n, "Probability parameter", theta); T_theta_ref theta_ref = theta; + const auto& n_arr = as_array_or_scalar(n); check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0, 1.0); @@ -47,38 +48,22 @@ return_type_t bernoulli_lcdf(const T_n& n, const T_prob& theta) { return 0.0; } - T_partials_return P(0.0); operands_and_partials ops_partials(theta_ref); - scalar_seq_view n_vec(n); - const auto& log1m_theta = to_ref(log1m(theta_ref)); - scalar_seq_view theta_vec(log1m_theta); - scalar_seq_view partials_vec(exp(-log1m_theta)); - size_t max_size_seq_view = max_size(n, theta); - // Explicit return for extreme values // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < stan::math::size(n); i++) { - if (n_vec.val(i) < 0) { - return ops_partials.build(NEGATIVE_INFTY); - } + if (sum(n_arr < 0)) { + return ops_partials.build(NEGATIVE_INFTY); } - for (size_t i = 0; i < max_size_seq_view; i++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (n_vec.val(i) >= 1) { - continue; - } - - P += theta_vec.val(i); + const auto& theta_arr = as_value_column_array_or_scalar(theta_ref); + const auto& log1m_theta = select(theta_arr == 1, 0.0, log1m(theta_arr)); - if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] -= partials_vec.val(i); - } + if (!is_constant_all::value) { + ops_partials.edge1_.partials_ = select(n_arr == 0, -exp(-log1m_theta), 0.0); } - return ops_partials.build(P); + return ops_partials.build(sum(select(n_arr == 0, log1m_theta, 0.0))); } } // namespace math From 0edeb5fe326eb5253143cf1193f58f1c156139bc Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 10 Jul 2022 20:51:18 +0300 Subject: [PATCH 07/30] Update doc --- stan/math/prim/fun/select.hpp | 83 ++++++++++++++++++++++++++++------- 1 file changed, 68 insertions(+), 15 deletions(-) diff --git a/stan/math/prim/fun/select.hpp b/stan/math/prim/fun/select.hpp index f23af49bdf1..a84c2a4f160 100644 --- a/stan/math/prim/fun/select.hpp +++ b/stan/math/prim/fun/select.hpp @@ -8,13 +8,9 @@ namespace math { /** * Return the second argument if the first argument is true - * and otherwise return the second argument. + * and otherwise return the third argument. * - *

This is just a convenience method to provide a function - * with the same behavior as the built-in ternary operator. - * In general, this function behaves as if defined by - * - *

select(c, y1, y0) = c ? y1 : y0. + * select(c, y1, y0) = c ? y1 : y0. * * @tparam T_true type of the true argument * @tparam T_false type of the false argument @@ -28,17 +24,37 @@ inline auto select(const bool c, const T_true y_true, const T_false y_false) { return c ? y_true : y_false; } +/** + * Return the second argument if the first argument is true + * and otherwise return the third argument. Overload for use with two Eigen + * objects. + * + * @tparam T_true type of the true argument + * @tparam T_false type of the false argument + * @param c Boolean condition value. + * @param y_true Value to return if condition is true. + * @param y_false Value to return if condition is false. + */ template * = nullptr> inline auto select(const bool c, const T_true y_true, const T_false y_false) { - using return_scalar_t = return_type_t; - if (c) { - return y_true.template cast().eval(); - } - - return y_false.template cast().eval(); + return y_true.binaryExpr(y_false, [&](auto&& x, auto&& y) { + return c ? x : y; + }).eval(); } +/** + * Return the second Eigen argument if the first argument is true + * and otherwise return the second Eigen argument. Overload for use with one + * scalar and one Eigen object. If chosen, the scalar is returned as an Eigen + * object of the same size and type as the provided argument. + * + * @tparam T_true type of the true argument + * @tparam T_false type of the false argument + * @param c Boolean condition value. + * @param y_true Value to return if condition is true. + * @param y_false Value to return if condition is false. + */ template * = nullptr, require_stan_scalar_t* = nullptr> @@ -51,18 +67,43 @@ inline plain_type_t select(const bool c, const T_true& y_true, return y_true.unaryExpr([&](auto&& y){ return y_false; }); } - +/** + * Return the second Eigen argument if the first argument is true + * and otherwise return the second Eigen argument. Overload for use with one + * scalar and one Eigen object. If chosen, the scalar is returned as an Eigen + * object of the same size and type as the provided argument. + * + * @tparam T_true type of the true argument + * @tparam T_false type of the false argument + * @param c Boolean condition value. + * @param y_true Value to return if condition is true. + * @param y_false Value to return if condition is false. + */ template * = nullptr, require_eigen_t* = nullptr> -inline plain_type_t select(const bool c, const T_true y_true, const T_false y_false) { +inline plain_type_t select(const bool c, const T_true y_true, + const T_false y_false) { if (c) { - return y_false.unaryExpr([&](auto&& y){ return y_true; }); + return y_false.unaryExpr([&](auto&& y){ return y_true; }); } return y_false; } +/** + * Return the second argument if the first argument is true + * and otherwise return the third argument. Overload for use with an Eigen + * object of booleans, and two scalars. The chosen scalar is returned as an + * Eigen object of the same dimension as the input Eigen argument + * + * @tparam T_bool type of Eigen boolean object + * @tparam T_true type of the true argument + * @tparam T_false type of the false argument + * @param c Eigen object of boolean condition values. + * @param y_true Value to return if condition is true. + * @param y_false Value to return if condition is false. + */ template * = nullptr, require_all_stan_scalar_t* = nullptr> @@ -73,6 +114,18 @@ inline auto select(const T_bool c, const T_true y_true, }).eval(); } +/** + * Return the second argument if the first argument is true + * and otherwise return the third argument. Overload for use with an Eigen + * object of booleans, and at least one Eigen object as input. + * + * @tparam T_bool type of Eigen boolean object + * @tparam T_true type of the true argument + * @tparam T_false type of the false argument + * @param c Eigen object of boolean condition values. + * @param y_true Value to return if condition is true. + * @param y_false Value to return if condition is false. + */ template * = nullptr, require_any_eigen_array_t* = nullptr> From f03c35f0dbde313d3f615177fd96dcafa5b1695c Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Sun, 10 Jul 2022 13:59:28 -0400 Subject: [PATCH 08/30] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/select.hpp | 27 ++++++++++----------------- 1 file changed, 10 insertions(+), 17 deletions(-) diff --git a/stan/math/prim/fun/select.hpp b/stan/math/prim/fun/select.hpp index a84c2a4f160..4dff54a7204 100644 --- a/stan/math/prim/fun/select.hpp +++ b/stan/math/prim/fun/select.hpp @@ -38,9 +38,9 @@ inline auto select(const bool c, const T_true y_true, const T_false y_false) { template * = nullptr> inline auto select(const bool c, const T_true y_true, const T_false y_false) { - return y_true.binaryExpr(y_false, [&](auto&& x, auto&& y) { - return c ? x : y; - }).eval(); + return y_true + .binaryExpr(y_false, [&](auto&& x, auto&& y) { return c ? x : y; }) + .eval(); } /** @@ -55,16 +55,15 @@ inline auto select(const bool c, const T_true y_true, const T_false y_false) { * @param y_true Value to return if condition is true. * @param y_false Value to return if condition is false. */ -template * = nullptr, +template * = nullptr, require_stan_scalar_t* = nullptr> inline plain_type_t select(const bool c, const T_true& y_true, - const T_false& y_false) { + const T_false& y_false) { if (c) { return y_true; } - return y_true.unaryExpr([&](auto&& y){ return y_false; }); + return y_true.unaryExpr([&](auto&& y) { return y_false; }); } /** @@ -85,7 +84,7 @@ template select(const bool c, const T_true y_true, const T_false y_false) { if (c) { - return y_false.unaryExpr([&](auto&& y){ return y_true; }); + return y_false.unaryExpr([&](auto&& y) { return y_true; }); } return y_false; @@ -107,11 +106,8 @@ inline plain_type_t select(const bool c, const T_true y_true, template * = nullptr, require_all_stan_scalar_t* = nullptr> -inline auto select(const T_bool c, const T_true y_true, - const T_false y_false) { - return c.unaryExpr([&](bool cond){ - return cond ? y_true : y_false; - }).eval(); +inline auto select(const T_bool c, const T_true y_true, const T_false y_false) { + return c.unaryExpr([&](bool cond) { return cond ? y_true : y_false; }).eval(); } /** @@ -129,13 +125,10 @@ inline auto select(const T_bool c, const T_true y_true, template * = nullptr, require_any_eigen_array_t* = nullptr> -inline auto select(const T_bool c, - const T_true y_true, - const T_false y_false) { +inline auto select(const T_bool c, const T_true y_true, const T_false y_false) { return c.select(y_true, y_false).eval(); } - } // namespace math } // namespace stan From 2b3cf01633ca808eccd05081443806bfbba77bd3 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 10 Jul 2022 21:28:22 +0300 Subject: [PATCH 09/30] Fix headers --- stan/math/prim/prob/bernoulli_lccdf.hpp | 1 - stan/math/prim/prob/bernoulli_lcdf.hpp | 1 - 2 files changed, 2 deletions(-) diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index b54e813f511..1c701de9bed 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -33,7 +33,6 @@ template * = nullptr> return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { - using T_partials_return = partials_return_t; using T_theta_ref = ref_type_t; using std::log; static const char* function = "bernoulli_lccdf"; diff --git a/stan/math/prim/prob/bernoulli_lcdf.hpp b/stan/math/prim/prob/bernoulli_lcdf.hpp index 800a94b4a7d..a1e827e83bb 100644 --- a/stan/math/prim/prob/bernoulli_lcdf.hpp +++ b/stan/math/prim/prob/bernoulli_lcdf.hpp @@ -33,7 +33,6 @@ template * = nullptr> return_type_t bernoulli_lcdf(const T_n& n, const T_prob& theta) { - using T_partials_return = partials_return_t; using T_theta_ref = ref_type_t; using std::log; static const char* function = "bernoulli_lcdf"; From a371c2bf35902614c10ab3041fd8cc364c7fd966 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 11 Jul 2022 11:07:45 +0300 Subject: [PATCH 10/30] Fix return types --- stan/math/prim/fun/select.hpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/fun/select.hpp b/stan/math/prim/fun/select.hpp index 4dff54a7204..c9ada8a8cdc 100644 --- a/stan/math/prim/fun/select.hpp +++ b/stan/math/prim/fun/select.hpp @@ -55,9 +55,12 @@ inline auto select(const bool c, const T_true y_true, const T_false y_false) { * @param y_true Value to return if condition is true. * @param y_false Value to return if condition is false. */ -template * = nullptr, +template , + plain_type_t>, + require_eigen_t* = nullptr, require_stan_scalar_t* = nullptr> -inline plain_type_t select(const bool c, const T_true& y_true, +inline ReturnT select(const bool c, const T_true& y_true, const T_false& y_false) { if (c) { return y_true; @@ -79,9 +82,11 @@ inline plain_type_t select(const bool c, const T_true& y_true, * @param y_false Value to return if condition is false. */ template , + plain_type_t>, require_stan_scalar_t* = nullptr, require_eigen_t* = nullptr> -inline plain_type_t select(const bool c, const T_true y_true, +inline ReturnT select(const bool c, const T_true y_true, const T_false y_false) { if (c) { return y_false.unaryExpr([&](auto&& y) { return y_true; }); From 2bade772120edf6a0925b2b95f5c4877a23f7cb3 Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Mon, 11 Jul 2022 04:08:42 -0400 Subject: [PATCH 11/30] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/select.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/select.hpp b/stan/math/prim/fun/select.hpp index c9ada8a8cdc..b6f57fc7ae2 100644 --- a/stan/math/prim/fun/select.hpp +++ b/stan/math/prim/fun/select.hpp @@ -61,7 +61,7 @@ template * = nullptr, require_stan_scalar_t* = nullptr> inline ReturnT select(const bool c, const T_true& y_true, - const T_false& y_false) { + const T_false& y_false) { if (c) { return y_true; } @@ -87,7 +87,7 @@ template * = nullptr, require_eigen_t* = nullptr> inline ReturnT select(const bool c, const T_true y_true, - const T_false y_false) { + const T_false y_false) { if (c) { return y_false.unaryExpr([&](auto&& y) { return y_true; }); } From 409c971da5e626e482a9d89664fe8731dfca0e94 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 11 Jul 2022 16:01:01 +0300 Subject: [PATCH 12/30] Ignore Eigen deprecation warning in expression tests --- Jenkinsfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Jenkinsfile b/Jenkinsfile index f3d438552fd..47ddd51ae78 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -416,7 +416,7 @@ pipeline { unstash 'MathSetup' script { sh "echo O=0 > make/local" - sh "echo CXX=${CLANG_CXX} -Werror >> make/local" + sh "echo CXX=${CLANG_CXX} -Werror -Wno-deprecated-declarations >> make/local" sh "python ./test/code_generator_test.py" sh "python ./test/signature_parser_test.py" sh "python ./test/statement_types_test.py" From 6220f1468dea24b1c082259285f04dc3c366b864 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 12 Jul 2022 09:35:36 +0300 Subject: [PATCH 13/30] Update compiler flags --- Jenkinsfile | 2 +- make/compiler_flags | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Jenkinsfile b/Jenkinsfile index 47ddd51ae78..f3d438552fd 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -416,7 +416,7 @@ pipeline { unstash 'MathSetup' script { sh "echo O=0 > make/local" - sh "echo CXX=${CLANG_CXX} -Werror -Wno-deprecated-declarations >> make/local" + sh "echo CXX=${CLANG_CXX} -Werror >> make/local" sh "python ./test/code_generator_test.py" sh "python ./test/signature_parser_test.py" sh "python ./test/statement_types_test.py" diff --git a/make/compiler_flags b/make/compiler_flags index b59b75b7f82..0e3f9afadd7 100644 --- a/make/compiler_flags +++ b/make/compiler_flags @@ -196,7 +196,7 @@ endif CXXFLAGS_OS += -D_REENTRANT ## silence warnings occuring due to the TBB and Eigen libraries -CXXFLAGS_WARNINGS += -Wno-ignored-attributes +CXXFLAGS_WARNINGS += -Wno-ignored-attributes -Wno-deprecated-declarations ################################################################################ # Setup OpenCL From 2c11c42119abfb2fbabe6bc41dc7980725783b89 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 12 Jul 2022 14:33:16 +0300 Subject: [PATCH 14/30] Update includes --- stan/math/opencl/kernel_generator/select.hpp | 1 + stan/math/prim/prob/bernoulli_lccdf.hpp | 1 + stan/math/prim/prob/bernoulli_lcdf.hpp | 1 + 3 files changed, 3 insertions(+) diff --git a/stan/math/opencl/kernel_generator/select.hpp b/stan/math/opencl/kernel_generator/select.hpp index c2e6717ce91..24b25e98534 100644 --- a/stan/math/opencl/kernel_generator/select.hpp +++ b/stan/math/opencl/kernel_generator/select.hpp @@ -3,6 +3,7 @@ #ifdef STAN_OPENCL #include +#include #include #include #include diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index 1c701de9bed..ddf1f48e529 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/prob/bernoulli_lcdf.hpp b/stan/math/prim/prob/bernoulli_lcdf.hpp index a1e827e83bb..2a668f10d23 100644 --- a/stan/math/prim/prob/bernoulli_lcdf.hpp +++ b/stan/math/prim/prob/bernoulli_lcdf.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include From da1f7a825cf67555e39bf6ac5eaa07c236bc1fe6 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 14 Nov 2022 14:42:56 +0200 Subject: [PATCH 15/30] review comments --- stan/math/prim/fun/select.hpp | 27 ++++++++++++++++--------- stan/math/prim/prob/bernoulli_lccdf.hpp | 3 +-- 2 files changed, 18 insertions(+), 12 deletions(-) diff --git a/stan/math/prim/fun/select.hpp b/stan/math/prim/fun/select.hpp index b6f57fc7ae2..1990554cf19 100644 --- a/stan/math/prim/fun/select.hpp +++ b/stan/math/prim/fun/select.hpp @@ -36,11 +36,17 @@ inline auto select(const bool c, const T_true y_true, const T_false y_false) { * @param y_false Value to return if condition is false. */ template * = nullptr> -inline auto select(const bool c, const T_true y_true, const T_false y_false) { - return y_true - .binaryExpr(y_false, [&](auto&& x, auto&& y) { return c ? x : y; }) - .eval(); + typename T_return = return_type_t, + typename T_true_plain = promote_scalar_t>, + typename T_false_plain = promote_scalar_t>, + require_all_eigen_t* = nullptr, + require_all_same_t* = nullptr> +inline T_true_plain select(const bool c, const T_true y_true, const T_false y_false) { + if (c) { + return y_true; + } else { + return y_false; + } } /** @@ -64,9 +70,9 @@ inline ReturnT select(const bool c, const T_true& y_true, const T_false& y_false) { if (c) { return y_true; + } else { + return y_true.unaryExpr([&](auto&& y) { return y_false; }); } - - return y_true.unaryExpr([&](auto&& y) { return y_false; }); } /** @@ -90,9 +96,9 @@ inline ReturnT select(const bool c, const T_true y_true, const T_false y_false) { if (c) { return y_false.unaryExpr([&](auto&& y) { return y_true; }); + } else { + return y_false; } - - return y_false; } /** @@ -129,7 +135,8 @@ inline auto select(const T_bool c, const T_true y_true, const T_false y_false) { */ template * = nullptr, - require_any_eigen_array_t* = nullptr> + require_any_eigen_array_t* = nullptr, + require_any_stan_scalar_t* = nullptr> inline auto select(const T_bool c, const T_true y_true, const T_false y_false) { return c.select(y_true, y_false).eval(); } diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index ddf1f48e529..b22c5f64b45 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -54,8 +54,7 @@ return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { // The gradients are technically ill-defined, but treated as zero if (sum(n_arr < 0)) { return ops_partials.build(0.0); - } - if (sum(n_arr >= 1)) { + } else if (sum(n_arr >= 1)) { return ops_partials.build(NEGATIVE_INFTY); } From abfa3647dec40d48e60fa5ab2506a5743c3128ee Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Mon, 14 Nov 2022 07:44:19 -0500 Subject: [PATCH 16/30] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/fun/select.hpp | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/stan/math/prim/fun/select.hpp b/stan/math/prim/fun/select.hpp index 1990554cf19..9202a4e3ae0 100644 --- a/stan/math/prim/fun/select.hpp +++ b/stan/math/prim/fun/select.hpp @@ -35,13 +35,15 @@ inline auto select(const bool c, const T_true y_true, const T_false y_false) { * @param y_true Value to return if condition is true. * @param y_false Value to return if condition is false. */ -template , - typename T_true_plain = promote_scalar_t>, - typename T_false_plain = promote_scalar_t>, - require_all_eigen_t* = nullptr, - require_all_same_t* = nullptr> -inline T_true_plain select(const bool c, const T_true y_true, const T_false y_false) { +template < + typename T_true, typename T_false, + typename T_return = return_type_t, + typename T_true_plain = promote_scalar_t>, + typename T_false_plain = promote_scalar_t>, + require_all_eigen_t* = nullptr, + require_all_same_t* = nullptr> +inline T_true_plain select(const bool c, const T_true y_true, + const T_false y_false) { if (c) { return y_true; } else { From 9872fcffb9a7f53f341c7065fa5d4ce06cca46c9 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 11 Aug 2023 16:16:57 +0300 Subject: [PATCH 17/30] Tidy includes --- stan/math/prim/prob/bernoulli_cdf.hpp | 12 ++++-------- stan/math/prim/prob/bernoulli_lccdf.hpp | 22 ++++++++-------------- stan/math/prim/prob/bernoulli_lcdf.hpp | 17 ++++------------- 3 files changed, 16 insertions(+), 35 deletions(-) diff --git a/stan/math/prim/prob/bernoulli_cdf.hpp b/stan/math/prim/prob/bernoulli_cdf.hpp index a2e77b3d5a1..44184399611 100644 --- a/stan/math/prim/prob/bernoulli_cdf.hpp +++ b/stan/math/prim/prob/bernoulli_cdf.hpp @@ -3,13 +3,10 @@ #include #include -#include +#include #include -#include -#include #include #include -#include #include namespace stan { @@ -38,8 +35,8 @@ return_type_t bernoulli_cdf(const T_n& n, const T_prob& theta) { "Probability parameter", theta); T_theta_ref theta_ref = theta; const auto& n_arr = as_array_or_scalar(n); - check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0, - 1.0); + const auto& theta_arr = as_value_column_array_or_scalar(theta_ref); + check_bounded(function, "Probability parameter", theta_arr, 0.0, 1.0); if (size_zero(n, theta)) { return 1.0; @@ -49,10 +46,9 @@ return_type_t bernoulli_cdf(const T_n& n, const T_prob& theta) { // Explicit return for extreme values // The gradients are technically ill-defined, but treated as zero - if (sum(n_arr < 0)) { + if (any(n_arr < 0)) { return ops_partials.build(0.0); } - const auto& theta_arr = as_value_column_array_or_scalar(theta_ref); const auto& log1m_theta = select(theta_arr == 1, 0.0, log1m(theta_arr)); const auto& P1 = select(n_arr == 0, log1m_theta, 0.0); diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index c14959b971a..b8c6f32fee2 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -3,17 +3,12 @@ #include #include -#include +#include #include #include -#include -#include #include -#include #include -#include #include -#include namespace stan { namespace math { @@ -35,14 +30,13 @@ template * = nullptr> return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { using T_theta_ref = ref_type_t; - using std::log; static const char* function = "bernoulli_lccdf"; check_consistent_sizes(function, "Random variable", n, "Probability parameter", theta); T_theta_ref theta_ref = theta; const auto& n_arr = as_array_or_scalar(n); - check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0, - 1.0); + const auto& theta_arr = as_value_column_array_or_scalar(theta_ref); + check_bounded(function, "Probability parameter", theta_arr, 0.0, 1.0); if (size_zero(n, theta)) { return 0.0; @@ -52,19 +46,19 @@ return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { // Explicit return for extreme values // The gradients are technically ill-defined, but treated as zero - if (sum(n_arr < 0)) { + if (any(n_arr < 0)) { return ops_partials.build(0.0); - } else if (sum(n_arr >= 1)) { + } else if (any(n_arr >= 1)) { return ops_partials.build(NEGATIVE_INFTY); } - const auto& theta_arr = as_value_column_array_or_scalar(theta_ref); + const auto& theta_broadcast = select(true, theta_arr, n_arr); if (!is_constant_all::value) { - partials<0>(ops_partials) = select(true, inv(theta_arr), n_arr); + partials<0>(ops_partials) = inv(theta_broadcast); } - return ops_partials.build(sum(select(true, log(theta_arr), n_arr))); + return ops_partials.build(sum(log(theta_broadcast))); } } // namespace math diff --git a/stan/math/prim/prob/bernoulli_lcdf.hpp b/stan/math/prim/prob/bernoulli_lcdf.hpp index e29b7acc0a2..384c083cd73 100644 --- a/stan/math/prim/prob/bernoulli_lcdf.hpp +++ b/stan/math/prim/prob/bernoulli_lcdf.hpp @@ -3,17 +3,10 @@ #include #include -#include -#include -#include -#include -#include +#include #include -#include #include -#include #include -#include namespace stan { namespace math { @@ -35,14 +28,13 @@ template * = nullptr> return_type_t bernoulli_lcdf(const T_n& n, const T_prob& theta) { using T_theta_ref = ref_type_t; - using std::log; static const char* function = "bernoulli_lcdf"; check_consistent_sizes(function, "Random variable", n, "Probability parameter", theta); T_theta_ref theta_ref = theta; const auto& n_arr = as_array_or_scalar(n); - check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0, - 1.0); + const auto& theta_arr = as_value_column_array_or_scalar(theta_ref); + check_bounded(function, "Probability parameter", theta_arr, 0.0, 1.0); if (size_zero(n, theta)) { return 0.0; @@ -52,11 +44,10 @@ return_type_t bernoulli_lcdf(const T_n& n, const T_prob& theta) { // Explicit return for extreme values // The gradients are technically ill-defined, but treated as zero - if (sum(n_arr < 0)) { + if (any(n_arr < 0)) { return ops_partials.build(NEGATIVE_INFTY); } - const auto& theta_arr = as_value_column_array_or_scalar(theta_ref); const auto& log1m_theta = select(theta_arr == 1, 0.0, log1m(theta_arr)); if (!is_constant_all::value) { From bd11e545b452a2041916168a0e43297c0f1ed648 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 11 Aug 2023 16:24:41 +0300 Subject: [PATCH 18/30] Missing headers --- make/compiler_flags | 2 +- stan/math/prim/prob/bernoulli_lccdf.hpp | 1 + stan/math/prim/prob/bernoulli_lcdf.hpp | 1 + 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/make/compiler_flags b/make/compiler_flags index 5e0e08d1534..1552381f232 100644 --- a/make/compiler_flags +++ b/make/compiler_flags @@ -196,7 +196,7 @@ endif CXXFLAGS_OS += -D_REENTRANT ## silence warnings occuring due to the TBB and Eigen libraries -CXXFLAGS_WARNINGS += -Wno-ignored-attributes -Wno-deprecated-declarations +CXXFLAGS_WARNINGS += -Wno-ignored-attributes ################################################################################ # Setup OpenCL diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index b8c6f32fee2..63f4a677be3 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/prob/bernoulli_lcdf.hpp b/stan/math/prim/prob/bernoulli_lcdf.hpp index 384c083cd73..9606955f9ce 100644 --- a/stan/math/prim/prob/bernoulli_lcdf.hpp +++ b/stan/math/prim/prob/bernoulli_lcdf.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include From 8f5cdb86a64339e118b96d4835a124a12be8af64 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 11 Aug 2023 16:29:01 +0300 Subject: [PATCH 19/30] Reduce unnecessary computation --- stan/math/prim/prob/bernoulli_lccdf.hpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index 63f4a677be3..e007120b309 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -53,13 +53,12 @@ return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { return ops_partials.build(NEGATIVE_INFTY); } - const auto& theta_broadcast = select(true, theta_arr, n_arr); - + // Use select() to broadcast theta values & gradients if necessary if (!is_constant_all::value) { - partials<0>(ops_partials) = inv(theta_broadcast); + partials<0>(ops_partials) = select(true, inv(theta_arr), n_arr); } - return ops_partials.build(sum(log(theta_broadcast))); + return ops_partials.build(sum(select(true, log(theta_arr), n_arr))); } } // namespace math From cabafd658f5969f1ca4b6695a6ea55c13c17d99c Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 16 Aug 2023 21:12:06 +0300 Subject: [PATCH 20/30] Remove select broadcast hack --- stan/math/prim/prob/bernoulli_lccdf.hpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index e007120b309..579118d73fc 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -53,12 +53,15 @@ return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { return ops_partials.build(NEGATIVE_INFTY); } - // Use select() to broadcast theta values & gradients if necessary + size_t theta_size = math::size(theta_arr); + size_t n_size = math::size(n_arr); + double broadcast_n = theta_size == n_size ? 1 : std::fmax(theta_size, n_size); + if (!is_constant_all::value) { - partials<0>(ops_partials) = select(true, inv(theta_arr), n_arr); + partials<0>(ops_partials) = inv(theta_arr) * broadcast_n; } - return ops_partials.build(sum(select(true, log(theta_arr), n_arr))); + return ops_partials.build(sum(log(theta_arr)) * broadcast_n); } } // namespace math From 47c4f9a0c1a71992fafcf66879d275ae646753a2 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 19 Aug 2023 17:38:51 +0300 Subject: [PATCH 21/30] Fix broadcast logic --- stan/math/prim/prob/bernoulli_lccdf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/prob/bernoulli_lccdf.hpp b/stan/math/prim/prob/bernoulli_lccdf.hpp index 579118d73fc..d6bc081c9d1 100644 --- a/stan/math/prim/prob/bernoulli_lccdf.hpp +++ b/stan/math/prim/prob/bernoulli_lccdf.hpp @@ -55,7 +55,7 @@ return_type_t bernoulli_lccdf(const T_n& n, const T_prob& theta) { size_t theta_size = math::size(theta_arr); size_t n_size = math::size(n_arr); - double broadcast_n = theta_size == n_size ? 1 : std::fmax(theta_size, n_size); + double broadcast_n = theta_size == n_size ? 1 : n_size; if (!is_constant_all::value) { partials<0>(ops_partials) = inv(theta_arr) * broadcast_n; From 7825fdf23d5157df30b82928e475f9f23523e331 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 28 Aug 2023 07:23:42 +0300 Subject: [PATCH 22/30] Cleanup numerical stability, constants, compound functions in kernels --- stan/math/opencl/kernels/device_functions/Phi.hpp | 4 ++-- .../math/opencl/kernels/device_functions/lbeta.hpp | 7 ++++--- .../kernels/device_functions/lgamma_stirling.hpp | 2 +- .../math/opencl/kernels/device_functions/logit.hpp | 2 +- .../opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp | 4 ++-- .../opencl/kernels/ordered_logistic_glm_lpmf.hpp | 14 ++------------ stan/math/opencl/kernels/ordered_logistic_lpmf.hpp | 14 ++------------ stan/math/opencl/kernels/tridiagonalization.hpp | 6 +++--- 8 files changed, 17 insertions(+), 36 deletions(-) diff --git a/stan/math/opencl/kernels/device_functions/Phi.hpp b/stan/math/opencl/kernels/device_functions/Phi.hpp index ba7ecb6a51f..6946fc6fd28 100644 --- a/stan/math/opencl/kernels/device_functions/Phi.hpp +++ b/stan/math/opencl/kernels/device_functions/Phi.hpp @@ -24,11 +24,11 @@ static const char* phi_device_function if (x < -37.5) { return 0; } else if (x < -5.0) { - return 0.5 * erfc(-1.0 / sqrt(2.0) * x); + return 0.5 * erfc(-M_SQRT1_2 * x); } else if (x > 8.25) { return 1; } else { - return 0.5 * (1.0 + erf(1.0 / sqrt(2.0) * x)); + return 0.5 * (1.0 + erf(M_SQRT1_2 * x)); } } // \cond diff --git a/stan/math/opencl/kernels/device_functions/lbeta.hpp b/stan/math/opencl/kernels/device_functions/lbeta.hpp index f0aeb8061fb..e5c377672d4 100644 --- a/stan/math/opencl/kernels/device_functions/lbeta.hpp +++ b/stan/math/opencl/kernels/device_functions/lbeta.hpp @@ -95,12 +95,13 @@ static const char* lbeta_device_function return lgamma(x) + lgamma(y) - lgamma(x + y); } double x_over_xy = x / (x + y); + double log_xpy = log(x + y); if (x < LGAMMA_STIRLING_DIFF_USEFUL) { // y large, x small double stirling_diff = lgamma_stirling_diff(y) - lgamma_stirling_diff(x + y); double stirling - = (y - 0.5) * log1p(-x_over_xy) + x * (1 - log(x + y)); + = (y - 0.5) * log1p(-x_over_xy) + x * (1 - log_xpy); return stirling + lgamma(x) + stirling_diff; } @@ -108,8 +109,8 @@ static const char* lbeta_device_function double stirling_diff = lgamma_stirling_diff(x) + lgamma_stirling_diff(y) - lgamma_stirling_diff(x + y); - double stirling = (x - 0.5) * log(x_over_xy) + y * log1p(-x_over_xy) - + 0.5 * log(2.0 * M_PI) - 0.5 * log(y); + double stirling = (x - 0.5) * (log(x) - log_xpy) + y * log1p(-x_over_xy) + + 0.5 * (M_LN2 + log(M_PI)) - 0.5 * log(y); return stirling + stirling_diff; } // \cond diff --git a/stan/math/opencl/kernels/device_functions/lgamma_stirling.hpp b/stan/math/opencl/kernels/device_functions/lgamma_stirling.hpp index 19fc87b968f..3934fe11600 100644 --- a/stan/math/opencl/kernels/device_functions/lgamma_stirling.hpp +++ b/stan/math/opencl/kernels/device_functions/lgamma_stirling.hpp @@ -28,7 +28,7 @@ static const char* lgamma_stirling_device_function * @return Stirling's approximation to lgamma(x). */ double lgamma_stirling(double x) { - return 0.5 * log(2.0 * M_PI) + (x - 0.5) * log(x) - x; + return 0.5 * (M_LN2 + log(M_PI)) + (x - 0.5) * log(x) - x; } // \cond ) "\n#endif\n"; // NOLINT diff --git a/stan/math/opencl/kernels/device_functions/logit.hpp b/stan/math/opencl/kernels/device_functions/logit.hpp index 8ec1bf8850b..b0b4524214e 100644 --- a/stan/math/opencl/kernels/device_functions/logit.hpp +++ b/stan/math/opencl/kernels/device_functions/logit.hpp @@ -49,7 +49,7 @@ static const char* logit_device_function * @param x argument * @return log odds of argument */ - double logit(double x) { return log(x / (1 - x)); } + double logit(double x) { return log(x) - log1m(x); } // \cond ) "\n#endif\n"; // NOLINT // \endcond diff --git a/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp index f79b4449f69..b52284c0faa 100644 --- a/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp +++ b/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp @@ -92,9 +92,9 @@ static const char* neg_binomial_2_log_glm_kernel_code = STRINGIFY( double log_phi = log(phi); double logsumexp_theta_logphi; if (theta > log_phi) { - logsumexp_theta_logphi = theta + log1p(exp(log_phi - theta)); + logsumexp_theta_logphi = theta + log1p_exp(log_phi - theta); } else { - logsumexp_theta_logphi = log_phi + log1p(exp(theta - log_phi)); + logsumexp_theta_logphi = log_phi + log1p_exp(theta - log_phi); } double y_plus_phi = y + phi; if (need_logp1) { diff --git a/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp b/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp index 74341ff76cc..2f07e2dae93 100644 --- a/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp +++ b/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp @@ -87,20 +87,10 @@ static const char* ordered_logistic_glm_kernel_code = STRINGIFY( if (need_location_derivative || need_cuts_derivative) { double exp_cuts_diff = exp(cut_y2 - cut_y1); - if (cut2 > 0) { - double exp_m_cut2 = exp(-cut2); - d1 = exp_m_cut2 / (1 + exp_m_cut2); - } else { - d1 = 1 / (1 + exp(cut2)); - } + d1 = inv_logit(-cut2); d1 -= exp_cuts_diff / (exp_cuts_diff - 1); d2 = 1 / (1 - exp_cuts_diff); - if (cut1 > 0) { - double exp_m_cut1 = exp(-cut1); - d2 -= exp_m_cut1 / (1 + exp_m_cut1); - } else { - d2 -= 1 / (1 + exp(cut1)); - } + d2 -= inv_logit(-cut1); if (need_location_derivative) { location_derivative[gid] = d1 - d2; diff --git a/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp b/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp index b06a6d47a46..89dda0d287b 100644 --- a/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp +++ b/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp @@ -83,20 +83,10 @@ static const char* ordered_logistic_kernel_code = STRINGIFY( if (need_lambda_derivative || need_cuts_derivative) { double exp_cuts_diff = exp(cut_y2 - cut_y1); - if (cut2 > 0) { - double exp_m_cut2 = exp(-cut2); - d1 = exp_m_cut2 / (1 + exp_m_cut2); - } else { - d1 = 1 / (1 + exp(cut2)); - } + d1 = inv_logit(-cut2); d1 -= exp_cuts_diff / (exp_cuts_diff - 1); d2 = 1 / (1 - exp_cuts_diff); - if (cut1 > 0) { - double exp_m_cut1 = exp(-cut1); - d2 -= exp_m_cut1 / (1 + exp_m_cut1); - } else { - d2 -= 1 / (1 + exp(cut1)); - } + d2 -= inv_logit(-cut1); if (need_lambda_derivative) { lambda_derivative[gid] = d1 - d2; diff --git a/stan/math/opencl/kernels/tridiagonalization.hpp b/stan/math/opencl/kernels/tridiagonalization.hpp index b0bf3d3288a..0e0c11ca5ed 100644 --- a/stan/math/opencl/kernels/tridiagonalization.hpp +++ b/stan/math/opencl/kernels/tridiagonalization.hpp @@ -84,7 +84,7 @@ static const char* tridiagonalization_householder_kernel_code = STRINGIFY( q = q_local[0]; alpha = q_local[1]; if (q != 0) { - double multi = sqrt(2.) / q; + double multi = M_SQRT2 / q; // normalize the Householder vector for (int i = lid + 1; i < P_span; i += lsize) { P[P_start + i] *= multi; @@ -92,7 +92,7 @@ static const char* tridiagonalization_householder_kernel_code = STRINGIFY( } if (gid == 0) { P[P_rows * (k + j + 1) + k + j] - = P[P_rows * (k + j) + k + j + 1] * q / sqrt(2.) + alpha; + = P[P_rows * (k + j) + k + j + 1] * q / M_SQRT2 + alpha; } } // \cond @@ -291,7 +291,7 @@ static const char* tridiagonalization_v_step_3_kernel_code = STRINGIFY( v[i] -= acc * u[i]; } if (gid == 0) { - P[P_rows * (k + j + 1) + k + j] -= *q / sqrt(2.) * u[0]; + P[P_rows * (k + j + 1) + k + j] -= *q / M_SQRT2 * u[0]; } } // \cond From 71653a3141caf154238de6c4959ac9db64719dea Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 28 Aug 2023 00:27:40 -0400 Subject: [PATCH 23/30] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/opencl/kernels/device_functions/lbeta.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stan/math/opencl/kernels/device_functions/lbeta.hpp b/stan/math/opencl/kernels/device_functions/lbeta.hpp index e5c377672d4..0e8e4d139d3 100644 --- a/stan/math/opencl/kernels/device_functions/lbeta.hpp +++ b/stan/math/opencl/kernels/device_functions/lbeta.hpp @@ -109,7 +109,8 @@ static const char* lbeta_device_function double stirling_diff = lgamma_stirling_diff(x) + lgamma_stirling_diff(y) - lgamma_stirling_diff(x + y); - double stirling = (x - 0.5) * (log(x) - log_xpy) + y * log1p(-x_over_xy) + double stirling = (x - 0.5) * (log(x) - log_xpy) + + y * log1p(-x_over_xy) + 0.5 * (M_LN2 + log(M_PI)) - 0.5 * log(y); return stirling + stirling_diff; } From 7e4c5912551ed0586b337523985bd2b767844706 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 28 Aug 2023 10:11:12 +0300 Subject: [PATCH 24/30] Fix kernel inclusion --- stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp | 2 ++ stan/math/opencl/kernels/ordered_logistic_lpmf.hpp | 2 ++ 2 files changed, 4 insertions(+) diff --git a/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp b/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp index 2f07e2dae93..3b1727b1aa9 100644 --- a/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp +++ b/stan/math/opencl/kernels/ordered_logistic_glm_lpmf.hpp @@ -5,6 +5,7 @@ #include #include #include +#include namespace stan { namespace math { @@ -171,6 +172,7 @@ const kernel_cl ordered_logistic_glm("ordered_logistic_glm", {log1p_exp_device_function, log1m_exp_device_function, + inv_logit_device_function, ordered_logistic_glm_kernel_code}, {{"REDUCTION_STEP_SIZE", 4}, {"LOCAL_SIZE_", 64}}); diff --git a/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp b/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp index 89dda0d287b..a77f9e64688 100644 --- a/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp +++ b/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp @@ -5,6 +5,7 @@ #include #include #include +#include namespace stan { namespace math { @@ -165,6 +166,7 @@ const kernel_cl ordered_logistic("ordered_logistic", {log1p_exp_device_function, log1m_exp_device_function, + inv_logit_device_function, ordered_logistic_kernel_code}, {{"REDUCTION_STEP_SIZE", 4}, {"LOCAL_SIZE_", 64}}); From febea5f092078187eea18208f68c6d788d102b3e Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 28 Aug 2023 03:12:38 -0400 Subject: [PATCH 25/30] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/opencl/kernels/ordered_logistic_lpmf.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp b/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp index a77f9e64688..7b5dbddc169 100644 --- a/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp +++ b/stan/math/opencl/kernels/ordered_logistic_lpmf.hpp @@ -166,8 +166,7 @@ const kernel_cl ordered_logistic("ordered_logistic", {log1p_exp_device_function, log1m_exp_device_function, - inv_logit_device_function, - ordered_logistic_kernel_code}, + inv_logit_device_function, ordered_logistic_kernel_code}, {{"REDUCTION_STEP_SIZE", 4}, {"LOCAL_SIZE_", 64}}); } // namespace opencl_kernels From 0a2a902c9d4081b11e98ad34b9175e14dd2c60a1 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 28 Aug 2023 11:44:24 +0300 Subject: [PATCH 26/30] Missing kernel inclide --- stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp index b52284c0faa..a87f20b4d47 100644 --- a/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp +++ b/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp @@ -4,6 +4,7 @@ #include #include +#include namespace stan { namespace math { @@ -197,6 +198,7 @@ const kernel_cl neg_binomial_2_log_glm("neg_binomial_2_log_glm", {digamma_device_function, + log1p_exp_device_function, neg_binomial_2_log_glm_kernel_code}, {{"REDUCTION_STEP_SIZE", 4}, {"LOCAL_SIZE_", 64}}); From 9c219fe20fedcb202857040eb5f1c49e04d8ac57 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 28 Aug 2023 04:45:32 -0400 Subject: [PATCH 27/30] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp b/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp index a87f20b4d47..0585df47654 100644 --- a/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp +++ b/stan/math/opencl/kernels/neg_binomial_2_log_glm_lpmf.hpp @@ -197,8 +197,7 @@ const kernel_cl neg_binomial_2_log_glm("neg_binomial_2_log_glm", - {digamma_device_function, - log1p_exp_device_function, + {digamma_device_function, log1p_exp_device_function, neg_binomial_2_log_glm_kernel_code}, {{"REDUCTION_STEP_SIZE", 4}, {"LOCAL_SIZE_", 64}}); From 97c3d7f15007d6df1ba47c8c5218f7f1ddc2d8a9 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 28 Aug 2023 14:14:47 +0300 Subject: [PATCH 28/30] logit include --- stan/math/opencl/kernel_generator/elt_function_cl.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stan/math/opencl/kernel_generator/elt_function_cl.hpp b/stan/math/opencl/kernel_generator/elt_function_cl.hpp index 8c2c66ac226..87388162f18 100644 --- a/stan/math/opencl/kernel_generator/elt_function_cl.hpp +++ b/stan/math/opencl/kernel_generator/elt_function_cl.hpp @@ -307,7 +307,8 @@ ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_square, opencl_kernels::inv_square_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_logit, opencl_kernels::inv_logit_device_function) -ADD_UNARY_FUNCTION_WITH_INCLUDES(logit, opencl_kernels::logit_device_function) +ADD_UNARY_FUNCTION_WITH_INCLUDES(logit, opencl_kernels::logit_device_function, + opencl_kernels::log1m_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi, opencl_kernels::phi_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi_approx, opencl_kernels::inv_logit_device_function, From 0aae3775eff72438335b6d86b39b4d81e6f4cc97 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Mon, 28 Aug 2023 20:10:12 +0300 Subject: [PATCH 29/30] Fix test failures --- stan/math/opencl/kernel_generator/elt_function_cl.hpp | 4 ++-- stan/math/opencl/kernels/device_functions/inv_logit.hpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/stan/math/opencl/kernel_generator/elt_function_cl.hpp b/stan/math/opencl/kernel_generator/elt_function_cl.hpp index 87388162f18..a2851709e98 100644 --- a/stan/math/opencl/kernel_generator/elt_function_cl.hpp +++ b/stan/math/opencl/kernel_generator/elt_function_cl.hpp @@ -307,8 +307,8 @@ ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_square, opencl_kernels::inv_square_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_logit, opencl_kernels::inv_logit_device_function) -ADD_UNARY_FUNCTION_WITH_INCLUDES(logit, opencl_kernels::logit_device_function, - opencl_kernels::log1m_device_function) +ADD_UNARY_FUNCTION_WITH_INCLUDES(logit, opencl_kernels::log1m_device_function, + opencl_kernels::logit_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi, opencl_kernels::phi_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi_approx, opencl_kernels::inv_logit_device_function, diff --git a/stan/math/opencl/kernels/device_functions/inv_logit.hpp b/stan/math/opencl/kernels/device_functions/inv_logit.hpp index 9cb9d56cc27..34c526e4fae 100644 --- a/stan/math/opencl/kernels/device_functions/inv_logit.hpp +++ b/stan/math/opencl/kernels/device_functions/inv_logit.hpp @@ -56,7 +56,7 @@ static const char* inv_logit_device_function */ double inv_logit(double x) { if (x < 0) { - if (x < log(2.2204460492503131E-16)) { + if (x < log(DBL_EPSILON)) { return exp(x); } return exp(x) / (1 + exp(x)); From 8c010a0c345728706de1f24fc9203792d8613f5a Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 28 Aug 2023 13:11:15 -0400 Subject: [PATCH 30/30] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/opencl/kernel_generator/elt_function_cl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/opencl/kernel_generator/elt_function_cl.hpp b/stan/math/opencl/kernel_generator/elt_function_cl.hpp index a2851709e98..e7dedc1c315 100644 --- a/stan/math/opencl/kernel_generator/elt_function_cl.hpp +++ b/stan/math/opencl/kernel_generator/elt_function_cl.hpp @@ -308,7 +308,7 @@ ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_square, ADD_UNARY_FUNCTION_WITH_INCLUDES(inv_logit, opencl_kernels::inv_logit_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(logit, opencl_kernels::log1m_device_function, - opencl_kernels::logit_device_function) + opencl_kernels::logit_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi, opencl_kernels::phi_device_function) ADD_UNARY_FUNCTION_WITH_INCLUDES(Phi_approx, opencl_kernels::inv_logit_device_function,