Class representing the infinite dimensional Markov chain Monte Carlo sampler. More...

#include <InfiniteDimensionalMCMCSampler.h>

Public Member Functions
	InfiniteDimensionalMCMCSampler (const BaseEnvironment &env, InfiniteDimensionalMeasureBase &prior, InfiniteDimensionalLikelihoodBase &llhd, InfiniteDimensionalMCMCSamplerOptions *ov)
	Construct an infinite dimensional MCMC chain with given `prior`, likelihood `llhd`, and options `ov`. More...

	~InfiniteDimensionalMCMCSampler ()
	Destructor. More...

void	step ()
	Do one iteration of the Markov chain. More...

double	llhd_val () const
	Get the current value of the llhd. More...

double	acc_prob ()
	Returns current acceptance probability. More...

double	avg_acc_prob ()
	Returns current average acceptance probability. More...

unsigned int	iteration () const
	Returns the current iteration number. More...

SharedPtr < InfiniteDimensionalMCMCSampler > ::Type	clone_and_reset () const
	Returns a pointer to new sampler, with all the moments reset. More...

Private Member Functions
void	_propose ()

void	_metropolis_hastings ()

void	_update_moments ()

void	_write_state ()

hid_t	_create_scalar_dataset (const std::string &name)

void	_append_scalar_dataset (hid_t dset, double data)

Private Attributes
unsigned int	_iteration

double	_llhd_val

double	_acc_prob

double	_avg_acc_prob

InfiniteDimensionalMeasureBase &	prior

InfiniteDimensionalLikelihoodBase &	llhd

InfiniteDimensionalMCMCSamplerOptions *	m_ov

const BaseEnvironment &	m_env

SharedPtr< FunctionBase >::Type	current_physical_state

SharedPtr< FunctionBase >::Type	proposed_physical_state

SharedPtr< FunctionBase >::Type	current_physical_mean

SharedPtr< FunctionBase >::Type	current_physical_var

SharedPtr< FunctionBase >::Type	_delta

SharedPtr< FunctionBase >::Type	_M2

hid_t	_outfile

hid_t	_acc_dset

hid_t	_avg_acc_dset

hid_t	_neg_log_llhd_dset

hid_t	_L2_norm_samples_dset

hid_t	_L2_norm_mean_dset

hid_t	_L2_norm_var_dset

Detailed Description

Class representing the infinite dimensional Markov chain Monte Carlo sampler.

Definition at line 51 of file InfiniteDimensionalMCMCSampler.h.

Constructor & Destructor Documentation

QUESO::InfiniteDimensionalMCMCSampler::InfiniteDimensionalMCMCSampler	(	const BaseEnvironment &	env,
		InfiniteDimensionalMeasureBase &	prior,
		InfiniteDimensionalLikelihoodBase &	llhd,
		InfiniteDimensionalMCMCSamplerOptions *	ov
	)

Construct an infinite dimensional MCMC chain with given prior, likelihood llhd, and options ov.

Definition at line 42 of file InfiniteDimensionalMCMCSampler.C.

References _acc_dset, _acc_prob, _avg_acc_dset, _avg_acc_prob, _create_scalar_dataset(), _delta, _iteration, _L2_norm_mean_dset, _L2_norm_samples_dset, _L2_norm_var_dset, _llhd_val, _M2, _neg_log_llhd_dset, _outfile, QUESO::CheckFilePath(), current_physical_mean, current_physical_state, current_physical_var, QUESO::InfiniteDimensionalLikelihoodBase::evaluate(), QUESO::InfiniteDimensionalMCMCSamplerOptions::m_dataOutputDirName, QUESO::InfiniteDimensionalMCMCSamplerOptions::m_dataOutputFileName, m_env, m_ov, QUESO::InfiniteDimensionalMCMCSamplerOptions::m_prefix, QUESO::BaseEnvironment::optionsInputFileName(), and QUESO::BaseEnvironment::subDisplayFile().

