From 2b8ae140373ce37a0c6086324ae569c5277f8651 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 8 Jul 2022 06:11:07 +0300 Subject: [PATCH 01/21] 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/21] 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/21] 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/21] [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/21] 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/21] 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/21] 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/21] [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/21] 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/21] 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/21] [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/21] 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/21] 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/21] 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/21] 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/21] [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/21] 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/21] 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/21] 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/21] 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/21] 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;