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