From a73eabe594bde42d15c922f05012ca30b8dd98b1 Mon Sep 17 00:00:00 2001 From: Ben Date: Fri, 11 Dec 2020 09:50:16 -0500 Subject: [PATCH 01/30] Initial commit to make matrix/vector constrains varmat friendly (Issue #2101) --- stan/math/prim/fun/exp.hpp | 7 +- stan/math/prim/fun/log1m.hpp | 3 +- stan/math/prim/fun/read_corr_L.hpp | 6 +- stan/math/prim/fun/square.hpp | 3 +- stan/math/prim/fun/tanh.hpp | 3 +- stan/math/prim/functor/apply_scalar_unary.hpp | 2 +- stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/exp.hpp | 19 ++ stan/math/rev/fun/log1m.hpp | 19 ++ stan/math/rev/fun/read_corr_L.hpp | 173 ++++++++++++++++++ stan/math/rev/fun/square.hpp | 19 ++ stan/math/rev/fun/tanh.hpp | 22 ++- .../unit/math/mix/fun/corr_constrain_test.cpp | 75 ++++++++ test/unit/math/mix/fun/exp_test.cpp | 18 ++ test/unit/math/mix/fun/log1m_test.cpp | 18 ++ test/unit/math/mix/fun/read_corr_L_test.cpp | 53 ++++++ test/unit/math/mix/fun/square_test.cpp | 18 ++ test/unit/math/mix/fun/tanh_test.cpp | 18 ++ 18 files changed, 466 insertions(+), 11 deletions(-) create mode 100644 stan/math/rev/fun/read_corr_L.hpp create mode 100644 test/unit/math/mix/fun/corr_constrain_test.cpp create mode 100644 test/unit/math/mix/fun/read_corr_L_test.cpp diff --git a/stan/math/prim/fun/exp.hpp b/stan/math/prim/fun/exp.hpp index 6766280f6af..29111a16069 100644 --- a/stan/math/prim/fun/exp.hpp +++ b/stan/math/prim/fun/exp.hpp @@ -40,9 +40,10 @@ struct exp_fun { * @return Elementwise application of exponentiation to the argument. */ template < - typename Container, - require_not_container_st* = nullptr, - require_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> + typename Container, + require_not_container_st* = nullptr, + require_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr, + require_not_var_matrix_t* = nullptr> inline auto exp(const Container& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/log1m.hpp b/stan/math/prim/fun/log1m.hpp index a24eafef36e..3c2892601ad 100644 --- a/stan/math/prim/fun/log1m.hpp +++ b/stan/math/prim/fun/log1m.hpp @@ -67,7 +67,8 @@ struct log1m_fun { * @param x container * @return Natural log of 1 minus each value in x. */ -template +template * = nullptr> inline auto log1m(const T& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/read_corr_L.hpp b/stan/math/prim/fun/read_corr_L.hpp index 659bc872b0a..51fa8238fb8 100644 --- a/stan/math/prim/fun/read_corr_L.hpp +++ b/stan/math/prim/fun/read_corr_L.hpp @@ -34,9 +34,9 @@ namespace math { * canonical partial correlations. */ template * = nullptr> -Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> read_corr_L( - const T& CPCs, // on (-1, 1) - size_t K) { +Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> read_corr_L(const T& CPCs, // on (-1, 1) + size_t K) { + using T_scalar = value_type_t; if (K == 0) { return {}; diff --git a/stan/math/prim/fun/square.hpp b/stan/math/prim/fun/square.hpp index 96e0426ebd2..92162e0c790 100644 --- a/stan/math/prim/fun/square.hpp +++ b/stan/math/prim/fun/square.hpp @@ -48,7 +48,8 @@ struct square_fun { * @return Each value in x squared. */ template * = nullptr> + require_not_container_st* = nullptr, + require_not_var_matrix_t* = nullptr> inline auto square(const Container& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/tanh.hpp b/stan/math/prim/fun/tanh.hpp index 2b93b6a437d..ec2ab88b684 100644 --- a/stan/math/prim/fun/tanh.hpp +++ b/stan/math/prim/fun/tanh.hpp @@ -36,7 +36,8 @@ struct tanh_fun { * @return Hyperbolic tangent of each value in x. */ template * = nullptr> + require_not_container_st* = nullptr, + require_not_var_matrix_t* = nullptr> inline auto tanh(const Container& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/functor/apply_scalar_unary.hpp b/stan/math/prim/functor/apply_scalar_unary.hpp index 6a249dd92cd..d21470d6e9a 100644 --- a/stan/math/prim/functor/apply_scalar_unary.hpp +++ b/stan/math/prim/functor/apply_scalar_unary.hpp @@ -51,7 +51,7 @@ struct apply_scalar_unary> { /** * Type of underlying scalar for the matrix type T. */ - using scalar_t = typename Eigen::internal::traits::Scalar; + using scalar_t = value_type_t; /** * Return the result of applying the function defined by the diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index 4b9de3df0ba..739dc064a24 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -126,6 +126,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/exp.hpp b/stan/math/rev/fun/exp.hpp index 4a14a4811b1..d6611158140 100644 --- a/stan/math/rev/fun/exp.hpp +++ b/stan/math/rev/fun/exp.hpp @@ -57,6 +57,25 @@ inline std::complex exp(const std::complex& z) { return internal::complex_exp(z); } +/** + * Return the exponentiation of the elements of x + * + * @tparam T type of x + * @param x argument + * @return elementwise exponentiation of x + */ +template * = nullptr> +inline auto exp(const T& x) { + T res = x.val().array().exp().matrix(); + + reverse_pass_callback([x, res]() mutable { + x.adj().noalias() += (res.val().array() * res.adj().array()).matrix(); + }); + + return res; +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/fun/log1m.hpp b/stan/math/rev/fun/log1m.hpp index d8c4d3abbd0..2fdc4ec4213 100644 --- a/stan/math/rev/fun/log1m.hpp +++ b/stan/math/rev/fun/log1m.hpp @@ -28,6 +28,25 @@ class log1m_vari : public op_v_vari { */ inline var log1m(const var& a) { return var(new internal::log1m_vari(a.vi_)); } +/** + * Return the elementwise log of 1 - x + * + * @tparam T type of x + * @param x argument + * @return elementwise log of 1 - x + */ +template * = nullptr> +inline auto log1m(const T& x) { + T res = stan::math::log1m(x.val()); + + reverse_pass_callback([x, res]() mutable { + x.adj().noalias() += (res.adj().array() / (x.val().array() - 1.0)).matrix(); + }); + + return res; +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/fun/read_corr_L.hpp b/stan/math/rev/fun/read_corr_L.hpp new file mode 100644 index 00000000000..3064f4022a4 --- /dev/null +++ b/stan/math/rev/fun/read_corr_L.hpp @@ -0,0 +1,173 @@ +#ifndef STAN_MATH_REV_FUN_READ_CORR_L_HPP +#define STAN_MATH_REV_FUN_READ_CORR_L_HPP + +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return the Cholesky factor of the correlation matrix of the + * specified dimensionality corresponding to the specified + * canonical partial correlations. + * + * It is generally better to work with the Cholesky factor rather + * than the correlation matrix itself when the determinant, + * inverse, etc. of the correlation matrix is needed for some + * statistical calculation. + * + *

See read_corr_matrix(Array, size_t, T) + * for more information. + * + * @tparam T type of input (must be a `var_value` where `T` + * inherits from EigenBase) + * @param CPCs The (K choose 2) canonical partial correlations in + * (-1, 1). + * @param K Dimensionality of correlation matrix. + * @return Cholesky factor of correlation matrix for specified + * canonical partial correlations. + */ +template * = nullptr> +auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) + using ret_type = var_value; + + if (K == 0) { + return ret_type(Eigen::MatrixXd()); + } + if (K == 1) { + return ret_type(Eigen::MatrixXd::Identity(1, 1)); + } + + using std::sqrt; + Eigen::ArrayXd acc(K - 1); + acc.setOnes(); + // Cholesky factor of correlation matrix + Eigen::MatrixXd L = Eigen::MatrixXd::Zero(K, K); + + size_t position = 0; + size_t pull = K - 1; + + L(0, 0) = 1.0; + L.col(0).tail(pull) = CPCs.val().head(pull); + acc.tail(pull) = 1.0 - CPCs.val().head(pull).array().square(); + for (size_t i = 1; i < (K - 1); i++) { + position += pull; + pull = K - 1 - i; + + const auto& cpc_seg = CPCs.val().array().segment(position, pull); + L(i, i) = sqrt(acc(i - 1)); + L.col(i).tail(pull) = cpc_seg * acc.tail(pull).sqrt(); + acc.tail(pull) = acc.tail(pull) * (1.0 - cpc_seg.square()); + } + + L(K - 1, K - 1) = sqrt(acc(K - 2)); + + ret_type L_res = L; + auto arena_acc = to_arena(acc); + + reverse_pass_callback([CPCs, arena_acc, K, L_res]() mutable { + Eigen::ArrayXd acc_val = arena_acc; + Eigen::ArrayXd acc_adj = Eigen::ArrayXd::Zero(K - 1); + + acc_adj(K - 2) += 0.5 * L_res.adj()(K - 1, K - 1) / L_res.val()(K - 1, K - 1); + + int position = CPCs.size() - 1; + for (size_t i = K - 2; i > 0; --i) { + int pull = K - 1 - i; + + const auto& cpc_seg_val = CPCs.val().array().segment(position, pull); + auto cpc_seg_adj = CPCs.adj().array().segment(position, pull); + acc_val.tail(pull) /= 1.0 - cpc_seg_val.square(); + cpc_seg_adj -= 2.0 * acc_adj.tail(pull) * acc_val.tail(pull) * cpc_seg_val; + acc_adj.tail(pull) = (acc_adj.tail(pull).array() * (1.0 - cpc_seg_val.square())).eval(); + cpc_seg_adj += L_res.adj().array().col(i).tail(pull) * acc_val.tail(pull).sqrt(); + acc_adj.tail(pull) += 0.5 * L_res.adj().array().col(i).tail(pull) * cpc_seg_val / acc_val.tail(pull).sqrt(); + acc_adj(i - 1) += 0.5 * L_res.adj()(i, i) / L_res.val()(i, i); + + position -= pull + 1; + } + + int pull = K - 1; + CPCs.adj().array().head(pull) -= 2.0 * acc_adj.tail(pull) * CPCs.val().array().head(pull); + CPCs.adj().head(pull) += L_res.adj().col(0).tail(pull); + }); + + return L_res; +} + +/** + * Return the Cholesky factor of the correlation matrix of the + * specified dimensionality corresponding to the specified + * canonical partial correlations, incrementing the specified + * scalar reference with the log absolute determinant of the + * Jacobian of the transformation. + * + *

The implementation is Ben Goodrich's Cholesky + * factor-based approach to the C-vine method of: + * + *

  • Daniel Lewandowski, Dorota Kurowicka, and Harry Joe, + * Generating random correlation matrices based on vines and + * extended onion method Journal of Multivariate Analysis 100 + * (2009) 1989–2001
