From 0f0ae2ee68b0550b47e741cf9c87977b7c51f6f2 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Wed, 18 Dec 2024 19:09:57 +0000 Subject: [PATCH 1/5] Updated docs for mcmc methods to indicate if they use sensitivities and make sure all mention single/multi chain --- pints/_mcmc/__init__.py | 9 ++++++--- pints/_mcmc/_adaptive_covariance.py | 2 +- pints/_mcmc/_differential_evolution.py | 2 +- pints/_mcmc/_dram_ac.py | 2 +- pints/_mcmc/_dream.py | 2 +- pints/_mcmc/_emcee_hammer.py | 2 ++ pints/_mcmc/_haario_ac.py | 2 +- pints/_mcmc/_haario_bardenet_ac.py | 2 +- pints/_mcmc/_hamiltonian.py | 7 +++++-- pints/_mcmc/_mala.py | 8 ++++++-- pints/_mcmc/_metropolis.py | 2 +- pints/_mcmc/_monomial_gamma_hamiltonian.py | 8 ++++++-- pints/_mcmc/_nuts.py | 8 ++++++-- pints/_mcmc/_population.py | 2 +- pints/_mcmc/_rao_blackwell_ac.py | 2 +- pints/_mcmc/_relativistic.py | 8 ++++++-- pints/_mcmc/_slice_doubling.py | 2 +- pints/_mcmc/_slice_rank_shrinking.py | 8 ++++++-- pints/_mcmc/_slice_stepout.py | 2 +- 19 files changed, 54 insertions(+), 26 deletions(-) diff --git a/pints/_mcmc/__init__.py b/pints/_mcmc/__init__.py index ff566bf54..cb8d5cf3e 100644 --- a/pints/_mcmc/__init__.py +++ b/pints/_mcmc/__init__.py @@ -735,6 +735,8 @@ def run(self): if reply is not None: # Unpack reply into position, evaluation, and status + # Note that evaluation is (fy, dfy) with sensitivities + # enabled. y, fy, accepted = reply # Inverse transform to model space if transformation is @@ -1034,9 +1036,10 @@ def set_log_pdf_storage(self, store_in_memory=False): Store :class:`LogPDF` evaluations in memory as they are generated. By default, evaluations of the :class:`LogPDF` are not stored. This - method can be used to enable storage of the evaluations for the - accepted samples. - After running, evaluations can be obtained using :meth:`evaluations()`. + method can be used to enable (or disable) storage for the accepted + samples, in memory. + + After running, evaluations can be obtained using :meth:`log_pdfs()`. """ self._evaluations_in_memory = bool(store_in_memory) diff --git a/pints/_mcmc/_adaptive_covariance.py b/pints/_mcmc/_adaptive_covariance.py index 49d1c9545..c895daf06 100644 --- a/pints/_mcmc/_adaptive_covariance.py +++ b/pints/_mcmc/_adaptive_covariance.py @@ -23,7 +23,7 @@ class AdaptiveCovarianceMC(pints.SingleChainMCMC): A basic implementation is provided for each, which extending methods can choose to override. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, does not use sensitivities. """ def __init__(self, x0, sigma0=None): diff --git a/pints/_mcmc/_differential_evolution.py b/pints/_mcmc/_differential_evolution.py index 8f8f4e9cb..feb87c673 100644 --- a/pints/_mcmc/_differential_evolution.py +++ b/pints/_mcmc/_differential_evolution.py @@ -29,7 +29,7 @@ class DifferentialEvolutionMCMC(pints.MultiChainMCMC): If ``x_proposed / x[i,r] > u ~ U(0,1)``, then ``x[i+1,r] = x_proposed``; otherwise, ``x[i+1,r] = x[i]``. - Extends :class:`MultiChainMCMC`. + Extends :class:`MultiChainMCMC`, does not use sensitivities. .. note:: This sampler requires a number of chains :math:`n \ge 3`, and diff --git a/pints/_mcmc/_dram_ac.py b/pints/_mcmc/_dram_ac.py index 2d7d2b83e..e3d43144c 100644 --- a/pints/_mcmc/_dram_ac.py +++ b/pints/_mcmc/_dram_ac.py @@ -50,7 +50,7 @@ class DramACMC(pints.AdaptiveCovarianceMC): are then adapted as ``[scale_1, scale_2] * sigma``, where the scale factors are set using ``set_sigma_scale``. - *Extends:* :class:`GlobalAdaptiveCovarianceMC` + *Extends:* :class:`GlobalAdaptiveCovarianceMC`, does not use sensitivities. References ---------- diff --git a/pints/_mcmc/_dream.py b/pints/_mcmc/_dream.py index 99e713c24..619476db0 100644 --- a/pints/_mcmc/_dream.py +++ b/pints/_mcmc/_dream.py @@ -53,7 +53,7 @@ class DreamMCMC(pints.MultiChainMCMC): Here b > 0, b* > 0, 1 >= p_g >= 0, 1 >= CR >= 0. - Extends :class:`MultiChainMCMC`. + Extends :class:`MultiChainMCMC`, does not use sensitivities. References ---------- diff --git a/pints/_mcmc/_emcee_hammer.py b/pints/_mcmc/_emcee_hammer.py index 9cc7c0985..111ed6e30 100644 --- a/pints/_mcmc/_emcee_hammer.py +++ b/pints/_mcmc/_emcee_hammer.py @@ -34,6 +34,8 @@ class EmceeHammerMCMC(pints.MultiChainMCMC): ``1 / sqrt(z)`` if ``z`` is in ``[1 / a, a]`` or to 0, otherwise (where ``a`` is a parameter with default value ``2``). + Extends :class:`MultiChainMCMC`, does not use sensitivities. + References ---------- .. [1] "emcee: The MCMC Hammer", Daniel Foreman-Mackey, David W. Hogg, diff --git a/pints/_mcmc/_haario_ac.py b/pints/_mcmc/_haario_ac.py index f2770f15a..33bafb84c 100644 --- a/pints/_mcmc/_haario_ac.py +++ b/pints/_mcmc/_haario_ac.py @@ -44,7 +44,7 @@ class HaarioACMC(pints.AdaptiveCovarianceMC): log lambda = log lambda + gamma (alpha - self._target_acceptance) gamma = adaptation_count^-eta - Extends :class:`AdaptiveCovarianceMC`. + Extends :class:`AdaptiveCovarianceMC`, does not use sensitivities. References ---------- diff --git a/pints/_mcmc/_haario_bardenet_ac.py b/pints/_mcmc/_haario_bardenet_ac.py index 520bb74a7..07fde6622 100644 --- a/pints/_mcmc/_haario_bardenet_ac.py +++ b/pints/_mcmc/_haario_bardenet_ac.py @@ -42,7 +42,7 @@ class HaarioBardenetACMC(pints.AdaptiveCovarianceMC): log lambda = log lambda + gamma (alpha - self._target_acceptance) gamma = adaptation_count^-eta - Extends :class:`AdaptiveCovarianceMC`. + Extends :class:`AdaptiveCovarianceMC`, does not use sensitivities. References ---------- diff --git a/pints/_mcmc/_hamiltonian.py b/pints/_mcmc/_hamiltonian.py index 9ed5dbb97..5840c3dca 100644 --- a/pints/_mcmc/_hamiltonian.py +++ b/pints/_mcmc/_hamiltonian.py @@ -46,7 +46,7 @@ class HamiltonianMCMC(pints.SingleChainMCMC): In particular, the algorithm we implement follows eqs. (4.14)-(4.16) in [1]_, since we allow different epsilon according to dimension. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, this method uses sensitivities. References ---------- @@ -190,7 +190,10 @@ def name(self): return 'Hamiltonian Monte Carlo' def needs_sensitivities(self): - """ See :meth:`pints.MCMCSampler.needs_sensitivities()`. """ + """ + This method requires sensitivities. + + See :meth:`pints.MCMCSampler.needs_sensitivities()`. """ return True def scaled_epsilon(self): diff --git a/pints/_mcmc/_mala.py b/pints/_mcmc/_mala.py index 0c05c3ada..f893ef16d 100644 --- a/pints/_mcmc/_mala.py +++ b/pints/_mcmc/_mala.py @@ -67,7 +67,7 @@ class MALAMCMC(pints.SingleChainMCMC): where :math:`\epsilon' = \epsilon \sqrt{M}` is given by the initial value of ``sigma0``. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, this method uses sensitivities. References ---------- @@ -181,7 +181,11 @@ def name(self): return 'Metropolis-Adjusted Langevin Algorithm (MALA)' def needs_sensitivities(self): - """ See :meth:`pints.MCMCSampler.needs_sensitivities()`. """ + """ + This method requires sensitivities. + + See :meth:`pints.MCMCSampler.needs_sensitivities()`. + """ return True def n_hyper_parameters(self): diff --git a/pints/_mcmc/_metropolis.py b/pints/_mcmc/_metropolis.py index 032e86d40..50871cc57 100644 --- a/pints/_mcmc/_metropolis.py +++ b/pints/_mcmc/_metropolis.py @@ -24,7 +24,7 @@ class MetropolisRandomWalkMCMC(pints.SingleChainMCMC): here Sigma is the covariance matrix of the proposal. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, does not use sensitivities. References ---------- diff --git a/pints/_mcmc/_monomial_gamma_hamiltonian.py b/pints/_mcmc/_monomial_gamma_hamiltonian.py index d5e497481..cc0a6a429 100644 --- a/pints/_mcmc/_monomial_gamma_hamiltonian.py +++ b/pints/_mcmc/_monomial_gamma_hamiltonian.py @@ -63,7 +63,7 @@ class MonomialGammaHamiltonianMCMC(pints.SingleChainMCMC): In particular, the algorithm we implement follows eqs. (4.14)-(4.16) in [2]_, since we allow different epsilon according to dimension. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, this method uses sensitivities. References ---------- @@ -297,7 +297,11 @@ def name(self): return 'Monomial-Gamma Hamiltonian Monte Carlo' def needs_sensitivities(self): - """ See :meth:`pints.MCMCSampler.needs_sensitivities()`. """ + """ + This method requires sensitivities. + + See also :meth:`pints.MCMCSampler.needs_sensitivities()`. + """ return True def n_hyper_parameters(self): diff --git a/pints/_mcmc/_nuts.py b/pints/_mcmc/_nuts.py index b665b3469..9e5e71102 100644 --- a/pints/_mcmc/_nuts.py +++ b/pints/_mcmc/_nuts.py @@ -499,7 +499,7 @@ class NoUTurnMCMC(pints.SingleChainMCMC): Note: This sampler is only supported on Python versions 3.3 and newer. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, this method uses sensitivities. References ---------- @@ -664,7 +664,11 @@ def name(self): return 'No-U-Turn MCMC' def needs_sensitivities(self): - """ See :meth:`pints.MCMCSampler.needs_sensitivities()`. """ + """ + This method requires sensitivities. + + See :meth:`pints.MCMCSampler.needs_sensitivities()`. + """ return True def number_adaption_steps(self): diff --git a/pints/_mcmc/_population.py b/pints/_mcmc/_population.py index 53eafdc7b..be5a8d5d6 100644 --- a/pints/_mcmc/_population.py +++ b/pints/_mcmc/_population.py @@ -39,7 +39,7 @@ class PopulationMCMC(pints.SingleChainMCMC): ``delta_T = 1 / num_temperatures``, and the chain with ``T_i = 0`` is the one whose target distribution we want to sample. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, does not use sensitivities. References ---------- diff --git a/pints/_mcmc/_rao_blackwell_ac.py b/pints/_mcmc/_rao_blackwell_ac.py index 4b6b90bec..69876c690 100644 --- a/pints/_mcmc/_rao_blackwell_ac.py +++ b/pints/_mcmc/_rao_blackwell_ac.py @@ -37,7 +37,7 @@ class RaoBlackwellACMC(pints.AdaptiveCovarianceMC): Y_t+1 ~ N(theta_t, lambda * sigma0) rather than Y_t+1 ~ N(theta_t, sigma0) - Extends :class:`AdaptiveCovarianceMC`. + Extends :class:`AdaptiveCovarianceMC`, does not use sensitivities. References ---------- diff --git a/pints/_mcmc/_relativistic.py b/pints/_mcmc/_relativistic.py index adff93555..2c4b73441 100644 --- a/pints/_mcmc/_relativistic.py +++ b/pints/_mcmc/_relativistic.py @@ -57,7 +57,7 @@ class RelativisticMCMC(pints.SingleChainMCMC): In particular, the algorithm we implement follows eqs. in section 2.1 of [1]_. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, this method uses sensitivities. References ---------- @@ -325,7 +325,11 @@ def name(self): return 'Relativistic MCMC' def needs_sensitivities(self): - """ See :meth:`pints.MCMCSampler.needs_sensitivities()`. """ + """ + This method requires sensitivities. + + See :meth:`pints.MCMCSampler.needs_sensitivities()`. + """ return True def _sample_momentum(self): diff --git a/pints/_mcmc/_slice_doubling.py b/pints/_mcmc/_slice_doubling.py index c55bb0bc8..62565ede6 100644 --- a/pints/_mcmc/_slice_doubling.py +++ b/pints/_mcmc/_slice_doubling.py @@ -109,7 +109,7 @@ class SliceDoublingMCMC(pints.SingleChainMCMC): :math:`\epsilon \sim \text{exp}(1)` and define the slice as :math:`S = {x : z < g(x)}`. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, does not use sensitivities. References ---------- diff --git a/pints/_mcmc/_slice_rank_shrinking.py b/pints/_mcmc/_slice_rank_shrinking.py index 1716251f9..dfd52c2b8 100644 --- a/pints/_mcmc/_slice_rank_shrinking.py +++ b/pints/_mcmc/_slice_rank_shrinking.py @@ -63,7 +63,7 @@ class SliceRankShrinkingMCMC(pints.SingleChainMCMC): :math:`\epsilon \sim \text{exp}(1)` and define the slice as :math:`S = {x : z < log f(x)}`. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, this method uses sensitivities. References ---------- @@ -156,7 +156,11 @@ def name(self): return 'Slice Sampling - Covariance-Adaptive: Rank Shrinking.' def needs_sensitivities(self): - """ See :meth:`pints.MCMCSampler.needs_sensitivities()`. """ + """ + This method requires sensitivities. + + See :meth:`pints.MCMCSampler.needs_sensitivities()`. + """ return True def n_hyper_parameters(self): diff --git a/pints/_mcmc/_slice_stepout.py b/pints/_mcmc/_slice_stepout.py index ed6e8fa63..e675910b9 100644 --- a/pints/_mcmc/_slice_stepout.py +++ b/pints/_mcmc/_slice_stepout.py @@ -120,7 +120,7 @@ class SliceStepoutMCMC(pints.SingleChainMCMC): :math:`\epsilon \sim \text{exp}(1)` and define the slice as :math:`S = {x : z < g(x)}`. - Extends :class:`SingleChainMCMC`. + Extends :class:`SingleChainMCMC`, does not use sensitivities. References ---------- From 3d16d1a69b8bd341619fc2c127114512ae3230b7 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Wed, 18 Dec 2024 20:23:28 +0000 Subject: [PATCH 2/5] Added tests for #1691. --- pints/tests/test_mcmc_controller.py | 54 ++++++++++++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/pints/tests/test_mcmc_controller.py b/pints/tests/test_mcmc_controller.py index 98220805b..8d2a462e7 100755 --- a/pints/tests/test_mcmc_controller.py +++ b/pints/tests/test_mcmc_controller.py @@ -690,6 +690,26 @@ def test_log_pdf_storage_in_memory_single(self): likelihoods = [self.log_likelihood(x) for x in chain] self.assertTrue(np.all(evals[i] == likelihoods)) + # Test with a sensitivity-using method + mcmc = pints.MCMCController( + self.log_posterior, n_chains, xs, method=pints.MALAMCMC) + mcmc.set_max_iterations(n_iterations) + mcmc.set_log_to_screen(False) + mcmc.set_log_pdf_storage(True) + chains = mcmc.run() + evals = mcmc.log_pdfs() + self.assertEqual(len(evals.shape), 3) + self.assertEqual(evals.shape[0], n_chains) + self.assertEqual(evals.shape[1], n_iterations) + self.assertEqual(evals.shape[2], 3) + for i, chain in enumerate(chains): + posteriors = [self.log_posterior(x) for x in chain] + self.assertTrue(np.all(evals[i, :, 0] == posteriors)) + likelihoods = [self.log_likelihood(x) for x in chain] + self.assertTrue(np.all(evals[i, :, 1] == likelihoods)) + priors = [self.log_prior(x) for x in chain] + self.assertTrue(np.all(evals[i, :, 2] == priors)) + # Test disabling again mcmc = pints.MCMCController(self.log_posterior, n_chains, xs) mcmc.set_max_iterations(n_iterations) @@ -706,7 +726,9 @@ def test_log_pdf_storage_in_memory_multi(self): x0 = np.array(self.real_parameters) * 1.05 x1 = np.array(self.real_parameters) * 1.15 x2 = np.array(self.real_parameters) * 0.95 - xs = [x0, x1, x2] + x3 = np.array(self.real_parameters) * 0.95 + x4 = np.array(self.real_parameters) * 0.95 + xs = [x0, x1, x2, x3, x4] n_chains = len(xs) n_iterations = 100 meth = pints.DifferentialEvolutionMCMC @@ -737,6 +759,16 @@ def test_log_pdf_storage_in_memory_multi(self): priors = [self.log_prior(x) for x in chain] self.assertTrue(np.all(evals[i, :, 2] == priors)) + # Test with a sensitivity-using method + # We don't have any of these! + mcmc = pints.MCMCController( + self.log_posterior, n_chains, xs, + method=FakeMultiChainSamplerWithSensitivities) + mcmc.set_max_iterations(5) + mcmc.set_log_to_screen(False) + mcmc.set_log_pdf_storage(True) + chains = mcmc.run() + # Test with a loglikelihood mcmc = pints.MCMCController( self.log_likelihood, n_chains, xs, method=meth) @@ -1651,6 +1683,26 @@ def tell(self, fx): return None if x is None else (x, fx, np.array([True] * self._n)) +class FakeMultiChainSamplerWithSensitivities(pints.MultiChainMCMC): + """ + Fake sampler that pretends to be a multi-chain method that uses + sensitivities. At the moment (2024-12-18) we don't have these in PINTS, but + we need to check that code potentially handling this type of sampler exists + or raises not-implemented errors. + """ + def ask(self): + self._xs = self._x0 + 1e-3 * np.random.uniform( + size=(self._n_chains, self._n_parameters)) + return self._xs + + def current_log_pdfs(self): + return self._fxs + + def tell(self, fxs): + self._fxs = fxs + return self._xs, self._fxs, [True] * self._n_chains + + class TestMCMCControllerSingleChainStorage(unittest.TestCase): """ Tests storage of samples and evaluations to disk, running with a From c4c76a6005c162e55d7fafc0a49f4d0ed7c5e47e Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Wed, 18 Dec 2024 23:44:15 +0000 Subject: [PATCH 3/5] Fixed test for 1691 multi-chain --- pints/tests/test_mcmc_controller.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pints/tests/test_mcmc_controller.py b/pints/tests/test_mcmc_controller.py index 8d2a462e7..c2446fdc8 100755 --- a/pints/tests/test_mcmc_controller.py +++ b/pints/tests/test_mcmc_controller.py @@ -768,6 +768,11 @@ def test_log_pdf_storage_in_memory_multi(self): mcmc.set_log_to_screen(False) mcmc.set_log_pdf_storage(True) chains = mcmc.run() + evals = mcmc.log_pdfs() + self.assertEqual(len(evals.shape), 3) + self.assertEqual(evals.shape[0], n_chains) + self.assertEqual(evals.shape[1], n_iterations) + self.assertEqual(evals.shape[2], 3) # Test with a loglikelihood mcmc = pints.MCMCController( From 0a077db62755901f8b7ac10ff8952e156d1932fc Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Wed, 18 Dec 2024 23:56:50 +0000 Subject: [PATCH 4/5] Fixed test for 1691 multi-chain. Again. --- pints/tests/test_mcmc_controller.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/pints/tests/test_mcmc_controller.py b/pints/tests/test_mcmc_controller.py index c2446fdc8..5c213a9e8 100755 --- a/pints/tests/test_mcmc_controller.py +++ b/pints/tests/test_mcmc_controller.py @@ -759,12 +759,12 @@ def test_log_pdf_storage_in_memory_multi(self): priors = [self.log_prior(x) for x in chain] self.assertTrue(np.all(evals[i, :, 2] == priors)) - # Test with a sensitivity-using method + # Test with a sensitivity-using multi-chain method # We don't have any of these! mcmc = pints.MCMCController( self.log_posterior, n_chains, xs, method=FakeMultiChainSamplerWithSensitivities) - mcmc.set_max_iterations(5) + mcmc.set_max_iterations(n_iterations) mcmc.set_log_to_screen(False) mcmc.set_log_pdf_storage(True) chains = mcmc.run() @@ -773,6 +773,13 @@ def test_log_pdf_storage_in_memory_multi(self): self.assertEqual(evals.shape[0], n_chains) self.assertEqual(evals.shape[1], n_iterations) self.assertEqual(evals.shape[2], 3) + for i, chain in enumerate(chains): + posteriors = [self.log_posterior(x) for x in chain] + self.assertTrue(np.all(evals[i, :, 0] == posteriors)) + likelihoods = [self.log_likelihood(x) for x in chain] + self.assertTrue(np.all(evals[i, :, 1] == likelihoods)) + priors = [self.log_prior(x) for x in chain] + self.assertTrue(np.all(evals[i, :, 2] == priors)) # Test with a loglikelihood mcmc = pints.MCMCController( @@ -1707,6 +1714,9 @@ def tell(self, fxs): self._fxs = fxs return self._xs, self._fxs, [True] * self._n_chains + def needs_sensitivities(self): + return True + class TestMCMCControllerSingleChainStorage(unittest.TestCase): """ From 8db1e239116c91434c0758ca1f22e72ee388c7c5 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Wed, 18 Dec 2024 23:57:04 +0000 Subject: [PATCH 5/5] Fixed #1691. --- pints/_mcmc/__init__.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/pints/_mcmc/__init__.py b/pints/_mcmc/__init__.py index cb8d5cf3e..5c137d5bc 100644 --- a/pints/_mcmc/__init__.py +++ b/pints/_mcmc/__init__.py @@ -756,7 +756,10 @@ def run(self): if store_evaluations: # If accepted, update log_pdf and prior for logging if accepted: - current_logpdf[i] = fy + if self._needs_sensitivities: + current_logpdf[i] = fy[0] + else: + current_logpdf[i] = fy if prior is not None: current_prior[i] = prior(y) @@ -820,7 +823,10 @@ def run(self): # Check if accepted, if so, update log_pdf and # prior to be logged if accepted[i]: - current_logpdf[i] = fys[i] + if self._needs_sensitivities: + current_logpdf[i] = fys[i][0] + else: + current_logpdf[i] = fys[i] if prior is not None: current_prior[i] = prior(ys[i])