1 /*
2 SPDX-FileCopyrightText: 2014-2017 Max Planck Society.
3 All rights reserved.
4
5 SPDX-License-Identifier: BSD-3-Clause
6 */
7
8 /* Created by Edgar Klenske <edgar.klenske@tuebingen.mpg.de>
9 *
10 * Provides the test cases for the Gaussian Process functionality.
11 *
12 */
13
14 #include <gtest/gtest.h>
15 #include <vector>
16 #include <iostream>
17 #include "math_tools.h"
18 #include "gaussian_process.h"
19 #include "covariance_functions.h"
20
21 #include <fstream>
22
23 class GPTest : public ::testing::Test
24 {
25 public:
GPTest()26 GPTest(): random_vector_(11), location_vector_(11), hyper_parameters_(4), extra_parameters_(1)
27 {
28 random_vector_ << -0.1799, -1.4215, -0.2774, 2.6056, 0.6471, -0.4366,
29 1.3820, 0.4340, 0.8970, -0.7286, -1.7046;
30 location_vector_ << 0, 0.1000, 0.2000, 0.3000, 0.4000, 0.5000, 0.6000,
31 0.7000, 0.8000, 0.9000, 1.0000;
32 hyper_parameters_ << 1, 2, 1, 2;
33 extra_parameters_ << 5;
34
35 covariance_function_ = covariance_functions::PeriodicSquareExponential(hyper_parameters_);
36 covariance_function_.setExtraParameters(extra_parameters_);
37
38 gp_ = GP(covariance_function_);
39 }
40 GP gp_;
41 Eigen::VectorXd random_vector_;
42 Eigen::VectorXd location_vector_;
43 Eigen::VectorXd hyper_parameters_;
44 Eigen::VectorXd extra_parameters_;
45 covariance_functions::PeriodicSquareExponential covariance_function_;
46 };
47
48 // This test is based on Matlab computations
TEST_F(GPTest,drawSample_prior_test)49 TEST_F(GPTest, drawSample_prior_test)
50 {
51 Eigen::VectorXd sample = gp_.drawSample(location_vector_, random_vector_);
52 Eigen::VectorXd expected_sample(11);
53 expected_sample << -1.8799, -2.2659, -2.6541, -3.0406, -3.4214, -3.7926, -4.1503, -4.4907, -4.8101, -5.1052, -5.3726;
54 for (int i = 0; i < expected_sample.rows(); i++)
55 {
56 EXPECT_NEAR(sample(i), expected_sample(i), 2e-1);
57 }
58 }
59
60 // This test is based on statistical expectations (mean)
TEST_F(GPTest,drawSamples_prior_mean_test)61 TEST_F(GPTest, drawSamples_prior_mean_test)
62 {
63
64 hyper_parameters_ << 1, 1, 1, 1; // use of smaller hypers needs less samples
65
66 covariance_function_ =
67 covariance_functions::PeriodicSquareExponential(hyper_parameters_);
68 gp_ = GP(covariance_function_);
69
70 int N = 10000; // number of samples to draw
71 location_vector_ = Eigen::VectorXd(1);
72 location_vector_ << 1;
73 Eigen::MatrixXd sample_collection(location_vector_.rows(), N);
74 sample_collection.setZero(); // just in case
75
76 for (int i = 0; i < N; i++)
77 {
78 Eigen::VectorXd sample = gp_.drawSample(location_vector_);
79 sample_collection.col(i) = sample;
80 }
81 Eigen::VectorXd sample_mean = sample_collection.rowwise().mean();
82
83 for (int i = 0; i < sample_mean.rows(); i++)
84 {
85 EXPECT_NEAR(0, sample_mean(i), 1e-1);
86 }
87 }
88
89 // This test is based on statistical expectations (covariance)
TEST_F(GPTest,drawSamples_prior_covariance_test)90 TEST_F(GPTest, drawSamples_prior_covariance_test)
91 {
92
93 hyper_parameters_ << 1, 1, 1, 1; // use of smaller hypers needs less samples
94
95 covariance_function_ =
96 covariance_functions::PeriodicSquareExponential(hyper_parameters_);
97 gp_ = GP(covariance_function_);
98
99 int N = 20000; // number of samples to draw
100 location_vector_ = Eigen::VectorXd(1);
101 location_vector_ << 1;
102 Eigen::MatrixXd sample_collection(location_vector_.rows(), N);
103 sample_collection.setZero(); // just in case
104
105 for (int i = 0; i < N; i++)
106 {
107 Eigen::VectorXd sample = gp_.drawSample(location_vector_);
108 sample_collection.col(i) = sample;
109 }
110 Eigen::MatrixXd sample_cov = sample_collection * sample_collection.transpose() / N;
111
112 Eigen::MatrixXd expected_cov = covariance_function_.evaluate(location_vector_, location_vector_);
113
114 for (int i = 0; i < sample_cov.rows(); i++)
115 {
116 for (int k = 0; k < sample_cov.cols(); k++)
117 {
118 EXPECT_NEAR(expected_cov(i, k), sample_cov(i, k), 2e-1);
119 }
120 }
121 }
122
123
TEST_F(GPTest,setCovarianceFunction)124 TEST_F(GPTest, setCovarianceFunction)
125 {
126 Eigen::VectorXd hyperparams(6);
127 hyperparams << 0.1, 15, 25, 15, 5000, 700;
128
129 GP instance_gp;
130 EXPECT_TRUE(instance_gp.setCovarianceFunction(covariance_functions::PeriodicSquareExponential(hyperparams.segment(1,4))));
131
132 GP instance_gp2 = GP(covariance_functions::PeriodicSquareExponential(Eigen::VectorXd::Zero(4)));
133 instance_gp2.setHyperParameters(hyperparams);
134
135 for (int i = 1; i < 5; i++) // first element is different (non set)
136 {
137 EXPECT_NEAR(instance_gp.getHyperParameters()[i], instance_gp2.getHyperParameters()[i], 1e-8);
138 }
139 }
140
141
142
TEST_F(GPTest,setCovarianceFunction_notworking_after_inference)143 TEST_F(GPTest, setCovarianceFunction_notworking_after_inference)
144 {
145 Eigen::VectorXd hyperparams(5);
146 hyperparams << 0.1, 15, 700, 25, 5000;
147
148 GP instance_gp;
149 EXPECT_TRUE(instance_gp.setCovarianceFunction(covariance_functions::PeriodicSquareExponential(hyperparams.tail(4))));
150
151 int N = 250;
152 Eigen::VectorXd location =
153 400 * math_tools::generate_uniform_random_matrix_0_1(N, 1)
154 - 200 * Eigen::MatrixXd::Ones(N, 1);
155
156 Eigen::VectorXd output_from_converged_hyperparams = instance_gp.drawSample(location);
157
158 instance_gp.infer(location, output_from_converged_hyperparams);
159 EXPECT_FALSE(instance_gp.setCovarianceFunction(covariance_functions::PeriodicSquareExponential(hyperparams.tail(4))));
160 }
161
162 // to be moved to some other file
TEST_F(GPTest,periodic_covariance_function_test)163 TEST_F(GPTest, periodic_covariance_function_test)
164 {
165 covariance_functions::PeriodicSquareExponential u;
166 EXPECT_EQ(u.getParameterCount(), 4);
167
168 GP instance_gp = GP(covariance_functions::PeriodicSquareExponential());
169 ASSERT_EQ(instance_gp.getHyperParameters().size(), 6);
170 instance_gp.setHyperParameters(Eigen::VectorXd::Zero(6)); // should not assert
171 }
172
173
174
TEST_F(GPTest,infer_prediction_clear_test)175 TEST_F(GPTest, infer_prediction_clear_test)
176 {
177 Eigen::VectorXd data_loc(1);
178 data_loc << 1;
179 Eigen::VectorXd data_out(1);
180 data_out << 1;
181 gp_.infer(data_loc, data_out);
182
183 Eigen::VectorXd prediction_location(2);
184 prediction_location << 1, 2;
185
186 Eigen::VectorXd prediction = gp_.predict(prediction_location);
187
188 EXPECT_NEAR(prediction(0), 1, 1e-6);
189 EXPECT_FALSE(std::abs(prediction(1) - 1) < 1e-6);
190
191 gp_.clearData();
192
193 prediction = gp_.predict(prediction_location);
194
195 EXPECT_NEAR(prediction(0), 0, 1e-6);
196 EXPECT_NEAR(prediction(1), 0, 1e-6);
197 }
198
TEST_F(GPTest,squareDistanceTest)199 TEST_F(GPTest, squareDistanceTest)
200 {
201 Eigen::MatrixXd a(4, 3);
202 Eigen::MatrixXd b(4, 5);
203 Eigen::MatrixXd c(3, 4);
204
205 a << 3, 5, 5, 4, 6, 6, 3, 2, 3, 1, 0, 3;
206
207 b << 1, 4, 5, 6, 7, 3, 4, 5, 6, 7, 0, 2, 4, 20, 2, 2, 3, -2, -2, 2;
208
209 c << 1, 2, 3, 4, 4, 5, 6, 7, 6, 7, 8, 9;
210
211
212 Eigen::MatrixXd sqdistc(4, 4);
213 Eigen::MatrixXd sqdistab(3, 5);
214
215 // Computed by Matlab
216 sqdistc << 0, 3, 12, 27, 3, 0, 3, 12, 12, 3, 0, 3, 27, 12, 3, 0;
217
218 sqdistab << 15, 6, 15, 311, 27, 33, 14, 9, 329, 9, 35, 6, 27, 315, 7;
219
220 // Test argument order
221 EXPECT_EQ(math_tools::squareDistance(a, b), math_tools::squareDistance
222 (b, a).transpose());
223
224 // Test that two identical Matrices give the same result
225 // (wether they are the same object or not)
226 EXPECT_EQ(math_tools::squareDistance(a, Eigen::MatrixXd(a)),
227 math_tools::squareDistance(a, a));
228 EXPECT_EQ(math_tools::squareDistance(a), math_tools::squareDistance(
229 a, a));
230
231 // Test that the implementation gives the same result as the Matlab
232 // implementation
233 EXPECT_EQ(math_tools::squareDistance(c, c), sqdistc);
234 EXPECT_EQ(math_tools::squareDistance(a, b), sqdistab);
235 }
236
TEST_F(GPTest,CovarianceTest2)237 TEST_F(GPTest, CovarianceTest2)
238 {
239 Eigen::VectorXd hyperParams(4), extraParams(1);
240 hyperParams << 1, 2, 1, 2;
241 hyperParams = hyperParams.array().log();
242 extraParams << 500;
243 extraParams = extraParams.array().log();
244
245 Eigen::VectorXd locations(5), X(3), Y(3);
246 locations << 0, 50, 100, 150, 200;
247 X << 0, 100, 200;
248 Y << 1, -1, 1;
249
250 covariance_functions::PeriodicSquareExponential covFunc(hyperParams);
251 covFunc.setExtraParameters(extraParams);
252 GP gp(covFunc);
253
254 Eigen::MatrixXd kxx_matlab(5, 5);
255 Eigen::MatrixXd kxX_matlab(5, 3);
256 Eigen::MatrixXd kXX_matlab(3, 3);
257
258 kxx_matlab << 8.0000, 3.3046, 2.0043, 1.0803, 0.6553,
259 3.3046, 8.0000, 3.3046, 2.0043, 1.0803,
260 2.0043, 3.3046, 8.0000, 3.3046, 2.0043,
261 1.0803, 2.0043, 3.3046, 8.0000, 3.3046,
262 0.6553, 1.0803, 2.0043, 3.3046, 8.0000;
263
264 kxX_matlab << 8.0000, 2.0043, 0.6553, 3.3046, 3.3046, 1.0803, 2.0043, 8.0000, 2.0043, 1.0803, 3.3046, 3.3046, 0.6553, 2.0043, 8.0000;
265
266 kXX_matlab << 8.0000, 2.0043, 0.6553, 2.0043, 8.0000, 2.0043, 0.6553, 2.0043, 8.0000;
267
268 Eigen::MatrixXd kxx = covFunc.evaluate(locations, locations);
269 Eigen::MatrixXd kxX = covFunc.evaluate(locations, X);
270 Eigen::MatrixXd kXX = covFunc.evaluate(X, X);
271
272 for (int col = 0; col < kxx.cols(); col++)
273 {
274 for (int row = 0; row < kxx.rows(); row++)
275 {
276 EXPECT_NEAR(kxx(row, col), kxx_matlab(row, col), 0.003);
277 }
278 }
279
280 for (int col = 0; col < kxX.cols(); col++)
281 {
282 for (int row = 0; row < kxX.rows(); row++)
283 {
284 EXPECT_NEAR(kxX(row, col), kxX_matlab(row, col), 0.003);
285 }
286 }
287
288 for (int col = 0; col < kXX.cols(); col++)
289 {
290 for (int row = 0; row < kXX.rows(); row++)
291 {
292 EXPECT_NEAR(kXX(row, col), kXX_matlab(row, col), 0.003);
293 }
294 }
295 }
296
TEST_F(GPTest,CovarianceTest3)297 TEST_F(GPTest, CovarianceTest3)
298 {
299 Eigen::Matrix<double, 6, 1> hyperParams;
300 hyperParams << 10, 1, 1, 1, 100, 1;
301 hyperParams = hyperParams.array().log();
302
303 Eigen::VectorXd periodLength(1);
304 periodLength << std::log(80);
305
306 Eigen::VectorXd locations(5), X(3), Y(3);
307 locations << 0, 50, 100, 150, 200;
308 X << 0, 100, 200;
309
310 covariance_functions::PeriodicSquareExponential2 covFunc(hyperParams);
311 covFunc.setExtraParameters(periodLength);
312
313 Eigen::MatrixXd kxx_matlab(5, 5);
314 kxx_matlab <<
315 3.00000, 1.06389, 0.97441, 1.07075, 0.27067,
316 1.06389, 3.00000, 1.06389, 0.97441, 1.07075,
317 0.97441, 1.06389, 3.00000, 1.06389, 0.97441,
318 1.07075, 0.97441, 1.06389, 3.00000, 1.06389,
319 0.27067, 1.07075, 0.97441, 1.06389, 3.00000;
320
321 Eigen::MatrixXd kxX_matlab(5, 3);
322 kxX_matlab <<
323 3.00000, 0.97441, 0.27067, 1.06389, 1.06389,
324 1.07075, 0.97441, 3.00000, 0.97441, 1.07075,
325 1.06389, 1.06389, 0.27067, 0.97441, 3.00000;
326
327 Eigen::MatrixXd kXX_matlab(3, 3);
328 kXX_matlab <<
329 3.00000, 0.97441, 0.27067,
330 0.97441, 3.00000, 0.97441,
331 0.27067, 0.97441, 3.00000;
332
333 Eigen::MatrixXd kxx = covFunc.evaluate(locations, locations);
334 Eigen::MatrixXd kxX = covFunc.evaluate(locations, X);
335 Eigen::MatrixXd kXX = covFunc.evaluate(X, X);
336
337 for (int col = 0; col < kxx.cols(); col++)
338 {
339 for (int row = 0; row < kxx.rows(); row++)
340 {
341 EXPECT_NEAR(kxx(row, col), kxx_matlab(row, col), 0.01);
342 }
343 }
344
345 for (int col = 0; col < kxX.cols(); col++)
346 {
347 for (int row = 0; row < kxX.rows(); row++)
348 {
349 EXPECT_NEAR(kxX(row, col), kxX_matlab(row, col), 0.01);
350 }
351 }
352
353 for (int col = 0; col < kXX.cols(); col++)
354 {
355 for (int row = 0; row < kXX.rows(); row++)
356 {
357 EXPECT_NEAR(kXX(row, col), kXX_matlab(row, col), 0.01);
358 }
359 }
360 }
361
main(int argc,char ** argv)362 int main(int argc, char** argv)
363 {
364 ::testing::InitGoogleTest(&argc, argv);
365 return RUN_ALL_TESTS();
366 }
367