1 #ifndef STAN_MCMC_HMC_UNIFORM_BASE_STATIC_UNIFORM_HPP
2 #define STAN_MCMC_HMC_UNIFORM_BASE_STATIC_UNIFORM_HPP
3 
4 #include <stan/callbacks/logger.hpp>
5 #include <stan/mcmc/hmc/base_hmc.hpp>
6 #include <stan/mcmc/hmc/hamiltonians/ps_point.hpp>
7 #include <boost/random/uniform_int_distribution.hpp>
8 #include <cmath>
9 #include <limits>
10 #include <string>
11 #include <vector>
12 
13 namespace stan {
14 
15 namespace mcmc {
16 /**
17  * Hamiltonian Monte Carlo implementation that uniformly samples
18  * from trajectories with a static integration time
19  */
20 template <class Model, template <class, class> class Hamiltonian,
21           template <class> class Integrator, class BaseRNG>
22 class base_static_uniform
23     : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
24  public:
base_static_uniform(const Model & model,BaseRNG & rng)25   base_static_uniform(const Model& model, BaseRNG& rng)
26       : base_hmc<Model, Hamiltonian, Integrator, BaseRNG>(model, rng),
27         T_(1),
28         energy_(0) {
29     update_L_();
30   }
31 
~base_static_uniform()32   ~base_static_uniform() {}
33 
transition(sample & init_sample,callbacks::logger & logger)34   sample transition(sample& init_sample, callbacks::logger& logger) {
35     this->sample_stepsize();
36 
37     this->seed(init_sample.cont_params());
38 
39     this->hamiltonian_.sample_p(this->z_, this->rand_int_);
40     this->hamiltonian_.init(this->z_, logger);
41 
42     ps_point z_init(this->z_);
43     double H0 = this->hamiltonian_.H(this->z_);
44 
45     ps_point z_sample(this->z_);
46     double sum_prob = 1;
47     double sum_metro_prob = 1;
48 
49     boost::random::uniform_int_distribution<> uniform(0, L_ - 1);
50     int Lp = uniform(this->rand_int_);
51 
52     for (int l = 0; l < Lp; ++l) {
53       this->integrator_.evolve(this->z_, this->hamiltonian_, -this->epsilon_,
54                                logger);
55 
56       double h = this->hamiltonian_.H(this->z_);
57       if (std::isnan(h))
58         h = std::numeric_limits<double>::infinity();
59 
60       double prob = std::exp(H0 - h);
61       sum_prob += prob;
62       sum_metro_prob += prob > 1 ? 1 : prob;
63 
64       if (this->rand_uniform_() < prob / sum_prob)
65         z_sample = this->z_;
66     }
67 
68     this->z_.ps_point::operator=(z_init);
69 
70     for (int l = 0; l < L_ - 1 - Lp; ++l) {
71       this->integrator_.evolve(this->z_, this->hamiltonian_, this->epsilon_,
72                                logger);
73 
74       double h = this->hamiltonian_.H(this->z_);
75       if (std::isnan(h))
76         h = std::numeric_limits<double>::infinity();
77 
78       double prob = std::exp(H0 - h);
79       sum_prob += prob;
80       sum_metro_prob += prob > 1 ? 1 : prob;
81 
82       if (this->rand_uniform_() < prob / sum_prob)
83         z_sample = this->z_;
84     }
85 
86     double accept_prob = sum_metro_prob / static_cast<double>(L_);
87 
88     this->z_.ps_point::operator=(z_sample);
89     this->energy_ = this->hamiltonian_.H(this->z_);
90     return sample(this->z_.q, -this->hamiltonian_.V(this->z_), accept_prob);
91   }
92 
get_sampler_param_names(std::vector<std::string> & names)93   void get_sampler_param_names(std::vector<std::string>& names) {
94     names.push_back("stepsize__");
95     names.push_back("int_time__");
96     names.push_back("energy__");
97   }
98 
get_sampler_params(std::vector<double> & values)99   void get_sampler_params(std::vector<double>& values) {
100     values.push_back(this->epsilon_);
101     values.push_back(this->T_);
102     values.push_back(this->energy_);
103   }
104 
set_nominal_stepsize_and_T(const double e,const double t)105   void set_nominal_stepsize_and_T(const double e, const double t) {
106     if (e > 0 && t > 0) {
107       this->nom_epsilon_ = e;
108       T_ = t;
109       update_L_();
110     }
111   }
112 
set_nominal_stepsize_and_L(const double e,const int l)113   void set_nominal_stepsize_and_L(const double e, const int l) {
114     if (e > 0 && l > 0) {
115       this->nom_epsilon_ = e;
116       L_ = l;
117       T_ = this->nom_epsilon_ * L_;
118     }
119   }
120 
set_T(const double t)121   void set_T(const double t) {
122     if (t > 0) {
123       T_ = t;
124       update_L_();
125     }
126   }
127 
set_nominal_stepsize(const double e)128   void set_nominal_stepsize(const double e) {
129     if (e > 0) {
130       this->nom_epsilon_ = e;
131       update_L_();
132     }
133   }
134 
get_T()135   double get_T() { return this->T_; }
136 
get_L()137   int get_L() { return this->L_; }
138 
139  protected:
140   double T_;
141   int L_;
142   double energy_;
143 
update_L_()144   void update_L_() {
145     L_ = static_cast<int>(T_ / this->nom_epsilon_);
146     L_ = L_ < 1 ? 1 : L_;
147   }
148 };
149 }  // namespace mcmc
150 }  // namespace stan
151 #endif
152