Automatic Differentiation
 
Loading...
Searching...
No Matches
cholesky_factor_constrain.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_CONSTRAINT_CHOLESKY_FACTOR_CONSTRAIN_HPP
2#define STAN_MATH_PRIM_CONSTRAINT_CHOLESKY_FACTOR_CONSTRAIN_HPP
3
9#include <cmath>
10#include <stdexcept>
11#include <vector>
12
13namespace stan {
14namespace math {
15
28template <typename T, require_eigen_col_vector_t<T>* = nullptr>
29inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
30cholesky_factor_constrain(const T& x, int M, int N) {
31 using std::exp;
32 using T_scalar = value_type_t<T>;
33 check_greater_or_equal("cholesky_factor_constrain",
34 "num rows (must be greater or equal to num cols)", M,
35 N);
36 check_size_match("cholesky_factor_constrain", "constrain size", x.size(),
37 "((N * (N + 1)) / 2 + (M - N) * N)",
38 ((N * (N + 1)) / 2 + (M - N) * N));
39 Eigen::Matrix<T_scalar, Eigen::Dynamic, Eigen::Dynamic> y(M, N);
40 int pos = 0;
41
42 const auto& x_ref = to_ref(x);
43 for (int m = 0; m < N; ++m) {
44 y.row(m).head(m) = x_ref.segment(pos, m);
45 pos += m;
46 y.coeffRef(m, m) = exp(x_ref.coeff(pos++));
47 y.row(m).tail(N - m - 1).setZero();
48 }
49
50 for (int m = N; m < M; ++m) {
51 y.row(m) = x_ref.segment(pos, N);
52 pos += N;
53 }
54 return y;
55}
56
73template <typename T, require_eigen_vector_t<T>* = nullptr>
74inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
75cholesky_factor_constrain(const T& x, int M, int N, return_type_t<T>& lp) {
76 check_size_match("cholesky_factor_constrain", "x.size()", x.size(),
77 "((N * (N + 1)) / 2 + (M - N) * N)",
78 ((N * (N + 1)) / 2 + (M - N) * N));
79 int pos = 0;
80 const auto& x_ref = to_ref(x);
81 for (int n = 0; n < N; ++n) {
82 pos += n;
83 lp += x_ref.coeff(pos++);
84 }
85 return cholesky_factor_constrain(x_ref, M, N);
86}
87
108template <bool Jacobian, typename T, require_not_std_vector_t<T>* = nullptr>
109inline auto cholesky_factor_constrain(const T& x, int M, int N,
110 return_type_t<T>& lp) {
111 if (Jacobian) {
112 return cholesky_factor_constrain(x, M, N, lp);
113 } else {
114 return cholesky_factor_constrain(x, M, N);
115 }
116}
117
138template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr>
139inline auto cholesky_factor_constrain(const T& x, int M, int N,
140 return_type_t<T>& lp) {
141 return apply_vector_unary<T>::apply(x, [&lp, M, N](auto&& v) {
142 return cholesky_factor_constrain<Jacobian>(v, M, N, lp);
143 });
144}
145
146} // namespace math
147} // namespace stan
148
149#endif
typename value_type< T >::type value_type_t
Helper function for accessing underlying type.
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
void check_greater_or_equal(const char *function, const char *name, const T_y &y, const T_low &low, Idxs... idxs)
Throw an exception if y is not greater or equal than low.
Eigen::Matrix< value_type_t< T >, Eigen::Dynamic, Eigen::Dynamic > cholesky_factor_constrain(const T &x, int M, int N)
Return the Cholesky factor of the specified size read from the specified vector.
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 ...