1 #include <stan/math/prim.hpp>
2 #include <test/unit/util.hpp>
3 #include <gtest/gtest.h>
4 #include <cmath>
5 #include <limits>
6 #include <string>
7 #include <vector>
8
TEST(MathPrimMat,vec_double_gp_exponential_cov1)9 TEST(MathPrimMat, vec_double_gp_exponential_cov1) {
10 double sigma = 0.2;
11 double l = 5.0;
12
13 std::vector<double> x = {-2, -1, -0.5};
14
15 Eigen::MatrixXd cov;
16 cov = stan::math::gp_exponential_cov(x, sigma, l);
17 for (int i = 0; i < 3; i++)
18 for (int j = 0; j < 3; j++)
19 EXPECT_FLOAT_EQ(
20 sigma * sigma * std::exp(-stan::math::distance(x[i], x[j]) / (l)),
21 cov(i, j))
22 << "index: (" << i << ", " << j << ")";
23 }
24
TEST(MathPrimMat,vec_double_double_gp_exponential_cov1)25 TEST(MathPrimMat, vec_double_double_gp_exponential_cov1) {
26 double sigma = 0.2;
27 double l = 5.0;
28
29 std::vector<double> x1 = {-2, -1, -0.5};
30
31 std::vector<double> x2 = {3, -4, -0.7};
32
33 Eigen::MatrixXd cov;
34 cov = stan::math::gp_exponential_cov(x1, x2, sigma, l);
35 for (int i = 0; i < 3; i++)
36 for (int j = 0; j < 3; j++)
37 EXPECT_FLOAT_EQ(
38 sigma * sigma
39 * std::exp(-1.0 * stan::math::distance(x1[i], x2[j]) / l),
40 cov(i, j))
41 << "index: (" << i << ", " << j << ")";
42 }
43
TEST(MathPrimMat,vec_eigen_gp_exponential_cov1)44 TEST(MathPrimMat, vec_eigen_gp_exponential_cov1) {
45 double l = 0.2;
46 double sigma = 0.3;
47
48 std::vector<Eigen::Matrix<double, -1, 1>> x1(3);
49 for (size_t i = 0; i < x1.size(); ++i) {
50 x1[i].resize(3, 1);
51 x1[i] << 1 * i, 2 * i, 3 * i;
52 }
53
54 Eigen::MatrixXd cov;
55 cov = stan::math::gp_exponential_cov(x1, sigma, l);
56 for (int i = 0; i < 3; i++) {
57 for (int j = 0; j < 3; j++) {
58 EXPECT_FLOAT_EQ(
59 sigma * sigma
60 * std::exp(-1.0 * stan::math::distance(x1[i], x1[j]) / l),
61 cov(i, j))
62 << "index: (" << i << ", " << j << ")";
63 }
64 }
65 }
66
TEST(MathPrimMat,vec_eigen_vec_eigen_gp_exponential_cov1)67 TEST(MathPrimMat, vec_eigen_vec_eigen_gp_exponential_cov1) {
68 double sigma = 0.3;
69 double l = 2.3;
70
71 std::vector<Eigen::Matrix<double, -1, 1>> x1(3);
72 for (size_t i = 0; i < x1.size(); ++i) {
73 x1[i].resize(3, 1);
74 x1[i] << 1 * i, 2 * i, 3 * i;
75 }
76
77 std::vector<Eigen::Matrix<double, -1, 1>> x2(3);
78 for (size_t i = 0; i < x2.size(); ++i) {
79 x2[i].resize(3, 1);
80 x2[i] << 2 * i, 3 * i, 4 * i;
81 }
82
83 Eigen::MatrixXd cov;
84 cov = stan::math::gp_exponential_cov(x1, x2, sigma, l);
85 for (int i = 0; i < 3; i++) {
86 for (int j = 0; j < 3; j++) {
87 EXPECT_FLOAT_EQ(
88 sigma * sigma
89 * std::exp(-1.0
90 * stan::math::sqrt(
91 stan::math::squared_distance(x1[i], x2[j]))
92 / l),
93 cov(i, j))
94 << "index: (" << i << ", " << j << ")";
95 }
96 }
97
98 cov = stan::math::gp_exponential_cov(x2, x1, sigma, l);
99 for (int i = 0; i < 3; i++) {
100 for (int j = 0; j < 3; j++) {
101 EXPECT_FLOAT_EQ(
102 sigma * sigma
103 * std::exp(-1.0 * stan::math::distance(x2[i], x1[j]) / l),
104 cov(i, j))
105 << "index: (" << i << ", " << j << ")";
106 }
107 }
108 }
109
TEST(MathPrimMat,vec_eigen_vec_eigen_ard_gp_exponential_cov1)110 TEST(MathPrimMat, vec_eigen_vec_eigen_ard_gp_exponential_cov1) {
111 double sigma = 0.3;
112 double temp;
113
114 std::vector<double> l = {1, 2, 3};
115
116 std::vector<Eigen::Matrix<double, -1, 1>> x1(3);
117 for (size_t i = 0; i < x1.size(); ++i) {
118 x1[i].resize(3, 1);
119 x1[i] << 1 * i, 2 * i, 3 * i;
120 }
121
122 std::vector<Eigen::Matrix<double, -1, 1>> x2(3);
123 for (size_t i = 0; i < x2.size(); ++i) {
124 x2[i].resize(3, 1);
125 x2[i] << 2 * i, 3 * i, 4 * i;
126 }
127
128 Eigen::MatrixXd cov;
129 cov = stan::math::gp_exponential_cov(x1, x2, sigma, l);
130 std::vector<Eigen::Matrix<double, -1, 1>> x1_new
131 = stan::math::divide_columns(x1, l);
132 std::vector<Eigen::Matrix<double, -1, 1>> x2_new
133 = stan::math::divide_columns(x2, l);
134
135 for (int i = 0; i < 3; i++) {
136 for (int j = 0; j < 3; j++) {
137 EXPECT_FLOAT_EQ(
138 sigma * sigma
139 * std::exp(-1.0 * stan::math::distance(x1_new[i], x2_new[j])),
140 cov(i, j))
141 << "index: (" << i << ", " << j << ")";
142 }
143 }
144 cov = stan::math::gp_exponential_cov(x2, x1, sigma, l);
145 for (int i = 0; i < 3; i++) {
146 for (int j = 0; j < 3; j++) {
147 EXPECT_FLOAT_EQ(
148 sigma * sigma * std::exp(-stan::math::distance(x2_new[i], x1_new[j])),
149 cov(i, j))
150 << "index: (" << i << ", " << j << ")";
151 }
152 }
153 }
154
TEST(MathPrimMat,domain_err_training_sig_l_gp_exp_cov)155 TEST(MathPrimMat, domain_err_training_sig_l_gp_exp_cov) {
156 double sigma = .2;
157 double l = 7;
158 double sigma_bad = -1.0;
159 double l_bad = -1.0;
160
161 std::vector<double> l_vec = {1, 2, 3};
162
163 std::vector<double> l_vec_bad = {1, -1, 2};
164
165 std::vector<double> x = {-2, -1, -0.5};
166
167 std::vector<Eigen::Matrix<double, -1, 1>> x_2(3);
168 for (size_t i = 0; i < x_2.size(); ++i) {
169 x_2[i].resize(3, 1);
170 x_2[i] << 1, 2, 3;
171 }
172
173 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x, sigma, l_bad),
174 std::domain_error, " length scale");
175
176 EXPECT_THROW(stan::math::gp_exponential_cov(x, sigma_bad, l_bad),
177 std::domain_error);
178
179 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x, sigma, l_bad),
180 std::domain_error, " length scale");
181 EXPECT_THROW(stan::math::gp_exponential_cov(x, sigma_bad, l_bad),
182 std::domain_error);
183
184 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, sigma, l_bad),
185 std::domain_error, " length scale");
186 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, sigma_bad, l_bad),
187 std::domain_error);
188
189 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, sigma, l_vec_bad),
190 std::domain_error, " length scale");
191 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, sigma_bad, l_vec_bad),
192 std::domain_error);
193
194 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, sigma_bad, l_vec),
195 std::domain_error, " magnitude");
196
197 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, x_2, sigma, l_bad),
198 std::domain_error, " length scale");
199 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, x_2, sigma_bad, l),
200 std::domain_error, " magnitude");
201 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, x_2, sigma_bad, l_bad),
202 std::domain_error);
203
204 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, x_2, sigma, l_bad),
205 std::domain_error, " length scale");
206 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, x_2, sigma_bad, l_vec),
207 std::domain_error, " magnitude");
208 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, x_2, sigma_bad, l_bad),
209 std::domain_error);
210
211 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, x_2, sigma, l_vec_bad),
212 std::domain_error, " length scale");
213 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, x_2, sigma_bad, l_vec_bad),
214 std::domain_error);
215 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, x_2, sigma_bad, l_vec),
216 std::domain_error, " magnitude");
217
218 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, x_2, sigma, l_vec_bad),
219 std::domain_error, " length scale");
220 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, x_2, sigma_bad, l_vec_bad),
221 std::domain_error);
222 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, x_2, sigma_bad, l_vec),
223 std::domain_error, " magnitude");
224
225 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, x_2, sigma, l_vec_bad),
226 std::domain_error, " length scale");
227 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, x_2, sigma_bad, l_vec_bad),
228 std::domain_error);
229 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x_2, x_2, sigma_bad, l_vec),
230 std::domain_error, " magnitude");
231 }
232
TEST(MathPrimMat,nan_error_training_sig_l_gp_exp_cov)233 TEST(MathPrimMat, nan_error_training_sig_l_gp_exp_cov) {
234 double sigma = 0.2;
235 double l = 5;
236
237 std::vector<double> x = {-2, -1, -0.5};
238
239 std::vector<Eigen::Matrix<double, -1, 1>> x_2(3);
240 for (size_t i = 0; i < x_2.size(); ++i) {
241 x_2[i].resize(3, 1);
242 x_2[i] << 1, 2, 3;
243 }
244
245 std::vector<double> l_vec = {1, 2, 3};
246
247 std::vector<double> l_vec_bad = {1, -1, 3};
248
249 std::vector<double> x_bad(x);
250 x_bad[1] = std::numeric_limits<double>::quiet_NaN();
251
252 std::vector<Eigen::Matrix<double, -1, 1>> x_bad_2(x_2);
253 x_bad_2[1](1) = std::numeric_limits<double>::quiet_NaN();
254
255 double sigma_bad = std::numeric_limits<double>::quiet_NaN();
256 double l_bad = std::numeric_limits<double>::quiet_NaN();
257
258 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x, sigma, l_bad),
259 std::domain_error, " length scale");
260 EXPECT_THROW_MSG(stan::math::gp_exponential_cov(x, sigma_bad, l),
261 std::domain_error, " magnitude");
262 EXPECT_THROW(stan::math::gp_exponential_cov(x, sigma_bad, l_bad),
263 std::domain_error);
264
265 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad, l, sigma),
266 std::domain_error);
267 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad, sigma_bad, l),
268 std::domain_error);
269 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad, sigma_bad, l_bad),
270 std::domain_error);
271 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad, sigma, l_bad),
272 std::domain_error);
273
274 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, sigma, l_bad),
275 std::domain_error);
276 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, sigma_bad, l),
277 std::domain_error);
278 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, sigma_bad, l_bad),
279 std::domain_error);
280
281 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad_2, sigma, l_bad),
282 std::domain_error);
283 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad_2, sigma_bad, l),
284 std::domain_error);
285 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad_2, sigma, l),
286 std::domain_error);
287 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad_2, sigma_bad, l_bad),
288 std::domain_error);
289
290 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, sigma, l_vec_bad),
291 std::domain_error);
292 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, sigma_bad, l_vec),
293 std::domain_error);
294 EXPECT_THROW(stan::math::gp_exponential_cov(x_2, sigma_bad, l_vec_bad),
295 std::domain_error);
296
297 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad_2, sigma, l_vec_bad),
298 std::domain_error);
299 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad_2, sigma_bad, l_vec),
300 std::domain_error);
301 EXPECT_THROW(stan::math::gp_exponential_cov(x_bad_2, sigma_bad, l_vec_bad),
302 std::domain_error);
303 }
304
TEST(MathPrimMat,check_dim_mismatch_gp_exp_cov)305 TEST(MathPrimMat, check_dim_mismatch_gp_exp_cov) {
306 double sig = 1.0;
307 double l = 1.0;
308
309 std::vector<Eigen::Matrix<double, -1, 1>> x(2);
310 x[0].resize(2, 1);
311 x[0] << 1, 2;
312 x[1].resize(3, 1);
313 x[1] << 1, 2, 3;
314
315 EXPECT_THROW(stan::math::gp_exponential_cov(x, sig, l),
316 std::invalid_argument);
317
318 std::vector<Eigen::Matrix<double, -1, 1>> x1(2);
319 x1[0].resize(2, 1);
320 x1[0] << 1, 2;
321 x1[1].resize(2, 1);
322 x1[1] << 1, 2;
323
324 std::vector<Eigen::Matrix<double, -1, 1>> x2(3);
325 x2[0].resize(2, 1);
326 x2[0] << 1, 2;
327 x2[1].resize(3, 1);
328 x2[1] << 1, 2, 3;
329
330 EXPECT_THROW(stan::math::gp_exponential_cov(x1, x2, sig, l),
331 std::invalid_argument);
332 }
333
TEST(MathPrimMat,zero_size_gp_exp_cov)334 TEST(MathPrimMat, zero_size_gp_exp_cov) {
335 double sigma = 0.2;
336
337 std::vector<double> l(0);
338
339 std::vector<Eigen::Matrix<double, -1, 1>> x(0);
340
341 Eigen::MatrixXd cov;
342 Eigen::MatrixXd cov2;
343 EXPECT_NO_THROW(cov = stan::math::gp_exponential_cov(x, sigma, l));
344 EXPECT_NO_THROW(cov2 = stan::math::gp_exponential_cov(x, x, sigma, l));
345 EXPECT_EQ(0, cov.rows());
346 EXPECT_EQ(0, cov.cols());
347 }
348
TEST(MathPrimMat,calculations_gp_exp_cov)349 TEST(MathPrimMat, calculations_gp_exp_cov) {
350 double sigma = 1.0;
351 double l = 1.0;
352
353 std::vector<double> x1 = {1, 2};
354 std::vector<double> x2 = {2, 3};
355
356 Eigen::MatrixXd cov;
357 EXPECT_NO_THROW(cov = stan::math::gp_exponential_cov(x1, sigma, l));
358 ASSERT_FLOAT_EQ(1.0, cov(0, 0));
359 ASSERT_FLOAT_EQ(1.0, cov(1, 1));
360 ASSERT_FLOAT_EQ(exp(-1), cov(1, 0));
361 ASSERT_FLOAT_EQ(exp(-1), cov(0, 1));
362 EXPECT_NO_THROW(cov = stan::math::gp_exponential_cov(x1, x2, sigma, l));
363 ASSERT_FLOAT_EQ(exp(-1), cov(0, 0));
364 ASSERT_FLOAT_EQ(exp(-1), cov(1, 1));
365 ASSERT_FLOAT_EQ(1.0, cov(1, 0));
366 ASSERT_FLOAT_EQ(exp(-2), cov(0, 1));
367 }
368
TEST(MathPrimMat,calculations_ard_gp_exp_cov)369 TEST(MathPrimMat, calculations_ard_gp_exp_cov) {
370 double sigma = 1.0;
371
372 std::vector<double> l = {1, 2};
373
374 std::vector<Eigen::Matrix<double, -1, 1>> x(2);
375 x[0].resize(2, 1);
376 x[0] << 1, 1;
377 x[1].resize(2, 1);
378 x[1] << 2, 4;
379
380 Eigen::MatrixXd cov;
381 Eigen::MatrixXd cov2;
382 cov = stan::math::gp_exponential_cov(x, sigma, l);
383 cov2 = stan::math::gp_exponential_cov(x, x, sigma, l);
384
385 EXPECT_FLOAT_EQ(1.0, cov(0, 0));
386 EXPECT_FLOAT_EQ(1.0, cov2(0, 0));
387 EXPECT_FLOAT_EQ(exp(-sqrt(1 + 9.0 / 4.0)), cov(1, 0));
388 EXPECT_FLOAT_EQ(exp(-sqrt(1 + 9.0 / 4.0)), cov2(1, 0));
389 EXPECT_FLOAT_EQ(exp(-sqrt(1 + 9.0 / 4.0)), cov(0, 1));
390 EXPECT_FLOAT_EQ(exp(-sqrt(1 + 9.0 / 4.0)), cov2(0, 1));
391 EXPECT_FLOAT_EQ(1.0, cov(1, 1));
392 EXPECT_FLOAT_EQ(1.0, cov2(1, 1));
393 }
394