diff --git a/stan/math/prim/fun/corr_matrix_constrain.hpp b/stan/math/prim/fun/corr_matrix_constrain.hpp index 0bb81772eb6..5b7e6f941cf 100644 --- a/stan/math/prim/fun/corr_matrix_constrain.hpp +++ b/stan/math/prim/fun/corr_matrix_constrain.hpp @@ -65,7 +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 +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/prim/fun/exp.hpp b/stan/math/prim/fun/exp.hpp index 6766280f6af..384d8d3c7d4 100644 --- a/stan/math/prim/fun/exp.hpp +++ b/stan/math/prim/fun/exp.hpp @@ -42,7 +42,8 @@ struct exp_fun { template < typename Container, require_not_container_st* = nullptr, - require_not_nonscalar_prim_or_rev_kernel_expression_t* = 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/inv_logit.hpp b/stan/math/prim/fun/inv_logit.hpp index ee6f08a0118..7c8355f9ffe 100644 --- a/stan/math/prim/fun/inv_logit.hpp +++ b/stan/math/prim/fun/inv_logit.hpp @@ -81,8 +81,9 @@ struct inv_logit_fun { * @param x container * @return Inverse logit applied to each value in x. */ -template * = 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/log.hpp b/stan/math/prim/fun/log.hpp index 94252b27dee..ba256709b95 100644 --- a/stan/math/prim/fun/log.hpp +++ b/stan/math/prim/fun/log.hpp @@ -45,6 +45,7 @@ struct log_fun { template < 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/log1m.hpp b/stan/math/prim/fun/log1m.hpp index c8607501535..2b81dbd1edb 100644 --- a/stan/math/prim/fun/log1m.hpp +++ b/stan/math/prim/fun/log1m.hpp @@ -67,8 +67,9 @@ struct log1m_fun { * @param x container * @return Natural log of 1 minus each value in x. */ -template * = 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/prim/fun/multiply_log.hpp b/stan/math/prim/fun/multiply_log.hpp index 664fb545327..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); } @@ -69,7 +65,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> +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/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/prim/fun/square.hpp b/stan/math/prim/fun/square.hpp index 96e0426ebd2..abb05be9cce 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 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/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 eaf819cfc2c..8d51d0c55e2 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -29,12 +29,17 @@ #include #include #include +#include +#include #include #include #include #include #include +#include +#include #include +#include #include #include #include @@ -127,6 +132,10 @@ #include #include #include +#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..f65c32c4770 --- /dev/null +++ b/stan/math/rev/fun/cholesky_corr_constrain.hpp @@ -0,0 +1,153 @@ +#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 { + +/** + * 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; + 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 = Eigen::MatrixXd::Zero(K, 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); + 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); + k++; + 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)); + 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); + k--; + z.adj().coeffRef(k) += x.adj().coeffRef(i, 0); + } + }); + + return x; +} + +/** + * 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) { + 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 = Eigen::MatrixXd::Zero(K, K); + + 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); + 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); + k++; + 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)); + + 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); + k--; + z.adj().coeffRef(k) += x.adj().coeffRef(i, 0); + } + }); + + return x; +} + +} // namespace math +} // namespace stan +#endif 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..93b626c7ae5 --- /dev/null +++ b/stan/math/rev/fun/cholesky_factor_constrain.hpp @@ -0,0 +1,119 @@ +#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 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 + * @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 = 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)); + pos++; + } + + 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) { + 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); + } + }); + + 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 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[out] lp Log density 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); + 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(); + pos++; + } + }); + + return cholesky_factor_constrain(x, M, N); +} + +} // namespace math +} // namespace stan + +#endif 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..f28d7ed5e01 --- /dev/null +++ b/stan/math/rev/fun/corr_matrix_constrain.hpp @@ -0,0 +1,80 @@ +#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 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 + * 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 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. + */ +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/stan/math/rev/fun/cov_matrix_constrain.hpp b/stan/math/rev/fun/cov_matrix_constrain.hpp new file mode 100644 index 00000000000..f0e6eca08f5 --- /dev/null +++ b/stan/math/rev/fun/cov_matrix_constrain.hpp @@ -0,0 +1,124 @@ +#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 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. + * @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; + + 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)); + i++; + } + + 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) { + 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); + } + }); + + 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 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 + * @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; + + 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 pos = 0; + for (Eigen::Index m = 0; m < K; ++m) { + 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); + 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)); + } + + lp += lp_val; + + 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 pos = x.size(); + for (int m = K - 1; m >= 0; --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); + } + }); + + return multiply_lower_tri_self_transpose(L); +} + +} // namespace math +} // namespace stan + +#endif 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..29ef3a2b0fe --- /dev/null +++ b/stan/math/rev/fun/cov_matrix_constrain_lkj.hpp @@ -0,0 +1,77 @@ +#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 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. + * @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 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. + * @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 4a14a4811b1..53389bbb33e 100644 --- a/stan/math/rev/fun/exp.hpp +++ b/stan/math/rev/fun/exp.hpp @@ -57,6 +57,21 @@ 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) { + 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 } // namespace stan #endif diff --git a/stan/math/rev/fun/inv_logit.hpp b/stan/math/rev/fun/inv_logit.hpp index 827f29f8249..e506178d695 100644 --- a/stan/math/rev/fun/inv_logit.hpp +++ b/stan/math/rev/fun/inv_logit.hpp @@ -33,6 +33,21 @@ 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) { + 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 } // namespace stan #endif diff --git a/stan/math/rev/fun/log.hpp b/stan/math/rev/fun/log.hpp index 89fcd529539..6eaf4ed9d0c 100644 --- a/stan/math/rev/fun/log.hpp +++ b/stan/math/rev/fun/log.hpp @@ -64,6 +64,21 @@ 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) { + 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 } // namespace stan #endif diff --git a/stan/math/rev/fun/log1m.hpp b/stan/math/rev/fun/log1m.hpp index d8c4d3abbd0..70306cbafe9 100644 --- a/stan/math/rev/fun/log1m.hpp +++ b/stan/math/rev/fun/log1m.hpp @@ -28,6 +28,21 @@ 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) { + 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 } // namespace stan #endif diff --git a/stan/math/rev/fun/multiply_log.hpp b/stan/math/rev/fun/multiply_log.hpp index 8f8546af335..5b8850d0f66 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 @@ -19,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 { @@ -38,24 +31,14 @@ 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 { public: 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_; - } - } + void chain() { bvi_->adj_ += adj_ * ad_ / bvi_->val_; } }; } // namespace internal @@ -64,8 +47,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. @@ -91,8 +73,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. @@ -105,6 +86,149 @@ 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) { + 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( + 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(); + }); + } 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(); + }); + } 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 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) { + using std::log; + + if (!is_constant::value && !is_constant::value) { + arena_t> arena_a = a; + var arena_b = 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(); + }); + } 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)); + }); + } else { + arena_t> arena_a = value_of(a); + var 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() + += (res.adj().array() * arena_a.array()).sum() / arena_b.val(); + }); + } +} + +/** + * 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) { + if (!is_constant::value && !is_constant::value) { + var 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 { + 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; + 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() + += (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(); + }); + } +} + } // namespace math } // namespace stan #endif 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_L.hpp b/stan/math/rev/fun/read_corr_L.hpp new file mode 100644 index 00000000000..dcfdde602d3 --- /dev/null +++ b/stan/math/rev/fun/read_corr_L.hpp @@ -0,0 +1,179 @@ +#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 vector (must be a `var_value` where `S` + * 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; + arena_t arena_acc = Eigen::ArrayXd::Ones(K - 1); + // Cholesky factor of correlation matrix + arena_t L = Eigen::MatrixXd::Zero(K, K); + + 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++) { + pos += pull; + pull = K - 1 - i; + + 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()); + } + + L.coeffRef(K - 1, K - 1) = sqrt(arena_acc.coeff(K - 2)); + + 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.coeffRef(K - 2) += 0.5 * L_res.adj().coeff(K - 1, K - 1) + / L_res.val().coeff(K - 1, K - 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(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; + 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.coeffRef(i - 1) + += 0.5 * L_res.adj().coeff(i, i) / L_res.val().coeff(i, i); + + pos -= 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 ret_type(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 vector (must be a `var_value` where `S` + * 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().coeff(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().coeffRef(pos) -= 2.0 * acc_adj * (K - k - 1) + * CPCs.val().coeff(pos) + / (1.0 - square(CPCs.val().coeff(pos))); + pos--; + } + } + }); + + return read_corr_L(CPCs, K); +} + +} // namespace math +} // namespace stan +#endif 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..a5d8281740d --- /dev/null +++ b/stan/math/rev/fun/read_corr_matrix.hpp @@ -0,0 +1,68 @@ +#ifndef STAN_MATH_REV_FUN_READ_CORR_MATRIX_HPP +#define STAN_MATH_REV_FUN_READ_CORR_MATRIX_HPP + +#include +#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 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 + * 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 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. + * @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_L.hpp b/stan/math/rev/fun/read_cov_L.hpp new file mode 100644 index 00000000000..0c294742932 --- /dev/null +++ b/stan/math/rev/fun/read_cov_L.hpp @@ -0,0 +1,57 @@ +#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 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 + * @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() += (res.adj().cwiseProduct(corr_L.val())).rowwise().sum(); + + sds.adj() += (K * log_prob.adj() / sds.val().array()).matrix(); + }); + + return res; +} + +} // 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..3885441079d --- /dev/null +++ b/stan/math/rev/fun/read_cov_matrix.hpp @@ -0,0 +1,65 @@ +#ifndef STAN_MATH_REV_FUN_READ_COV_MATRIX_HPP +#define STAN_MATH_REV_FUN_READ_COV_MATRIX_HPP + +#include +#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 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 + * @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 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) + */ +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() += (prod.adj().cwiseProduct(corr_L.val())).rowwise().sum(); + }); + + return multiply_lower_tri_self_transpose(prod); +} + +} // namespace math +} // namespace stan + +#endif 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/square.hpp b/stan/math/rev/fun/square.hpp index 7de37294524..b502aa8d807 100644 --- a/stan/math/rev/fun/square.hpp +++ b/stan/math/rev/fun/square.hpp @@ -44,6 +44,22 @@ 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) { + 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 } // namespace stan #endif diff --git a/stan/math/rev/fun/tanh.hpp b/stan/math/rev/fun/tanh.hpp index f6feecf1fd8..39eb1bc9579 100644 --- a/stan/math/rev/fun/tanh.hpp +++ b/stan/math/rev/fun/tanh.hpp @@ -2,6 +2,7 @@ #define STAN_MATH_REV_FUN_TANH_HPP #include +#include #include #include #include @@ -46,12 +47,11 @@ inline var tanh(const var& a) { } /** - * Return the tanhgent of a radian-scaled variable (cmath). + * Return the hyperbolic tangent of elements of a * - * - * @tparam Varmat a `var_value` with inner Eigen type - * @param a Variable for radians of angle. - * @return Tangent of variable. + * @tparam T type of a + * @param a argument + * @return elementwise hyperbolic tangent of a */ template * = nullptr> inline auto tanh(const VarMat& a) { 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..eadb61e5468 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,53 @@ 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); +} 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..0350fb4465b --- /dev/null +++ b/test/unit/math/mix/fun/cholesky_factor_constrain_test.cpp @@ -0,0 +1,53 @@ +#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); +} 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..ff8e9f29b37 --- /dev/null +++ b/test/unit/math/mix/fun/corr_constrain_test.cpp @@ -0,0 +1,73 @@ +#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/corr_matrix_constrain_test.cpp b/test/unit/math/mix/fun/corr_matrix_constrain_test.cpp index 25d01fc4027..10dbb8990b5 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,55 @@ 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 new file mode 100644 index 00000000000..42fd1668ac9 --- /dev/null +++ b/test/unit/math/mix/fun/cov_matrix_constrain_lkj_test.cpp @@ -0,0 +1,53 @@ +#include + +TEST(mathMixMatFun, cov_matrix_constrain_lkj) { + 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); +} 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..e8785c3b91c 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); +} diff --git a/test/unit/math/mix/fun/exp_test.cpp b/test/unit/math/mix/fun/exp_test.cpp index 1a75ed232e8..ab326f1c540 100644 --- a/test/unit/math/mix/fun/exp_test.cpp +++ b/test/unit/math/mix/fun/exp_test.cpp @@ -9,4 +9,10 @@ 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); + + std::vector com_args = stan::test::internal::common_nonzero_args(); + std::vector args{-2.6, -0.5, 0.5, 1.5}; + + stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(com_args)); + stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(args)); } diff --git a/test/unit/math/mix/fun/inv_logit_test.cpp b/test/unit/math/mix/fun/inv_logit_test.cpp index ca0ae131503..8fee66f7eeb 100644 --- a/test/unit/math/mix/fun/inv_logit_test.cpp +++ b/test/unit/math/mix/fun/inv_logit_test.cpp @@ -5,4 +5,10 @@ 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); + + std::vector com_args = stan::test::internal::common_nonzero_args(); + std::vector args{-2.6, -0.5, 0.5, 1.5}; + + stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(com_args)); + stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(args)); } diff --git a/test/unit/math/mix/fun/log1m_test.cpp b/test/unit/math/mix/fun/log1m_test.cpp index d52ba182a9c..288cb52c6e0 100644 --- a/test/unit/math/mix/fun/log1m_test.cpp +++ b/test/unit/math/mix/fun/log1m_test.cpp @@ -4,4 +4,10 @@ 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); + + std::vector com_args = stan::test::internal::common_nonzero_args(); + std::vector args{-2.6, -0.5, 0.5, 1.5}; + + stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(com_args)); + stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(args)); } diff --git a/test/unit/math/mix/fun/log_test.cpp b/test/unit/math/mix/fun/log_test.cpp index 21ff2c8fa03..f6b4e5b8871 100644 --- a/test/unit/math/mix/fun/log_test.cpp +++ b/test/unit/math/mix/fun/log_test.cpp @@ -24,4 +24,10 @@ 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) + + std::vector com_args = stan::test::internal::common_nonzero_args(); + std::vector args{0.1, 2.5, 5.5}; + + stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(com_args)); + stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(args)); } 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); } 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..bbd4450a64c --- /dev/null +++ b/test/unit/math/mix/fun/read_corr_L_test.cpp @@ -0,0 +1,51 @@ +#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/read_corr_matrix_test.cpp b/test/unit/math/mix/fun/read_corr_matrix_test.cpp new file mode 100644 index 00000000000..e50a766bbf3 --- /dev/null +++ b/test/unit/math/mix/fun/read_corr_matrix_test.cpp @@ -0,0 +1,51 @@ +#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_L_test.cpp b/test/unit/math/mix/fun/read_cov_L_test.cpp new file mode 100644 index 00000000000..321aa521998 --- /dev/null +++ b/test/unit/math/mix/fun/read_cov_L_test.cpp @@ -0,0 +1,42 @@ +#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); +} 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..8370f6e50d5 --- /dev/null +++ b/test/unit/math/mix/fun/read_cov_matrix_test.cpp @@ -0,0 +1,65 @@ +#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); +} diff --git a/test/unit/math/mix/fun/square_test.cpp b/test/unit/math/mix/fun/square_test.cpp index 2ef89b26009..5577bcda7a0 100644 --- a/test/unit/math/mix/fun/square_test.cpp +++ b/test/unit/math/mix/fun/square_test.cpp @@ -5,4 +5,10 @@ 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); + + std::vector com_args = stan::test::internal::common_nonzero_args(); + std::vector args{-2.6, -0.5, 0.5, 1.5}; + + stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(com_args)); + stan::test::expect_ad_vector_matvar(f, stan::math::to_vector(args)); } diff --git a/test/unit/math/mix/fun/tanh_test.cpp b/test/unit/math/mix/fun/tanh_test.cpp index c0472d2b24e..be8e0326528 100644 --- a/test/unit/math/mix/fun/tanh_test.cpp +++ b/test/unit/math/mix/fun/tanh_test.cpp @@ -19,7 +19,7 @@ TEST(mathMixMatFun, tanh_varmat) { return tanh(x1); }; std::vector com_args = common_nonzero_args(); - std::vector args{-2.6, -2, -1.2, -0.5, 0.5, 1.5}; + std::vector args{-2.6, -0.5, 0.5, 1.5}; auto all_args = vec_concat(com_args, args); Eigen::VectorXd A(all_args.size()); for (int i = 0; i < all_args.size(); ++i) { diff --git a/test/unit/math/test_ad.hpp b/test/unit/math/test_ad.hpp index f1af2c7caf3..66463162bd6 100644 --- a/test/unit/math/test_ad.hpp +++ b/test/unit/math/test_ad.hpp @@ -2110,7 +2110,7 @@ void expect_ad_matvar(const ad_tolerances& tols, const F& f, */ template * = nullptr> -void expect_ad_vector_matvar(const F& f, EigVec& x) { +void expect_ad_vector_matvar(const F& f, const EigVec& x) { ad_tolerances tols; expect_ad_vector_matvar(tols, f, x); }