Loading [MathJax]/extensions/TeX/AMSsymbols.js
Automatic Differentiation
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Loading...
Searching...
No Matches
sum_to_zero_constrain.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP
2#define STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP
3
9#include <cmath>
10
11namespace stan {
12namespace math {
13
38template <typename Vec, require_eigen_col_vector_t<Vec>* = nullptr,
39 require_not_st_var<Vec>* = nullptr>
41 const auto N = y.size();
42
43 plain_type_t<Vec> z = Eigen::VectorXd::Zero(N + 1);
44 if (unlikely(N == 0)) {
45 return z;
46 }
47
48 auto&& y_ref = to_ref(y);
49
50 value_type_t<Vec> sum_w(0);
51 for (int i = N; i > 0; --i) {
52 double n = static_cast<double>(i);
53 auto w = y_ref.coeff(i - 1) * inv_sqrt(n * (n + 1));
54 sum_w += w;
55
56 z.coeffRef(i - 1) += sum_w;
57 z.coeffRef(i) -= w * n;
58 }
59
60 return z;
61}
62
73template <typename Mat, require_eigen_matrix_dynamic_t<Mat>* = nullptr,
74 require_not_st_var<Mat>* = nullptr>
76 const auto N = x.rows();
77 const auto M = x.cols();
78
79 plain_type_t<Mat> Z = Eigen::MatrixXd::Zero(N + 1, M + 1);
80 if (unlikely(N == 0 || M == 0)) {
81 return Z;
82 }
83 auto&& x_ref = to_ref(x);
84
85 Eigen::Matrix<value_type_t<Mat>, -1, 1> beta = Eigen::VectorXd::Zero(N);
86
87 for (int j = M - 1; j >= 0; --j) {
88 value_type_t<Mat> ax_previous(0);
89
90 double a_j = inv_sqrt((j + 1.0) * (j + 2.0));
91 double b_j = (j + 1.0) * a_j;
92
93 for (int i = N - 1; i >= 0; --i) {
94 double a_i = inv_sqrt((i + 1.0) * (i + 2.0));
95 double b_i = (i + 1.0) * a_i;
96
97 auto b_i_x = b_i * x_ref.coeff(i, j) - ax_previous;
98
99 Z.coeffRef(i, j) = (b_j * b_i_x) - beta.coeff(i);
100 beta.coeffRef(i) += a_j * b_i_x;
101
102 Z.coeffRef(N, j) -= Z.coeff(i, j);
103 Z.coeffRef(i, M) -= Z.coeff(i, j);
104
105 ax_previous += a_i * x_ref.coeff(i, j);
106 }
107 Z.coeffRef(N, M) -= Z.coeff(N, j);
108 }
109
110 return Z;
111}
112
125template <typename T, typename Lp, require_not_st_var<T>* = nullptr>
126inline plain_type_t<T> sum_to_zero_constrain(const T& y, Lp& lp) {
127 return sum_to_zero_constrain(y);
128}
129
140template <typename T, require_std_vector_t<T>* = nullptr>
141inline auto sum_to_zero_constrain(const T& y) {
143 y, [](auto&& v) { return sum_to_zero_constrain(v); });
144}
145
158template <typename T, typename Lp, require_std_vector_t<T>* = nullptr,
159 require_convertible_t<return_type_t<T>, Lp>* = nullptr>
160inline auto sum_to_zero_constrain(const T& y, Lp& lp) {
162 y, [](auto&& v) { return sum_to_zero_constrain(v); });
163}
164
177template <bool Jacobian, typename T, typename Lp>
178inline plain_type_t<T> sum_to_zero_constrain(const T& y, Lp& lp) {
179 return sum_to_zero_constrain(y);
180}
181
182} // namespace math
183} // namespace stan
184
185#endif
#define unlikely(x)
typename value_type< T >::type value_type_t
Helper function for accessing underlying type.
ref_type_t< T && > to_ref(T &&a)
This evaluates expensive Eigen expressions.
Definition to_ref.hpp:17
plain_type_t< Vec > sum_to_zero_constrain(const Vec &y)
Return a vector with sum zero corresponding to the specified free vector.
fvar< T > beta(const fvar< T > &x1, const fvar< T > &x2)
Return fvar with the beta function applied to the specified arguments and its gradient.
Definition beta.hpp:51
fvar< T > inv_sqrt(const fvar< T > &x)
Definition inv_sqrt.hpp:14
typename plain_type< T >::type plain_type_t
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...