Automatic Differentiation
 
Loading...
Searching...
No Matches
gp_matern52_cov.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_GP_MATERN52_COV_HPP
2#define STAN_MATH_PRIM_FUN_GP_MATERN52_COV_HPP
3
12#include <cmath>
13#include <vector>
14
15namespace stan {
16namespace math {
17
38template <typename T_x, typename T_s, typename T_l>
39inline typename Eigen::Matrix<return_type_t<T_x, T_s, T_l>, Eigen::Dynamic,
40 Eigen::Dynamic>
41gp_matern52_cov(const std::vector<T_x> &x, const T_s &sigma,
42 const T_l &length_scale) {
43 using std::exp;
44 using std::sqrt;
45
46 size_t x_size = x.size();
47 Eigen::Matrix<return_type_t<T_x, T_s, T_l>, Eigen::Dynamic, Eigen::Dynamic>
48 cov(x_size, x_size);
49
50 if (x_size == 0) {
51 return cov;
52 }
53
54 for (size_t n = 0; n < x_size; ++n) {
55 check_not_nan("gp_matern52_cov", "x", x[n]);
56 }
57
58 check_positive_finite("gp_matern52_cov", "magnitude", sigma);
59 check_positive_finite("gp_matern52_cov", "length scale", length_scale);
60
61 T_s sigma_sq = square(sigma);
62 T_l root_5_inv_l = sqrt(5.0) / length_scale;
63 T_l inv_l_sq_5_3 = 5.0 / (3.0 * square(length_scale));
64
65 size_t block_size = 10;
66 for (size_t jb = 0; jb < x_size; jb += block_size) {
67 for (size_t ib = jb; ib < x_size; ib += block_size) {
68 size_t j_end = std::min(x_size, jb + block_size);
69 for (size_t j = jb; j < j_end; ++j) {
70 cov.coeffRef(j, j) = sigma_sq;
71 size_t i_end = std::min(x_size, ib + block_size);
72 for (size_t i = std::max(ib, j + 1); i < i_end; ++i) {
73 return_type_t<T_x> sq_distance = squared_distance(x[i], x[j]);
74 return_type_t<T_x> dist = sqrt(sq_distance);
75 cov.coeffRef(j, i) = cov.coeffRef(i, j)
76 = sigma_sq
77 * (1.0 + root_5_inv_l * dist + inv_l_sq_5_3 * sq_distance)
78 * exp(-root_5_inv_l * dist);
79 }
80 }
81 }
82 }
83 return cov;
84}
85
109template <typename T_x, typename T_s, typename T_l>
110inline typename Eigen::Matrix<return_type_t<T_x, T_s, T_l>, Eigen::Dynamic,
111 Eigen::Dynamic>
112gp_matern52_cov(const std::vector<Eigen::Matrix<T_x, Eigen::Dynamic, 1>> &x,
113 const T_s &sigma, const std::vector<T_l> &length_scale) {
114 using std::exp;
115
116 size_t x_size = x.size();
117 Eigen::Matrix<return_type_t<T_x, T_s, T_l>, Eigen::Dynamic, Eigen::Dynamic>
118 cov(x_size, x_size);
119 if (x_size == 0) {
120 return cov;
121 }
122
123 size_t l_size = length_scale.size();
124 for (size_t n = 0; n < x_size; ++n) {
125 check_not_nan("gp_matern52_cov", "x", x[n]);
126 }
127
128 check_positive_finite("gp_matern52_cov", "magnitude", sigma);
129 check_positive_finite("gp_matern52_cov", "length scale", length_scale);
130 check_size_match("gp_matern52_cov", "x dimension", x[0].size(),
131 "number of length scales", l_size);
132
133 T_s sigma_sq = square(sigma);
134 double root_5 = sqrt(5.0);
135 double five_thirds = 5.0 / 3.0;
136
137 std::vector<Eigen::Matrix<return_type_t<T_x, T_l>, -1, 1>> x_new
138 = divide_columns(x, length_scale);
139 size_t block_size = 10;
140
141 for (size_t jb = 0; jb < x_size; jb += block_size) {
142 for (size_t ib = jb; ib < x_size; ib += block_size) {
143 size_t j_end = std::min(x_size, jb + block_size);
144 for (size_t j = jb; j < j_end; ++j) {
145 cov.coeffRef(j, j) = sigma_sq;
146 size_t i_end = std::min(x_size, ib + block_size);
147 for (size_t i = std::max(ib, j + 1); i < i_end; ++i) {
148 return_type_t<T_x, T_l> sq_distance
149 = squared_distance(x_new[i], x_new[j]);
150 return_type_t<T_x, T_l> dist = sqrt(sq_distance);
151 cov.coeffRef(j, i) = cov.coeffRef(i, j)
152 = sigma_sq * (1.0 + root_5 * dist + five_thirds * sq_distance)
153 * exp(-root_5 * dist);
154 }
155 }
156 }
157 }
158 return cov;
159}
160
185template <typename T_x1, typename T_x2, typename T_s, typename T_l>
186inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>,
187 Eigen::Dynamic, Eigen::Dynamic>
188gp_matern52_cov(const std::vector<T_x1> &x1, const std::vector<T_x2> &x2,
189 const T_s &sigma, const T_l &length_scale) {
190 using std::exp;
191
192 size_t x1_size = x1.size();
193 size_t x2_size = x2.size();
194 Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>, Eigen::Dynamic,
195 Eigen::Dynamic>
196 cov(x1_size, x2_size);
197
198 if (x1_size == 0 || x2_size == 0) {
199 return cov;
200 }
201
202 for (size_t n = 0; n < x1_size; ++n) {
203 check_not_nan("gp_matern52_cov", "x1", x1[n]);
204 }
205 for (size_t n = 0; n < x2_size; ++n) {
206 check_not_nan("gp_matern52_cov", "x1", x2[n]);
207 }
208
209 check_positive_finite("gp_matern52_cov", "magnitude", sigma);
210 check_positive_finite("gp_matern52_cov", "length scale", length_scale);
211
212 T_s sigma_sq = square(sigma);
213 T_l root_5_inv_l = sqrt(5.0) / length_scale;
214 T_l inv_l_sq_5_3 = 5.0 / (3.0 * square(length_scale));
215
216 size_t block_size = 10;
217
218 for (size_t ib = 0; ib < x1_size; ib += block_size) {
219 for (size_t jb = 0; jb < x2_size; jb += block_size) {
220 size_t j_end = std::min(x2_size, jb + block_size);
221 for (size_t j = jb; j < j_end; ++j) {
222 size_t i_end = std::min(x1_size, ib + block_size);
223 for (size_t i = ib; i < i_end; ++i) {
224 return_type_t<T_x1, T_x2> sq_distance
225 = squared_distance(x1[i], x2[j]);
226 return_type_t<T_x1, T_x2> dist = sqrt(sq_distance);
227 cov.coeffRef(i, j)
228 = sigma_sq
229 * (1.0 + root_5_inv_l * dist + inv_l_sq_5_3 * sq_distance)
230 * exp(-root_5_inv_l * dist);
231 }
232 }
233 }
234 }
235 return cov;
236}
237
267template <typename T_x1, typename T_x2, typename T_s, typename T_l>
268inline typename Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>,
269 Eigen::Dynamic, Eigen::Dynamic>
270gp_matern52_cov(const std::vector<Eigen::Matrix<T_x1, Eigen::Dynamic, 1>> &x1,
271 const std::vector<Eigen::Matrix<T_x2, Eigen::Dynamic, 1>> &x2,
272 const T_s &sigma, const std::vector<T_l> &length_scale) {
273 using std::exp;
274
275 size_t x1_size = x1.size();
276 size_t x2_size = x2.size();
277 Eigen::Matrix<return_type_t<T_x1, T_x2, T_s, T_l>, Eigen::Dynamic,
278 Eigen::Dynamic>
279 cov(x1_size, x2_size);
280
281 if (x1_size == 0 || x2_size == 0) {
282 return cov;
283 }
284
285 size_t l_size = length_scale.size();
286
287 for (size_t n = 0; n < x1_size; ++n) {
288 check_not_nan("gp_matern52_cov", "x1", x1[n]);
289 }
290 for (size_t n = 0; n < x2_size; ++n) {
291 check_not_nan("gp_matern52_cov", "x1", x2[n]);
292 }
293
294 check_positive_finite("gp_matern52_cov", "magnitude", sigma);
295 check_positive_finite("gp_matern52_cov", "length scale", length_scale);
296 check_size_match("gp_matern52_cov", "x dimension", x1[0].size(),
297 "number of length scales", l_size);
298 check_size_match("gp_matern52_cov", "x dimension", x2[0].size(),
299 "number of length scales", l_size);
300
301 T_s sigma_sq = square(sigma);
302 double root_5 = sqrt(5.0);
303 double five_thirds = 5.0 / 3.0;
304
305 std::vector<Eigen::Matrix<return_type_t<T_x1, T_l>, -1, 1>> x1_new
306 = divide_columns(x1, length_scale);
307 std::vector<Eigen::Matrix<return_type_t<T_x2, T_l>, -1, 1>> x2_new
308 = divide_columns(x2, length_scale);
309 size_t block_size = 10;
310
311 for (size_t ib = 0; ib < x1_size; ib += block_size) {
312 for (size_t jb = 0; jb < x2_size; jb += block_size) {
313 size_t j_end = std::min(x2_size, jb + block_size);
314 for (size_t j = jb; j < j_end; ++j) {
315 size_t i_end = std::min(x1_size, ib + block_size);
316 for (size_t i = ib; i < i_end; ++i) {
318 = squared_distance(x1_new[i], x2_new[j]);
319 return_type_t<T_x1, T_x2, T_l> dist = sqrt(sq_distance);
320 cov.coeffRef(i, j)
321 = sigma_sq * (1.0 + root_5 * dist + five_thirds * sq_distance)
322 * exp(-root_5 * dist);
323 }
324 }
325 }
326 }
327 return cov;
328}
329} // namespace math
330} // namespace stan
331#endif
void divide_columns(matrix_cl< T1 > &A, const matrix_cl< T2 > &B)
Divides each column of a matrix by a vector.
matrix_cl< return_type_t< T1, T2, T3 > > gp_matern52_cov(const T1 &x, const T2 sigma, const T3 length_scale)
Matern 5/2 kernel on the GPU.
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
int64_t size(const T &m)
Returns the size (number of the elements) of a matrix_cl or var_value<matrix_cl<T>>.
Definition size.hpp:19
fvar< T > sqrt(const fvar< T > &x)
Definition sqrt.hpp:17
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
auto squared_distance(const T_a &a, const T_b &b)
Returns the squared distance.
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.
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
fvar< T > square(const fvar< T > &x)
Definition square.hpp:12
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 ...