From 44f12df6cf2bc58403ccb422bb53fecb09023e60 Mon Sep 17 00:00:00 2001 From: michascant <89426143+MichaScant@users.noreply.github.com> Date: Mon, 25 Nov 2024 12:14:28 -0500 Subject: [PATCH 1/7] first commit --- stan/math/fwd/fun.hpp | 1 + stan/math/fwd/fun/log_add_exp.hpp | 123 ++++++++++++++++++++ stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/log_add_exp.hpp | 105 +++++++++++++++++ stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/log_add_exp.hpp | 88 ++++++++++++++ test/unit/math/mix/fun/log_add_exp_test.cpp | 71 +++++++++++ 7 files changed, 390 insertions(+) create mode 100644 stan/math/fwd/fun/log_add_exp.hpp create mode 100644 stan/math/prim/fun/log_add_exp.hpp create mode 100644 stan/math/rev/fun/log_add_exp.hpp create mode 100644 test/unit/math/mix/fun/log_add_exp_test.cpp diff --git a/stan/math/fwd/fun.hpp b/stan/math/fwd/fun.hpp index b56232a0121..a100117751c 100644 --- a/stan/math/fwd/fun.hpp +++ b/stan/math/fwd/fun.hpp @@ -80,6 +80,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/fwd/fun/log_add_exp.hpp b/stan/math/fwd/fun/log_add_exp.hpp new file mode 100644 index 00000000000..99cd0284bf7 --- /dev/null +++ b/stan/math/fwd/fun/log_add_exp.hpp @@ -0,0 +1,123 @@ +#ifndef STAN_MATH_FWD_FUN_LOG_ADD_EXP_HPP +#define STAN_MATH_FWD_FUN_LOG_ADD_EXP_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +// Overload for fvar and fvar +template +inline fvar log_add_exp(const fvar& x1, const fvar& x2) { + auto val = stan::math::log_add_exp(x1.val_, x2.val_); + + auto exp_x1 = stan::math::exp(x1.val_); + auto exp_x2 = stan::math::exp(x2.val_); + auto sum_exp = exp_x1 + exp_x2; + + auto grad1 = exp_x1 / sum_exp; + auto grad2 = exp_x2 / sum_exp; + + return fvar(val, x1.d_ * grad1 + x2.d_ * grad2); +} + +template +inline fvar log_add_exp(const fvar& x1, double x2) { + if (x1.val_ == NEGATIVE_INFTY) { + return fvar(x2, 0.0); // log_add_exp(-∞, b) = b + } + return log_add_exp(x2, x1); +} + +template +inline fvar log_add_exp(double x1, const fvar& x2) { + if (x2.val_ == NEGATIVE_INFTY) { + return fvar(x1, 0.0); // log_add_exp(a, -∞) = a + } + auto val = stan::math::log_add_exp(x1, x2.val_); + auto exp_x2 = stan::math::exp(x2.val_); + auto grad = exp_x2 / (stan::math::exp(x1) + exp_x2); + return fvar(val, x2.d_ * grad); +} + +// Overload for matrices of fvar +template +inline fvar log_add_exp(const Eigen::Matrix, -1, -1>& a, + const Eigen::Matrix, -1, -1>& b) { + using fvar_mat_type = Eigen::Matrix, -1, -1>; + fvar_mat_type result(a.rows(), a.cols()); + + // Check for empty inputs + if (a.size() == 0 || b.size() == 0) { + throw std::invalid_argument("Input containers must not be empty."); + } + + // Check for NaN + if (a.array().isNaN().any() || b.array().isNaN().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } + + // Check for infinity + if (a.array().isInf().any() || b.array().isInf().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } + + // Apply the log_add_exp operation directly + for (int i = 0; i < a.rows(); ++i) { + for (int j = 0; j < a.cols(); ++j) { + result(i, j) = stan::math::log_add_exp(a(i, j), b(i, j)); + } + } + + return result; // Return the result matrix +} + +// Specialization for nested fvar types +template +inline auto log_add_exp(const Eigen::Matrix>, -1, -1>& a, + const Eigen::Matrix>, -1, -1>& b) { + using nested_fvar_mat_type = Eigen::Matrix>, -1, -1>; + nested_fvar_mat_type result(a.rows(), a.cols()); + + // Check for empty inputs + if (a.size() == 0 || b.size() == 0) { + throw std::invalid_argument("Input containers must not be empty."); + } + + // Check for NaN + if (a.array().isNaN().any() || b.array().isNaN().any()) { + result.setConstant(stan::math::fvar>(std::numeric_limits::quiet_NaN())); + return result; + } + + // Check for infinity + if (a.array().isInf().any() || b.array().isInf().any()) { + result.setConstant(stan::math::fvar>(std::numeric_limits::quiet_NaN())); + return result; + } + + // Implement the logic for log_add_exp for nested fvar types + for (int i = 0; i < a.rows(); ++i) { + for (int j = 0; j < a.cols(); ++j) { + auto inner_a = a(i, j); + auto inner_b = b(i, j); + result(i, j) = stan::math::log_add_exp(inner_a, inner_b); + } + } + + return result; // Return the result matrix +} + +} +} + +#endif diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 21403200b41..1885331725e 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -188,6 +188,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/log_add_exp.hpp b/stan/math/prim/fun/log_add_exp.hpp new file mode 100644 index 00000000000..574c5dde1b1 --- /dev/null +++ b/stan/math/prim/fun/log_add_exp.hpp @@ -0,0 +1,105 @@ +#ifndef STAN_MATH_PRIM_FUN_LOG_ADD_EXP_HPP +#define STAN_MATH_PRIM_FUN_LOG_ADD_EXP_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Calculates the elementwise sum of exponentials without overflow. + * + * \f$\log (\exp(a) + \exp(b)) = m + \log(\exp(a-m) + \exp(b-m))\f$, + * + * where \f$m = max(a, b)\f$. + * + * @tparam T1 type of the first variable + * @tparam T2 type of the second variable + * @param a the first variable + * @param b the second variable + */ + +template * = nullptr, + require_all_stan_scalar_t* = nullptr> +inline return_type_t log_add_exp(const T2& a, const T1& b) { + if (a == NEGATIVE_INFTY) { + return b; + } + if (b == NEGATIVE_INFTY) { + return a; + } + if (a == INFTY || b == INFTY) { + return INFTY; + } + + const double max_val = std::max(a, b); + return max_val + std::log(std::exp(a - max_val) + std::exp(b - max_val)); +} + +/** + * Calculates the element-wise log sum of exponentials for two containers. + * For vectors a and b, computes log(exp(a[i]) + exp(b[i])) for each element i. + * If sizes don't match, uses the smaller size. + * + * @tparam T1 type of first container + * @tparam T2 type of second container + * @param a First input container + * @param b Second input container + * @return Container with element-wise log_add_exp results + */ +template * = nullptr> +inline auto log_add_exp(const T& a, const T& b) { + if (a.size() != b.size()) { + throw std::invalid_argument("Binary function: size of x (" + std::to_string(a.size()) + ") and size of y (" + std::to_string(b.size()) + ") must match in size"); + } + + const size_t min_size = std::min(a.size(), b.size()); + using return_t = return_type_t; + + std::vector result(min_size); + + for (size_t i = 0; i < min_size; ++i) { + if (a[i] == NEGATIVE_INFTY) { + result[i] = b[i]; // log_add_exp(-∞, b) = b + } else if (b[i] == NEGATIVE_INFTY) { + result[i] = a[i]; // log_add_exp(a, -∞) = a + } else if (a[i] == INFTY || b[i] == INFTY) { + result[i] = INFTY; // log_add_exp(∞, b) = ∞ + } else { + // Log-add-exp trick + const double max_val = std::max(a[i], b[i]); + result[i] = max_val + std::log(std::exp(a[i] - max_val) + std::exp(b[i] - max_val)); + } + } + + return result; +} + +/** + * Enables the vectorized application of the log_add_exp function, + * when the first and/or second arguments are containers. + * + * @tparam T1 + * @tparam T2 + * @param a + * @param b + * @return auto + */ +template * = nullptr> +inline auto log_add_exp(const T1& a, const T2& b) { + return apply_scalar_binary( + a, b, [](const auto& c, const auto& d) { return log_sum_exp(c, d); }); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index 929e88aa4c3..edc1d829a16 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -116,6 +116,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/log_add_exp.hpp b/stan/math/rev/fun/log_add_exp.hpp new file mode 100644 index 00000000000..eba7006e388 --- /dev/null +++ b/stan/math/rev/fun/log_add_exp.hpp @@ -0,0 +1,88 @@ +#ifndef STAN_MATH_REV_FUN_LOG_ADD_EXP_HPP +#define STAN_MATH_REV_FUN_LOG_ADD_EXP_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace stan { +namespace math { +namespace internal { + +class log_add_exp_vv_vari : public op_vv_vari { + public: + log_add_exp_vv_vari(vari* avi, vari* bvi) + : op_vv_vari(log_add_exp(avi->val_, bvi->val_), avi, bvi) {} + void chain() { + double exp_a = std::exp(avi_->val_); + double exp_b = std::exp(bvi_->val_); + double sum_exp = exp_a + exp_b; + + avi_->adj_ += adj_ * (exp_a / sum_exp); + bvi_->adj_ += adj_ * (exp_b / sum_exp); + } +}; + +class log_add_exp_vd_vari : public op_vd_vari { + public: + log_add_exp_vd_vari(vari* avi, double b) + : op_vd_vari(log_add_exp(avi->val_, b), avi, b) {} + void chain() { + if (val_ == NEGATIVE_INFTY) { + avi_->adj_ += adj_; + } else { + double exp_a = std::exp(avi_->val_); + avi_->adj_ += adj_ * (exp_a / (exp_a + std::exp(bd_))); + } + } +}; + +} // namespace internal + +/** + * Returns the element-wise log sum of exponentials. + */ +inline var log_add_exp(const var& a, const var& b) { + return var(new internal::log_add_exp_vv_vari(a.vi_, b.vi_)); +} + +/** + * Returns the log sum of exponentials. + */ +inline var log_add_exp(const var& a, double b) { + return var(new internal::log_add_exp_vd_vari(a.vi_, b)); +} + +/** + * Returns the element-wise log sum of exponentials. + */ +inline var log_add_exp(double a, const var& b) { + return var(new internal::log_add_exp_vd_vari(b.vi_, a)); +} + +/** + * Returns element-wise log sum of exponentials for Eigen types. + * + * @tparam T A type inheriting from EigenBase with var scalar type + * @param x First input + * @param y Second input + */ +template * = nullptr> +inline T log_add_exp(const T& x, const T& y) { + return apply_scalar_binary( + x, y, [](const auto& a, const auto& b) { return log_add_exp(a, b); }); +} + +} // namespace math +} // namespace stan + +#endif \ No newline at end of file diff --git a/test/unit/math/mix/fun/log_add_exp_test.cpp b/test/unit/math/mix/fun/log_add_exp_test.cpp new file mode 100644 index 00000000000..de2c1e834dc --- /dev/null +++ b/test/unit/math/mix/fun/log_add_exp_test.cpp @@ -0,0 +1,71 @@ +#include +#include +#include +#include + +TEST(MathMixMatFun, logAddExp) { + auto f = [](const auto& x, const auto& y) { + return stan::math::log_add_exp(x, y); + }; + // Test with finite values + Eigen::VectorXd x1(2); + x1 << 2.0, 1.0; + Eigen::VectorXd y1(2); + y1 << 3.0, 2.0; + stan::test::expect_ad(f, x1, y1); + + // Test with negative infinity + + stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); + stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); + + // Test with infinity + stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); +} + +TEST(MathMixMatFun, log_add_exp_elementwise_values) { + auto f = [](const auto& x, const auto& y) { + return stan::math::log_add_exp(x, y); + }; + + Eigen::VectorXd x1(2); + x1 << 2.0, 1.0; + Eigen::VectorXd y1(2); + y1 << 3.0, 2.0; + stan::test::expect_ad(f, x1, y1); + + Eigen::VectorXd x2(2); + x2 << 0.5, -1.0; + Eigen::VectorXd y2(2); + y2 << 1.0, 2.0; + stan::test::expect_ad(f, x2, y2); + + // Test with infinity + Eigen::VectorXd x3(2); + x3 << std::numeric_limits::infinity(), 1.0; + Eigen::VectorXd y3(2); + y3 << 2.0, std::numeric_limits::infinity(); + stan::test::expect_ad(f, x3, y3); +} + +TEST(MathMixMatFun, log_add_exp_edge_cases) { + auto f = [](const auto& x, const auto& y) { + return stan::math::log_add_exp(x, y); + }; + + stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); + stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); + stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); +} + +TEST(MathMixMatFun, log_add_exp_mismatched_sizes) { + auto f = [](const auto& x, const auto& y) { + return stan::math::log_add_exp(x, y); + }; + + std::vector x{1.0, 2.0}; + std::vector y{1.0, 2.0, 3.0}; + + stan::test::expect_ad(f, x, y); + stan::test::expect_ad(f, y, x); +} From 85c4a94bd48d4689bcfa21336343c6fe8d2dd0dd Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 25 Nov 2024 12:38:57 -0500 Subject: [PATCH 2/7] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/fwd/fun/log_add_exp.hpp | 165 ++++++++++---------- stan/math/prim/fun/log_add_exp.hpp | 60 +++---- stan/math/rev/fun/log_add_exp.hpp | 2 +- test/unit/math/mix/fun/log_add_exp_test.cpp | 4 +- 4 files changed, 120 insertions(+), 111 deletions(-) diff --git a/stan/math/fwd/fun/log_add_exp.hpp b/stan/math/fwd/fun/log_add_exp.hpp index 99cd0284bf7..f6c0827b1a0 100644 --- a/stan/math/fwd/fun/log_add_exp.hpp +++ b/stan/math/fwd/fun/log_add_exp.hpp @@ -16,108 +16,113 @@ namespace math { // Overload for fvar and fvar template inline fvar log_add_exp(const fvar& x1, const fvar& x2) { - auto val = stan::math::log_add_exp(x1.val_, x2.val_); - - auto exp_x1 = stan::math::exp(x1.val_); - auto exp_x2 = stan::math::exp(x2.val_); - auto sum_exp = exp_x1 + exp_x2; - - auto grad1 = exp_x1 / sum_exp; - auto grad2 = exp_x2 / sum_exp; - - return fvar(val, x1.d_ * grad1 + x2.d_ * grad2); + auto val = stan::math::log_add_exp(x1.val_, x2.val_); + + auto exp_x1 = stan::math::exp(x1.val_); + auto exp_x2 = stan::math::exp(x2.val_); + auto sum_exp = exp_x1 + exp_x2; + + auto grad1 = exp_x1 / sum_exp; + auto grad2 = exp_x2 / sum_exp; + + return fvar(val, x1.d_ * grad1 + x2.d_ * grad2); } template inline fvar log_add_exp(const fvar& x1, double x2) { - if (x1.val_ == NEGATIVE_INFTY) { - return fvar(x2, 0.0); // log_add_exp(-∞, b) = b - } - return log_add_exp(x2, x1); + if (x1.val_ == NEGATIVE_INFTY) { + return fvar(x2, 0.0); // log_add_exp(-∞, b) = b + } + return log_add_exp(x2, x1); } template inline fvar log_add_exp(double x1, const fvar& x2) { - if (x2.val_ == NEGATIVE_INFTY) { - return fvar(x1, 0.0); // log_add_exp(a, -∞) = a - } - auto val = stan::math::log_add_exp(x1, x2.val_); - auto exp_x2 = stan::math::exp(x2.val_); - auto grad = exp_x2 / (stan::math::exp(x1) + exp_x2); - return fvar(val, x2.d_ * grad); + if (x2.val_ == NEGATIVE_INFTY) { + return fvar(x1, 0.0); // log_add_exp(a, -∞) = a + } + auto val = stan::math::log_add_exp(x1, x2.val_); + auto exp_x2 = stan::math::exp(x2.val_); + auto grad = exp_x2 / (stan::math::exp(x1) + exp_x2); + return fvar(val, x2.d_ * grad); } // Overload for matrices of fvar template inline fvar log_add_exp(const Eigen::Matrix, -1, -1>& a, - const Eigen::Matrix, -1, -1>& b) { - using fvar_mat_type = Eigen::Matrix, -1, -1>; - fvar_mat_type result(a.rows(), a.cols()); - - // Check for empty inputs - if (a.size() == 0 || b.size() == 0) { - throw std::invalid_argument("Input containers must not be empty."); - } - - // Check for NaN - if (a.array().isNaN().any() || b.array().isNaN().any()) { - result.setConstant(fvar(std::numeric_limits::quiet_NaN())); - return result; - } - - // Check for infinity - if (a.array().isInf().any() || b.array().isInf().any()) { - result.setConstant(fvar(std::numeric_limits::quiet_NaN())); - return result; + const Eigen::Matrix, -1, -1>& b) { + using fvar_mat_type = Eigen::Matrix, -1, -1>; + fvar_mat_type result(a.rows(), a.cols()); + + // Check for empty inputs + if (a.size() == 0 || b.size() == 0) { + throw std::invalid_argument("Input containers must not be empty."); + } + + // Check for NaN + if (a.array().isNaN().any() || b.array().isNaN().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } + + // Check for infinity + if (a.array().isInf().any() || b.array().isInf().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } + + // Apply the log_add_exp operation directly + for (int i = 0; i < a.rows(); ++i) { + for (int j = 0; j < a.cols(); ++j) { + result(i, j) = stan::math::log_add_exp(a(i, j), b(i, j)); } + } - // Apply the log_add_exp operation directly - for (int i = 0; i < a.rows(); ++i) { - for (int j = 0; j < a.cols(); ++j) { - result(i, j) = stan::math::log_add_exp(a(i, j), b(i, j)); - } - } - - return result; // Return the result matrix + return result; // Return the result matrix } // Specialization for nested fvar types template -inline auto log_add_exp(const Eigen::Matrix>, -1, -1>& a, - const Eigen::Matrix>, -1, -1>& b) { - using nested_fvar_mat_type = Eigen::Matrix>, -1, -1>; - nested_fvar_mat_type result(a.rows(), a.cols()); - - // Check for empty inputs - if (a.size() == 0 || b.size() == 0) { - throw std::invalid_argument("Input containers must not be empty."); - } - - // Check for NaN - if (a.array().isNaN().any() || b.array().isNaN().any()) { - result.setConstant(stan::math::fvar>(std::numeric_limits::quiet_NaN())); - return result; - } - - // Check for infinity - if (a.array().isInf().any() || b.array().isInf().any()) { - result.setConstant(stan::math::fvar>(std::numeric_limits::quiet_NaN())); - return result; - } - - // Implement the logic for log_add_exp for nested fvar types - for (int i = 0; i < a.rows(); ++i) { - for (int j = 0; j < a.cols(); ++j) { - auto inner_a = a(i, j); - auto inner_b = b(i, j); - result(i, j) = stan::math::log_add_exp(inner_a, inner_b); - } +inline auto log_add_exp( + const Eigen::Matrix>, -1, -1>& a, + const Eigen::Matrix>, -1, -1>& + b) { + using nested_fvar_mat_type + = Eigen::Matrix>, -1, -1>; + nested_fvar_mat_type result(a.rows(), a.cols()); + + // Check for empty inputs + if (a.size() == 0 || b.size() == 0) { + throw std::invalid_argument("Input containers must not be empty."); + } + + // Check for NaN + if (a.array().isNaN().any() || b.array().isNaN().any()) { + result.setConstant(stan::math::fvar>( + std::numeric_limits::quiet_NaN())); + return result; + } + + // Check for infinity + if (a.array().isInf().any() || b.array().isInf().any()) { + result.setConstant(stan::math::fvar>( + std::numeric_limits::quiet_NaN())); + return result; + } + + // Implement the logic for log_add_exp for nested fvar types + for (int i = 0; i < a.rows(); ++i) { + for (int j = 0; j < a.cols(); ++j) { + auto inner_a = a(i, j); + auto inner_b = b(i, j); + result(i, j) = stan::math::log_add_exp(inner_a, inner_b); } + } - return result; // Return the result matrix + return result; // Return the result matrix } -} -} +} // namespace math +} // namespace stan #endif diff --git a/stan/math/prim/fun/log_add_exp.hpp b/stan/math/prim/fun/log_add_exp.hpp index 574c5dde1b1..7117d876179 100644 --- a/stan/math/prim/fun/log_add_exp.hpp +++ b/stan/math/prim/fun/log_add_exp.hpp @@ -26,8 +26,7 @@ namespace math { * @param b the second variable */ -template * = nullptr, +template * = nullptr, require_all_stan_scalar_t* = nullptr> inline return_type_t log_add_exp(const T2& a, const T1& b) { if (a == NEGATIVE_INFTY) { @@ -39,7 +38,7 @@ inline return_type_t log_add_exp(const T2& a, const T1& b) { if (a == INFTY || b == INFTY) { return INFTY; } - + const double max_val = std::max(a, b); return max_val + std::log(std::exp(a - max_val) + std::exp(b - max_val)); } @@ -57,41 +56,46 @@ inline return_type_t log_add_exp(const T2& a, const T1& b) { */ template * = nullptr> inline auto log_add_exp(const T& a, const T& b) { - if (a.size() != b.size()) { - throw std::invalid_argument("Binary function: size of x (" + std::to_string(a.size()) + ") and size of y (" + std::to_string(b.size()) + ") must match in size"); - } + if (a.size() != b.size()) { + throw std::invalid_argument("Binary function: size of x (" + + std::to_string(a.size()) + ") and size of y (" + + std::to_string(b.size()) + + ") must match in size"); + } - const size_t min_size = std::min(a.size(), b.size()); - using return_t = return_type_t; + const size_t min_size = std::min(a.size(), b.size()); + using return_t = return_type_t; - std::vector result(min_size); + std::vector result(min_size); - for (size_t i = 0; i < min_size; ++i) { - if (a[i] == NEGATIVE_INFTY) { - result[i] = b[i]; // log_add_exp(-∞, b) = b - } else if (b[i] == NEGATIVE_INFTY) { - result[i] = a[i]; // log_add_exp(a, -∞) = a - } else if (a[i] == INFTY || b[i] == INFTY) { - result[i] = INFTY; // log_add_exp(∞, b) = ∞ - } else { - // Log-add-exp trick - const double max_val = std::max(a[i], b[i]); - result[i] = max_val + std::log(std::exp(a[i] - max_val) + std::exp(b[i] - max_val)); - } + for (size_t i = 0; i < min_size; ++i) { + if (a[i] == NEGATIVE_INFTY) { + result[i] = b[i]; // log_add_exp(-∞, b) = b + } else if (b[i] == NEGATIVE_INFTY) { + result[i] = a[i]; // log_add_exp(a, -∞) = a + } else if (a[i] == INFTY || b[i] == INFTY) { + result[i] = INFTY; // log_add_exp(∞, b) = ∞ + } else { + // Log-add-exp trick + const double max_val = std::max(a[i], b[i]); + result[i] + = max_val + + std::log(std::exp(a[i] - max_val) + std::exp(b[i] - max_val)); } + } - return result; + return result; } /** * Enables the vectorized application of the log_add_exp function, * when the first and/or second arguments are containers. - * - * @tparam T1 - * @tparam T2 - * @param a - * @param b - * @return auto + * + * @tparam T1 + * @tparam T2 + * @param a + * @param b + * @return auto */ template * = nullptr> inline auto log_add_exp(const T1& a, const T2& b) { diff --git a/stan/math/rev/fun/log_add_exp.hpp b/stan/math/rev/fun/log_add_exp.hpp index eba7006e388..c0663e064a2 100644 --- a/stan/math/rev/fun/log_add_exp.hpp +++ b/stan/math/rev/fun/log_add_exp.hpp @@ -85,4 +85,4 @@ inline T log_add_exp(const T& x, const T& y) { } // namespace math } // namespace stan -#endif \ No newline at end of file +#endif \ No newline at end of file diff --git a/test/unit/math/mix/fun/log_add_exp_test.cpp b/test/unit/math/mix/fun/log_add_exp_test.cpp index de2c1e834dc..d311aa40034 100644 --- a/test/unit/math/mix/fun/log_add_exp_test.cpp +++ b/test/unit/math/mix/fun/log_add_exp_test.cpp @@ -15,7 +15,7 @@ TEST(MathMixMatFun, logAddExp) { stan::test::expect_ad(f, x1, y1); // Test with negative infinity - + stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); @@ -65,7 +65,7 @@ TEST(MathMixMatFun, log_add_exp_mismatched_sizes) { std::vector x{1.0, 2.0}; std::vector y{1.0, 2.0, 3.0}; - + stan::test::expect_ad(f, x, y); stan::test::expect_ad(f, y, x); } From 6d39d20ba08aa5a1e8cd80a662b5c64f043d45c7 Mon Sep 17 00:00:00 2001 From: michascant <89426143+MichaScant@users.noreply.github.com> Date: Mon, 25 Nov 2024 12:44:11 -0500 Subject: [PATCH 3/7] correct indentation in a few areas for improved readability --- stan/math/fwd/fun/log_add_exp.hpp | 1 + stan/math/prim/fun/log_add_exp.hpp | 27 ++++--- stan/math/rev/fun/log_add_exp.hpp | 51 ++++++------ test/unit/math/mix/fun/log_add_exp_test.cpp | 88 ++++++++++----------- 4 files changed, 85 insertions(+), 82 deletions(-) diff --git a/stan/math/fwd/fun/log_add_exp.hpp b/stan/math/fwd/fun/log_add_exp.hpp index 99cd0284bf7..31e2ba0543e 100644 --- a/stan/math/fwd/fun/log_add_exp.hpp +++ b/stan/math/fwd/fun/log_add_exp.hpp @@ -51,6 +51,7 @@ inline fvar log_add_exp(double x1, const fvar& x2) { template inline fvar log_add_exp(const Eigen::Matrix, -1, -1>& a, const Eigen::Matrix, -1, -1>& b) { + using fvar_mat_type = Eigen::Matrix, -1, -1>; fvar_mat_type result(a.rows(), a.cols()); diff --git a/stan/math/prim/fun/log_add_exp.hpp b/stan/math/prim/fun/log_add_exp.hpp index 574c5dde1b1..e47cfca16f3 100644 --- a/stan/math/prim/fun/log_add_exp.hpp +++ b/stan/math/prim/fun/log_add_exp.hpp @@ -27,21 +27,21 @@ namespace math { */ template * = nullptr, - require_all_stan_scalar_t* = nullptr> + require_all_not_st_var* = nullptr, + require_all_stan_scalar_t* = nullptr> inline return_type_t log_add_exp(const T2& a, const T1& b) { - if (a == NEGATIVE_INFTY) { + if (a == NEGATIVE_INFTY) { return b; - } - if (b == NEGATIVE_INFTY) { + } + if (b == NEGATIVE_INFTY) { return a; - } - if (a == INFTY || b == INFTY) { + } + if (a == INFTY || b == INFTY) { return INFTY; - } - - const double max_val = std::max(a, b); - return max_val + std::log(std::exp(a - max_val) + std::exp(b - max_val)); + } + + const double max_val = std::max(a, b); + return max_val + std::log(std::exp(a - max_val) + std::exp(b - max_val)); } /** @@ -95,8 +95,9 @@ inline auto log_add_exp(const T& a, const T& b) { */ template * = nullptr> inline auto log_add_exp(const T1& a, const T2& b) { - return apply_scalar_binary( - a, b, [](const auto& c, const auto& d) { return log_sum_exp(c, d); }); + return apply_scalar_binary( + a, b, [](const auto& c, const auto& d) { return log_sum_exp(c, d); } + ); } } // namespace math diff --git a/stan/math/rev/fun/log_add_exp.hpp b/stan/math/rev/fun/log_add_exp.hpp index eba7006e388..38954e8f6b0 100644 --- a/stan/math/rev/fun/log_add_exp.hpp +++ b/stan/math/rev/fun/log_add_exp.hpp @@ -19,31 +19,31 @@ namespace math { namespace internal { class log_add_exp_vv_vari : public op_vv_vari { - public: - log_add_exp_vv_vari(vari* avi, vari* bvi) - : op_vv_vari(log_add_exp(avi->val_, bvi->val_), avi, bvi) {} - void chain() { - double exp_a = std::exp(avi_->val_); - double exp_b = std::exp(bvi_->val_); - double sum_exp = exp_a + exp_b; + public: + log_add_exp_vv_vari(vari* avi, vari* bvi) + : op_vv_vari(log_add_exp(avi->val_, bvi->val_), avi, bvi) {} + void chain() { + double exp_a = std::exp(avi_->val_); + double exp_b = std::exp(bvi_->val_); + double sum_exp = exp_a + exp_b; - avi_->adj_ += adj_ * (exp_a / sum_exp); - bvi_->adj_ += adj_ * (exp_b / sum_exp); - } + avi_->adj_ += adj_ * (exp_a / sum_exp); + bvi_->adj_ += adj_ * (exp_b / sum_exp); + } }; class log_add_exp_vd_vari : public op_vd_vari { - public: - log_add_exp_vd_vari(vari* avi, double b) - : op_vd_vari(log_add_exp(avi->val_, b), avi, b) {} - void chain() { - if (val_ == NEGATIVE_INFTY) { - avi_->adj_ += adj_; - } else { - double exp_a = std::exp(avi_->val_); - avi_->adj_ += adj_ * (exp_a / (exp_a + std::exp(bd_))); + public: + log_add_exp_vd_vari(vari* avi, double b) + : op_vd_vari(log_add_exp(avi->val_, b), avi, b) {} + void chain() { + if (val_ == NEGATIVE_INFTY) { + avi_->adj_ += adj_; + } else { + double exp_a = std::exp(avi_->val_); + avi_->adj_ += adj_ * (exp_a / (exp_a + std::exp(bd_))); + } } - } }; } // namespace internal @@ -52,21 +52,21 @@ class log_add_exp_vd_vari : public op_vd_vari { * Returns the element-wise log sum of exponentials. */ inline var log_add_exp(const var& a, const var& b) { - return var(new internal::log_add_exp_vv_vari(a.vi_, b.vi_)); + return var(new internal::log_add_exp_vv_vari(a.vi_, b.vi_)); } /** * Returns the log sum of exponentials. */ inline var log_add_exp(const var& a, double b) { - return var(new internal::log_add_exp_vd_vari(a.vi_, b)); + return var(new internal::log_add_exp_vd_vari(a.vi_, b)); } /** * Returns the element-wise log sum of exponentials. */ inline var log_add_exp(double a, const var& b) { - return var(new internal::log_add_exp_vd_vari(b.vi_, a)); + return var(new internal::log_add_exp_vd_vari(b.vi_, a)); } /** @@ -78,8 +78,9 @@ inline var log_add_exp(double a, const var& b) { */ template * = nullptr> inline T log_add_exp(const T& x, const T& y) { - return apply_scalar_binary( - x, y, [](const auto& a, const auto& b) { return log_add_exp(a, b); }); + return apply_scalar_binary( + x, y, [](const auto& a, const auto& b) { return log_add_exp(a, b); } + ); } } // namespace math diff --git a/test/unit/math/mix/fun/log_add_exp_test.cpp b/test/unit/math/mix/fun/log_add_exp_test.cpp index de2c1e834dc..6db2f3f8067 100644 --- a/test/unit/math/mix/fun/log_add_exp_test.cpp +++ b/test/unit/math/mix/fun/log_add_exp_test.cpp @@ -4,68 +4,68 @@ #include TEST(MathMixMatFun, logAddExp) { - auto f = [](const auto& x, const auto& y) { + auto f = [](const auto& x, const auto& y) { return stan::math::log_add_exp(x, y); - }; - // Test with finite values - Eigen::VectorXd x1(2); - x1 << 2.0, 1.0; - Eigen::VectorXd y1(2); - y1 << 3.0, 2.0; - stan::test::expect_ad(f, x1, y1); + }; + // Test with finite values + Eigen::VectorXd x1(2); + x1 << 2.0, 1.0; + Eigen::VectorXd y1(2); + y1 << 3.0, 2.0; + stan::test::expect_ad(f, x1, y1); - // Test with negative infinity - - stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); - stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); + // Test with negative infinity - // Test with infinity - stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); + stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); + stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); + + // Test with infinity + stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); } TEST(MathMixMatFun, log_add_exp_elementwise_values) { - auto f = [](const auto& x, const auto& y) { + auto f = [](const auto& x, const auto& y) { return stan::math::log_add_exp(x, y); - }; + }; - Eigen::VectorXd x1(2); - x1 << 2.0, 1.0; - Eigen::VectorXd y1(2); - y1 << 3.0, 2.0; - stan::test::expect_ad(f, x1, y1); + Eigen::VectorXd x1(2); + x1 << 2.0, 1.0; + Eigen::VectorXd y1(2); + y1 << 3.0, 2.0; + stan::test::expect_ad(f, x1, y1); - Eigen::VectorXd x2(2); - x2 << 0.5, -1.0; - Eigen::VectorXd y2(2); - y2 << 1.0, 2.0; - stan::test::expect_ad(f, x2, y2); + Eigen::VectorXd x2(2); + x2 << 0.5, -1.0; + Eigen::VectorXd y2(2); + y2 << 1.0, 2.0; + stan::test::expect_ad(f, x2, y2); - // Test with infinity - Eigen::VectorXd x3(2); - x3 << std::numeric_limits::infinity(), 1.0; - Eigen::VectorXd y3(2); - y3 << 2.0, std::numeric_limits::infinity(); - stan::test::expect_ad(f, x3, y3); + // Test with infinity + Eigen::VectorXd x3(2); + x3 << std::numeric_limits::infinity(), 1.0; + Eigen::VectorXd y3(2); + y3 << 2.0, std::numeric_limits::infinity(); + stan::test::expect_ad(f, x3, y3); } TEST(MathMixMatFun, log_add_exp_edge_cases) { - auto f = [](const auto& x, const auto& y) { + auto f = [](const auto& x, const auto& y) { return stan::math::log_add_exp(x, y); - }; + }; - stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); - stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); - stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); + stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); + stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); + stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); } TEST(MathMixMatFun, log_add_exp_mismatched_sizes) { - auto f = [](const auto& x, const auto& y) { + auto f = [](const auto& x, const auto& y) { return stan::math::log_add_exp(x, y); - }; + }; + + std::vector x{1.0, 2.0}; + std::vector y{1.0, 2.0, 3.0}; - std::vector x{1.0, 2.0}; - std::vector y{1.0, 2.0, 3.0}; - - stan::test::expect_ad(f, x, y); - stan::test::expect_ad(f, y, x); + stan::test::expect_ad(f, x, y); + stan::test::expect_ad(f, y, x); } From 0ef959683efda3bb6869f04bbf86e8a7a23a8683 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 25 Nov 2024 12:47:50 -0500 Subject: [PATCH 4/7] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/fwd/fun/log_add_exp.hpp | 36 ++++----- stan/math/prim/fun/log_add_exp.hpp | 26 +++---- stan/math/rev/fun/log_add_exp.hpp | 51 ++++++------- test/unit/math/mix/fun/log_add_exp_test.cpp | 84 ++++++++++----------- 4 files changed, 97 insertions(+), 100 deletions(-) diff --git a/stan/math/fwd/fun/log_add_exp.hpp b/stan/math/fwd/fun/log_add_exp.hpp index 26fcff8068c..f6c0827b1a0 100644 --- a/stan/math/fwd/fun/log_add_exp.hpp +++ b/stan/math/fwd/fun/log_add_exp.hpp @@ -49,36 +49,36 @@ inline fvar log_add_exp(double x1, const fvar& x2) { // Overload for matrices of fvar template -inline fvar log_add_exp(const Eigen::Matrix, -1, -1>& a, const Eigen::Matrix, -1, -1>& b) { - - using fvar_mat_type = Eigen::Matrix, -1, -1>; - fvar_mat_type result(a.rows(), a.cols()); +inline fvar log_add_exp(const Eigen::Matrix, -1, -1>& a, + const Eigen::Matrix, -1, -1>& b) { + using fvar_mat_type = Eigen::Matrix, -1, -1>; + fvar_mat_type result(a.rows(), a.cols()); - // Check for empty inputs - if (a.size() == 0 || b.size() == 0) { + // Check for empty inputs + if (a.size() == 0 || b.size() == 0) { throw std::invalid_argument("Input containers must not be empty."); - } + } - // Check for NaN - if (a.array().isNaN().any() || b.array().isNaN().any()) { + // Check for NaN + if (a.array().isNaN().any() || b.array().isNaN().any()) { result.setConstant(fvar(std::numeric_limits::quiet_NaN())); return result; - } + } - // Check for infinity - if (a.array().isInf().any() || b.array().isInf().any()) { + // Check for infinity + if (a.array().isInf().any() || b.array().isInf().any()) { result.setConstant(fvar(std::numeric_limits::quiet_NaN())); return result; - } + } - // Apply the log_add_exp operation directly - for (int i = 0; i < a.rows(); ++i) { + // Apply the log_add_exp operation directly + for (int i = 0; i < a.rows(); ++i) { for (int j = 0; j < a.cols(); ++j) { - result(i, j) = stan::math::log_add_exp(a(i, j), b(i, j)); - } + result(i, j) = stan::math::log_add_exp(a(i, j), b(i, j)); } + } - return result; // Return the result matrix + return result; // Return the result matrix } // Specialization for nested fvar types diff --git a/stan/math/prim/fun/log_add_exp.hpp b/stan/math/prim/fun/log_add_exp.hpp index 82aa9a83fff..7117d876179 100644 --- a/stan/math/prim/fun/log_add_exp.hpp +++ b/stan/math/prim/fun/log_add_exp.hpp @@ -26,22 +26,21 @@ namespace math { * @param b the second variable */ -template * = nullptr, - require_all_stan_scalar_t* = nullptr> +template * = nullptr, + require_all_stan_scalar_t* = nullptr> inline return_type_t log_add_exp(const T2& a, const T1& b) { - if (a == NEGATIVE_INFTY) { + if (a == NEGATIVE_INFTY) { return b; - } - if (b == NEGATIVE_INFTY) { + } + if (b == NEGATIVE_INFTY) { return a; - } - if (a == INFTY || b == INFTY) { + } + if (a == INFTY || b == INFTY) { return INFTY; - } + } - const double max_val = std::max(a, b); - return max_val + std::log(std::exp(a - max_val) + std::exp(b - max_val)); + const double max_val = std::max(a, b); + return max_val + std::log(std::exp(a - max_val) + std::exp(b - max_val)); } /** @@ -100,9 +99,8 @@ inline auto log_add_exp(const T& a, const T& b) { */ template * = nullptr> inline auto log_add_exp(const T1& a, const T2& b) { - return apply_scalar_binary( - a, b, [](const auto& c, const auto& d) { return log_sum_exp(c, d); } - ); + return apply_scalar_binary( + a, b, [](const auto& c, const auto& d) { return log_sum_exp(c, d); }); } } // namespace math diff --git a/stan/math/rev/fun/log_add_exp.hpp b/stan/math/rev/fun/log_add_exp.hpp index be870405e69..c0663e064a2 100644 --- a/stan/math/rev/fun/log_add_exp.hpp +++ b/stan/math/rev/fun/log_add_exp.hpp @@ -19,31 +19,31 @@ namespace math { namespace internal { class log_add_exp_vv_vari : public op_vv_vari { - public: - log_add_exp_vv_vari(vari* avi, vari* bvi) - : op_vv_vari(log_add_exp(avi->val_, bvi->val_), avi, bvi) {} - void chain() { - double exp_a = std::exp(avi_->val_); - double exp_b = std::exp(bvi_->val_); - double sum_exp = exp_a + exp_b; + public: + log_add_exp_vv_vari(vari* avi, vari* bvi) + : op_vv_vari(log_add_exp(avi->val_, bvi->val_), avi, bvi) {} + void chain() { + double exp_a = std::exp(avi_->val_); + double exp_b = std::exp(bvi_->val_); + double sum_exp = exp_a + exp_b; - avi_->adj_ += adj_ * (exp_a / sum_exp); - bvi_->adj_ += adj_ * (exp_b / sum_exp); - } + avi_->adj_ += adj_ * (exp_a / sum_exp); + bvi_->adj_ += adj_ * (exp_b / sum_exp); + } }; class log_add_exp_vd_vari : public op_vd_vari { - public: - log_add_exp_vd_vari(vari* avi, double b) - : op_vd_vari(log_add_exp(avi->val_, b), avi, b) {} - void chain() { - if (val_ == NEGATIVE_INFTY) { - avi_->adj_ += adj_; - } else { - double exp_a = std::exp(avi_->val_); - avi_->adj_ += adj_ * (exp_a / (exp_a + std::exp(bd_))); - } + public: + log_add_exp_vd_vari(vari* avi, double b) + : op_vd_vari(log_add_exp(avi->val_, b), avi, b) {} + void chain() { + if (val_ == NEGATIVE_INFTY) { + avi_->adj_ += adj_; + } else { + double exp_a = std::exp(avi_->val_); + avi_->adj_ += adj_ * (exp_a / (exp_a + std::exp(bd_))); } + } }; } // namespace internal @@ -52,21 +52,21 @@ class log_add_exp_vd_vari : public op_vd_vari { * Returns the element-wise log sum of exponentials. */ inline var log_add_exp(const var& a, const var& b) { - return var(new internal::log_add_exp_vv_vari(a.vi_, b.vi_)); + return var(new internal::log_add_exp_vv_vari(a.vi_, b.vi_)); } /** * Returns the log sum of exponentials. */ inline var log_add_exp(const var& a, double b) { - return var(new internal::log_add_exp_vd_vari(a.vi_, b)); + return var(new internal::log_add_exp_vd_vari(a.vi_, b)); } /** * Returns the element-wise log sum of exponentials. */ inline var log_add_exp(double a, const var& b) { - return var(new internal::log_add_exp_vd_vari(b.vi_, a)); + return var(new internal::log_add_exp_vd_vari(b.vi_, a)); } /** @@ -78,9 +78,8 @@ inline var log_add_exp(double a, const var& b) { */ template * = nullptr> inline T log_add_exp(const T& x, const T& y) { - return apply_scalar_binary( - x, y, [](const auto& a, const auto& b) { return log_add_exp(a, b); } - ); + return apply_scalar_binary( + x, y, [](const auto& a, const auto& b) { return log_add_exp(a, b); }); } } // namespace math diff --git a/test/unit/math/mix/fun/log_add_exp_test.cpp b/test/unit/math/mix/fun/log_add_exp_test.cpp index 6db2f3f8067..d311aa40034 100644 --- a/test/unit/math/mix/fun/log_add_exp_test.cpp +++ b/test/unit/math/mix/fun/log_add_exp_test.cpp @@ -4,68 +4,68 @@ #include TEST(MathMixMatFun, logAddExp) { - auto f = [](const auto& x, const auto& y) { + auto f = [](const auto& x, const auto& y) { return stan::math::log_add_exp(x, y); - }; - // Test with finite values - Eigen::VectorXd x1(2); - x1 << 2.0, 1.0; - Eigen::VectorXd y1(2); - y1 << 3.0, 2.0; - stan::test::expect_ad(f, x1, y1); + }; + // Test with finite values + Eigen::VectorXd x1(2); + x1 << 2.0, 1.0; + Eigen::VectorXd y1(2); + y1 << 3.0, 2.0; + stan::test::expect_ad(f, x1, y1); - // Test with negative infinity + // Test with negative infinity - stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); - stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); + stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); + stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); - // Test with infinity - stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); + // Test with infinity + stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); } TEST(MathMixMatFun, log_add_exp_elementwise_values) { - auto f = [](const auto& x, const auto& y) { + auto f = [](const auto& x, const auto& y) { return stan::math::log_add_exp(x, y); - }; + }; - Eigen::VectorXd x1(2); - x1 << 2.0, 1.0; - Eigen::VectorXd y1(2); - y1 << 3.0, 2.0; - stan::test::expect_ad(f, x1, y1); + Eigen::VectorXd x1(2); + x1 << 2.0, 1.0; + Eigen::VectorXd y1(2); + y1 << 3.0, 2.0; + stan::test::expect_ad(f, x1, y1); - Eigen::VectorXd x2(2); - x2 << 0.5, -1.0; - Eigen::VectorXd y2(2); - y2 << 1.0, 2.0; - stan::test::expect_ad(f, x2, y2); + Eigen::VectorXd x2(2); + x2 << 0.5, -1.0; + Eigen::VectorXd y2(2); + y2 << 1.0, 2.0; + stan::test::expect_ad(f, x2, y2); - // Test with infinity - Eigen::VectorXd x3(2); - x3 << std::numeric_limits::infinity(), 1.0; - Eigen::VectorXd y3(2); - y3 << 2.0, std::numeric_limits::infinity(); - stan::test::expect_ad(f, x3, y3); + // Test with infinity + Eigen::VectorXd x3(2); + x3 << std::numeric_limits::infinity(), 1.0; + Eigen::VectorXd y3(2); + y3 << 2.0, std::numeric_limits::infinity(); + stan::test::expect_ad(f, x3, y3); } TEST(MathMixMatFun, log_add_exp_edge_cases) { - auto f = [](const auto& x, const auto& y) { + auto f = [](const auto& x, const auto& y) { return stan::math::log_add_exp(x, y); - }; + }; - stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); - stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); - stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); + stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); + stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); + stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); } TEST(MathMixMatFun, log_add_exp_mismatched_sizes) { - auto f = [](const auto& x, const auto& y) { + auto f = [](const auto& x, const auto& y) { return stan::math::log_add_exp(x, y); - }; + }; - std::vector x{1.0, 2.0}; - std::vector y{1.0, 2.0, 3.0}; + std::vector x{1.0, 2.0}; + std::vector y{1.0, 2.0, 3.0}; - stan::test::expect_ad(f, x, y); - stan::test::expect_ad(f, y, x); + stan::test::expect_ad(f, x, y); + stan::test::expect_ad(f, y, x); } From 1cfdfae32484e10ca50381b9e503145209cb11df Mon Sep 17 00:00:00 2001 From: michascant <89426143+MichaScant@users.noreply.github.com> Date: Mon, 25 Nov 2024 13:02:02 -0500 Subject: [PATCH 5/7] added trailing whitespace --- stan/math/rev/fun/log_add_exp.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/log_add_exp.hpp b/stan/math/rev/fun/log_add_exp.hpp index c0663e064a2..cfe7711b796 100644 --- a/stan/math/rev/fun/log_add_exp.hpp +++ b/stan/math/rev/fun/log_add_exp.hpp @@ -85,4 +85,4 @@ inline T log_add_exp(const T& x, const T& y) { } // namespace math } // namespace stan -#endif \ No newline at end of file +#endif From 1595abab9df0b6309cec13c86d22f4c29985be1e Mon Sep 17 00:00:00 2001 From: michascant <89426143+MichaScant@users.noreply.github.com> Date: Thu, 28 Nov 2024 14:21:48 -0500 Subject: [PATCH 6/7] incorperating log_sum_exp into log_add_exp, added functionality for matrices --- stan/math/fwd/fun/log_add_exp.hpp | 80 +++++++++++----- stan/math/prim/fun/log_add_exp.hpp | 101 ++++++++++++++------ test/unit/math/mix/fun/log_add_exp_test.cpp | 39 +++++--- 3 files changed, 157 insertions(+), 63 deletions(-) diff --git a/stan/math/fwd/fun/log_add_exp.hpp b/stan/math/fwd/fun/log_add_exp.hpp index f6c0827b1a0..57c8e68cf6f 100644 --- a/stan/math/fwd/fun/log_add_exp.hpp +++ b/stan/math/fwd/fun/log_add_exp.hpp @@ -49,36 +49,68 @@ inline fvar log_add_exp(double x1, const fvar& x2) { // Overload for matrices of fvar template -inline fvar log_add_exp(const Eigen::Matrix, -1, -1>& a, - const Eigen::Matrix, -1, -1>& b) { - using fvar_mat_type = Eigen::Matrix, -1, -1>; - fvar_mat_type result(a.rows(), a.cols()); +inline Eigen::Matrix, -1, -1> log_add_exp(const Eigen::Matrix, -1, -1>& a, + const Eigen::Matrix, -1, -1>& b) { + using fvar_mat_type = Eigen::Matrix, -1, -1>; + fvar_mat_type result(a.rows(), a.cols()); + + // Check for empty inputs + if (a.size() == 0 || b.size() == 0) { + throw std::invalid_argument("Input containers must not be empty."); + } - // Check for empty inputs - if (a.size() == 0 || b.size() == 0) { - throw std::invalid_argument("Input containers must not be empty."); - } + // Check for NaN + if (a.array().isNaN().any() || b.array().isNaN().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } - // Check for NaN - if (a.array().isNaN().any() || b.array().isNaN().any()) { - result.setConstant(fvar(std::numeric_limits::quiet_NaN())); - return result; - } + // Check for infinity + if (a.array().isInf().any() || b.array().isInf().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } - // Check for infinity - if (a.array().isInf().any() || b.array().isInf().any()) { - result.setConstant(fvar(std::numeric_limits::quiet_NaN())); - return result; - } + // Apply the log_add_exp operation directly + for (int i = 0; i < a.rows(); ++i) { + for (int j = 0; j < a.cols(); ++j) { + result(i, j) = stan::math::log_add_exp(a(i, j), b(i, j)); + } + } - // Apply the log_add_exp operation directly - for (int i = 0; i < a.rows(); ++i) { - for (int j = 0; j < a.cols(); ++j) { - result(i, j) = stan::math::log_add_exp(a(i, j), b(i, j)); + return result; // Return the result matrix +} + +// Overload for Eigen vectors +template +inline Eigen::Matrix, -1, 1> log_add_exp(const Eigen::Matrix, -1, 1>& a, + const Eigen::Matrix, -1, 1>& b) { + using fvar_vec_type = Eigen::Matrix, -1, 1>; + fvar_vec_type result(a.rows()); + + // Check for empty inputs + if (a.size() == 0 || b.size() == 0) { + throw std::invalid_argument("Input containers must not be empty."); } - } - return result; // Return the result matrix + // Check for NaN + if (a.array().isNaN().any() || b.array().isNaN().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } + + // Check for infinity + if (a.array().isInf().any() || b.array().isInf().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } + + // Apply the log_add_exp operation directly + for (int i = 0; i < a.rows(); ++i) { + result(i) = stan::math::log_add_exp(a(i), b(i)); + } + + return result; // Return the result vector } // Specialization for nested fvar types diff --git a/stan/math/prim/fun/log_add_exp.hpp b/stan/math/prim/fun/log_add_exp.hpp index 7117d876179..eeaf2d49612 100644 --- a/stan/math/prim/fun/log_add_exp.hpp +++ b/stan/math/prim/fun/log_add_exp.hpp @@ -5,10 +5,13 @@ #include #include #include +#include #include #include #include #include +#include +#include namespace stan { namespace math { @@ -29,6 +32,7 @@ namespace math { template * = nullptr, require_all_stan_scalar_t* = nullptr> inline return_type_t log_add_exp(const T2& a, const T1& b) { + if (a == NEGATIVE_INFTY) { return b; } @@ -56,35 +60,66 @@ inline return_type_t log_add_exp(const T2& a, const T1& b) { */ template * = nullptr> inline auto log_add_exp(const T& a, const T& b) { - if (a.size() != b.size()) { - throw std::invalid_argument("Binary function: size of x (" - + std::to_string(a.size()) + ") and size of y (" - + std::to_string(b.size()) - + ") must match in size"); - } - const size_t min_size = std::min(a.size(), b.size()); - using return_t = return_type_t; + // Check if sizes are compatible + if constexpr (stan::is_eigen::value) { + // Check if both matrices/vectors have the same dimensions + stan::math::check_matching_dims("log_add_exp", "a", a, "b", b); + + // Determine the number of rows and columns for the result + size_t rows = a.rows(); + size_t cols = b.cols(); + using return_t = return_type_t; + + Eigen::Matrix result(rows, cols); + + // Iterate over each element + for (size_t i = 0; i < rows; ++i) { + for (size_t j = 0; j < cols; ++j) { + double a_val = (a.cols() == 1) ? a(i, 0) : a(i, j); // Handle column vector or matrix + double b_val = (b.rows() == 1) ? b(0, j) : b(i, j); // Handle row vector or matrix + + if (a_val == NEGATIVE_INFTY) { + result(i, j) = b_val; + } else if (b_val == NEGATIVE_INFTY) { + result(i, j) = a_val; + } else if (a_val == INFTY || b_val == INFTY) { + result(i, j) = INFTY; + } else { + result(i, j) = log_sum_exp(a_val, b_val); + } + } + } + + return result; + } else if constexpr (std::is_same_v>) { + // Handle std::vector + if (a.size() != b.size()) { + throw std::invalid_argument("Sizes of x and y must match."); + } - std::vector result(min_size); + using return_t = return_type_t; + std::vector result(a.size()); - for (size_t i = 0; i < min_size; ++i) { - if (a[i] == NEGATIVE_INFTY) { - result[i] = b[i]; // log_add_exp(-∞, b) = b - } else if (b[i] == NEGATIVE_INFTY) { - result[i] = a[i]; // log_add_exp(a, -∞) = a - } else if (a[i] == INFTY || b[i] == INFTY) { - result[i] = INFTY; // log_add_exp(∞, b) = ∞ + for (size_t i = 0; i < a.size(); ++i) { + double a_val = a[i]; + double b_val = b[i]; + + if (a_val == NEGATIVE_INFTY) { + result[i] = b_val; + } else if (b_val == NEGATIVE_INFTY) { + result[i] = a_val; + } else if (a_val == INFTY || b_val == INFTY) { + result[i] = INFTY; + } else { + result[i] = log_sum_exp(a_val, b_val); + } + } + + return result; } else { - // Log-add-exp trick - const double max_val = std::max(a[i], b[i]); - result[i] - = max_val - + std::log(std::exp(a[i] - max_val) + std::exp(b[i] - max_val)); + throw std::invalid_argument("Unsupported container type."); } - } - - return result; } /** @@ -99,11 +134,23 @@ inline auto log_add_exp(const T& a, const T& b) { */ template * = nullptr> inline auto log_add_exp(const T1& a, const T2& b) { - return apply_scalar_binary( - a, b, [](const auto& c, const auto& d) { return log_sum_exp(c, d); }); + // Check if both are Eigen/vectors + if constexpr (stan::is_eigen::value && stan::is_eigen::value) { + // Check if both matrices/vectors have the same dimensions + stan::math::check_matching_dims("log_add_exp", "a", a, "b", b); + } else { + // Check if sizes are compatible for other types + if (a.size() != b.size()) { + throw std::invalid_argument("Sizes of x and y must match or be compatible."); + } + } + + // If dimensions are verified to match, apply the operation + return apply_scalar_binary( + a, b, [](const auto& c, const auto& d) { return log_add_exp(c, d); }); } } // namespace math -} // namespace stan +} // nfamespace stan #endif diff --git a/test/unit/math/mix/fun/log_add_exp_test.cpp b/test/unit/math/mix/fun/log_add_exp_test.cpp index d311aa40034..024b2f331e0 100644 --- a/test/unit/math/mix/fun/log_add_exp_test.cpp +++ b/test/unit/math/mix/fun/log_add_exp_test.cpp @@ -1,7 +1,6 @@ #include #include #include -#include TEST(MathMixMatFun, logAddExp) { auto f = [](const auto& x, const auto& y) { @@ -14,7 +13,7 @@ TEST(MathMixMatFun, logAddExp) { y1 << 3.0, 2.0; stan::test::expect_ad(f, x1, y1); - // Test with negative infinity + //Test with negative infinity stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); @@ -45,19 +44,13 @@ TEST(MathMixMatFun, log_add_exp_elementwise_values) { x3 << std::numeric_limits::infinity(), 1.0; Eigen::VectorXd y3(2); y3 << 2.0, std::numeric_limits::infinity(); - stan::test::expect_ad(f, x3, y3); -} - -TEST(MathMixMatFun, log_add_exp_edge_cases) { - auto f = [](const auto& x, const auto& y) { - return stan::math::log_add_exp(x, y); - }; - stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); - stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); - stan::test::expect_ad(f, stan::math::INFTY, stan::math::INFTY); + Eigen::VectorXd result = f(x3, y3); + EXPECT_TRUE(std::isinf(result[0])); // Expect infinity for the first element + EXPECT_TRUE(std::isinf(result[1])); } + TEST(MathMixMatFun, log_add_exp_mismatched_sizes) { auto f = [](const auto& x, const auto& y) { return stan::math::log_add_exp(x, y); @@ -69,3 +62,25 @@ TEST(MathMixMatFun, log_add_exp_mismatched_sizes) { stan::test::expect_ad(f, x, y); stan::test::expect_ad(f, y, x); } + +TEST(MathMixMatFun, log_add_exp_container_tests) { + auto f = [](const auto& x, const auto& y) { + return stan::math::log_add_exp(x, y); + }; + + // Test with Eigen::MatrixXd + Eigen::MatrixXd x_row(1, 2); + x_row << 2.0, 1.0; + Eigen::MatrixXd y_row(1, 2); + y_row << 3.0, 2.0; + + stan::test::expect_ad(f, x_row, y_row); + + // Additional tests with mismatched sizes + Eigen::MatrixXd x_mismatch(2, 1); + x_mismatch << 0.5, -1.0; + Eigen::MatrixXd y_mismatch(1, 3); + y_mismatch << 1.0, 2.0, 3.0; + + EXPECT_THROW(stan::math::log_add_exp(x_mismatch, y_mismatch), std::invalid_argument); +} From 88b378faefd4dbcf251a64b283dd712d9423e425 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 28 Nov 2024 14:46:41 -0500 Subject: [PATCH 7/7] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/fwd/fun/log_add_exp.hpp | 100 +++++++------- stan/math/prim/fun/log_add_exp.hpp | 141 ++++++++++---------- test/unit/math/mix/fun/log_add_exp_test.cpp | 12 +- 3 files changed, 129 insertions(+), 124 deletions(-) diff --git a/stan/math/fwd/fun/log_add_exp.hpp b/stan/math/fwd/fun/log_add_exp.hpp index 57c8e68cf6f..2a446a301cb 100644 --- a/stan/math/fwd/fun/log_add_exp.hpp +++ b/stan/math/fwd/fun/log_add_exp.hpp @@ -49,68 +49,70 @@ inline fvar log_add_exp(double x1, const fvar& x2) { // Overload for matrices of fvar template -inline Eigen::Matrix, -1, -1> log_add_exp(const Eigen::Matrix, -1, -1>& a, - const Eigen::Matrix, -1, -1>& b) { - using fvar_mat_type = Eigen::Matrix, -1, -1>; - fvar_mat_type result(a.rows(), a.cols()); - - // Check for empty inputs - if (a.size() == 0 || b.size() == 0) { - throw std::invalid_argument("Input containers must not be empty."); - } +inline Eigen::Matrix, -1, -1> log_add_exp( + const Eigen::Matrix, -1, -1>& a, + const Eigen::Matrix, -1, -1>& b) { + using fvar_mat_type = Eigen::Matrix, -1, -1>; + fvar_mat_type result(a.rows(), a.cols()); - // Check for NaN - if (a.array().isNaN().any() || b.array().isNaN().any()) { - result.setConstant(fvar(std::numeric_limits::quiet_NaN())); - return result; - } + // Check for empty inputs + if (a.size() == 0 || b.size() == 0) { + throw std::invalid_argument("Input containers must not be empty."); + } - // Check for infinity - if (a.array().isInf().any() || b.array().isInf().any()) { - result.setConstant(fvar(std::numeric_limits::quiet_NaN())); - return result; - } + // Check for NaN + if (a.array().isNaN().any() || b.array().isNaN().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } + + // Check for infinity + if (a.array().isInf().any() || b.array().isInf().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } - // Apply the log_add_exp operation directly - for (int i = 0; i < a.rows(); ++i) { - for (int j = 0; j < a.cols(); ++j) { - result(i, j) = stan::math::log_add_exp(a(i, j), b(i, j)); - } + // Apply the log_add_exp operation directly + for (int i = 0; i < a.rows(); ++i) { + for (int j = 0; j < a.cols(); ++j) { + result(i, j) = stan::math::log_add_exp(a(i, j), b(i, j)); } + } - return result; // Return the result matrix + return result; // Return the result matrix } // Overload for Eigen vectors template -inline Eigen::Matrix, -1, 1> log_add_exp(const Eigen::Matrix, -1, 1>& a, - const Eigen::Matrix, -1, 1>& b) { - using fvar_vec_type = Eigen::Matrix, -1, 1>; - fvar_vec_type result(a.rows()); - - // Check for empty inputs - if (a.size() == 0 || b.size() == 0) { - throw std::invalid_argument("Input containers must not be empty."); - } +inline Eigen::Matrix, -1, 1> log_add_exp( + const Eigen::Matrix, -1, 1>& a, + const Eigen::Matrix, -1, 1>& b) { + using fvar_vec_type = Eigen::Matrix, -1, 1>; + fvar_vec_type result(a.rows()); - // Check for NaN - if (a.array().isNaN().any() || b.array().isNaN().any()) { - result.setConstant(fvar(std::numeric_limits::quiet_NaN())); - return result; - } + // Check for empty inputs + if (a.size() == 0 || b.size() == 0) { + throw std::invalid_argument("Input containers must not be empty."); + } - // Check for infinity - if (a.array().isInf().any() || b.array().isInf().any()) { - result.setConstant(fvar(std::numeric_limits::quiet_NaN())); - return result; - } + // Check for NaN + if (a.array().isNaN().any() || b.array().isNaN().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } - // Apply the log_add_exp operation directly - for (int i = 0; i < a.rows(); ++i) { - result(i) = stan::math::log_add_exp(a(i), b(i)); - } + // Check for infinity + if (a.array().isInf().any() || b.array().isInf().any()) { + result.setConstant(fvar(std::numeric_limits::quiet_NaN())); + return result; + } + + // Apply the log_add_exp operation directly + for (int i = 0; i < a.rows(); ++i) { + result(i) = stan::math::log_add_exp(a(i), b(i)); + } - return result; // Return the result vector + return result; // Return the result vector } // Specialization for nested fvar types diff --git a/stan/math/prim/fun/log_add_exp.hpp b/stan/math/prim/fun/log_add_exp.hpp index eeaf2d49612..d41aa4b485c 100644 --- a/stan/math/prim/fun/log_add_exp.hpp +++ b/stan/math/prim/fun/log_add_exp.hpp @@ -32,7 +32,6 @@ namespace math { template * = nullptr, require_all_stan_scalar_t* = nullptr> inline return_type_t log_add_exp(const T2& a, const T1& b) { - if (a == NEGATIVE_INFTY) { return b; } @@ -60,66 +59,69 @@ inline return_type_t log_add_exp(const T2& a, const T1& b) { */ template * = nullptr> inline auto log_add_exp(const T& a, const T& b) { - - // Check if sizes are compatible - if constexpr (stan::is_eigen::value) { - // Check if both matrices/vectors have the same dimensions - stan::math::check_matching_dims("log_add_exp", "a", a, "b", b); - - // Determine the number of rows and columns for the result - size_t rows = a.rows(); - size_t cols = b.cols(); - using return_t = return_type_t; - - Eigen::Matrix result(rows, cols); - - // Iterate over each element - for (size_t i = 0; i < rows; ++i) { - for (size_t j = 0; j < cols; ++j) { - double a_val = (a.cols() == 1) ? a(i, 0) : a(i, j); // Handle column vector or matrix - double b_val = (b.rows() == 1) ? b(0, j) : b(i, j); // Handle row vector or matrix - - if (a_val == NEGATIVE_INFTY) { - result(i, j) = b_val; - } else if (b_val == NEGATIVE_INFTY) { - result(i, j) = a_val; - } else if (a_val == INFTY || b_val == INFTY) { - result(i, j) = INFTY; - } else { - result(i, j) = log_sum_exp(a_val, b_val); - } - } - } - - return result; - } else if constexpr (std::is_same_v>) { - // Handle std::vector - if (a.size() != b.size()) { - throw std::invalid_argument("Sizes of x and y must match."); + // Check if sizes are compatible + if constexpr (stan::is_eigen::value) { + // Check if both matrices/vectors have the same dimensions + stan::math::check_matching_dims("log_add_exp", "a", a, "b", b); + + // Determine the number of rows and columns for the result + size_t rows = a.rows(); + size_t cols = b.cols(); + using return_t = return_type_t; + + Eigen::Matrix result(rows, cols); + + // Iterate over each element + for (size_t i = 0; i < rows; ++i) { + for (size_t j = 0; j < cols; ++j) { + double a_val = (a.cols() == 1) + ? a(i, 0) + : a(i, j); // Handle column vector or matrix + double b_val = (b.rows() == 1) + ? b(0, j) + : b(i, j); // Handle row vector or matrix + + if (a_val == NEGATIVE_INFTY) { + result(i, j) = b_val; + } else if (b_val == NEGATIVE_INFTY) { + result(i, j) = a_val; + } else if (a_val == INFTY || b_val == INFTY) { + result(i, j) = INFTY; + } else { + result(i, j) = log_sum_exp(a_val, b_val); } + } + } - using return_t = return_type_t; - std::vector result(a.size()); - - for (size_t i = 0; i < a.size(); ++i) { - double a_val = a[i]; - double b_val = b[i]; - - if (a_val == NEGATIVE_INFTY) { - result[i] = b_val; - } else if (b_val == NEGATIVE_INFTY) { - result[i] = a_val; - } else if (a_val == INFTY || b_val == INFTY) { - result[i] = INFTY; - } else { - result[i] = log_sum_exp(a_val, b_val); - } - } + return result; + } else if constexpr (std::is_same_v>) { + // Handle std::vector + if (a.size() != b.size()) { + throw std::invalid_argument("Sizes of x and y must match."); + } - return result; - } else { - throw std::invalid_argument("Unsupported container type."); + using return_t = return_type_t; + std::vector result(a.size()); + + for (size_t i = 0; i < a.size(); ++i) { + double a_val = a[i]; + double b_val = b[i]; + + if (a_val == NEGATIVE_INFTY) { + result[i] = b_val; + } else if (b_val == NEGATIVE_INFTY) { + result[i] = a_val; + } else if (a_val == INFTY || b_val == INFTY) { + result[i] = INFTY; + } else { + result[i] = log_sum_exp(a_val, b_val); + } } + + return result; + } else { + throw std::invalid_argument("Unsupported container type."); + } } /** @@ -134,23 +136,24 @@ inline auto log_add_exp(const T& a, const T& b) { */ template * = nullptr> inline auto log_add_exp(const T1& a, const T2& b) { - // Check if both are Eigen/vectors - if constexpr (stan::is_eigen::value && stan::is_eigen::value) { - // Check if both matrices/vectors have the same dimensions - stan::math::check_matching_dims("log_add_exp", "a", a, "b", b); - } else { - // Check if sizes are compatible for other types - if (a.size() != b.size()) { - throw std::invalid_argument("Sizes of x and y must match or be compatible."); - } + // Check if both are Eigen/vectors + if constexpr (stan::is_eigen::value && stan::is_eigen::value) { + // Check if both matrices/vectors have the same dimensions + stan::math::check_matching_dims("log_add_exp", "a", a, "b", b); + } else { + // Check if sizes are compatible for other types + if (a.size() != b.size()) { + throw std::invalid_argument( + "Sizes of x and y must match or be compatible."); } + } - // If dimensions are verified to match, apply the operation - return apply_scalar_binary( + // If dimensions are verified to match, apply the operation + return apply_scalar_binary( a, b, [](const auto& c, const auto& d) { return log_add_exp(c, d); }); } } // namespace math -} // nfamespace stan +} // namespace stan #endif diff --git a/test/unit/math/mix/fun/log_add_exp_test.cpp b/test/unit/math/mix/fun/log_add_exp_test.cpp index 024b2f331e0..cd58b217b83 100644 --- a/test/unit/math/mix/fun/log_add_exp_test.cpp +++ b/test/unit/math/mix/fun/log_add_exp_test.cpp @@ -13,7 +13,7 @@ TEST(MathMixMatFun, logAddExp) { y1 << 3.0, 2.0; stan::test::expect_ad(f, x1, y1); - //Test with negative infinity + // Test with negative infinity stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 1.0); stan::test::expect_ad(f, 1.0, stan::math::NEGATIVE_INFTY); @@ -50,7 +50,6 @@ TEST(MathMixMatFun, log_add_exp_elementwise_values) { EXPECT_TRUE(std::isinf(result[1])); } - TEST(MathMixMatFun, log_add_exp_mismatched_sizes) { auto f = [](const auto& x, const auto& y) { return stan::math::log_add_exp(x, y); @@ -73,14 +72,15 @@ TEST(MathMixMatFun, log_add_exp_container_tests) { x_row << 2.0, 1.0; Eigen::MatrixXd y_row(1, 2); y_row << 3.0, 2.0; - + stan::test::expect_ad(f, x_row, y_row); - + // Additional tests with mismatched sizes Eigen::MatrixXd x_mismatch(2, 1); x_mismatch << 0.5, -1.0; Eigen::MatrixXd y_mismatch(1, 3); y_mismatch << 1.0, 2.0, 3.0; - - EXPECT_THROW(stan::math::log_add_exp(x_mismatch, y_mismatch), std::invalid_argument); + + EXPECT_THROW(stan::math::log_add_exp(x_mismatch, y_mismatch), + std::invalid_argument); }