queso-0.53.0
InfiniteDimensionalMCMCSampler.C
Go to the documentation of this file.
1 //--------------------------------------------------------------------------
2 //
3 // QUESO - a library to support the Quantification of Uncertainty
4 // for Estimation, Simulation and Optimization
5 //
6 // Copyright (C) 2008,2009,2010,2011,2012,2013 The PECOS Development Team
7 //
8 // This library is free software; you can redistribute it and/or
9 // modify it under the terms of the Version 2.1 GNU Lesser General
10 // Public License as published by the Free Software Foundation.
11 //
12 // This library is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 // Lesser General Public License for more details.
16 //
17 // You should have received a copy of the GNU Lesser General Public
18 // License along with this library; if not, write to the Free Software
19 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
20 // Boston, MA 02110-1301 USA
21 //
22 //-----------------------------------------------------------------------el-
23 
24 #include <queso/Defines.h>
25 
26 #ifdef QUESO_HAVE_HDF5
27 
28 #include <cmath>
29 #include <algorithm>
30 #include <boost/math/constants/constants.hpp>
31 
32 // QUESO includes
33 #include <queso/Miscellaneous.h>
34 #include <queso/InfiniteDimensionalMeasureBase.h>
35 
36 #include <queso/InfiniteDimensionalLikelihoodBase.h>
37 #include <queso/InfiniteDimensionalMCMCSampler.h>
38 #include <queso/InfiniteDimensionalMCMCSamplerOptions.h>
39 
40 namespace QUESO {
41 
42 InfiniteDimensionalMCMCSampler::InfiniteDimensionalMCMCSampler(
43  const BaseEnvironment& env,
44  InfiniteDimensionalMeasureBase & prior,
45  InfiniteDimensionalLikelihoodBase & llhd,
46  InfiniteDimensionalMCMCSamplerOptions * ov)
47  : prior(prior),
48  llhd(llhd),
49  m_env(env),
50  current_physical_state(prior.draw()),
51  proposed_physical_state(prior.draw()),
52  current_physical_mean(prior.draw()),
53  current_physical_var(prior.draw()),
54  _delta(prior.draw()),
55  _M2(prior.draw())
56 {
57  if (ov != NULL) {
58  this->m_ov = ov;
59  }
60 
61 #ifdef QUESO_MEMORY_DEBUGGING
62  std::cout << "Entering InfiniteDimensionalMCMCSampler class" << std::endl;
63 #endif
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()
68  << std::endl;
69  }
70 
71 #ifdef QUESO_MEMORY_DEBUGGING
72  std::cout << "In InfiniteDimensionalMCMCSamplerOptions,"
73  << " finished scanning options" << std::endl;
74 #endif
75 
76  // Verify parent directory exists (for cases when a user specifies a
77  // relative path for the desired output file).
78  // Only subprocess of subrank 0 creates the output file
79  if ((this->m_env).subRank() == 0) {
80  int irtrn = CheckFilePath((this->m_ov->m_dataOutputDirName +
81  (this->m_env).subIdString() +
82  "/test.txt").c_str());
83  queso_require_greater_equal_msg(irtrn, 0, "unable to verify output path");
84  }
85 
86  // Only subprocess of subrank 0 creates the output file
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);
92  }
93 
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");
100 
101  // Ensure that we created the path and the output files
102  ((this->m_env).fullComm()).Barrier();
103 
104  this->_iteration = 0;
105  this->_acc_prob = 0.0;
106  this->_avg_acc_prob = 0.0;
107 
108  // FIXME: Should use the QUESO rng here instead
109  r = gsl_rng_alloc(gsl_rng_taus2);
110 
111  // Zero out these guys. There's probably a better way of doing this than
112  // calling prior.draw() at the start, but creation of a Sampler object
113  // should only be done a O(1) times anyway.
114  this->_delta->zero();
115  this->_M2->zero();
116  this->current_physical_mean->zero();
117  this->current_physical_var->zero();
118 
119  // LibMeshFunction & p = libmesh_cast_ref<LibMeshFunction &>(*(this->current_physical_state));
120 
121  // Seed the sampler at the truth
122  // for (unsigned int ii = 0; ii < 513; ii++) {
123  // p.equation_systems->get_system<ExplicitSystem>("Function").solution->set(ii, std::cos(2.0 * boost::math::constants::pi<double>() * ii / 512.0));
124  // }
125  // p.equation_systems->get_system<ExplicitSystem>("Function").solution->close();
126  // std::string seedname("seedfn");
127  // p.save_function(seedname);
128 
129  // Initialise cost function
130  this->_llhd_val = this->llhd.evaluate(*(this->current_physical_state));
131 
132  std::cout.setf(std::ios::fixed, std::ios::floatfield);
133  std::cout.precision(15);
134 }
135 
136 InfiniteDimensionalMCMCSampler::~InfiniteDimensionalMCMCSampler()
137 {
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);
146  }
147 }
148 
149 void InfiniteDimensionalMCMCSampler::_propose()
150 {
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);
153 
154  boost::shared_ptr<FunctionBase> p(prior.draw());
155 
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);
159 }
160 
161 void InfiniteDimensionalMCMCSampler::_metropolis_hastings()
162 {
163  // Downcast since we know we're dealing with a libmesh function
164  // LibMeshFunction & p = static_cast<LibMeshFunction &>(*(this->proposed_physical_state));
165  // LibMeshFunction & q = static_cast<LibMeshFunction &>(*(this->current_physical_state));
166  // Evaluate the likelihood at the proposed state
167  double proposed_llhd = this->llhd.evaluate(*(this->proposed_physical_state));
168  // double current_llhd = this->llhd.evaluate(q);
169  // std::cout << "llhd of current is: " << current_llhd << std::endl;
170  // std::cout << "llhd of proposal is: " << proposed_llhd << std::endl;
171 
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);
175  if (rand < alpha) {
176  // Accept
177  // std::cout << "accepted" << std::endl;
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;
182  }
183  else {
184  // std::cout << "rejected" << std::endl;
185  this->_acc_prob = 0.0;
186  }
187 }
188 
189 void InfiniteDimensionalMCMCSampler::_update_moments()
190 {
191  // Increment the current iteration number and update the running mean and
192  // variance.
193  this->_iteration += 1;
194 
195  // Make _delta contain the difference from the mean
196  this->_delta->zero();
197  this->_delta->add(1.0, *(this->current_physical_state));
198  this->_delta->add(-1.0, *(this->current_physical_mean));
199 
200  // Scale _delta to update the mean field
201  this->current_physical_mean->add(1.0 / this->iteration(), *(this->_delta));
202 
203  // Update running sum-of-squares
204  boost::shared_ptr<FunctionBase> temp_ptr(this->_delta->zero_clone());
205  // LibMeshFunction & temp = static_cast<LibMeshFunction &>(*temp_ptr);
206 
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);
211 
212  if (this->iteration() > 1) {
213  this->current_physical_var->zero();
214  this->current_physical_var->add(1.0 / (this->iteration() - 1), *(this->_M2));
215  }
216 
217  // Update acceptance rate
218  double delta_acc = this->acc_prob() - this->avg_acc_prob();
219  this->_avg_acc_prob += delta_acc / this->iteration();
220 }
221 
222 void InfiniteDimensionalMCMCSampler::step()
223 {
224  this->_propose();
225  this->_metropolis_hastings();
226  this->_update_moments();
227 
228  // We never save the 0th iteration
229  if (this->_iteration % this->m_ov->m_save_freq == 0) {
230  this->_write_state();
231  }
232 }
233 
234 hid_t InfiniteDimensionalMCMCSampler::_create_scalar_dataset(const std::string & name)
235 {
236  // Only subprocess with rank 0 manipulates the output file
237  if ((this->m_env).subRank() == 0) {
238  // Create a 1D dataspace. Unlimited size. Initially set to 0. We will
239  // extend it later
240  const int ndims = 1;
241  hsize_t dims[ndims] = {0}; // dataset dimensions at creation
242  hsize_t maxdims[ndims] = {H5S_UNLIMITED};
243 
244  hid_t file_space = H5Screate_simple(ndims, dims, maxdims);
245 
246  // Create dataset creation property list. Unlimited datasets must be
247  // chunked. Choosing the chunk size is an issue, here we set it to 1
248  hsize_t chunk_dims[ndims] = {1};
249  hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
250 
251  H5Pset_layout(plist, H5D_CHUNKED);
252  H5Pset_chunk(plist, ndims, chunk_dims);
253 
254  // Create the dataset
255  hid_t dset = H5Dcreate(this->_outfile, name.c_str(), H5T_NATIVE_DOUBLE,
256  file_space, H5P_DEFAULT, plist, H5P_DEFAULT);
257 
258  // We don't need the property list anymore. We also don't need the file
259  // dataspace anymore because we'll extend it later, making this one
260  // invalild anyway.
261  H5Pclose(plist);
262  H5Sclose(file_space);
263 
264  return dset;
265  }
266 
267  hid_t dummy;
268  return dummy;
269 }
270 
271 void InfiniteDimensionalMCMCSampler::_append_scalar_dataset(hid_t dset, double data)
272 {
273  // Only subprocess with rank 0 manipulates the output file
274  if ((this->m_env).subRank() == 0) {
275  int err;
276  // Create a memory dataspace for data to append
277  const int ndims = 1;
278  hsize_t dims[ndims] = { 1 }; // Only writing one double
279  hid_t mem_space = H5Screate_simple(ndims, dims, NULL);
280 
281  // Extend the dataset
282  // Set dims to be the *new* dimension of the extended dataset
283  dims[0] = this->_iteration / this->m_ov->m_save_freq;
284  err = H5Dset_extent(dset, dims);
285 
286  // Select hyperslab on file dataset
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};
290 
291  err = H5Sselect_hyperslab(file_space, H5S_SELECT_SET, start, NULL, count, NULL);
292 
293  // hsize_t size[1];
294  // size[0] = this->_iteration / this->m_ov->m_save_freq;
295 
296  // Write the data
297  H5Dwrite(dset, H5T_NATIVE_DOUBLE, mem_space, file_space, H5P_DEFAULT, &data);
298 
299  // Close a bunch of stuff
300  H5Sclose(file_space);
301  H5Sclose(mem_space);
302  }
303 }
304 
305 void InfiniteDimensionalMCMCSampler::_write_state()
306 {
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());
313 
314  // Now to write the fields. It appears to be a pain in the arse to write a
315  // method to spit this out to HDF5 format. Also, this won't scale to
316  // non-uniform finite element meshes. Therefore, I'm going to spit out an
317  // ExodusII file for each sample, average and variance. Got disk space?
318  std::stringstream basename;
319  basename << this->m_ov->m_dataOutputDirName;
320  basename << (this->m_env).subIdString();
321  basename << "/"; // Sigh
322 
323  std::ostringstream curr_iter;
324  curr_iter << this->iteration();
325 
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());
329 
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());
333 
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());
337 }
338 
339 boost::shared_ptr<InfiniteDimensionalMCMCSampler> InfiniteDimensionalMCMCSampler::clone_and_reset() const
340 {
341  // Set up a clone
342  boost::shared_ptr<InfiniteDimensionalMCMCSampler> clone(new InfiniteDimensionalMCMCSampler(this->m_env, this->prior, this->llhd, this->m_ov));
343 
344  // Copy the state.
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;
350 
351  return clone;
352 }
353 
354 double InfiniteDimensionalMCMCSampler::acc_prob()
355 {
356  return this->_acc_prob;
357 }
358 
359 double InfiniteDimensionalMCMCSampler::avg_acc_prob()
360 {
361  return this->_avg_acc_prob;
362 }
363 
364 double InfiniteDimensionalMCMCSampler::llhd_val() const
365 {
366  return this->_llhd_val;
367 }
368 
369 unsigned int InfiniteDimensionalMCMCSampler::iteration() const
370 {
371  return this->_iteration;
372 }
373 
374 } // End namespace QUESO
375 
376 #endif // QUESO_HAVE_HDF5
int CheckFilePath(const char *path)
#define queso_require_greater_equal_msg(expr1, expr2, msg)
Definition: asserts.h:90

Generated on Thu Jun 11 2015 13:52:31 for queso-0.53.0 by  doxygen 1.8.5