1 #include <stan/mcmc/hmc/integrators/expl_leapfrog.hpp>
2 #include <gtest/gtest.h>
3
4 #include <sstream>
5 #include <stan/callbacks/stream_logger.hpp>
6 #include <test/test-models/good/mcmc/hmc/integrators/command.hpp>
7
8 #include <stan/io/dump.hpp>
9
10 #include <stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp>
11 #include <stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp>
12 #include <boost/random/additive_combine.hpp> // L'Ecuyer RNG
13 #include <test/unit/util.hpp>
14
15 // namespace
16 //************************************************************
17
18 typedef boost::ecuyer1988 rng_t;
19
20 class McmcHmcIntegratorsExplLeapfrogF : public testing::Test {
21 public:
McmcHmcIntegratorsExplLeapfrogF()22 McmcHmcIntegratorsExplLeapfrogF()
23 : logger(debug, info, warn, error, fatal),
24 unit_e_integrator(),
25 diag_e_integrator() {}
26
SetUp()27 void SetUp() {
28 static const std::string DATA("mu <- 0.0\ny <- 0\n");
29 std::stringstream data_stream(DATA);
30 // setup hamiltonian
31 stan::io::dump data_var_context(data_stream);
32
33 model = new command_model_namespace::command_model(data_var_context);
34 debug.str("");
35 info.str("");
36 warn.str("");
37 error.str("");
38 fatal.str("");
39 }
40
TearDown()41 void TearDown() { delete (model); }
42
43 std::stringstream debug, info, warn, error, fatal;
44 stan::callbacks::stream_logger logger;
45
46 // integrator under test
47 stan::mcmc::expl_leapfrog<
48 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t> >
49 unit_e_integrator;
50
51 stan::mcmc::expl_leapfrog<
52 stan::mcmc::diag_e_metric<command_model_namespace::command_model, rng_t> >
53 diag_e_integrator;
54
55 // model
56 command_model_namespace::command_model *model;
57
58 // output
59 std::stringstream output;
60 };
61
TEST_F(McmcHmcIntegratorsExplLeapfrogF,begin_update_p)62 TEST_F(McmcHmcIntegratorsExplLeapfrogF, begin_update_p) {
63 // setup z
64 stan::mcmc::unit_e_point z(1);
65 z.V = 1.99974742955684;
66 z.q(0) = 1.99987371079118;
67 z.p(0) = -1.58612292129732;
68 z.g(0) = 1.99987371079118;
69 EXPECT_NEAR(z.V, 1.99974742955684, 1e-15);
70 EXPECT_NEAR(z.q(0), 1.99987371079118, 1e-15);
71 EXPECT_NEAR(z.p(0), -1.58612292129732, 1e-15);
72 EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
73
74 // setup hamiltonian
75 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
76 hamiltonian(*model);
77
78 // setup epsilon
79 double epsilon = 0.1;
80
81 unit_e_integrator.begin_update_p(z, hamiltonian, 0.5 * epsilon, logger);
82 EXPECT_NEAR(z.V, 1.99974742955684, 5e-14);
83 EXPECT_NEAR(z.q(0), 1.99987371079118, 5e-14);
84 EXPECT_NEAR(z.p(0), -1.68611660683688, 5e-14);
85 EXPECT_NEAR(z.g(0), 1.99987371079118, 5e-14);
86
87 EXPECT_EQ("", debug.str());
88 EXPECT_EQ("", info.str());
89 EXPECT_EQ("", warn.str());
90 EXPECT_EQ("", error.str());
91 EXPECT_EQ("", fatal.str());
92 }
93
TEST_F(McmcHmcIntegratorsExplLeapfrogF,update_q)94 TEST_F(McmcHmcIntegratorsExplLeapfrogF, update_q) {
95 // setup z
96 stan::mcmc::unit_e_point z(1);
97 z.V = 1.99974742955684;
98 z.q(0) = 1.99987371079118;
99 z.p(0) = -1.68611660683688;
100 z.g(0) = 1.99987371079118;
101 EXPECT_NEAR(z.V, 1.99974742955684, 1e-15);
102 EXPECT_NEAR(z.q(0), 1.99987371079118, 1e-15);
103 EXPECT_NEAR(z.p(0), -1.68611660683688, 1e-15);
104 EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
105
106 // setup hamiltonian
107 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
108 hamiltonian(*model);
109
110 // setup epsilon
111 double epsilon = 0.1;
112
113 unit_e_integrator.update_q(z, hamiltonian, epsilon, logger);
114 EXPECT_NEAR(z.V, 1.6767603480819471, 5e-14);
115 EXPECT_NEAR(z.q(0), 1.83126205010749, 5e-14);
116 EXPECT_NEAR(z.p(0), -1.68611660683688, 5e-14);
117 EXPECT_NEAR(z.g(0), 1.8312620501074919, 5e-14);
118
119 EXPECT_EQ("", debug.str());
120 EXPECT_EQ("", info.str());
121 EXPECT_EQ("", warn.str());
122 EXPECT_EQ("", error.str());
123 EXPECT_EQ("", fatal.str());
124 }
125
TEST_F(McmcHmcIntegratorsExplLeapfrogF,end_update_p)126 TEST_F(McmcHmcIntegratorsExplLeapfrogF, end_update_p) {
127 // setup z
128 stan::mcmc::unit_e_point z(1);
129 z.V = 1.39887860643153;
130 z.q(0) = 1.67264975797776;
131 z.p(0) = -1.68611660683688;
132 z.g(0) = 1.67264975797776;
133 EXPECT_NEAR(z.V, 1.39887860643153, 1e-15);
134 EXPECT_NEAR(z.q(0), 1.67264975797776, 1e-15);
135 EXPECT_NEAR(z.p(0), -1.68611660683688, 1e-15);
136 EXPECT_NEAR(z.g(0), 1.67264975797776, 1e-15);
137
138 // setup hamiltonian
139 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
140 hamiltonian(*model);
141
142 // setup epsilon
143 double epsilon = 0.1;
144
145 unit_e_integrator.end_update_p(z, hamiltonian, 0.5 * epsilon, logger);
146 EXPECT_NEAR(z.V, 1.39887860643153, 5e-14);
147 EXPECT_NEAR(z.q(0), 1.67264975797776, 5e-14);
148 EXPECT_NEAR(z.p(0), -1.76974909473577, 5e-14);
149 EXPECT_NEAR(z.g(0), 1.67264975797776, 5e-14);
150
151 EXPECT_EQ("", debug.str());
152 EXPECT_EQ("", info.str());
153 EXPECT_EQ("", warn.str());
154 EXPECT_EQ("", error.str());
155 EXPECT_EQ("", fatal.str());
156 }
157
TEST_F(McmcHmcIntegratorsExplLeapfrogF,evolve_1)158 TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_1) {
159 // setup z
160 stan::mcmc::unit_e_point z(1);
161 z.V = 1.99974742955684;
162 z.q(0) = 1.99987371079118;
163 z.p(0) = -1.58612292129732;
164 z.g(0) = 1.99987371079118;
165 EXPECT_NEAR(z.V, 1.99974742955684, 1e-15);
166 EXPECT_NEAR(z.q(0), 1.99987371079118, 1e-15);
167 EXPECT_NEAR(z.p(0), -1.58612292129732, 1e-15);
168 EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
169
170 // setup hamiltonian
171 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
172 hamiltonian(*model);
173
174 // setup epsilon
175 double epsilon = 0.1;
176
177 unit_e_integrator.evolve(z, hamiltonian, epsilon, logger);
178 EXPECT_NEAR(z.V, 1.67676034808195, 5e-14);
179 EXPECT_NEAR(z.q(0), 1.83126205010749, 5e-14);
180 EXPECT_NEAR(z.p(0), -1.77767970934226, 5e-14);
181 EXPECT_NEAR(z.g(0), 1.83126205010749, 5e-14);
182
183 EXPECT_EQ("", debug.str());
184 EXPECT_EQ("", info.str());
185 EXPECT_EQ("", warn.str());
186 EXPECT_EQ("", error.str());
187 EXPECT_EQ("", fatal.str());
188 }
189
TEST_F(McmcHmcIntegratorsExplLeapfrogF,evolve_2)190 TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_2) {
191 // setup z
192 stan::mcmc::unit_e_point z(1);
193 z.V = 1.99974742955684;
194 z.q(0) = 1.99987371079118;
195 z.p(0) = 0.531143888645192;
196 z.g(0) = 1.99987371079118;
197 EXPECT_NEAR(z.V, 1.99974742955684, 1e-15);
198 EXPECT_NEAR(z.q(0), 1.99987371079118, 1e-15);
199 EXPECT_NEAR(z.p(0), 0.531143888645192, 1e-15);
200 EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
201
202 // setup hamiltonian
203 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
204 hamiltonian(*model);
205
206 // setup epsilon
207 double epsilon = 0.2;
208
209 unit_e_integrator.evolve(z, hamiltonian, epsilon, logger);
210 EXPECT_NEAR(z.V, 2.13439496506688, 5e-14);
211 EXPECT_NEAR(z.q(0), 2.06610501430439, 5e-14);
212 EXPECT_NEAR(z.p(0), 0.124546016135635, 5e-14);
213 EXPECT_NEAR(z.g(0), 2.06610501430439, 5e-14);
214
215 EXPECT_EQ("", debug.str());
216 EXPECT_EQ("", info.str());
217 EXPECT_EQ("", warn.str());
218 EXPECT_EQ("", error.str());
219 EXPECT_EQ("", fatal.str());
220 }
221
TEST_F(McmcHmcIntegratorsExplLeapfrogF,evolve_3)222 TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_3) {
223 // setup z
224 stan::mcmc::unit_e_point z(1);
225 z.V = 1.99974742955684;
226 z.q(0) = 1.99987371079118;
227 z.p(0) = 0.531143888645192;
228 z.g(0) = 1.99987371079118;
229 EXPECT_NEAR(z.V, 1.99974742955684, 1e-15);
230 EXPECT_NEAR(z.q(0), 1.99987371079118, 1e-15);
231 EXPECT_NEAR(z.p(0), 0.531143888645192, 1e-15);
232 EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
233
234 // setup hamiltonian
235 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
236 hamiltonian(*model);
237
238 // setup epsilon
239 double epsilon = 0.2;
240
241 unit_e_integrator.evolve(z, hamiltonian, epsilon, logger);
242 EXPECT_NEAR(z.V, 2.13439496506688, 5e-14);
243 EXPECT_NEAR(z.q(0), 2.06610501430439, 5e-14);
244 EXPECT_NEAR(z.p(0), 0.124546016135635, 5e-14);
245 EXPECT_NEAR(z.g(0), 2.06610501430439, 5e-14);
246
247 EXPECT_EQ("", debug.str());
248 EXPECT_EQ("", info.str());
249 EXPECT_EQ("", warn.str());
250 EXPECT_EQ("", error.str());
251 EXPECT_EQ("", fatal.str());
252 }
253
TEST_F(McmcHmcIntegratorsExplLeapfrogF,evolve_4)254 TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_4) {
255 // setup z
256 stan::mcmc::unit_e_point z(1);
257 z.V = 1.99974742955684;
258 z.q(0) = 1.99987371079118;
259 z.p(0) = -1.01150787313287;
260 z.g(0) = 1.99987371079118;
261 EXPECT_NEAR(z.V, 1.99974742955684, 1e-15);
262 EXPECT_NEAR(z.q(0), 1.99987371079118, 1e-15);
263 EXPECT_NEAR(z.p(0), -1.01150787313287, 1e-15);
264 EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
265
266 // setup hamiltonian
267 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
268 hamiltonian(*model);
269
270 // setup epsilon
271 double epsilon = 0.4;
272
273 unit_e_integrator.evolve(z, hamiltonian, epsilon, logger);
274 EXPECT_NEAR(z.V, 1.03001529319458, 5e-14);
275 EXPECT_NEAR(z.q(0), 1.43528066467474, 5e-14);
276 EXPECT_NEAR(z.p(0), -1.69853874822605, 5e-14);
277 EXPECT_NEAR(z.g(0), 1.43528066467474, 5e-14);
278
279 EXPECT_EQ("", debug.str());
280 EXPECT_EQ("", info.str());
281 EXPECT_EQ("", warn.str());
282 EXPECT_EQ("", error.str());
283 EXPECT_EQ("", fatal.str());
284 }
285
TEST_F(McmcHmcIntegratorsExplLeapfrogF,evolve_5)286 TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_5) {
287 // setup z
288 stan::mcmc::unit_e_point z(1);
289 z.V = 1.99974742955684;
290 z.q(0) = 1.99987371079118;
291 z.p(0) = -0.141638464197442;
292 z.g(0) = 1.99987371079118;
293 EXPECT_NEAR(z.V, 1.99974742955684, 1e-15);
294 EXPECT_NEAR(z.q(0), 1.99987371079118, 1e-15);
295 EXPECT_NEAR(z.p(0), -0.141638464197442, 1e-15);
296 EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
297
298 // setup hamiltonian
299 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
300 hamiltonian(*model);
301
302 // setup epsilon
303 double epsilon = 0.8;
304
305 unit_e_integrator.evolve(z, hamiltonian, epsilon, logger);
306 EXPECT_NEAR(z.V, 0.777009958583946, 5e-14);
307 EXPECT_NEAR(z.q(0), 1.24660335198005, 5e-14);
308 EXPECT_NEAR(z.p(0), -1.44022928930593, 5e-14);
309 EXPECT_NEAR(z.g(0), 1.24660335198005, 5e-14);
310
311 EXPECT_EQ("", debug.str());
312 EXPECT_EQ("", info.str());
313 EXPECT_EQ("", warn.str());
314 EXPECT_EQ("", error.str());
315 EXPECT_EQ("", fatal.str());
316 }
TEST_F(McmcHmcIntegratorsExplLeapfrogF,evolve_6)317 TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_6) {
318 // setup z
319 stan::mcmc::unit_e_point z(1);
320 z.V = 1.99974742955684;
321 z.q(0) = 1.99987371079118;
322 z.p(0) = 0.0942249427134016;
323 z.g(0) = 1.99987371079118;
324 EXPECT_NEAR(z.V, 1.99974742955684, 1e-15);
325 EXPECT_NEAR(z.q(0), 1.99987371079118, 1e-15);
326 EXPECT_NEAR(z.p(0), 0.0942249427134016, 1e-15);
327 EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
328
329 // setup hamiltonian
330 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
331 hamiltonian(*model);
332
333 // setup epsilon
334 double epsilon = 1.6;
335
336 unit_e_integrator.evolve(z, hamiltonian, epsilon, logger);
337 EXPECT_NEAR(z.V, 0.0837242558054816, 5e-14);
338 EXPECT_NEAR(z.q(0), -0.409204730680088, 5e-14);
339 EXPECT_NEAR(z.p(0), -1.17831024137547, 5e-14);
340 EXPECT_NEAR(z.g(0), -0.409204730680088, 5e-14);
341
342 EXPECT_EQ("", debug.str());
343 EXPECT_EQ("", info.str());
344 EXPECT_EQ("", warn.str());
345 EXPECT_EQ("", error.str());
346 EXPECT_EQ("", fatal.str());
347 }
348
TEST_F(McmcHmcIntegratorsExplLeapfrogF,evolve_7)349 TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_7) {
350 // setup z
351 stan::mcmc::unit_e_point z(1);
352 z.V = 1.99974742955684;
353 z.q(0) = 1.99987371079118;
354 z.p(0) = 1.01936184962275;
355 z.g(0) = 1.99987371079118;
356 EXPECT_NEAR(z.V, 1.99974742955684, 1e-15);
357 EXPECT_NEAR(z.q(0), 1.99987371079118, 1e-15);
358 EXPECT_NEAR(z.p(0), 1.01936184962275, 1e-15);
359 EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
360
361 // setup hamiltonian
362 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
363 hamiltonian(*model);
364
365 // setup epsilon
366 double epsilon = 3.2;
367
368 unit_e_integrator.evolve(z, hamiltonian, epsilon, logger);
369 EXPECT_NEAR(z.V, 12.3878614837537, 5e-13);
370 EXPECT_NEAR(z.q(0), -4.97752176966686, 5e-14);
371 EXPECT_NEAR(z.p(0), 5.78359874382383, 5e-14);
372 EXPECT_NEAR(z.g(0), -4.97752176966686, 5e-14);
373
374 EXPECT_EQ("", debug.str());
375 EXPECT_EQ("", info.str());
376 EXPECT_EQ("", warn.str());
377 EXPECT_EQ("", error.str());
378 EXPECT_EQ("", fatal.str());
379 }
380
TEST_F(McmcHmcIntegratorsExplLeapfrogF,evolve_8)381 TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_8) {
382 // setup z
383 stan::mcmc::unit_e_point z(1);
384 z.V = 1.99974742955684;
385 z.q(0) = 1.99987371079118;
386 z.p(0) = -2.73131279771964;
387 z.g(0) = 1.99987371079118;
388 EXPECT_NEAR(z.V, 1.99974742955684, 1e-15);
389 EXPECT_NEAR(z.q(0), 1.99987371079118, 1e-15);
390 EXPECT_NEAR(z.p(0), -2.73131279771964, 1e-15);
391 EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
392
393 // setup hamiltonian
394 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
395 hamiltonian(*model);
396
397 // setup epsilon
398 double epsilon = -1;
399
400 unit_e_integrator.evolve(z, hamiltonian, epsilon, logger);
401 EXPECT_NEAR(z.V, 6.96111198693627, 5e-14);
402 EXPECT_NEAR(z.q(0), 3.73124965311523, 5e-14);
403 EXPECT_NEAR(z.p(0), 0.134248884233563, 5e-14);
404 EXPECT_NEAR(z.g(0), 3.73124965311523, 5e-14);
405
406 EXPECT_EQ("", debug.str());
407 EXPECT_EQ("", info.str());
408 EXPECT_EQ("", warn.str());
409 EXPECT_EQ("", error.str());
410 EXPECT_EQ("", fatal.str());
411 }
412
TEST_F(McmcHmcIntegratorsExplLeapfrogF,evolve_9)413 TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_9) {
414 // setup z
415 stan::mcmc::diag_e_point z(1);
416 z.V = 0.807684865121721;
417 z.q(0) = 1.27097196280777;
418 z.p(0) = -0.159996782671291;
419 z.g(0) = 1.27097196280777;
420 z.inv_e_metric_(0) = 0.733184698671436;
421 EXPECT_NEAR(z.V, 0.807684865121721, 1e-15);
422 EXPECT_NEAR(z.q(0), 1.27097196280777, 1e-15);
423 EXPECT_NEAR(z.p(0), -0.159996782671291, 1e-15);
424 EXPECT_NEAR(z.g(0), 1.27097196280777, 1e-15);
425
426 // setup hamiltonian
427 stan::mcmc::diag_e_metric<command_model_namespace::command_model, rng_t>
428 hamiltonian(*model);
429
430 // setup epsilon
431 double epsilon = 2.40769920051673;
432
433 diag_e_integrator.evolve(z, hamiltonian, epsilon, logger);
434 EXPECT_NEAR(z.V, 1.46626604258356, 5e-14);
435 EXPECT_NEAR(z.q(0), -1.71246374711032, 5e-14);
436 EXPECT_NEAR(z.p(0), 0.371492925378682, 5e-14);
437 EXPECT_NEAR(z.g(0), -1.71246374711032, 5e-14);
438
439 EXPECT_EQ("", debug.str());
440 EXPECT_EQ("", info.str());
441 EXPECT_EQ("", warn.str());
442 EXPECT_EQ("", error.str());
443 EXPECT_EQ("", fatal.str());
444 }
445
TEST_F(McmcHmcIntegratorsExplLeapfrogF,streams)446 TEST_F(McmcHmcIntegratorsExplLeapfrogF, streams) {
447 stan::test::capture_std_streams();
448
449 typedef stan::mcmc::expl_leapfrog<
450 stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t> >
451 integrator;
452
453 EXPECT_NO_THROW(integrator i);
454
455 stan::test::reset_std_streams();
456 EXPECT_EQ("", stan::test::cout_ss.str());
457 EXPECT_EQ("", stan::test::cerr_ss.str());
458 }
459