From 8330f931339e34b8646f318e1a066b7cadc652ec Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 30 Jul 2024 12:15:28 -0400 Subject: [PATCH 01/22] Add sum_to_zero constraint, free, and check --- stan/math/prim/constraint.hpp | 2 + .../prim/constraint/sum_to_zero_constrain.hpp | 106 ++++++++++++++++++ .../math/prim/constraint/sum_to_zero_free.hpp | 52 +++++++++ stan/math/prim/err.hpp | 1 + stan/math/prim/err/check_sum_to_zero.hpp | 67 +++++++++++ .../constraint/sum_to_zero_constrain_test.cpp | 58 ++++++++++ .../constraint/sum_to_zero_transform_test.cpp | 62 ++++++++++ .../math/prim/err/check_sum_to_zero_test.cpp | 93 +++++++++++++++ 8 files changed, 441 insertions(+) create mode 100644 stan/math/prim/constraint/sum_to_zero_constrain.hpp create mode 100644 stan/math/prim/constraint/sum_to_zero_free.hpp create mode 100644 stan/math/prim/err/check_sum_to_zero.hpp create mode 100644 test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp create mode 100644 test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp create mode 100644 test/unit/math/prim/err/check_sum_to_zero_test.cpp diff --git a/stan/math/prim/constraint.hpp b/stan/math/prim/constraint.hpp index 169fd8f2e31..5c494c83214 100644 --- a/stan/math/prim/constraint.hpp +++ b/stan/math/prim/constraint.hpp @@ -37,6 +37,8 @@ #include #include #include +#include +#include #include #include #include diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp new file mode 100644 index 00000000000..f7e612012fc --- /dev/null +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -0,0 +1,106 @@ +#ifndef STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP +#define STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP + +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return a vector with sum zero corresponding to the specified + * free vector. + * + * The sum-to-zero transform is defined such that the first K-1 + * elements are unconstrained and the last element is the negative + * sum of those elements. + * + * @tparam Vec type of the vector + * @param y Free vector input of dimensionality K - 1. + * @return Zero-sum vector of dimensionality K. + */ +template * = nullptr> +inline plain_type_t sum_to_zero_constrain(const Vec& y) { + using T = value_type_t; + + int Km1 = y.size(); + plain_type_t x(Km1 + 1); + // copy the first Km1 elements + x.head(Km1) = y; + // set the last element to -sum(y) + x.coeffRef(Km1) = -sum(y); + return x; +} + +/** + * Return a vector with sum zero corresponding to the specified + * free vector. + * + * The sum-to-zero transform is defined such that the first K-1 + * elements are unconstrained and the last element is the negative + * sum of those elements. This is a linear transform, with no + * Jacobian. + * + * @tparam Vec type of the vector + * @param y Free vector input of dimensionality K - 1. + * @param lp unused + * @return Zero-sum vector of dimensionality K. + */ +template * = nullptr> +inline plain_type_t sum_to_zero_constrain(const Vec& y, + value_type_t& lp) { + return sum_to_zero_constrain(y); +} + +/** + * Return a vector with sum zero corresponding to the specified + * free vector. + * + * The sum-to-zero transform is defined such that the first K-1 + * elements are unconstrained and the last element is the negative + * sum of those elements. This is a linear transform, with no + * Jacobian. + * + * @tparam Jacobian unused + * @tparam Vec A type inheriting from `Eigen::DenseBase` or a `var_value` with + * inner type inheriting from `Eigen::DenseBase` with compile time dynamic rows + * and 1 column + * @param[in] y free vector + * @param[in, out] lp unused + * @return Zero-sum vector of dimensionality one greater than `y` + */ +template * = nullptr> +inline plain_type_t sum_to_zero_constrain(const Vec& y, + return_type_t& lp) { + return sum_to_zero_constrain(y); +} + +/** + * Return a vector with sum zero corresponding to the specified + * free vector. + * + * The sum-to-zero transform is defined such that the first K-1 + * elements are unconstrained and the last element is the negative + * sum of those elements. This is a linear transform, with no + * Jacobian. + * + * @tparam Jacobian unused + * @tparam Vec A standard vector with inner type inheriting from + * `Eigen::DenseBase` or a `var_value` with inner type inheriting from + * `Eigen::DenseBase` with compile time dynamic rows and 1 column + * @param[in] y free vector + * @param[in, out] lp unused + * @return Zero-sum vectors of dimensionality one greater than `y` + */ +template * = nullptr> +inline auto sum_to_zero_constrain(const T& y, return_type_t& lp) { + return apply_vector_unary::apply( + y, [](auto&& v) { return sum_to_zero_constrain(v); }); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp new file mode 100644 index 00000000000..b34615257fe --- /dev/null +++ b/stan/math/prim/constraint/sum_to_zero_free.hpp @@ -0,0 +1,52 @@ +#ifndef STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_FREE_HPP +#define STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_FREE_HPP + +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return an unconstrained vector. + * + * The sum-to-zero transform is defined such that the first K-1 + * elements are unconstrained and the last element is the negative + * sum of those elements. + * + * @tparam ColVec a column vector type + * @param x Vector of length K. + * @return Free vector of length (K-1). + * @throw std::domain_error if x does not sum to zero + */ +template * = nullptr> +inline plain_type_t sum_to_zero_free(const Vec& x) { + const auto& x_ref = to_ref(x); + check_sum_to_zero("stan::math::sum_to_zero_free", "sum_to_zero variable", + x_ref); + if (x_ref.size() == 0) { + return plain_type_t(0); + } + return x_ref.head(x_ref.size() - 1); +} + +/** + * Overload of `sum_to_zero_free()` to untransform each Eigen vector + * in a standard vector. + * @tparam T A standard vector with with a `value_type` which inherits from + * `Eigen::MatrixBase` with compile time rows or columns equal to 1. + * @param x The standard vector to untransform. + */ +template * = nullptr> +auto sum_to_zero_free(const T& x) { + return apply_vector_unary::apply( + x, [](auto&& v) { return sum_to_zero_free(v); }); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/prim/err.hpp b/stan/math/prim/err.hpp index c78344c95d8..72d6c733e39 100644 --- a/stan/math/prim/err.hpp +++ b/stan/math/prim/err.hpp @@ -41,6 +41,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/err/check_sum_to_zero.hpp b/stan/math/prim/err/check_sum_to_zero.hpp new file mode 100644 index 00000000000..ef7b4fc2472 --- /dev/null +++ b/stan/math/prim/err/check_sum_to_zero.hpp @@ -0,0 +1,67 @@ +#ifndef STAN_MATH_PRIM_ERR_CHECK_SUM_TO_ZERO_HPP +#define STAN_MATH_PRIM_ERR_CHECK_SUM_TO_ZERO_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Throw an exception if the specified vector does not sum to 0. + * This function tests that the sum is within the tolerance specified by + * `CONSTRAINT_TOLERANCE`. + * This function only accepts Eigen vectors, statically + * typed vectors, not general matrices with 1 column. + * @tparam T A type inheriting from `Eigen::EigenBase` + * @param function Function name (for error messages) + * @param name Variable name (for error messages) + * @param theta Vector to test + * @throw `std::domain_error` if the vector does not sum to zero + */ +template * = nullptr> +void check_sum_to_zero(const char* function, const char* name, const T& theta) { + using std::fabs; + auto&& theta_ref = to_ref(value_of_rec(theta)); + if (!(fabs(theta_ref.sum()) <= CONSTRAINT_TOLERANCE)) { + [&]() STAN_COLD_PATH { + std::stringstream msg; + scalar_type_t sum = theta_ref.sum(); + msg << "does not sum to zero."; + msg.precision(10); + msg << " sum(" << name << ") = " << sum << ", but should be "; + std::string msg_str(msg.str()); + throw_domain_error(function, name, 0.0, msg_str.c_str()); + }(); + } +} + +/** + * Throw an exception if any vector in a standard vector does not sum to 0. + * This function tests that the sum is within the tolerance specified by + * `CONSTRAINT_TOLERANCE`. + * @tparam T A standard vector with inner type inheriting from + * `Eigen::EigenBase` + * @param function Function name (for error messages) + * @param name Variable name (for error messages) + * @param theta Vector to test. + * @throw `std::domain_error` if the vector does not sum to zero + */ +template * = nullptr> +void check_sum_to_zero(const char* function, const char* name, const T& theta) { + for (size_t i = 0; i < theta.size(); ++i) { + check_sum_to_zero(function, internal::make_iter_name(name, i).c_str(), + theta[i]); + } +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp b/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp new file mode 100644 index 00000000000..6423ec5cb77 --- /dev/null +++ b/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp @@ -0,0 +1,58 @@ +#include + +namespace sum_to_zero_constrain_test { +template +T g1(const T& x) { + stan::scalar_type_t lp = 0; + return stan::math::sum_to_zero_constrain(x, lp); +} +template +T g2(const T& x) { + stan::scalar_type_t lp = 0; + return stan::math::sum_to_zero_constrain(x, lp); +} +template +typename stan::scalar_type::type g3(const T& x) { + stan::scalar_type_t lp = 0; + stan::math::sum_to_zero_constrain(x, lp); + return lp; +} + +template +void expect_sum_to_zero_transform(const T& x) { + auto f1 = [](const auto& x) { return g1(x); }; + auto f2 = [](const auto& x) { return g2(x); }; + auto f3 = [](const auto& x) { return g3(x); }; + stan::test::expect_ad(f1, x); +// stan::test::expect_ad_matvar(f1, x); + stan::test::expect_ad(f2, x); +// stan::test::expect_ad_matvar(f2, x); + stan::test::expect_ad(f3, x); +// stan::test::expect_ad_matvar(f3, x); +} +} // namespace sum_to_zero_constrain_test + +TEST(MathMixMatFun, sum_to_zeroTransform) { + Eigen::VectorXd v0(0); + sum_to_zero_constrain_test::expect_sum_to_zero_transform(v0); + + Eigen::VectorXd v1(1); + v1 << 1; + sum_to_zero_constrain_test::expect_sum_to_zero_transform(v1); + + Eigen::VectorXd v2(2); + v2 << 3, -1; + sum_to_zero_constrain_test::expect_sum_to_zero_transform(v2); + + Eigen::VectorXd v3(3); + v3 << 2, 3, -1; + sum_to_zero_constrain_test::expect_sum_to_zero_transform(v3); + + Eigen::VectorXd v4(4); + v4 << 2, -1, 0, -1.1; + sum_to_zero_constrain_test::expect_sum_to_zero_transform(v4); + + Eigen::VectorXd v5(5); + v5 << 1, -3, 2, 0, -1; + sum_to_zero_constrain_test::expect_sum_to_zero_transform(v5); +} diff --git a/test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp b/test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp new file mode 100644 index 00000000000..d3114e0a663 --- /dev/null +++ b/test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp @@ -0,0 +1,62 @@ +#include +#include +#include + +TEST(prob_transform, sum_to_zero_rt0) { + using Eigen::Dynamic; + using Eigen::Matrix; + double lp = 0; + Matrix x(4); + x << 0.0, 0.0, 0.0, 0.0; + std::vector> x_vec{x, x, x}; + std::vector> y_vec + = stan::math::sum_to_zero_constrain(x_vec, lp); + for (auto&& y_i : y_vec) { + EXPECT_MATRIX_FLOAT_EQ(Eigen::VectorXd::Zero(5), y_i); + } + std::vector> xrt + = stan::math::sum_to_zero_free(y_vec); + EXPECT_EQ(x.size() + 1, y_vec[2].size()); + for (auto&& x_i : xrt) { + EXPECT_EQ(x.size(), x_i.size()); + for (int i = 0; i < x.size(); ++i) { + EXPECT_NEAR(x[i], x_i[i], 1E-10); + } + } +} +TEST(prob_transform, sum_to_zero_rt) { + using Eigen::Dynamic; + using Eigen::Matrix; + Matrix x(3); + x << 1.0, -1.0, 2.0; + Matrix y = stan::math::sum_to_zero_constrain(x); + EXPECT_NO_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y)); + Matrix xrt = stan::math::sum_to_zero_free(y); + EXPECT_EQ(x.size() + 1, y.size()); + EXPECT_EQ(x.size(), xrt.size()); + for (int i = 0; i < x.size(); ++i) { + EXPECT_FLOAT_EQ(x[i], xrt[i]); + } +} +TEST(prob_transform, sum_to_zero_match) { + using Eigen::Dynamic; + using Eigen::Matrix; + Matrix x(3); + x << 1.0, -1.0, 2.0; + double lp = 0; + Matrix y = stan::math::sum_to_zero_constrain(x); + Matrix y2 = stan::math::sum_to_zero_constrain(x, lp); + + EXPECT_EQ(4, y.size()); + EXPECT_EQ(4, y2.size()); + for (int i = 0; i < x.size(); ++i) + EXPECT_FLOAT_EQ(y[i], y2[i]); +} + +TEST(prob_transform, sum_to_zero_f_exception) { + using Eigen::Dynamic; + using Eigen::Matrix; + Matrix y(2); + y << 0.5, -0.55; + EXPECT_THROW(stan::math::sum_to_zero_free(y), std::domain_error); +} diff --git a/test/unit/math/prim/err/check_sum_to_zero_test.cpp b/test/unit/math/prim/err/check_sum_to_zero_test.cpp new file mode 100644 index 00000000000..e637f8294ca --- /dev/null +++ b/test/unit/math/prim/err/check_sum_to_zero_test.cpp @@ -0,0 +1,93 @@ +#include +#include +#include +#include +#include + +TEST(ErrorHandlingMatrix, checkSumToZero_edges) { + Eigen::Matrix zero(0); + + EXPECT_NO_THROW( + stan::math::check_sum_to_zero("checkSumToZero", "zero", zero)); + + Eigen::Matrix y_vec(1); + y_vec << 0.0; + + EXPECT_NO_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y_vec)); + + y_vec[0] = 0.1; + EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y_vec), + std::domain_error); +} + +TEST(ErrorHandlingMatrix, checkSumToZero) { + Eigen::Matrix y_vec(2); + std::vector> y{y_vec, y_vec, y_vec}; + for (auto& y_i : y) { + y_i << 0.5, -0.5; + } + + EXPECT_NO_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y)); + + for (auto& y_i : y) { + y_i[0] = 0.55; + } + EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y), + std::domain_error); +} + +TEST(ErrorHandlingMatrix, checkSumToZero_message_sum) { + std::string message; + Eigen::Matrix y_vec(100); + y_vec.setZero(); + std::vector> y{y_vec, y_vec, y_vec}; + for (auto& y_i : y) { + y_i[13] = 0.1; + } + + try { + stan::math::check_sum_to_zero("checkSumToZero", "y", y); + FAIL() << "should have thrown"; + } catch (std::domain_error& e) { + message = e.what(); + } catch (...) { + FAIL() << "threw the wrong error"; + } + + EXPECT_TRUE(std::string::npos != message.find(" y[1] does not sum to zero")) + << message; + + EXPECT_TRUE(std::string::npos != message.find("sum(y[1]) = 0.1")) << message; +} + + +TEST(ErrorHandlingMatrix, checkSumToZero_nan) { + Eigen::Matrix y_vec(2); + std::vector> y{y_vec, y_vec, y_vec}; + constexpr double nan = std::numeric_limits::quiet_NaN(); + for (auto& y_i : y) { + y_i << nan, 0.5; + } + + EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y), + std::domain_error); + + for (auto& y_i : y) { + y_i[1] = 0.55; + } + EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y), + std::domain_error); + + for (auto& y_i : y) { + y_i[0] = 0.5; + y_i[1] = nan; + } + EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y), + std::domain_error); + + for (auto& y_i : y) { + y_i[0] = nan; + } + EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y), + std::domain_error); +} From 45649359a910d4e0ff26fb3eb423f578e31d12e2 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 30 Jul 2024 12:20:07 -0400 Subject: [PATCH 02/22] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- .../unit/math/mix/constraint/sum_to_zero_constrain_test.cpp | 6 +++--- test/unit/math/prim/err/check_sum_to_zero_test.cpp | 1 - 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp b/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp index 6423ec5cb77..cec8c2cd068 100644 --- a/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp +++ b/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp @@ -24,11 +24,11 @@ void expect_sum_to_zero_transform(const T& x) { auto f2 = [](const auto& x) { return g2(x); }; auto f3 = [](const auto& x) { return g3(x); }; stan::test::expect_ad(f1, x); -// stan::test::expect_ad_matvar(f1, x); + // stan::test::expect_ad_matvar(f1, x); stan::test::expect_ad(f2, x); -// stan::test::expect_ad_matvar(f2, x); + // stan::test::expect_ad_matvar(f2, x); stan::test::expect_ad(f3, x); -// stan::test::expect_ad_matvar(f3, x); + // stan::test::expect_ad_matvar(f3, x); } } // namespace sum_to_zero_constrain_test diff --git a/test/unit/math/prim/err/check_sum_to_zero_test.cpp b/test/unit/math/prim/err/check_sum_to_zero_test.cpp index e637f8294ca..713ecfe33a6 100644 --- a/test/unit/math/prim/err/check_sum_to_zero_test.cpp +++ b/test/unit/math/prim/err/check_sum_to_zero_test.cpp @@ -60,7 +60,6 @@ TEST(ErrorHandlingMatrix, checkSumToZero_message_sum) { EXPECT_TRUE(std::string::npos != message.find("sum(y[1]) = 0.1")) << message; } - TEST(ErrorHandlingMatrix, checkSumToZero_nan) { Eigen::Matrix y_vec(2); std::vector> y{y_vec, y_vec, y_vec}; From c75a7939a036525e4f4deca66d93bdb330964978 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 30 Jul 2024 14:16:28 -0400 Subject: [PATCH 03/22] Header check fix --- stan/math/prim/constraint/sum_to_zero_constrain.hpp | 3 +-- stan/math/prim/constraint/sum_to_zero_free.hpp | 1 + 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index f7e612012fc..277ae99544e 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include namespace stan { @@ -23,8 +24,6 @@ namespace math { */ template * = nullptr> inline plain_type_t sum_to_zero_constrain(const Vec& y) { - using T = value_type_t; - int Km1 = y.size(); plain_type_t x(Km1 + 1); // copy the first Km1 elements diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp index b34615257fe..1dec7fea53c 100644 --- a/stan/math/prim/constraint/sum_to_zero_free.hpp +++ b/stan/math/prim/constraint/sum_to_zero_free.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include namespace stan { From 69ffa68bc872c5d36369f34b79272bcf5bdfd5fb Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 30 Jul 2024 15:30:31 -0400 Subject: [PATCH 04/22] Try rev overload --- .../prim/constraint/sum_to_zero_constrain.hpp | 6 +- stan/math/rev/constraint.hpp | 1 + .../rev/constraint/sum_to_zero_constrain.hpp | 79 +++++++++++++++++++ 3 files changed, 84 insertions(+), 2 deletions(-) create mode 100644 stan/math/rev/constraint/sum_to_zero_constrain.hpp diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index 277ae99544e..92fa4a2daa0 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -22,7 +22,8 @@ namespace math { * @param y Free vector input of dimensionality K - 1. * @return Zero-sum vector of dimensionality K. */ -template * = nullptr> +template * = nullptr, + require_not_st_var* = nullptr> inline plain_type_t sum_to_zero_constrain(const Vec& y) { int Km1 = y.size(); plain_type_t x(Km1 + 1); @@ -47,7 +48,8 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y) { * @param lp unused * @return Zero-sum vector of dimensionality K. */ -template * = nullptr> +template * = nullptr, + require_not_st_var* = nullptr> inline plain_type_t sum_to_zero_constrain(const Vec& y, value_type_t& lp) { return sum_to_zero_constrain(y); diff --git a/stan/math/rev/constraint.hpp b/stan/math/rev/constraint.hpp index 3ab4a6f72c7..ac036a31f52 100644 --- a/stan/math/rev/constraint.hpp +++ b/stan/math/rev/constraint.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp new file mode 100644 index 00000000000..05aa581c412 --- /dev/null +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -0,0 +1,79 @@ +#ifndef STAN_MATH_REV_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP +#define STAN_MATH_REV_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return a vector with sum zero corresponding to the specified + * free vector. + * + * The sum-to-zero transform is defined such that the first K-1 + * elements are unconstrained and the last element is the negative + * sum of those elements. + * + * @tparam T type of the vector + * @param y Free vector input of dimensionality K - 1. + * @return Zero-sum vector of dimensionality K. + */ +template * = nullptr> +inline auto sum_to_zero_constrain(const T& y) { + using ret_type = plain_type_t; + + size_t N = y.size(); + Eigen::VectorXd x_val(N + 1); + + arena_t arena_y = y; + + if (unlikely(N == 0)) { + x_val << 0; + return ret_type(x_val); + } + + x_val.head(N) = y.val(); + arena_t arena_x = x_val; + + var x_N = -sum(y); + + arena_x.coeffRef(N) = x_N; + + + reverse_pass_callback([arena_y, arena_x, x_N, N]() mutable { + arena_y.adj() += arena_x.adj().head(N); + x_N.adj() += arena_x.adj().coeff(N); + }); + + return ret_type(arena_x); +} + +/** + * Return a vector with sum zero corresponding to the specified + * free vector. + * + * The sum-to-zero transform is defined such that the first K-1 + * elements are unconstrained and the last element is the negative + * sum of those elements. This is a linear transform, with no + * Jacobian. + * + * @tparam Vec type of the vector + * @param y Free vector input of dimensionality K - 1. + * @param lp unused + * @return Zero-sum vector of dimensionality K. + */ +template * = nullptr> +auto sum_to_zero_constrain(const T& y, scalar_type_t& lp) { + return sum_to_zero_constrain(y); +} + +} // namespace math +} // namespace stan +#endif From fb97a2433441284398b40a2f9634933d6e678e67 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 30 Jul 2024 15:31:33 -0400 Subject: [PATCH 05/22] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/rev/constraint/sum_to_zero_constrain.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index 05aa581c412..10cbcd295de 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -46,10 +46,9 @@ inline auto sum_to_zero_constrain(const T& y) { arena_x.coeffRef(N) = x_N; - reverse_pass_callback([arena_y, arena_x, x_N, N]() mutable { arena_y.adj() += arena_x.adj().head(N); - x_N.adj() += arena_x.adj().coeff(N); + x_N.adj() += arena_x.adj().coeff(N); }); return ret_type(arena_x); From 3fb6ff49948fa92efcae782bbd904b7922809d72 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Tue, 30 Jul 2024 16:33:19 -0400 Subject: [PATCH 06/22] fix bug in reverse pass for sum to zero constraint --- .../prim/constraint/sum_to_zero_constrain.hpp | 7 +++-- .../rev/constraint/sum_to_zero_constrain.hpp | 28 +++++++------------ .../constraint/sum_to_zero_constrain_test.cpp | 6 ++-- 3 files changed, 17 insertions(+), 24 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index 92fa4a2daa0..212aed2790a 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -25,12 +25,13 @@ namespace math { template * = nullptr, require_not_st_var* = nullptr> inline plain_type_t sum_to_zero_constrain(const Vec& y) { - int Km1 = y.size(); + const auto Km1 = y.size(); plain_type_t x(Km1 + 1); // copy the first Km1 elements - x.head(Km1) = y; + auto&& y_ref = to_ref(y); + x.head(Km1) = y_ref; // set the last element to -sum(y) - x.coeffRef(Km1) = -sum(y); + x.coeffRef(Km1) = -sum(y_ref); return x; } diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index 10cbcd295de..01878004691 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -29,28 +29,20 @@ template * = nullptr> inline auto sum_to_zero_constrain(const T& y) { using ret_type = plain_type_t; - size_t N = y.size(); - Eigen::VectorXd x_val(N + 1); - - arena_t arena_y = y; - + const auto N = y.size(); if (unlikely(N == 0)) { - x_val << 0; - return ret_type(x_val); + return ret_type(Eigen::VectorXd{{0}}); } - - x_val.head(N) = y.val(); + Eigen::VectorXd x_val = Eigen::VectorXd::Zero(N + 1); + auto arena_y = to_arena(y); + double x_sum = -sum(arena_y.val()); + x_val.head(N) = arena_y.val(); + x_val(N) = x_sum; arena_t arena_x = x_val; - - var x_N = -sum(y); - - arena_x.coeffRef(N) = x_N; - - reverse_pass_callback([arena_y, arena_x, x_N, N]() mutable { - arena_y.adj() += arena_x.adj().head(N); - x_N.adj() += arena_x.adj().coeff(N); + reverse_pass_callback([arena_y, arena_x, x_sum, N]() mutable { + arena_y.adj().array() -= arena_x.adj_op()(N); + arena_y.adj() += arena_x.adj_op().head(N); }); - return ret_type(arena_x); } diff --git a/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp b/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp index cec8c2cd068..ae4ce5516db 100644 --- a/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp +++ b/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp @@ -24,11 +24,11 @@ void expect_sum_to_zero_transform(const T& x) { auto f2 = [](const auto& x) { return g2(x); }; auto f3 = [](const auto& x) { return g3(x); }; stan::test::expect_ad(f1, x); - // stan::test::expect_ad_matvar(f1, x); + stan::test::expect_ad_matvar(f1, x); stan::test::expect_ad(f2, x); - // stan::test::expect_ad_matvar(f2, x); + stan::test::expect_ad_matvar(f2, x); stan::test::expect_ad(f3, x); - // stan::test::expect_ad_matvar(f3, x); + stan::test::expect_ad_matvar(f3, x); } } // namespace sum_to_zero_constrain_test From 40d0e4b7bf95d51f9159e1a0a8ffbaf4d1be32b2 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 30 Jul 2024 16:34:20 -0400 Subject: [PATCH 07/22] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- .../unit/math/mix/constraint/sum_to_zero_constrain_test.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp b/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp index ae4ce5516db..a42db039fcd 100644 --- a/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp +++ b/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp @@ -24,11 +24,11 @@ void expect_sum_to_zero_transform(const T& x) { auto f2 = [](const auto& x) { return g2(x); }; auto f3 = [](const auto& x) { return g3(x); }; stan::test::expect_ad(f1, x); - stan::test::expect_ad_matvar(f1, x); + stan::test::expect_ad_matvar(f1, x); stan::test::expect_ad(f2, x); - stan::test::expect_ad_matvar(f2, x); + stan::test::expect_ad_matvar(f2, x); stan::test::expect_ad(f3, x); - stan::test::expect_ad_matvar(f3, x); + stan::test::expect_ad_matvar(f3, x); } } // namespace sum_to_zero_constrain_test From cf9b012cd76f132a758ecea4e3aecce16d3bf96e Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 31 Jul 2024 12:03:36 -0400 Subject: [PATCH 08/22] Changes per Steve's comments --- stan/math/prim/constraint/sum_to_zero_constrain.hpp | 3 +++ stan/math/prim/constraint/sum_to_zero_free.hpp | 4 +--- stan/math/prim/err/check_sum_to_zero.hpp | 7 ++++++- stan/math/rev/constraint/sum_to_zero_constrain.hpp | 6 +++--- test/unit/math/prim/err/check_sum_to_zero_test.cpp | 4 ++-- 5 files changed, 15 insertions(+), 9 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index 212aed2790a..528e979061c 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -26,6 +26,9 @@ template * = nullptr, require_not_st_var* = nullptr> inline plain_type_t sum_to_zero_constrain(const Vec& y) { const auto Km1 = y.size(); + if (unlikely(Km1 == 0)) { + return plain_type_t(Eigen::VectorXd{{0}}); + } plain_type_t x(Km1 + 1); // copy the first Km1 elements auto&& y_ref = to_ref(y); diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp index 1dec7fea53c..ff1a7390711 100644 --- a/stan/math/prim/constraint/sum_to_zero_free.hpp +++ b/stan/math/prim/constraint/sum_to_zero_free.hpp @@ -28,9 +28,7 @@ inline plain_type_t sum_to_zero_free(const Vec& x) { const auto& x_ref = to_ref(x); check_sum_to_zero("stan::math::sum_to_zero_free", "sum_to_zero variable", x_ref); - if (x_ref.size() == 0) { - return plain_type_t(0); - } + return x_ref.head(x_ref.size() - 1); } diff --git a/stan/math/prim/err/check_sum_to_zero.hpp b/stan/math/prim/err/check_sum_to_zero.hpp index ef7b4fc2472..b384f5cc778 100644 --- a/stan/math/prim/err/check_sum_to_zero.hpp +++ b/stan/math/prim/err/check_sum_to_zero.hpp @@ -24,13 +24,17 @@ namespace math { * @param function Function name (for error messages) * @param name Variable name (for error messages) * @param theta Vector to test + * @throw `std::invalid_argument` if `theta` is a 0-vector * @throw `std::domain_error` if the vector does not sum to zero */ template * = nullptr> void check_sum_to_zero(const char* function, const char* name, const T& theta) { using std::fabs; + // the size-zero case is technically a valid sum-to-zero vector, + // but it cannot be unconstrained to anything + check_nonzero_size(function, name, theta); auto&& theta_ref = to_ref(value_of_rec(theta)); - if (!(fabs(theta_ref.sum()) <= CONSTRAINT_TOLERANCE)) { + if (unlikely(!(fabs(theta_ref.sum()) <= CONSTRAINT_TOLERANCE))) { [&]() STAN_COLD_PATH { std::stringstream msg; scalar_type_t sum = theta_ref.sum(); @@ -52,6 +56,7 @@ void check_sum_to_zero(const char* function, const char* name, const T& theta) { * @param function Function name (for error messages) * @param name Variable name (for error messages) * @param theta Vector to test. + * @throw `std::invalid_argument` if `theta` is a 0-vector * @throw `std::domain_error` if the vector does not sum to zero */ template * = nullptr> diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index 01878004691..fe664c660c4 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -31,7 +31,7 @@ inline auto sum_to_zero_constrain(const T& y) { const auto N = y.size(); if (unlikely(N == 0)) { - return ret_type(Eigen::VectorXd{{0}}); + return arena_t(Eigen::VectorXd{{0}}); } Eigen::VectorXd x_val = Eigen::VectorXd::Zero(N + 1); auto arena_y = to_arena(y); @@ -43,7 +43,7 @@ inline auto sum_to_zero_constrain(const T& y) { arena_y.adj().array() -= arena_x.adj_op()(N); arena_y.adj() += arena_x.adj_op().head(N); }); - return ret_type(arena_x); + return arena_x; } /** @@ -61,7 +61,7 @@ inline auto sum_to_zero_constrain(const T& y) { * @return Zero-sum vector of dimensionality K. */ template * = nullptr> -auto sum_to_zero_constrain(const T& y, scalar_type_t& lp) { +inline auto sum_to_zero_constrain(const T& y, scalar_type_t& lp) { return sum_to_zero_constrain(y); } diff --git a/test/unit/math/prim/err/check_sum_to_zero_test.cpp b/test/unit/math/prim/err/check_sum_to_zero_test.cpp index 713ecfe33a6..1e6ad716a81 100644 --- a/test/unit/math/prim/err/check_sum_to_zero_test.cpp +++ b/test/unit/math/prim/err/check_sum_to_zero_test.cpp @@ -7,8 +7,8 @@ TEST(ErrorHandlingMatrix, checkSumToZero_edges) { Eigen::Matrix zero(0); - EXPECT_NO_THROW( - stan::math::check_sum_to_zero("checkSumToZero", "zero", zero)); + EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "zero", zero), + std::invalid_argument); Eigen::Matrix y_vec(1); y_vec << 0.0; From ed4a64d3c6404724f94391b3a4e7b7aa53453f25 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 1 Aug 2024 12:49:24 -0400 Subject: [PATCH 09/22] inv ilr constrain prim impl --- .../prim/constraint/sum_to_zero_constrain.hpp | 23 +++++++++++-------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index 528e979061c..a383ae60d20 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -14,9 +14,8 @@ namespace math { * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined such that the first K-1 - * elements are unconstrained and the last element is the negative - * sum of those elements. + * The sum-to-zero transform is defined using the inverse of the + * isometric log ratio transform (ILR) * * @tparam Vec type of the vector * @param y Free vector input of dimensionality K - 1. @@ -25,16 +24,20 @@ namespace math { template * = nullptr, require_not_st_var* = nullptr> inline plain_type_t sum_to_zero_constrain(const Vec& y) { - const auto Km1 = y.size(); - if (unlikely(Km1 == 0)) { + const auto N = y.size() + 1; + if (unlikely(N == 1)) { return plain_type_t(Eigen::VectorXd{{0}}); } - plain_type_t x(Km1 + 1); - // copy the first Km1 elements + auto&& y_ref = to_ref(y); - x.head(Km1) = y_ref; - // set the last element to -sum(y) - x.coeffRef(Km1) = -sum(y_ref); + + // straightforward port of Stan code, may be optimized + Eigen::VectorXd ns = linspaced_vector(N - 1, 1, N - 1); + auto w = y_ref.array() / sqrt(ns.array() * (ns.array() + 1)); + + plain_type_t x = append_row(reverse(cumulative_sum(reverse(w))), 0) + - append_row(0, ns.array() * w.array()); + return x; } From aa41da41e11da8300546d71c1523c5ff56b03212 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 1 Aug 2024 15:10:53 -0400 Subject: [PATCH 10/22] Single loop --- .../prim/constraint/sum_to_zero_constrain.hpp | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index a383ae60d20..27f2259cde8 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -25,18 +25,24 @@ template * = nullptr, require_not_st_var* = nullptr> inline plain_type_t sum_to_zero_constrain(const Vec& y) { const auto N = y.size() + 1; + + plain_type_t x = Eigen::VectorXd::Zero(N); if (unlikely(N == 1)) { - return plain_type_t(Eigen::VectorXd{{0}}); + return x; } auto&& y_ref = to_ref(y); - // straightforward port of Stan code, may be optimized - Eigen::VectorXd ns = linspaced_vector(N - 1, 1, N - 1); - auto w = y_ref.array() / sqrt(ns.array() * (ns.array() + 1)); + typename plain_type_t::Scalar cumsum(0); + for (int i = N - 2; i >= 0; --i) { + double n = i + 1; + auto w = y_ref(i) * inv_sqrt(n * (n + 1)); + cumsum += w; + + x.coeffRef(N - 2 - i) += cumsum; - plain_type_t x = append_row(reverse(cumulative_sum(reverse(w))), 0) - - append_row(0, ns.array() * w.array()); + x.coeffRef(i + 1) -= w * n; + } return x; } From ba179cb83f2520955271f74b3587fed27bea72f8 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 2 Aug 2024 10:42:32 -0400 Subject: [PATCH 11/22] off-by-1 fixes --- .../prim/constraint/sum_to_zero_constrain.hpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index 27f2259cde8..c8c2e7a4c17 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -24,23 +24,22 @@ namespace math { template * = nullptr, require_not_st_var* = nullptr> inline plain_type_t sum_to_zero_constrain(const Vec& y) { - const auto N = y.size() + 1; + const auto N = y.size(); - plain_type_t x = Eigen::VectorXd::Zero(N); - if (unlikely(N == 1)) { + plain_type_t x = Eigen::VectorXd::Zero(N+1); + if (unlikely(N == 0)) { return x; } auto&& y_ref = to_ref(y); - typename plain_type_t::Scalar cumsum(0); - for (int i = N - 2; i >= 0; --i) { + typename plain_type_t::Scalar total(0); + for (int i = N - 1; i >= 0; --i) { double n = i + 1; auto w = y_ref(i) * inv_sqrt(n * (n + 1)); - cumsum += w; - - x.coeffRef(N - 2 - i) += cumsum; + total += w; + x.coeffRef(i) += total; x.coeffRef(i + 1) -= w * n; } From 4197bfd809f6021212f9b785cc7cb5e9ada8d196 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 2 Aug 2024 11:10:29 -0400 Subject: [PATCH 12/22] Names --- stan/math/prim/constraint/sum_to_zero_constrain.hpp | 10 +++++----- stan/math/rev/constraint/sum_to_zero_constrain.hpp | 13 ++++--------- 2 files changed, 9 insertions(+), 14 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index c8c2e7a4c17..0aa84f26b16 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -26,9 +26,9 @@ template * = nullptr, inline plain_type_t sum_to_zero_constrain(const Vec& y) { const auto N = y.size(); - plain_type_t x = Eigen::VectorXd::Zero(N+1); + plain_type_t z = Eigen::VectorXd::Zero(N+1); if (unlikely(N == 0)) { - return x; + return z; } auto&& y_ref = to_ref(y); @@ -39,11 +39,11 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y) { auto w = y_ref(i) * inv_sqrt(n * (n + 1)); total += w; - x.coeffRef(i) += total; - x.coeffRef(i + 1) -= w * n; + z.coeffRef(i) += total; + z.coeffRef(i + 1) -= w * n; } - return x; + return z; } /** diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index fe664c660c4..3e8433d71ee 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -28,20 +28,15 @@ namespace math { template * = nullptr> inline auto sum_to_zero_constrain(const T& y) { using ret_type = plain_type_t; - const auto N = y.size(); if (unlikely(N == 0)) { return arena_t(Eigen::VectorXd{{0}}); } - Eigen::VectorXd x_val = Eigen::VectorXd::Zero(N + 1); auto arena_y = to_arena(y); - double x_sum = -sum(arena_y.val()); - x_val.head(N) = arena_y.val(); - x_val(N) = x_sum; - arena_t arena_x = x_val; - reverse_pass_callback([arena_y, arena_x, x_sum, N]() mutable { - arena_y.adj().array() -= arena_x.adj_op()(N); - arena_y.adj() += arena_x.adj_op().head(N); + arena_t arena_x = sum_to_zero_constrain(arena_y.val()); + + reverse_pass_callback([arena_y, arena_x]() mutable { + // TODO }); return arena_x; } From 004c7fe8a0815213405b0d1725f1a52d0d36f9d9 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 2 Aug 2024 11:43:32 -0400 Subject: [PATCH 13/22] Free impl --- .../math/prim/constraint/sum_to_zero_free.hpp | 35 ++++++++++++++----- .../rev/constraint/sum_to_zero_constrain.hpp | 1 + .../constraint/sum_to_zero_transform_test.cpp | 2 ++ 3 files changed, 29 insertions(+), 9 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp index ff1a7390711..9a8106cf092 100644 --- a/stan/math/prim/constraint/sum_to_zero_free.hpp +++ b/stan/math/prim/constraint/sum_to_zero_free.hpp @@ -19,17 +19,34 @@ namespace math { * sum of those elements. * * @tparam ColVec a column vector type - * @param x Vector of length K. + * @param z Vector of length K. * @return Free vector of length (K-1). - * @throw std::domain_error if x does not sum to zero + * @throw std::domain_error if z does not sum to zero */ template * = nullptr> -inline plain_type_t sum_to_zero_free(const Vec& x) { - const auto& x_ref = to_ref(x); +inline plain_type_t sum_to_zero_free(const Vec& z) { + const auto& z_ref = to_ref(z); check_sum_to_zero("stan::math::sum_to_zero_free", "sum_to_zero variable", - x_ref); + z_ref); - return x_ref.head(x_ref.size() - 1); + const auto N = z.size() - 1; + + plain_type_t y = Eigen::VectorXd::Zero(N); + if (unlikely(N == 0)) { + return y; + } + + y.coeffRef(N - 1) = -z_ref(N) * sqrt(N * (N + 1)) / N; + typename plain_type_t::Scalar total(0); + + for (int i = N - 2; i >= 0; --i) { + double n = i + 1; + auto w = y(i + 1) / sqrt((n + 1) * (n + 2)); + total += w; + y.coeffRef(i) = (total - z_ref(i + 1)) * sqrt(n * (n + 1)) / n; + } + + return y; } /** @@ -37,12 +54,12 @@ inline plain_type_t sum_to_zero_free(const Vec& x) { * in a standard vector. * @tparam T A standard vector with with a `value_type` which inherits from * `Eigen::MatrixBase` with compile time rows or columns equal to 1. - * @param x The standard vector to untransform. + * @param z The standard vector to untransform. */ template * = nullptr> -auto sum_to_zero_free(const T& x) { +auto sum_to_zero_free(const T& z) { return apply_vector_unary::apply( - x, [](auto&& v) { return sum_to_zero_free(v); }); + z, [](auto&& v) { return sum_to_zero_free(v); }); } } // namespace math diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index 3e8433d71ee..b546f7a49cb 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include diff --git a/test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp b/test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp index d3114e0a663..3c09aa1b768 100644 --- a/test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp +++ b/test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp @@ -11,6 +11,8 @@ TEST(prob_transform, sum_to_zero_rt0) { std::vector> x_vec{x, x, x}; std::vector> y_vec = stan::math::sum_to_zero_constrain(x_vec, lp); + EXPECT_NO_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y_vec)); + for (auto&& y_i : y_vec) { EXPECT_MATRIX_FLOAT_EQ(Eigen::VectorXd::Zero(5), y_i); } From e132643a4c20d86983b1f48a5da072cf357deac0 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Sat, 3 Aug 2024 17:39:08 -0400 Subject: [PATCH 14/22] Basic rev impl --- .../prim/constraint/sum_to_zero_constrain.hpp | 14 +++---- .../math/prim/constraint/sum_to_zero_free.hpp | 7 ++-- .../rev/constraint/sum_to_zero_constrain.hpp | 37 ++++++++++++++++--- 3 files changed, 42 insertions(+), 16 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index 0aa84f26b16..03e0d5b8ab6 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -33,14 +33,14 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y) { auto&& y_ref = to_ref(y); - typename plain_type_t::Scalar total(0); - for (int i = N - 1; i >= 0; --i) { - double n = i + 1; - auto w = y_ref(i) * inv_sqrt(n * (n + 1)); - total += w; + typename plain_type_t::Scalar sum_w(0); + for (int i = N; i > 0; --i) { + double n = i; + auto w = y_ref(i-1) * inv_sqrt(n * (n + 1)); + sum_w += w; - z.coeffRef(i) += total; - z.coeffRef(i + 1) -= w * n; + z.coeffRef(i-1) += sum_w; + z.coeffRef(i) -= w * n; } return z; diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp index 9a8106cf092..5c601a90297 100644 --- a/stan/math/prim/constraint/sum_to_zero_free.hpp +++ b/stan/math/prim/constraint/sum_to_zero_free.hpp @@ -37,13 +37,14 @@ inline plain_type_t sum_to_zero_free(const Vec& z) { } y.coeffRef(N - 1) = -z_ref(N) * sqrt(N * (N + 1)) / N; - typename plain_type_t::Scalar total(0); + + typename plain_type_t::Scalar sum_w(0); for (int i = N - 2; i >= 0; --i) { double n = i + 1; auto w = y(i + 1) / sqrt((n + 1) * (n + 2)); - total += w; - y.coeffRef(i) = (total - z_ref(i + 1)) * sqrt(n * (n + 1)) / n; + sum_w += w; + y.coeffRef(i) = (sum_w - z_ref(i + 1)) * sqrt(n * (n + 1)) / n; } return y; diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index b546f7a49cb..c6cba13c51b 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -7,6 +7,8 @@ #include #include #include +#include +#include #include #include #include @@ -29,17 +31,40 @@ namespace math { template * = nullptr> inline auto sum_to_zero_constrain(const T& y) { using ret_type = plain_type_t; - const auto N = y.size(); - if (unlikely(N == 0)) { + if (unlikely(y.size() == 0)) { return arena_t(Eigen::VectorXd{{0}}); } auto arena_y = to_arena(y); - arena_t arena_x = sum_to_zero_constrain(arena_y.val()); + arena_t arena_z = sum_to_zero_constrain(arena_y.val()); - reverse_pass_callback([arena_y, arena_x]() mutable { - // TODO + reverse_pass_callback([arena_y, arena_z]() mutable { + const auto N = arena_y.size(); + + Eigen::VectorXd ns = linspaced_vector(N, 1, N ); + + /* + u3.adj += z.adj[1:N-1] + v.adj -= z.adj[2:N] + w.adj += v.adj .* ns + u2.adj += reverse(u3.adj) + u1.adj += cumulative_sum(u2.adj) + w.adj += reverse(u1.adj) + y.adj += w.adj ./ sqrt(ns .* (ns + 1)) + */ + + Eigen::VectorXd u_adj = arena_z.adj_op().head(N).eval(); + Eigen::VectorXd v_adj = -arena_z.adj_op().tail(N).eval(); + + Eigen::VectorXd w_from_v = v_adj.array() * ns.array(); + + Eigen::VectorXd w_from_u = cumulative_sum(u_adj); + + Eigen::VectorXd w_adj = (w_from_v.array() + w_from_u.array()); + + arena_y.adj() += (w_adj.array() / sqrt(ns.array() * (ns.array() + 1))).matrix(); }); - return arena_x; + + return arena_z; } /** From 9cf06f12c91b69647ab14b6055685f4167f16ac7 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Sat, 3 Aug 2024 17:48:19 -0400 Subject: [PATCH 15/22] Loop version --- .../rev/constraint/sum_to_zero_constrain.hpp | 27 +++++++------------ 1 file changed, 9 insertions(+), 18 deletions(-) diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index c6cba13c51b..e27d99039aa 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -40,28 +40,19 @@ inline auto sum_to_zero_constrain(const T& y) { reverse_pass_callback([arena_y, arena_z]() mutable { const auto N = arena_y.size(); - Eigen::VectorXd ns = linspaced_vector(N, 1, N ); + double sum_u_adj = 0; + for (int i = 0; i < N; ++i) { + double n = i + 1; - /* - u3.adj += z.adj[1:N-1] - v.adj -= z.adj[2:N] - w.adj += v.adj .* ns - u2.adj += reverse(u3.adj) - u1.adj += cumulative_sum(u2.adj) - w.adj += reverse(u1.adj) - y.adj += w.adj ./ sqrt(ns .* (ns + 1)) - */ + double u_adj = arena_z.adj()(i); + sum_u_adj += u_adj; - Eigen::VectorXd u_adj = arena_z.adj_op().head(N).eval(); - Eigen::VectorXd v_adj = -arena_z.adj_op().tail(N).eval(); + double v_adj = -arena_z.adj()(i + 1); - Eigen::VectorXd w_from_v = v_adj.array() * ns.array(); + double w = (v_adj * n) + sum_u_adj; - Eigen::VectorXd w_from_u = cumulative_sum(u_adj); - - Eigen::VectorXd w_adj = (w_from_v.array() + w_from_u.array()); - - arena_y.adj() += (w_adj.array() / sqrt(ns.array() * (ns.array() + 1))).matrix(); + arena_y.adj()(i) += w / sqrt(n * (n + 1)); + } }); return arena_z; From 86fed93a6dcf52a183cb217f42766e9070e9327b Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Sat, 3 Aug 2024 17:59:37 -0400 Subject: [PATCH 16/22] Clarifying --- stan/math/rev/constraint/sum_to_zero_constrain.hpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index e27d99039aa..b722d6b4f72 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -44,14 +44,15 @@ inline auto sum_to_zero_constrain(const T& y) { for (int i = 0; i < N; ++i) { double n = i + 1; - double u_adj = arena_z.adj()(i); - sum_u_adj += u_adj; + // adjoint of the reverse cumulative sum computed in the forward mode + sum_u_adj += arena_z.adj()(i); - double v_adj = -arena_z.adj()(i + 1); + // adjoint of the offset subtraction + double v_adj = -arena_z.adj()(i + 1) * n; - double w = (v_adj * n) + sum_u_adj; + double w_adj = v_adj + sum_u_adj; - arena_y.adj()(i) += w / sqrt(n * (n + 1)); + arena_y.adj()(i) += w_adj / sqrt(n * (n + 1)); } }); From 25460e922b95aac0976c94c1f59946782d218399 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Sat, 3 Aug 2024 18:02:25 -0400 Subject: [PATCH 17/22] Docstrings --- .../prim/constraint/sum_to_zero_constrain.hpp | 56 +++++++++++++------ .../rev/constraint/sum_to_zero_constrain.hpp | 25 ++++++--- 2 files changed, 57 insertions(+), 24 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index 03e0d5b8ab6..c7baa81f306 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -14,8 +14,15 @@ namespace math { * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined using the inverse of the - * isometric log ratio transform (ILR) + * The sum-to-zero transform is defined using a modified version of the + * the inverse of the isometric log ratio transform (ILR). + * See: + * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; + * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for + * compositional data analysis", Mathematical Geology, 35 (3): 279–300, + * doi:10.1023/A:1023818214614, S2CID 122844634 + * + * This is a linear transform, with no Jacobian. * * @tparam Vec type of the vector * @param y Free vector input of dimensionality K - 1. @@ -26,7 +33,7 @@ template * = nullptr, inline plain_type_t sum_to_zero_constrain(const Vec& y) { const auto N = y.size(); - plain_type_t z = Eigen::VectorXd::Zero(N+1); + plain_type_t z = Eigen::VectorXd::Zero(N + 1); if (unlikely(N == 0)) { return z; } @@ -36,10 +43,10 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y) { typename plain_type_t::Scalar sum_w(0); for (int i = N; i > 0; --i) { double n = i; - auto w = y_ref(i-1) * inv_sqrt(n * (n + 1)); + auto w = y_ref(i - 1) * inv_sqrt(n * (n + 1)); sum_w += w; - z.coeffRef(i-1) += sum_w; + z.coeffRef(i - 1) += sum_w; z.coeffRef(i) -= w * n; } @@ -50,10 +57,15 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y) { * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined such that the first K-1 - * elements are unconstrained and the last element is the negative - * sum of those elements. This is a linear transform, with no - * Jacobian. + * The sum-to-zero transform is defined using a modified version of the + * the inverse of the isometric log ratio transform (ILR). + * See: + * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; + * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for + * compositional data analysis", Mathematical Geology, 35 (3): 279–300, + * doi:10.1023/A:1023818214614, S2CID 122844634 + * + * This is a linear transform, with no Jacobian. * * @tparam Vec type of the vector * @param y Free vector input of dimensionality K - 1. @@ -71,10 +83,15 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y, * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined such that the first K-1 - * elements are unconstrained and the last element is the negative - * sum of those elements. This is a linear transform, with no - * Jacobian. + * The sum-to-zero transform is defined using a modified version of the + * the inverse of the isometric log ratio transform (ILR). + * See: + * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; + * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for + * compositional data analysis", Mathematical Geology, 35 (3): 279–300, + * doi:10.1023/A:1023818214614, S2CID 122844634 + * + * This is a linear transform, with no Jacobian. * * @tparam Jacobian unused * @tparam Vec A type inheriting from `Eigen::DenseBase` or a `var_value` with @@ -94,10 +111,15 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y, * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined such that the first K-1 - * elements are unconstrained and the last element is the negative - * sum of those elements. This is a linear transform, with no - * Jacobian. + * The sum-to-zero transform is defined using a modified version of the + * the inverse of the isometric log ratio transform (ILR). + * See: + * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; + * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for + * compositional data analysis", Mathematical Geology, 35 (3): 279–300, + * doi:10.1023/A:1023818214614, S2CID 122844634 + * + * This is a linear transform, with no Jacobian. * * @tparam Jacobian unused * @tparam Vec A standard vector with inner type inheriting from diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index b722d6b4f72..f23c0b650a1 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -20,9 +20,15 @@ namespace math { * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined such that the first K-1 - * elements are unconstrained and the last element is the negative - * sum of those elements. + * The sum-to-zero transform is defined using a modified version of the + * the inverse of the isometric log ratio transform (ILR). + * See: + * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; + * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for + * compositional data analysis", Mathematical Geology, 35 (3): 279–300, + * doi:10.1023/A:1023818214614, S2CID 122844634 + * + * This is a linear transform, with no Jacobian. * * @tparam T type of the vector * @param y Free vector input of dimensionality K - 1. @@ -63,10 +69,15 @@ inline auto sum_to_zero_constrain(const T& y) { * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined such that the first K-1 - * elements are unconstrained and the last element is the negative - * sum of those elements. This is a linear transform, with no - * Jacobian. + * The sum-to-zero transform is defined using a modified version of the + * the inverse of the isometric log ratio transform (ILR). + * See: + * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; + * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for + * compositional data analysis", Mathematical Geology, 35 (3): 279–300, + * doi:10.1023/A:1023818214614, S2CID 122844634 + * + * This is a linear transform, with no Jacobian. * * @tparam Vec type of the vector * @param y Free vector input of dimensionality K - 1. From edd826efdabaed32e7883013da89135faf191d03 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Mon, 5 Aug 2024 09:29:18 -0400 Subject: [PATCH 18/22] Comment fixes --- stan/math/prim/constraint/sum_to_zero_constrain.hpp | 6 +++--- stan/math/prim/constraint/sum_to_zero_free.hpp | 10 +++++++--- stan/math/rev/constraint/sum_to_zero_constrain.hpp | 6 +++--- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index c7baa81f306..02c7e52697c 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -83,7 +83,7 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y, * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined using a modified version of the + * The sum-to-zero transform is defined using a modified version of * the inverse of the isometric log ratio transform (ILR). * See: * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; @@ -111,7 +111,7 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y, * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined using a modified version of the + * The sum-to-zero transform is defined using a modified version of * the inverse of the isometric log ratio transform (ILR). * See: * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; @@ -119,7 +119,7 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y, * compositional data analysis", Mathematical Geology, 35 (3): 279–300, * doi:10.1023/A:1023818214614, S2CID 122844634 * - * This is a linear transform, with no Jacobian. + * This is a linear transform, with no Jacobian. * * @tparam Jacobian unused * @tparam Vec A standard vector with inner type inheriting from diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp index 5c601a90297..672fb3a3919 100644 --- a/stan/math/prim/constraint/sum_to_zero_free.hpp +++ b/stan/math/prim/constraint/sum_to_zero_free.hpp @@ -14,9 +14,13 @@ namespace math { /** * Return an unconstrained vector. * - * The sum-to-zero transform is defined such that the first K-1 - * elements are unconstrained and the last element is the negative - * sum of those elements. + * The sum-to-zero transform is defined using a modified version of the + * isometric log ratio transform (ILR). + * See: + * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; + * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for + * compositional data analysis", Mathematical Geology, 35 (3): 279–300, + * doi:10.1023/A:1023818214614, S2CID 122844634 * * @tparam ColVec a column vector type * @param z Vector of length K. diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index f23c0b650a1..91ed0d3504f 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -20,7 +20,7 @@ namespace math { * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined using a modified version of the + * The sum-to-zero transform is defined using a modified version of * the inverse of the isometric log ratio transform (ILR). * See: * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; @@ -69,7 +69,7 @@ inline auto sum_to_zero_constrain(const T& y) { * Return a vector with sum zero corresponding to the specified * free vector. * - * The sum-to-zero transform is defined using a modified version of the + * The sum-to-zero transform is defined using a modified version of * the inverse of the isometric log ratio transform (ILR). * See: * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria; @@ -77,7 +77,7 @@ inline auto sum_to_zero_constrain(const T& y) { * compositional data analysis", Mathematical Geology, 35 (3): 279–300, * doi:10.1023/A:1023818214614, S2CID 122844634 * - * This is a linear transform, with no Jacobian. + * This is a linear transform, with no Jacobian. * * @tparam Vec type of the vector * @param y Free vector input of dimensionality K - 1. From 8be67c6511cbbbf454b7253faf2fe145399dfdf6 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Mon, 5 Aug 2024 10:12:01 -0400 Subject: [PATCH 19/22] Header fixes --- stan/math/prim/constraint/sum_to_zero_constrain.hpp | 3 ++- stan/math/prim/constraint/sum_to_zero_free.hpp | 1 + stan/math/rev/constraint/sum_to_zero_constrain.hpp | 4 +--- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index 02c7e52697c..e67e482e2ae 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -3,7 +3,8 @@ #include #include -#include +#include +#include #include #include diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp index 672fb3a3919..28a63162511 100644 --- a/stan/math/prim/constraint/sum_to_zero_free.hpp +++ b/stan/math/prim/constraint/sum_to_zero_free.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index 91ed0d3504f..b6743d153ca 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -4,11 +4,9 @@ #include #include #include -#include #include +#include #include -#include -#include #include #include #include From c64bace1ec515aac5af7703760089647a1bb8673 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Mon, 5 Aug 2024 14:25:17 -0400 Subject: [PATCH 20/22] Second citation from @spinkney --- .../prim/constraint/sum_to_zero_constrain.hpp | 24 +++++++++++++++++++ .../math/prim/constraint/sum_to_zero_free.hpp | 6 +++++ .../rev/constraint/sum_to_zero_constrain.hpp | 12 ++++++++++ 3 files changed, 42 insertions(+) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index e67e482e2ae..31bc3db4121 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -23,6 +23,12 @@ namespace math { * compositional data analysis", Mathematical Geology, 35 (3): 279–300, * doi:10.1023/A:1023818214614, S2CID 122844634 * + * This implementation is closer to the description of the same using "pivot + * coordinates" in + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. + * In: Applied Compositional Data Analysis. Springer Series in Statistics. + * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * * This is a linear transform, with no Jacobian. * * @tparam Vec type of the vector @@ -66,6 +72,12 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y) { * compositional data analysis", Mathematical Geology, 35 (3): 279–300, * doi:10.1023/A:1023818214614, S2CID 122844634 * + * This implementation is closer to the description of the same using "pivot + * coordinates" in + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. + * In: Applied Compositional Data Analysis. Springer Series in Statistics. + * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * * This is a linear transform, with no Jacobian. * * @tparam Vec type of the vector @@ -92,6 +104,12 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y, * compositional data analysis", Mathematical Geology, 35 (3): 279–300, * doi:10.1023/A:1023818214614, S2CID 122844634 * + * This implementation is closer to the description of the same using "pivot + * coordinates" in + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. + * In: Applied Compositional Data Analysis. Springer Series in Statistics. + * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * * This is a linear transform, with no Jacobian. * * @tparam Jacobian unused @@ -120,6 +138,12 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y, * compositional data analysis", Mathematical Geology, 35 (3): 279–300, * doi:10.1023/A:1023818214614, S2CID 122844634 * + * This implementation is closer to the description of the same using "pivot + * coordinates" in + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. + * In: Applied Compositional Data Analysis. Springer Series in Statistics. + * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * * This is a linear transform, with no Jacobian. * * @tparam Jacobian unused diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp index 28a63162511..11a5962243b 100644 --- a/stan/math/prim/constraint/sum_to_zero_free.hpp +++ b/stan/math/prim/constraint/sum_to_zero_free.hpp @@ -23,6 +23,12 @@ namespace math { * compositional data analysis", Mathematical Geology, 35 (3): 279–300, * doi:10.1023/A:1023818214614, S2CID 122844634 * + * This implementation is closer to the description of the same using "pivot + * coordinates" in + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. + * In: Applied Compositional Data Analysis. Springer Series in Statistics. + * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * * @tparam ColVec a column vector type * @param z Vector of length K. * @return Free vector of length (K-1). diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index b6743d153ca..7a404f3f325 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -26,6 +26,12 @@ namespace math { * compositional data analysis", Mathematical Geology, 35 (3): 279–300, * doi:10.1023/A:1023818214614, S2CID 122844634 * + * This implementation is closer to the description of the same using "pivot + * coordinates" in + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. + * In: Applied Compositional Data Analysis. Springer Series in Statistics. + * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * * This is a linear transform, with no Jacobian. * * @tparam T type of the vector @@ -75,6 +81,12 @@ inline auto sum_to_zero_constrain(const T& y) { * compositional data analysis", Mathematical Geology, 35 (3): 279–300, * doi:10.1023/A:1023818214614, S2CID 122844634 * + * This implementation is closer to the description of the same using "pivot + * coordinates" in + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. + * In: Applied Compositional Data Analysis. Springer Series in Statistics. + * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * * This is a linear transform, with no Jacobian. * * @tparam Vec type of the vector From 51ac4c61cc3b42f703a539dbcd996abadcd1e38e Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 5 Aug 2024 14:26:22 -0400 Subject: [PATCH 21/22] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- .../prim/constraint/sum_to_zero_constrain.hpp | 24 +++++++++---------- .../math/prim/constraint/sum_to_zero_free.hpp | 6 ++--- .../rev/constraint/sum_to_zero_constrain.hpp | 12 +++++----- 3 files changed, 21 insertions(+), 21 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index 31bc3db4121..d56a1e9ce52 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -25,9 +25,9 @@ namespace math { * * This implementation is closer to the description of the same using "pivot * coordinates" in - * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. - * In: Applied Compositional Data Analysis. Springer Series in Statistics. - * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of + * Compositional Data. In: Applied Compositional Data Analysis. Springer Series + * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 * * This is a linear transform, with no Jacobian. * @@ -74,9 +74,9 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y) { * * This implementation is closer to the description of the same using "pivot * coordinates" in - * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. - * In: Applied Compositional Data Analysis. Springer Series in Statistics. - * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of + * Compositional Data. In: Applied Compositional Data Analysis. Springer Series + * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 * * This is a linear transform, with no Jacobian. * @@ -106,9 +106,9 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y, * * This implementation is closer to the description of the same using "pivot * coordinates" in - * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. - * In: Applied Compositional Data Analysis. Springer Series in Statistics. - * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of + * Compositional Data. In: Applied Compositional Data Analysis. Springer Series + * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 * * This is a linear transform, with no Jacobian. * @@ -140,9 +140,9 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y, * * This implementation is closer to the description of the same using "pivot * coordinates" in - * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. - * In: Applied Compositional Data Analysis. Springer Series in Statistics. - * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of + * Compositional Data. In: Applied Compositional Data Analysis. Springer Series + * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 * * This is a linear transform, with no Jacobian. * diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp index 11a5962243b..19164e7ff3b 100644 --- a/stan/math/prim/constraint/sum_to_zero_free.hpp +++ b/stan/math/prim/constraint/sum_to_zero_free.hpp @@ -25,9 +25,9 @@ namespace math { * * This implementation is closer to the description of the same using "pivot * coordinates" in - * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. - * In: Applied Compositional Data Analysis. Springer Series in Statistics. - * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of + * Compositional Data. In: Applied Compositional Data Analysis. Springer Series + * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 * * @tparam ColVec a column vector type * @param z Vector of length K. diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index 7a404f3f325..dd56944acb8 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -28,9 +28,9 @@ namespace math { * * This implementation is closer to the description of the same using "pivot * coordinates" in - * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. - * In: Applied Compositional Data Analysis. Springer Series in Statistics. - * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of + * Compositional Data. In: Applied Compositional Data Analysis. Springer Series + * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 * * This is a linear transform, with no Jacobian. * @@ -83,9 +83,9 @@ inline auto sum_to_zero_constrain(const T& y) { * * This implementation is closer to the description of the same using "pivot * coordinates" in - * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of Compositional Data. - * In: Applied Compositional Data Analysis. Springer Series in Statistics. - * Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 + * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of + * Compositional Data. In: Applied Compositional Data Analysis. Springer Series + * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3 * * This is a linear transform, with no Jacobian. * From 5ae902d4c9461d9f1e1043606338883dcc090752 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 6 Aug 2024 16:50:57 -0400 Subject: [PATCH 22/22] Changes per @SteveBronder --- stan/math/prim/constraint/sum_to_zero_constrain.hpp | 4 ++-- stan/math/prim/constraint/sum_to_zero_free.hpp | 4 ++-- stan/math/rev/constraint/sum_to_zero_constrain.hpp | 10 +++++----- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp index d56a1e9ce52..c3979720c11 100644 --- a/stan/math/prim/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp @@ -47,9 +47,9 @@ inline plain_type_t sum_to_zero_constrain(const Vec& y) { auto&& y_ref = to_ref(y); - typename plain_type_t::Scalar sum_w(0); + value_type_t sum_w(0); for (int i = N; i > 0; --i) { - double n = i; + double n = static_cast(i); auto w = y_ref(i - 1) * inv_sqrt(n * (n + 1)); sum_w += w; diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp index 19164e7ff3b..650970f992d 100644 --- a/stan/math/prim/constraint/sum_to_zero_free.hpp +++ b/stan/math/prim/constraint/sum_to_zero_free.hpp @@ -49,10 +49,10 @@ inline plain_type_t sum_to_zero_free(const Vec& z) { y.coeffRef(N - 1) = -z_ref(N) * sqrt(N * (N + 1)) / N; - typename plain_type_t::Scalar sum_w(0); + value_type_t sum_w(0); for (int i = N - 2; i >= 0; --i) { - double n = i + 1; + double n = static_cast(i + 1); auto w = y(i + 1) / sqrt((n + 1) * (n + 2)); sum_w += w; y.coeffRef(i) = (sum_w - z_ref(i + 1)) * sqrt(n * (n + 1)) / n; diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp index dd56944acb8..38248fe1413 100644 --- a/stan/math/rev/constraint/sum_to_zero_constrain.hpp +++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp @@ -39,12 +39,12 @@ namespace math { * @return Zero-sum vector of dimensionality K. */ template * = nullptr> -inline auto sum_to_zero_constrain(const T& y) { +inline auto sum_to_zero_constrain(T&& y) { using ret_type = plain_type_t; if (unlikely(y.size() == 0)) { return arena_t(Eigen::VectorXd{{0}}); } - auto arena_y = to_arena(y); + auto arena_y = to_arena(std::forward(y)); arena_t arena_z = sum_to_zero_constrain(arena_y.val()); reverse_pass_callback([arena_y, arena_z]() mutable { @@ -52,7 +52,7 @@ inline auto sum_to_zero_constrain(const T& y) { double sum_u_adj = 0; for (int i = 0; i < N; ++i) { - double n = i + 1; + double n = static_cast(i + 1); // adjoint of the reverse cumulative sum computed in the forward mode sum_u_adj += arena_z.adj()(i); @@ -95,8 +95,8 @@ inline auto sum_to_zero_constrain(const T& y) { * @return Zero-sum vector of dimensionality K. */ template * = nullptr> -inline auto sum_to_zero_constrain(const T& y, scalar_type_t& lp) { - return sum_to_zero_constrain(y); +inline auto sum_to_zero_constrain(T&& y, scalar_type_t& lp) { + return sum_to_zero_constrain(std::forward(y)); } } // namespace math