From 08909203cda72f41593cbd8888bdeffa471d1665 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Thu, 26 May 2022 13:21:55 -0400 Subject: [PATCH 1/6] adds internals needed for supporting vector returning lpdfs --- .../prim/functor/operands_and_partials.hpp | 3 +- stan/math/prim/prob/normal_log.hpp | 4 +- stan/math/prim/prob/normal_lpdf.hpp | 141 ++++++++++++++++-- .../rev/functor/operands_and_partials.hpp | 62 ++++++++ test/unit/math/prim/prob/normal_log_test.cpp | 20 ++- 5 files changed, 211 insertions(+), 19 deletions(-) diff --git a/stan/math/prim/functor/operands_and_partials.hpp b/stan/math/prim/functor/operands_and_partials.hpp index 3d2f70b6fb4..aea2481467c 100644 --- a/stan/math/prim/functor/operands_and_partials.hpp +++ b/stan/math/prim/functor/operands_and_partials.hpp @@ -139,7 +139,8 @@ class operands_and_partials { * @param value the return value of the function we are compressing * @return the value with its derivative */ - inline double build(double value) const noexcept { return value; } + template + inline auto build(T&& value) const noexcept { return std::forward(value); } // These will always be 0 size base template instantiations (above). internal::ops_partials_edge> edge1_; diff --git a/stan/math/prim/prob/normal_log.hpp b/stan/math/prim/prob/normal_log.hpp index 962d6657f08..9458314173b 100644 --- a/stan/math/prim/prob/normal_log.hpp +++ b/stan/math/prim/prob/normal_log.hpp @@ -32,7 +32,7 @@ template inline return_type_t normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) { - return normal_lpdf(y, mu, sigma); + return normal_lpdf(y, mu, sigma); } /** \ingroup prob_dists @@ -42,7 +42,7 @@ template inline return_type_t normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) { - return normal_lpdf(y, mu, sigma); + return normal_lpdf(y, mu, sigma); } } // namespace math diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp index e5a1be8265a..25517a9fcf5 100644 --- a/stan/math/prim/prob/normal_lpdf.hpp +++ b/stan/math/prim/prob/normal_lpdf.hpp @@ -19,6 +19,123 @@ namespace stan { namespace math { +enum class ProbReturnType {Scalar, Vector}; + +template +struct prob_broadcaster; + +template +struct prob_broadcaster> { + T ret_; + template * = nullptr> + prob_broadcaster(EigArr&& x) : ret_(sum(std::forward(x))) {} + + template * = nullptr> + prob_broadcaster(Scalar&& x) : ret_(x) {} + + template * = nullptr> + inline auto operator=(EigArr&& x) { + ret_ = sum(x); + return *this; + } + + template * = nullptr> + inline auto operator=(Scalar x) { + ret_ = x; + return *this; + } + + template * = nullptr> + inline auto operator+=(EigArr&& x) { + ret_ += sum(x); + return *this; + } + + template * = nullptr> + inline auto operator+=(Scalar&& x) { + ret_ += x; + return *this; + } + + template * = nullptr> + inline auto operator-=(EigArr&& x) { + ret_ -= sum(x); + return *this; + } + + template * = nullptr> + inline auto operator-=(Scalar&& x) { + ret_ -= x; + return *this; + } + inline auto ret() noexcept { + return ret_; + } + template + static auto zero(T1&& /* */) { + return T(0); + } + +}; + +template +struct prob_broadcaster> { + T ret_; + template * = nullptr> + prob_broadcaster(EigArr&& x) : ret_(std::forward(x)) {} + + template * = nullptr> + inline auto operator=(EigArr&& x) { + ret_ = sum(x); + return *this; + } + + template * = nullptr> + inline auto operator=(Scalar x) { + ret_ = Eigen::Array, -1, 1>::Constant(x, ret_.size()); + return *this; + } + + template * = nullptr> + inline auto operator+=(EigArr&& x) { + ret_ += x; + return *this; + } + + template * = nullptr> + inline auto operator+=(Scalar&& x) { + ret_ += x; + return *this; + } + + template * = nullptr> + inline auto operator-=(EigArr&& x) { + ret_ -= x; + return *this; + } + + template * = nullptr> + inline auto operator-=(Scalar&& x) { + ret_ -= x; + return *this; + } + + inline auto&& ret() noexcept { + return std::move(ret_); + } + + template + static auto zero(T1&& size) { + return Eigen::Array, -1, 1>::Constant(0, size).eval(); + } + +}; + + + +template +using prob_return_t = prob_broadcaster, Eigen::Array, -1, 1>>>; + /** \ingroup prob_dists * The log of the normal density for the specified scalar(s) given * the specified mean(s) and deviation(s). y, mu, or sigma can @@ -38,10 +155,10 @@ namespace math { * @return The log of the product of the densities. * @throw std::domain_error if the scale is not positive. */ -template * = nullptr> -inline return_type_t normal_lpdf(const T_y& y, +inline auto normal_lpdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { using T_partials_return = partials_return_t; @@ -62,12 +179,13 @@ inline return_type_t normal_lpdf(const T_y& y, check_not_nan(function, "Random variable", y_val); check_finite(function, "Location parameter", mu_val); check_positive(function, "Scale parameter", sigma_val); - + using ret_t = prob_return_t; + const size_t N = max_size(y, mu, sigma); if (size_zero(y, mu, sigma)) { - return 0.0; + return ret_t::zero(N); } if (!include_summand::value) { - return 0.0; + return ret_t::zero(N); } operands_and_partials ops_partials( @@ -79,13 +197,16 @@ inline return_type_t normal_lpdf(const T_y& y, const auto& y_scaled_sq = to_ref_if::value>(y_scaled * y_scaled); - size_t N = max_size(y, mu, sigma); - T_partials_return logp = -0.5 * sum(y_scaled_sq); + prob_return_t logp = -0.5 * y_scaled_sq; if (include_summand::value) { logp += NEG_LOG_SQRT_TWO_PI * N; } if (include_summand::value) { - logp -= sum(log(sigma_val)) * N / size(sigma); + if (RetType == ProbReturnType::Scalar) { + logp -= sum(log(sigma_val)) * N / size(sigma); + } else { + logp -= log(sigma_val); + } } if (!is_constant_all::value) { @@ -103,11 +224,11 @@ inline return_type_t normal_lpdf(const T_y& y, ops_partials.edge2_.partials_ = std::move(scaled_diff); } } - return ops_partials.build(logp); + return ops_partials.build(logp.ret()); } template -inline return_type_t normal_lpdf(const T_y& y, +inline auto normal_lpdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { return normal_lpdf(y, mu, sigma); diff --git a/stan/math/rev/functor/operands_and_partials.hpp b/stan/math/rev/functor/operands_and_partials.hpp index 908d5803174..e43d1d275c1 100644 --- a/stan/math/rev/functor/operands_and_partials.hpp +++ b/stan/math/rev/functor/operands_and_partials.hpp @@ -74,6 +74,41 @@ inline void update_adjoints(StdVec1& x, const Vec2& y, const vari& z) { } } +template * = nullptr> +inline void update_adjoints(var_value& x, const T2& y, const T3& z) { + x.adj() += z.adj() * y; +} + +template * = nullptr, + require_not_var_matrix_t* = nullptr, + require_arithmetic_t* = nullptr, + require_eigen_t* = nullptr> +inline void update_adjoints(Scalar1 x, Scalar2 y, const T3& z) noexcept { + x.adj() += sum(z.adj() * y); +} +template * = nullptr, + require_st_arithmetic* = nullptr, + require_eigen_t* = nullptr> +inline void update_adjoints(Matrix1& x, const Matrix2& y, const T3& z) { + x.adj().array() += z.adj() * y.array(); +} + +template * = nullptr, typename T3, require_eigen_t* = nullptr> +inline constexpr void update_adjoints(Arith&& /* x */, Alt&& /* y */, + const T3& /* z */) noexcept {} + +template * = nullptr, + require_st_arithmetic* = nullptr, + require_eigen_t* = nullptr> +inline void update_adjoints(StdVec1& x, const Vec2& y, const T3& z) { + for (size_t i = 0; i < x.size(); ++i) { + update_adjoints(x[i], y[i], z[i]); + } +} + } // namespace internal /** \ingroup type_trait @@ -169,6 +204,33 @@ class operands_and_partials { } }); } + template * = nullptr> + auto build(T&& value) { + arena_t> ret = value; + reverse_pass_callback([ret, operand1 = edge1_.operand(), partial1 = edge1_.partial(), + operand2 = edge2_.operand(), partial2 = edge2_.partial(), + operand3 = edge3_.operand(), partial3 = edge3_.partial(), + operand4 = edge4_.operand(), partial4 = edge4_.partial(), + operand5 = edge5_.operand(), + partial5 = edge5_.partial()]() mutable { + if (!is_constant::value) { + internal::update_adjoints(operand1, partial1, vi); + } + if (!is_constant::value) { + internal::update_adjoints(operand2, partial2, vi); + } + if (!is_constant::value) { + internal::update_adjoints(operand3, partial3, vi); + } + if (!is_constant::value) { + internal::update_adjoints(operand4, partial4, vi); + } + if (!is_constant::value) { + internal::update_adjoints(operand5, partial5, vi); + } + }); + return plain_type_t>(ret); + } }; namespace internal { diff --git a/test/unit/math/prim/prob/normal_log_test.cpp b/test/unit/math/prim/prob/normal_log_test.cpp index fc4906fdc7b..b31fb0a100d 100644 --- a/test/unit/math/prim/prob/normal_log_test.cpp +++ b/test/unit/math/prim/prob/normal_log_test.cpp @@ -13,12 +13,20 @@ TEST(ProbNormal, log_matches_lpdf) { EXPECT_FLOAT_EQ((stan::math::normal_lpdf(y, mu, sigma)), (stan::math::normal_log(y, mu, sigma))); EXPECT_FLOAT_EQ( - (stan::math::normal_lpdf(y, mu, sigma)), - (stan::math::normal_log(y, mu, sigma))); + (stan::math::normal_lpdf(y, mu, sigma)), + (stan::math::normal_log(y, mu, sigma))); EXPECT_FLOAT_EQ( - (stan::math::normal_lpdf(y, mu, sigma)), - (stan::math::normal_log(y, mu, sigma))); + (stan::math::normal_lpdf(y, mu, sigma)), + (stan::math::normal_log(y, mu, sigma))); EXPECT_FLOAT_EQ( - (stan::math::normal_lpdf(y, mu, sigma)), - (stan::math::normal_log(y, mu, sigma))); + (stan::math::normal_lpdf(y, mu, sigma)), + (stan::math::normal_log(y, mu, sigma))); +} + +TEST(ProbNormal, test_vlpdf) { + Eigen::Matrix Y = Eigen::Matrix::Random(5); + Eigen::Matrix Mu = Eigen::Matrix::Random(5); + Eigen::Matrix Sigma = stan::math::abs(Eigen::Matrix::Random(5)); + Eigen::Matrix A = stan::math::normal_lpdf(Y, Mu, Sigma); + } From d888a0f411e921667396acea2dd456b5f70b9ff0 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Thu, 26 May 2022 13:39:09 -0400 Subject: [PATCH 2/6] get vector working for prim and reverse mode --- stan/math/prim/prob/normal_log.hpp | 4 ++-- stan/math/prim/prob/normal_lpdf.hpp | 19 ++++++++++--------- .../rev/functor/operands_and_partials.hpp | 12 ++++++------ test/unit/math/rev/prob/normal_log_test.cpp | 10 ++++++++++ 4 files changed, 28 insertions(+), 17 deletions(-) diff --git a/stan/math/prim/prob/normal_log.hpp b/stan/math/prim/prob/normal_log.hpp index 9458314173b..a420e805da8 100644 --- a/stan/math/prim/prob/normal_log.hpp +++ b/stan/math/prim/prob/normal_log.hpp @@ -29,7 +29,7 @@ namespace math { * @tparam T_loc Type of location parameter. */ template -inline return_type_t normal_log(const T_y& y, +inline auto normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) { return normal_lpdf(y, mu, sigma); @@ -39,7 +39,7 @@ inline return_type_t normal_log(const T_y& y, * @deprecated use normal_lpdf */ template -inline return_type_t normal_log(const T_y& y, +inline auto normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) { return normal_lpdf(y, mu, sigma); diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp index 25517a9fcf5..238d192c332 100644 --- a/stan/math/prim/prob/normal_lpdf.hpp +++ b/stan/math/prim/prob/normal_lpdf.hpp @@ -69,11 +69,12 @@ struct prob_broadcaster> { return *this; } inline auto ret() noexcept { + static_assert(!is_var::value, "NOOO"); return ret_; } - template - static auto zero(T1&& /* */) { - return T(0); + template + static auto zero(int /* */) { + return return_type_t(0); } }; @@ -86,7 +87,7 @@ struct prob_broadcaster> { template * = nullptr> inline auto operator=(EigArr&& x) { - ret_ = sum(x); + ret_ = x; return *this; } @@ -124,9 +125,9 @@ struct prob_broadcaster> { return std::move(ret_); } - template - static auto zero(T1&& size) { - return Eigen::Array, -1, 1>::Constant(0, size).eval(); + template + static auto zero(int size) { + return Eigen::Array, -1, 1>::Constant(0, size).eval(); } }; @@ -182,10 +183,10 @@ inline auto normal_lpdf(const T_y& y, using ret_t = prob_return_t; const size_t N = max_size(y, mu, sigma); if (size_zero(y, mu, sigma)) { - return ret_t::zero(N); + return ret_t::template zero(N); } if (!include_summand::value) { - return ret_t::zero(N); + return ret_t::template zero(N); } operands_and_partials ops_partials( diff --git a/stan/math/rev/functor/operands_and_partials.hpp b/stan/math/rev/functor/operands_and_partials.hpp index e43d1d275c1..4b948888ab5 100644 --- a/stan/math/rev/functor/operands_and_partials.hpp +++ b/stan/math/rev/functor/operands_and_partials.hpp @@ -87,7 +87,7 @@ template * = nullptr, require_st_arithmetic* = nullptr, require_eigen_t* = nullptr> @@ -214,19 +214,19 @@ class operands_and_partials { operand5 = edge5_.operand(), partial5 = edge5_.partial()]() mutable { if (!is_constant::value) { - internal::update_adjoints(operand1, partial1, vi); + internal::update_adjoints(operand1, partial1, ret); } if (!is_constant::value) { - internal::update_adjoints(operand2, partial2, vi); + internal::update_adjoints(operand2, partial2, ret); } if (!is_constant::value) { - internal::update_adjoints(operand3, partial3, vi); + internal::update_adjoints(operand3, partial3, ret); } if (!is_constant::value) { - internal::update_adjoints(operand4, partial4, vi); + internal::update_adjoints(operand4, partial4, ret); } if (!is_constant::value) { - internal::update_adjoints(operand5, partial5, vi); + internal::update_adjoints(operand5, partial5, ret); } }); return plain_type_t>(ret); diff --git a/test/unit/math/rev/prob/normal_log_test.cpp b/test/unit/math/rev/prob/normal_log_test.cpp index da5d17a6f3b..47480943e7b 100644 --- a/test/unit/math/rev/prob/normal_log_test.cpp +++ b/test/unit/math/rev/prob/normal_log_test.cpp @@ -21,3 +21,13 @@ TEST(ProbDistributionsNormal, intVsDouble) { EXPECT_FLOAT_EQ(lp1adj, lp2adj); } } + + +TEST(ProbNormal, test_vlpdf) { + using stan::math::var; + Eigen::Matrix Y = Eigen::Matrix::Random(5); + Eigen::Matrix Mu = Eigen::Matrix::Random(5); + Eigen::Matrix Sigma = stan::math::abs(Eigen::Matrix::Random(5)); + Eigen::Matrix A = stan::math::normal_lpdf(Y, Mu, Sigma); + +} From 793b4f8fd3cf8d7a59f3cacf50518bd6a25ffaf9 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Fri, 3 Jun 2022 13:27:18 -0400 Subject: [PATCH 3/6] prototype reducer to allow distributions to return vectors instead of scalars --- .../fwd/functor/operands_and_partials.hpp | 37 +++ .../prim/functor/operands_and_partials.hpp | 1 + stan/math/prim/functor/prob_reducer.hpp | 247 ++++++++++++++++++ stan/math/prim/prob/normal_lpdf.hpp | 127 +-------- .../rev/functor/operands_and_partials.hpp | 35 ++- .../math/fwd/prob/normal_vectorized_test.cpp | 11 + test/unit/math/mix/prob/normal_test.cpp | 17 +- .../math/mix/prob/test_distribution_ad.hpp | 172 ++++++++++++ test/unit/math/test_ad.hpp | 1 + 9 files changed, 526 insertions(+), 122 deletions(-) create mode 100644 stan/math/prim/functor/prob_reducer.hpp create mode 100644 test/unit/math/fwd/prob/normal_vectorized_test.cpp create mode 100644 test/unit/math/mix/prob/test_distribution_ad.hpp diff --git a/stan/math/fwd/functor/operands_and_partials.hpp b/stan/math/fwd/functor/operands_and_partials.hpp index bf8e087d521..d7e78c57f9a 100644 --- a/stan/math/fwd/functor/operands_and_partials.hpp +++ b/stan/math/fwd/functor/operands_and_partials.hpp @@ -27,6 +27,7 @@ class ops_partials_edge> { const Op& operand_; Dx dx() { return this->partials_[0] * this->operand_.d_; } + Dx dx_v() { return this->partials_[0] * this->operand_.d_; } }; } // namespace internal @@ -106,6 +107,13 @@ class operands_and_partials> { = edge1_.dx() + edge2_.dx() + edge3_.dx() + edge4_.dx() + edge5_.dx(); return T_return_type(value, deriv); } + + template * = nullptr> + auto build(EigVec&& value) { + Eigen::Array, -1, 1> ret(value.template cast>()); + ret.d() = (edge1_.dx_v() + edge2_.dx_v() + edge3_.dx_v() + edge4_.dx_v() + edge5_.dx_v()); + return ret; + } }; namespace internal { @@ -135,6 +143,14 @@ class ops_partials_edge>> { } return derivative; } + Eigen::Array dx_v() { + Eigen::Array derivative(this->operands_.size()); + for (size_t i = 0; i < this->operands_.size(); ++i) { + derivative[i] = this->partials_[i] * this->operands_[i].d_; + } + return derivative; + } + }; template @@ -161,6 +177,10 @@ class ops_partials_edge, R, C>> { } return derivative; } + + Eigen::Array dx_v() { + return this->partials_.array() * this->operands_.d().array(); + } }; // Multivariate; vectors of eigen types @@ -191,6 +211,14 @@ class ops_partials_edge, R, C>>> { } return derivative; } + + Eigen::Array dx_v() { + Eigen::Array derivative(this->operands_.size()); + for (size_t i = 0; i < this->operands_.size(); ++i) { + derivative[i] = this->partials_vec_[i].dot(this->operands_[i].d()); + } + return derivative; + } }; template @@ -220,6 +248,15 @@ class ops_partials_edge>>> { } return derivative; } + Eigen::Array dx_v() { + Eigen::Array derivative = Eigen::Array::Zero(this->operands_.size()); + for (size_t i = 0; i < this->operands_.size(); ++i) { + for (int j = 0; j < this->operands_[i].size(); ++j) { + derivative[i] = this->partials_vec_[i][j] * this->operands_[i][j].d_; + } + } + return derivative; + } }; } // namespace internal diff --git a/stan/math/prim/functor/operands_and_partials.hpp b/stan/math/prim/functor/operands_and_partials.hpp index aea2481467c..cd256a7211c 100644 --- a/stan/math/prim/functor/operands_and_partials.hpp +++ b/stan/math/prim/functor/operands_and_partials.hpp @@ -62,6 +62,7 @@ struct ops_partials_edge> { * expression returning zero. */ static constexpr double dx() noexcept { return 0.0; } + static constexpr double dx_v() noexcept { return 0.0; } /** * Return the size of the operand for the edge. For doubles this is a compile * time expression returning zero. diff --git a/stan/math/prim/functor/prob_reducer.hpp b/stan/math/prim/functor/prob_reducer.hpp new file mode 100644 index 00000000000..eb27b43b4f2 --- /dev/null +++ b/stan/math/prim/functor/prob_reducer.hpp @@ -0,0 +1,247 @@ +#ifndef STAN_MATH_PRIM_FUNCTOR_PROB_REDUCER_HPP +#define STAN_MATH_PRIM_FUNCTOR_PROB_REDUCER_HPP + + +#include + +namespace stan { +namespace math { + +/** + * Used by distributions to decide whether return type shoudl be Scalar or Vector. + */ +enum class ProbReturnType {Scalar, Vector}; + +/** + * For scalars performs summations and is a no-op reducer for eigen vectors. + * Used in the probability distributions for scalar or vector return types. + */ +template +struct prob_reducer; + +/** + * For scalars performs summations when given eigen types. + * @tparam A stan scalar type + */ +template +struct prob_reducer> { + T ret_; // Underlying return type + + /** + * Construct from an Eigen type + * @tparam EigArr A type inheriting from `Eigen::EigenBase` + * @param x will be summed and passed into `ret_`. + */ + template * = nullptr> + prob_reducer(EigArr&& x) : ret_(sum(std::forward(x))) {} + + /** + * Construct from an Eigen type while ignoring size argument passed. + * @tparam EigArr A type inheriting from `Eigen::EigenBase` + * @tparam Tossed An integral type + * @param x will be summed and passed into `ret_`. + */ + template * = nullptr, + require_integral_t* = nullptr> + prob_reducer(EigArr&& x, Tossed&& /* */) : ret_(sum(std::forward(x))) {} + + /** + * Construct from a scalar type. + * @tparam Scalar a scalar + * @param x passed to `ret_`. + */ + template * = nullptr> + prob_reducer(Scalar&& x) : ret_(x) {} + + /** + * Construct from an Eigen type while ignoring size argument passed. + * @tparam Scalar a scalar + * @tparam Tossed an integral type + * @param x will be summed and inserted into `ret_`. + */ + template * = nullptr> + prob_reducer(Scalar&& x, Tossed&& /* */) : ret_(x) {} + + /** + * Perform summation and then assignment + * @tparam EigArr A type inheriting from `Eigen::EigenBase` + */ + template * = nullptr> + inline auto operator=(EigArr&& x) { + ret_ = sum(x); + return *this; + } + + /** + * Assignment + * @tparam Scalar A scalar type + */ + template * = nullptr> + inline auto operator=(Scalar x) { + ret_ = x; + return *this; + } + + /** + * Perform summation and then `+=` + * @tparam EigArr A type inheriting from `Eigen::EigenBase` + * @param x Eigen object to be summed. + */ + template * = nullptr> + inline auto operator+=(EigArr&& x) { + ret_ += sum(x); + return *this; + } + + template * = nullptr> + inline auto operator+=(Scalar&& x) { + ret_ += x; + return *this; + } + + /** + * Perform summation and then `-=` + * @tparam EigArr A type inheriting from `Eigen::EigenBase` + * @param x Eigen object to be summed. + */ + template * = nullptr> + inline auto operator-=(EigArr&& x) { + ret_ -= sum(x); + return *this; + } + + template * = nullptr> + inline auto operator-=(Scalar&& x) { + ret_ -= x; + return *this; + } + + /** + * Return the underlying scalar return type. + */ + inline auto ret() noexcept { + return ret_; + } + + /** + * Return a zero value, used when distribution has special cases that + * immedietly return zero. + * @tparam Types types to deduce the overall return type of the function. + */ + template + static auto zero(int /* */) { + return return_type_t(0); + } + +}; + +template +struct prob_reducer> { + T ret_; + + /** + * Construct from an Eigen type + * @tparam EigArr A type inheriting from `Eigen::EigenBase` + * @param x will be forwarded into `ret_`. + */ + template * = nullptr> + prob_reducer(EigArr&& x) : ret_(std::forward(x)) {} + + /** + * Construct from an Eigen type while ignoring size argument passed. + * @tparam EigArr A type inheriting from `Eigen::EigenBase` + * @tparam Tossed An integral type + * @param x will be forwarded to `ret_`. + */ + template * = nullptr, require_integral_t* = nullptr> + prob_reducer(EigArr&& x, Size /* x */) : ret_(std::forward(x)) {} + + /** + * Construct from a scalar type. + * @tparam Scalar a scalar + * @tparam Size An integral type + * @param x passed to `ret_` along with size to fill with a base value. + * @param n The size `ret_` should be + */ + template * = nullptr, + require_integral_t* = nullptr> + prob_reducer(Scalar constant, Size n) : ret_(T::Constant(n, constant)) {} + + /** + * assignment + * @tparam EigArr A type inheriting from `Eigen::EigenBase` + */ + template * = nullptr> + inline auto operator=(EigArr&& x) { + ret_ = std::forward(x); + return *this; + } + + /** + * assignm scalar by propogating value over `ret_` + * @tparam Scalar a stan scalar + * @param x The value to fill `ret_` with. + */ + template * = nullptr> + inline auto operator=(Scalar x) { + ret_ = Eigen::Array, -1, 1>::Constant(x, ret_.size()); + return *this; + } + + template * = nullptr> + inline auto operator+=(EigArr&& x) { + ret_ += std::forward(x); + return *this; + } + + template * = nullptr> + inline auto operator+=(Scalar&& x) { + ret_ += x; + return *this; + } + + template * = nullptr> + inline auto operator-=(EigArr&& x) { + ret_ -= std::forward(x); + return *this; + } + + template * = nullptr> + inline auto operator-=(Scalar&& x) { + ret_ -= x; + return *this; + } + + /** + * Return the underlying scalar return type. Passed the underlying by + * moving it which can cause `ret_` to be uninitialized after. + */ + inline auto&& ret() noexcept { + return std::move(ret_); + } + + /** + * Return a zero value, used when distribution has special cases that + * immedietly return zero. + * @tparam Types types to deduce the overall return type of the function. + * @param size The size of the array to return + */ + template + static auto zero(int size) { + return Eigen::Array, -1, 1>::Constant(0, size).eval(); + } + +}; + +/** + * Generate a reducer with correct return type. + * @tparam ReturnType Either Scalar or Vector. + * @tparam Types A parameter pack of types to deduce the underlying scalar type from + */ +template +using prob_return_t = prob_reducer, Eigen::Array, -1, 1>>>; + +} +} + +#endif diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp index 238d192c332..f79db64ada7 100644 --- a/stan/math/prim/prob/normal_lpdf.hpp +++ b/stan/math/prim/prob/normal_lpdf.hpp @@ -14,129 +14,12 @@ #include #include #include +#include #include namespace stan { namespace math { -enum class ProbReturnType {Scalar, Vector}; - -template -struct prob_broadcaster; - -template -struct prob_broadcaster> { - T ret_; - template * = nullptr> - prob_broadcaster(EigArr&& x) : ret_(sum(std::forward(x))) {} - - template * = nullptr> - prob_broadcaster(Scalar&& x) : ret_(x) {} - - template * = nullptr> - inline auto operator=(EigArr&& x) { - ret_ = sum(x); - return *this; - } - - template * = nullptr> - inline auto operator=(Scalar x) { - ret_ = x; - return *this; - } - - template * = nullptr> - inline auto operator+=(EigArr&& x) { - ret_ += sum(x); - return *this; - } - - template * = nullptr> - inline auto operator+=(Scalar&& x) { - ret_ += x; - return *this; - } - - template * = nullptr> - inline auto operator-=(EigArr&& x) { - ret_ -= sum(x); - return *this; - } - - template * = nullptr> - inline auto operator-=(Scalar&& x) { - ret_ -= x; - return *this; - } - inline auto ret() noexcept { - static_assert(!is_var::value, "NOOO"); - return ret_; - } - template - static auto zero(int /* */) { - return return_type_t(0); - } - -}; - -template -struct prob_broadcaster> { - T ret_; - template * = nullptr> - prob_broadcaster(EigArr&& x) : ret_(std::forward(x)) {} - - template * = nullptr> - inline auto operator=(EigArr&& x) { - ret_ = x; - return *this; - } - - template * = nullptr> - inline auto operator=(Scalar x) { - ret_ = Eigen::Array, -1, 1>::Constant(x, ret_.size()); - return *this; - } - - template * = nullptr> - inline auto operator+=(EigArr&& x) { - ret_ += x; - return *this; - } - - template * = nullptr> - inline auto operator+=(Scalar&& x) { - ret_ += x; - return *this; - } - - template * = nullptr> - inline auto operator-=(EigArr&& x) { - ret_ -= x; - return *this; - } - - template * = nullptr> - inline auto operator-=(Scalar&& x) { - ret_ -= x; - return *this; - } - - inline auto&& ret() noexcept { - return std::move(ret_); - } - - template - static auto zero(int size) { - return Eigen::Array, -1, 1>::Constant(0, size).eval(); - } - -}; - - - -template -using prob_return_t = prob_broadcaster, Eigen::Array, -1, 1>>>; - /** \ingroup prob_dists * The log of the normal density for the specified scalar(s) given * the specified mean(s) and deviation(s). y, mu, or sigma can @@ -198,9 +81,13 @@ inline auto normal_lpdf(const T_y& y, const auto& y_scaled_sq = to_ref_if::value>(y_scaled * y_scaled); - prob_return_t logp = -0.5 * y_scaled_sq; + prob_return_t logp(-0.5 * y_scaled_sq, N); if (include_summand::value) { - logp += NEG_LOG_SQRT_TWO_PI * N; + if (RetType == ProbReturnType::Scalar) { + logp += NEG_LOG_SQRT_TWO_PI * N; + } else { + logp += NEG_LOG_SQRT_TWO_PI; + } } if (include_summand::value) { if (RetType == ProbReturnType::Scalar) { diff --git a/stan/math/rev/functor/operands_and_partials.hpp b/stan/math/rev/functor/operands_and_partials.hpp index 4b948888ab5..dca93b8f261 100644 --- a/stan/math/rev/functor/operands_and_partials.hpp +++ b/stan/math/rev/functor/operands_and_partials.hpp @@ -48,12 +48,25 @@ inline void update_adjoints(var_value& x, const T2& y, const vari& z) { x.adj() += z.adj() * y; } +template * = nullptr> +inline void update_adjoints(var_value& x, const T2& y, const var& z) { + x.adj() += z.adj() * y; +} + template * = nullptr, require_not_var_matrix_t* = nullptr, require_arithmetic_t* = nullptr> inline void update_adjoints(Scalar1 x, Scalar2 y, const vari& z) noexcept { x.adj() += z.adj() * y; } +template * = nullptr, + require_not_var_matrix_t* = nullptr, + require_arithmetic_t* = nullptr> +inline void update_adjoints(Scalar1 x, Scalar2 y, const var& z) noexcept { + x.adj() += z.adj() * y; +} + template * = nullptr, require_st_arithmetic* = nullptr> @@ -61,10 +74,21 @@ inline void update_adjoints(Matrix1& x, const Matrix2& y, const vari& z) { x.adj().array() += z.adj() * y.array(); } +template * = nullptr, + require_st_arithmetic* = nullptr> +inline void update_adjoints(Matrix1& x, const Matrix2& y, const var& z) { + x.adj().array() += z.adj() * y.array(); +} + template * = nullptr> inline constexpr void update_adjoints(Arith&& /* x */, Alt&& /* y */, const vari& /* z */) noexcept {} +template * = nullptr> +inline constexpr void update_adjoints(Arith&& /* x */, Alt&& /* y */, + const var& /* z */) noexcept {} + template * = nullptr, require_st_arithmetic* = nullptr> @@ -74,6 +98,15 @@ inline void update_adjoints(StdVec1& x, const Vec2& y, const vari& z) { } } +template * = nullptr, + require_st_arithmetic* = nullptr> +inline void update_adjoints(StdVec1& x, const Vec2& y, const var& z) { + for (size_t i = 0; i < x.size(); ++i) { + update_adjoints(x[i], y[i], z); + } +} + template * = nullptr> inline void update_adjoints(var_value& x, const T2& y, const T3& z) { @@ -95,7 +128,7 @@ inline void update_adjoints(Matrix1& x, const Matrix2& y, const T3& z) { x.adj().array() += z.adj() * y.array(); } -template * = nullptr, typename T3, require_eigen_t* = nullptr> +template * = nullptr, require_eigen_t* = nullptr> inline constexpr void update_adjoints(Arith&& /* x */, Alt&& /* y */, const T3& /* z */) noexcept {} diff --git a/test/unit/math/fwd/prob/normal_vectorized_test.cpp b/test/unit/math/fwd/prob/normal_vectorized_test.cpp new file mode 100644 index 00000000000..c116c215fe2 --- /dev/null +++ b/test/unit/math/fwd/prob/normal_vectorized_test.cpp @@ -0,0 +1,11 @@ +#include +#include + +TEST(ProbNormal, fvar_test_vlpdf) { + using stan::math::fvar; + Eigen::Matrix, -1, 1> Y = Eigen::Matrix::Random(5); + Eigen::Matrix, -1, 1> Mu = Eigen::Matrix::Random(5); + Eigen::Matrix, -1, 1> Sigma = stan::math::abs(Eigen::Matrix::Random(5)); + Eigen::Matrix, -1, 1> A = stan::math::normal_lpdf(Y, Mu, Sigma); + +} diff --git a/test/unit/math/mix/prob/normal_test.cpp b/test/unit/math/mix/prob/normal_test.cpp index 81136b202d2..3a7cd9f986b 100644 --- a/test/unit/math/mix/prob/normal_test.cpp +++ b/test/unit/math/mix/prob/normal_test.cpp @@ -1,5 +1,5 @@ #include -#include +#include TEST(mathMixScalFun, normal_lpdf) { auto f = [](const double mu, const double sigma) { @@ -10,3 +10,18 @@ TEST(mathMixScalFun, normal_lpdf) { stan::test::expect_ad(f(0, 1), 0.0); stan::test::expect_ad(f(0, 1), 1.7); } + +TEST(mathMixScalFun, vnormal_lpdf) { + auto f = [](const auto& y, const auto& mu, const auto& sigma) { + return stan::math::normal_lpdf(y, mu, sigma); + }; + Eigen::VectorXd y = Eigen::VectorXd::Random(5); + //y << 0, 0; + Eigen::VectorXd mu = Eigen::VectorXd::Random(5); + //mu << 0, 0; + Eigen::VectorXd sigma = stan::math::abs(Eigen::VectorXd::Random(5)); + //sigma << 1, 1; + stan::test::expect_ad_distribution(f, y, mu, sigma); + //stan::test::expect_ad(f, y, mu, sigma); + +} diff --git a/test/unit/math/mix/prob/test_distribution_ad.hpp b/test/unit/math/mix/prob/test_distribution_ad.hpp new file mode 100644 index 00000000000..c9676cf0ebc --- /dev/null +++ b/test/unit/math/mix/prob/test_distribution_ad.hpp @@ -0,0 +1,172 @@ +#ifndef TEST_UNIT_MATH_TEST_DISTRIBUTION_AD_HPP +#define TEST_UNIT_MATH_TEST_DISTRIBUTION_AD_HPP + +#include +#include +namespace stan { + namespace test { + template * = nullptr> + Cast do_cast(T&& x) { + return x; + } + template * = nullptr> + Cast do_cast(T&& x) { + return Cast(x.data(), x.data() + x.size()); + } + template * = nullptr> + Cast do_cast(T&& x) { + return x(0); + } + + using unary_base_set = std::tuple, + std::tuple>, + std::tuple>>; + + using binary_base_set = std::tuple, + std::tuple>, + std::tuple>, + std::tuple, double>, + std::tuple, std::vector>, + std::tuple, Eigen::Matrix>, + std::tuple, double>, + std::tuple, std::vector>, + std::tuple, Eigen::Matrix>>; + + + using ternary_base_sets = std::tuple, + std::tuple>, + std::tuple>, + std::tuple, double>, + std::tuple, std::vector>, + std::tuple, Eigen::Matrix>, + std::tuple, double>, + std::tuple, std::vector>, + std::tuple, Eigen::Matrix>, + std::tuple, double, double>, + std::tuple, double, std::vector>, + std::tuple, double, Eigen::Matrix>, + std::tuple, std::vector, double>, + std::tuple, std::vector, std::vector>, + std::tuple, std::vector, Eigen::Matrix>, + std::tuple, Eigen::Matrix, double>, + std::tuple, Eigen::Matrix, std::vector>, + std::tuple, Eigen::Matrix, Eigen::Matrix>, + std::tuple, double, double>, + std::tuple, double, std::vector>, + std::tuple, double, Eigen::Matrix>, + std::tuple, std::vector, double>, + std::tuple, std::vector, std::vector>, + std::tuple, std::vector, Eigen::Matrix>, + std::tuple, Eigen::Matrix, double>, + std::tuple, Eigen::Matrix, std::vector>, + std::tuple, Eigen::Matrix, Eigen::Matrix>>; + + +using quad_base_set = std::tuple, + std::tuple>, + std::tuple>, + std::tuple, double>, + std::tuple, std::vector>, + std::tuple, Eigen::Matrix>, + std::tuple, double>, + std::tuple, std::vector>, + std::tuple, Eigen::Matrix>, + std::tuple, double, double>, + std::tuple, double, std::vector>, + std::tuple, double, Eigen::Matrix>, + std::tuple, std::vector, double>, + std::tuple, std::vector, std::vector>, + std::tuple, std::vector, Eigen::Matrix>, + std::tuple, Eigen::Matrix, double>, + std::tuple, Eigen::Matrix, std::vector>, + std::tuple, Eigen::Matrix, Eigen::Matrix>, + std::tuple, double, double>, + std::tuple, double, std::vector>, + std::tuple, double, Eigen::Matrix>, + std::tuple, std::vector, double>, + std::tuple, std::vector, std::vector>, + std::tuple, std::vector, Eigen::Matrix>, + std::tuple, Eigen::Matrix, double>, + std::tuple, Eigen::Matrix, std::vector>, + std::tuple, Eigen::Matrix, Eigen::Matrix>, + std::tuple, double, double, double>, + std::tuple, double, double, std::vector>, + std::tuple, double, double, Eigen::Matrix>, + std::tuple, double, std::vector, double>, + std::tuple, double, std::vector, std::vector>, + std::tuple, double, std::vector, Eigen::Matrix>, + std::tuple, double, Eigen::Matrix, double>, + std::tuple, double, Eigen::Matrix, std::vector>, + std::tuple, double, Eigen::Matrix, Eigen::Matrix>, + std::tuple, std::vector, double, double>, + std::tuple, std::vector, double, std::vector>, + std::tuple, std::vector, double, Eigen::Matrix>, + std::tuple, std::vector, std::vector, double>, + std::tuple, std::vector, std::vector, std::vector>, + std::tuple, std::vector, std::vector, Eigen::Matrix>, + std::tuple, std::vector, Eigen::Matrix, double>, + std::tuple, std::vector, Eigen::Matrix, std::vector>, + std::tuple, std::vector, Eigen::Matrix, Eigen::Matrix>, + std::tuple, Eigen::Matrix, double, double>, + std::tuple, Eigen::Matrix, double, std::vector>, + std::tuple, Eigen::Matrix, double, Eigen::Matrix>, + std::tuple, Eigen::Matrix, std::vector, double>, + std::tuple, Eigen::Matrix, std::vector, std::vector>, + std::tuple, Eigen::Matrix, std::vector, Eigen::Matrix>, + std::tuple, Eigen::Matrix, Eigen::Matrix, double>, + std::tuple, Eigen::Matrix, Eigen::Matrix, std::vector>, + std::tuple, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>, + std::tuple, double, double, double>, + std::tuple, double, double, std::vector>, + std::tuple, double, double, Eigen::Matrix>, + std::tuple, double, std::vector, double>, + std::tuple, double, std::vector, std::vector>, + std::tuple, double, std::vector, Eigen::Matrix>, + std::tuple, double, Eigen::Matrix, double>, + std::tuple, double, Eigen::Matrix, std::vector>, + std::tuple, double, Eigen::Matrix, Eigen::Matrix>, + std::tuple, std::vector, double, double>, + std::tuple, std::vector, double, std::vector>, + std::tuple, std::vector, double, Eigen::Matrix>, + std::tuple, std::vector, std::vector, double>, + std::tuple, std::vector, std::vector, std::vector>, + std::tuple, std::vector, std::vector, Eigen::Matrix>, + std::tuple, std::vector, Eigen::Matrix, double>, + std::tuple, std::vector, Eigen::Matrix, std::vector>, + std::tuple, std::vector, Eigen::Matrix, Eigen::Matrix>, + std::tuple, Eigen::Matrix, double, double>, + std::tuple, Eigen::Matrix, double, std::vector>, + std::tuple, Eigen::Matrix, double, Eigen::Matrix>, + std::tuple, Eigen::Matrix, std::vector, double>, + std::tuple, Eigen::Matrix, std::vector, std::vector>, + std::tuple, Eigen::Matrix, std::vector, Eigen::Matrix>, + std::tuple, Eigen::Matrix, Eigen::Matrix, double>, + std::tuple, Eigen::Matrix, Eigen::Matrix, std::vector>, + std::tuple, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>>; + + + +template +using base_set_t = std::conditional_t>>>; + + + + template + void expect_ad_distribution(F&& f, const EigVecs&... args) { + // Need to test all combinations of scalars, Eigen vectors, and std vectors + using base_set = base_set_t; + // TODO Write a version of for each that just takes in a template + // and not an object + stan::math::for_each([&f, &args...](auto&& type_tuple) { + stan::math::apply([&f, &args...](auto&&... types) { + stan::test::expect_ad(ad_tolerances{}, f, do_cast>(args)...); + }, type_tuple); + }, base_set{}); + } + + } +} +#endif diff --git a/test/unit/math/test_ad.hpp b/test/unit/math/test_ad.hpp index 5a8fff9ec1d..a21de8354ec 100644 --- a/test/unit/math/test_ad.hpp +++ b/test/unit/math/test_ad.hpp @@ -2150,6 +2150,7 @@ std::vector square_test_matrices(int low, int high) { return xs; } + } // namespace test } // namespace stan #endif From f7b6c3dda3342d976148400edd48052cb8a8fb37 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Fri, 3 Jun 2022 13:29:00 -0400 Subject: [PATCH 4/6] clang format --- .../fwd/functor/operands_and_partials.hpp | 9 +- stan/math/prim/functor/prob_reducer.hpp | 51 ++- stan/math/prim/prob/normal_lpdf.hpp | 11 +- test/unit/math/mix/prob/normal_test.cpp | 12 +- .../math/mix/prob/test_distribution_ad.hpp | 424 +++++++++++------- 5 files changed, 309 insertions(+), 198 deletions(-) diff --git a/stan/math/fwd/functor/operands_and_partials.hpp b/stan/math/fwd/functor/operands_and_partials.hpp index d7e78c57f9a..e8d53d4a4a6 100644 --- a/stan/math/fwd/functor/operands_and_partials.hpp +++ b/stan/math/fwd/functor/operands_and_partials.hpp @@ -111,7 +111,8 @@ class operands_and_partials> { template * = nullptr> auto build(EigVec&& value) { Eigen::Array, -1, 1> ret(value.template cast>()); - ret.d() = (edge1_.dx_v() + edge2_.dx_v() + edge3_.dx_v() + edge4_.dx_v() + edge5_.dx_v()); + ret.d() = (edge1_.dx_v() + edge2_.dx_v() + edge3_.dx_v() + edge4_.dx_v() + + edge5_.dx_v()); return ret; } }; @@ -150,7 +151,6 @@ class ops_partials_edge>> { } return derivative; } - }; template @@ -215,7 +215,7 @@ class ops_partials_edge, R, C>>> { Eigen::Array dx_v() { Eigen::Array derivative(this->operands_.size()); for (size_t i = 0; i < this->operands_.size(); ++i) { - derivative[i] = this->partials_vec_[i].dot(this->operands_[i].d()); + derivative[i] = this->partials_vec_[i].dot(this->operands_[i].d()); } return derivative; } @@ -249,7 +249,8 @@ class ops_partials_edge>>> { return derivative; } Eigen::Array dx_v() { - Eigen::Array derivative = Eigen::Array::Zero(this->operands_.size()); + Eigen::Array derivative + = Eigen::Array::Zero(this->operands_.size()); for (size_t i = 0; i < this->operands_.size(); ++i) { for (int j = 0; j < this->operands_[i].size(); ++j) { derivative[i] = this->partials_vec_[i][j] * this->operands_[i][j].d_; diff --git a/stan/math/prim/functor/prob_reducer.hpp b/stan/math/prim/functor/prob_reducer.hpp index eb27b43b4f2..ef9bd5f897a 100644 --- a/stan/math/prim/functor/prob_reducer.hpp +++ b/stan/math/prim/functor/prob_reducer.hpp @@ -1,16 +1,16 @@ #ifndef STAN_MATH_PRIM_FUNCTOR_PROB_REDUCER_HPP #define STAN_MATH_PRIM_FUNCTOR_PROB_REDUCER_HPP - #include namespace stan { namespace math { /** - * Used by distributions to decide whether return type shoudl be Scalar or Vector. + * Used by distributions to decide whether return type shoudl be Scalar or + * Vector. */ -enum class ProbReturnType {Scalar, Vector}; +enum class ProbReturnType { Scalar, Vector }; /** * For scalars performs summations and is a no-op reducer for eigen vectors. @@ -25,7 +25,7 @@ struct prob_reducer; */ template struct prob_reducer> { - T ret_; // Underlying return type + T ret_; // Underlying return type /** * Construct from an Eigen type @@ -41,9 +41,11 @@ struct prob_reducer> { * @tparam Tossed An integral type * @param x will be summed and passed into `ret_`. */ - template * = nullptr, - require_integral_t* = nullptr> - prob_reducer(EigArr&& x, Tossed&& /* */) : ret_(sum(std::forward(x))) {} + template * = nullptr, + require_integral_t* = nullptr> + prob_reducer(EigArr&& x, Tossed&& /* */) + : ret_(sum(std::forward(x))) {} /** * Construct from a scalar type. @@ -59,7 +61,8 @@ struct prob_reducer> { * @tparam Tossed an integral type * @param x will be summed and inserted into `ret_`. */ - template * = nullptr> + template * = nullptr> prob_reducer(Scalar&& x, Tossed&& /* */) : ret_(x) {} /** @@ -119,9 +122,7 @@ struct prob_reducer> { /** * Return the underlying scalar return type. */ - inline auto ret() noexcept { - return ret_; - } + inline auto ret() noexcept { return ret_; } /** * Return a zero value, used when distribution has special cases that @@ -132,7 +133,6 @@ struct prob_reducer> { static auto zero(int /* */) { return return_type_t(0); } - }; template @@ -153,7 +153,8 @@ struct prob_reducer> { * @tparam Tossed An integral type * @param x will be forwarded to `ret_`. */ - template * = nullptr, require_integral_t* = nullptr> + template * = nullptr, + require_integral_t* = nullptr> prob_reducer(EigArr&& x, Size /* x */) : ret_(std::forward(x)) {} /** @@ -163,8 +164,9 @@ struct prob_reducer> { * @param x passed to `ret_` along with size to fill with a base value. * @param n The size `ret_` should be */ - template * = nullptr, - require_integral_t* = nullptr> + template * = nullptr, + require_integral_t* = nullptr> prob_reducer(Scalar constant, Size n) : ret_(T::Constant(n, constant)) {} /** @@ -216,9 +218,7 @@ struct prob_reducer> { * Return the underlying scalar return type. Passed the underlying by * moving it which can cause `ret_` to be uninitialized after. */ - inline auto&& ret() noexcept { - return std::move(ret_); - } + inline auto&& ret() noexcept { return std::move(ret_); } /** * Return a zero value, used when distribution has special cases that @@ -228,20 +228,23 @@ struct prob_reducer> { */ template static auto zero(int size) { - return Eigen::Array, -1, 1>::Constant(0, size).eval(); + return Eigen::Array, -1, 1>::Constant(0, size) + .eval(); } - }; /** * Generate a reducer with correct return type. * @tparam ReturnType Either Scalar or Vector. - * @tparam Types A parameter pack of types to deduce the underlying scalar type from + * @tparam Types A parameter pack of types to deduce the underlying scalar type + * from */ template -using prob_return_t = prob_reducer, Eigen::Array, -1, 1>>>; +using prob_return_t = prob_reducer, + Eigen::Array, -1, 1>>>; -} -} +} // namespace math +} // namespace stan #endif diff --git a/stan/math/prim/prob/normal_lpdf.hpp b/stan/math/prim/prob/normal_lpdf.hpp index f79db64ada7..f5c60e4db25 100644 --- a/stan/math/prim/prob/normal_lpdf.hpp +++ b/stan/math/prim/prob/normal_lpdf.hpp @@ -39,12 +39,11 @@ namespace math { * @return The log of the product of the densities. * @throw std::domain_error if the scale is not positive. */ -template * = nullptr> -inline auto normal_lpdf(const T_y& y, - const T_loc& mu, - const T_scale& sigma) { +inline auto normal_lpdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { using T_partials_return = partials_return_t; using T_y_ref = ref_type_if_t::value, T_y>; using T_mu_ref = ref_type_if_t::value, T_loc>; @@ -116,9 +115,7 @@ inline auto normal_lpdf(const T_y& y, } template -inline auto normal_lpdf(const T_y& y, - const T_loc& mu, - const T_scale& sigma) { +inline auto normal_lpdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { return normal_lpdf(y, mu, sigma); } diff --git a/test/unit/math/mix/prob/normal_test.cpp b/test/unit/math/mix/prob/normal_test.cpp index 3a7cd9f986b..95fd189ff08 100644 --- a/test/unit/math/mix/prob/normal_test.cpp +++ b/test/unit/math/mix/prob/normal_test.cpp @@ -13,15 +13,15 @@ TEST(mathMixScalFun, normal_lpdf) { TEST(mathMixScalFun, vnormal_lpdf) { auto f = [](const auto& y, const auto& mu, const auto& sigma) { - return stan::math::normal_lpdf(y, mu, sigma); + return stan::math::normal_lpdf( + y, mu, sigma); }; Eigen::VectorXd y = Eigen::VectorXd::Random(5); - //y << 0, 0; + // y << 0, 0; Eigen::VectorXd mu = Eigen::VectorXd::Random(5); - //mu << 0, 0; + // mu << 0, 0; Eigen::VectorXd sigma = stan::math::abs(Eigen::VectorXd::Random(5)); - //sigma << 1, 1; + // sigma << 1, 1; stan::test::expect_ad_distribution(f, y, mu, sigma); - //stan::test::expect_ad(f, y, mu, sigma); - + // stan::test::expect_ad(f, y, mu, sigma); } diff --git a/test/unit/math/mix/prob/test_distribution_ad.hpp b/test/unit/math/mix/prob/test_distribution_ad.hpp index c9676cf0ebc..70fc16fc5d0 100644 --- a/test/unit/math/mix/prob/test_distribution_ad.hpp +++ b/test/unit/math/mix/prob/test_distribution_ad.hpp @@ -4,169 +4,279 @@ #include #include namespace stan { - namespace test { - template * = nullptr> - Cast do_cast(T&& x) { - return x; - } - template * = nullptr> - Cast do_cast(T&& x) { - return Cast(x.data(), x.data() + x.size()); - } - template * = nullptr> - Cast do_cast(T&& x) { - return x(0); - } - - using unary_base_set = std::tuple, - std::tuple>, - std::tuple>>; - - using binary_base_set = std::tuple, - std::tuple>, - std::tuple>, - std::tuple, double>, - std::tuple, std::vector>, - std::tuple, Eigen::Matrix>, - std::tuple, double>, - std::tuple, std::vector>, - std::tuple, Eigen::Matrix>>; - - - using ternary_base_sets = std::tuple, - std::tuple>, - std::tuple>, - std::tuple, double>, - std::tuple, std::vector>, - std::tuple, Eigen::Matrix>, - std::tuple, double>, - std::tuple, std::vector>, - std::tuple, Eigen::Matrix>, - std::tuple, double, double>, - std::tuple, double, std::vector>, - std::tuple, double, Eigen::Matrix>, - std::tuple, std::vector, double>, - std::tuple, std::vector, std::vector>, - std::tuple, std::vector, Eigen::Matrix>, - std::tuple, Eigen::Matrix, double>, - std::tuple, Eigen::Matrix, std::vector>, - std::tuple, Eigen::Matrix, Eigen::Matrix>, - std::tuple, double, double>, - std::tuple, double, std::vector>, - std::tuple, double, Eigen::Matrix>, - std::tuple, std::vector, double>, - std::tuple, std::vector, std::vector>, - std::tuple, std::vector, Eigen::Matrix>, - std::tuple, Eigen::Matrix, double>, - std::tuple, Eigen::Matrix, std::vector>, - std::tuple, Eigen::Matrix, Eigen::Matrix>>; +namespace test { +template * = nullptr> +Cast do_cast(T&& x) { + return x; +} +template * = nullptr> +Cast do_cast(T&& x) { + return Cast(x.data(), x.data() + x.size()); +} +template * = nullptr> +Cast do_cast(T&& x) { + return x(0); +} +using unary_base_set + = std::tuple, std::tuple>, + std::tuple>>; -using quad_base_set = std::tuple, - std::tuple>, - std::tuple>, - std::tuple, double>, - std::tuple, std::vector>, - std::tuple, Eigen::Matrix>, - std::tuple, double>, - std::tuple, std::vector>, - std::tuple, Eigen::Matrix>, - std::tuple, double, double>, - std::tuple, double, std::vector>, - std::tuple, double, Eigen::Matrix>, - std::tuple, std::vector, double>, - std::tuple, std::vector, std::vector>, - std::tuple, std::vector, Eigen::Matrix>, - std::tuple, Eigen::Matrix, double>, - std::tuple, Eigen::Matrix, std::vector>, - std::tuple, Eigen::Matrix, Eigen::Matrix>, - std::tuple, double, double>, - std::tuple, double, std::vector>, - std::tuple, double, Eigen::Matrix>, - std::tuple, std::vector, double>, - std::tuple, std::vector, std::vector>, - std::tuple, std::vector, Eigen::Matrix>, - std::tuple, Eigen::Matrix, double>, - std::tuple, Eigen::Matrix, std::vector>, - std::tuple, Eigen::Matrix, Eigen::Matrix>, - std::tuple, double, double, double>, - std::tuple, double, double, std::vector>, - std::tuple, double, double, Eigen::Matrix>, - std::tuple, double, std::vector, double>, - std::tuple, double, std::vector, std::vector>, - std::tuple, double, std::vector, Eigen::Matrix>, - std::tuple, double, Eigen::Matrix, double>, - std::tuple, double, Eigen::Matrix, std::vector>, - std::tuple, double, Eigen::Matrix, Eigen::Matrix>, - std::tuple, std::vector, double, double>, - std::tuple, std::vector, double, std::vector>, - std::tuple, std::vector, double, Eigen::Matrix>, - std::tuple, std::vector, std::vector, double>, - std::tuple, std::vector, std::vector, std::vector>, - std::tuple, std::vector, std::vector, Eigen::Matrix>, - std::tuple, std::vector, Eigen::Matrix, double>, - std::tuple, std::vector, Eigen::Matrix, std::vector>, - std::tuple, std::vector, Eigen::Matrix, Eigen::Matrix>, - std::tuple, Eigen::Matrix, double, double>, - std::tuple, Eigen::Matrix, double, std::vector>, - std::tuple, Eigen::Matrix, double, Eigen::Matrix>, - std::tuple, Eigen::Matrix, std::vector, double>, - std::tuple, Eigen::Matrix, std::vector, std::vector>, - std::tuple, Eigen::Matrix, std::vector, Eigen::Matrix>, - std::tuple, Eigen::Matrix, Eigen::Matrix, double>, - std::tuple, Eigen::Matrix, Eigen::Matrix, std::vector>, - std::tuple, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>, - std::tuple, double, double, double>, - std::tuple, double, double, std::vector>, - std::tuple, double, double, Eigen::Matrix>, - std::tuple, double, std::vector, double>, - std::tuple, double, std::vector, std::vector>, - std::tuple, double, std::vector, Eigen::Matrix>, - std::tuple, double, Eigen::Matrix, double>, - std::tuple, double, Eigen::Matrix, std::vector>, - std::tuple, double, Eigen::Matrix, Eigen::Matrix>, - std::tuple, std::vector, double, double>, - std::tuple, std::vector, double, std::vector>, - std::tuple, std::vector, double, Eigen::Matrix>, - std::tuple, std::vector, std::vector, double>, - std::tuple, std::vector, std::vector, std::vector>, - std::tuple, std::vector, std::vector, Eigen::Matrix>, - std::tuple, std::vector, Eigen::Matrix, double>, - std::tuple, std::vector, Eigen::Matrix, std::vector>, - std::tuple, std::vector, Eigen::Matrix, Eigen::Matrix>, - std::tuple, Eigen::Matrix, double, double>, - std::tuple, Eigen::Matrix, double, std::vector>, - std::tuple, Eigen::Matrix, double, Eigen::Matrix>, - std::tuple, Eigen::Matrix, std::vector, double>, - std::tuple, Eigen::Matrix, std::vector, std::vector>, - std::tuple, Eigen::Matrix, std::vector, Eigen::Matrix>, - std::tuple, Eigen::Matrix, Eigen::Matrix, double>, - std::tuple, Eigen::Matrix, Eigen::Matrix, std::vector>, - std::tuple, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>>; +using binary_base_set = std::tuple< + std::tuple, std::tuple>, + std::tuple>, + std::tuple, double>, + std::tuple, std::vector>, + std::tuple, Eigen::Matrix>, + std::tuple, double>, + std::tuple, std::vector>, + std::tuple, + Eigen::Matrix>>; +using ternary_base_sets = std::tuple< + std::tuple, + std::tuple>, + std::tuple>, + std::tuple, double>, + std::tuple, std::vector>, + std::tuple, + Eigen::Matrix>, + std::tuple, double>, + std::tuple, + std::vector>, + std::tuple, + Eigen::Matrix>, + std::tuple, double, double>, + std::tuple, double, std::vector>, + std::tuple, double, + Eigen::Matrix>, + std::tuple, std::vector, double>, + std::tuple, std::vector, std::vector>, + std::tuple, std::vector, + Eigen::Matrix>, + std::tuple, Eigen::Matrix, + double>, + std::tuple, Eigen::Matrix, + std::vector>, + std::tuple, Eigen::Matrix, + Eigen::Matrix>, + std::tuple, double, double>, + std::tuple, double, + std::vector>, + std::tuple, double, + Eigen::Matrix>, + std::tuple, std::vector, + double>, + std::tuple, std::vector, + std::vector>, + std::tuple, std::vector, + Eigen::Matrix>, + std::tuple, + Eigen::Matrix, double>, + std::tuple, + Eigen::Matrix, std::vector>, + std::tuple, + Eigen::Matrix, + Eigen::Matrix>>; +using quad_base_set = std::tuple< + std::tuple, + std::tuple>, + std::tuple>, + std::tuple, double>, + std::tuple, std::vector>, + std::tuple, + Eigen::Matrix>, + std::tuple, + double>, + std::tuple, + std::vector>, + std::tuple, + Eigen::Matrix>, + std::tuple, double, double>, + std::tuple, double, std::vector>, + std::tuple, double, + Eigen::Matrix>, + std::tuple, std::vector, double>, + std::tuple, std::vector, + std::vector>, + std::tuple, std::vector, + Eigen::Matrix>, + std::tuple, + Eigen::Matrix, double>, + std::tuple, + Eigen::Matrix, std::vector>, + std::tuple, + Eigen::Matrix, + Eigen::Matrix>, + std::tuple, double, + double>, + std::tuple, double, + std::vector>, + std::tuple, double, + Eigen::Matrix>, + std::tuple, + std::vector, double>, + std::tuple, + std::vector, std::vector>, + std::tuple, + std::vector, Eigen::Matrix>, + std::tuple, + Eigen::Matrix, double>, + std::tuple, + Eigen::Matrix, std::vector>, + std::tuple, + Eigen::Matrix, + Eigen::Matrix>, + std::tuple, double, double, double>, + std::tuple, double, double, std::vector>, + std::tuple, double, double, + Eigen::Matrix>, + std::tuple, double, std::vector, double>, + std::tuple, double, std::vector, + std::vector>, + std::tuple, double, std::vector, + Eigen::Matrix>, + std::tuple, double, + Eigen::Matrix, double>, + std::tuple, double, + Eigen::Matrix, std::vector>, + std::tuple, double, + Eigen::Matrix, + Eigen::Matrix>, + std::tuple, std::vector, double, double>, + std::tuple, std::vector, double, + std::vector>, + std::tuple, std::vector, double, + Eigen::Matrix>, + std::tuple, std::vector, std::vector, + double>, + std::tuple, std::vector, std::vector, + std::vector>, + std::tuple, std::vector, std::vector, + Eigen::Matrix>, + std::tuple, std::vector, + Eigen::Matrix, double>, + std::tuple, std::vector, + Eigen::Matrix, std::vector>, + std::tuple, std::vector, + Eigen::Matrix, + Eigen::Matrix>, + std::tuple, Eigen::Matrix, + double, double>, + std::tuple, Eigen::Matrix, + double, std::vector>, + std::tuple, Eigen::Matrix, + double, Eigen::Matrix>, + std::tuple, Eigen::Matrix, + std::vector, double>, + std::tuple, Eigen::Matrix, + std::vector, std::vector>, + std::tuple, Eigen::Matrix, + std::vector, Eigen::Matrix>, + std::tuple, Eigen::Matrix, + Eigen::Matrix, double>, + std::tuple, Eigen::Matrix, + Eigen::Matrix, std::vector>, + std::tuple, Eigen::Matrix, + Eigen::Matrix, + Eigen::Matrix>, + std::tuple, double, double, + double>, + std::tuple, double, double, + std::vector>, + std::tuple, double, double, + Eigen::Matrix>, + std::tuple, double, + std::vector, double>, + std::tuple, double, + std::vector, std::vector>, + std::tuple, double, + std::vector, Eigen::Matrix>, + std::tuple, double, + Eigen::Matrix, double>, + std::tuple, double, + Eigen::Matrix, std::vector>, + std::tuple, double, + Eigen::Matrix, + Eigen::Matrix>, + std::tuple, std::vector, + double, double>, + std::tuple, std::vector, + double, std::vector>, + std::tuple, std::vector, + double, Eigen::Matrix>, + std::tuple, std::vector, + std::vector, double>, + std::tuple, std::vector, + std::vector, std::vector>, + std::tuple, std::vector, + std::vector, Eigen::Matrix>, + std::tuple, std::vector, + Eigen::Matrix, double>, + std::tuple, std::vector, + Eigen::Matrix, std::vector>, + std::tuple, std::vector, + Eigen::Matrix, + Eigen::Matrix>, + std::tuple, + Eigen::Matrix, double, double>, + std::tuple, + Eigen::Matrix, double, + std::vector>, + std::tuple, + Eigen::Matrix, double, + Eigen::Matrix>, + std::tuple, + Eigen::Matrix, std::vector, + double>, + std::tuple, + Eigen::Matrix, std::vector, + std::vector>, + std::tuple, + Eigen::Matrix, std::vector, + Eigen::Matrix>, + std::tuple, + Eigen::Matrix, + Eigen::Matrix, double>, + std::tuple, + Eigen::Matrix, + Eigen::Matrix, std::vector>, + std::tuple, + Eigen::Matrix, + Eigen::Matrix, + Eigen::Matrix>>; template -using base_set_t = std::conditional_t>>>; - - +using base_set_t = std::conditional_t< + sizeof_args == 1, unary_base_set, + std::conditional_t< + sizeof_args == 2, binary_base_set, + std::conditional_t>>>; - template - void expect_ad_distribution(F&& f, const EigVecs&... args) { - // Need to test all combinations of scalars, Eigen vectors, and std vectors - using base_set = base_set_t; - // TODO Write a version of for each that just takes in a template - // and not an object - stan::math::for_each([&f, &args...](auto&& type_tuple) { - stan::math::apply([&f, &args...](auto&&... types) { - stan::test::expect_ad(ad_tolerances{}, f, do_cast>(args)...); - }, type_tuple); - }, base_set{}); - } - - } +template +void expect_ad_distribution(F&& f, const EigVecs&... args) { + // Need to test all combinations of scalars, Eigen vectors, and std vectors + using base_set = base_set_t; + // Would be nice to have version of for_each that just took in a template + stan::math::for_each( + [&f, &args...](auto&& type_tuple) { + stan::math::apply( + [&f, &args...](auto&&... types) { + stan::test::expect_ad( + ad_tolerances{}, f, + do_cast>(args)...); + }, + type_tuple); + }, + base_set{}); } + +} // namespace test +} // namespace stan #endif From 5a0c081c3613c2075f17af2e0065824eae432fbd Mon Sep 17 00:00:00 2001 From: Stan BuildBot Date: Fri, 3 Jun 2022 13:43:59 -0400 Subject: [PATCH 5/6] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- .../prim/functor/operands_and_partials.hpp | 4 +- stan/math/prim/prob/normal_log.hpp | 8 +-- .../rev/functor/operands_and_partials.hpp | 56 ++++++++++--------- .../math/fwd/prob/normal_vectorized_test.cpp | 14 +++-- test/unit/math/prim/prob/normal_log_test.cpp | 23 ++++---- test/unit/math/rev/prob/normal_log_test.cpp | 9 +-- test/unit/math/test_ad.hpp | 1 - 7 files changed, 60 insertions(+), 55 deletions(-) diff --git a/stan/math/prim/functor/operands_and_partials.hpp b/stan/math/prim/functor/operands_and_partials.hpp index cd256a7211c..68907d5758b 100644 --- a/stan/math/prim/functor/operands_and_partials.hpp +++ b/stan/math/prim/functor/operands_and_partials.hpp @@ -141,7 +141,9 @@ class operands_and_partials { * @return the value with its derivative */ template - inline auto build(T&& value) const noexcept { return std::forward(value); } + inline auto build(T&& value) const noexcept { + return std::forward(value); + } // These will always be 0 size base template instantiations (above). internal::ops_partials_edge> edge1_; diff --git a/stan/math/prim/prob/normal_log.hpp b/stan/math/prim/prob/normal_log.hpp index a420e805da8..bb16d728ae9 100644 --- a/stan/math/prim/prob/normal_log.hpp +++ b/stan/math/prim/prob/normal_log.hpp @@ -29,9 +29,7 @@ namespace math { * @tparam T_loc Type of location parameter. */ template -inline auto normal_log(const T_y& y, - const T_loc& mu, - const T_scale& sigma) { +inline auto normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) { return normal_lpdf(y, mu, sigma); } @@ -39,9 +37,7 @@ inline auto normal_log(const T_y& y, * @deprecated use normal_lpdf */ template -inline auto normal_log(const T_y& y, - const T_loc& mu, - const T_scale& sigma) { +inline auto normal_log(const T_y& y, const T_loc& mu, const T_scale& sigma) { return normal_lpdf(y, mu, sigma); } diff --git a/stan/math/rev/functor/operands_and_partials.hpp b/stan/math/rev/functor/operands_and_partials.hpp index dca93b8f261..fec51b1725c 100644 --- a/stan/math/rev/functor/operands_and_partials.hpp +++ b/stan/math/rev/functor/operands_and_partials.hpp @@ -107,13 +107,15 @@ inline void update_adjoints(StdVec1& x, const Vec2& y, const var& z) { } } -template * = nullptr> +template < + typename T1, typename T2, typename T3, + require_all_kernel_expressions_and_none_scalar_t* = nullptr> inline void update_adjoints(var_value& x, const T2& y, const T3& z) { x.adj() += z.adj() * y; } -template * = nullptr, +template * = nullptr, require_not_var_matrix_t* = nullptr, require_arithmetic_t* = nullptr, require_eigen_t* = nullptr> @@ -128,7 +130,9 @@ inline void update_adjoints(Matrix1& x, const Matrix2& y, const T3& z) { x.adj().array() += z.adj() * y.array(); } -template * = nullptr, require_eigen_t* = nullptr> +template * = nullptr, + require_eigen_t* = nullptr> inline constexpr void update_adjoints(Arith&& /* x */, Alt&& /* y */, const T3& /* z */) noexcept {} @@ -240,28 +244,28 @@ class operands_and_partials { template * = nullptr> auto build(T&& value) { arena_t> ret = value; - reverse_pass_callback([ret, operand1 = edge1_.operand(), partial1 = edge1_.partial(), - operand2 = edge2_.operand(), partial2 = edge2_.partial(), - operand3 = edge3_.operand(), partial3 = edge3_.partial(), - operand4 = edge4_.operand(), partial4 = edge4_.partial(), - operand5 = edge5_.operand(), - partial5 = edge5_.partial()]() mutable { - if (!is_constant::value) { - internal::update_adjoints(operand1, partial1, ret); - } - if (!is_constant::value) { - internal::update_adjoints(operand2, partial2, ret); - } - if (!is_constant::value) { - internal::update_adjoints(operand3, partial3, ret); - } - if (!is_constant::value) { - internal::update_adjoints(operand4, partial4, ret); - } - if (!is_constant::value) { - internal::update_adjoints(operand5, partial5, ret); - } - }); + reverse_pass_callback( + [ret, operand1 = edge1_.operand(), partial1 = edge1_.partial(), + operand2 = edge2_.operand(), partial2 = edge2_.partial(), + operand3 = edge3_.operand(), partial3 = edge3_.partial(), + operand4 = edge4_.operand(), partial4 = edge4_.partial(), + operand5 = edge5_.operand(), partial5 = edge5_.partial()]() mutable { + if (!is_constant::value) { + internal::update_adjoints(operand1, partial1, ret); + } + if (!is_constant::value) { + internal::update_adjoints(operand2, partial2, ret); + } + if (!is_constant::value) { + internal::update_adjoints(operand3, partial3, ret); + } + if (!is_constant::value) { + internal::update_adjoints(operand4, partial4, ret); + } + if (!is_constant::value) { + internal::update_adjoints(operand5, partial5, ret); + } + }); return plain_type_t>(ret); } }; diff --git a/test/unit/math/fwd/prob/normal_vectorized_test.cpp b/test/unit/math/fwd/prob/normal_vectorized_test.cpp index c116c215fe2..b10f2fc833b 100644 --- a/test/unit/math/fwd/prob/normal_vectorized_test.cpp +++ b/test/unit/math/fwd/prob/normal_vectorized_test.cpp @@ -3,9 +3,13 @@ TEST(ProbNormal, fvar_test_vlpdf) { using stan::math::fvar; - Eigen::Matrix, -1, 1> Y = Eigen::Matrix::Random(5); - Eigen::Matrix, -1, 1> Mu = Eigen::Matrix::Random(5); - Eigen::Matrix, -1, 1> Sigma = stan::math::abs(Eigen::Matrix::Random(5)); - Eigen::Matrix, -1, 1> A = stan::math::normal_lpdf(Y, Mu, Sigma); - + Eigen::Matrix, -1, 1> Y + = Eigen::Matrix::Random(5); + Eigen::Matrix, -1, 1> Mu + = Eigen::Matrix::Random(5); + Eigen::Matrix, -1, 1> Sigma + = stan::math::abs(Eigen::Matrix::Random(5)); + Eigen::Matrix, -1, 1> A + = stan::math::normal_lpdf( + Y, Mu, Sigma); } diff --git a/test/unit/math/prim/prob/normal_log_test.cpp b/test/unit/math/prim/prob/normal_log_test.cpp index b31fb0a100d..e637d1055b1 100644 --- a/test/unit/math/prim/prob/normal_log_test.cpp +++ b/test/unit/math/prim/prob/normal_log_test.cpp @@ -12,21 +12,20 @@ TEST(ProbNormal, log_matches_lpdf) { (stan::math::normal_log(y, mu, sigma))); EXPECT_FLOAT_EQ((stan::math::normal_lpdf(y, mu, sigma)), (stan::math::normal_log(y, mu, sigma))); - EXPECT_FLOAT_EQ( - (stan::math::normal_lpdf(y, mu, sigma)), - (stan::math::normal_log(y, mu, sigma))); - EXPECT_FLOAT_EQ( - (stan::math::normal_lpdf(y, mu, sigma)), - (stan::math::normal_log(y, mu, sigma))); - EXPECT_FLOAT_EQ( - (stan::math::normal_lpdf(y, mu, sigma)), - (stan::math::normal_log(y, mu, sigma))); + EXPECT_FLOAT_EQ((stan::math::normal_lpdf(y, mu, sigma)), + (stan::math::normal_log(y, mu, sigma))); + EXPECT_FLOAT_EQ((stan::math::normal_lpdf(y, mu, sigma)), + (stan::math::normal_log(y, mu, sigma))); + EXPECT_FLOAT_EQ((stan::math::normal_lpdf(y, mu, sigma)), + (stan::math::normal_log(y, mu, sigma))); } TEST(ProbNormal, test_vlpdf) { Eigen::Matrix Y = Eigen::Matrix::Random(5); Eigen::Matrix Mu = Eigen::Matrix::Random(5); - Eigen::Matrix Sigma = stan::math::abs(Eigen::Matrix::Random(5)); - Eigen::Matrix A = stan::math::normal_lpdf(Y, Mu, Sigma); - + Eigen::Matrix Sigma + = stan::math::abs(Eigen::Matrix::Random(5)); + Eigen::Matrix A + = stan::math::normal_lpdf( + Y, Mu, Sigma); } diff --git a/test/unit/math/rev/prob/normal_log_test.cpp b/test/unit/math/rev/prob/normal_log_test.cpp index 47480943e7b..5b5e10ca647 100644 --- a/test/unit/math/rev/prob/normal_log_test.cpp +++ b/test/unit/math/rev/prob/normal_log_test.cpp @@ -22,12 +22,13 @@ TEST(ProbDistributionsNormal, intVsDouble) { } } - TEST(ProbNormal, test_vlpdf) { using stan::math::var; Eigen::Matrix Y = Eigen::Matrix::Random(5); Eigen::Matrix Mu = Eigen::Matrix::Random(5); - Eigen::Matrix Sigma = stan::math::abs(Eigen::Matrix::Random(5)); - Eigen::Matrix A = stan::math::normal_lpdf(Y, Mu, Sigma); - + Eigen::Matrix Sigma + = stan::math::abs(Eigen::Matrix::Random(5)); + Eigen::Matrix A + = stan::math::normal_lpdf( + Y, Mu, Sigma); } diff --git a/test/unit/math/test_ad.hpp b/test/unit/math/test_ad.hpp index a21de8354ec..5a8fff9ec1d 100644 --- a/test/unit/math/test_ad.hpp +++ b/test/unit/math/test_ad.hpp @@ -2150,7 +2150,6 @@ std::vector square_test_matrices(int low, int high) { return xs; } - } // namespace test } // namespace stan #endif From bb1918ccab2a6257c5aeccb1c983f8041b0207e0 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Fri, 3 Jun 2022 16:31:42 -0400 Subject: [PATCH 6/6] remove unneeded constructors from prob_reducer --- stan/math/prim/functor/prob_reducer.hpp | 26 +------------------------ test/unit/math/mix/prob/normal_test.cpp | 4 ---- 2 files changed, 1 insertion(+), 29 deletions(-) diff --git a/stan/math/prim/functor/prob_reducer.hpp b/stan/math/prim/functor/prob_reducer.hpp index ef9bd5f897a..4509b7b85c8 100644 --- a/stan/math/prim/functor/prob_reducer.hpp +++ b/stan/math/prim/functor/prob_reducer.hpp @@ -27,14 +27,6 @@ template struct prob_reducer> { T ret_; // Underlying return type - /** - * Construct from an Eigen type - * @tparam EigArr A type inheriting from `Eigen::EigenBase` - * @param x will be summed and passed into `ret_`. - */ - template * = nullptr> - prob_reducer(EigArr&& x) : ret_(sum(std::forward(x))) {} - /** * Construct from an Eigen type while ignoring size argument passed. * @tparam EigArr A type inheriting from `Eigen::EigenBase` @@ -47,14 +39,6 @@ struct prob_reducer> { prob_reducer(EigArr&& x, Tossed&& /* */) : ret_(sum(std::forward(x))) {} - /** - * Construct from a scalar type. - * @tparam Scalar a scalar - * @param x passed to `ret_`. - */ - template * = nullptr> - prob_reducer(Scalar&& x) : ret_(x) {} - /** * Construct from an Eigen type while ignoring size argument passed. * @tparam Scalar a scalar @@ -139,18 +123,10 @@ template struct prob_reducer> { T ret_; - /** - * Construct from an Eigen type - * @tparam EigArr A type inheriting from `Eigen::EigenBase` - * @param x will be forwarded into `ret_`. - */ - template * = nullptr> - prob_reducer(EigArr&& x) : ret_(std::forward(x)) {} - /** * Construct from an Eigen type while ignoring size argument passed. * @tparam EigArr A type inheriting from `Eigen::EigenBase` - * @tparam Tossed An integral type + * @tparam Size An integral type * @param x will be forwarded to `ret_`. */ template * = nullptr, diff --git a/test/unit/math/mix/prob/normal_test.cpp b/test/unit/math/mix/prob/normal_test.cpp index 95fd189ff08..156f498e3ab 100644 --- a/test/unit/math/mix/prob/normal_test.cpp +++ b/test/unit/math/mix/prob/normal_test.cpp @@ -17,11 +17,7 @@ TEST(mathMixScalFun, vnormal_lpdf) { y, mu, sigma); }; Eigen::VectorXd y = Eigen::VectorXd::Random(5); - // y << 0, 0; Eigen::VectorXd mu = Eigen::VectorXd::Random(5); - // mu << 0, 0; Eigen::VectorXd sigma = stan::math::abs(Eigen::VectorXd::Random(5)); - // sigma << 1, 1; stan::test::expect_ad_distribution(f, y, mu, sigma); - // stan::test::expect_ad(f, y, mu, sigma); }