Automatic Differentiation
 
Loading...
Searching...
No Matches
cov_matrix_constrain.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_CONSTRAINT_COV_MATRIX_CONSTRAIN_HPP
2#define STAN_MATH_PRIM_CONSTRAINT_COV_MATRIX_CONSTRAIN_HPP
3
11#include <cmath>
12
13namespace stan {
14namespace math {
15
30template <typename T, require_eigen_col_vector_t<T>* = nullptr>
31inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
32cov_matrix_constrain(const T& x, Eigen::Index K) {
33 using Eigen::Dynamic;
34 using Eigen::Matrix;
35 using std::exp;
36
37 Matrix<value_type_t<T>, Dynamic, Dynamic> L(K, K);
38 check_size_match("cov_matrix_constrain", "x.size()", x.size(),
39 "K + (K choose 2)", (K * (K + 1)) / 2);
40 const auto& x_ref = to_ref(x);
41 int i = 0;
42 for (Eigen::Index m = 0; m < K; ++m) {
43 L.row(m).head(m) = x_ref.segment(i, m);
44 i += m;
45 L.coeffRef(m, m) = exp(x_ref.coeff(i++));
46 L.row(m).tail(K - m - 1).setZero();
47 }
49}
50
63template <typename T, require_eigen_col_vector_t<T>* = nullptr>
64inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
65cov_matrix_constrain(const T& x, Eigen::Index K, return_type_t<T>& lp) {
66 using Eigen::Dynamic;
67 using Eigen::Matrix;
68 using std::exp;
69 using std::log;
70 check_size_match("cov_matrix_constrain", "x.size()", x.size(),
71 "K + (K choose 2)", (K * (K + 1)) / 2);
72 Matrix<value_type_t<T>, Dynamic, Dynamic> L(K, K);
73 const auto& x_ref = to_ref(x);
74 int i = 0;
75 for (Eigen::Index m = 0; m < K; ++m) {
76 L.row(m).head(m) = x_ref.segment(i, m);
77 i += m;
78 L.coeffRef(m, m) = exp(x_ref.coeff(i++));
79 L.row(m).tail(K - m - 1).setZero();
80 }
81 // Jacobian for complete transform, including exp() above
82 lp += (K * LOG_TWO); // needless constant; want propto
83 for (Eigen::Index k = 0; k < K; ++k) {
84 lp += (K - k + 1) * log(L.coeff(k, k)); // only +1 because index from 0
85 }
87}
88
107template <bool Jacobian, typename T, require_not_std_vector_t<T>* = nullptr>
108inline auto cov_matrix_constrain(const T& x, Eigen::Index K,
109 return_type_t<T>& lp) {
110 if (Jacobian) {
111 return cov_matrix_constrain(x, K, lp);
112 } else {
113 return cov_matrix_constrain(x, K);
114 }
115}
116
135template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr>
136inline auto cov_matrix_constrain(const T& x, Eigen::Index K,
137 return_type_t<T>& lp) {
138 return apply_vector_unary<T>::apply(x, [&lp, K](auto&& v) {
139 return cov_matrix_constrain<Jacobian>(v, K, lp);
140 });
141}
142
143} // namespace math
144} // namespace stan
145
146#endif
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
Eigen::Matrix< value_type_t< EigMat >, EigMat::RowsAtCompileTime, EigMat::RowsAtCompileTime > multiply_lower_tri_self_transpose(const EigMat &m)
Eigen::Matrix< value_type_t< T >, Eigen::Dynamic, Eigen::Dynamic > cov_matrix_constrain(const T &x, Eigen::Index K)
Return the symmetric, positive-definite matrix of dimensions K by K resulting from transforming the s...
fvar< T > log(const fvar< T > &x)
Definition log.hpp:15
static constexpr double LOG_TWO
The natural logarithm of 2, .
Definition constants.hpp:80
ref_type_t< T && > to_ref(T &&a)
This evaluates expensive Eigen expressions.
Definition to_ref.hpp:17
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Check if the provided sizes match.
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...