24 #include <queso/Defines.h>
26 #ifdef QUESO_HAVE_HDF5
30 #include <boost/math/constants/constants.hpp>
33 #include <queso/Miscellaneous.h>
34 #include <queso/InfiniteDimensionalMeasureBase.h>
36 #include <queso/InfiniteDimensionalLikelihoodBase.h>
37 #include <queso/InfiniteDimensionalMCMCSampler.h>
38 #include <queso/InfiniteDimensionalMCMCSamplerOptions.h>
42 InfiniteDimensionalMCMCSampler::InfiniteDimensionalMCMCSampler(
43 const BaseEnvironment& env,
44 InfiniteDimensionalMeasureBase & prior,
45 InfiniteDimensionalLikelihoodBase & llhd,
46 InfiniteDimensionalMCMCSamplerOptions * ov)
50 current_physical_state(prior.draw()),
51 proposed_physical_state(prior.draw()),
52 current_physical_mean(prior.draw()),
53 current_physical_var(prior.draw()),
61 #ifdef QUESO_MEMORY_DEBUGGING
62 std::cout <<
"Entering InfiniteDimensionalMCMCSampler class" << std::endl;
64 if (m_env.subDisplayFile()) {
65 *m_env.subDisplayFile() <<
"Entering InfiniteDimensionalMCMCSampler::constructor()"
66 <<
": prefix = " << this->m_ov->m_prefix
67 <<
", m_env.optionsInputFileName() = " << this->m_env.optionsInputFileName()
71 #ifdef QUESO_MEMORY_DEBUGGING
72 std::cout <<
"In InfiniteDimensionalMCMCSamplerOptions,"
73 <<
" finished scanning options" << std::endl;
79 if ((this->m_env).subRank() == 0) {
81 (this->m_env).subIdString() +
82 "/test.txt").c_str());
87 if ((this->m_env).subRank() == 0) {
88 this->_outfile = H5Fcreate((this->m_ov->m_dataOutputDirName +
89 (this->m_env).subIdString() +
"/" +
90 this->m_ov->m_dataOutputFileName).c_str(),
91 H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
94 this->_acc_dset = this->_create_scalar_dataset(
"acc");
95 this->_avg_acc_dset = this->_create_scalar_dataset(
"avg_acc");
96 this->_neg_log_llhd_dset = this->_create_scalar_dataset(
"neg_log_llhd");
97 this->_L2_norm_samples_dset = this->_create_scalar_dataset(
"L2_norm_samples");
98 this->_L2_norm_mean_dset = this->_create_scalar_dataset(
"L2_norm_mean");
99 this->_L2_norm_var_dset = this->_create_scalar_dataset(
"L2_norm_var");
102 ((this->m_env).fullComm()).Barrier();
104 this->_iteration = 0;
105 this->_acc_prob = 0.0;
106 this->_avg_acc_prob = 0.0;
109 r = gsl_rng_alloc(gsl_rng_taus2);
114 this->_delta->zero();
116 this->current_physical_mean->zero();
117 this->current_physical_var->zero();
130 this->_llhd_val = this->llhd.evaluate(*(this->current_physical_state));
132 std::cout.setf(std::ios::fixed, std::ios::floatfield);
133 std::cout.precision(15);
136 InfiniteDimensionalMCMCSampler::~InfiniteDimensionalMCMCSampler()
138 if ((this->m_env).subRank() == 0) {
139 H5Dclose(this->_acc_dset);
140 H5Dclose(this->_avg_acc_dset);
141 H5Dclose(this->_neg_log_llhd_dset);
142 H5Dclose(this->_L2_norm_samples_dset);
143 H5Dclose(this->_L2_norm_mean_dset);
144 H5Dclose(this->_L2_norm_var_dset);
145 H5Fclose(this->_outfile);
149 void InfiniteDimensionalMCMCSampler::_propose()
151 const double rwmh_step_sq = (this->m_ov->m_rwmh_step * this->m_ov->m_rwmh_step);
152 const double coeff = std::sqrt(1.0 - rwmh_step_sq);
154 boost::shared_ptr<FunctionBase> p(prior.draw());
156 this->proposed_physical_state->zero();
157 this->proposed_physical_state->add(coeff, *(this->current_physical_state));
158 this->proposed_physical_state->add(this->m_ov->m_rwmh_step, *p);
161 void InfiniteDimensionalMCMCSampler::_metropolis_hastings()
167 double proposed_llhd = this->llhd.evaluate(*(this->proposed_physical_state));
172 double diff = this->_llhd_val - proposed_llhd;
173 double alpha = std::min(1.0, std::exp(diff));
174 double rand = gsl_rng_uniform(this->r);
178 this->current_physical_state.swap(this->proposed_physical_state);
179 this->proposed_physical_state->zero();
180 this->_acc_prob = alpha;
181 this->_llhd_val = proposed_llhd;
185 this->_acc_prob = 0.0;
189 void InfiniteDimensionalMCMCSampler::_update_moments()
193 this->_iteration += 1;
196 this->_delta->zero();
197 this->_delta->add(1.0, *(this->current_physical_state));
198 this->_delta->add(-1.0, *(this->current_physical_mean));
201 this->current_physical_mean->add(1.0 / this->iteration(), *(this->_delta));
204 boost::shared_ptr<FunctionBase> temp_ptr(this->_delta->zero_clone());
207 temp_ptr->pointwise_mult(*(this->_delta), *(this->current_physical_state));
208 this->_M2->add(1.0, *temp_ptr);
209 temp_ptr->pointwise_mult(*(this->_delta), *(this->current_physical_mean));
210 this->_M2->add(-1.0, *temp_ptr);
212 if (this->iteration() > 1) {
213 this->current_physical_var->zero();
214 this->current_physical_var->add(1.0 / (this->iteration() - 1), *(this->_M2));
218 double delta_acc = this->acc_prob() - this->avg_acc_prob();
219 this->_avg_acc_prob += delta_acc / this->iteration();
222 void InfiniteDimensionalMCMCSampler::step()
225 this->_metropolis_hastings();
226 this->_update_moments();
229 if (this->_iteration % this->m_ov->m_save_freq == 0) {
230 this->_write_state();
234 hid_t InfiniteDimensionalMCMCSampler::_create_scalar_dataset(
const std::string & name)
237 if ((this->m_env).subRank() == 0) {
241 hsize_t dims[ndims] = {0};
242 hsize_t maxdims[ndims] = {H5S_UNLIMITED};
244 hid_t file_space = H5Screate_simple(ndims, dims, maxdims);
248 hsize_t chunk_dims[ndims] = {1};
249 hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
251 H5Pset_layout(plist, H5D_CHUNKED);
252 H5Pset_chunk(plist, ndims, chunk_dims);
255 hid_t dset = H5Dcreate(this->_outfile, name.c_str(), H5T_NATIVE_DOUBLE,
256 file_space, H5P_DEFAULT, plist, H5P_DEFAULT);
262 H5Sclose(file_space);
271 void InfiniteDimensionalMCMCSampler::_append_scalar_dataset(hid_t dset,
double data)
274 if ((this->m_env).subRank() == 0) {
278 hsize_t dims[ndims] = { 1 };
279 hid_t mem_space = H5Screate_simple(ndims, dims, NULL);
283 dims[0] = this->_iteration / this->m_ov->m_save_freq;
284 err = H5Dset_extent(dset, dims);
287 hid_t file_space = H5Dget_space(dset);
288 hsize_t start[1] = {(this->_iteration / this->m_ov->m_save_freq) - 1};
289 hsize_t count[1] = {1};
291 err = H5Sselect_hyperslab(file_space, H5S_SELECT_SET, start, NULL, count, NULL);
297 H5Dwrite(dset, H5T_NATIVE_DOUBLE, mem_space, file_space, H5P_DEFAULT, &data);
300 H5Sclose(file_space);
305 void InfiniteDimensionalMCMCSampler::_write_state()
307 this->_append_scalar_dataset(this->_acc_dset, this->acc_prob());
308 this->_append_scalar_dataset(this->_avg_acc_dset, this->avg_acc_prob());
309 this->_append_scalar_dataset(this->_neg_log_llhd_dset, this->_llhd_val);
310 this->_append_scalar_dataset(this->_L2_norm_samples_dset, this->current_physical_state->L2_norm());
311 this->_append_scalar_dataset(this->_L2_norm_mean_dset, this->current_physical_mean->L2_norm());
312 this->_append_scalar_dataset(this->_L2_norm_var_dset, this->current_physical_var->L2_norm());
318 std::stringstream basename;
319 basename << this->m_ov->m_dataOutputDirName;
320 basename << (this->m_env).subIdString();
323 std::ostringstream curr_iter;
324 curr_iter << this->iteration();
326 std::string sample_name(basename.str() +
"sample.e-s.");
327 sample_name += curr_iter.str();
328 this->current_physical_state->save_function(sample_name, this->iteration());
330 std::string mean_name(basename.str() +
"mean.e-s.");
331 mean_name += curr_iter.str();
332 this->current_physical_mean->save_function(mean_name, this->iteration());
334 std::string var_name(basename.str() +
"var.e-s.");
335 var_name += curr_iter.str();
336 this->current_physical_var->save_function(var_name, this->iteration());
339 boost::shared_ptr<InfiniteDimensionalMCMCSampler> InfiniteDimensionalMCMCSampler::clone_and_reset()
const
342 boost::shared_ptr<InfiniteDimensionalMCMCSampler> clone(
new InfiniteDimensionalMCMCSampler(this->m_env, this->prior, this->llhd, this->m_ov));
345 clone->current_physical_state = this->current_physical_state;
346 clone->proposed_physical_state = this->proposed_physical_state;
347 clone->_llhd_val = this->_llhd_val;
348 clone->_acc_prob = this->_acc_prob;
349 clone->_avg_acc_prob = 0.0;
354 double InfiniteDimensionalMCMCSampler::acc_prob()
356 return this->_acc_prob;
359 double InfiniteDimensionalMCMCSampler::avg_acc_prob()
361 return this->_avg_acc_prob;
364 double InfiniteDimensionalMCMCSampler::llhd_val()
const
366 return this->_llhd_val;
369 unsigned int InfiniteDimensionalMCMCSampler::iteration()
const
371 return this->_iteration;
376 #endif // QUESO_HAVE_HDF5
int CheckFilePath(const char *path)
#define queso_require_greater_equal_msg(expr1, expr2, msg)