Skip to content

Commit

Permalink
Ports SequentialTemperingMCMC.py to v4
Browse files Browse the repository at this point in the history
  • Loading branch information
dimtsap committed Nov 25, 2022
1 parent e758876 commit 1b479bd
Show file tree
Hide file tree
Showing 2 changed files with 24 additions and 15 deletions.
28 changes: 16 additions & 12 deletions src/UQpy/sampling/tempering_mcmc/SequentialTemperingMCMC.py
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,8 @@ def __init__(self, pdf_intermediate=None, log_pdf_intermediate=None, args_pdf_in
# Initialize input distributions
self.evaluate_log_reference, self.seed = self._preprocess_reference(dist_=distribution_reference,
seed_=seed, nsamples=nsamples,
dimension=self.__dimension)
dimension=self.__dimension,
random_state=self.random_state)

# Initialize flag that indicates whether default proposal is to be used (default proposal defined adaptively
# during run)
Expand All @@ -95,12 +96,12 @@ def __init__(self, pdf_intermediate=None, log_pdf_intermediate=None, args_pdf_in

# Call the run function
if nsamples is not None:
self._run(nsamples=nsamples)
self.run(nsamples=nsamples)
else:
raise ValueError('UQpy: a value for "nsamples" must be specified ')

@beartype
def _run(self, nsamples: PositiveInteger = None):
def run(self, nsamples: PositiveInteger = None):

self.logger.info('TMCMC Start')

Expand Down Expand Up @@ -175,16 +176,17 @@ def _run(self, nsamples: PositiveInteger = None):
for i in range(self.n_resamples):

# Resampling from previous tempering level
lead_index = int(np.random.choice(pts_index, p=weight_probabilities))
lead_index = int(self.random_state.choice(pts_index, p=weight_probabilities))
lead = points[lead_index]

# Defining the default proposal
if self.proposal_given_flag is False:
self.proposal = MultivariateNormal(lead, cov=sigma_matrix)

# Single MH-MCMC step
x = MetropolisHastings(dimension=self.__dimension, log_pdf_target=mcmc_log_pdf_target, seed=lead,
nsamples=1, nchains=1, nburn=self.resampling_burn_length, proposal=self.proposal,
x = MetropolisHastings(dimension=self.__dimension, log_pdf_target=mcmc_log_pdf_target, seed=list(lead),
nsamples=1, n_chains=1, burn_length=self.resampling_burn_length,
proposal=self.proposal, random_state=self.random_state,
proposal_is_symmetric=self.proposal_is_symmetric)

# Setting the generated sample in the array
Expand All @@ -199,12 +201,14 @@ def _run(self, nsamples: PositiveInteger = None):
self.logger.info('Begin MCMC')
mcmc_seed = self._mcmc_seed_generator(resampled_pts=points[0:self.n_resamples, :],
arr_length=self.n_resamples,
seed_length=self.__n_chains)
seed_length=self.__n_chains,
random_state=self.random_state)

y = copy.deepcopy(self.sampler)
self.update_target_and_seed(y, mcmc_seed, mcmc_log_pdf_target)
y = self.sampler.__copy__(log_pdf_target=mcmc_log_pdf_target, seed=mcmc_seed,
nsamples_per_chain=self.n_samples_per_chain, concat_chains=True)
nsamples_per_chain=self.n_samples_per_chain,
concatenate_chains=True, random_state=self.random_state)
points[self.n_resamples:, :] = y.samples

if self.save_intermediate_samples is True:
Expand Down Expand Up @@ -288,7 +292,7 @@ def _find_temper_param(temper_param_prev, samples, q_func, n, iter_lim=1000, ite
raise RuntimeError('UQpy: unable to find tempering exponent due to nonconvergence')
return temper_param_trial

def _preprocess_reference(self, dist_, seed_=None, nsamples=None, dimension=None):
def _preprocess_reference(self, dist_, seed_=None, nsamples=None, dimension=None, random_state=None):
"""
Preprocess the target pdf inputs.
Expand Down Expand Up @@ -316,7 +320,7 @@ def _preprocess_reference(self, dist_, seed_=None, nsamples=None, dimension=None
if not (isinstance(dist_, Distribution)):
raise TypeError('UQpy: A UQpy.Distribution object must be provided.')
evaluate_log_pdf = (lambda x: dist_.log_pdf(x))
seed_values = dist_.rvs(nsamples=nsamples)
seed_values = dist_.rvs(nsamples=nsamples, random_state=random_state)
elif seed_ is not None:
if seed_.shape[0] != nsamples or seed_.shape[1] != dimension:
raise TypeError('UQpy: the seed values should be a numpy array of size (nsamples, dimension)')
Expand All @@ -328,7 +332,7 @@ def _preprocess_reference(self, dist_, seed_=None, nsamples=None, dimension=None
return evaluate_log_pdf, seed_values

@staticmethod
def _mcmc_seed_generator(resampled_pts, arr_length, seed_length):
def _mcmc_seed_generator(resampled_pts, arr_length, seed_length, random_state):
"""
Generates the seed from the resampled samples for the mcmc step
Expand All @@ -346,6 +350,6 @@ def _mcmc_seed_generator(resampled_pts, arr_length, seed_length):
* evaluate_log_pdf (callable): Callable that computes the log of the target density function (the prior)
"""
index_arr = np.arange(arr_length)
seed_indices = np.random.choice(index_arr, size=seed_length, replace=False)
seed_indices = random_state.choice(index_arr, size=seed_length, replace=False)
mcmc_seed = resampled_pts[seed_indices, :]
return mcmc_seed
11 changes: 8 additions & 3 deletions tests/unit_tests/sampling/test_tempering.py
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,12 @@ def likelihood(x, b):

def test_sequential():
prior = JointIndependent(marginals=[Uniform(loc=-2.0, scale=4.0), Uniform(loc=-2.0, scale=4.0)])
test = SequentialTemperingMCMC(dimension=2, nsamples=100, pdf_intermediate=likelihood,
distribution_reference=prior, nchains=20, save_intermediate_samples=True,
percentage_resampling=10, mcmc_class=MH, verbose=False, random_state=960242069)
sampler = MetropolisHastings(dimension=2, n_chains=20)
test = SequentialTemperingMCMC(pdf_intermediate=likelihood,
distribution_reference=prior,
save_intermediate_samples=True,
percentage_resampling=10,
random_state=960242069,
sampler=sampler,
nsamples=100)
assert np.round(test.evidence, 4) == 0.0656

0 comments on commit 1b479bd

Please sign in to comment.