base_leapfrog.hpp
Go to the documentation of this file.
1 #ifndef STAN_MCMC_HMC_INTEGRATORS_BASE_LEAPFROG_HPP
2 #define STAN_MCMC_HMC_INTEGRATORS_BASE_LEAPFROG_HPP
3 
6 #include <iostream>
7 #include <iomanip>
8 
9 namespace stan {
10  namespace mcmc {
11 
12  template <class Hamiltonian>
13  class base_leapfrog : public base_integrator<Hamiltonian> {
14  public:
16  : base_integrator<Hamiltonian>() {}
17 
18  void evolve(typename Hamiltonian::PointType& z,
19  Hamiltonian& hamiltonian,
20  const double epsilon,
21  callbacks::logger& logger) {
22  begin_update_p(z, hamiltonian, 0.5 * epsilon,
23  logger);
24  update_q(z, hamiltonian, epsilon,
25  logger);
26  end_update_p(z, hamiltonian, 0.5 * epsilon,
27  logger);
28  }
29 
30  void
31  verbose_evolve(typename Hamiltonian::PointType& z,
32  Hamiltonian& hamiltonian,
33  const double epsilon,
34  callbacks::logger& logger) {
35  std::stringstream msg;
36  msg.precision(6);
37 
38  int width = 14;
39  int nColumn = 4;
40 
41  msg << "Verbose Hamiltonian Evolution, Step Size = " << epsilon << ":";
42  logger.info(msg);
43 
44  msg.str("");
45  msg << " " << std::setw(nColumn * width)
46  << std::setfill('-')
47  << "" << std::setfill(' ');
48  logger.info(msg);
49 
50  msg.str("");
51  msg << " "
52  << std::setw(width) << std::left << "Poisson"
53  << std::setw(width) << std::left << "Initial"
54  << std::setw(width) << std::left << "Current"
55  << std::setw(width) << std::left << "DeltaH";
56  logger.info(msg);
57 
58  msg.str("");
59  msg << " "
60  << std::setw(width) << std::left << "Operator"
61  << std::setw(width) << std::left << "Hamiltonian"
62  << std::setw(width) << std::left << "Hamiltonian"
63  << std::setw(width) << std::left << "/ Stepsize^{2}";
64  logger.info(msg);
65 
66  msg.str("");
67  msg << " " << std::setw(nColumn * width)
68  << std::setfill('-')
69  << "" << std::setfill(' ');
70  logger.info(msg);
71 
72  double H0 = hamiltonian.H(z);
73 
74  begin_update_p(z, hamiltonian, 0.5 * epsilon,
75  logger);
76 
77  double H1 = hamiltonian.H(z);
78 
79  msg.str("");
80  msg << " "
81  << std::setw(width) << std::left << "hat{V}/2"
82  << std::setw(width) << std::left << H0
83  << std::setw(width) << std::left << H1
84  << std::setw(width) << std::left << (H1 - H0) / (epsilon * epsilon);
85  logger.info(msg);
86 
87  update_q(z, hamiltonian, epsilon, logger);
88 
89  double H2 = hamiltonian.H(z);
90 
91  msg.str("");
92  msg << " "
93  << std::setw(width) << std::left << "hat{T}"
94  << std::setw(width) << std::left << H0
95  << std::setw(width) << std::left << H2
96  << std::setw(width) << std::left << (H2 - H0) / (epsilon * epsilon);
97  logger.info(msg);
98 
99  end_update_p(z, hamiltonian, 0.5 * epsilon, logger);
100 
101  double H3 = hamiltonian.H(z);
102 
103  msg.str("");
104  msg << " "
105  << std::setw(width) << std::left << "hat{V}/2"
106  << std::setw(width) << std::left << H0
107  << std::setw(width) << std::left << H3
108  << std::setw(width) << std::left << (H3 - H0) / (epsilon * epsilon);
109  logger.info(msg);
110 
111  msg.str("");
112  msg << " "
113  << std::setw(nColumn * width)
114  << std::setfill('-')
115  << ""
116  << std::setfill(' ');
117  logger.info(msg);
118  }
119 
120  virtual
121  void begin_update_p(typename Hamiltonian::PointType& z,
122  Hamiltonian& hamiltonian, double epsilon,
123  callbacks::logger& logger) = 0;
124 
125  virtual
126  void update_q(typename Hamiltonian::PointType& z,
127  Hamiltonian& hamiltonian, double epsilon,
128  callbacks::logger& logger) = 0;
129  virtual
130  void end_update_p(typename Hamiltonian::PointType& z,
131  Hamiltonian& hamiltonian, double epsilon,
132  callbacks::logger& logger) = 0;
133  };
134 
135  } // mcmc
136 } // stan
137 #endif
#define H1(NAME, LABEL, N, X1, X2)
void verbose_evolve(typename Hamiltonian::PointType &z, Hamiltonian &hamiltonian, const double epsilon, callbacks::logger &logger)
virtual void update_q(typename Hamiltonian::PointType &z, Hamiltonian &hamiltonian, double epsilon, callbacks::logger &logger)=0
#define H2(NAME, LABEL, NX, X1, X2, NY, Y1, Y2)
virtual void begin_update_p(typename Hamiltonian::PointType &z, Hamiltonian &hamiltonian, double epsilon, callbacks::logger &logger)=0
virtual void end_update_p(typename Hamiltonian::PointType &z, Hamiltonian &hamiltonian, double epsilon, callbacks::logger &logger)=0
z
Definition: test.py:28
virtual void info(const std::string &message)
Definition: logger.hpp:47
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
void evolve(typename Hamiltonian::PointType &z, Hamiltonian &hamiltonian, const double epsilon, callbacks::logger &logger)
double epsilon