Referenced by clone_and_reset().

   : prior(prior),
     llhd(llhd),
     m_env(env),
     current_physical_state(prior.draw()),
     proposed_physical_state(prior.draw()),
     current_physical_mean(prior.draw()),
     current_physical_var(prior.draw()),
     _delta(prior.draw()),
     _M2(prior.draw())
 {
   if (ov != NULL) {
     this->m_ov = ov;
   }
 
 #ifdef QUESO_MEMORY_DEBUGGING
   std::cout << "Entering InfiniteDimensionalMCMCSampler class" << std::endl;
 #endif
   if (m_env.subDisplayFile()) {
     *m_env.subDisplayFile() << "Entering InfiniteDimensionalMCMCSampler::constructor()"
                             << ": prefix = " << this->m_ov->m_prefix
                             << ", m_env.optionsInputFileName() = " << this->m_env.optionsInputFileName()
                             << std::endl;
   }
 
 #ifdef QUESO_MEMORY_DEBUGGING
   std::cout << "In InfiniteDimensionalMCMCSamplerOptions,"
             << " finished scanning options" << std::endl;
 #endif
 
   // Verify parent directory exists (for cases when a user specifies a
   // relative path for the desired output file).
   // Only subprocess of subrank 0 creates the output file
   if ((this->m_env).subRank() == 0) {
     int irtrn = CheckFilePath((this->m_ov->m_dataOutputDirName +
                                (this->m_env).subIdString() +
                                "/test.txt").c_str());
     queso_require_greater_equal_msg(irtrn, 0, "unable to verify output path");
   }
 
   // Only subprocess of subrank 0 creates the output file
   if ((this->m_env).subRank() == 0) {
     this->_outfile = H5Fcreate((this->m_ov->m_dataOutputDirName +
           (this->m_env).subIdString() + "/" +
           this->m_ov->m_dataOutputFileName).c_str(),
         H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
   }
 
   this->_acc_dset = this->_create_scalar_dataset("acc");
   this->_avg_acc_dset = this->_create_scalar_dataset("avg_acc");
   this->_neg_log_llhd_dset = this->_create_scalar_dataset("neg_log_llhd");
   this->_L2_norm_samples_dset = this->_create_scalar_dataset("L2_norm_samples");
   this->_L2_norm_mean_dset = this->_create_scalar_dataset("L2_norm_mean");
   this->_L2_norm_var_dset = this->_create_scalar_dataset("L2_norm_var");
 
   // Ensure that we created the path and the output files
   ((this->m_env).fullComm()).Barrier();
 
   this->_iteration = 0;
   this->_acc_prob = 0.0;
   this->_avg_acc_prob = 0.0;
 
   // Zero out these guys.  There's probably a better way of doing this than
   // calling prior.draw() at the start, but creation of a Sampler object
   // should only be done a O(1) times anyway.
   this->_delta->zero();
   this->_M2->zero();
   this->current_physical_mean->zero();
   this->current_physical_var->zero();
 
   // LibMeshFunction & p = libmesh_cast_ref<LibMeshFunction &>(*(this->current_physical_state));
 
   // Seed the sampler at the truth
   // for (unsigned int ii = 0; ii < 513; ii++) {
   //   p.equation_systems->get_system<ExplicitSystem>("Function").solution->set(ii, std::cos(2.0 * boost::math::constants::pi<double>() * ii / 512.0));
   // }
   // p.equation_systems->get_system<ExplicitSystem>("Function").solution->close();
   // std::string seedname("seedfn");
   // p.save_function(seedname);
 
   // Initialise cost function
   this->_llhd_val = this->llhd.evaluate(*(this->current_physical_state));
 
   std::cout.setf(std::ios::fixed, std::ios::floatfield);
   std::cout.precision(15);
 }

QUESO::InfiniteDimensionalMCMCSampler::~InfiniteDimensionalMCMCSampler ( )

Destructor.

Definition at line 133 of file InfiniteDimensionalMCMCSampler.C.

References _acc_dset, _avg_acc_dset, _L2_norm_mean_dset, _L2_norm_samples_dset, _L2_norm_var_dset, _neg_log_llhd_dset, _outfile, and m_env.

 {
   if ((this->m_env).subRank() == 0) {
     H5Dclose(this->_acc_dset);
     H5Dclose(this->_avg_acc_dset);
     H5Dclose(this->_neg_log_llhd_dset);
     H5Dclose(this->_L2_norm_samples_dset);
     H5Dclose(this->_L2_norm_mean_dset);
     H5Dclose(this->_L2_norm_var_dset);
     H5Fclose(this->_outfile);
   }
 }

Member Function Documentation

void QUESO::InfiniteDimensionalMCMCSampler::_append_scalar_dataset	(	hid_t	dset,
		double	data
	)

private

Definition at line 268 of file InfiniteDimensionalMCMCSampler.C.

References _iteration, m_env, m_ov, and QUESO::InfiniteDimensionalMCMCSamplerOptions::m_save_freq.

Referenced by _write_state().

 {
   // Only subprocess with rank 0 manipulates the output file
   if ((this->m_env).subRank() == 0) {
     int err;
     // Create a memory dataspace for data to append
     const int ndims = 1;
     hsize_t dims[ndims] = { 1 };  // Only writing one double
     hid_t mem_space = H5Screate_simple(ndims, dims, NULL);
 
     // Extend the dataset
     // Set dims to be the *new* dimension of the extended dataset
     dims[0] = this->_iteration / this->m_ov->m_save_freq;
     err = H5Dset_extent(dset, dims);
     queso_require_greater_equal_msg(err, 0, "H5DSet_extent(dset, dims) failed");
 
     // Select hyperslab on file dataset
     hid_t file_space = H5Dget_space(dset);
     hsize_t start[1] = {(this->_iteration / this->m_ov->m_save_freq) - 1};
     hsize_t count[1] = {1};
 
     err = H5Sselect_hyperslab(file_space, H5S_SELECT_SET, start, NULL, count, NULL);
     queso_require_greater_equal_msg(err, 0, "H5Sselect_hyperslab(file_space, H5S_SELECT_SET, start, NULL, count, NULL) failed");
 
     // hsize_t      size[1];
     // size[0]   = this->_iteration / this->m_ov->m_save_freq;
 
     // Write the data
     H5Dwrite(dset, H5T_NATIVE_DOUBLE, mem_space, file_space, H5P_DEFAULT, &data);
 
     // Close a bunch of stuff
     H5Sclose(file_space);
     H5Sclose(mem_space);
   }
 }

hid_t QUESO::InfiniteDimensionalMCMCSampler::_create_scalar_dataset ( const std::string & name )

private

Definition at line 231 of file InfiniteDimensionalMCMCSampler.C.

References _outfile, and m_env.

Referenced by InfiniteDimensionalMCMCSampler().

 {
   // Only subprocess with rank 0 manipulates the output file
   if ((this->m_env).subRank() == 0) {
     // Create a 1D dataspace.  Unlimited size.  Initially set to 0.  We will
     // extend it later
     const int ndims = 1;
     hsize_t dims[ndims] = {0};  // dataset dimensions at creation
     hsize_t maxdims[ndims] = {H5S_UNLIMITED};
 
     hid_t file_space = H5Screate_simple(ndims, dims, maxdims);
 
     // Create dataset creation property list.  Unlimited datasets must be
     // chunked.  Choosing the chunk size is an issue, here we set it to 1
     hsize_t chunk_dims[ndims] = {1};
     hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
 
     H5Pset_layout(plist, H5D_CHUNKED);
     H5Pset_chunk(plist, ndims, chunk_dims);
 
     // Create the dataset
     hid_t dset = H5Dcreate(this->_outfile, name.c_str(), H5T_NATIVE_DOUBLE,
         file_space, H5P_DEFAULT, plist, H5P_DEFAULT);
 
     // We don't need the property list anymore.  We also don't need the file
     // dataspace anymore because we'll extend it later, making this one
     // invalild anyway.
     H5Pclose(plist);
     H5Sclose(file_space);
 
     return dset;
   }
 
   hid_t dummy = -1;
   return dummy;
 }

void QUESO::InfiniteDimensionalMCMCSampler::_metropolis_hastings ( )

private

Definition at line 158 of file InfiniteDimensionalMCMCSampler.C.

References _acc_prob, _llhd_val, current_physical_state, QUESO::InfiniteDimensionalLikelihoodBase::evaluate(), llhd, m_env, proposed_physical_state, QUESO::BaseEnvironment::rngObject(), and QUESO::RngBase::uniformSample().

Referenced by step().

 {
   // Downcast since we know we're dealing with a libmesh function
   // LibMeshFunction & p = static_cast<LibMeshFunction &>(*(this->proposed_physical_state));
   // LibMeshFunction & q = static_cast<LibMeshFunction &>(*(this->current_physical_state));
   // Evaluate the likelihood at the proposed state
   double proposed_llhd = this->llhd.evaluate(*(this->proposed_physical_state));
   // double current_llhd = this->llhd.evaluate(q);
   // std::cout << "llhd of  current is: " << current_llhd << std::endl;
   // std::cout << "llhd of proposal is: " << proposed_llhd << std::endl;
 
   double diff = this->_llhd_val - proposed_llhd;
   double alpha = std::min(1.0, std::exp(diff));
   double rand = this->m_env.rngObject()->uniformSample();
   if (rand < alpha) {
     // Accept
     // std::cout << "accepted" << std::endl;
     this->current_physical_state.swap(this->proposed_physical_state);
     this->proposed_physical_state->zero();
     this->_acc_prob = alpha;
     this->_llhd_val = proposed_llhd;
   }
   else {
     // std::cout << "rejected" << std::endl;
     this->_acc_prob = 0.0;
   }
 }

void QUESO::InfiniteDimensionalMCMCSampler::_propose ( )

private

Make a proposal from the prior using a standard random walk

Definition at line 146 of file InfiniteDimensionalMCMCSampler.C.

References current_physical_state, QUESO::InfiniteDimensionalMeasureBase::draw(), m_ov, QUESO::InfiniteDimensionalMCMCSamplerOptions::m_rwmh_step, prior, and proposed_physical_state.

Referenced by step().

 {
   const double rwmh_step_sq = (this->m_ov->m_rwmh_step * this->m_ov->m_rwmh_step);
   const double coeff = std::sqrt(1.0 - rwmh_step_sq);
 
   SharedPtr<FunctionBase>::Type p(prior.draw());
 
   this->proposed_physical_state->zero();
   this->proposed_physical_state->add(coeff, *(this->current_physical_state));
   this->proposed_physical_state->add(this->m_ov->m_rwmh_step, *p);
 }

void QUESO::InfiniteDimensionalMCMCSampler::_update_moments ( )

private

Definition at line 186 of file InfiniteDimensionalMCMCSampler.C.

References _avg_acc_prob, _delta, _iteration, _M2, acc_prob(), avg_acc_prob(), current_physical_mean, current_physical_state, current_physical_var, and iteration().

Referenced by step().

 {
   // Increment the current iteration number and update the running mean and
   // variance.
   this->_iteration += 1;
 
   // Make _delta contain the difference from the mean
   this->_delta->zero();
   this->_delta->add(1.0, *(this->current_physical_state));
   this->_delta->add(-1.0, *(this->current_physical_mean));
 
   // Scale _delta to update the mean field
   this->current_physical_mean->add(1.0 / this->iteration(), *(this->_delta));
 
   // Update running sum-of-squares
   SharedPtr<FunctionBase>::Type temp_ptr(this->_delta->zero_clone());
   // LibMeshFunction & temp = static_cast<LibMeshFunction &>(*temp_ptr);
 
   temp_ptr->pointwise_mult(*(this->_delta), *(this->current_physical_state));
   this->_M2->add(1.0, *temp_ptr);
   temp_ptr->pointwise_mult(*(this->_delta), *(this->current_physical_mean));
   this->_M2->add(-1.0, *temp_ptr);
 
   if (this->iteration() > 1) {
     this->current_physical_var->zero();
     this->current_physical_var->add(1.0 / (this->iteration() - 1), *(this->_M2));
   }
 
   // Update acceptance rate
   double delta_acc = this->acc_prob() - this->avg_acc_prob();
   this->_avg_acc_prob += delta_acc / this->iteration();
 }

void QUESO::InfiniteDimensionalMCMCSampler::_write_state ( )

private

Definition at line 304 of file InfiniteDimensionalMCMCSampler.C.

References _acc_dset, _append_scalar_dataset(), _avg_acc_dset, _L2_norm_mean_dset, _L2_norm_samples_dset, _L2_norm_var_dset, _llhd_val, _neg_log_llhd_dset, acc_prob(), avg_acc_prob(), current_physical_mean, current_physical_state, current_physical_var, iteration(), QUESO::InfiniteDimensionalMCMCSamplerOptions::m_dataOutputDirName, m_env, and m_ov.

Referenced by step().

 {
   this->_append_scalar_dataset(this->_acc_dset, this->acc_prob());
   this->_append_scalar_dataset(this->_avg_acc_dset, this->avg_acc_prob());
   this->_append_scalar_dataset(this->_neg_log_llhd_dset, this->_llhd_val);
   this->_append_scalar_dataset(this->_L2_norm_samples_dset, this->current_physical_state->L2_norm());
   this->_append_scalar_dataset(this->_L2_norm_mean_dset, this->current_physical_mean->L2_norm());
   this->_append_scalar_dataset(this->_L2_norm_var_dset, this->current_physical_var->L2_norm());
 
   // Now to write the fields.  It appears to be a pain in the arse to write a
   // method to spit this out to HDF5 format.  Also, this won't scale to
   // non-uniform finite element meshes.  Therefore, I'm going to spit out an
   // ExodusII file for each sample, average and variance.  Got disk space?
   std::stringstream basename;
   basename << this->m_ov->m_dataOutputDirName;
   basename << (this->m_env).subIdString();
   basename << "/";  // Sigh
 
   std::ostringstream curr_iter;
   curr_iter << this->iteration();
 
   std::string sample_name(basename.str() + "sample.e-s.");
   sample_name += curr_iter.str();
   this->current_physical_state->save_function(sample_name, this->iteration());
 
   std::string mean_name(basename.str() + "mean.e-s.");
   mean_name += curr_iter.str();
   this->current_physical_mean->save_function(mean_name, this->iteration());
 
   std::string var_name(basename.str() + "var.e-s.");
   var_name += curr_iter.str();
   this->current_physical_var->save_function(var_name, this->iteration());
 }

double QUESO::InfiniteDimensionalMCMCSampler::acc_prob ( )

Returns current acceptance probability.

Definition at line 353 of file InfiniteDimensionalMCMCSampler.C.

References _acc_prob.

Referenced by _update_moments(), and _write_state().

 {
   return this->_acc_prob;
 }

double QUESO::InfiniteDimensionalMCMCSampler::avg_acc_prob ( )

Returns current average acceptance probability.

Definition at line 358 of file InfiniteDimensionalMCMCSampler.C.

References _avg_acc_prob.

Referenced by _update_moments(), and _write_state().

 {
   return this->_avg_acc_prob;
 }

SharedPtr< InfiniteDimensionalMCMCSampler >::Type QUESO::InfiniteDimensionalMCMCSampler::clone_and_reset ( ) const

Returns a pointer to new sampler, with all the moments reset.

Definition at line 338 of file InfiniteDimensionalMCMCSampler.C.

References _acc_prob, _llhd_val, current_physical_state, InfiniteDimensionalMCMCSampler(), llhd, m_env, m_ov, prior, and proposed_physical_state.

 {
   // Set up a clone
   SharedPtr<InfiniteDimensionalMCMCSampler>::Type clone(new InfiniteDimensionalMCMCSampler(this->m_env, this->prior, this->llhd, this->m_ov));
 
   // Copy the state.
   clone->current_physical_state = this->current_physical_state;
   clone->proposed_physical_state = this->proposed_physical_state;
   clone->_llhd_val = this->_llhd_val;
   clone->_acc_prob = this->_acc_prob;
   clone->_avg_acc_prob = 0.0;
 
   return clone;
 }

unsigned int QUESO::InfiniteDimensionalMCMCSampler::iteration ( ) const

Returns the current iteration number.

Definition at line 368 of file InfiniteDimensionalMCMCSampler.C.

References _iteration.

Referenced by _update_moments(), and _write_state().

 {
   return this->_iteration;
 }

double QUESO::InfiniteDimensionalMCMCSampler::llhd_val ( ) const

Get the current value of the llhd.

Definition at line 363 of file InfiniteDimensionalMCMCSampler.C.

References _llhd_val.

 {
   return this->_llhd_val;
 }

void QUESO::InfiniteDimensionalMCMCSampler::step ( )

Do one iteration of the Markov chain.

Definition at line 219 of file InfiniteDimensionalMCMCSampler.C.

References _iteration, _metropolis_hastings(), _propose(), _update_moments(), _write_state(), m_ov, and QUESO::InfiniteDimensionalMCMCSamplerOptions::m_save_freq.

 {
   this->_propose();
   this->_metropolis_hastings();
   this->_update_moments();
 
   // We never save the 0th iteration
   if (this->_iteration % this->m_ov->m_save_freq == 0) {
     this->_write_state();
   }
 }