1 #include "MUQ/Approximation/GaussianProcesses/CovarianceKernels.h"
2 #include "MUQ/Approximation/GaussianProcesses/StateSpaceGP.h"
3 #include "MUQ/Approximation/GaussianProcesses/GaussianProcess.h"
4
5 #include "MUQ/Utilities/Exceptions.h"
6
7 #include <gtest/gtest.h>
8
9 #include <memory>
10 #include <random>
11 #include <iostream>
12
13
14 using namespace muq::Approximation;
15
16
TEST(Approximation_GP,MaternStateSpace)17 TEST(Approximation_GP, MaternStateSpace)
18 {
19
20 const double sigma2 = 1.0;
21 const double length = 0.15;
22
23 const double nu = 3.0/2.0;
24
25 MaternKernel kernel(1, sigma2, length, nu);
26
27 ZeroMean mu(1,1);
28 StateSpaceGP gp(mu, kernel);
29
30 EXPECT_EQ(nu+0.5, gp.stateDim);
31
32 // draw a random sample from the SDE model
33 Eigen::VectorXd obsTimes = Eigen::VectorXd::LinSpaced(100, 0, 1);
34 Eigen::MatrixXd realization = gp.Sample(obsTimes);
35
36 }
37
TEST(Approximation_GP,StateSpace_DistributionIntegration)38 TEST(Approximation_GP, StateSpace_DistributionIntegration)
39 {
40
41 const double sigma2 = 1.0;
42 const double length = 0.15;
43
44 const double nu = 3.0/2.0;
45
46 MaternKernel kernel(1, sigma2, length, nu);
47
48 std::shared_ptr<muq::Modeling::LinearSDE> sde; // The underying SDE
49 std::shared_ptr<muq::Modeling::LinearOperator> H; // Observation matrix to go from SDE state to actual GP variable
50 Eigen::MatrixXd pinf; // Steady state covariance matrix
51
52 std::tie(sde, H, pinf) = kernel.GetStateSpace();
53
54 Eigen::MatrixXd p0 = Eigen::MatrixXd::Identity(sde->stateDim, sde->stateDim);
55 Eigen::MatrixXd mu0 = Eigen::VectorXd::Ones(sde->stateDim);
56
57 // Integrate the SDE for a while
58 Eigen::MatrixXd muT, pT;
59 std::tie(muT, pT) = sde->EvolveDistribution(mu0,p0,30);
60
61 EXPECT_NEAR(0.0, muT(0), 1e-14);
62 EXPECT_NEAR(0.0, muT(1), 1e-14);
63
64 EXPECT_NEAR(pinf(0,0), pT(0,0), 1e-11);
65 EXPECT_NEAR(pinf(0,1), pT(0,1), 1e-11);
66 EXPECT_NEAR(pinf(1,0), pT(1,0), 1e-11);
67 EXPECT_NEAR(pinf(1,1), pT(1,1), 1e-11);
68 }
69
TEST(Approximation_GP,StateSpace_DistributionIntegration2)70 TEST(Approximation_GP, StateSpace_DistributionIntegration2)
71 {
72
73 const double sigma2 = 1.0;
74 const double length = 0.15;
75
76 const double nu = 3.0/2.0;
77
78 MaternKernel kernel(1, sigma2, length, nu);
79
80 std::shared_ptr<muq::Modeling::LinearSDE> sde; // The underying SDE
81 std::shared_ptr<muq::Modeling::LinearOperator> H; // Observation matrix to go from SDE state to actual GP variable
82 Eigen::MatrixXd pinf; // Steady state covariance matrix
83
84 std::tie(sde, H, pinf) = kernel.GetStateSpace();
85
86 Eigen::MatrixXd p0 = Eigen::MatrixXd::Identity(sde->stateDim, sde->stateDim);
87 Eigen::VectorXd mu0 = Eigen::VectorXd::Ones(sde->stateDim);
88 std::pair<Eigen::VectorXd, Eigen::MatrixXd> dist0 = std::make_pair(mu0,p0);
89
90 // Integrate the SDE for a while
91 Eigen::MatrixXd muT, pT;
92 std::tie(muT, pT) = sde->EvolveDistribution(dist0,30);
93
94 EXPECT_NEAR(0.0, muT(0), 1e-14);
95 EXPECT_NEAR(0.0, muT(1), 1e-14);
96
97 EXPECT_NEAR(pinf(0,0), pT(0,0), 1e-11);
98 EXPECT_NEAR(pinf(0,1), pT(0,1), 1e-11);
99 EXPECT_NEAR(pinf(1,0), pT(1,0), 1e-11);
100 EXPECT_NEAR(pinf(1,1), pT(1,1), 1e-11);
101 }
102
TEST(Approximation_GP,PeriodicStateSpace)103 TEST(Approximation_GP, PeriodicStateSpace)
104 {
105
106 const double sigma2 = 1.0;
107 const double length = 0.6;
108 const double period = 0.25;
109 const double periodN = 50; // how many steps per period
110
111 PeriodicKernel kernel(1, sigma2, length, period);
112
113 boost::property_tree::ptree options;
114 options.put("PeriodicKernel.StateSpace.NumTerms",8);
115 options.put("SDE.dt", 5e-5);
116
117 ZeroMean mu(1,1);
118 StateSpaceGP gp(mu, kernel, options);
119
120 // draw a random sample from the SDE model
121 Eigen::VectorXd obsTimes = Eigen::VectorXd::LinSpaced(5*periodN+1, 0, 5*period);
122
123 Eigen::MatrixXd realization = gp.Sample(obsTimes);
124
125 // Make sure the sample is periodic
126 for(int i=0; i<obsTimes.size()-periodN-1; ++i)
127 EXPECT_NEAR(realization(0,i), realization(0,i+periodN), 1e-1);
128
129 }
130
TEST(Approximation_GP,SumStateSpace)131 TEST(Approximation_GP, SumStateSpace)
132 {
133
134 MaternKernel kernel1(1, 3.0, 1.0, 7.0/2.0);
135 MaternKernel kernel2(1, 0.15, 0.25, 1.0/2.0);
136
137 auto kernel = kernel1 + kernel2;
138
139 boost::property_tree::ptree options;
140 options.put("PeriodicKernel.StateSpace.NumTerms",8);
141 options.put("SDE.dt", 5e-5);
142
143 ZeroMean mu(1,1);
144 StateSpaceGP gp(mu, kernel, options);
145
146 // draw a random sample from the SDE model
147 Eigen::VectorXd obsTimes = Eigen::VectorXd::LinSpaced(100, 0, 4);
148
149 Eigen::MatrixXd realization = gp.Sample(obsTimes);
150 }
151
TEST(Approximation_GP,ProductStateSpace)152 TEST(Approximation_GP, ProductStateSpace)
153 {
154
155 const double sigma2 = 1.0;
156 const double length = 0.6;
157 const double nu = 3.0/2.0;
158 const double period = 0.25;
159 const double periodN = 50; // how many steps per period
160
161 PeriodicKernel kernel1(1, sigma2, 0.8, period);
162 MaternKernel kernel2(1, sigma2, 2.0, nu);
163
164 auto kernel12 = kernel1*kernel2;
165 auto kernel21 = kernel2*kernel1;
166 auto kernel22 = kernel2*kernel2;
167
168 boost::property_tree::ptree options;
169 options.put("PeriodicKernel.StateSpace.NumTerms",7);
170 options.put("SDE.dt", 1e-4);
171
172 ZeroMean mu(1,1);
173 StateSpaceGP gp1(mu, kernel1, options);
174 StateSpaceGP gp2(mu, kernel2, options);
175 EXPECT_EQ(int(nu+0.5), gp2.stateDim);
176
177 StateSpaceGP gp12(mu, kernel12, options);
178 StateSpaceGP gp21(mu, kernel21, options);
179
180 EXPECT_THROW(kernel22.GetStateSpace(options), muq::NotImplementedError);
181 EXPECT_EQ(gp1.stateDim *gp2.stateDim, gp12.stateDim);
182
183
184 // draw a random sample from the SDE model
185 Eigen::VectorXd obsTimes = Eigen::VectorXd::LinSpaced(5*periodN+1, 0, 5*period);
186
187 Eigen::MatrixXd realization1 = gp1.Sample(obsTimes);
188 Eigen::MatrixXd realization2 = gp2.Sample(obsTimes);
189 Eigen::MatrixXd realization12 = gp12.Sample(obsTimes);
190
191 }
192
TEST(Approximation_GP,StateSpacePredict_Interpolation)193 TEST(Approximation_GP, StateSpacePredict_Interpolation)
194 {
195
196 const double sigma2 = 1.0;
197 const double length = 0.3;
198 const double nu = 3.0/2.0;
199
200 MaternKernel kernel(1, sigma2, length, nu);
201
202 boost::property_tree::ptree options;
203 options.put("SDE.dt", 1e-5);
204
205 ZeroMean mu(1,1);
206 StateSpaceGP gp1(mu, kernel, options);
207 GaussianProcess gp2(mu, kernel);
208
209 const int numEvals = 100;
210 Eigen::MatrixXd evalPts(1,numEvals);
211 evalPts.row(0) = Eigen::VectorXd::LinSpaced(numEvals, 0, 2);
212
213 // Make a prediction about the prior
214 Eigen::MatrixXd predMu, predCov;
215 std::tie(predMu, predCov) = gp1.Predict(evalPts, GaussianProcess::DiagonalCov);
216
217 Eigen::MatrixXd predMu2, predCov2;
218 std::tie(predMu2, predCov2) = gp2.Predict(evalPts, GaussianProcess::DiagonalCov);
219
220 // Condition both GPs with some data
221 Eigen::MatrixXd obsLoc = 0.5*Eigen::MatrixXd::Ones(1,1);
222 Eigen::MatrixXd obsData = 0.25*Eigen::MatrixXd::Ones(1,1);
223 const double obsVar = 1e-2;
224
225 gp1.Condition(obsLoc, obsData, obsVar);
226 gp2.Condition(obsLoc, obsData, obsVar);
227
228 // add another observation
229 obsLoc(0) = 1.5;
230 gp1.Condition(obsLoc, obsData, obsVar);
231 gp2.Condition(obsLoc, obsData, obsVar);
232
233 // Make a posterior prediction
234 std::tie(predMu, predCov) = gp1.Predict(evalPts, GaussianProcess::DiagonalCov);
235 std::tie(predMu2, predCov2) = gp2.Predict(evalPts, GaussianProcess::DiagonalCov);
236
237 const double meanTol = 1e-2;
238 const double covTol = 1e-2;
239 for(int j=0; j<predMu.size(); ++j){
240 EXPECT_NEAR(predMu2(0,j), predMu(0,j), meanTol);
241 EXPECT_NEAR(predCov2(0,j), predCov(0,j), covTol);
242 }
243 }
244
TEST(Approximation_GP,StateSpacePredict_ExtrapolateRight)245 TEST(Approximation_GP, StateSpacePredict_ExtrapolateRight)
246 {
247
248 const double sigma2 = 1.0;
249 const double length = 0.3;
250 const double nu = 3.0/2.0;
251
252 MaternKernel kernel(1, sigma2, length, nu);
253
254 boost::property_tree::ptree options;
255 options.put("SDE.dt", 1e-5);
256
257 ZeroMean mu(1,1);
258 StateSpaceGP gp1(mu, kernel, options);
259 GaussianProcess gp2(mu, kernel);
260
261 const int numEvals = 100;
262 Eigen::MatrixXd evalPts(1,numEvals);
263 evalPts.row(0) = Eigen::VectorXd::LinSpaced(numEvals, 1.6, 3);
264
265 // Make a prediction about the prior
266 Eigen::MatrixXd predMu, predCov;
267 std::tie(predMu, predCov) = gp1.Predict(evalPts, GaussianProcess::DiagonalCov);
268
269 Eigen::MatrixXd predMu2, predCov2;
270 std::tie(predMu2, predCov2) = gp2.Predict(evalPts, GaussianProcess::DiagonalCov);
271
272 // Condition both GPs with some data
273 Eigen::MatrixXd obsLoc = 0.5*Eigen::MatrixXd::Ones(1,1);
274 Eigen::MatrixXd obsData = 0.25*Eigen::MatrixXd::Ones(1,1);
275 const double obsVar = 1e-2;
276
277 gp1.Condition(obsLoc, obsData, obsVar);
278 gp2.Condition(obsLoc, obsData, obsVar);
279
280 // add another observation
281 obsLoc(0) = 1.5;
282 gp1.Condition(obsLoc, obsData, obsVar);
283 gp2.Condition(obsLoc, obsData, obsVar);
284
285 // Make a posterior prediction
286 std::tie(predMu, predCov) = gp1.Predict(evalPts, GaussianProcess::DiagonalCov);
287 std::tie(predMu2, predCov2) = gp2.Predict(evalPts, GaussianProcess::DiagonalCov);
288
289 const double meanTol = 1e-2;
290 const double covTol = 1e-2;
291 for(int j=0; j<predMu.size(); ++j){
292 EXPECT_NEAR(predMu2(0,j), predMu(0,j), meanTol);
293 EXPECT_NEAR(predCov2(0,j), predCov(0,j), covTol);
294 }
295
296 }
297
TEST(Approximation_GP,StateSpacePredict_ExtrapolateLeft)298 TEST(Approximation_GP, StateSpacePredict_ExtrapolateLeft)
299 {
300
301 const double sigma2 = 1.0;
302 const double length = 0.3;
303 const double nu = 3.0/2.0;
304
305 MaternKernel kernel(1, sigma2, length, nu);
306
307 boost::property_tree::ptree options;
308 options.put("SDE.dt", 1e-5);
309
310 ZeroMean mu(1,1);
311 StateSpaceGP gp1(mu, kernel, options);
312 GaussianProcess gp2(mu, kernel);
313
314 const int numEvals = 100;
315 Eigen::MatrixXd evalPts(1,numEvals);
316 evalPts.row(0) = Eigen::VectorXd::LinSpaced(numEvals, -1, 0.2);
317
318 // Make a prediction about the prior
319 Eigen::MatrixXd predMu, predCov;
320 std::tie(predMu, predCov) = gp1.Predict(evalPts, GaussianProcess::DiagonalCov);
321
322 Eigen::MatrixXd predMu2, predCov2;
323 std::tie(predMu2, predCov2) = gp2.Predict(evalPts, GaussianProcess::DiagonalCov);
324
325 // Condition both GPs with some data
326 Eigen::MatrixXd obsLoc = 0.5*Eigen::MatrixXd::Ones(1,1);
327 Eigen::MatrixXd obsData = 0.25*Eigen::MatrixXd::Ones(1,1);
328 const double obsVar = 1e-2;
329
330 gp1.Condition(obsLoc, obsData, obsVar);
331 gp2.Condition(obsLoc, obsData, obsVar);
332
333 // add another observation
334 obsLoc(0) = 1.5;
335 gp1.Condition(obsLoc, obsData, obsVar);
336 gp2.Condition(obsLoc, obsData, obsVar);
337
338 // Make a posterior prediction
339 std::tie(predMu, predCov) = gp1.Predict(evalPts, GaussianProcess::DiagonalCov);
340 std::tie(predMu2, predCov2) = gp2.Predict(evalPts, GaussianProcess::DiagonalCov);
341
342 const double meanTol = 1e-2;
343 const double covTol = 1e-2;
344 for(int j=0; j<predMu.size(); ++j){
345 EXPECT_NEAR(predMu2(0,j), predMu(0,j), meanTol);
346 EXPECT_NEAR(predCov2(0,j), predCov(0,j), covTol);
347 }
348 }
349