+ * + * @tparam T type of input (must be a `var_value` where `T` + * inherits from EigenBase) + * @param CPCs The (K choose 2) canonical partial correlations in + * (-1, 1). + * @param K Dimensionality of correlation matrix. + * @param log_prob Reference to variable to increment with the log + * Jacobian determinant. + * @return Cholesky factor of correlation matrix for specified + * partial correlations. + */ +template * = nullptr, + require_stan_scalar_t* = nullptr> +auto read_corr_L(const T1& CPCs, size_t K, T2& log_prob) { + using ret_val_type = Eigen::MatrixXd; + using ret_type = var_value; + + if (K == 0) { + return ret_type(ret_val_type()); + } + if (K == 1) { + return ret_type(Eigen::MatrixXd::Identity(1, 1)); + } + + size_t pos = 0; + double acc_val = 0; + // no need to abs() because this Jacobian determinant + // is strictly positive (and triangular) + // see inverse of Jacobian in equation 11 of LKJ paper + for (size_t k = 1; k <= (K - 2); k++) { + for (size_t i = k + 1; i <= K; i++) { + acc_val += (K - k - 1) * log1m(square(CPCs.val()(pos))); + pos++; + } + } + + log_prob += 0.5 * acc_val; + + reverse_pass_callback([K, CPCs, log_prob]() mutable { + double acc_adj = log_prob.adj() * 0.5; + + size_t pos = CPCs.size() - 2; + for (size_t k = K - 2; k >= 1; --k) { + for(size_t i = K; i >= k + 1; --i) { + CPCs.adj()(pos) -= 2.0 * acc_adj * (K - k - 1) * CPCs.val()(pos) / (1.0 - square(CPCs.val()(pos))); + pos--; + } + } + }); + + return read_corr_L(CPCs, K); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/rev/fun/square.hpp b/stan/math/rev/fun/square.hpp index 7de37294524..83a0ba6f80d 100644 --- a/stan/math/rev/fun/square.hpp +++ b/stan/math/rev/fun/square.hpp @@ -44,6 +44,25 @@ inline var square(const var& x) { return var(new internal::square_vari(x.vi_)); } +/** + * Return the elementwise square of x + * + * @tparam T type of x + * @param x argument + * @return elementwise square of x + */ +template * = nullptr> +inline auto square(const T& x) { + T res = (x.val().array() * x.val().array()).matrix(); + + reverse_pass_callback([x, res]() mutable { + x.adj().noalias() += (2.0 * x.val().array() * res.adj().array()).matrix(); + }); + + return res; +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/fun/tanh.hpp b/stan/math/rev/fun/tanh.hpp index ef9763989e2..7769b9f49a8 100644 --- a/stan/math/rev/fun/tanh.hpp +++ b/stan/math/rev/fun/tanh.hpp @@ -2,9 +2,9 @@ #define STAN_MATH_REV_FUN_TANH_HPP #include +#include #include #include -#include #include #include @@ -61,6 +61,26 @@ inline std::complex tanh(const std::complex& z) { return stan::math::internal::complex_tanh(z); } +/** + * Return the hyperbolic tangent of of x + * + * @tparam T type of x + * @param x argument + * @return elementwise hyperbolic tangent of x + */ +template * = nullptr> +inline auto tanh(const T& x) { + T res = stan::math::tanh(x.val()); + + reverse_pass_callback([x, res]() mutable { + auto cosh = stan::math::cosh(x.val()); + x.adj().noalias() += (res.adj().array() / (cosh.array() * cosh.array())).matrix(); + }); + + return res; +} + } // namespace math } // namespace stan #endif diff --git a/test/unit/math/mix/fun/corr_constrain_test.cpp b/test/unit/math/mix/fun/corr_constrain_test.cpp new file mode 100644 index 00000000000..57ac6689a18 --- /dev/null +++ b/test/unit/math/mix/fun/corr_constrain_test.cpp @@ -0,0 +1,75 @@ +#include + +TEST(mathMixMatFun, corr_constrain) { + auto f = [](const auto& x1) { + return stan::math::corr_constrain(x1); + }; + + std::vector x0 = { -1.0, 2.0, 3.0 }; + Eigen::VectorXd x1(3); + x1 << -1.0, -2.0, 0.3; + Eigen::RowVectorXd x2(3); + x2 << -1.0, -2.0, 0.3; + stan::test::expect_ad(f, x0); + stan::test::expect_ad(f, x1); + stan::test::expect_ad(f, x2); + stan::test::expect_ad_matvar(f, x1); + stan::test::expect_ad_matvar(f, x2); + + std::vector x4; + Eigen::VectorXd x5(0); + Eigen::RowVectorXd x6(0); + + stan::test::expect_ad(f, x4); + stan::test::expect_ad(f, x5); + stan::test::expect_ad(f, x6); + + stan::test::expect_ad_matvar(f, x5); + stan::test::expect_ad_matvar(f, x6); +} + +TEST(mathMixMatFun, corr_constrain_lp) { + auto f1 = [](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + return stan::math::corr_constrain(x1, lp); + }; + + auto f2 = [](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + stan::math::corr_constrain(x1, lp); + return lp; + }; + + std::vector x0 = { -1.0, 2.0, 3.0 }; + Eigen::VectorXd x1(3); + x1 << -1.0, -2.0, 0.3; + Eigen::RowVectorXd x2(3); + x2 << -1.0, -2.0, 0.3; + stan::test::expect_ad(f1, x0); + stan::test::expect_ad(f1, x1); + stan::test::expect_ad(f1, x2); + stan::test::expect_ad_matvar(f1, x1); + stan::test::expect_ad_matvar(f1, x2); + + stan::test::expect_ad(f2, x0); + stan::test::expect_ad(f2, x1); + stan::test::expect_ad(f2, x2); + stan::test::expect_ad_matvar(f2, x1); + stan::test::expect_ad_matvar(f2, x2); + + std::vector x4; + Eigen::VectorXd x5(0); + Eigen::RowVectorXd x6(0); + + stan::test::expect_ad(f1, x4); + stan::test::expect_ad(f1, x5); + stan::test::expect_ad(f1, x6); + stan::test::expect_ad_matvar(f1, x5); + stan::test::expect_ad_matvar(f1, x6); + + stan::test::expect_ad(f2, x4); + stan::test::expect_ad(f2, x5); + stan::test::expect_ad(f2, x6); + stan::test::expect_ad_matvar(f2, x5); + stan::test::expect_ad_matvar(f2, x6); +} diff --git a/test/unit/math/mix/fun/exp_test.cpp b/test/unit/math/mix/fun/exp_test.cpp index 1a75ed232e8..28a19e80f46 100644 --- a/test/unit/math/mix/fun/exp_test.cpp +++ b/test/unit/math/mix/fun/exp_test.cpp @@ -9,4 +9,22 @@ TEST(mathMixMatFun, exp) { stan::test::expect_unary_vectorized(f, -15.2, -10, -0.5, 0.5, 1, 1.0, 1.3, 5, 10); stan::test::expect_complex_common(f); + + Eigen::VectorXd x1(3); + x1 << -1.0, 2.0, 3.0; + Eigen::RowVectorXd x2(3); + x2 << -1.0, 2.0, 3.0; + Eigen::MatrixXd x3(2, 3); + x3 << -1.0, 2.0, 3.0, 4.0, 5.0, 6.0; + stan::test::expect_ad_matvar(f, x1); + stan::test::expect_ad_matvar(f, x2); + stan::test::expect_ad_matvar(f, x3); + + Eigen::VectorXd x4(0); + Eigen::RowVectorXd x5(0); + Eigen::MatrixXd x6(0, 0); + + stan::test::expect_ad_matvar(f, x4); + stan::test::expect_ad_matvar(f, x5); + stan::test::expect_ad_matvar(f, x6); } diff --git a/test/unit/math/mix/fun/log1m_test.cpp b/test/unit/math/mix/fun/log1m_test.cpp index d52ba182a9c..a5c6e3542d1 100644 --- a/test/unit/math/mix/fun/log1m_test.cpp +++ b/test/unit/math/mix/fun/log1m_test.cpp @@ -4,4 +4,22 @@ TEST(mathMixMatFun, log1m) { auto f = [](const auto& x1) { return stan::math::log1m(x1); }; stan::test::expect_common_unary_vectorized(f); stan::test::expect_unary_vectorized(f, -21.5, -21, -1, 0.9, 3); + + Eigen::VectorXd x1(3); + x1 << -1.0, -2.0, 0.3; + Eigen::RowVectorXd x2(3); + x2 << -1.0, -2.0, 0.3; + Eigen::MatrixXd x3(2, 3); + x3 << -1.0, -2.0, 0.3, 0.4, -5.0, -6.0; + stan::test::expect_ad_matvar(f, x1); + stan::test::expect_ad_matvar(f, x2); + stan::test::expect_ad_matvar(f, x3); + + Eigen::VectorXd x4(0); + Eigen::RowVectorXd x5(0); + Eigen::MatrixXd x6(0, 0); + + stan::test::expect_ad_matvar(f, x4); + stan::test::expect_ad_matvar(f, x5); + stan::test::expect_ad_matvar(f, x6); } diff --git a/test/unit/math/mix/fun/read_corr_L_test.cpp b/test/unit/math/mix/fun/read_corr_L_test.cpp new file mode 100644 index 00000000000..c7b83e946c6 --- /dev/null +++ b/test/unit/math/mix/fun/read_corr_L_test.cpp @@ -0,0 +1,53 @@ +#include + +/*TEST(mathMixMatFun, read_corr_L) { + auto f = [](int K) { + return [K](const auto& x1) { + return stan::math::read_corr_L(x1, K); + }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + stan::test::expect_ad(f(4), x1); + stan::test::expect_ad(f(3), x2); + stan::test::expect_ad_matvar(f(4), x1); + stan::test::expect_ad_matvar(f(3), x2); + }*/ + +TEST(mathMixMatFun, read_corr_L_lp) { + auto f1 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + return stan::math::read_corr_L(x1, K, lp); + }; + }; + + auto f2 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + stan::math::read_corr_L(x1, K, lp); + return lp; + }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.0, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + Eigen::VectorXd x3(1); + x3 << 0.1; + stan::test::expect_ad(f1(4), x1); + stan::test::expect_ad(f1(3), x2); + + stan::test::expect_ad_matvar(f1(4), x1); + stan::test::expect_ad_matvar(f1(3), x2); + + stan::test::expect_ad(f2(4), x1); + stan::test::expect_ad(f2(3), x2); + + stan::test::expect_ad_matvar(f2(4), x1); + stan::test::expect_ad_matvar(f2(3), x2); +} diff --git a/test/unit/math/mix/fun/square_test.cpp b/test/unit/math/mix/fun/square_test.cpp index 2ef89b26009..fc8ac70efaa 100644 --- a/test/unit/math/mix/fun/square_test.cpp +++ b/test/unit/math/mix/fun/square_test.cpp @@ -5,4 +5,22 @@ TEST(mathMixMatFun, square) { stan::test::expect_common_unary_vectorized(f); stan::test::expect_unary_vectorized(f, -2.6, -1.0, -0.5, -0.2, 0.5, 1.3, 3, 5, 1e5); + + Eigen::VectorXd x1(3); + x1 << -1.0, 2.0, 3.0; + Eigen::RowVectorXd x2(3); + x2 << -1.0, 2.0, 3.0; + Eigen::MatrixXd x3(2, 3); + x3 << -1.0, 2.0, 3.0, 4.0, 5.0, 6.0; + stan::test::expect_ad_matvar(f, x1); + stan::test::expect_ad_matvar(f, x2); + stan::test::expect_ad_matvar(f, x3); + + Eigen::VectorXd x4(0); + Eigen::RowVectorXd x5(0); + Eigen::MatrixXd x6(0, 0); + + stan::test::expect_ad_matvar(f, x4); + stan::test::expect_ad_matvar(f, x5); + stan::test::expect_ad_matvar(f, x6); } diff --git a/test/unit/math/mix/fun/tanh_test.cpp b/test/unit/math/mix/fun/tanh_test.cpp index 78bca8cef01..4d857c41a5f 100644 --- a/test/unit/math/mix/fun/tanh_test.cpp +++ b/test/unit/math/mix/fun/tanh_test.cpp @@ -8,4 +8,22 @@ TEST(mathMixMatFun, tanh) { stan::test::expect_common_nonzero_unary_vectorized(f); stan::test::expect_unary_vectorized(f, -2.6, -2, -1.2, -0.5, 0.5, 1.5); stan::test::expect_complex_common(f); + + Eigen::VectorXd x1(3); + x1 << -1.0, 2.0, 3.0; + Eigen::RowVectorXd x2(3); + x2 << -1.0, 2.0, 3.0; + Eigen::MatrixXd x3(2, 3); + x3 << -1.0, 2.0, 3.0, 4.0, 5.0, 6.0; + stan::test::expect_ad_matvar(f, x1); + stan::test::expect_ad_matvar(f, x2); + stan::test::expect_ad_matvar(f, x3); + + Eigen::VectorXd x4(0); + Eigen::RowVectorXd x5(0); + Eigen::MatrixXd x6(0, 0); + + stan::test::expect_ad_matvar(f, x4); + stan::test::expect_ad_matvar(f, x5); + stan::test::expect_ad_matvar(f, x6); } From fe825bd4f93ce0d42f9ed75cf895357c65329b71 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 11 Dec 2020 14:59:52 +0000 Subject: [PATCH 02/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/exp.hpp | 8 ++-- stan/math/prim/fun/log1m.hpp | 3 +- stan/math/prim/fun/read_corr_L.hpp | 6 +-- stan/math/prim/fun/square.hpp | 2 +- stan/math/prim/fun/tanh.hpp | 2 +- stan/math/rev/fun/exp.hpp | 5 +-- stan/math/rev/fun/log1m.hpp | 5 +-- stan/math/rev/fun/read_corr_L.hpp | 40 +++++++++++-------- stan/math/rev/fun/square.hpp | 5 +-- stan/math/rev/fun/tanh.hpp | 8 ++-- .../unit/math/mix/fun/corr_constrain_test.cpp | 8 ++-- 11 files changed, 46 insertions(+), 46 deletions(-) diff --git a/stan/math/prim/fun/exp.hpp b/stan/math/prim/fun/exp.hpp index 29111a16069..384d8d3c7d4 100644 --- a/stan/math/prim/fun/exp.hpp +++ b/stan/math/prim/fun/exp.hpp @@ -40,10 +40,10 @@ struct exp_fun { * @return Elementwise application of exponentiation to the argument. */ template < - typename Container, - require_not_container_st* = nullptr, - require_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr, - require_not_var_matrix_t* = nullptr> + typename Container, + require_not_container_st* = nullptr, + require_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr, + require_not_var_matrix_t* = nullptr> inline auto exp(const Container& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/log1m.hpp b/stan/math/prim/fun/log1m.hpp index 3c2892601ad..924ba7a8615 100644 --- a/stan/math/prim/fun/log1m.hpp +++ b/stan/math/prim/fun/log1m.hpp @@ -67,8 +67,7 @@ struct log1m_fun { * @param x container * @return Natural log of 1 minus each value in x. */ -template * = nullptr> +template * = nullptr> inline auto log1m(const T& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/read_corr_L.hpp b/stan/math/prim/fun/read_corr_L.hpp index 51fa8238fb8..659bc872b0a 100644 --- a/stan/math/prim/fun/read_corr_L.hpp +++ b/stan/math/prim/fun/read_corr_L.hpp @@ -34,9 +34,9 @@ namespace math { * canonical partial correlations. */ template * = nullptr> -Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> read_corr_L(const T& CPCs, // on (-1, 1) - size_t K) { - +Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> read_corr_L( + const T& CPCs, // on (-1, 1) + size_t K) { using T_scalar = value_type_t; if (K == 0) { return {}; diff --git a/stan/math/prim/fun/square.hpp b/stan/math/prim/fun/square.hpp index 92162e0c790..abb05be9cce 100644 --- a/stan/math/prim/fun/square.hpp +++ b/stan/math/prim/fun/square.hpp @@ -49,7 +49,7 @@ struct square_fun { */ template * = nullptr, - require_not_var_matrix_t* = nullptr> + require_not_var_matrix_t* = nullptr> inline auto square(const Container& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/tanh.hpp b/stan/math/prim/fun/tanh.hpp index ec2ab88b684..e0a2233d180 100644 --- a/stan/math/prim/fun/tanh.hpp +++ b/stan/math/prim/fun/tanh.hpp @@ -37,7 +37,7 @@ struct tanh_fun { */ template * = nullptr, - require_not_var_matrix_t* = nullptr> + require_not_var_matrix_t* = nullptr> inline auto tanh(const Container& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/rev/fun/exp.hpp b/stan/math/rev/fun/exp.hpp index d6611158140..41000fc2ebd 100644 --- a/stan/math/rev/fun/exp.hpp +++ b/stan/math/rev/fun/exp.hpp @@ -64,11 +64,10 @@ inline std::complex exp(const std::complex& z) { * @param x argument * @return elementwise exponentiation of x */ -template * = nullptr> +template * = nullptr> inline auto exp(const T& x) { T res = x.val().array().exp().matrix(); - + reverse_pass_callback([x, res]() mutable { x.adj().noalias() += (res.val().array() * res.adj().array()).matrix(); }); diff --git a/stan/math/rev/fun/log1m.hpp b/stan/math/rev/fun/log1m.hpp index 2fdc4ec4213..8e9523bd933 100644 --- a/stan/math/rev/fun/log1m.hpp +++ b/stan/math/rev/fun/log1m.hpp @@ -35,11 +35,10 @@ inline var log1m(const var& a) { return var(new internal::log1m_vari(a.vi_)); } * @param x argument * @return elementwise log of 1 - x */ -template * = nullptr> +template * = nullptr> inline auto log1m(const T& x) { T res = stan::math::log1m(x.val()); - + reverse_pass_callback([x, res]() mutable { x.adj().noalias() += (res.adj().array() / (x.val().array() - 1.0)).matrix(); }); diff --git a/stan/math/rev/fun/read_corr_L.hpp b/stan/math/rev/fun/read_corr_L.hpp index 3064f4022a4..8c0aada980f 100644 --- a/stan/math/rev/fun/read_corr_L.hpp +++ b/stan/math/rev/fun/read_corr_L.hpp @@ -33,7 +33,7 @@ namespace math { * canonical partial correlations. */ template * = nullptr> -auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) +auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) using ret_type = var_value; if (K == 0) { @@ -69,13 +69,14 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) ret_type L_res = L; auto arena_acc = to_arena(acc); - + reverse_pass_callback([CPCs, arena_acc, K, L_res]() mutable { Eigen::ArrayXd acc_val = arena_acc; Eigen::ArrayXd acc_adj = Eigen::ArrayXd::Zero(K - 1); - acc_adj(K - 2) += 0.5 * L_res.adj()(K - 1, K - 1) / L_res.val()(K - 1, K - 1); - + acc_adj(K - 2) + += 0.5 * L_res.adj()(K - 1, K - 1) / L_res.val()(K - 1, K - 1); + int position = CPCs.size() - 1; for (size_t i = K - 2; i > 0; --i) { int pull = K - 1 - i; @@ -83,17 +84,22 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) const auto& cpc_seg_val = CPCs.val().array().segment(position, pull); auto cpc_seg_adj = CPCs.adj().array().segment(position, pull); acc_val.tail(pull) /= 1.0 - cpc_seg_val.square(); - cpc_seg_adj -= 2.0 * acc_adj.tail(pull) * acc_val.tail(pull) * cpc_seg_val; - acc_adj.tail(pull) = (acc_adj.tail(pull).array() * (1.0 - cpc_seg_val.square())).eval(); - cpc_seg_adj += L_res.adj().array().col(i).tail(pull) * acc_val.tail(pull).sqrt(); - acc_adj.tail(pull) += 0.5 * L_res.adj().array().col(i).tail(pull) * cpc_seg_val / acc_val.tail(pull).sqrt(); + cpc_seg_adj + -= 2.0 * acc_adj.tail(pull) * acc_val.tail(pull) * cpc_seg_val; + acc_adj.tail(pull) + = (acc_adj.tail(pull).array() * (1.0 - cpc_seg_val.square())).eval(); + cpc_seg_adj + += L_res.adj().array().col(i).tail(pull) * acc_val.tail(pull).sqrt(); + acc_adj.tail(pull) += 0.5 * L_res.adj().array().col(i).tail(pull) + * cpc_seg_val / acc_val.tail(pull).sqrt(); acc_adj(i - 1) += 0.5 * L_res.adj()(i, i) / L_res.val()(i, i); - + position -= pull + 1; } int pull = K - 1; - CPCs.adj().array().head(pull) -= 2.0 * acc_adj.tail(pull) * CPCs.val().array().head(pull); + CPCs.adj().array().head(pull) + -= 2.0 * acc_adj.tail(pull) * CPCs.val().array().head(pull); CPCs.adj().head(pull) += L_res.adj().col(0).tail(pull); }); @@ -125,9 +131,8 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) * @return Cholesky factor of correlation matrix for specified * partial correlations. */ -template * = nullptr, - require_stan_scalar_t* = nullptr> +template * = nullptr, + require_stan_scalar_t* = nullptr> auto read_corr_L(const T1& CPCs, size_t K, T2& log_prob) { using ret_val_type = Eigen::MatrixXd; using ret_type = var_value; @@ -158,13 +163,14 @@ auto read_corr_L(const T1& CPCs, size_t K, T2& log_prob) { size_t pos = CPCs.size() - 2; for (size_t k = K - 2; k >= 1; --k) { - for(size_t i = K; i >= k + 1; --i) { - CPCs.adj()(pos) -= 2.0 * acc_adj * (K - k - 1) * CPCs.val()(pos) / (1.0 - square(CPCs.val()(pos))); - pos--; + for (size_t i = K; i >= k + 1; --i) { + CPCs.adj()(pos) -= 2.0 * acc_adj * (K - k - 1) * CPCs.val()(pos) + / (1.0 - square(CPCs.val()(pos))); + pos--; } } }); - + return read_corr_L(CPCs, K); } diff --git a/stan/math/rev/fun/square.hpp b/stan/math/rev/fun/square.hpp index 83a0ba6f80d..0ce53d01654 100644 --- a/stan/math/rev/fun/square.hpp +++ b/stan/math/rev/fun/square.hpp @@ -51,11 +51,10 @@ inline var square(const var& x) { * @param x argument * @return elementwise square of x */ -template * = nullptr> +template * = nullptr> inline auto square(const T& x) { T res = (x.val().array() * x.val().array()).matrix(); - + reverse_pass_callback([x, res]() mutable { x.adj().noalias() += (2.0 * x.val().array() * res.adj().array()).matrix(); }); diff --git a/stan/math/rev/fun/tanh.hpp b/stan/math/rev/fun/tanh.hpp index 7769b9f49a8..369def13b43 100644 --- a/stan/math/rev/fun/tanh.hpp +++ b/stan/math/rev/fun/tanh.hpp @@ -68,14 +68,14 @@ inline std::complex tanh(const std::complex& z) { * @param x argument * @return elementwise hyperbolic tangent of x */ -template * = nullptr> +template * = nullptr> inline auto tanh(const T& x) { T res = stan::math::tanh(x.val()); - + reverse_pass_callback([x, res]() mutable { auto cosh = stan::math::cosh(x.val()); - x.adj().noalias() += (res.adj().array() / (cosh.array() * cosh.array())).matrix(); + x.adj().noalias() + += (res.adj().array() / (cosh.array() * cosh.array())).matrix(); }); return res; diff --git a/test/unit/math/mix/fun/corr_constrain_test.cpp b/test/unit/math/mix/fun/corr_constrain_test.cpp index 57ac6689a18..ff8e9f29b37 100644 --- a/test/unit/math/mix/fun/corr_constrain_test.cpp +++ b/test/unit/math/mix/fun/corr_constrain_test.cpp @@ -1,11 +1,9 @@ #include TEST(mathMixMatFun, corr_constrain) { - auto f = [](const auto& x1) { - return stan::math::corr_constrain(x1); - }; + auto f = [](const auto& x1) { return stan::math::corr_constrain(x1); }; - std::vector x0 = { -1.0, 2.0, 3.0 }; + std::vector x0 = {-1.0, 2.0, 3.0}; Eigen::VectorXd x1(3); x1 << -1.0, -2.0, 0.3; Eigen::RowVectorXd x2(3); @@ -40,7 +38,7 @@ TEST(mathMixMatFun, corr_constrain_lp) { return lp; }; - std::vector x0 = { -1.0, 2.0, 3.0 }; + std::vector x0 = {-1.0, 2.0, 3.0}; Eigen::VectorXd x1(3); x1 << -1.0, -2.0, 0.3; Eigen::RowVectorXd x2(3); From c5b9b4e12859b364c845d8d4555681e3a0c0d14a Mon Sep 17 00:00:00 2001 From: Ben Date: Sat, 12 Dec 2020 19:57:14 -0500 Subject: [PATCH 03/30] Added varmat support for `read_cov_L` (Issue #2101) --- stan/math/prim/fun/read_cov_L.hpp | 2 +- stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/read_cov_L.hpp | 55 +++++++++++++++++++++ test/unit/math/mix/fun/read_corr_L_test.cpp | 4 +- test/unit/math/mix/fun/read_cov_L_test.cpp | 43 ++++++++++++++++ 5 files changed, 102 insertions(+), 3 deletions(-) create mode 100644 stan/math/rev/fun/read_cov_L.hpp create mode 100644 test/unit/math/mix/fun/read_cov_L_test.cpp diff --git a/stan/math/prim/fun/read_cov_L.hpp b/stan/math/prim/fun/read_cov_L.hpp index 6dce25df5d4..dbc0734639d 100644 --- a/stan/math/prim/fun/read_cov_L.hpp +++ b/stan/math/prim/fun/read_cov_L.hpp @@ -27,7 +27,7 @@ namespace math { template * = nullptr, require_vt_same* = nullptr> -Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> read_cov_L( +Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> read_cov_L( const T_CPCs& CPCs, const T_sds& sds, value_type_t& log_prob) { size_t K = sds.rows(); // adjust due to transformation from correlations to covariances diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index 739dc064a24..fac8263c21d 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -127,6 +127,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/read_cov_L.hpp b/stan/math/rev/fun/read_cov_L.hpp new file mode 100644 index 00000000000..470693a1abb --- /dev/null +++ b/stan/math/rev/fun/read_cov_L.hpp @@ -0,0 +1,55 @@ +#ifndef STAN_MATH_REV_FUN_READ_COV_L_HPP +#define STAN_MATH_REV_FUN_READ_COV_L_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * This is the function that should be called prior to evaluating + * the density of any elliptical distribution + * + * @tparam T_CPCs type of \c T_CPCs vector (must be derived from \c + * Eigen::ArrayBase and have one compile-time dimension equal to 1) + * @tparam T_sds type of \c sds vector (must be derived from \c Eigen::ArrayBase + * and have one compile-time dimension equal to 1) + * @param CPCs on (-1, 1) + * @param sds on (0, inf) + * @param log_prob the log probability value to increment with the Jacobian + * @return Cholesky factor of covariance matrix for specified + * partial correlations. + */ +template * = nullptr, + require_vt_same* = nullptr> +inline auto read_cov_L(const T_CPCs& CPCs, const T_sds& sds, scalar_type_t& log_prob) { + size_t K = sds.rows(); + // adjust due to transformation from correlations to covariances + log_prob += (sum(log(sds.val())) + LOG_TWO) * K; + + auto corr_L = read_corr_L(CPCs, K, log_prob); + var_value res = sds.val().matrix().asDiagonal() * corr_L.val(); + + reverse_pass_callback([CPCs, sds, corr_L, log_prob, res]() mutable { + size_t K = sds.size(); + + corr_L.adj() += sds.val().matrix().asDiagonal() * res.adj(); + sds.adj() += rows_dot_product(res.adj(), corr_L.val()); + + sds.adj() += (K * log_prob.adj() / sds.val().array()).matrix(); + }); + + return res; +} + +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/mix/fun/read_corr_L_test.cpp b/test/unit/math/mix/fun/read_corr_L_test.cpp index c7b83e946c6..990217cd820 100644 --- a/test/unit/math/mix/fun/read_corr_L_test.cpp +++ b/test/unit/math/mix/fun/read_corr_L_test.cpp @@ -1,6 +1,6 @@ #include -/*TEST(mathMixMatFun, read_corr_L) { +TEST(mathMixMatFun, read_corr_L) { auto f = [](int K) { return [K](const auto& x1) { return stan::math::read_corr_L(x1, K); @@ -15,7 +15,7 @@ stan::test::expect_ad(f(3), x2); stan::test::expect_ad_matvar(f(4), x1); stan::test::expect_ad_matvar(f(3), x2); - }*/ +} TEST(mathMixMatFun, read_corr_L_lp) { auto f1 = [](int K) { diff --git a/test/unit/math/mix/fun/read_cov_L_test.cpp b/test/unit/math/mix/fun/read_cov_L_test.cpp new file mode 100644 index 00000000000..3e6e4808aaa --- /dev/null +++ b/test/unit/math/mix/fun/read_cov_L_test.cpp @@ -0,0 +1,43 @@ +#include + +TEST(mathMixMatFun, read_cov_L_lp) { + auto f1 = [](int K) { + Eigen::VectorXd rx2 = + (Eigen::VectorXd::Random(K).array() + 2.0).matrix(); + return [K, rx2](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + std::decay_t x2 = stan::math::add(x1.head(K), rx2); + return stan::math::read_cov_L(x1, x2, lp); + }; + }; + + auto f2 = [](int K) { + Eigen::VectorXd rx2 = + (Eigen::VectorXd::Random(K).array() * 0.0 + 2.0).matrix(); + return [K, rx2](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + auto x2 = stan::math::eval(stan::math::add(x1.head(K), rx2)); + stan::math::read_cov_L(x1, x2, lp); + return lp; + }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.0, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + Eigen::VectorXd x3(1); + x3 << 0.1; + + stan::test::expect_ad(f1(4), x1); + stan::test::expect_ad(f1(3), x2); + + stan::test::expect_ad_matvar(f1(4), x1); + stan::test::expect_ad_matvar(f1(3), x2); + + stan::test::expect_ad(f2(4), x1); + stan::test::expect_ad(f2(3), x2); + + stan::test::expect_ad_matvar(f2(4), x1); + stan::test::expect_ad_matvar(f2(3), x2); +} From 50757f2d22be3c1f1b8e9c1b75d1c6679e267130 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 13 Dec 2020 01:04:03 +0000 Subject: [PATCH 04/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/read_cov_L.hpp | 8 +++++--- test/unit/math/mix/fun/read_corr_L_test.cpp | 4 +--- test/unit/math/mix/fun/read_cov_L_test.cpp | 7 +++---- 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/stan/math/rev/fun/read_cov_L.hpp b/stan/math/rev/fun/read_cov_L.hpp index 470693a1abb..f8d64cdb4b8 100644 --- a/stan/math/rev/fun/read_cov_L.hpp +++ b/stan/math/rev/fun/read_cov_L.hpp @@ -29,13 +29,15 @@ namespace math { template * = nullptr, require_vt_same* = nullptr> -inline auto read_cov_L(const T_CPCs& CPCs, const T_sds& sds, scalar_type_t& log_prob) { +inline auto read_cov_L(const T_CPCs& CPCs, const T_sds& sds, + scalar_type_t& log_prob) { size_t K = sds.rows(); // adjust due to transformation from correlations to covariances log_prob += (sum(log(sds.val())) + LOG_TWO) * K; auto corr_L = read_corr_L(CPCs, K, log_prob); - var_value res = sds.val().matrix().asDiagonal() * corr_L.val(); + var_value res + = sds.val().matrix().asDiagonal() * corr_L.val(); reverse_pass_callback([CPCs, sds, corr_L, log_prob, res]() mutable { size_t K = sds.size(); @@ -45,7 +47,7 @@ inline auto read_cov_L(const T_CPCs& CPCs, const T_sds& sds, scalar_type_t> lp = 0.0; std::decay_t x2 = stan::math::add(x1.head(K), rx2); @@ -12,8 +11,8 @@ TEST(mathMixMatFun, read_cov_L_lp) { }; auto f2 = [](int K) { - Eigen::VectorXd rx2 = - (Eigen::VectorXd::Random(K).array() * 0.0 + 2.0).matrix(); + Eigen::VectorXd rx2 + = (Eigen::VectorXd::Random(K).array() * 0.0 + 2.0).matrix(); return [K, rx2](const auto& x1) { stan::scalar_type_t> lp = 0.0; auto x2 = stan::math::eval(stan::math::add(x1.head(K), rx2)); From 3015e1bccdaa5836ab0859cb231a9e6f2e0474c9 Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 13 Dec 2020 10:29:54 -0500 Subject: [PATCH 05/30] Added `read_corr_matrix` and `read_cov_matrix` (Issue #2101) --- stan/math/rev/fun.hpp | 2 + stan/math/rev/fun/read_corr_matrix.hpp | 67 +++++++++++++++++++ stan/math/rev/fun/read_cov_matrix.hpp | 64 ++++++++++++++++++ .../math/mix/fun/read_corr_matrix_test.cpp | 53 +++++++++++++++ .../math/mix/fun/read_cov_matrix_test.cpp | 67 +++++++++++++++++++ 5 files changed, 253 insertions(+) create mode 100644 stan/math/rev/fun/read_corr_matrix.hpp create mode 100644 stan/math/rev/fun/read_cov_matrix.hpp create mode 100644 test/unit/math/mix/fun/read_corr_matrix_test.cpp create mode 100644 test/unit/math/mix/fun/read_cov_matrix_test.cpp diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index fac8263c21d..eefb1f0bc57 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -127,7 +127,9 @@ #include #include #include +#include #include +#include #include #include #include diff --git a/stan/math/rev/fun/read_corr_matrix.hpp b/stan/math/rev/fun/read_corr_matrix.hpp new file mode 100644 index 00000000000..5a1835ff9bc --- /dev/null +++ b/stan/math/rev/fun/read_corr_matrix.hpp @@ -0,0 +1,67 @@ +#ifndef STAN_MATH_REV_FUN_READ_CORR_MATRIX_HPP +#define STAN_MATH_REV_FUN_READ_CORR_MATRIX_HPP + +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return the correlation matrix of the specified dimensionality + * corresponding to the specified canonical partial correlations. + * + *

See read_corr_matrix(Array, size_t, T) + * for more information. + * + * @tparam T_CPCs type of the array (must be derived from \c Eigen::ArrayBase + * and have one compile-time dimension equal to 1) + * @param CPCs The (K choose 2) canonical partial correlations in (-1, 1). + * @param K Dimensionality of correlation matrix. + * @return Cholesky factor of correlation matrix for specified + * canonical partial correlations. + */ +template * = nullptr> +inline var_value +read_corr_matrix(const T_CPCs& CPCs, size_t K) { + if (K == 0) { + return Eigen::MatrixXd(); + } + + return multiply_lower_tri_self_transpose(read_corr_L(CPCs, K)); +} + +/** + * Return the correlation matrix of the specified dimensionality + * corresponding to the specified canonical partial correlations, + * incrementing the specified scalar reference with the log + * absolute determinant of the Jacobian of the transformation. + * + * It is usually preferable to utilize the version that returns + * the Cholesky factor of the correlation matrix rather than the + * correlation matrix itself in statistical calculations. + * + * @tparam T_CPCs type of the array (must be derived from \c Eigen::ArrayBase + * and have one compile-time dimension equal to 1) + * @param CPCs The (K choose 2) canonical partial correlations in + * (-1, 1). + * @param K Dimensionality of correlation matrix. + * @param log_prob Reference to variable to increment with the log + * Jacobian determinant. + * @return Correlation matrix for specified partial correlations. + */ +template * = nullptr> +inline var_value +read_corr_matrix(const T_CPCs& CPCs, size_t K, scalar_type_t& log_prob) { + if (K == 0) { + return Eigen::MatrixXd(); + } + + return multiply_lower_tri_self_transpose(read_corr_L(CPCs, K, log_prob)); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun/read_cov_matrix.hpp b/stan/math/rev/fun/read_cov_matrix.hpp new file mode 100644 index 00000000000..06f9e75bbb6 --- /dev/null +++ b/stan/math/rev/fun/read_cov_matrix.hpp @@ -0,0 +1,64 @@ +#ifndef STAN_MATH_REV_FUN_READ_COV_MATRIX_HPP +#define STAN_MATH_REV_FUN_READ_COV_MATRIX_HPP + +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * A generally worse alternative to call prior to evaluating the + * density of an elliptical distribution + * + * @tparam T_CPCs type of \c T_CPCs vector (must be derived from \c + * Eigen::ArrayBase and have one compile-time dimension equal to 1) + * @tparam T_sds type of \c sds vector (must be derived from \c Eigen::ArrayBase + * and have one compile-time dimension equal to 1) + * @param CPCs on (-1, 1) + * @param sds on (0, inf) + * @param log_prob the log probability value to increment with the Jacobian + * @return Covariance matrix for specified partial correlations. + */ +template * = nullptr> +var_value +read_cov_matrix(const T_CPCs& CPCs, const T_sds& sds, + scalar_type_t& log_prob) { + return multiply_lower_tri_self_transpose(read_cov_L(CPCs, sds, log_prob)); +} + +/** + * Builds a covariance matrix from CPCs and standard deviations + * + * @tparam T_CPCs type of \c T_CPCs vector (must be derived from \c + * Eigen::ArrayBase and have one compile-time dimension equal to 1) + * @tparam T_sds type of \c sds vector (must be derived from \c Eigen::ArrayBase + * and have one compile-time dimension equal to 1) + * @param CPCs in (-1, 1) + * @param sds in (0, inf) + */ +template * = nullptr> +inline var_value +read_cov_matrix(const T_CPCs& CPCs, const T_sds& sds) { + size_t K = sds.rows(); + + var_value corr_L = read_corr_L(CPCs, K); + var_value prod = sds.val().matrix().asDiagonal() * corr_L.val(); + + reverse_pass_callback([sds, corr_L, prod]() mutable { + corr_L.adj() += sds.val().matrix().asDiagonal() * prod.adj(); + sds.adj() += rows_dot_product(prod.adj(), corr_L.val()); + }); + + return multiply_lower_tri_self_transpose(prod); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/mix/fun/read_corr_matrix_test.cpp b/test/unit/math/mix/fun/read_corr_matrix_test.cpp new file mode 100644 index 00000000000..5cf0e653f02 --- /dev/null +++ b/test/unit/math/mix/fun/read_corr_matrix_test.cpp @@ -0,0 +1,53 @@ +#include + +TEST(mathMixMatFun, read_corr_matrix) { + auto f = [](int K) { + return [K](const auto& x1) { + return stan::math::read_corr_matrix(x1, K); + }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + stan::test::expect_ad(f(4), x1); + stan::test::expect_ad(f(3), x2); + stan::test::expect_ad_matvar(f(4), x1); + stan::test::expect_ad_matvar(f(3), x2); +} + +TEST(mathMixMatFun, read_corr_matrix_lp) { + auto f1 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + return stan::math::read_corr_matrix(x1, K, lp); + }; + }; + + auto f2 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + stan::math::read_corr_matrix(x1, K, lp); + return lp; + }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.0, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + Eigen::VectorXd x3(1); + x3 << 0.1; + stan::test::expect_ad(f1(4), x1); + stan::test::expect_ad(f1(3), x2); + + stan::test::expect_ad_matvar(f1(4), x1); + stan::test::expect_ad_matvar(f1(3), x2); + + stan::test::expect_ad(f2(4), x1); + stan::test::expect_ad(f2(3), x2); + + stan::test::expect_ad_matvar(f2(4), x1); + stan::test::expect_ad_matvar(f2(3), x2); +} diff --git a/test/unit/math/mix/fun/read_cov_matrix_test.cpp b/test/unit/math/mix/fun/read_cov_matrix_test.cpp new file mode 100644 index 00000000000..27d272cd4cd --- /dev/null +++ b/test/unit/math/mix/fun/read_cov_matrix_test.cpp @@ -0,0 +1,67 @@ +#include + +TEST(mathMixMatFun, read_cov_matrix) { + auto f = [](int K) { + Eigen::VectorXd rx2 = + (Eigen::VectorXd::Random(K).array() + 2.0).matrix(); + return [K, rx2](const auto& x1) { + std::decay_t x2 = stan::math::add(x1.head(K), rx2); + return stan::math::read_cov_matrix(x1, x2); + }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.0, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + Eigen::VectorXd x3(1); + x3 << 0.1; + + stan::test::expect_ad(f(4), x1); + stan::test::expect_ad(f(3), x2); + + stan::test::expect_ad_matvar(f(4), x1); + stan::test::expect_ad_matvar(f(3), x2); +} + +TEST(mathMixMatFun, read_cov_matrix_lp) { + auto f1 = [](int K) { + Eigen::VectorXd rx2 = + (Eigen::VectorXd::Random(K).array() + 2.0).matrix(); + return [K, rx2](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + std::decay_t x2 = stan::math::add(x1.head(K), rx2); + return stan::math::read_cov_matrix(x1, x2, lp); + }; + }; + + auto f2 = [](int K) { + Eigen::VectorXd rx2 = + (Eigen::VectorXd::Random(K).array() * 0.0 + 2.0).matrix(); + return [K, rx2](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + auto x2 = stan::math::eval(stan::math::add(x1.head(K), rx2)); + stan::math::read_cov_matrix(x1, x2, lp); + return lp; + }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.0, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + Eigen::VectorXd x3(1); + x3 << 0.1; + + stan::test::expect_ad(f1(4), x1); + stan::test::expect_ad(f1(3), x2); + + stan::test::expect_ad_matvar(f1(4), x1); + stan::test::expect_ad_matvar(f1(3), x2); + + stan::test::expect_ad(f2(4), x1); + stan::test::expect_ad(f2(3), x2); + + stan::test::expect_ad_matvar(f2(4), x1); + stan::test::expect_ad_matvar(f2(3), x2); +} From b527fb487d04ba9b719ae7a0d3e6a4ba9574f341 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 13 Dec 2020 15:40:19 +0000 Subject: [PATCH 06/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/read_corr_matrix.hpp | 8 ++++---- stan/math/rev/fun/read_cov_matrix.hpp | 12 ++++++------ test/unit/math/mix/fun/read_corr_matrix_test.cpp | 4 +--- test/unit/math/mix/fun/read_cov_matrix_test.cpp | 10 ++++------ 4 files changed, 15 insertions(+), 19 deletions(-) diff --git a/stan/math/rev/fun/read_corr_matrix.hpp b/stan/math/rev/fun/read_corr_matrix.hpp index 5a1835ff9bc..5e30198ced4 100644 --- a/stan/math/rev/fun/read_corr_matrix.hpp +++ b/stan/math/rev/fun/read_corr_matrix.hpp @@ -23,8 +23,8 @@ namespace math { * canonical partial correlations. */ template * = nullptr> -inline var_value -read_corr_matrix(const T_CPCs& CPCs, size_t K) { +inline var_value read_corr_matrix(const T_CPCs& CPCs, + size_t K) { if (K == 0) { return Eigen::MatrixXd(); } @@ -52,8 +52,8 @@ read_corr_matrix(const T_CPCs& CPCs, size_t K) { * @return Correlation matrix for specified partial correlations. */ template * = nullptr> -inline var_value -read_corr_matrix(const T_CPCs& CPCs, size_t K, scalar_type_t& log_prob) { +inline var_value read_corr_matrix( + const T_CPCs& CPCs, size_t K, scalar_type_t& log_prob) { if (K == 0) { return Eigen::MatrixXd(); } diff --git a/stan/math/rev/fun/read_cov_matrix.hpp b/stan/math/rev/fun/read_cov_matrix.hpp index 06f9e75bbb6..9cdc47a3745 100644 --- a/stan/math/rev/fun/read_cov_matrix.hpp +++ b/stan/math/rev/fun/read_cov_matrix.hpp @@ -25,9 +25,8 @@ namespace math { */ template * = nullptr> -var_value -read_cov_matrix(const T_CPCs& CPCs, const T_sds& sds, - scalar_type_t& log_prob) { +var_value read_cov_matrix(const T_CPCs& CPCs, const T_sds& sds, + scalar_type_t& log_prob) { return multiply_lower_tri_self_transpose(read_cov_L(CPCs, sds, log_prob)); } @@ -43,12 +42,13 @@ read_cov_matrix(const T_CPCs& CPCs, const T_sds& sds, */ template * = nullptr> -inline var_value -read_cov_matrix(const T_CPCs& CPCs, const T_sds& sds) { +inline var_value read_cov_matrix(const T_CPCs& CPCs, + const T_sds& sds) { size_t K = sds.rows(); var_value corr_L = read_corr_L(CPCs, K); - var_value prod = sds.val().matrix().asDiagonal() * corr_L.val(); + var_value prod + = sds.val().matrix().asDiagonal() * corr_L.val(); reverse_pass_callback([sds, corr_L, prod]() mutable { corr_L.adj() += sds.val().matrix().asDiagonal() * prod.adj(); diff --git a/test/unit/math/mix/fun/read_corr_matrix_test.cpp b/test/unit/math/mix/fun/read_corr_matrix_test.cpp index 5cf0e653f02..e50a766bbf3 100644 --- a/test/unit/math/mix/fun/read_corr_matrix_test.cpp +++ b/test/unit/math/mix/fun/read_corr_matrix_test.cpp @@ -2,9 +2,7 @@ TEST(mathMixMatFun, read_corr_matrix) { auto f = [](int K) { - return [K](const auto& x1) { - return stan::math::read_corr_matrix(x1, K); - }; + return [K](const auto& x1) { return stan::math::read_corr_matrix(x1, K); }; }; Eigen::VectorXd x1(6); diff --git a/test/unit/math/mix/fun/read_cov_matrix_test.cpp b/test/unit/math/mix/fun/read_cov_matrix_test.cpp index 27d272cd4cd..8370f6e50d5 100644 --- a/test/unit/math/mix/fun/read_cov_matrix_test.cpp +++ b/test/unit/math/mix/fun/read_cov_matrix_test.cpp @@ -2,8 +2,7 @@ TEST(mathMixMatFun, read_cov_matrix) { auto f = [](int K) { - Eigen::VectorXd rx2 = - (Eigen::VectorXd::Random(K).array() + 2.0).matrix(); + Eigen::VectorXd rx2 = (Eigen::VectorXd::Random(K).array() + 2.0).matrix(); return [K, rx2](const auto& x1) { std::decay_t x2 = stan::math::add(x1.head(K), rx2); return stan::math::read_cov_matrix(x1, x2); @@ -26,8 +25,7 @@ TEST(mathMixMatFun, read_cov_matrix) { TEST(mathMixMatFun, read_cov_matrix_lp) { auto f1 = [](int K) { - Eigen::VectorXd rx2 = - (Eigen::VectorXd::Random(K).array() + 2.0).matrix(); + Eigen::VectorXd rx2 = (Eigen::VectorXd::Random(K).array() + 2.0).matrix(); return [K, rx2](const auto& x1) { stan::scalar_type_t> lp = 0.0; std::decay_t x2 = stan::math::add(x1.head(K), rx2); @@ -36,8 +34,8 @@ TEST(mathMixMatFun, read_cov_matrix_lp) { }; auto f2 = [](int K) { - Eigen::VectorXd rx2 = - (Eigen::VectorXd::Random(K).array() * 0.0 + 2.0).matrix(); + Eigen::VectorXd rx2 + = (Eigen::VectorXd::Random(K).array() * 0.0 + 2.0).matrix(); return [K, rx2](const auto& x1) { stan::scalar_type_t> lp = 0.0; auto x2 = stan::math::eval(stan::math::add(x1.head(K), rx2)); From 82be9cc24e570fb22b56c3908715b6d975634e2a Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 13 Dec 2020 10:48:33 -0500 Subject: [PATCH 07/30] Added `varmat` implementation of cov_matrix_constrain_lkj (Issue #2101) --- stan/math/rev/fun.hpp | 1 + .../math/rev/fun/cov_matrix_constrain_lkj.hpp | 76 +++++++++++++++++++ stan/math/rev/fun/exp.hpp | 2 +- stan/math/rev/fun/tanh.hpp | 2 +- .../mix/fun/cov_matrix_constrain_lkj_test.cpp | 51 +++++++++++++ 5 files changed, 130 insertions(+), 2 deletions(-) create mode 100644 stan/math/rev/fun/cov_matrix_constrain_lkj.hpp create mode 100644 test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index eefb1f0bc57..298a3ef7f21 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -35,6 +35,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp b/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp new file mode 100644 index 00000000000..2e55c2c8f70 --- /dev/null +++ b/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp @@ -0,0 +1,76 @@ +#ifndef STAN_MATH_REV_FUN_COV_MATRIX_CONSTRAIN_LKJ_HPP +#define STAN_MATH_REV_FUN_COV_MATRIX_CONSTRAIN_LKJ_HPP + +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return the covariance matrix of the specified dimensionality + * derived from constraining the specified vector of unconstrained + * values. The input vector must be of length \f$k \choose 2 + + * k\f$. The first \f$k \choose 2\f$ values in the input + * represent unconstrained (partial) correlations and the last + * \f$k\f$ are unconstrained standard deviations of the dimensions. + * + *

The transform scales the correlation matrix transform defined + * in corr_matrix_constrain(Matrix, size_t) + * with the constrained deviations. + * + * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and + * have one compile-time dimension equal to 1) + * @param x Input vector of unconstrained partial correlations and + * standard deviations. + * @param k Dimensionality of returned covariance matrix. + * @return Covariance matrix derived from the unconstrained partial + * correlations and deviations. + */ +template * = nullptr> +var_value cov_matrix_constrain_lkj(const T& x, size_t k) { + size_t k_choose_2 = (k * (k - 1)) / 2; + return read_cov_matrix(corr_constrain(x.head(k_choose_2)), + positive_constrain(x.tail(k))); +} + +/** + * Return the covariance matrix of the specified dimensionality + * derived from constraining the specified vector of unconstrained + * values and increment the specified log probability reference + * with the log absolute Jacobian determinant. + * + *

The transform is defined as for + * cov_matrix_constrain(Matrix, size_t). + * + *

The log absolute Jacobian determinant is derived by + * composing the log absolute Jacobian determinant for the + * underlying correlation matrix as defined in + * cov_matrix_constrain(Matrix, size_t, T&) with + * the Jacobian of the transform of the correlation matrix + * into a covariance matrix by scaling by standard deviations. + * + * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and + * have one compile-time dimension equal to 1) + * @param x Input vector of unconstrained partial correlations and + * standard deviations. + * @param k Dimensionality of returned covariance matrix. + * @param lp Log probability reference to increment. + * @return Covariance matrix derived from the unconstrained partial + * correlations and deviations. + */ +template * = nullptr> +var_value cov_matrix_constrain_lkj(const T& x, size_t k, scalar_type_t& lp) { + size_t k_choose_2 = (k * (k - 1)) / 2; + return read_cov_matrix(corr_constrain(x.head(k_choose_2)), + positive_constrain(x.tail(k)), lp); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun/exp.hpp b/stan/math/rev/fun/exp.hpp index 41000fc2ebd..c08477f9567 100644 --- a/stan/math/rev/fun/exp.hpp +++ b/stan/math/rev/fun/exp.hpp @@ -66,7 +66,7 @@ inline std::complex exp(const std::complex& z) { */ template * = nullptr> inline auto exp(const T& x) { - T res = x.val().array().exp().matrix(); + plain_type_t res = x.val().array().exp().matrix(); reverse_pass_callback([x, res]() mutable { x.adj().noalias() += (res.val().array() * res.adj().array()).matrix(); diff --git a/stan/math/rev/fun/tanh.hpp b/stan/math/rev/fun/tanh.hpp index 369def13b43..0e545ef3d8d 100644 --- a/stan/math/rev/fun/tanh.hpp +++ b/stan/math/rev/fun/tanh.hpp @@ -70,7 +70,7 @@ inline std::complex tanh(const std::complex& z) { */ template * = nullptr> inline auto tanh(const T& x) { - T res = stan::math::tanh(x.val()); + plain_type_t res = stan::math::tanh(x.val()); reverse_pass_callback([x, res]() mutable { auto cosh = stan::math::cosh(x.val()); diff --git a/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp b/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp new file mode 100644 index 00000000000..c570fb41ae7 --- /dev/null +++ b/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp @@ -0,0 +1,51 @@ +#include + +TEST(mathMixMatFun, read_corr_L) { + auto f = [](int K) { + return [K](const auto& x1) { return stan::math::cov_matrix_constrain_lkj(x1, K); }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + stan::test::expect_ad(f(4), x1); + stan::test::expect_ad(f(3), x2); + stan::test::expect_ad_matvar(f(4), x1); + stan::test::expect_ad_matvar(f(3), x2); +} + +TEST(mathMixMatFun, cov_matrix_constrain_lkj_lp) { + auto f1 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + return stan::math::cov_matrix_constrain_lkj(x1, K, lp); + }; + }; + + auto f2 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + stan::math::cov_matrix_constrain_lkj(x1, K, lp); + return lp; + }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.0, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + Eigen::VectorXd x3(1); + x3 << 0.1; + stan::test::expect_ad(f1(4), x1); + stan::test::expect_ad(f1(3), x2); + + stan::test::expect_ad_matvar(f1(4), x1); + stan::test::expect_ad_matvar(f1(3), x2); + + stan::test::expect_ad(f2(4), x1); + stan::test::expect_ad(f2(3), x2); + + stan::test::expect_ad_matvar(f2(4), x1); + stan::test::expect_ad_matvar(f2(3), x2); +} From 6125831e9a51c120886e67bed88aa7f344e088e6 Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 13 Dec 2020 10:57:59 -0500 Subject: [PATCH 08/30] Added varmat implementatin of corr_matrix_constrain (Issue #2101) --- stan/math/prim/fun/corr_matrix_constrain.hpp | 3 +- stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/corr_matrix_constrain.hpp | 82 +++++++++++++++++++ .../mix/fun/corr_matrix_constrain_test.cpp | 50 +++++++++++ .../mix/fun/cov_matrix_constrain_lkj_test.cpp | 2 +- 5 files changed, 136 insertions(+), 2 deletions(-) create mode 100644 stan/math/rev/fun/corr_matrix_constrain.hpp diff --git a/stan/math/prim/fun/corr_matrix_constrain.hpp b/stan/math/prim/fun/corr_matrix_constrain.hpp index 0bb81772eb6..39c2fdc7fb7 100644 --- a/stan/math/prim/fun/corr_matrix_constrain.hpp +++ b/stan/math/prim/fun/corr_matrix_constrain.hpp @@ -65,7 +65,8 @@ corr_matrix_constrain(const T& x, Eigen::Index k) { * @param k Dimensionality of returned correlation matrix. * @param lp Log probability reference to increment. */ -template +template * = nullptr> Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> corr_matrix_constrain(const T& x, Eigen::Index k, value_type_t& lp) { Eigen::Index k_choose_2 = (k * (k - 1)) / 2; diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index 298a3ef7f21..d7fddb8106d 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/corr_matrix_constrain.hpp b/stan/math/rev/fun/corr_matrix_constrain.hpp new file mode 100644 index 00000000000..7e6615dc26a --- /dev/null +++ b/stan/math/rev/fun/corr_matrix_constrain.hpp @@ -0,0 +1,82 @@ +#ifndef STAN_MATH_REV_FUN_CORR_MATRIX_CONSTRAIN_HPP +#define STAN_MATH_REV_FUN_CORR_MATRIX_CONSTRAIN_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return the correlation matrix of the specified dimensionality + * derived from the specified vector of unconstrained values. The + * input vector must be of length \f${k \choose 2} = + * \frac{k(k-1)}{2}\f$. The values in the input vector represent + * unconstrained (partial) correlations among the dimensions. + * + *

The transform based on partial correlations is as specified + * in + * + *

  • Lewandowski, Daniel, Dorota Kurowicka, and Harry + * Joe. 2009. Generating random correlation matrices based on + * vines and extended onion method. Journal of Multivariate + * Analysis 100:1989–-2001.
+ * + *

The free vector entries are first constrained to be + * valid correlation values using corr_constrain(T). + * + * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and + * have one compile-time dimension equal to 1) + * @param x Vector of unconstrained partial correlations. + * @param k Dimensionality of returned correlation matrix. + * @throw std::invalid_argument if x is not a valid correlation + * matrix. + */ +template * = nullptr> +var_value +corr_matrix_constrain(const T& x, Eigen::Index k) { + Eigen::Index k_choose_2 = (k * (k - 1)) / 2; + check_size_match("cov_matrix_constrain", "x.size()", x.size(), "k_choose_2", + k_choose_2); + return read_corr_matrix(corr_constrain(x), k); +} + +/** + * Return the correlation matrix of the specified dimensionality + * derived from the specified vector of unconstrained values. The + * input vector must be of length \f${k \choose 2} = + * \frac{k(k-1)}{2}\f$. The values in the input vector represent + * unconstrained (partial) correlations among the dimensions. + * + *

The transform is as specified for + * corr_matrix_constrain(Matrix, size_t); the + * paper it cites also defines the Jacobians for correlation inputs, + * which are composed with the correlation constrained Jacobians + * defined in corr_constrain(T, double) for + * this function. + * + * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and + * have one compile-time dimension equal to 1) + * @param x Vector of unconstrained partial correlations. + * @param k Dimensionality of returned correlation matrix. + * @param lp Log probability reference to increment. + */ +template * = nullptr> +var_value +corr_matrix_constrain(const T& x, Eigen::Index k, scalar_type_t& lp) { + Eigen::Index k_choose_2 = (k * (k - 1)) / 2; + check_size_match("cov_matrix_constrain", "x.size()", x.size(), "k_choose_2", + k_choose_2); + return read_corr_matrix(corr_constrain(x, lp), k, lp); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/mix/fun/corr_matrix_constrain_test.cpp b/test/unit/math/mix/fun/corr_matrix_constrain_test.cpp index 25d01fc4027..bf862d1a169 100644 --- a/test/unit/math/mix/fun/corr_matrix_constrain_test.cpp +++ b/test/unit/math/mix/fun/corr_matrix_constrain_test.cpp @@ -65,3 +65,53 @@ TEST(MathMixMatFun, corr_matrixTransform) { v6 << 1, 2, -3, 1.5, 0.2, 2; corr_matrix_constrain_test::expect_corr_matrix_transform(v6); } + +TEST(mathMixMatFun, corr_matrix_constrain) { + auto f = [](int K) { + return [K](const auto& x1) { return stan::math::corr_matrix_constrain(x1, K); }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + stan::test::expect_ad(f(4), x1); + stan::test::expect_ad(f(3), x2); + stan::test::expect_ad_matvar(f(4), x1); + stan::test::expect_ad_matvar(f(3), x2); +} + +TEST(mathMixMatFun, corr_matrix_constrain_lp) { + auto f1 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + return stan::math::corr_matrix_constrain(x1, K, lp); + }; + }; + + auto f2 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + stan::math::corr_matrix_constrain(x1, K, lp); + return lp; + }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.0, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + Eigen::VectorXd x3(1); + x3 << 0.1; + stan::test::expect_ad(f1(4), x1); + stan::test::expect_ad(f1(3), x2); + + stan::test::expect_ad_matvar(f1(4), x1); + stan::test::expect_ad_matvar(f1(3), x2); + + stan::test::expect_ad(f2(4), x1); + stan::test::expect_ad(f2(3), x2); + + stan::test::expect_ad_matvar(f2(4), x1); + stan::test::expect_ad_matvar(f2(3), x2); +} diff --git a/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp b/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp index c570fb41ae7..6c6c0f17611 100644 --- a/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp +++ b/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp @@ -1,6 +1,6 @@ #include -TEST(mathMixMatFun, read_corr_L) { +TEST(mathMixMatFun, cov_matrix_constrain_lkj) { auto f = [](int K) { return [K](const auto& x1) { return stan::math::cov_matrix_constrain_lkj(x1, K); }; }; From e118b1a4075d7e16001f43a52e2960cfe05dea98 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 13 Dec 2020 16:03:57 +0000 Subject: [PATCH 09/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/corr_matrix_constrain.hpp | 3 +-- stan/math/rev/fun/corr_matrix_constrain.hpp | 10 ++++------ stan/math/rev/fun/cov_matrix_constrain_lkj.hpp | 3 ++- test/unit/math/mix/fun/corr_matrix_constrain_test.cpp | 4 +++- .../math/mix/fun/cov_matrix_constrain_lkj_test.cpp | 4 +++- 5 files changed, 13 insertions(+), 11 deletions(-) diff --git a/stan/math/prim/fun/corr_matrix_constrain.hpp b/stan/math/prim/fun/corr_matrix_constrain.hpp index 39c2fdc7fb7..5b7e6f941cf 100644 --- a/stan/math/prim/fun/corr_matrix_constrain.hpp +++ b/stan/math/prim/fun/corr_matrix_constrain.hpp @@ -65,8 +65,7 @@ corr_matrix_constrain(const T& x, Eigen::Index k) { * @param k Dimensionality of returned correlation matrix. * @param lp Log probability reference to increment. */ -template * = nullptr> +template * = nullptr> Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> corr_matrix_constrain(const T& x, Eigen::Index k, value_type_t& lp) { Eigen::Index k_choose_2 = (k * (k - 1)) / 2; diff --git a/stan/math/rev/fun/corr_matrix_constrain.hpp b/stan/math/rev/fun/corr_matrix_constrain.hpp index 7e6615dc26a..d8b1d9c11a7 100644 --- a/stan/math/rev/fun/corr_matrix_constrain.hpp +++ b/stan/math/rev/fun/corr_matrix_constrain.hpp @@ -38,8 +38,7 @@ namespace math { * matrix. */ template * = nullptr> -var_value -corr_matrix_constrain(const T& x, Eigen::Index k) { +var_value corr_matrix_constrain(const T& x, Eigen::Index k) { Eigen::Index k_choose_2 = (k * (k - 1)) / 2; check_size_match("cov_matrix_constrain", "x.size()", x.size(), "k_choose_2", k_choose_2); @@ -66,10 +65,9 @@ corr_matrix_constrain(const T& x, Eigen::Index k) { * @param k Dimensionality of returned correlation matrix. * @param lp Log probability reference to increment. */ -template * = nullptr> -var_value -corr_matrix_constrain(const T& x, Eigen::Index k, scalar_type_t& lp) { +template * = nullptr> +var_value corr_matrix_constrain(const T& x, Eigen::Index k, + scalar_type_t& lp) { Eigen::Index k_choose_2 = (k * (k - 1)) / 2; check_size_match("cov_matrix_constrain", "x.size()", x.size(), "k_choose_2", k_choose_2); diff --git a/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp b/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp index 2e55c2c8f70..5f0473ede40 100644 --- a/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp +++ b/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp @@ -64,7 +64,8 @@ var_value cov_matrix_constrain_lkj(const T& x, size_t k) { * correlations and deviations. */ template * = nullptr> -var_value cov_matrix_constrain_lkj(const T& x, size_t k, scalar_type_t& lp) { +var_value cov_matrix_constrain_lkj(const T& x, size_t k, + scalar_type_t& lp) { size_t k_choose_2 = (k * (k - 1)) / 2; return read_cov_matrix(corr_constrain(x.head(k_choose_2)), positive_constrain(x.tail(k)), lp); diff --git a/test/unit/math/mix/fun/corr_matrix_constrain_test.cpp b/test/unit/math/mix/fun/corr_matrix_constrain_test.cpp index bf862d1a169..10dbb8990b5 100644 --- a/test/unit/math/mix/fun/corr_matrix_constrain_test.cpp +++ b/test/unit/math/mix/fun/corr_matrix_constrain_test.cpp @@ -68,7 +68,9 @@ TEST(MathMixMatFun, corr_matrixTransform) { TEST(mathMixMatFun, corr_matrix_constrain) { auto f = [](int K) { - return [K](const auto& x1) { return stan::math::corr_matrix_constrain(x1, K); }; + return [K](const auto& x1) { + return stan::math::corr_matrix_constrain(x1, K); + }; }; Eigen::VectorXd x1(6); diff --git a/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp b/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp index 6c6c0f17611..42fd1668ac9 100644 --- a/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp +++ b/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp @@ -2,7 +2,9 @@ TEST(mathMixMatFun, cov_matrix_constrain_lkj) { auto f = [](int K) { - return [K](const auto& x1) { return stan::math::cov_matrix_constrain_lkj(x1, K); }; + return [K](const auto& x1) { + return stan::math::cov_matrix_constrain_lkj(x1, K); + }; }; Eigen::VectorXd x1(6); From 77e5deb78b5282fcfdf3cf8ef0bff62a0a23800f Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 13 Dec 2020 12:13:57 -0500 Subject: [PATCH 10/30] Added varmat implementation of cov_matrix_constrain (Issue #2101) --- stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/cov_matrix_constrain.hpp | 121 ++++++++++++++++++ .../mix/fun/cov_matrix_constrain_test.cpp | 49 +++++++ 3 files changed, 171 insertions(+) create mode 100644 stan/math/rev/fun/cov_matrix_constrain.hpp diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index d7fddb8106d..e35d1398b50 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -35,6 +35,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/cov_matrix_constrain.hpp b/stan/math/rev/fun/cov_matrix_constrain.hpp new file mode 100644 index 00000000000..fb21c3f8c49 --- /dev/null +++ b/stan/math/rev/fun/cov_matrix_constrain.hpp @@ -0,0 +1,121 @@ +#ifndef STAN_MATH_REV_FUN_COV_MATRIX_CONSTRAIN_HPP +#define STAN_MATH_REV_FUN_COV_MATRIX_CONSTRAIN_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return the symmetric, positive-definite matrix of dimensions K + * by K resulting from transforming the specified finite vector of + * size K plus (K choose 2). + * + *

See cov_matrix_free() for the inverse transform. + * + * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and + * have one compile-time dimension equal to 1) + * @param x The vector to convert to a covariance matrix. + * @param K The number of rows and columns of the resulting + * covariance matrix. + * @throws std::invalid_argument if (x.size() != K + (K choose 2)). + */ +template * = nullptr> +var_value +cov_matrix_constrain(const T& x, Eigen::Index K) { + using std::exp; + + arena_t L_val(K, K); + check_size_match("cov_matrix_constrain", "x.size()", x.size(), + "K + (K choose 2)", (K * (K + 1)) / 2); + int i = 0; + for (Eigen::Index m = 0; m < K; ++m) { + L_val.row(m).head(m) = x.val().segment(i, m); + i += m; + L_val.coeffRef(m, m) = exp(x.val().coeff(i++)); + L_val.row(m).tail(K - m - 1).setZero(); + } + + var_value L = L_val; + + reverse_pass_callback([x, L]() mutable { + Eigen::Index K = L.rows(); + int i = x.size(); + for(int m = K - 1; m >= 0; --m) { + x.adj()(--i) += L.adj().coeff(m, m) * L.val().coeff(m, m); + i -= m; + x.adj().segment(i, m) += L.adj().row(m).head(m); + } + }); + + return multiply_lower_tri_self_transpose(L); +} + +/** + * Return the symmetric, positive-definite matrix of dimensions K + * by K resulting from transforming the specified finite vector of + * size K plus (K choose 2). + * + *

See cov_matrix_free() for the inverse transform. + * + * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and + * have one compile-time dimension equal to 1) + * @param x The vector to convert to a covariance matrix. + * @param K The dimensions of the resulting covariance matrix. + * @param lp Reference + * @throws std::domain_error if (x.size() != K + (K choose 2)). + */ +template * = nullptr> +var_value +cov_matrix_constrain(const T& x, Eigen::Index K, scalar_type_t& lp) { + using std::exp; + using std::log; + + arena_t L_val(K, K); + check_size_match("cov_matrix_constrain", "x.size()", x.size(), + "K + (K choose 2)", (K * (K + 1)) / 2); + int i = 0; + for (Eigen::Index m = 0; m < K; ++m) { + L_val.row(m).head(m) = x.val().segment(i, m); + i += m; + L_val.coeffRef(m, m) = exp(x.val().coeff(i++)); + L_val.row(m).tail(K - m - 1).setZero(); + } + + // Jacobian for complete transform, including exp() above + lp += (K * LOG_TWO); // needless constant; want propto + for (Eigen::Index k = 0; k < K; ++k) { + // only +1 because index from 0 + lp += (K - k + 1) * log(L_val.coeff(k, k)); + } + + var_value L = L_val; + + reverse_pass_callback([x, L, lp]() mutable { + Eigen::Index K = L.rows(); + for (Eigen::Index k = 0; k < K; ++k) { + L.adj().coeffRef(k, k) += (K - k + 1) * lp.adj() / L.val().coeff(k, k); + } + int i = x.size(); + for(int m = K - 1; m >= 0; --m) { + x.adj()(--i) += L.adj().coeff(m, m) * L.val().coeff(m, m); + i -= m; + x.adj().segment(i, m) += L.adj().row(m).head(m); + } + }); + + return multiply_lower_tri_self_transpose(L); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp b/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp index c9da8670c5d..17f88cd7c3a 100644 --- a/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp +++ b/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp @@ -74,3 +74,52 @@ TEST(MathMixMatFun, cov_matrixTransform) { v10 << 1, 2, -3, 1.7, 9.8, -12.2, 0.4, 0.2, 1.2, 2.7; cov_matrix_constrain_test::expect_cov_matrix_transform(v10); } + +TEST(mathMixMatFun, cov_matrix_constrain) { + auto f = [](int K) { + return [K](const auto& x1) { return stan::math::cov_matrix_constrain(x1, K); }; + }; + + Eigen::VectorXd x1(10); + x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3, -0.1, -0.2, -0.3, -0.4; + Eigen::VectorXd x2(6); + x2 << -0.3, 0.2, -0.99, 0.1, 0.2, 0.3; + stan::test::expect_ad(f(4), x1); + stan::test::expect_ad(f(3), x2); + stan::test::expect_ad_matvar(f(4), x1); + stan::test::expect_ad_matvar(f(3), x2); +} + +TEST(mathMixMatFun, cov_matrix_constrain_lp) { + auto f1 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + return stan::math::cov_matrix_constrain(x1, K, lp); + }; + }; + + auto f2 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + stan::math::cov_matrix_constrain(x1, K, lp); + return lp; + }; + }; + + Eigen::VectorXd x1(10); + x1 << -0.9, 0.0, 0.99, 0.1, 0.2, 0.3, -0.1, -0.2, -0.3, -0.4; + Eigen::VectorXd x2(6); + x2 << -0.3, 0.2, -0.99, 0.1, 0.2, 0.3; + stan::test::expect_ad(f1(4), x1); + stan::test::expect_ad(f1(3), x2); + + stan::test::expect_ad_matvar(f1(4), x1); + stan::test::expect_ad_matvar(f1(3), x2); + + stan::test::expect_ad(f2(4), x1); + stan::test::expect_ad(f2(3), x2); + + stan::test::expect_ad_matvar(f2(4), x1); + stan::test::expect_ad_matvar(f2(3), x2); +} + From 974c6d67dffaf016ebd11aba1a577b65b671f20d Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 13 Dec 2020 17:16:19 +0000 Subject: [PATCH 11/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/cov_matrix_constrain.hpp | 13 ++++++------- .../unit/math/mix/fun/cov_matrix_constrain_test.cpp | 4 ++-- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/stan/math/rev/fun/cov_matrix_constrain.hpp b/stan/math/rev/fun/cov_matrix_constrain.hpp index fb21c3f8c49..aff3135426e 100644 --- a/stan/math/rev/fun/cov_matrix_constrain.hpp +++ b/stan/math/rev/fun/cov_matrix_constrain.hpp @@ -29,8 +29,7 @@ namespace math { * @throws std::invalid_argument if (x.size() != K + (K choose 2)). */ template * = nullptr> -var_value -cov_matrix_constrain(const T& x, Eigen::Index K) { +var_value cov_matrix_constrain(const T& x, Eigen::Index K) { using std::exp; arena_t L_val(K, K); @@ -49,13 +48,13 @@ cov_matrix_constrain(const T& x, Eigen::Index K) { reverse_pass_callback([x, L]() mutable { Eigen::Index K = L.rows(); int i = x.size(); - for(int m = K - 1; m >= 0; --m) { + for (int m = K - 1; m >= 0; --m) { x.adj()(--i) += L.adj().coeff(m, m) * L.val().coeff(m, m); i -= m; x.adj().segment(i, m) += L.adj().row(m).head(m); } }); - + return multiply_lower_tri_self_transpose(L); } @@ -74,8 +73,8 @@ cov_matrix_constrain(const T& x, Eigen::Index K) { * @throws std::domain_error if (x.size() != K + (K choose 2)). */ template * = nullptr> -var_value -cov_matrix_constrain(const T& x, Eigen::Index K, scalar_type_t& lp) { +var_value cov_matrix_constrain(const T& x, Eigen::Index K, + scalar_type_t& lp) { using std::exp; using std::log; @@ -105,7 +104,7 @@ cov_matrix_constrain(const T& x, Eigen::Index K, scalar_type_t& lp) { L.adj().coeffRef(k, k) += (K - k + 1) * lp.adj() / L.val().coeff(k, k); } int i = x.size(); - for(int m = K - 1; m >= 0; --m) { + for (int m = K - 1; m >= 0; --m) { x.adj()(--i) += L.adj().coeff(m, m) * L.val().coeff(m, m); i -= m; x.adj().segment(i, m) += L.adj().row(m).head(m); diff --git a/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp b/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp index 17f88cd7c3a..e8785c3b91c 100644 --- a/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp +++ b/test/unit/math/mix/fun/cov_matrix_constrain_test.cpp @@ -77,7 +77,8 @@ TEST(MathMixMatFun, cov_matrixTransform) { TEST(mathMixMatFun, cov_matrix_constrain) { auto f = [](int K) { - return [K](const auto& x1) { return stan::math::cov_matrix_constrain(x1, K); }; + return + [K](const auto& x1) { return stan::math::cov_matrix_constrain(x1, K); }; }; Eigen::VectorXd x1(10); @@ -122,4 +123,3 @@ TEST(mathMixMatFun, cov_matrix_constrain_lp) { stan::test::expect_ad_matvar(f2(4), x1); stan::test::expect_ad_matvar(f2(3), x2); } - From a6ae96f031228f337faa59cb8dffdc0b6d7bbf83 Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 13 Dec 2020 13:26:42 -0500 Subject: [PATCH 12/30] Added varmat implementation of cholesky_factor_constrain (Issue #2101) --- stan/math/rev/fun.hpp | 1 + .../rev/fun/cholesky_factor_constrain.hpp | 117 ++++++++++++++++++ stan/math/rev/fun/cov_matrix_constrain.hpp | 6 +- .../fun/cholesky_factor_constrain_test.cpp | 52 ++++++++ 4 files changed, 174 insertions(+), 2 deletions(-) create mode 100644 stan/math/rev/fun/cholesky_factor_constrain.hpp create mode 100644 test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index e35d1398b50..01b9d7bbe90 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/cholesky_factor_constrain.hpp b/stan/math/rev/fun/cholesky_factor_constrain.hpp new file mode 100644 index 00000000000..2ca1fa8111f --- /dev/null +++ b/stan/math/rev/fun/cholesky_factor_constrain.hpp @@ -0,0 +1,117 @@ +#ifndef STAN_MATH_REV_FUN_CHOLESKY_FACTOR_CONSTRAIN_HPP +#define STAN_MATH_REV_FUN_CHOLESKY_FACTOR_CONSTRAIN_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Return the Cholesky factor of the specified size read from the + * specified vector. A total of (N choose 2) + N + (M - N) * N + * elements are required to read an M by N Cholesky factor. + * + * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and + * have one compile-time dimension equal to 1) + * @param x Vector of unconstrained values + * @param M number of rows + * @param N number of columns + * @return Cholesky factor + */ +template * = nullptr> +var_value +cholesky_factor_constrain(const T& x, int M, int N) { + using std::exp; + using T_scalar = value_type_t; + check_greater_or_equal("cholesky_factor_constrain", + "num rows (must be greater or equal to num cols)", M, + N); + check_size_match("cholesky_factor_constrain", "x.size()", x.size(), + "((N * (N + 1)) / 2 + (M - N) * N)", + ((N * (N + 1)) / 2 + (M - N) * N)); + arena_t y_val(M, N); + + int pos = 0; + for (int m = 0; m < N; ++m) { + y_val.row(m).head(m) = x.val().segment(pos, m); + pos += m; + y_val.coeffRef(m, m) = exp(x.val().coeff(pos++)); + y_val.row(m).tail(N - m - 1).setZero(); + } + + for (int m = N; m < M; ++m) { + y_val.row(m) = x.val().segment(pos, N); + pos += N; + } + + var_value y = y_val; + + reverse_pass_callback([x, M, N, y]() mutable { + int pos = x.size(); + for (int m = M - 1; m >= N; --m) { + pos -= N; + x.adj().segment(pos, N) += y.adj().row(m); + } + + for (int m = N - 1; m >= 0; --m) { + x.adj().coeffRef(--pos) += y.adj().coeff(m, m) * y.val().coeff(m, m); + pos -= m; + x.adj().segment(pos, m) += y.adj().row(m).head(m); + } + }); + + return y; +} + +/** + * Return the Cholesky factor of the specified size read from the + * specified vector and increment the specified log probability + * reference with the log Jacobian adjustment of the transform. A total + * of (N choose 2) + N + N * (M - N) free parameters are required to read + * an M by N Cholesky factor. + * + * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and + * have one compile-time dimension equal to 1) + * @param x Vector of unconstrained values + * @param M number of rows + * @param N number of columns + * @param lp Log probability that is incremented with the log Jacobian + * @return Cholesky factor + */ +template * = nullptr> +var_value +cholesky_factor_constrain(const T& x, int M, int N, scalar_type_t& lp) { + check_size_match("cholesky_factor_constrain", "x.size()", x.size(), + "((N * (N + 1)) / 2 + (M - N) * N)", + ((N * (N + 1)) / 2 + (M - N) * N)); + int pos = 0; + double lp_val = 0.0; + for (int n = 0; n < N; ++n) { + pos += n; + lp_val += x.val().coeff(pos++); + } + lp += lp_val; + + reverse_pass_callback([x, N, lp]() mutable { + int pos = 0; + for (int n = 0; n < N; ++n) { + pos += n; + x.adj().coeffRef(pos++) += lp.adj(); + } + }); + + return cholesky_factor_constrain(x, M, N); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun/cov_matrix_constrain.hpp b/stan/math/rev/fun/cov_matrix_constrain.hpp index fb21c3f8c49..e4d139f33d6 100644 --- a/stan/math/rev/fun/cov_matrix_constrain.hpp +++ b/stan/math/rev/fun/cov_matrix_constrain.hpp @@ -91,12 +91,14 @@ cov_matrix_constrain(const T& x, Eigen::Index K, scalar_type_t& lp) { } // Jacobian for complete transform, including exp() above - lp += (K * LOG_TWO); // needless constant; want propto + double lp_val = (K * LOG_TWO); // needless constant; want propto for (Eigen::Index k = 0; k < K; ++k) { // only +1 because index from 0 - lp += (K - k + 1) * log(L_val.coeff(k, k)); + lp_val += (K - k + 1) * log(L_val.coeff(k, k)); } + lp += lp_val; + var_value L = L_val; reverse_pass_callback([x, L, lp]() mutable { diff --git a/test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp b/test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp new file mode 100644 index 00000000000..f39f8ecac2e --- /dev/null +++ b/test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp @@ -0,0 +1,52 @@ +#include + +TEST(mathMixMatFun, cholesky_factor_constrain) { + auto f = [](int M, int N) { + return [M, N](const auto& x1) { + return stan::math::cholesky_factor_constrain(x1, M, N); + }; + }; + + Eigen::VectorXd x1(14); + x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3, -0.1, -0.2, -0.3, -0.4, -0.4, -0.5, -0.6, -0.7; + Eigen::VectorXd x2(12); + x2 << -0.3, 0.2, -0.99, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, -0.9; + stan::test::expect_ad(f(5, 4), x1); + stan::test::expect_ad(f(5, 3), x2); + stan::test::expect_ad_matvar(f(5, 4), x1); + stan::test::expect_ad_matvar(f(5, 3), x2); +} + +TEST(mathMixMatFun, cholesky_factor_constrain_lp) { + auto f1 = [](int M, int N) { + return [M, N](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + return stan::math::cholesky_factor_constrain(x1, M, N, lp); + }; + }; + + auto f2 = [](int M, int N) { + return [M, N](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + stan::math::cholesky_factor_constrain(x1, M, N, lp); + return lp; + }; + }; + + Eigen::VectorXd x1(14); + x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3, -0.1, -0.2, -0.3, -0.4, -0.4, -0.5, -0.6, -0.7; + Eigen::VectorXd x2(12); + x2 << -0.3, 0.2, -0.99, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, -0.9; + stan::test::expect_ad(f1(5, 4), x1); + stan::test::expect_ad(f1(5, 3), x2); + + stan::test::expect_ad_matvar(f1(5, 4), x1); + stan::test::expect_ad_matvar(f1(5, 3), x2); + + stan::test::expect_ad(f2(5, 4), x1); + stan::test::expect_ad(f2(5, 3), x2); + + stan::test::expect_ad_matvar(f2(5, 4), x1); + stan::test::expect_ad_matvar(f2(5, 3), x2); +} + From e89d629979cc05d9cf583bc0a8c47b0e57d4af7c Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 13 Dec 2020 18:27:46 +0000 Subject: [PATCH 13/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/cholesky_factor_constrain.hpp | 7 +++---- test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp | 7 ++++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/stan/math/rev/fun/cholesky_factor_constrain.hpp b/stan/math/rev/fun/cholesky_factor_constrain.hpp index 2ca1fa8111f..4ddabe20739 100644 --- a/stan/math/rev/fun/cholesky_factor_constrain.hpp +++ b/stan/math/rev/fun/cholesky_factor_constrain.hpp @@ -27,8 +27,7 @@ namespace math { * @return Cholesky factor */ template * = nullptr> -var_value -cholesky_factor_constrain(const T& x, int M, int N) { +var_value cholesky_factor_constrain(const T& x, int M, int N) { using std::exp; using T_scalar = value_type_t; check_greater_or_equal("cholesky_factor_constrain", @@ -87,8 +86,8 @@ cholesky_factor_constrain(const T& x, int M, int N) { * @return Cholesky factor */ template * = nullptr> -var_value -cholesky_factor_constrain(const T& x, int M, int N, scalar_type_t& lp) { +var_value cholesky_factor_constrain(const T& x, int M, int N, + scalar_type_t& lp) { check_size_match("cholesky_factor_constrain", "x.size()", x.size(), "((N * (N + 1)) / 2 + (M - N) * N)", ((N * (N + 1)) / 2 + (M - N) * N)); diff --git a/test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp b/test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp index f39f8ecac2e..0350fb4465b 100644 --- a/test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp +++ b/test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp @@ -8,7 +8,8 @@ TEST(mathMixMatFun, cholesky_factor_constrain) { }; Eigen::VectorXd x1(14); - x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3, -0.1, -0.2, -0.3, -0.4, -0.4, -0.5, -0.6, -0.7; + x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3, -0.1, -0.2, -0.3, -0.4, -0.4, -0.5, + -0.6, -0.7; Eigen::VectorXd x2(12); x2 << -0.3, 0.2, -0.99, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, -0.9; stan::test::expect_ad(f(5, 4), x1); @@ -34,7 +35,8 @@ TEST(mathMixMatFun, cholesky_factor_constrain_lp) { }; Eigen::VectorXd x1(14); - x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3, -0.1, -0.2, -0.3, -0.4, -0.4, -0.5, -0.6, -0.7; + x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3, -0.1, -0.2, -0.3, -0.4, -0.4, -0.5, + -0.6, -0.7; Eigen::VectorXd x2(12); x2 << -0.3, 0.2, -0.99, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, -0.9; stan::test::expect_ad(f1(5, 4), x1); @@ -49,4 +51,3 @@ TEST(mathMixMatFun, cholesky_factor_constrain_lp) { stan::test::expect_ad_matvar(f2(5, 4), x1); stan::test::expect_ad_matvar(f2(5, 3), x2); } - From 24e8aa13cf5e7d64a655816c68876050b0e97337 Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 13 Dec 2020 15:35:49 -0500 Subject: [PATCH 14/30] Add varmat implementation of cholesky_corr_constrain (Issue #2101) --- stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/cholesky_corr_constrain.hpp | 127 ++++++++++++++++++ .../mix/fun/cholesky_corr_constrain_test.cpp | 48 +++++++ 3 files changed, 176 insertions(+) create mode 100644 stan/math/rev/fun/cholesky_corr_constrain.hpp diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index 01b9d7bbe90..975a3135b6a 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/cholesky_corr_constrain.hpp b/stan/math/rev/fun/cholesky_corr_constrain.hpp new file mode 100644 index 00000000000..9ba4ed80e74 --- /dev/null +++ b/stan/math/rev/fun/cholesky_corr_constrain.hpp @@ -0,0 +1,127 @@ +#ifndef STAN_MATH_REV_FUN_CHOLESKY_CORR_CONSTRAIN_HPP +#define STAN_MATH_REV_FUN_CHOLESKY_CORR_CONSTRAIN_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +template * = nullptr> +var_value +cholesky_corr_constrain(const T& y, int K) { + using std::sqrt; + int k_choose_2 = (K * (K - 1)) / 2; + check_size_match("cholesky_corr_constrain", "y.size()", y.size(), + "k_choose_2", k_choose_2); + + if (K == 0) { + return Eigen::MatrixXd(); + } + + var_value z = corr_constrain(y); + arena_t x_val(K, K); + + x_val.setZero(); + x_val.coeffRef(0, 0) = 1; + int k = 0; + for (int i = 1; i < K; ++i) { + x_val.coeffRef(i, 0) = z.val().coeff(k++); + double sum_sqs = square(x_val.coeff(i, 0)); + for (int j = 1; j < i; ++j) { + x_val.coeffRef(i, j) = z.val().coeff(k++) * sqrt(1.0 - sum_sqs); + sum_sqs += square(x_val.coeff(i, j)); + } + x_val.coeffRef(i, i) = sqrt(1.0 - sum_sqs); + } + + var_value x = x_val; + + reverse_pass_callback([z, K, x]() mutable { + size_t k = z.size(); + for (int i = K - 1; i > 0; --i) { + double sum_sqs_val = 1.0 - square(x.val().coeffRef(i, i)); + double sum_sqs_adj = -0.5 * x.adj().coeff(i, i) / x.val().coeff(i, i); + for(int j = i - 1; j > 0; --j) { + x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); + sum_sqs_val -= square(x.val().coeff(i, j)); + sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * + z.val().coeff(--k) / sqrt(1.0 - sum_sqs_val); + z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); + } + x.adj().coeffRef(i, 0) += 2 * sum_sqs_adj * x.val().coeff(i, 0); + z.adj().coeffRef(--k) += x.adj().coeffRef(i, 0); + } + }); + + return x; +} + +// FIXME to match above after debugged +template * = nullptr> +var_value +cholesky_corr_constrain(const T& y, int K, scalar_type_t& lp) { + using Eigen::Dynamic; + using Eigen::Matrix; + using std::sqrt; + int k_choose_2 = (K * (K - 1)) / 2; + check_size_match("cholesky_corr_constrain", "y.size()", y.size(), + "k_choose_2", k_choose_2); + + if (K == 0) { + return Eigen::MatrixXd(); + } + + var_value z = corr_constrain(y, lp); + arena_t x_val(K, K); + + x_val.setZero(); + x_val.coeffRef(0, 0) = 1; + int k = 0; + double lp_val = 0.0; + for (int i = 1; i < K; ++i) { + x_val.coeffRef(i, 0) = z.val().coeff(k++); + double sum_sqs = square(x_val.coeff(i, 0)); + for (int j = 1; j < i; ++j) { + lp_val += 0.5 * log1m(sum_sqs); + x_val.coeffRef(i, j) = z.val().coeff(k++) * sqrt(1.0 - sum_sqs); + sum_sqs += square(x_val.coeff(i, j)); + } + x_val.coeffRef(i, i) = sqrt(1.0 - sum_sqs); + } + + lp += lp_val; + var_value x = x_val; + + reverse_pass_callback([z, K, x, lp]() mutable { + size_t k = z.size(); + for (int i = K - 1; i > 0; --i) { + double sum_sqs_val = 1.0 - square(x.val().coeffRef(i, i)); + double sum_sqs_adj = -0.5 * x.adj().coeff(i, i) / x.val().coeff(i, i); + for(int j = i - 1; j > 0; --j) { + x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); + sum_sqs_val -= square(x.val().coeff(i, j)); + + sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * + z.val().coeff(--k) / sqrt(1.0 - sum_sqs_val); + z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); + sum_sqs_adj -= 0.5 * lp.adj() / (1 - sum_sqs_val); + } + x.adj().coeffRef(i, 0) += 2 * sum_sqs_adj * x.val().coeff(i, 0); + z.adj().coeffRef(--k) += x.adj().coeffRef(i, 0); + } + }); + + return x; +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/mix/fun/cholesky_corr_constrain_test.cpp b/test/unit/math/mix/fun/cholesky_corr_constrain_test.cpp index 089b52256c0..87f1678cf8d 100644 --- a/test/unit/math/mix/fun/cholesky_corr_constrain_test.cpp +++ b/test/unit/math/mix/fun/cholesky_corr_constrain_test.cpp @@ -65,3 +65,51 @@ TEST(MathMixMatFun, cholesky_corrTransform) { v6 << 1, 2, -3, 1.5, 0.2, 2; cholesky_corr_constrain_test::expect_cholesky_corr_transform(v6); } + +TEST(mathMixMatFun, cholesky_corr_constrain) { + auto f = [](int K) { + return [K](const auto& x1) { return stan::math::cholesky_corr_constrain(x1, K); }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.2, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + stan::test::expect_ad(f(4), x1); + stan::test::expect_ad(f(3), x2); + stan::test::expect_ad_matvar(f(4), x1); + stan::test::expect_ad_matvar(f(3), x2); +} + +TEST(mathMixMatFun, cholesky_corr_constrain_lp) { + auto f1 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + return stan::math::cholesky_corr_constrain(x1, K, lp); + }; + }; + + auto f2 = [](int K) { + return [K](const auto& x1) { + stan::scalar_type_t> lp = 0.0; + stan::math::cholesky_corr_constrain(x1, K, lp); + return lp; + }; + }; + + Eigen::VectorXd x1(6); + x1 << -0.9, 0.0, 0.99, 0.1, 0.2, 0.3; + Eigen::VectorXd x2(3); + x2 << -0.3, 0.2, -0.99; + stan::test::expect_ad(f1(4), x1); + stan::test::expect_ad(f1(3), x2); + + stan::test::expect_ad_matvar(f1(4), x1); + stan::test::expect_ad_matvar(f1(3), x2); + + stan::test::expect_ad(f2(4), x1); + stan::test::expect_ad(f2(3), x2); + + stan::test::expect_ad_matvar(f2(4), x1); + stan::test::expect_ad_matvar(f2(3), x2); +} From 640978177bb7b46ea493e48d546aca1ba59a3af3 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 13 Dec 2020 20:36:59 +0000 Subject: [PATCH 15/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/cholesky_corr_constrain.hpp | 40 +++++++++---------- .../mix/fun/cholesky_corr_constrain_test.cpp | 4 +- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/stan/math/rev/fun/cholesky_corr_constrain.hpp b/stan/math/rev/fun/cholesky_corr_constrain.hpp index 9ba4ed80e74..010d6dd09d5 100644 --- a/stan/math/rev/fun/cholesky_corr_constrain.hpp +++ b/stan/math/rev/fun/cholesky_corr_constrain.hpp @@ -14,8 +14,7 @@ namespace stan { namespace math { template * = nullptr> -var_value -cholesky_corr_constrain(const T& y, int K) { +var_value cholesky_corr_constrain(const T& y, int K) { using std::sqrt; int k_choose_2 = (K * (K - 1)) / 2; check_size_match("cholesky_corr_constrain", "y.size()", y.size(), @@ -48,26 +47,25 @@ cholesky_corr_constrain(const T& y, int K) { for (int i = K - 1; i > 0; --i) { double sum_sqs_val = 1.0 - square(x.val().coeffRef(i, i)); double sum_sqs_adj = -0.5 * x.adj().coeff(i, i) / x.val().coeff(i, i); - for(int j = i - 1; j > 0; --j) { - x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); - sum_sqs_val -= square(x.val().coeff(i, j)); - sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * - z.val().coeff(--k) / sqrt(1.0 - sum_sqs_val); - z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); + for (int j = i - 1; j > 0; --j) { + x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); + sum_sqs_val -= square(x.val().coeff(i, j)); + sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * z.val().coeff(--k) + / sqrt(1.0 - sum_sqs_val); + z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); } x.adj().coeffRef(i, 0) += 2 * sum_sqs_adj * x.val().coeff(i, 0); z.adj().coeffRef(--k) += x.adj().coeffRef(i, 0); } }); - + return x; } // FIXME to match above after debugged -template * = nullptr> -var_value -cholesky_corr_constrain(const T& y, int K, scalar_type_t& lp) { +template * = nullptr> +var_value cholesky_corr_constrain(const T& y, int K, + scalar_type_t& lp) { using Eigen::Dynamic; using Eigen::Matrix; using std::sqrt; @@ -105,14 +103,14 @@ cholesky_corr_constrain(const T& y, int K, scalar_type_t& lp) { for (int i = K - 1; i > 0; --i) { double sum_sqs_val = 1.0 - square(x.val().coeffRef(i, i)); double sum_sqs_adj = -0.5 * x.adj().coeff(i, i) / x.val().coeff(i, i); - for(int j = i - 1; j > 0; --j) { - x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); - sum_sqs_val -= square(x.val().coeff(i, j)); - - sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * - z.val().coeff(--k) / sqrt(1.0 - sum_sqs_val); - z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); - sum_sqs_adj -= 0.5 * lp.adj() / (1 - sum_sqs_val); + for (int j = i - 1; j > 0; --j) { + x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); + sum_sqs_val -= square(x.val().coeff(i, j)); + + sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * z.val().coeff(--k) + / sqrt(1.0 - sum_sqs_val); + z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); + sum_sqs_adj -= 0.5 * lp.adj() / (1 - sum_sqs_val); } x.adj().coeffRef(i, 0) += 2 * sum_sqs_adj * x.val().coeff(i, 0); z.adj().coeffRef(--k) += x.adj().coeffRef(i, 0); diff --git a/test/unit/math/mix/fun/cholesky_corr_constrain_test.cpp b/test/unit/math/mix/fun/cholesky_corr_constrain_test.cpp index 87f1678cf8d..eadb61e5468 100644 --- a/test/unit/math/mix/fun/cholesky_corr_constrain_test.cpp +++ b/test/unit/math/mix/fun/cholesky_corr_constrain_test.cpp @@ -68,7 +68,9 @@ TEST(MathMixMatFun, cholesky_corrTransform) { TEST(mathMixMatFun, cholesky_corr_constrain) { auto f = [](int K) { - return [K](const auto& x1) { return stan::math::cholesky_corr_constrain(x1, K); }; + return [K](const auto& x1) { + return stan::math::cholesky_corr_constrain(x1, K); + }; }; Eigen::VectorXd x1(6); From 2ca1cde478ce57545a60009516bd7e98548e2781 Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 13 Dec 2020 17:21:10 -0500 Subject: [PATCH 16/30] Added varmat vectorized versions of log, multiply_log, inv_logit (Issue #2101) --- stan/math/prim/fun/inv_logit.hpp | 3 +- stan/math/prim/fun/log.hpp | 7 +-- stan/math/prim/fun/multiply_log.hpp | 4 +- stan/math/rev/fun/inv_logit.hpp | 18 +++++++ stan/math/rev/fun/log.hpp | 18 +++++++ stan/math/rev/fun/multiply_log.hpp | 16 +++++++ test/unit/math/mix/fun/inv_logit_test.cpp | 18 +++++++ test/unit/math/mix/fun/log_test.cpp | 18 +++++++ test/unit/math/mix/fun/multiply_log_test.cpp | 49 ++++++++++++++++++++ 9 files changed, 146 insertions(+), 5 deletions(-) diff --git a/stan/math/prim/fun/inv_logit.hpp b/stan/math/prim/fun/inv_logit.hpp index e9e565cc840..e5b5b2f3bd4 100644 --- a/stan/math/prim/fun/inv_logit.hpp +++ b/stan/math/prim/fun/inv_logit.hpp @@ -81,7 +81,8 @@ struct inv_logit_fun { * @param x container * @return Inverse logit applied to each value in x. */ -template +template * = nullptr> inline auto inv_logit(const T& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/log.hpp b/stan/math/prim/fun/log.hpp index 94252b27dee..a38790f303f 100644 --- a/stan/math/prim/fun/log.hpp +++ b/stan/math/prim/fun/log.hpp @@ -43,9 +43,10 @@ struct log_fun { * @return Elementwise application of natural log to the argument. */ template < - typename Container, - require_not_container_st* = nullptr, - require_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> + typename Container, + require_not_container_st* = nullptr, + require_not_var_matrix_t* = nullptr, + require_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> inline auto log(const Container& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/multiply_log.hpp b/stan/math/prim/fun/multiply_log.hpp index 664fb545327..74d060de4c3 100644 --- a/stan/math/prim/fun/multiply_log.hpp +++ b/stan/math/prim/fun/multiply_log.hpp @@ -69,7 +69,9 @@ inline return_type_t multiply_log(const T_a a, const T_b b) { * @param b Second input * @return multiply_log function applied to the two inputs. */ -template * = nullptr> +template * = nullptr, + require_all_not_var_matrix_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { return apply_scalar_binary( a, b, [&](const auto& c, const auto& d) { return multiply_log(c, d); }); diff --git a/stan/math/rev/fun/inv_logit.hpp b/stan/math/rev/fun/inv_logit.hpp index 827f29f8249..dff7eaecc7a 100644 --- a/stan/math/rev/fun/inv_logit.hpp +++ b/stan/math/rev/fun/inv_logit.hpp @@ -33,6 +33,24 @@ inline var inv_logit(const var& a) { return var(new internal::inv_logit_vari(a.vi_)); } +/** + * Return the inverse logit of the elements of x + * + * @tparam T type of x + * @param x argument + * @return elementwise inverse logit of x + */ +template * = nullptr> +inline auto inv_logit(const T& x) { + plain_type_t res = inv_logit(x.val()); + + reverse_pass_callback([x, res]() mutable { + x.adj().noalias() += (res.adj().array() * res.val().array() * (1.0 - res.val().array())).matrix(); + }); + + return res; +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/fun/log.hpp b/stan/math/rev/fun/log.hpp index 89fcd529539..d23aa017c39 100644 --- a/stan/math/rev/fun/log.hpp +++ b/stan/math/rev/fun/log.hpp @@ -64,6 +64,24 @@ inline std::complex log(const std::complex& z) { return internal::complex_log(z); } +/** + * Return the natural log of the elements of x + * + * @tparam T type of x + * @param x argument + * @return elementwise natural log of x + */ +template * = nullptr> +inline auto log(const T& x) { + plain_type_t res = x.val().array().log().matrix(); + + reverse_pass_callback([x, res]() mutable { + x.adj().noalias() += (res.adj().array() / x.val().array()).matrix(); + }); + + return res; +} + } // namespace math } // namespace stan #endif diff --git a/stan/math/rev/fun/multiply_log.hpp b/stan/math/rev/fun/multiply_log.hpp index 8f8546af335..c6d01dd88fa 100644 --- a/stan/math/rev/fun/multiply_log.hpp +++ b/stan/math/rev/fun/multiply_log.hpp @@ -4,6 +4,8 @@ #include #include #include +#include +#include #include #include #include @@ -105,6 +107,20 @@ inline var multiply_log(double a, const var& b) { return var(new internal::multiply_log_dv_vari(a, b.vi_)); } +template* = nullptr, + require_any_var_matrix_t* = nullptr> +inline auto multiply_log(const T1& a, const T2& b) { + return elt_multiply(a, log(b)); +} + +template* = nullptr, + require_any_stan_scalar_t* = nullptr> +inline auto multiply_log(const T1& a, const T2& b) { + return multiply(a, log(b)); +} + } // namespace math } // namespace stan #endif diff --git a/test/unit/math/mix/fun/inv_logit_test.cpp b/test/unit/math/mix/fun/inv_logit_test.cpp index ca0ae131503..f765cda4580 100644 --- a/test/unit/math/mix/fun/inv_logit_test.cpp +++ b/test/unit/math/mix/fun/inv_logit_test.cpp @@ -5,4 +5,22 @@ TEST(mathMixMatFun, invLogit) { stan::test::expect_common_unary_vectorized(f); stan::test::expect_unary_vectorized(f, -2.6, -2, -1.2, -0.2, 0.5, 1, 1.3, 1.5, 3); + + Eigen::VectorXd x1(3); + x1 << -1.0, 2.0, 3.0; + Eigen::RowVectorXd x2(3); + x2 << -1.0, 2.0, 3.0; + Eigen::MatrixXd x3(2, 3); + x3 << -1.0, 2.0, 3.0, 4.0, 5.0, 6.0; + stan::test::expect_ad_matvar(f, x1); + stan::test::expect_ad_matvar(f, x2); + stan::test::expect_ad_matvar(f, x3); + + Eigen::VectorXd x4(0); + Eigen::RowVectorXd x5(0); + Eigen::MatrixXd x6(0, 0); + + stan::test::expect_ad_matvar(f, x4); + stan::test::expect_ad_matvar(f, x5); + stan::test::expect_ad_matvar(f, x6); } diff --git a/test/unit/math/mix/fun/log_test.cpp b/test/unit/math/mix/fun/log_test.cpp index 21ff2c8fa03..ca5143ff3b2 100644 --- a/test/unit/math/mix/fun/log_test.cpp +++ b/test/unit/math/mix/fun/log_test.cpp @@ -24,4 +24,22 @@ TEST(mathMixMatFun, log) { stan::test::expect_ad(f, std::complex{-0.0, -2.1}); stan::test::expect_ad(f, std::complex{2.1, -0.0}); // (negative real and zero imaginary illegal) + + Eigen::VectorXd x1(3); + x1 << 0.1, 2.0, 3.0; + Eigen::RowVectorXd x2(3); + x2 << 0.1, 2.0, 3.0; + Eigen::MatrixXd x3(2, 3); + x3 << 0.1, 2.0, 3.0, 4.0, 5.0, 6.0; + stan::test::expect_ad_matvar(f, x1); + stan::test::expect_ad_matvar(f, x2); + stan::test::expect_ad_matvar(f, x3); + + Eigen::VectorXd x4(0); + Eigen::RowVectorXd x5(0); + Eigen::MatrixXd x6(0, 0); + + stan::test::expect_ad_matvar(f, x4); + stan::test::expect_ad_matvar(f, x5); + stan::test::expect_ad_matvar(f, x6); } diff --git a/test/unit/math/mix/fun/multiply_log_test.cpp b/test/unit/math/mix/fun/multiply_log_test.cpp index d928e171f03..61f4b9db250 100644 --- a/test/unit/math/mix/fun/multiply_log_test.cpp +++ b/test/unit/math/mix/fun/multiply_log_test.cpp @@ -29,4 +29,53 @@ TEST(mathMixScalFun, multiplyLog_vec) { Eigen::VectorXd in2(2); in2 << 0.5, 3.4; stan::test::expect_ad_vectorized_binary(f, in1, in2); + + Eigen::VectorXd x1(3); + x1 << 1.0, 2.0, 3.0; + Eigen::RowVectorXd x2(3); + x2 << 1.0, 2.0, 3.0; + Eigen::MatrixXd x3(2, 3); + x3 << 1.0, 2.0, 3.0, 4.0, 5.0, 6.0; + + stan::test::expect_ad(f, x1, x1); + stan::test::expect_ad(f, x1, 2.0); + stan::test::expect_ad(f, 3.0, x1); + stan::test::expect_ad(f, x2, x2); + stan::test::expect_ad(f, x2, 2.5); + stan::test::expect_ad(f, 3.5, x2); + stan::test::expect_ad(f, x3, x3); + stan::test::expect_ad(f, x3, 4.0); + stan::test::expect_ad(f, 5.0, x3); + stan::test::expect_ad_matvar(f, x1, x1); + stan::test::expect_ad_matvar(f, x1, 2.0); + stan::test::expect_ad_matvar(f, 3.0, x1); + stan::test::expect_ad_matvar(f, x2, x2); + stan::test::expect_ad_matvar(f, x2, 2.5); + stan::test::expect_ad_matvar(f, 3.5, x2); + stan::test::expect_ad_matvar(f, x3, x3); + stan::test::expect_ad_matvar(f, x3, 4.0); + stan::test::expect_ad_matvar(f, 5.0, x3); + + Eigen::VectorXd x4(0); + Eigen::RowVectorXd x5(0); + Eigen::MatrixXd x6(0, 0); + + stan::test::expect_ad(f, x4, x4); + stan::test::expect_ad(f, x4, 2.0); + stan::test::expect_ad(f, 3.0, x4); + stan::test::expect_ad(f, x5, x5); + stan::test::expect_ad(f, x5, 2.5); + stan::test::expect_ad(f, 3.5, x5); + stan::test::expect_ad(f, x6, x6); + stan::test::expect_ad(f, x6, 4.0); + stan::test::expect_ad(f, 5.0, x6); + stan::test::expect_ad_matvar(f, x4, x4); + stan::test::expect_ad_matvar(f, x4, 2.0); + stan::test::expect_ad_matvar(f, 3.0, x4); + stan::test::expect_ad_matvar(f, x5, x5); + stan::test::expect_ad_matvar(f, x5, 2.5); + stan::test::expect_ad_matvar(f, 3.5, x5); + stan::test::expect_ad_matvar(f, x6, x6); + stan::test::expect_ad_matvar(f, x6, 4.0); + stan::test::expect_ad_matvar(f, 5.0, x6); } From 965a63dc3b850349d135954384d0e8975bc11e49 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 13 Dec 2020 22:26:48 +0000 Subject: [PATCH 17/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/inv_logit.hpp | 3 +-- stan/math/prim/fun/log.hpp | 8 ++++---- stan/math/prim/fun/multiply_log.hpp | 5 ++--- stan/math/rev/fun/inv_logit.hpp | 4 +++- stan/math/rev/fun/multiply_log.hpp | 10 ++++------ 5 files changed, 14 insertions(+), 16 deletions(-) diff --git a/stan/math/prim/fun/inv_logit.hpp b/stan/math/prim/fun/inv_logit.hpp index e5b5b2f3bd4..1899cce4390 100644 --- a/stan/math/prim/fun/inv_logit.hpp +++ b/stan/math/prim/fun/inv_logit.hpp @@ -81,8 +81,7 @@ struct inv_logit_fun { * @param x container * @return Inverse logit applied to each value in x. */ -template * = nullptr> +template * = nullptr> inline auto inv_logit(const T& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/log.hpp b/stan/math/prim/fun/log.hpp index a38790f303f..ba256709b95 100644 --- a/stan/math/prim/fun/log.hpp +++ b/stan/math/prim/fun/log.hpp @@ -43,10 +43,10 @@ struct log_fun { * @return Elementwise application of natural log to the argument. */ template < - typename Container, - require_not_container_st* = nullptr, - require_not_var_matrix_t* = nullptr, - require_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> + typename Container, + require_not_container_st* = nullptr, + require_not_var_matrix_t* = nullptr, + require_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> inline auto log(const Container& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/multiply_log.hpp b/stan/math/prim/fun/multiply_log.hpp index 74d060de4c3..731f7ddbbba 100644 --- a/stan/math/prim/fun/multiply_log.hpp +++ b/stan/math/prim/fun/multiply_log.hpp @@ -69,9 +69,8 @@ inline return_type_t multiply_log(const T_a a, const T_b b) { * @param b Second input * @return multiply_log function applied to the two inputs. */ -template * = nullptr, - require_all_not_var_matrix_t* = nullptr> +template * = nullptr, + require_all_not_var_matrix_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { return apply_scalar_binary( a, b, [&](const auto& c, const auto& d) { return multiply_log(c, d); }); diff --git a/stan/math/rev/fun/inv_logit.hpp b/stan/math/rev/fun/inv_logit.hpp index dff7eaecc7a..03bbd635938 100644 --- a/stan/math/rev/fun/inv_logit.hpp +++ b/stan/math/rev/fun/inv_logit.hpp @@ -45,7 +45,9 @@ inline auto inv_logit(const T& x) { plain_type_t res = inv_logit(x.val()); reverse_pass_callback([x, res]() mutable { - x.adj().noalias() += (res.adj().array() * res.val().array() * (1.0 - res.val().array())).matrix(); + x.adj().noalias() + += (res.adj().array() * res.val().array() * (1.0 - res.val().array())) + .matrix(); }); return res; diff --git a/stan/math/rev/fun/multiply_log.hpp b/stan/math/rev/fun/multiply_log.hpp index c6d01dd88fa..fc51f2c156e 100644 --- a/stan/math/rev/fun/multiply_log.hpp +++ b/stan/math/rev/fun/multiply_log.hpp @@ -107,16 +107,14 @@ inline var multiply_log(double a, const var& b) { return var(new internal::multiply_log_dv_vari(a, b.vi_)); } -template* = nullptr, - require_any_var_matrix_t* = nullptr> +template * = nullptr, + require_any_var_matrix_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { return elt_multiply(a, log(b)); } -template* = nullptr, - require_any_stan_scalar_t* = nullptr> +template * = nullptr, + require_any_stan_scalar_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { return multiply(a, log(b)); } From c3f84147601b0e76c7a681e449b84c601ac6d8d3 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 16 Dec 2020 20:08:35 -0500 Subject: [PATCH 18/30] Cleaning up constrain functions (Issue #2101) --- stan/math/prim/fun/read_corr_L.hpp | 3 +- stan/math/rev/fun/cholesky_corr_constrain.hpp | 30 +++++++++++++--- .../rev/fun/cholesky_factor_constrain.hpp | 13 ++++--- stan/math/rev/fun/corr_matrix_constrain.hpp | 8 ++--- stan/math/rev/fun/cov_matrix_constrain.hpp | 18 +++++----- .../math/rev/fun/cov_matrix_constrain_lkj.hpp | 8 ++--- stan/math/rev/fun/read_corr_L.hpp | 36 +++++++++---------- stan/math/rev/fun/read_corr_matrix.hpp | 8 ++--- stan/math/rev/fun/read_cov_L.hpp | 8 ++--- stan/math/rev/fun/read_cov_matrix.hpp | 16 ++++----- 10 files changed, 82 insertions(+), 66 deletions(-) diff --git a/stan/math/prim/fun/read_corr_L.hpp b/stan/math/prim/fun/read_corr_L.hpp index 659bc872b0a..7093d9f4874 100644 --- a/stan/math/prim/fun/read_corr_L.hpp +++ b/stan/math/prim/fun/read_corr_L.hpp @@ -33,7 +33,8 @@ namespace math { * @return Cholesky factor of correlation matrix for specified * canonical partial correlations. */ -template * = nullptr> +template * = nullptr> Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> read_corr_L( const T& CPCs, // on (-1, 1) size_t K) { diff --git a/stan/math/rev/fun/cholesky_corr_constrain.hpp b/stan/math/rev/fun/cholesky_corr_constrain.hpp index 010d6dd09d5..124ba13b920 100644 --- a/stan/math/rev/fun/cholesky_corr_constrain.hpp +++ b/stan/math/rev/fun/cholesky_corr_constrain.hpp @@ -13,6 +13,17 @@ namespace stan { namespace math { +/** + * Return the Cholesky factor of the correlation matrix of the sepcified + * size read from the unconstrained vector `y`. A total of K choose 2 + * elements are required to build a K by K Cholesky factor. + * + * @tparam T type of input vector (must be a `var_value` where `S` + * inherits from EigenBase) + * @param y Vector of unconstrained values + * @param K number of rows + * @return Cholesky factor of correlation matrix + */ template * = nullptr> var_value cholesky_corr_constrain(const T& y, int K) { using std::sqrt; @@ -25,9 +36,8 @@ var_value cholesky_corr_constrain(const T& y, int K) { } var_value z = corr_constrain(y); - arena_t x_val(K, K); + arena_t x_val = Eigen::MatrixXd::Zero(K, K); - x_val.setZero(); x_val.coeffRef(0, 0) = 1; int k = 0; for (int i = 1; i < K; ++i) { @@ -62,7 +72,18 @@ var_value cholesky_corr_constrain(const T& y, int K) { return x; } -// FIXME to match above after debugged +/** + * Return the Cholesky factor of the correlation matrix of the sepcified + * size read from the unconstrained vector `y`. A total of K choose 2 + * elements are required to build a K by K Cholesky factor. + * + * @tparam T type of input vector (must be a `var_value` where `S` + * inherits from EigenBase) + * @param y Vector of unconstrained values + * @param K number of rows + * @param[out] lp Log density that is incremented with the log Jacobian + * @return Cholesky factor of correlation matrix + */ template * = nullptr> var_value cholesky_corr_constrain(const T& y, int K, scalar_type_t& lp) { @@ -78,9 +99,8 @@ var_value cholesky_corr_constrain(const T& y, int K, } var_value z = corr_constrain(y, lp); - arena_t x_val(K, K); + arena_t x_val = Eigen::MatrixXd::Zero(K, K); - x_val.setZero(); x_val.coeffRef(0, 0) = 1; int k = 0; double lp_val = 0.0; diff --git a/stan/math/rev/fun/cholesky_factor_constrain.hpp b/stan/math/rev/fun/cholesky_factor_constrain.hpp index 4ddabe20739..d1e70f10083 100644 --- a/stan/math/rev/fun/cholesky_factor_constrain.hpp +++ b/stan/math/rev/fun/cholesky_factor_constrain.hpp @@ -19,8 +19,8 @@ namespace math { * specified vector. A total of (N choose 2) + N + (M - N) * N * elements are required to read an M by N Cholesky factor. * - * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and - * have one compile-time dimension equal to 1) + * @tparam T type of input vector (must be a `var_value` where `S` + * inherits from EigenBase) * @param x Vector of unconstrained values * @param M number of rows * @param N number of columns @@ -36,14 +36,13 @@ var_value cholesky_factor_constrain(const T& x, int M, int N) { check_size_match("cholesky_factor_constrain", "x.size()", x.size(), "((N * (N + 1)) / 2 + (M - N) * N)", ((N * (N + 1)) / 2 + (M - N) * N)); - arena_t y_val(M, N); + arena_t y_val = Eigen::MatrixXd::Zero(M, N); int pos = 0; for (int m = 0; m < N; ++m) { y_val.row(m).head(m) = x.val().segment(pos, m); pos += m; y_val.coeffRef(m, m) = exp(x.val().coeff(pos++)); - y_val.row(m).tail(N - m - 1).setZero(); } for (int m = N; m < M; ++m) { @@ -77,12 +76,12 @@ var_value cholesky_factor_constrain(const T& x, int M, int N) { * of (N choose 2) + N + N * (M - N) free parameters are required to read * an M by N Cholesky factor. * - * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and - * have one compile-time dimension equal to 1) + * @tparam T type of input vector (must be a `var_value` where `S` + * inherits from EigenBase) * @param x Vector of unconstrained values * @param M number of rows * @param N number of columns - * @param lp Log probability that is incremented with the log Jacobian + * @param[out] lp Log density that is incremented with the log Jacobian * @return Cholesky factor */ template * = nullptr> diff --git a/stan/math/rev/fun/corr_matrix_constrain.hpp b/stan/math/rev/fun/corr_matrix_constrain.hpp index d8b1d9c11a7..f28d7ed5e01 100644 --- a/stan/math/rev/fun/corr_matrix_constrain.hpp +++ b/stan/math/rev/fun/corr_matrix_constrain.hpp @@ -30,8 +30,8 @@ namespace math { *

The free vector entries are first constrained to be * valid correlation values using corr_constrain(T). * - * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and - * have one compile-time dimension equal to 1) + * @tparam T type of input vector (must be a `var_value` where `S` + * inherits from EigenBase) * @param x Vector of unconstrained partial correlations. * @param k Dimensionality of returned correlation matrix. * @throw std::invalid_argument if x is not a valid correlation @@ -59,8 +59,8 @@ var_value corr_matrix_constrain(const T& x, Eigen::Index k) { * defined in corr_constrain(T, double) for * this function. * - * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and - * have one compile-time dimension equal to 1) + * @tparam T type of input vector (must be a `var_value` where `S` + * inherits from EigenBase) * @param x Vector of unconstrained partial correlations. * @param k Dimensionality of returned correlation matrix. * @param lp Log probability reference to increment. diff --git a/stan/math/rev/fun/cov_matrix_constrain.hpp b/stan/math/rev/fun/cov_matrix_constrain.hpp index 727fd9e9c48..e8f4ccf6abf 100644 --- a/stan/math/rev/fun/cov_matrix_constrain.hpp +++ b/stan/math/rev/fun/cov_matrix_constrain.hpp @@ -21,8 +21,8 @@ namespace math { * *

See cov_matrix_free() for the inverse transform. * - * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and - * have one compile-time dimension equal to 1) + * @tparam T type of input vector (must be a `var_value` where `S` + * inherits from EigenBase) * @param x The vector to convert to a covariance matrix. * @param K The number of rows and columns of the resulting * covariance matrix. @@ -32,15 +32,14 @@ template * = nullptr> var_value cov_matrix_constrain(const T& x, Eigen::Index K) { using std::exp; - arena_t L_val(K, K); check_size_match("cov_matrix_constrain", "x.size()", x.size(), "K + (K choose 2)", (K * (K + 1)) / 2); + arena_t L_val = Eigen::MatrixXd::Zero(K, K); int i = 0; for (Eigen::Index m = 0; m < K; ++m) { L_val.row(m).head(m) = x.val().segment(i, m); i += m; L_val.coeffRef(m, m) = exp(x.val().coeff(i++)); - L_val.row(m).tail(K - m - 1).setZero(); } var_value L = L_val; @@ -49,7 +48,7 @@ var_value cov_matrix_constrain(const T& x, Eigen::Index K) { Eigen::Index K = L.rows(); int i = x.size(); for (int m = K - 1; m >= 0; --m) { - x.adj()(--i) += L.adj().coeff(m, m) * L.val().coeff(m, m); + x.adj().coeffRef(--i) += L.adj().coeff(m, m) * L.val().coeff(m, m); i -= m; x.adj().segment(i, m) += L.adj().row(m).head(m); } @@ -65,8 +64,8 @@ var_value cov_matrix_constrain(const T& x, Eigen::Index K) { * *

See cov_matrix_free() for the inverse transform. * - * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and - * have one compile-time dimension equal to 1) + * @tparam T type of input vector (must be a `var_value` where `S` + * inherits from EigenBase) * @param x The vector to convert to a covariance matrix. * @param K The dimensions of the resulting covariance matrix. * @param lp Reference @@ -78,15 +77,14 @@ var_value cov_matrix_constrain(const T& x, Eigen::Index K, using std::exp; using std::log; - arena_t L_val(K, K); check_size_match("cov_matrix_constrain", "x.size()", x.size(), "K + (K choose 2)", (K * (K + 1)) / 2); + arena_t L_val = Eigen::MatrixXd::Zero(K, K); int i = 0; for (Eigen::Index m = 0; m < K; ++m) { L_val.row(m).head(m) = x.val().segment(i, m); i += m; L_val.coeffRef(m, m) = exp(x.val().coeff(i++)); - L_val.row(m).tail(K - m - 1).setZero(); } // Jacobian for complete transform, including exp() above @@ -107,7 +105,7 @@ var_value cov_matrix_constrain(const T& x, Eigen::Index K, } int i = x.size(); for (int m = K - 1; m >= 0; --m) { - x.adj()(--i) += L.adj().coeff(m, m) * L.val().coeff(m, m); + x.adj().coeffRef(--i) += L.adj().coeff(m, m) * L.val().coeff(m, m); i -= m; x.adj().segment(i, m) += L.adj().row(m).head(m); } diff --git a/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp b/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp index 5f0473ede40..29ef3a2b0fe 100644 --- a/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp +++ b/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp @@ -23,8 +23,8 @@ namespace math { * in corr_matrix_constrain(Matrix, size_t) * with the constrained deviations. * - * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and - * have one compile-time dimension equal to 1) + * @tparam T type of input vector (must be a `var_value` where `S` + * inherits from EigenBase) * @param x Input vector of unconstrained partial correlations and * standard deviations. * @param k Dimensionality of returned covariance matrix. @@ -54,8 +54,8 @@ var_value cov_matrix_constrain_lkj(const T& x, size_t k) { * the Jacobian of the transform of the correlation matrix * into a covariance matrix by scaling by standard deviations. * - * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and - * have one compile-time dimension equal to 1) + * @tparam T type of input vector (must be a `var_value` where `S` + * inherits from EigenBase) * @param x Input vector of unconstrained partial correlations and * standard deviations. * @param k Dimensionality of returned covariance matrix. diff --git a/stan/math/rev/fun/read_corr_L.hpp b/stan/math/rev/fun/read_corr_L.hpp index 8c0aada980f..6ffa34bd2e1 100644 --- a/stan/math/rev/fun/read_corr_L.hpp +++ b/stan/math/rev/fun/read_corr_L.hpp @@ -24,7 +24,7 @@ namespace math { *

See read_corr_matrix(Array, size_t, T) * for more information. * - * @tparam T type of input (must be a `var_value` where `T` + * @tparam T type of input vector (must be a `var_value` where `S` * inherits from EigenBase) * @param CPCs The (K choose 2) canonical partial correlations in * (-1, 1). @@ -44,38 +44,36 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) } using std::sqrt; - Eigen::ArrayXd acc(K - 1); - acc.setOnes(); + arena_t arena_acc = Eigen::ArrayXd::Ones(K - 1); // Cholesky factor of correlation matrix - Eigen::MatrixXd L = Eigen::MatrixXd::Zero(K, K); + arena_t L = Eigen::MatrixXd::Zero(K, K); size_t position = 0; size_t pull = K - 1; L(0, 0) = 1.0; L.col(0).tail(pull) = CPCs.val().head(pull); - acc.tail(pull) = 1.0 - CPCs.val().head(pull).array().square(); + arena_acc.tail(pull) = 1.0 - CPCs.val().head(pull).array().square(); for (size_t i = 1; i < (K - 1); i++) { position += pull; pull = K - 1 - i; const auto& cpc_seg = CPCs.val().array().segment(position, pull); - L(i, i) = sqrt(acc(i - 1)); - L.col(i).tail(pull) = cpc_seg * acc.tail(pull).sqrt(); - acc.tail(pull) = acc.tail(pull) * (1.0 - cpc_seg.square()); + L.coeffRef(i, i) = sqrt(arena_acc.coeff(i - 1)); + L.col(i).tail(pull) = cpc_seg * arena_acc.tail(pull).sqrt(); + arena_acc.tail(pull) = arena_acc.tail(pull) * (1.0 - cpc_seg.square()); } - L(K - 1, K - 1) = sqrt(acc(K - 2)); + L.coeffRef(K - 1, K - 1) = sqrt(arena_acc.coeff(K - 2)); - ret_type L_res = L; - auto arena_acc = to_arena(acc); + arena_t L_res = L; reverse_pass_callback([CPCs, arena_acc, K, L_res]() mutable { Eigen::ArrayXd acc_val = arena_acc; Eigen::ArrayXd acc_adj = Eigen::ArrayXd::Zero(K - 1); - acc_adj(K - 2) - += 0.5 * L_res.adj()(K - 1, K - 1) / L_res.val()(K - 1, K - 1); + acc_adj.coeffRef(K - 2) + += 0.5 * L_res.adj().coeff(K - 1, K - 1) / L_res.val().coeff(K - 1, K - 1); int position = CPCs.size() - 1; for (size_t i = K - 2; i > 0; --i) { @@ -92,7 +90,7 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) += L_res.adj().array().col(i).tail(pull) * acc_val.tail(pull).sqrt(); acc_adj.tail(pull) += 0.5 * L_res.adj().array().col(i).tail(pull) * cpc_seg_val / acc_val.tail(pull).sqrt(); - acc_adj(i - 1) += 0.5 * L_res.adj()(i, i) / L_res.val()(i, i); + acc_adj.coeffRef(i - 1) += 0.5 * L_res.adj().coeff(i, i) / L_res.val().coeff(i, i); position -= pull + 1; } @@ -103,7 +101,7 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) CPCs.adj().head(pull) += L_res.adj().col(0).tail(pull); }); - return L_res; + return ret_type(L_res); } /** @@ -121,7 +119,7 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) * extended onion method Journal of Multivariate Analysis 100 * (2009) 1989–2001 * - * @tparam T type of input (must be a `var_value` where `T` + * @tparam T type of input vector (must be a `var_value` where `S` * inherits from EigenBase) * @param CPCs The (K choose 2) canonical partial correlations in * (-1, 1). @@ -151,7 +149,7 @@ auto read_corr_L(const T1& CPCs, size_t K, T2& log_prob) { // see inverse of Jacobian in equation 11 of LKJ paper for (size_t k = 1; k <= (K - 2); k++) { for (size_t i = k + 1; i <= K; i++) { - acc_val += (K - k - 1) * log1m(square(CPCs.val()(pos))); + acc_val += (K - k - 1) * log1m(square(CPCs.val().coeff(pos))); pos++; } } @@ -164,8 +162,8 @@ auto read_corr_L(const T1& CPCs, size_t K, T2& log_prob) { size_t pos = CPCs.size() - 2; for (size_t k = K - 2; k >= 1; --k) { for (size_t i = K; i >= k + 1; --i) { - CPCs.adj()(pos) -= 2.0 * acc_adj * (K - k - 1) * CPCs.val()(pos) - / (1.0 - square(CPCs.val()(pos))); + CPCs.adj().coeffRef(pos) -= 2.0 * acc_adj * (K - k - 1) * + CPCs.val().coeff(pos) / (1.0 - square(CPCs.val().coeff(pos))); pos--; } } diff --git a/stan/math/rev/fun/read_corr_matrix.hpp b/stan/math/rev/fun/read_corr_matrix.hpp index 5e30198ced4..61ab3a05e65 100644 --- a/stan/math/rev/fun/read_corr_matrix.hpp +++ b/stan/math/rev/fun/read_corr_matrix.hpp @@ -15,8 +15,8 @@ namespace math { *

See read_corr_matrix(Array, size_t, T) * for more information. * - * @tparam T_CPCs type of the array (must be derived from \c Eigen::ArrayBase - * and have one compile-time dimension equal to 1) + * @tparam T_CPCs type of input vector (must be a `var_value` where `T` + * inherits from EigenBase) * @param CPCs The (K choose 2) canonical partial correlations in (-1, 1). * @param K Dimensionality of correlation matrix. * @return Cholesky factor of correlation matrix for specified @@ -42,8 +42,8 @@ inline var_value read_corr_matrix(const T_CPCs& CPCs, * the Cholesky factor of the correlation matrix rather than the * correlation matrix itself in statistical calculations. * - * @tparam T_CPCs type of the array (must be derived from \c Eigen::ArrayBase - * and have one compile-time dimension equal to 1) + * @tparam T_CPCs type of input vector (must be a `var_value` where `T` + * inherits from EigenBase) * @param CPCs The (K choose 2) canonical partial correlations in * (-1, 1). * @param K Dimensionality of correlation matrix. diff --git a/stan/math/rev/fun/read_cov_L.hpp b/stan/math/rev/fun/read_cov_L.hpp index f8d64cdb4b8..aec016967cb 100644 --- a/stan/math/rev/fun/read_cov_L.hpp +++ b/stan/math/rev/fun/read_cov_L.hpp @@ -16,10 +16,10 @@ namespace math { * This is the function that should be called prior to evaluating * the density of any elliptical distribution * - * @tparam T_CPCs type of \c T_CPCs vector (must be derived from \c - * Eigen::ArrayBase and have one compile-time dimension equal to 1) - * @tparam T_sds type of \c sds vector (must be derived from \c Eigen::ArrayBase - * and have one compile-time dimension equal to 1) + * @tparam T_CPCs type of CPCs vector (must be a `var_value` where `T` + * inherits from EigenBase) + * @tparam T_sds type of sds vector (must be a `var_value` where `T` + * inherits from EigenBase) * @param CPCs on (-1, 1) * @param sds on (0, inf) * @param log_prob the log probability value to increment with the Jacobian diff --git a/stan/math/rev/fun/read_cov_matrix.hpp b/stan/math/rev/fun/read_cov_matrix.hpp index 9cdc47a3745..2bc946f3ba1 100644 --- a/stan/math/rev/fun/read_cov_matrix.hpp +++ b/stan/math/rev/fun/read_cov_matrix.hpp @@ -14,10 +14,10 @@ namespace math { * A generally worse alternative to call prior to evaluating the * density of an elliptical distribution * - * @tparam T_CPCs type of \c T_CPCs vector (must be derived from \c - * Eigen::ArrayBase and have one compile-time dimension equal to 1) - * @tparam T_sds type of \c sds vector (must be derived from \c Eigen::ArrayBase - * and have one compile-time dimension equal to 1) + * @tparam T_CPCs type of CPCs vector (must be a `var_value` where `T` + * inherits from EigenBase) + * @tparam T_sds type of sds vector (must be a `var_value` where `T` + * inherits from EigenBase) * @param CPCs on (-1, 1) * @param sds on (0, inf) * @param log_prob the log probability value to increment with the Jacobian @@ -33,10 +33,10 @@ var_value read_cov_matrix(const T_CPCs& CPCs, const T_sds& sds, /** * Builds a covariance matrix from CPCs and standard deviations * - * @tparam T_CPCs type of \c T_CPCs vector (must be derived from \c - * Eigen::ArrayBase and have one compile-time dimension equal to 1) - * @tparam T_sds type of \c sds vector (must be derived from \c Eigen::ArrayBase - * and have one compile-time dimension equal to 1) + * @tparam T_CPCs type of CPCs vector (must be a `var_value` where `T` + * inherits from EigenBase) + * @tparam T_sds type of sds vector (must be a `var_value` where `T` + * inherits from EigenBase) * @param CPCs in (-1, 1) * @param sds in (0, inf) */ From 65fbb1c602b28a5410459e05655facea91294ced Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 17 Dec 2020 01:30:14 +0000 Subject: [PATCH 19/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/read_corr_L.hpp | 3 +-- stan/math/prim/fun/tanh.hpp | 5 ++--- stan/math/rev/fun/cholesky_corr_constrain.hpp | 2 +- stan/math/rev/fun/read_corr_L.hpp | 12 +++++++----- 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/stan/math/prim/fun/read_corr_L.hpp b/stan/math/prim/fun/read_corr_L.hpp index 7093d9f4874..659bc872b0a 100644 --- a/stan/math/prim/fun/read_corr_L.hpp +++ b/stan/math/prim/fun/read_corr_L.hpp @@ -33,8 +33,7 @@ namespace math { * @return Cholesky factor of correlation matrix for specified * canonical partial correlations. */ -template * = nullptr> +template * = nullptr> Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> read_corr_L( const T& CPCs, // on (-1, 1) size_t K) { diff --git a/stan/math/prim/fun/tanh.hpp b/stan/math/prim/fun/tanh.hpp index 4f86c371278..444ba696978 100644 --- a/stan/math/prim/fun/tanh.hpp +++ b/stan/math/prim/fun/tanh.hpp @@ -38,9 +38,8 @@ struct tanh_fun { template * = nullptr, require_not_var_matrix_t* = nullptr> - require_all_not_nonscalar_prim_or_rev_kernel_expression_t< - Container>* = nullptr> -inline auto tanh(const Container& x) { +require_all_not_nonscalar_prim_or_rev_kernel_expression_t< + Container>* = nullptr > inline auto tanh(const Container& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/rev/fun/cholesky_corr_constrain.hpp b/stan/math/rev/fun/cholesky_corr_constrain.hpp index 124ba13b920..104465204b7 100644 --- a/stan/math/rev/fun/cholesky_corr_constrain.hpp +++ b/stan/math/rev/fun/cholesky_corr_constrain.hpp @@ -81,7 +81,7 @@ var_value cholesky_corr_constrain(const T& y, int K) { * inherits from EigenBase) * @param y Vector of unconstrained values * @param K number of rows - * @param[out] lp Log density that is incremented with the log Jacobian + * @param[out] lp Log density that is incremented with the log Jacobian * @return Cholesky factor of correlation matrix */ template * = nullptr> diff --git a/stan/math/rev/fun/read_corr_L.hpp b/stan/math/rev/fun/read_corr_L.hpp index 6ffa34bd2e1..ac66e09b82a 100644 --- a/stan/math/rev/fun/read_corr_L.hpp +++ b/stan/math/rev/fun/read_corr_L.hpp @@ -72,8 +72,8 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) Eigen::ArrayXd acc_val = arena_acc; Eigen::ArrayXd acc_adj = Eigen::ArrayXd::Zero(K - 1); - acc_adj.coeffRef(K - 2) - += 0.5 * L_res.adj().coeff(K - 1, K - 1) / L_res.val().coeff(K - 1, K - 1); + acc_adj.coeffRef(K - 2) += 0.5 * L_res.adj().coeff(K - 1, K - 1) + / L_res.val().coeff(K - 1, K - 1); int position = CPCs.size() - 1; for (size_t i = K - 2; i > 0; --i) { @@ -90,7 +90,8 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) += L_res.adj().array().col(i).tail(pull) * acc_val.tail(pull).sqrt(); acc_adj.tail(pull) += 0.5 * L_res.adj().array().col(i).tail(pull) * cpc_seg_val / acc_val.tail(pull).sqrt(); - acc_adj.coeffRef(i - 1) += 0.5 * L_res.adj().coeff(i, i) / L_res.val().coeff(i, i); + acc_adj.coeffRef(i - 1) + += 0.5 * L_res.adj().coeff(i, i) / L_res.val().coeff(i, i); position -= pull + 1; } @@ -162,8 +163,9 @@ auto read_corr_L(const T1& CPCs, size_t K, T2& log_prob) { size_t pos = CPCs.size() - 2; for (size_t k = K - 2; k >= 1; --k) { for (size_t i = K; i >= k + 1; --i) { - CPCs.adj().coeffRef(pos) -= 2.0 * acc_adj * (K - k - 1) * - CPCs.val().coeff(pos) / (1.0 - square(CPCs.val().coeff(pos))); + CPCs.adj().coeffRef(pos) -= 2.0 * acc_adj * (K - k - 1) + * CPCs.val().coeff(pos) + / (1.0 - square(CPCs.val().coeff(pos))); pos--; } } From 44322d49edee809162892c63a713e58d4f58ffff Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 20 Dec 2020 19:37:19 +0000 Subject: [PATCH 20/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/multiply_log.hpp | 115 ++++++++++++++++------------- 1 file changed, 64 insertions(+), 51 deletions(-) diff --git a/stan/math/rev/fun/multiply_log.hpp b/stan/math/rev/fun/multiply_log.hpp index 9ea5a1a6530..b0e08896426 100644 --- a/stan/math/rev/fun/multiply_log.hpp +++ b/stan/math/rev/fun/multiply_log.hpp @@ -47,9 +47,7 @@ class multiply_log_dv_vari : public op_dv_vari { public: multiply_log_dv_vari(double a, vari* bvi) : op_dv_vari(multiply_log(a, bvi->val_), a, bvi) {} - void chain() { - bvi_->adj_ += adj_ * ad_ / bvi_->val_; - } + void chain() { bvi_->adj_ += adj_ * ad_ / bvi_->val_; } }; } // namespace internal @@ -102,95 +100,110 @@ inline var multiply_log(double a, const var& b) { template * = nullptr, require_any_var_matrix_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { - if(!is_constant::value && !is_constant::value) { + if (!is_constant::value && !is_constant::value) { arena_t> arena_a = a; arena_t> arena_b = b; - return make_callback_var((arena_a.val().array() * arena_b.val().array().log()).matrix(), - [arena_a, arena_b](const auto& res) mutable { - arena_a.adj().array() += res.adj().array() * arena_b.val().array().log(); - arena_b.adj().array() += res.adj().array() * arena_a.val().array() / arena_b.val().array(); - }); - } else if(!is_constant::value) { + return make_callback_var( + (arena_a.val().array() * arena_b.val().array().log()).matrix(), + [arena_a, arena_b](const auto& res) mutable { + arena_a.adj().array() + += res.adj().array() * arena_b.val().array().log(); + arena_b.adj().array() += res.adj().array() * arena_a.val().array() + / arena_b.val().array(); + }); + } else if (!is_constant::value) { arena_t> arena_a = a; auto arena_b_log = to_arena(value_of(b).array().log()); return make_callback_var((arena_a.val().array() * arena_b_log).matrix(), - [arena_a, arena_b_log](const auto& res) mutable { - arena_a.adj().array() += res.adj().array() * arena_b_log; - }); + [arena_a, arena_b_log](const auto& res) mutable { + arena_a.adj().array() + += res.adj().array() * arena_b_log; + }); } else { auto arena_a = to_arena(value_of(a)); arena_t> arena_b = b; - return make_callback_var((arena_a.array() * arena_b.val().array().log()).matrix(), - [arena_a, arena_b](const auto& res) mutable { - arena_b.adj().array() += res.adj().array() * arena_a.array() / arena_b.val().array(); - }); + return make_callback_var( + (arena_a.array() * arena_b.val().array().log()).matrix(), + [arena_a, arena_b](const auto& res) mutable { + arena_b.adj().array() + += res.adj().array() * arena_a.array() / arena_b.val().array(); + }); } } -template * = nullptr, +template * = nullptr, require_stan_scalar_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { using std::log; - - if(!is_constant::value && !is_constant::value) { + + if (!is_constant::value && !is_constant::value) { arena_t> arena_a = a; var arena_b = b; - return make_callback_var(arena_a.val() * log(arena_b.val()), - [arena_a, arena_b](const auto& res) mutable { - arena_a.adj() += res.adj() * log(arena_b.val()); - arena_b.adj() += (res.adj().array() * arena_a.val().array() / arena_b.val()).sum(); - }); - } else if(!is_constant::value) { + return make_callback_var( + arena_a.val() * log(arena_b.val()), + [arena_a, arena_b](const auto& res) mutable { + arena_a.adj() += res.adj() * log(arena_b.val()); + arena_b.adj() + += (res.adj().array() * arena_a.val().array() / arena_b.val()) + .sum(); + }); + } else if (!is_constant::value) { arena_t> arena_a = a; return make_callback_var(arena_a.val() * log(value_of(b)), - [arena_a, b](const auto& res) mutable { - arena_a.adj() += res.adj() * log(value_of(b)); - }); + [arena_a, b](const auto& res) mutable { + arena_a.adj() += res.adj() * log(value_of(b)); + }); } else { arena_t> arena_a = value_of(a); var arena_b = b; - return make_callback_var(arena_a * log(arena_b.val()), - [arena_a, arena_b](const auto& res) mutable { - arena_b.adj() += (res.adj().array() * arena_a.array() / arena_b.val()).sum(); - }); + return make_callback_var( + arena_a * log(arena_b.val()), + [arena_a, arena_b](const auto& res) mutable { + arena_b.adj() + += (res.adj().array() * arena_a.array() / arena_b.val()).sum(); + }); } } -template * = nullptr, - require_var_matrix_t* = nullptr> +template * = nullptr, + require_var_matrix_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { - if(!is_constant::value && !is_constant::value) { + if (!is_constant::value && !is_constant::value) { var arena_a = a; arena_t> arena_b = b; - return make_callback_var((arena_a.val() * arena_b.val().array().log()).matrix(), - [arena_a, arena_b](const auto& res) mutable { - arena_a.adj() += (res.adj().array() * arena_b.val().array().log()).sum(); - arena_b.adj().array() += res.adj().array() * arena_a.val() / arena_b.val().array(); - }); - } else if(!is_constant::value) { + return make_callback_var( + (arena_a.val() * arena_b.val().array().log()).matrix(), + [arena_a, arena_b](const auto& res) mutable { + arena_a.adj() + += (res.adj().array() * arena_b.val().array().log()).sum(); + arena_b.adj().array() + += res.adj().array() * arena_a.val() / arena_b.val().array(); + }); + } else if (!is_constant::value) { var arena_a = a; auto arena_b_log = to_arena(value_of(b).array().log()); return make_callback_var((arena_a.val() * arena_b_log).matrix(), - [arena_a, arena_b_log](const auto& res) mutable { - arena_a.adj() += (res.adj().array() * arena_b_log).sum(); - }); + [arena_a, arena_b_log](const auto& res) mutable { + arena_a.adj() + += (res.adj().array() * arena_b_log).sum(); + }); } else { arena_t> arena_b = b; - return make_callback_var((value_of(a) * arena_b.val().array().log()).matrix(), - [a, arena_b](const auto& res) mutable { - arena_b.adj().array() += res.adj().array() * value_of(a) / arena_b.val().array(); - }); + return make_callback_var( + (value_of(a) * arena_b.val().array().log()).matrix(), + [a, arena_b](const auto& res) mutable { + arena_b.adj().array() + += res.adj().array() * value_of(a) / arena_b.val().array(); + }); } } From e0b9d53eebe7c5289e1ffba4aa57650791e2e129 Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 20 Dec 2020 16:48:25 -0500 Subject: [PATCH 21/30] Fixed headers test and docs build (Issue #2101) --- stan/math/prim/fun/tanh.hpp | 1 - stan/math/rev/fun/multiply_lower_tri_self_transpose.hpp | 1 + stan/math/rev/fun/read_corr_matrix.hpp | 1 + stan/math/rev/fun/read_cov_matrix.hpp | 7 ++++--- stan/math/rev/fun/rows_dot_product.hpp | 1 + stan/math/rev/fun/tanh.hpp | 9 +++++---- 6 files changed, 12 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/fun/tanh.hpp b/stan/math/prim/fun/tanh.hpp index 2da63ff9014..884760e377f 100644 --- a/stan/math/prim/fun/tanh.hpp +++ b/stan/math/prim/fun/tanh.hpp @@ -70,7 +70,6 @@ namespace internal { template inline std::complex complex_tanh(const std::complex& z) { using std::exp; - using stan::math::operator/; auto exp_z = exp(z); auto exp_neg_z = exp(-z); return stan::math::internal::complex_divide(exp_z - exp_neg_z, diff --git a/stan/math/rev/fun/multiply_lower_tri_self_transpose.hpp b/stan/math/rev/fun/multiply_lower_tri_self_transpose.hpp index 5357ac1f173..5e4d6c42c42 100644 --- a/stan/math/rev/fun/multiply_lower_tri_self_transpose.hpp +++ b/stan/math/rev/fun/multiply_lower_tri_self_transpose.hpp @@ -8,6 +8,7 @@ #include #include #include +#include namespace stan { namespace math { diff --git a/stan/math/rev/fun/read_corr_matrix.hpp b/stan/math/rev/fun/read_corr_matrix.hpp index 61ab3a05e65..a5d8281740d 100644 --- a/stan/math/rev/fun/read_corr_matrix.hpp +++ b/stan/math/rev/fun/read_corr_matrix.hpp @@ -1,6 +1,7 @@ #ifndef STAN_MATH_REV_FUN_READ_CORR_MATRIX_HPP #define STAN_MATH_REV_FUN_READ_CORR_MATRIX_HPP +#include #include #include #include diff --git a/stan/math/rev/fun/read_cov_matrix.hpp b/stan/math/rev/fun/read_cov_matrix.hpp index 2bc946f3ba1..ceeee85f5d8 100644 --- a/stan/math/rev/fun/read_cov_matrix.hpp +++ b/stan/math/rev/fun/read_cov_matrix.hpp @@ -2,9 +2,10 @@ #define STAN_MATH_REV_FUN_READ_COV_MATRIX_HPP #include -#include -#include -#include +#include +#include +#include +#include #include namespace stan { diff --git a/stan/math/rev/fun/rows_dot_product.hpp b/stan/math/rev/fun/rows_dot_product.hpp index 0f9c834e21c..61a014c8356 100644 --- a/stan/math/rev/fun/rows_dot_product.hpp +++ b/stan/math/rev/fun/rows_dot_product.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include namespace stan { diff --git a/stan/math/rev/fun/tanh.hpp b/stan/math/rev/fun/tanh.hpp index 488b5dfd4c6..39eb1bc9579 100644 --- a/stan/math/rev/fun/tanh.hpp +++ b/stan/math/rev/fun/tanh.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -46,11 +47,11 @@ inline var tanh(const var& a) { } /** - * Return the hyperbolic tangent of of x + * Return the hyperbolic tangent of elements of a * - * @tparam T type of x - * @param x argument - * @return elementwise hyperbolic tangent of x + * @tparam T type of a + * @param a argument + * @return elementwise hyperbolic tangent of a */ template * = nullptr> inline auto tanh(const VarMat& a) { From 2754c66604acc44d50d247917ea9bf02bd882b39 Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 20 Dec 2020 16:55:36 -0500 Subject: [PATCH 22/30] Updated `multiply_log` docs (Issue #2101) --- stan/math/rev/fun/multiply_log.hpp | 36 ++++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/stan/math/rev/fun/multiply_log.hpp b/stan/math/rev/fun/multiply_log.hpp index b0e08896426..aa7126be9c8 100644 --- a/stan/math/rev/fun/multiply_log.hpp +++ b/stan/math/rev/fun/multiply_log.hpp @@ -56,8 +56,7 @@ class multiply_log_dv_vari : public op_dv_vari { * * When both a and b are 0, the value returned is 0. * The partial derivative with respect to a is log(b). - * The partial derivative with respect to b is a/b. When - * a and b are both 0, this is set to Inf. + * The partial derivative with respect to b is a/b. * * @param a First variable. * @param b Second variable. @@ -83,8 +82,7 @@ inline var multiply_log(const var& a, double b) { * Return the value of a*log(b). * * When both a and b are 0, the value returned is 0. - * The partial derivative with respect to b is a/b. When - * a and b are both 0, this is set to Inf. + * The partial derivative with respect to b is a/b. * * @param a First scalar. * @param b Second variable. @@ -97,6 +95,18 @@ inline var multiply_log(double a, const var& b) { return var(new internal::multiply_log_dv_vari(a, b.vi_)); } +/** + * Return the elementwise product `a * log(b)`. + * + * Both `T1` and `T2` are matrices, and one of `T1` or `T2` must be a + * `var_value` + * + * @tparam T1 Type of first argument + * @tparam T2 Type of second argument + * @param a First argument + * @param b Second argument + * @return elementwise product of `a` and `log(b)` + */ template * = nullptr, require_any_var_matrix_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { @@ -134,6 +144,15 @@ inline auto multiply_log(const T1& a, const T2& b) { } } +/** + * Return the product `a * log(b)`. + * + * @tparam T1 Type of matrix argument + * @tparam T2 Type of scalar argument + * @param a Matrix argument + * @param b Scalar argument + * @return Product of `a` and `log(b)` + */ template * = nullptr, require_stan_scalar_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { @@ -171,6 +190,15 @@ inline auto multiply_log(const T1& a, const T2& b) { } } +/** + * Return the product `a * log(b)`. + * + * @tparam T1 Type of scalar argument + * @tparam T2 Type of matrix argument + * @param a Scalar argument + * @param b Matrix argument + * @return Product of `a` and `log(b)` + */ template * = nullptr, require_var_matrix_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { From bc40d8d307217ff5b67f3c7fe8df7a3e0838c410 Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 20 Dec 2020 20:49:20 -0500 Subject: [PATCH 23/30] Added checks back to `multiply_log` (Issue #2101) --- stan/math/rev/fun/multiply_log.hpp | 191 +++++++++++++++++++++-------- 1 file changed, 142 insertions(+), 49 deletions(-) diff --git a/stan/math/rev/fun/multiply_log.hpp b/stan/math/rev/fun/multiply_log.hpp index aa7126be9c8..78221236bd3 100644 --- a/stan/math/rev/fun/multiply_log.hpp +++ b/stan/math/rev/fun/multiply_log.hpp @@ -26,7 +26,11 @@ class multiply_log_vv_vari : public op_vv_vari { bvi_->adj_ = NOT_A_NUMBER; } else { avi_->adj_ += adj_ * log(bvi_->val_); - bvi_->adj_ += adj_ * avi_->val_ / bvi_->val_; + if (bvi_->val_ == 0.0 && avi_->val_ == 0) { + bvi_->adj_ += adj_ * INFTY; + } else { + bvi_->adj_ += adj_ * avi_->val_ / bvi_->val_; + } } } }; @@ -47,7 +51,13 @@ class multiply_log_dv_vari : public op_dv_vari { public: multiply_log_dv_vari(double a, vari* bvi) : op_dv_vari(multiply_log(a, bvi->val_), a, bvi) {} - void chain() { bvi_->adj_ += adj_ * ad_ / bvi_->val_; } + void chain() { + if (bvi_->val_ == 0.0 && ad_ == 0.0) { + bvi_->adj_ += adj_ * INFTY; + } else { + bvi_->adj_ += adj_ * ad_ / bvi_->val_; + } + } }; } // namespace internal @@ -110,37 +120,65 @@ inline var multiply_log(double a, const var& b) { template * = nullptr, require_any_var_matrix_t* = nullptr> inline auto multiply_log(const T1& a, const T2& b) { + check_matching_dims("multiply_log", "a", a, "b", b); if (!is_constant::value && !is_constant::value) { arena_t> arena_a = a; arena_t> arena_b = b; - return make_callback_var( - (arena_a.val().array() * arena_b.val().array().log()).matrix(), - [arena_a, arena_b](const auto& res) mutable { - arena_a.adj().array() - += res.adj().array() * arena_b.val().array().log(); - arena_b.adj().array() += res.adj().array() * arena_a.val().array() - / arena_b.val().array(); - }); + return make_callback_var(multiply_log(arena_a.val(), arena_b.val()), + [arena_a, arena_b](const auto& res) mutable { + for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val().coeff(i, j), arena_b.val().coeff(i, j)))) { + arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; + arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + arena_a.adj().coeffRef(i, j) += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); + if (arena_b.val().coeff(i, j) == 0.0 && arena_a.val().coeff(i, j) == 0) { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * arena_a.val().coeff(i, j) / arena_b.val().coeff(i, j); + } + } + } + } + }); } else if (!is_constant::value) { arena_t> arena_a = a; - auto arena_b_log = to_arena(value_of(b).array().log()); + auto arena_b = to_arena(value_of(b)); - return make_callback_var((arena_a.val().array() * arena_b_log).matrix(), - [arena_a, arena_b_log](const auto& res) mutable { - arena_a.adj().array() - += res.adj().array() * arena_b_log; + return make_callback_var(multiply_log(arena_a.val(), arena_b), + [arena_a, arena_b](const auto& res) mutable { + for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val().coeff(i, j), arena_b.coeff(i, j)))) { + arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + arena_a.adj().coeffRef(i, j) += res.adj().coeff(i, j) * log(arena_b.coeff(i, j)); + } + } + } }); } else { auto arena_a = to_arena(value_of(a)); arena_t> arena_b = b; - return make_callback_var( - (arena_a.array() * arena_b.val().array().log()).matrix(), - [arena_a, arena_b](const auto& res) mutable { - arena_b.adj().array() - += res.adj().array() * arena_a.array() / arena_b.val().array(); - }); + return make_callback_var(multiply_log(arena_a, arena_b.val()), + [arena_a, arena_b](const auto& res) mutable { + for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val().coeff(i, j), arena_b.val().coeff(i, j)))) { + arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + if (arena_b.val().coeff(i, j) == 0.0 && arena_a.val().coeff(i, j) == 0) { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * arena_a.coeff(i, j) / arena_b.val().coeff(i, j); + } + } + } + } + }); } } @@ -162,30 +200,58 @@ inline auto multiply_log(const T1& a, const T2& b) { arena_t> arena_a = a; var arena_b = b; - return make_callback_var( - arena_a.val() * log(arena_b.val()), + return make_callback_var(multiply_log(arena_a.val(), arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - arena_a.adj() += res.adj() * log(arena_b.val()); - arena_b.adj() - += (res.adj().array() * arena_a.val().array() / arena_b.val()) - .sum(); + for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val().coeff(i, j), arena_b.val()))) { + arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; + arena_b.adj() = NOT_A_NUMBER; + } else { + arena_a.adj().coeffRef(i, j) += res.adj().coeff(i, j) * log(arena_b.val()); + if (arena_b.val() == 0.0 && arena_a.val().coeff(i, j) == 0) { + arena_b.adj() += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj() += res.adj().coeff(i, j) * arena_a.val().coeff(i, j) / arena_b.val(); + } + } + } + } }); } else if (!is_constant::value) { arena_t> arena_a = a; - return make_callback_var(arena_a.val() * log(value_of(b)), + return make_callback_var(multiply_log(arena_a.val(), value_of(b)), [arena_a, b](const auto& res) mutable { - arena_a.adj() += res.adj() * log(value_of(b)); + for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val().coeff(i, j), value_of(b)))) { + arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + arena_a.adj().coeffRef(i, j) += res.adj().coeff(i, j) * log(value_of(b)); + } + } + } }); } else { arena_t> arena_a = value_of(a); var arena_b = b; - return make_callback_var( - arena_a * log(arena_b.val()), + return make_callback_var(multiply_log(arena_a, arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - arena_b.adj() - += (res.adj().array() * arena_a.array() / arena_b.val()).sum(); + for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val().coeff(i, j), arena_b.val()))) { + arena_b.adj() = NOT_A_NUMBER; + } else { + if (arena_b.val() == 0.0 && arena_a.val().coeff(i, j) == 0) { + arena_b.adj() += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj() += res.adj().coeff(i, j) * arena_a.val().coeff(i, j) / arena_b.val(); + } + } + } + } }); } } @@ -206,31 +272,58 @@ inline auto multiply_log(const T1& a, const T2& b) { var arena_a = a; arena_t> arena_b = b; - return make_callback_var( - (arena_a.val() * arena_b.val().array().log()).matrix(), + return make_callback_var(multiply_log(arena_a.val(), arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - arena_a.adj() - += (res.adj().array() * arena_b.val().array().log()).sum(); - arena_b.adj().array() - += res.adj().array() * arena_a.val() / arena_b.val().array(); + for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val(), arena_b.val().coeff(i, j)))) { + arena_a.adj() = NOT_A_NUMBER; + arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + arena_a.adj() += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); + if (arena_b.val().coeff(i, j) == 0.0 && arena_a.val() == 0) { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * arena_a.val() / arena_b.val().coeff(i, j); + } + } + } + } }); } else if (!is_constant::value) { var arena_a = a; - auto arena_b_log = to_arena(value_of(b).array().log()); + auto arena_b = to_arena(value_of(b)); - return make_callback_var((arena_a.val() * arena_b_log).matrix(), - [arena_a, arena_b_log](const auto& res) mutable { - arena_a.adj() - += (res.adj().array() * arena_b_log).sum(); - }); + return make_callback_var(multiply_log(arena_a.val(), arena_b), + [arena_a, arena_b](const auto& res) mutable { + for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val(), arena_b.val().coeff(i, j)))) { + arena_a.adj() = NOT_A_NUMBER; + } else { + arena_a.adj() += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); + } + } + } + }); } else { arena_t> arena_b = b; - return make_callback_var( - (value_of(a) * arena_b.val().array().log()).matrix(), + return make_callback_var(multiply_log(value_of(a), arena_b.val()), [a, arena_b](const auto& res) mutable { - arena_b.adj().array() - += res.adj().array() * value_of(a) / arena_b.val().array(); + for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(value_of(a), arena_b.val().coeff(i, j)))) { + arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + if (arena_b.val().coeff(i, j) == 0.0 && value_of(a) == 0) { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * value_of(a) / arena_b.val().coeff(i, j); + } + } + } + } }); } } From 014e515bd69f6c4f387d352f4048251b1e95c5eb Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 21 Dec 2020 01:51:10 +0000 Subject: [PATCH 24/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/multiply_log.hpp | 296 ++++++++++++++++------------- 1 file changed, 166 insertions(+), 130 deletions(-) diff --git a/stan/math/rev/fun/multiply_log.hpp b/stan/math/rev/fun/multiply_log.hpp index 78221236bd3..ea691f2ad64 100644 --- a/stan/math/rev/fun/multiply_log.hpp +++ b/stan/math/rev/fun/multiply_log.hpp @@ -125,60 +125,74 @@ inline auto multiply_log(const T1& a, const T2& b) { arena_t> arena_a = a; arena_t> arena_b = b; - return make_callback_var(multiply_log(arena_a.val(), arena_b.val()), - [arena_a, arena_b](const auto& res) mutable { - for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val().coeff(i, j), arena_b.val().coeff(i, j)))) { - arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; - arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - arena_a.adj().coeffRef(i, j) += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); - if (arena_b.val().coeff(i, j) == 0.0 && arena_a.val().coeff(i, j) == 0) { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * arena_a.val().coeff(i, j) / arena_b.val().coeff(i, j); - } - } - } - } - }); + return make_callback_var( + multiply_log(arena_a.val(), arena_b.val()), + [arena_a, arena_b](const auto& res) mutable { + for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val().coeff(i, j), + arena_b.val().coeff(i, j)))) { + arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; + arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + arena_a.adj().coeffRef(i, j) + += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); + if (arena_b.val().coeff(i, j) == 0.0 + && arena_a.val().coeff(i, j) == 0) { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) + * arena_a.val().coeff(i, j) + / arena_b.val().coeff(i, j); + } + } + } + } + }); } else if (!is_constant::value) { arena_t> arena_a = a; auto arena_b = to_arena(value_of(b)); - return make_callback_var(multiply_log(arena_a.val(), arena_b), - [arena_a, arena_b](const auto& res) mutable { - for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val().coeff(i, j), arena_b.coeff(i, j)))) { - arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - arena_a.adj().coeffRef(i, j) += res.adj().coeff(i, j) * log(arena_b.coeff(i, j)); - } - } - } - }); + return make_callback_var( + multiply_log(arena_a.val(), arena_b), + [arena_a, arena_b](const auto& res) mutable { + for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val().coeff(i, j), + arena_b.coeff(i, j)))) { + arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + arena_a.adj().coeffRef(i, j) + += res.adj().coeff(i, j) * log(arena_b.coeff(i, j)); + } + } + } + }); } else { auto arena_a = to_arena(value_of(a)); arena_t> arena_b = b; - return make_callback_var(multiply_log(arena_a, arena_b.val()), - [arena_a, arena_b](const auto& res) mutable { - for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val().coeff(i, j), arena_b.val().coeff(i, j)))) { - arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - if (arena_b.val().coeff(i, j) == 0.0 && arena_a.val().coeff(i, j) == 0) { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * arena_a.coeff(i, j) / arena_b.val().coeff(i, j); - } - } - } - } - }); + return make_callback_var( + multiply_log(arena_a, arena_b.val()), + [arena_a, arena_b](const auto& res) mutable { + for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely(is_any_nan(arena_a.val().coeff(i, j), + arena_b.val().coeff(i, j)))) { + arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + if (arena_b.val().coeff(i, j) == 0.0 + && arena_a.val().coeff(i, j) == 0) { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) + * arena_a.coeff(i, j) + / arena_b.val().coeff(i, j); + } + } + } + } + }); } } @@ -200,58 +214,68 @@ inline auto multiply_log(const T1& a, const T2& b) { arena_t> arena_a = a; var arena_b = b; - return make_callback_var(multiply_log(arena_a.val(), arena_b.val()), + return make_callback_var( + multiply_log(arena_a.val(), arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val().coeff(i, j), arena_b.val()))) { - arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; - arena_b.adj() = NOT_A_NUMBER; - } else { - arena_a.adj().coeffRef(i, j) += res.adj().coeff(i, j) * log(arena_b.val()); - if (arena_b.val() == 0.0 && arena_a.val().coeff(i, j) == 0) { - arena_b.adj() += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj() += res.adj().coeff(i, j) * arena_a.val().coeff(i, j) / arena_b.val(); - } - } - } - } + for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely( + is_any_nan(arena_a.val().coeff(i, j), arena_b.val()))) { + arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; + arena_b.adj() = NOT_A_NUMBER; + } else { + arena_a.adj().coeffRef(i, j) + += res.adj().coeff(i, j) * log(arena_b.val()); + if (arena_b.val() == 0.0 && arena_a.val().coeff(i, j) == 0) { + arena_b.adj() += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj() += res.adj().coeff(i, j) + * arena_a.val().coeff(i, j) / arena_b.val(); + } + } + } + } }); } else if (!is_constant::value) { arena_t> arena_a = a; - return make_callback_var(multiply_log(arena_a.val(), value_of(b)), - [arena_a, b](const auto& res) mutable { - for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val().coeff(i, j), value_of(b)))) { - arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - arena_a.adj().coeffRef(i, j) += res.adj().coeff(i, j) * log(value_of(b)); - } - } - } - }); + return make_callback_var( + multiply_log(arena_a.val(), value_of(b)), + [arena_a, b](const auto& res) mutable { + for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely( + is_any_nan(arena_a.val().coeff(i, j), value_of(b)))) { + arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + arena_a.adj().coeffRef(i, j) + += res.adj().coeff(i, j) * log(value_of(b)); + } + } + } + }); } else { arena_t> arena_a = value_of(a); var arena_b = b; - return make_callback_var(multiply_log(arena_a, arena_b.val()), + return make_callback_var( + multiply_log(arena_a, arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val().coeff(i, j), arena_b.val()))) { - arena_b.adj() = NOT_A_NUMBER; - } else { - if (arena_b.val() == 0.0 && arena_a.val().coeff(i, j) == 0) { - arena_b.adj() += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj() += res.adj().coeff(i, j) * arena_a.val().coeff(i, j) / arena_b.val(); - } - } - } - } + for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely( + is_any_nan(arena_a.val().coeff(i, j), arena_b.val()))) { + arena_b.adj() = NOT_A_NUMBER; + } else { + if (arena_b.val() == 0.0 && arena_a.val().coeff(i, j) == 0) { + arena_b.adj() += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj() += res.adj().coeff(i, j) + * arena_a.val().coeff(i, j) / arena_b.val(); + } + } + } + } }); } } @@ -272,58 +296,70 @@ inline auto multiply_log(const T1& a, const T2& b) { var arena_a = a; arena_t> arena_b = b; - return make_callback_var(multiply_log(arena_a.val(), arena_b.val()), + return make_callback_var( + multiply_log(arena_a.val(), arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val(), arena_b.val().coeff(i, j)))) { - arena_a.adj() = NOT_A_NUMBER; - arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - arena_a.adj() += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); - if (arena_b.val().coeff(i, j) == 0.0 && arena_a.val() == 0) { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * arena_a.val() / arena_b.val().coeff(i, j); - } - } - } - } + for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely( + is_any_nan(arena_a.val(), arena_b.val().coeff(i, j)))) { + arena_a.adj() = NOT_A_NUMBER; + arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + arena_a.adj() + += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); + if (arena_b.val().coeff(i, j) == 0.0 && arena_a.val() == 0) { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) + * arena_a.val() + / arena_b.val().coeff(i, j); + } + } + } + } }); } else if (!is_constant::value) { var arena_a = a; auto arena_b = to_arena(value_of(b)); - return make_callback_var(multiply_log(arena_a.val(), arena_b), - [arena_a, arena_b](const auto& res) mutable { - for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val(), arena_b.val().coeff(i, j)))) { - arena_a.adj() = NOT_A_NUMBER; - } else { - arena_a.adj() += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); - } - } - } - }); + return make_callback_var( + multiply_log(arena_a.val(), arena_b), + [arena_a, arena_b](const auto& res) mutable { + for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely( + is_any_nan(arena_a.val(), arena_b.val().coeff(i, j)))) { + arena_a.adj() = NOT_A_NUMBER; + } else { + arena_a.adj() + += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); + } + } + } + }); } else { arena_t> arena_b = b; - return make_callback_var(multiply_log(value_of(a), arena_b.val()), + return make_callback_var( + multiply_log(value_of(a), arena_b.val()), [a, arena_b](const auto& res) mutable { - for(Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for(Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(value_of(a), arena_b.val().coeff(i, j)))) { - arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - if (arena_b.val().coeff(i, j) == 0.0 && value_of(a) == 0) { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * value_of(a) / arena_b.val().coeff(i, j); - } - } - } - } + for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { + for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { + if (unlikely( + is_any_nan(value_of(a), arena_b.val().coeff(i, j)))) { + arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; + } else { + if (arena_b.val().coeff(i, j) == 0.0 && value_of(a) == 0) { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; + } else { + arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) + * value_of(a) + / arena_b.val().coeff(i, j); + } + } + } + } }); } } From 37e5252d8870ae238695c0d4d056269f68544590 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 30 Dec 2020 14:50:45 -0500 Subject: [PATCH 25/30] Don't use noalias and do use make callback var (Issue #2101) --- stan/math/rev/fun/exp.hpp | 9 +++------ stan/math/rev/fun/inv_logit.hpp | 15 ++++++--------- stan/math/rev/fun/log.hpp | 11 ++++------- stan/math/rev/fun/log1m.hpp | 11 ++++------- stan/math/rev/fun/square.hpp | 11 ++++------- 5 files changed, 21 insertions(+), 36 deletions(-) diff --git a/stan/math/rev/fun/exp.hpp b/stan/math/rev/fun/exp.hpp index c08477f9567..e8799327e6c 100644 --- a/stan/math/rev/fun/exp.hpp +++ b/stan/math/rev/fun/exp.hpp @@ -66,13 +66,10 @@ inline std::complex exp(const std::complex& z) { */ template * = nullptr> inline auto exp(const T& x) { - plain_type_t res = x.val().array().exp().matrix(); - - reverse_pass_callback([x, res]() mutable { - x.adj().noalias() += (res.val().array() * res.adj().array()).matrix(); + return make_callback_var(x.val().array().exp().matrix(), + [x](const auto& vi) mutable { + x.adj() += (vi.val().array() * vi.adj().array()).matrix(); }); - - return res; } } // namespace math diff --git a/stan/math/rev/fun/inv_logit.hpp b/stan/math/rev/fun/inv_logit.hpp index 03bbd635938..cbf9109116a 100644 --- a/stan/math/rev/fun/inv_logit.hpp +++ b/stan/math/rev/fun/inv_logit.hpp @@ -42,15 +42,12 @@ inline var inv_logit(const var& a) { */ template * = nullptr> inline auto inv_logit(const T& x) { - plain_type_t res = inv_logit(x.val()); - - reverse_pass_callback([x, res]() mutable { - x.adj().noalias() - += (res.adj().array() * res.val().array() * (1.0 - res.val().array())) - .matrix(); - }); - - return res; + return make_callback_var(inv_logit(x.val()), + [x](const auto& vi) mutable { + x.adj() + += (vi.adj().array() * vi.val().array() * (1.0 - vi.val().array())) + .matrix(); + }); } } // namespace math diff --git a/stan/math/rev/fun/log.hpp b/stan/math/rev/fun/log.hpp index d23aa017c39..aaa8e3316ca 100644 --- a/stan/math/rev/fun/log.hpp +++ b/stan/math/rev/fun/log.hpp @@ -73,13 +73,10 @@ inline std::complex log(const std::complex& z) { */ template * = nullptr> inline auto log(const T& x) { - plain_type_t res = x.val().array().log().matrix(); - - reverse_pass_callback([x, res]() mutable { - x.adj().noalias() += (res.adj().array() / x.val().array()).matrix(); - }); - - return res; + return make_callback_var(x.val().array().log().matrix(), + [x](const auto& vi) mutable { + x.adj() += (vi.adj().array() / x.val().array()).matrix(); + }); } } // namespace math diff --git a/stan/math/rev/fun/log1m.hpp b/stan/math/rev/fun/log1m.hpp index 8e9523bd933..ee63c536fcf 100644 --- a/stan/math/rev/fun/log1m.hpp +++ b/stan/math/rev/fun/log1m.hpp @@ -37,13 +37,10 @@ inline var log1m(const var& a) { return var(new internal::log1m_vari(a.vi_)); } */ template * = nullptr> inline auto log1m(const T& x) { - T res = stan::math::log1m(x.val()); - - reverse_pass_callback([x, res]() mutable { - x.adj().noalias() += (res.adj().array() / (x.val().array() - 1.0)).matrix(); - }); - - return res; + return make_callback_var(stan::math::log1m(x.val()), + [x](const auto& vi) mutable { + x.adj() += (vi.adj().array() / (x.val().array() - 1.0)).matrix(); + }); } } // namespace math diff --git a/stan/math/rev/fun/square.hpp b/stan/math/rev/fun/square.hpp index 0ce53d01654..b33f51e1218 100644 --- a/stan/math/rev/fun/square.hpp +++ b/stan/math/rev/fun/square.hpp @@ -53,13 +53,10 @@ inline var square(const var& x) { */ template * = nullptr> inline auto square(const T& x) { - T res = (x.val().array() * x.val().array()).matrix(); - - reverse_pass_callback([x, res]() mutable { - x.adj().noalias() += (2.0 * x.val().array() * res.adj().array()).matrix(); - }); - - return res; + return make_callback_var((x.val().array() * x.val().array()).matrix(), + [x](const auto& vi) mutable { + x.adj() += (2.0 * x.val().array() * vi.adj().array()).matrix(); + }); } } // namespace math From 46d590b744a224dccff34b2d6461592e84fd0a57 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 30 Dec 2020 15:09:07 -0500 Subject: [PATCH 26/30] Cleaned up `multiply_log` (Issue #2101) --- stan/math/prim/fun/multiply_log.hpp | 12 +- stan/math/rev/fun/multiply_log.hpp | 195 +++++----------------------- 2 files changed, 37 insertions(+), 170 deletions(-) diff --git a/stan/math/prim/fun/multiply_log.hpp b/stan/math/prim/fun/multiply_log.hpp index 731f7ddbbba..63f1d792e97 100644 --- a/stan/math/prim/fun/multiply_log.hpp +++ b/stan/math/prim/fun/multiply_log.hpp @@ -20,26 +20,21 @@ namespace math { \mbox{multiply\_log}(x, y) = \begin{cases} 0 & \mbox{if } x=y=0\\ - x\ln y & \mbox{if } x, y\neq0 \\[6pt] - \textrm{NaN} & \mbox{if } x = \textrm{NaN or } y = \textrm{NaN} + x\ln y & \mbox{if } x, y\neq 0 \\[6pt] \end{cases} \f] \f[ \frac{\partial\, \mbox{multiply\_log}(x, y)}{\partial x} = \begin{cases} - \infty & \mbox{if } x=y=0\\ - \ln y & \mbox{if } x, y\neq 0 \\[6pt] - \textrm{NaN} & \mbox{if } x = \textrm{NaN or } y = \textrm{NaN} + \ln y \\[6pt] \end{cases} \f] \f[ \frac{\partial\, \mbox{multiply\_log}(x, y)}{\partial y} = \begin{cases} - \infty & \mbox{if } x=y=0\\ - \frac{x}{y} & \mbox{if } x, y\neq 0 \\[6pt] - \textrm{NaN} & \mbox{if } x = \textrm{NaN or } y = \textrm{NaN} + \frac{x}{y} \\[6pt] \end{cases} \f] * @@ -56,6 +51,7 @@ inline return_type_t multiply_log(const T_a a, const T_b b) { if (b == 0.0 && a == 0.0) { return 0.0; } + return a * log(b); } diff --git a/stan/math/rev/fun/multiply_log.hpp b/stan/math/rev/fun/multiply_log.hpp index ea691f2ad64..cbbd4cdcdab 100644 --- a/stan/math/rev/fun/multiply_log.hpp +++ b/stan/math/rev/fun/multiply_log.hpp @@ -21,17 +21,8 @@ class multiply_log_vv_vari : public op_vv_vari { : op_vv_vari(multiply_log(avi->val_, bvi->val_), avi, bvi) {} void chain() { using std::log; - if (unlikely(is_any_nan(avi_->val_, bvi_->val_))) { - avi_->adj_ = NOT_A_NUMBER; - bvi_->adj_ = NOT_A_NUMBER; - } else { - avi_->adj_ += adj_ * log(bvi_->val_); - if (bvi_->val_ == 0.0 && avi_->val_ == 0) { - bvi_->adj_ += adj_ * INFTY; - } else { - bvi_->adj_ += adj_ * avi_->val_ / bvi_->val_; - } - } + avi_->adj_ += adj_ * log(bvi_->val_); + bvi_->adj_ += adj_ * avi_->val_ / bvi_->val_; } }; class multiply_log_vd_vari : public op_vd_vari { @@ -40,11 +31,7 @@ class multiply_log_vd_vari : public op_vd_vari { : op_vd_vari(multiply_log(avi->val_, b), avi, b) {} void chain() { using std::log; - if (unlikely(is_any_nan(avi_->val_, bd_))) { - avi_->adj_ = NOT_A_NUMBER; - } else { - avi_->adj_ += adj_ * log(bd_); - } + avi_->adj_ += adj_ * log(bd_); } }; class multiply_log_dv_vari : public op_dv_vari { @@ -52,11 +39,7 @@ class multiply_log_dv_vari : public op_dv_vari { multiply_log_dv_vari(double a, vari* bvi) : op_dv_vari(multiply_log(a, bvi->val_), a, bvi) {} void chain() { - if (bvi_->val_ == 0.0 && ad_ == 0.0) { - bvi_->adj_ += adj_ * INFTY; - } else { - bvi_->adj_ += adj_ * ad_ / bvi_->val_; - } + bvi_->adj_ += adj_ * ad_ / bvi_->val_; } }; } // namespace internal @@ -128,70 +111,30 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(arena_a.val(), arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val().coeff(i, j), - arena_b.val().coeff(i, j)))) { - arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; - arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - arena_a.adj().coeffRef(i, j) - += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); - if (arena_b.val().coeff(i, j) == 0.0 - && arena_a.val().coeff(i, j) == 0) { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) - * arena_a.val().coeff(i, j) - / arena_b.val().coeff(i, j); - } - } - } - } + arena_a.adj().array() += res.adj().array() * arena_b.val().array().log(); + arena_b.adj().array() += res.adj().array() + * arena_a.val().array() + / arena_b.val().array(); }); } else if (!is_constant::value) { arena_t> arena_a = a; - auto arena_b = to_arena(value_of(b)); + arena_t> arena_b = value_of(b); return make_callback_var( multiply_log(arena_a.val(), arena_b), [arena_a, arena_b](const auto& res) mutable { - for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val().coeff(i, j), - arena_b.coeff(i, j)))) { - arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - arena_a.adj().coeffRef(i, j) - += res.adj().coeff(i, j) * log(arena_b.coeff(i, j)); - } - } - } + arena_a.adj().array() += res.adj().array() * arena_b.val().array().log(); }); } else { - auto arena_a = to_arena(value_of(a)); + arena_t> arena_a = value_of(a); arena_t> arena_b = b; return make_callback_var( multiply_log(arena_a, arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely(is_any_nan(arena_a.val().coeff(i, j), - arena_b.val().coeff(i, j)))) { - arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - if (arena_b.val().coeff(i, j) == 0.0 - && arena_a.val().coeff(i, j) == 0) { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) - * arena_a.coeff(i, j) - / arena_b.val().coeff(i, j); - } - } - } - } + arena_b.adj().array() += res.adj().array() + * arena_a.val().array() + / arena_b.val().array(); }); } } @@ -217,24 +160,10 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(arena_a.val(), arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely( - is_any_nan(arena_a.val().coeff(i, j), arena_b.val()))) { - arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; - arena_b.adj() = NOT_A_NUMBER; - } else { - arena_a.adj().coeffRef(i, j) - += res.adj().coeff(i, j) * log(arena_b.val()); - if (arena_b.val() == 0.0 && arena_a.val().coeff(i, j) == 0) { - arena_b.adj() += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj() += res.adj().coeff(i, j) - * arena_a.val().coeff(i, j) / arena_b.val(); - } - } - } - } + arena_a.adj().array() + += res.adj().array() * log(arena_b.val()); + arena_b.adj() += (res.adj().array() + * arena_a.val().array()).sum() / arena_b.val(); }); } else if (!is_constant::value) { arena_t> arena_a = a; @@ -242,17 +171,8 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(arena_a.val(), value_of(b)), [arena_a, b](const auto& res) mutable { - for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely( - is_any_nan(arena_a.val().coeff(i, j), value_of(b)))) { - arena_a.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - arena_a.adj().coeffRef(i, j) - += res.adj().coeff(i, j) * log(value_of(b)); - } - } - } + arena_a.adj().array() + += res.adj().array() * log(value_of(b)); }); } else { arena_t> arena_a = value_of(a); @@ -261,21 +181,8 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(arena_a, arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely( - is_any_nan(arena_a.val().coeff(i, j), arena_b.val()))) { - arena_b.adj() = NOT_A_NUMBER; - } else { - if (arena_b.val() == 0.0 && arena_a.val().coeff(i, j) == 0) { - arena_b.adj() += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj() += res.adj().coeff(i, j) - * arena_a.val().coeff(i, j) / arena_b.val(); - } - } - } - } + arena_b.adj() += (res.adj().array() + * arena_a.array()).sum() / arena_b.val(); }); } } @@ -299,44 +206,21 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(arena_a.val(), arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely( - is_any_nan(arena_a.val(), arena_b.val().coeff(i, j)))) { - arena_a.adj() = NOT_A_NUMBER; - arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - arena_a.adj() - += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); - if (arena_b.val().coeff(i, j) == 0.0 && arena_a.val() == 0) { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) - * arena_a.val() - / arena_b.val().coeff(i, j); - } - } - } - } + arena_a.adj() + += (res.adj().array() * arena_b.val().array().log()).sum(); + arena_b.adj().array() += arena_a.val() + * res.adj().array() + / arena_b.val().array(); }); } else if (!is_constant::value) { var arena_a = a; - auto arena_b = to_arena(value_of(b)); + arena_t> arena_b = value_of(b); return make_callback_var( multiply_log(arena_a.val(), arena_b), [arena_a, arena_b](const auto& res) mutable { - for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely( - is_any_nan(arena_a.val(), arena_b.val().coeff(i, j)))) { - arena_a.adj() = NOT_A_NUMBER; - } else { - arena_a.adj() - += res.adj().coeff(i, j) * log(arena_b.val().coeff(i, j)); - } - } - } + arena_a.adj() + += (res.adj().array() * arena_b.val().array().log()).sum(); }); } else { arena_t> arena_b = b; @@ -344,22 +228,9 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(value_of(a), arena_b.val()), [a, arena_b](const auto& res) mutable { - for (Eigen::Index j = 0; j < res.adj().cols(); ++j) { - for (Eigen::Index i = 0; i < res.adj().rows(); ++i) { - if (unlikely( - is_any_nan(value_of(a), arena_b.val().coeff(i, j)))) { - arena_b.adj().coeffRef(i, j) = NOT_A_NUMBER; - } else { - if (arena_b.val().coeff(i, j) == 0.0 && value_of(a) == 0) { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) * INFTY; - } else { - arena_b.adj().coeffRef(i, j) += res.adj().coeff(i, j) - * value_of(a) - / arena_b.val().coeff(i, j); - } - } - } - } + arena_b.adj().array() += value_of(a) + * res.adj().array() + / arena_b.val().array(); }); } } From f293d1c320183bed87526bc9615e5e021bf2dbe6 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 30 Dec 2020 20:15:33 +0000 Subject: [PATCH 27/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/inv_logit.hpp | 6 +-- stan/math/prim/fun/log1m.hpp | 6 +-- stan/math/rev/fun/exp.hpp | 8 +-- stan/math/rev/fun/inv_logit.hpp | 10 ++-- stan/math/rev/fun/log.hpp | 8 +-- stan/math/rev/fun/log1m.hpp | 8 +-- stan/math/rev/fun/multiply_log.hpp | 84 ++++++++++++++---------------- stan/math/rev/fun/square.hpp | 9 ++-- 8 files changed, 66 insertions(+), 73 deletions(-) diff --git a/stan/math/prim/fun/inv_logit.hpp b/stan/math/prim/fun/inv_logit.hpp index f3168b2eeea..7c8355f9ffe 100644 --- a/stan/math/prim/fun/inv_logit.hpp +++ b/stan/math/prim/fun/inv_logit.hpp @@ -81,9 +81,9 @@ struct inv_logit_fun { * @param x container * @return Inverse logit applied to each value in x. */ -template * = nullptr, - require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> +template < + typename T, require_not_var_matrix_t* = nullptr, + require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> inline auto inv_logit(const T& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/prim/fun/log1m.hpp b/stan/math/prim/fun/log1m.hpp index fcd8f744d53..2b81dbd1edb 100644 --- a/stan/math/prim/fun/log1m.hpp +++ b/stan/math/prim/fun/log1m.hpp @@ -67,9 +67,9 @@ struct log1m_fun { * @param x container * @return Natural log of 1 minus each value in x. */ -template * = nullptr, - require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> +template < + typename T, require_not_var_matrix_t* = nullptr, + require_all_not_nonscalar_prim_or_rev_kernel_expression_t* = nullptr> inline auto log1m(const T& x) { return apply_scalar_unary::apply(x); } diff --git a/stan/math/rev/fun/exp.hpp b/stan/math/rev/fun/exp.hpp index e8799327e6c..53389bbb33e 100644 --- a/stan/math/rev/fun/exp.hpp +++ b/stan/math/rev/fun/exp.hpp @@ -66,10 +66,10 @@ inline std::complex exp(const std::complex& z) { */ template * = nullptr> inline auto exp(const T& x) { - return make_callback_var(x.val().array().exp().matrix(), - [x](const auto& vi) mutable { - x.adj() += (vi.val().array() * vi.adj().array()).matrix(); - }); + return make_callback_var( + x.val().array().exp().matrix(), [x](const auto& vi) mutable { + x.adj() += (vi.val().array() * vi.adj().array()).matrix(); + }); } } // namespace math diff --git a/stan/math/rev/fun/inv_logit.hpp b/stan/math/rev/fun/inv_logit.hpp index cbf9109116a..e506178d695 100644 --- a/stan/math/rev/fun/inv_logit.hpp +++ b/stan/math/rev/fun/inv_logit.hpp @@ -42,12 +42,10 @@ inline var inv_logit(const var& a) { */ template * = nullptr> inline auto inv_logit(const T& x) { - return make_callback_var(inv_logit(x.val()), - [x](const auto& vi) mutable { - x.adj() - += (vi.adj().array() * vi.val().array() * (1.0 - vi.val().array())) - .matrix(); - }); + return make_callback_var(inv_logit(x.val()), [x](const auto& vi) mutable { + x.adj() += (vi.adj().array() * vi.val().array() * (1.0 - vi.val().array())) + .matrix(); + }); } } // namespace math diff --git a/stan/math/rev/fun/log.hpp b/stan/math/rev/fun/log.hpp index aaa8e3316ca..6eaf4ed9d0c 100644 --- a/stan/math/rev/fun/log.hpp +++ b/stan/math/rev/fun/log.hpp @@ -73,10 +73,10 @@ inline std::complex log(const std::complex& z) { */ template * = nullptr> inline auto log(const T& x) { - return make_callback_var(x.val().array().log().matrix(), - [x](const auto& vi) mutable { - x.adj() += (vi.adj().array() / x.val().array()).matrix(); - }); + return make_callback_var( + x.val().array().log().matrix(), [x](const auto& vi) mutable { + x.adj() += (vi.adj().array() / x.val().array()).matrix(); + }); } } // namespace math diff --git a/stan/math/rev/fun/log1m.hpp b/stan/math/rev/fun/log1m.hpp index ee63c536fcf..70306cbafe9 100644 --- a/stan/math/rev/fun/log1m.hpp +++ b/stan/math/rev/fun/log1m.hpp @@ -37,10 +37,10 @@ inline var log1m(const var& a) { return var(new internal::log1m_vari(a.vi_)); } */ template * = nullptr> inline auto log1m(const T& x) { - return make_callback_var(stan::math::log1m(x.val()), - [x](const auto& vi) mutable { - x.adj() += (vi.adj().array() / (x.val().array() - 1.0)).matrix(); - }); + return make_callback_var( + stan::math::log1m(x.val()), [x](const auto& vi) mutable { + x.adj() += (vi.adj().array() / (x.val().array() - 1.0)).matrix(); + }); } } // namespace math diff --git a/stan/math/rev/fun/multiply_log.hpp b/stan/math/rev/fun/multiply_log.hpp index cbbd4cdcdab..5b8850d0f66 100644 --- a/stan/math/rev/fun/multiply_log.hpp +++ b/stan/math/rev/fun/multiply_log.hpp @@ -38,9 +38,7 @@ class multiply_log_dv_vari : public op_dv_vari { public: multiply_log_dv_vari(double a, vari* bvi) : op_dv_vari(multiply_log(a, bvi->val_), a, bvi) {} - void chain() { - bvi_->adj_ += adj_ * ad_ / bvi_->val_; - } + void chain() { bvi_->adj_ += adj_ * ad_ / bvi_->val_; } }; } // namespace internal @@ -111,31 +109,31 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(arena_a.val(), arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - arena_a.adj().array() += res.adj().array() * arena_b.val().array().log(); - arena_b.adj().array() += res.adj().array() - * arena_a.val().array() - / arena_b.val().array(); + arena_a.adj().array() + += res.adj().array() * arena_b.val().array().log(); + arena_b.adj().array() += res.adj().array() * arena_a.val().array() + / arena_b.val().array(); }); } else if (!is_constant::value) { arena_t> arena_a = a; arena_t> arena_b = value_of(b); - return make_callback_var( - multiply_log(arena_a.val(), arena_b), - [arena_a, arena_b](const auto& res) mutable { - arena_a.adj().array() += res.adj().array() * arena_b.val().array().log(); - }); + return make_callback_var(multiply_log(arena_a.val(), arena_b), + [arena_a, arena_b](const auto& res) mutable { + arena_a.adj().array() + += res.adj().array() + * arena_b.val().array().log(); + }); } else { arena_t> arena_a = value_of(a); arena_t> arena_b = b; - return make_callback_var( - multiply_log(arena_a, arena_b.val()), - [arena_a, arena_b](const auto& res) mutable { - arena_b.adj().array() += res.adj().array() - * arena_a.val().array() - / arena_b.val().array(); - }); + return make_callback_var(multiply_log(arena_a, arena_b.val()), + [arena_a, arena_b](const auto& res) mutable { + arena_b.adj().array() += res.adj().array() + * arena_a.val().array() + / arena_b.val().array(); + }); } } @@ -160,20 +158,18 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(arena_a.val(), arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - arena_a.adj().array() - += res.adj().array() * log(arena_b.val()); - arena_b.adj() += (res.adj().array() - * arena_a.val().array()).sum() / arena_b.val(); + arena_a.adj().array() += res.adj().array() * log(arena_b.val()); + arena_b.adj() += (res.adj().array() * arena_a.val().array()).sum() + / arena_b.val(); }); } else if (!is_constant::value) { arena_t> arena_a = a; - return make_callback_var( - multiply_log(arena_a.val(), value_of(b)), - [arena_a, b](const auto& res) mutable { - arena_a.adj().array() - += res.adj().array() * log(value_of(b)); - }); + return make_callback_var(multiply_log(arena_a.val(), value_of(b)), + [arena_a, b](const auto& res) mutable { + arena_a.adj().array() + += res.adj().array() * log(value_of(b)); + }); } else { arena_t> arena_a = value_of(a); var arena_b = b; @@ -181,8 +177,8 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(arena_a, arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - arena_b.adj() += (res.adj().array() - * arena_a.array()).sum() / arena_b.val(); + arena_b.adj() + += (res.adj().array() * arena_a.array()).sum() / arena_b.val(); }); } } @@ -206,11 +202,10 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(arena_a.val(), arena_b.val()), [arena_a, arena_b](const auto& res) mutable { - arena_a.adj() - += (res.adj().array() * arena_b.val().array().log()).sum(); - arena_b.adj().array() += arena_a.val() - * res.adj().array() - / arena_b.val().array(); + arena_a.adj() + += (res.adj().array() * arena_b.val().array().log()).sum(); + arena_b.adj().array() + += arena_a.val() * res.adj().array() / arena_b.val().array(); }); } else if (!is_constant::value) { var arena_a = a; @@ -219,19 +214,18 @@ inline auto multiply_log(const T1& a, const T2& b) { return make_callback_var( multiply_log(arena_a.val(), arena_b), [arena_a, arena_b](const auto& res) mutable { - arena_a.adj() - += (res.adj().array() * arena_b.val().array().log()).sum(); + arena_a.adj() + += (res.adj().array() * arena_b.val().array().log()).sum(); }); } else { arena_t> arena_b = b; - return make_callback_var( - multiply_log(value_of(a), arena_b.val()), - [a, arena_b](const auto& res) mutable { - arena_b.adj().array() += value_of(a) - * res.adj().array() - / arena_b.val().array(); - }); + return make_callback_var(multiply_log(value_of(a), arena_b.val()), + [a, arena_b](const auto& res) mutable { + arena_b.adj().array() += value_of(a) + * res.adj().array() + / arena_b.val().array(); + }); } } diff --git a/stan/math/rev/fun/square.hpp b/stan/math/rev/fun/square.hpp index b33f51e1218..b502aa8d807 100644 --- a/stan/math/rev/fun/square.hpp +++ b/stan/math/rev/fun/square.hpp @@ -53,10 +53,11 @@ inline var square(const var& x) { */ template * = nullptr> inline auto square(const T& x) { - return make_callback_var((x.val().array() * x.val().array()).matrix(), - [x](const auto& vi) mutable { - x.adj() += (2.0 * x.val().array() * vi.adj().array()).matrix(); - }); + return make_callback_var( + (x.val().array() * x.val().array()).matrix(), + [x](const auto& vi) mutable { + x.adj() += (2.0 * x.val().array() * vi.adj().array()).matrix(); + }); } } // namespace math From 768a6652749b74170b0845b2a569e151206de617 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 30 Dec 2020 15:50:00 -0500 Subject: [PATCH 28/30] Edits to constrain functions (Issue #2101) --- stan/math/rev/fun/cov_matrix_constrain.hpp | 26 +++++++++++++--------- stan/math/rev/fun/read_corr_L.hpp | 14 ++++++------ stan/math/rev/fun/read_cov_L.hpp | 2 +- stan/math/rev/fun/read_cov_matrix.hpp | 2 +- 4 files changed, 24 insertions(+), 20 deletions(-) diff --git a/stan/math/rev/fun/cov_matrix_constrain.hpp b/stan/math/rev/fun/cov_matrix_constrain.hpp index e8f4ccf6abf..f0e6eca08f5 100644 --- a/stan/math/rev/fun/cov_matrix_constrain.hpp +++ b/stan/math/rev/fun/cov_matrix_constrain.hpp @@ -39,7 +39,8 @@ var_value cov_matrix_constrain(const T& x, Eigen::Index K) { for (Eigen::Index m = 0; m < K; ++m) { L_val.row(m).head(m) = x.val().segment(i, m); i += m; - L_val.coeffRef(m, m) = exp(x.val().coeff(i++)); + L_val.coeffRef(m, m) = exp(x.val().coeff(i)); + i++; } var_value L = L_val; @@ -48,7 +49,8 @@ var_value cov_matrix_constrain(const T& x, Eigen::Index K) { Eigen::Index K = L.rows(); int i = x.size(); for (int m = K - 1; m >= 0; --m) { - x.adj().coeffRef(--i) += L.adj().coeff(m, m) * L.val().coeff(m, m); + i--; + x.adj().coeffRef(i) += L.adj().coeff(m, m) * L.val().coeff(m, m); i -= m; x.adj().segment(i, m) += L.adj().row(m).head(m); } @@ -80,15 +82,16 @@ var_value cov_matrix_constrain(const T& x, Eigen::Index K, check_size_match("cov_matrix_constrain", "x.size()", x.size(), "K + (K choose 2)", (K * (K + 1)) / 2); arena_t L_val = Eigen::MatrixXd::Zero(K, K); - int i = 0; + int pos = 0; for (Eigen::Index m = 0; m < K; ++m) { - L_val.row(m).head(m) = x.val().segment(i, m); - i += m; - L_val.coeffRef(m, m) = exp(x.val().coeff(i++)); + L_val.row(m).head(m) = x.val().segment(pos, m); + pos += m; + L_val.coeffRef(m, m) = exp(x.val().coeff(pos)); + pos++; } // Jacobian for complete transform, including exp() above - double lp_val = (K * LOG_TWO); // needless constant; want propto + double lp_val = (K * LOG_TWO); for (Eigen::Index k = 0; k < K; ++k) { // only +1 because index from 0 lp_val += (K - k + 1) * log(L_val.coeff(k, k)); @@ -103,11 +106,12 @@ var_value cov_matrix_constrain(const T& x, Eigen::Index K, for (Eigen::Index k = 0; k < K; ++k) { L.adj().coeffRef(k, k) += (K - k + 1) * lp.adj() / L.val().coeff(k, k); } - int i = x.size(); + int pos = x.size(); for (int m = K - 1; m >= 0; --m) { - x.adj().coeffRef(--i) += L.adj().coeff(m, m) * L.val().coeff(m, m); - i -= m; - x.adj().segment(i, m) += L.adj().row(m).head(m); + pos--; + x.adj().coeffRef(pos) += L.adj().coeff(m, m) * L.val().coeff(m, m); + pos -= m; + x.adj().segment(pos, m) += L.adj().row(m).head(m); } }); diff --git a/stan/math/rev/fun/read_corr_L.hpp b/stan/math/rev/fun/read_corr_L.hpp index ac66e09b82a..dcfdde602d3 100644 --- a/stan/math/rev/fun/read_corr_L.hpp +++ b/stan/math/rev/fun/read_corr_L.hpp @@ -48,17 +48,17 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) // Cholesky factor of correlation matrix arena_t L = Eigen::MatrixXd::Zero(K, K); - size_t position = 0; + size_t pos = 0; size_t pull = K - 1; L(0, 0) = 1.0; L.col(0).tail(pull) = CPCs.val().head(pull); arena_acc.tail(pull) = 1.0 - CPCs.val().head(pull).array().square(); for (size_t i = 1; i < (K - 1); i++) { - position += pull; + pos += pull; pull = K - 1 - i; - const auto& cpc_seg = CPCs.val().array().segment(position, pull); + const auto& cpc_seg = CPCs.val().array().segment(pos, pull); L.coeffRef(i, i) = sqrt(arena_acc.coeff(i - 1)); L.col(i).tail(pull) = cpc_seg * arena_acc.tail(pull).sqrt(); arena_acc.tail(pull) = arena_acc.tail(pull) * (1.0 - cpc_seg.square()); @@ -75,12 +75,12 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) acc_adj.coeffRef(K - 2) += 0.5 * L_res.adj().coeff(K - 1, K - 1) / L_res.val().coeff(K - 1, K - 1); - int position = CPCs.size() - 1; + int pos = CPCs.size() - 1; for (size_t i = K - 2; i > 0; --i) { int pull = K - 1 - i; - const auto& cpc_seg_val = CPCs.val().array().segment(position, pull); - auto cpc_seg_adj = CPCs.adj().array().segment(position, pull); + const auto& cpc_seg_val = CPCs.val().array().segment(pos, pull); + auto cpc_seg_adj = CPCs.adj().array().segment(pos, pull); acc_val.tail(pull) /= 1.0 - cpc_seg_val.square(); cpc_seg_adj -= 2.0 * acc_adj.tail(pull) * acc_val.tail(pull) * cpc_seg_val; @@ -93,7 +93,7 @@ auto read_corr_L(const T& CPCs, size_t K) { // on (-1, 1) acc_adj.coeffRef(i - 1) += 0.5 * L_res.adj().coeff(i, i) / L_res.val().coeff(i, i); - position -= pull + 1; + pos -= pull + 1; } int pull = K - 1; diff --git a/stan/math/rev/fun/read_cov_L.hpp b/stan/math/rev/fun/read_cov_L.hpp index aec016967cb..0c294742932 100644 --- a/stan/math/rev/fun/read_cov_L.hpp +++ b/stan/math/rev/fun/read_cov_L.hpp @@ -43,7 +43,7 @@ inline auto read_cov_L(const T_CPCs& CPCs, const T_sds& sds, size_t K = sds.size(); corr_L.adj() += sds.val().matrix().asDiagonal() * res.adj(); - sds.adj() += rows_dot_product(res.adj(), corr_L.val()); + sds.adj() += (res.adj().cwiseProduct(corr_L.val())).rowwise().sum(); sds.adj() += (K * log_prob.adj() / sds.val().array()).matrix(); }); diff --git a/stan/math/rev/fun/read_cov_matrix.hpp b/stan/math/rev/fun/read_cov_matrix.hpp index ceeee85f5d8..3885441079d 100644 --- a/stan/math/rev/fun/read_cov_matrix.hpp +++ b/stan/math/rev/fun/read_cov_matrix.hpp @@ -53,7 +53,7 @@ inline var_value read_cov_matrix(const T_CPCs& CPCs, reverse_pass_callback([sds, corr_L, prod]() mutable { corr_L.adj() += sds.val().matrix().asDiagonal() * prod.adj(); - sds.adj() += rows_dot_product(prod.adj(), corr_L.val()); + sds.adj() += (prod.adj().cwiseProduct(corr_L.val())).rowwise().sum(); }); return multiply_lower_tri_self_transpose(prod); From 647bd401680b1b299a64a7377609a6e8d821e85b Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 30 Dec 2020 16:40:24 -0500 Subject: [PATCH 29/30] Updated indexing in constrain functions (Issue #2101) --- stan/math/rev/fun/cholesky_corr_constrain.hpp | 24 ++++++++++++------- .../rev/fun/cholesky_factor_constrain.hpp | 12 ++++++---- 2 files changed, 24 insertions(+), 12 deletions(-) diff --git a/stan/math/rev/fun/cholesky_corr_constrain.hpp b/stan/math/rev/fun/cholesky_corr_constrain.hpp index 104465204b7..4e624b7f894 100644 --- a/stan/math/rev/fun/cholesky_corr_constrain.hpp +++ b/stan/math/rev/fun/cholesky_corr_constrain.hpp @@ -41,10 +41,12 @@ var_value cholesky_corr_constrain(const T& y, int K) { x_val.coeffRef(0, 0) = 1; int k = 0; for (int i = 1; i < K; ++i) { - x_val.coeffRef(i, 0) = z.val().coeff(k++); + x_val.coeffRef(i, 0) = z.val().coeff(k); + k++; double sum_sqs = square(x_val.coeff(i, 0)); for (int j = 1; j < i; ++j) { - x_val.coeffRef(i, j) = z.val().coeff(k++) * sqrt(1.0 - sum_sqs); + x_val.coeffRef(i, j) = z.val().coeff(k) * sqrt(1.0 - sum_sqs); + k++; sum_sqs += square(x_val.coeff(i, j)); } x_val.coeffRef(i, i) = sqrt(1.0 - sum_sqs); @@ -60,12 +62,14 @@ var_value cholesky_corr_constrain(const T& y, int K) { for (int j = i - 1; j > 0; --j) { x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); sum_sqs_val -= square(x.val().coeff(i, j)); - sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * z.val().coeff(--k) + k--; + sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * z.val().coeff(k) / sqrt(1.0 - sum_sqs_val); z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); } x.adj().coeffRef(i, 0) += 2 * sum_sqs_adj * x.val().coeff(i, 0); - z.adj().coeffRef(--k) += x.adj().coeffRef(i, 0); + k--; + z.adj().coeffRef(k) += x.adj().coeffRef(i, 0); } }); @@ -105,11 +109,13 @@ var_value cholesky_corr_constrain(const T& y, int K, int k = 0; double lp_val = 0.0; for (int i = 1; i < K; ++i) { - x_val.coeffRef(i, 0) = z.val().coeff(k++); + x_val.coeffRef(i, 0) = z.val().coeff(k); + k++; double sum_sqs = square(x_val.coeff(i, 0)); for (int j = 1; j < i; ++j) { lp_val += 0.5 * log1m(sum_sqs); - x_val.coeffRef(i, j) = z.val().coeff(k++) * sqrt(1.0 - sum_sqs); + x_val.coeffRef(i, j) = z.val().coeff(k) * sqrt(1.0 - sum_sqs); + k++; sum_sqs += square(x_val.coeff(i, j)); } x_val.coeffRef(i, i) = sqrt(1.0 - sum_sqs); @@ -127,13 +133,15 @@ var_value cholesky_corr_constrain(const T& y, int K, x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); sum_sqs_val -= square(x.val().coeff(i, j)); - sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * z.val().coeff(--k) + k--; + sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * z.val().coeff(k) / sqrt(1.0 - sum_sqs_val); z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); sum_sqs_adj -= 0.5 * lp.adj() / (1 - sum_sqs_val); } x.adj().coeffRef(i, 0) += 2 * sum_sqs_adj * x.val().coeff(i, 0); - z.adj().coeffRef(--k) += x.adj().coeffRef(i, 0); + k--; + z.adj().coeffRef(k) += x.adj().coeffRef(i, 0); } }); diff --git a/stan/math/rev/fun/cholesky_factor_constrain.hpp b/stan/math/rev/fun/cholesky_factor_constrain.hpp index d1e70f10083..93b626c7ae5 100644 --- a/stan/math/rev/fun/cholesky_factor_constrain.hpp +++ b/stan/math/rev/fun/cholesky_factor_constrain.hpp @@ -42,7 +42,8 @@ var_value cholesky_factor_constrain(const T& x, int M, int N) { for (int m = 0; m < N; ++m) { y_val.row(m).head(m) = x.val().segment(pos, m); pos += m; - y_val.coeffRef(m, m) = exp(x.val().coeff(pos++)); + y_val.coeffRef(m, m) = exp(x.val().coeff(pos)); + pos++; } for (int m = N; m < M; ++m) { @@ -60,7 +61,8 @@ var_value cholesky_factor_constrain(const T& x, int M, int N) { } for (int m = N - 1; m >= 0; --m) { - x.adj().coeffRef(--pos) += y.adj().coeff(m, m) * y.val().coeff(m, m); + pos--; + x.adj().coeffRef(pos) += y.adj().coeff(m, m) * y.val().coeff(m, m); pos -= m; x.adj().segment(pos, m) += y.adj().row(m).head(m); } @@ -94,7 +96,8 @@ var_value cholesky_factor_constrain(const T& x, int M, int N, double lp_val = 0.0; for (int n = 0; n < N; ++n) { pos += n; - lp_val += x.val().coeff(pos++); + lp_val += x.val().coeff(pos); + pos++; } lp += lp_val; @@ -102,7 +105,8 @@ var_value cholesky_factor_constrain(const T& x, int M, int N, int pos = 0; for (int n = 0; n < N; ++n) { pos += n; - x.adj().coeffRef(pos++) += lp.adj(); + x.adj().coeffRef(pos) += lp.adj(); + pos++; } }); From 523100ddf38aae0189392d853bdae5f585ea6142 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 30 Dec 2020 21:41:21 +0000 Subject: [PATCH 30/30] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/cholesky_corr_constrain.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/rev/fun/cholesky_corr_constrain.hpp b/stan/math/rev/fun/cholesky_corr_constrain.hpp index 4e624b7f894..f65c32c4770 100644 --- a/stan/math/rev/fun/cholesky_corr_constrain.hpp +++ b/stan/math/rev/fun/cholesky_corr_constrain.hpp @@ -62,7 +62,7 @@ var_value cholesky_corr_constrain(const T& y, int K) { for (int j = i - 1; j > 0; --j) { x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); sum_sqs_val -= square(x.val().coeff(i, j)); - k--; + k--; sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * z.val().coeff(k) / sqrt(1.0 - sum_sqs_val); z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val); @@ -133,7 +133,7 @@ var_value cholesky_corr_constrain(const T& y, int K, x.adj().coeffRef(i, j) += 2 * sum_sqs_adj * x.val().coeff(i, j); sum_sqs_val -= square(x.val().coeff(i, j)); - k--; + k--; sum_sqs_adj += -0.5 * x.adj().coeffRef(i, j) * z.val().coeff(k) / sqrt(1.0 - sum_sqs_val); z.adj().coeffRef(k) += x.adj().coeffRef(i, j) * sqrt(1.0 - sum_sqs_val);