queso-0.57.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
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 <sstream>
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 
43  const BaseEnvironment& env,
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  // Zero out these guys. There's probably a better way of doing this than
109  // calling prior.draw() at the start, but creation of a Sampler object
110  // should only be done a O(1) times anyway.
111  this->_delta->zero();
112  this->_M2->zero();
113  this->current_physical_mean->zero();
114  this->current_physical_var->zero();
115 
116  // LibMeshFunction & p = libmesh_cast_ref<LibMeshFunction &>(*(this->current_physical_state));
117 
118  // Seed the sampler at the truth
119  // for (unsigned int ii = 0; ii < 513; ii++) {
120  // p.equation_systems->get_system<ExplicitSystem>("Function").solution->set(ii, std::cos(2.0 * boost::math::constants::pi<double>() * ii / 512.0));
121  // }
122  // p.equation_systems->get_system<ExplicitSystem>("Function").solution->close();
123  // std::string seedname("seedfn");
124  // p.save_function(seedname);
125 
126  // Initialise cost function
127  this->_llhd_val = this->llhd.evaluate(*(this->current_physical_state));
128 
129  std::cout.setf(std::ios::fixed, std::ios::floatfield);
130  std::cout.precision(15);
131 }
132 
134 {
135  if ((this->m_env).subRank() == 0) {
136  H5Dclose(this->_acc_dset);
137  H5Dclose(this->_avg_acc_dset);
138  H5Dclose(this->_neg_log_llhd_dset);
139  H5Dclose(this->_L2_norm_samples_dset);
140  H5Dclose(this->_L2_norm_mean_dset);
141  H5Dclose(this->_L2_norm_var_dset);
142  H5Fclose(this->_outfile);
143  }
144 }
145 
147 {
148  const double rwmh_step_sq = (this->m_ov->m_rwmh_step * this->m_ov->m_rwmh_step);
149  const double coeff = std::sqrt(1.0 - rwmh_step_sq);
150 
152 
153  this->proposed_physical_state->zero();
154  this->proposed_physical_state->add(coeff, *(this->current_physical_state));
155  this->proposed_physical_state->add(this->m_ov->m_rwmh_step, *p);
156 }
157 
159 {
160  // Downcast since we know we're dealing with a libmesh function
161  // LibMeshFunction & p = static_cast<LibMeshFunction &>(*(this->proposed_physical_state));
162  // LibMeshFunction & q = static_cast<LibMeshFunction &>(*(this->current_physical_state));
163  // Evaluate the likelihood at the proposed state
164  double proposed_llhd = this->llhd.evaluate(*(this->proposed_physical_state));
165  // double current_llhd = this->llhd.evaluate(q);
166  // std::cout << "llhd of current is: " << current_llhd << std::endl;
167  // std::cout << "llhd of proposal is: " << proposed_llhd << std::endl;
168 
169  double diff = this->_llhd_val - proposed_llhd;
170  double alpha = std::min(1.0, std::exp(diff));
171  double rand = this->m_env.rngObject()->uniformSample();
172  if (rand < alpha) {
173  // Accept
174  // std::cout << "accepted" << std::endl;
176  this->proposed_physical_state->zero();
177  this->_acc_prob = alpha;
178  this->_llhd_val = proposed_llhd;
179  }
180  else {
181  // std::cout << "rejected" << std::endl;
182  this->_acc_prob = 0.0;
183  }
184 }
185 
187 {
188  // Increment the current iteration number and update the running mean and
189  // variance.
190  this->_iteration += 1;
191 
192  // Make _delta contain the difference from the mean
193  this->_delta->zero();
194  this->_delta->add(1.0, *(this->current_physical_state));
195  this->_delta->add(-1.0, *(this->current_physical_mean));
196 
197  // Scale _delta to update the mean field
198  this->current_physical_mean->add(1.0 / this->iteration(), *(this->_delta));
199 
200  // Update running sum-of-squares
201  SharedPtr<FunctionBase>::Type temp_ptr(this->_delta->zero_clone());
202  // LibMeshFunction & temp = static_cast<LibMeshFunction &>(*temp_ptr);
203 
204  temp_ptr->pointwise_mult(*(this->_delta), *(this->current_physical_state));
205  this->_M2->add(1.0, *temp_ptr);
206  temp_ptr->pointwise_mult(*(this->_delta), *(this->current_physical_mean));
207  this->_M2->add(-1.0, *temp_ptr);
208 
209  if (this->iteration() > 1) {
210  this->current_physical_var->zero();
211  this->current_physical_var->add(1.0 / (this->iteration() - 1), *(this->_M2));
212  }
213 
214  // Update acceptance rate
215  double delta_acc = this->acc_prob() - this->avg_acc_prob();
216  this->_avg_acc_prob += delta_acc / this->iteration();
217 }
218 
220 {
221  this->_propose();
222  this->_metropolis_hastings();
223  this->_update_moments();
224 
225  // We never save the 0th iteration
226  if (this->_iteration % this->m_ov->m_save_freq == 0) {
227  this->_write_state();
228  }
229 }
230 
232 {
233  // Only subprocess with rank 0 manipulates the output file
234  if ((this->m_env).subRank() == 0) {
235  // Create a 1D dataspace. Unlimited size. Initially set to 0. We will
236  // extend it later
237  const int ndims = 1;
238  hsize_t dims[ndims] = {0}; // dataset dimensions at creation
239  hsize_t maxdims[ndims] = {H5S_UNLIMITED};
240 
241  hid_t file_space = H5Screate_simple(ndims, dims, maxdims);
242 
243  // Create dataset creation property list. Unlimited datasets must be
244  // chunked. Choosing the chunk size is an issue, here we set it to 1
245  hsize_t chunk_dims[ndims] = {1};
246  hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
247 
248  H5Pset_layout(plist, H5D_CHUNKED);
249  H5Pset_chunk(plist, ndims, chunk_dims);
250 
251  // Create the dataset
252  hid_t dset = H5Dcreate(this->_outfile, name.c_str(), H5T_NATIVE_DOUBLE,
253  file_space, H5P_DEFAULT, plist, H5P_DEFAULT);
254 
255  // We don't need the property list anymore. We also don't need the file
256  // dataspace anymore because we'll extend it later, making this one
257  // invalild anyway.
258  H5Pclose(plist);
259  H5Sclose(file_space);
260 
261  return dset;
262  }
263 
264  hid_t dummy = -1;
265  return dummy;
266 }
267 
269 {
270  // Only subprocess with rank 0 manipulates the output file
271  if ((this->m_env).subRank() == 0) {
272  int err;
273  // Create a memory dataspace for data to append
274  const int ndims = 1;
275  hsize_t dims[ndims] = { 1 }; // Only writing one double
276  hid_t mem_space = H5Screate_simple(ndims, dims, NULL);
277 
278  // Extend the dataset
279  // Set dims to be the *new* dimension of the extended dataset
280  dims[0] = this->_iteration / this->m_ov->m_save_freq;
281  err = H5Dset_extent(dset, dims);
282  queso_require_greater_equal_msg(err, 0, "H5DSet_extent(dset, dims) failed");
283 
284  // Select hyperslab on file dataset
285  hid_t file_space = H5Dget_space(dset);
286  hsize_t start[1] = {(this->_iteration / this->m_ov->m_save_freq) - 1};
287  hsize_t count[1] = {1};
288 
289  err = H5Sselect_hyperslab(file_space, H5S_SELECT_SET, start, NULL, count, NULL);
290  queso_require_greater_equal_msg(err, 0, "H5Sselect_hyperslab(file_space, H5S_SELECT_SET, start, NULL, count, NULL) failed");
291 
292  // hsize_t size[1];
293  // size[0] = this->_iteration / this->m_ov->m_save_freq;
294 
295  // Write the data
296  H5Dwrite(dset, H5T_NATIVE_DOUBLE, mem_space, file_space, H5P_DEFAULT, &data);
297 
298  // Close a bunch of stuff
299  H5Sclose(file_space);
300  H5Sclose(mem_space);
301  }
302 }
303 
305 {
306  this->_append_scalar_dataset(this->_acc_dset, this->acc_prob());
307  this->_append_scalar_dataset(this->_avg_acc_dset, this->avg_acc_prob());
310  this->_append_scalar_dataset(this->_L2_norm_mean_dset, this->current_physical_mean->L2_norm());
311  this->_append_scalar_dataset(this->_L2_norm_var_dset, this->current_physical_var->L2_norm());
312 
313  // Now to write the fields. It appears to be a pain in the arse to write a
314  // method to spit this out to HDF5 format. Also, this won't scale to
315  // non-uniform finite element meshes. Therefore, I'm going to spit out an
316  // ExodusII file for each sample, average and variance. Got disk space?
317  std::stringstream basename;
318  basename << this->m_ov->m_dataOutputDirName;
319  basename << (this->m_env).subIdString();
320  basename << "/"; // Sigh
321 
322  std::ostringstream curr_iter;
323  curr_iter << this->iteration();
324 
325  std::string sample_name(basename.str() + "sample.e-s.");
326  sample_name += curr_iter.str();
327  this->current_physical_state->save_function(sample_name, this->iteration());
328 
329  std::string mean_name(basename.str() + "mean.e-s.");
330  mean_name += curr_iter.str();
331  this->current_physical_mean->save_function(mean_name, this->iteration());
332 
333  std::string var_name(basename.str() + "var.e-s.");
334  var_name += curr_iter.str();
335  this->current_physical_var->save_function(var_name, this->iteration());
336 }
337 
339 {
340  // Set up a clone
342 
343  // Copy the state.
344  clone->current_physical_state = this->current_physical_state;
345  clone->proposed_physical_state = this->proposed_physical_state;
346  clone->_llhd_val = this->_llhd_val;
347  clone->_acc_prob = this->_acc_prob;
348  clone->_avg_acc_prob = 0.0;
349 
350  return clone;
351 }
352 
354 {
355  return this->_acc_prob;
356 }
357 
359 {
360  return this->_avg_acc_prob;
361 }
362 
364 {
365  return this->_llhd_val;
366 }
367 
369 {
370  return this->_iteration;
371 }
372 
373 } // End namespace QUESO
374 
375 #endif // QUESO_HAVE_HDF5
std::shared_ptr< T > Type
Definition: SharedPtr.h:69
This (virtual) class sets up the environment underlying the use of the QUESO library by an executable...
Definition: Environment.h:198
SharedPtr< FunctionBase >::Type current_physical_mean
std::string m_dataOutputFileName
Name of the HDF5 output file to store chain statistics.
InfiniteDimensionalMCMCSampler(const BaseEnvironment &env, InfiniteDimensionalMeasureBase &prior, InfiniteDimensionalLikelihoodBase &llhd, InfiniteDimensionalMCMCSamplerOptions *ov)
Construct an infinite dimensional MCMC chain with given prior, likelihood llhd, and options ov...
Abstract class representing the likelihood. Users must subclass this.
InfiniteDimensionalMCMCSamplerOptions * m_ov
double llhd_val() const
Get the current value of the llhd.
std::string m_dataOutputDirName
Name of the output dir to save infinite dimensional output files to.
std::string optionsInputFileName() const
Access to the attribute m_optionsInputFileName, which stores the name of the input file passed by the...
Definition: Environment.C:354
const RngBase * rngObject() const
Access to the RNG object.
Definition: Environment.C:471
SharedPtr< InfiniteDimensionalMCMCSampler >::Type clone_and_reset() const
Returns a pointer to new sampler, with all the moments reset.
hid_t _create_scalar_dataset(const std::string &name)
unsigned int iteration() const
Returns the current iteration number.
int CheckFilePath(const char *path)
double acc_prob()
Returns current acceptance probability.
Abstract base class for infinite dimensional measures.
unsigned int m_save_freq
The frequency at which to save the state of the chain.
void step()
Do one iteration of the Markov chain.
This class defines the options that specify the behaviour of the MCMC sampler.
SharedPtr< FunctionBase >::Type current_physical_state
double avg_acc_prob()
Returns current average acceptance probability.
SharedPtr< FunctionBase >::Type current_physical_var
std::string m_prefix
The prefix to look for in the input file.
virtual SharedPtr< FunctionBase >::Type draw()=0
Draw from the measure, and then return a shared pointer to the draw.
virtual double uniformSample() const =0
Samples a value from a uniform distribution.
SharedPtr< FunctionBase >::Type proposed_physical_state
virtual double evaluate(FunctionBase &flow)=0
Evaluate the likelihood functional at flow. Subclasses must implement this method.
std::ofstream * subDisplayFile() const
Access function for m_subDisplayFile (displays file on stream).
Definition: Environment.C:320

Generated on Tue Jun 5 2018 19:48:55 for queso-0.57.1 by  doxygen 1.8.5