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