Skip to content

Commit 0f1ad25

Browse files
authored
Merge pull request #3151 from stan-dev/broaden-constraint-lp-template
Broaden the template bound on the lp argument to constraints
2 parents 5dc041b + 5aad82a commit 0f1ad25

40 files changed

+768
-470
lines changed

stan/math/prim/constraint/cholesky_corr_constrain.hpp

+45-27
Original file line numberDiff line numberDiff line change
@@ -43,9 +43,11 @@ cholesky_corr_constrain(const EigVec& y, int K) {
4343
}
4444

4545
// FIXME to match above after debugged
46-
template <typename EigVec, require_eigen_vector_t<EigVec>* = nullptr>
46+
template <typename EigVec, typename Lp,
47+
require_eigen_vector_t<EigVec>* = nullptr,
48+
require_convertible_t<return_type_t<EigVec>, Lp>* = nullptr>
4749
inline Eigen::Matrix<value_type_t<EigVec>, Eigen::Dynamic, Eigen::Dynamic>
48-
cholesky_corr_constrain(const EigVec& y, int K, return_type_t<EigVec>& lp) {
50+
cholesky_corr_constrain(const EigVec& y, int K, Lp& lp) {
4951
using Eigen::Dynamic;
5052
using Eigen::Matrix;
5153
using std::sqrt;
@@ -75,28 +77,39 @@ cholesky_corr_constrain(const EigVec& y, int K, return_type_t<EigVec>& lp) {
7577
}
7678

7779
/**
78-
* Return The cholesky of a `KxK` correlation matrix. If the `Jacobian`
79-
* parameter is `true`, the log density accumulator is incremented with the log
80-
* absolute Jacobian determinant of the transform. All of the transforms are
81-
* specified with their Jacobians in the *Stan Reference Manual* chapter
82-
* Constraint Transforms.
83-
* @tparam Jacobian if `true`, increment log density accumulator with log
84-
* absolute Jacobian determinant of constraining transform
85-
* @tparam T A type inheriting from `Eigen::DenseBase` or a `var_value` with
86-
* inner type inheriting from `Eigen::DenseBase` with compile time dynamic rows
87-
* and 1 column
80+
* Return The cholesky of a `KxK` correlation matrix.
81+
* This overload handles looping over the elements of a standard vector.
82+
* @tparam T A standard vector with inner type inheriting from
83+
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
84+
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
85+
* @param y Linearly Serialized vector of size `(K * (K - 1))/2` holding the
86+
* column major order elements of the lower triangurlar
87+
* @param K The size of the matrix to return
88+
*/
89+
template <typename T, require_std_vector_t<T>* = nullptr>
90+
inline auto cholesky_corr_constrain(const T& y, int K) {
91+
return apply_vector_unary<T>::apply(
92+
y, [K](auto&& v) { return cholesky_corr_constrain(v, K); });
93+
}
94+
95+
/**
96+
* Return The cholesky of a `KxK` correlation matrix.
97+
* This overload handles looping over the elements of a standard vector.
98+
* @tparam T A standard vector with inner type inheriting from
99+
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
100+
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
101+
* @tparam Lp Scalar type for the lp argument. The scalar type of T should be
102+
* convertable to this.
88103
* @param y Linearly Serialized vector of size `(K * (K - 1))/2` holding the
89104
* column major order elements of the lower triangurlar
90105
* @param K The size of the matrix to return
91106
* @param[in,out] lp log density accumulator
92107
*/
93-
template <bool Jacobian, typename T, require_not_std_vector_t<T>* = nullptr>
94-
inline auto cholesky_corr_constrain(const T& y, int K, return_type_t<T>& lp) {
95-
if (Jacobian) {
96-
return cholesky_corr_constrain(y, K, lp);
97-
} else {
98-
return cholesky_corr_constrain(y, K);
99-
}
108+
template <typename T, typename Lp, require_std_vector_t<T>* = nullptr,
109+
require_convertible_t<return_type_t<T>, Lp>* = nullptr>
110+
inline auto cholesky_corr_constrain(const T& y, int K, Lp& lp) {
111+
return apply_vector_unary<T>::apply(
112+
y, [&lp, K](auto&& v) { return cholesky_corr_constrain(v, K, lp); });
100113
}
101114

102115
/**
@@ -107,19 +120,24 @@ inline auto cholesky_corr_constrain(const T& y, int K, return_type_t<T>& lp) {
107120
* Constraint Transforms.
108121
* @tparam Jacobian if `true`, increment log density accumulator with log
109122
* absolute Jacobian determinant of constraining transform
110-
* @tparam T A standard vector with inner type inheriting from
111-
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
112-
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
123+
* @tparam T A type inheriting from `Eigen::DenseBase` or a `var_value` with
124+
* inner type inheriting from `Eigen::DenseBase` with compile time dynamic rows
125+
* and 1 column, or a standard vector thereof
126+
* @tparam Lp A scalar type for the lp argument. The scalar type of T should be
127+
* convertable to this.
113128
* @param y Linearly Serialized vector of size `(K * (K - 1))/2` holding the
114129
* column major order elements of the lower triangurlar
115130
* @param K The size of the matrix to return
116131
* @param[in,out] lp log density accumulator
117132
*/
118-
template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr>
119-
inline auto cholesky_corr_constrain(const T& y, int K, return_type_t<T>& lp) {
120-
return apply_vector_unary<T>::apply(y, [&lp, K](auto&& v) {
121-
return cholesky_corr_constrain<Jacobian>(v, K, lp);
122-
});
133+
template <bool Jacobian, typename T, typename Lp,
134+
require_convertible_t<return_type_t<T>, Lp>* = nullptr>
135+
inline auto cholesky_corr_constrain(const T& y, int K, Lp& lp) {
136+
if constexpr (Jacobian) {
137+
return cholesky_corr_constrain(y, K, lp);
138+
} else {
139+
return cholesky_corr_constrain(y, K);
140+
}
123141
}
124142

125143
} // namespace math

stan/math/prim/constraint/cholesky_corr_free.hpp

+5-4
Original file line numberDiff line numberDiff line change
@@ -11,9 +11,10 @@ namespace stan {
1111
namespace math {
1212

1313
template <typename T, require_eigen_t<T>* = nullptr>
14-
auto cholesky_corr_free(const T& x) {
14+
inline auto cholesky_corr_free(const T& x) {
1515
using Eigen::Dynamic;
1616
using Eigen::Matrix;
17+
using std::sqrt;
1718

1819
check_square("cholesky_corr_free", "x", x);
1920
// should validate lower-triangular, unit lengths
@@ -24,9 +25,9 @@ auto cholesky_corr_free(const T& x) {
2425
int k = 0;
2526
for (int i = 1; i < x.rows(); ++i) {
2627
z.coeffRef(k++) = corr_free(x_ref.coeff(i, 0));
27-
double sum_sqs = square(x_ref.coeff(i, 0));
28+
auto sum_sqs = square(x_ref.coeff(i, 0));
2829
for (int j = 1; j < i; ++j) {
29-
z.coeffRef(k++) = corr_free(x_ref.coeff(i, j) / std::sqrt(1.0 - sum_sqs));
30+
z.coeffRef(k++) = corr_free(x_ref.coeff(i, j) / sqrt(1.0 - sum_sqs));
3031
sum_sqs += square(x_ref.coeff(i, j));
3132
}
3233
}
@@ -41,7 +42,7 @@ auto cholesky_corr_free(const T& x) {
4142
* @param x The standard vector to untransform.
4243
*/
4344
template <typename T, require_std_vector_t<T>* = nullptr>
44-
auto cholesky_corr_free(const T& x) {
45+
inline auto cholesky_corr_free(const T& x) {
4546
return apply_vector_unary<T>::apply(
4647
x, [](auto&& v) { return cholesky_corr_free(v); });
4748
}

stan/math/prim/constraint/cholesky_factor_constrain.hpp

+49-29
Original file line numberDiff line numberDiff line change
@@ -70,9 +70,10 @@ cholesky_factor_constrain(const T& x, int M, int N) {
7070
* determinant
7171
* @return Cholesky factor
7272
*/
73-
template <typename T, require_eigen_vector_t<T>* = nullptr>
73+
template <typename T, typename Lp, require_eigen_vector_t<T>* = nullptr,
74+
require_convertible_t<return_type_t<T>, Lp>* = nullptr>
7475
inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
75-
cholesky_factor_constrain(const T& x, int M, int N, return_type_t<T>& lp) {
76+
cholesky_factor_constrain(const T& x, int M, int N, Lp& lp) {
7677
check_size_match("cholesky_factor_constrain", "x.size()", x.size(),
7778
"((N * (N + 1)) / 2 + (M - N) * N)",
7879
((N * (N + 1)) / 2 + (M - N) * N));
@@ -88,31 +89,46 @@ cholesky_factor_constrain(const T& x, int M, int N, return_type_t<T>& lp) {
8889
/**
8990
* Return the Cholesky factor of the specified size read from the specified
9091
* vector. A total of (N choose 2) + N + N * (M - N) free parameters are
91-
* required to read an M by N Cholesky factor. If the `Jacobian` parameter is
92-
* `true`, the log density accumulator is incremented with the log absolute
93-
* Jacobian determinant of the transform. All of the transforms are specified
94-
* with their Jacobians in the *Stan Reference Manual* chapter Constraint
95-
* Transforms.
92+
* required to read an M by N Cholesky factor.
93+
* This overload handles looping over the elements of a standard vector.
9694
*
97-
* @tparam Jacobian if `true`, increment log density accumulator with log
98-
* absolute Jacobian determinant of constraining transform
99-
* @tparam T A type inheriting from `Eigen::DenseBase` or a `var_value` with
100-
* inner type inheriting from `Eigen::DenseBase` with compile time dynamic rows
101-
* and 1 column
95+
* @tparam T A standard vector with inner type inheriting from
96+
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
97+
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
98+
* @param x Vector of unconstrained values
99+
* @param M number of rows
100+
* @param N number of columns
101+
* @return Cholesky factor
102+
*/
103+
template <typename T, require_std_vector_t<T>* = nullptr>
104+
inline auto cholesky_factor_constrain(const T& x, int M, int N) {
105+
return apply_vector_unary<T>::apply(
106+
x, [M, N](auto&& v) { return cholesky_factor_constrain(v, M, N); });
107+
}
108+
109+
/**
110+
* Return the Cholesky factor of the specified size read from the specified
111+
* vector. A total of (N choose 2) + N + N * (M - N) free parameters are
112+
* required to read an M by N Cholesky factor.
113+
* This overload handles looping over the elements of a standard vector.
114+
*
115+
* @tparam T A standard vector with inner type inheriting from
116+
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
117+
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
118+
* @tparam Lp Scalar type for the lp argument. The scalar type of T should be
119+
* convertable to this.
102120
* @param x Vector of unconstrained values
103121
* @param M number of rows
104122
* @param N number of columns
105123
* @param[in,out] lp log density accumulator
106124
* @return Cholesky factor
107125
*/
108-
template <bool Jacobian, typename T, require_not_std_vector_t<T>* = nullptr>
109-
inline 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-
}
126+
template <typename T, typename Lp, require_std_vector_t<T>* = nullptr,
127+
require_convertible_t<return_type_t<T>, Lp>* = nullptr>
128+
inline auto cholesky_factor_constrain(const T& x, int M, int N, Lp& lp) {
129+
return apply_vector_unary<T>::apply(x, [&lp, M, N](auto&& v) {
130+
return cholesky_factor_constrain(v, M, N, lp);
131+
});
116132
}
117133

118134
/**
@@ -126,21 +142,25 @@ inline auto cholesky_factor_constrain(const T& x, int M, int N,
126142
*
127143
* @tparam Jacobian if `true`, increment log density accumulator with log
128144
* absolute Jacobian determinant of constraining transform
129-
* @tparam T A standard vector with inner type inheriting from
130-
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
131-
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
145+
* @tparam T A type inheriting from `Eigen::DenseBase` or a `var_value` with
146+
* inner type inheriting from `Eigen::DenseBase` with compile time dynamic rows
147+
* and 1 column
148+
* @tparam Lp A scalar type for the lp argument. The scalar type of T should be
149+
* convertable to this.
132150
* @param x Vector of unconstrained values
133151
* @param M number of rows
134152
* @param N number of columns
135153
* @param[in,out] lp log density accumulator
136154
* @return Cholesky factor
137155
*/
138-
template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr>
139-
inline 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-
});
156+
template <bool Jacobian, typename T, typename Lp,
157+
require_convertible_t<return_type_t<T>, Lp>* = nullptr>
158+
inline auto cholesky_factor_constrain(const T& x, int M, int N, Lp& lp) {
159+
if constexpr (Jacobian) {
160+
return cholesky_factor_constrain(x, M, N, lp);
161+
} else {
162+
return cholesky_factor_constrain(x, M, N);
163+
}
144164
}
145165

146166
} // namespace math

stan/math/prim/constraint/cholesky_factor_free.hpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,7 @@ Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, 1> cholesky_factor_free(
5656
* @param x The standard vector to untransform.
5757
*/
5858
template <typename T, require_std_vector_t<T>* = nullptr>
59-
auto cholesky_factor_free(const T& x) {
59+
inline auto cholesky_factor_free(const T& x) {
6060
return apply_vector_unary<T>::apply(
6161
x, [](auto&& v) { return cholesky_factor_free(v); });
6262
}

stan/math/prim/constraint/corr_constrain.hpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -66,7 +66,7 @@ inline auto corr_constrain(const T_x& x, T_lp& lp) {
6666
*/
6767
template <bool Jacobian, typename T_x, typename T_lp>
6868
inline auto corr_constrain(const T_x& x, T_lp& lp) {
69-
if (Jacobian) {
69+
if constexpr (Jacobian) {
7070
return corr_constrain(x, lp);
7171
} else {
7272
return corr_constrain(x);

stan/math/prim/constraint/corr_free.hpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ namespace math {
2525
* @return free scalar that transforms to the specified input
2626
*/
2727
template <typename T>
28-
inline T corr_free(const T& y) {
28+
inline plain_type_t<T> corr_free(const T& y) {
2929
check_bounded("lub_free", "Correlation variable", y, -1.0, 1.0);
3030
return atanh(y);
3131
}

stan/math/prim/constraint/corr_matrix_constrain.hpp

+53-32
Original file line numberDiff line numberDiff line change
@@ -61,13 +61,16 @@ corr_matrix_constrain(const T& x, Eigen::Index k) {
6161
*
6262
* @tparam T type of the vector (must be derived from \c Eigen::MatrixBase and
6363
* have one compile-time dimension equal to 1)
64+
* @tparam Lp A scalar type for the lp argument. The scalar type of T should be
65+
* convertable to this.
6466
* @param x Vector of unconstrained partial correlations.
6567
* @param k Dimensionality of returned correlation matrix.
6668
* @param lp Log probability reference to increment.
6769
*/
68-
template <typename T, require_eigen_col_vector_t<T>* = nullptr>
70+
template <typename T, typename Lp, require_eigen_col_vector_t<T>* = nullptr,
71+
require_convertible_t<return_type_t<T>, Lp>* = nullptr>
6972
inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
70-
corr_matrix_constrain(const T& x, Eigen::Index k, return_type_t<T>& lp) {
73+
corr_matrix_constrain(const T& x, Eigen::Index k, Lp& lp) {
7174
Eigen::Index k_choose_2 = (k * (k - 1)) / 2;
7275
check_size_match("cov_matrix_constrain", "x.size()", x.size(), "k_choose_2",
7376
k_choose_2);
@@ -79,28 +82,41 @@ corr_matrix_constrain(const T& x, Eigen::Index k, return_type_t<T>& lp) {
7982
* the specified vector of unconstrained values. The input vector must be of
8083
* length \f${k \choose 2} = \frac{k(k-1)}{2}\f$. The values in the input
8184
* vector represent unconstrained (partial) correlations among the dimensions.
82-
* If the `Jacobian` parameter is `true`, the log density accumulator is
83-
* incremented with the log absolute Jacobian determinant of the transform. All
84-
* of the transforms are specified with their Jacobians in the *Stan Reference
85-
* Manual* chapter Constraint Transforms.
85+
* This overload handles looping over the elements of a standard vector.
8686
*
87-
* @tparam Jacobian if `true`, increment log density accumulator with log
88-
* absolute Jacobian determinant of constraining transform
89-
* @tparam T A type inheriting from `Eigen::DenseBase` or a `var_value` with
90-
* inner type inheriting from `Eigen::DenseBase` with compile time dynamic rows
91-
* and 1 column
92-
* @param x Vector of unconstrained partial correlations
93-
* @param k Dimensionality of returned correlation matrix
94-
* @param[in,out] lp log density accumulator
87+
* @tparam T A standard vector with inner type inheriting from
88+
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
89+
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
90+
* @param y Vector of unconstrained partial correlations
91+
* @param K Dimensionality of returned correlation matrix
9592
*/
96-
template <bool Jacobian, typename T, require_not_std_vector_t<T>* = nullptr>
97-
inline auto corr_matrix_constrain(const T& x, Eigen::Index k,
98-
return_type_t<T>& lp) {
99-
if (Jacobian) {
100-
return corr_matrix_constrain(x, k, lp);
101-
} else {
102-
return corr_matrix_constrain(x, k);
103-
}
93+
template <typename T, require_std_vector_t<T>* = nullptr>
94+
inline auto corr_matrix_constrain(const T& y, int K) {
95+
return apply_vector_unary<T>::apply(
96+
y, [K](auto&& v) { return corr_matrix_constrain(v, K); });
97+
}
98+
99+
/**
100+
* Return the correlation matrix of the specified dimensionality derived from
101+
* the specified vector of unconstrained values. The input vector must be of
102+
* length \f${k \choose 2} = \frac{k(k-1)}{2}\f$. The values in the input
103+
* vector represent unconstrained (partial) correlations among the dimensions.
104+
* This overload handles looping over the elements of a standard vector.
105+
*
106+
* @tparam T A standard vector with inner type inheriting from
107+
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
108+
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
109+
* @tparam Lp Scalar type for the lp argument. The scalar type of T should be
110+
* convertable to this.
111+
* @param y Vector of unconstrained partial correlations
112+
* @param K Dimensionality of returned correlation matrix
113+
* @param[in, out] lp log density accumulator o
114+
*/
115+
template <typename T, typename Lp, require_std_vector_t<T>* = nullptr,
116+
require_convertible_t<return_type_t<T>, Lp>* = nullptr>
117+
inline auto corr_matrix_constrain(const T& y, int K, Lp& lp) {
118+
return apply_vector_unary<T>::apply(
119+
y, [&lp, K](auto&& v) { return corr_matrix_constrain(v, K, lp); });
104120
}
105121

106122
/**
@@ -115,18 +131,23 @@ inline auto corr_matrix_constrain(const T& x, Eigen::Index k,
115131
*
116132
* @tparam Jacobian if `true`, increment log density accumulator with log
117133
* absolute Jacobian determinant of constraining transform
118-
* @tparam T A standard vector with inner type inheriting from
119-
* `Eigen::DenseBase` or a `var_value` with inner type inheriting from
120-
* `Eigen::DenseBase` with compile time dynamic rows and 1 column
121-
* @param y Vector of unconstrained partial correlations
122-
* @param K Dimensionality of returned correlation matrix
134+
* @tparam T A type inheriting from `Eigen::DenseBase` or a `var_value` with
135+
* inner type inheriting from `Eigen::DenseBase` with compile time dynamic rows
136+
* and 1 column or standard vector thereof
137+
* @tparam Lp A scalar type for the lp argument. The scalar type of T should be
138+
* convertable to this.
139+
* @param x Vector of unconstrained partial correlations
140+
* @param k Dimensionality of returned correlation matrix
123141
* @param[in,out] lp log density accumulator
124142
*/
125-
template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr>
126-
inline auto corr_matrix_constrain(const T& y, int K, return_type_t<T>& lp) {
127-
return apply_vector_unary<T>::apply(y, [&lp, K](auto&& v) {
128-
return corr_matrix_constrain<Jacobian>(v, K, lp);
129-
});
143+
template <bool Jacobian, typename T, typename Lp,
144+
require_convertible_t<return_type_t<T>, Lp>* = nullptr>
145+
inline auto corr_matrix_constrain(const T& x, Eigen::Index k, Lp& lp) {
146+
if constexpr (Jacobian) {
147+
return corr_matrix_constrain(x, k, lp);
148+
} else {
149+
return corr_matrix_constrain(x, k);
150+
}
130151
}
131152

132153
} // namespace math

0 commit comments

Comments
 (0)