From aca9102cd0a02b9ec48747b34c51c0608ed7f538 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Thu, 20 May 2021 23:13:58 +0200 Subject: [PATCH 01/18] mcmc: tempereture: basic implementation [skip travis] --- cobaya/collection.py | 32 +++++++++++++++++++++++++++----- cobaya/conventions.py | 4 ++-- cobaya/samplers/mcmc/mcmc.py | 27 +++++++++++++++++++-------- cobaya/samplers/mcmc/mcmc.yaml | 2 ++ 4 files changed, 50 insertions(+), 15 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index ad7cef384..11cc036c2 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -102,7 +102,7 @@ class SampleCollection(BaseCollection): def __init__(self, model, output=None, cache_size=_default_cache_size, name=None, extension=None, file_name=None, resuming=False, load=False, - onload_skip=0, onload_thin=1): + weight_temperature=None, onload_skip=0, onload_thin=1): super().__init__(model, name) self.cache_size = cache_size # Create/load the main data frame and the tracking indices @@ -167,6 +167,7 @@ def __init__(self, model, output=None, cache_size=_default_cache_size, name=None self._out_delete() if not resuming and not load: self.reset() + self.weight_temperature = weight_temperature # Prepare fast numpy cache self._icol = {col: i for i, col in enumerate(self.columns)} self._cache_reset() @@ -390,7 +391,8 @@ def copy(self) -> 'SampleCollection': return self._copy() # Statistical computations - def mean(self, first=None, last=None, derived=False, pweight=False): + def mean(self, first=None, last=None, derived=False, pweight=False, + ignore_temperature=False): """ Returns the (weighted) mean of the parameters in the chain, between `first` (default 0) and `last` (default last obtained), @@ -404,12 +406,16 @@ def mean(self, first=None, last=None, derived=False, pweight=False): logps -= max(logps) weights = np.exp(logps) else: - weights = self[OutPar.weight][first:last].to_numpy() + if ignore_temperature: + weights = self[OutPar.weight][first:last].to_numpy() + else: + weights = self.tempered_weights(first, last) return np.average(self[list(self.sampled_params) + (list(self.derived_params) if derived else [])] [first:last].T, weights=weights, axis=-1) - def cov(self, first=None, last=None, derived=False, pweight=False): + def cov(self, first=None, last=None, derived=False, pweight=False, + ignore_temperature=False): """ Returns the (weighted) covariance matrix of the parameters in the chain, between `first` (default 0) and `last` (default last obtained), @@ -424,7 +430,10 @@ def cov(self, first=None, last=None, derived=False, pweight=False): weights = np.exp(logps) kwarg = "aweights" else: - weights = self[OutPar.weight][first:last].to_numpy() + if ignore_temperature: + weights = self[OutPar.weight][first:last].to_numpy() + else: + weights = self.tempered_weights(first, last) kwarg = "fweights" if np.allclose(np.round(weights), weights) else "aweights" weights_kwarg = {kwarg: weights} return np.atleast_2d(np.cov( @@ -432,6 +441,19 @@ def cov(self, first=None, last=None, derived=False, pweight=False): (list(self.derived_params) if derived else [])][first:last].T, **weights_kwarg)) + def tempered_weights(self, first=None, last=None): + """ + Returns the weights after applying the sampling temperature. + """ + if self.weight_temperature: + maxlogpost = -self.MAP()[OutPar.minuslogpost] + reweighting = np.exp( + (-self[OutPar.minuslogpost][first:last].to_numpy() - maxlogpost) * + (1 - 1 / self.weight_temperature)) + else: + reweighting = 1 + return self[OutPar.weight][first:last].to_numpy() * reweighting + def filtered_copy(self, where) -> 'SampleCollection': return self._copy(self.data[where].reset_index(drop=True)) diff --git a/cobaya/conventions.py b/cobaya/conventions.py index 78ae92840..90da9c379 100644 --- a/cobaya/conventions.py +++ b/cobaya/conventions.py @@ -40,9 +40,9 @@ def get_version(): # Names for the samples' fields internally and in the output class OutPar: weight = "weight" # sample weight - # log-posterior, or in general the total log-probability + # minus log-posterior, or in general the total minus log-probability minuslogpost = "minuslogpost" - minuslogprior = "minuslogprior" # log-prior + minuslogprior = "minuslogprior" # minus log-prior chi2 = "chi2" # chi^2 = -2 * loglike (not always normalized to be useful) diff --git a/cobaya/samplers/mcmc/mcmc.py b/cobaya/samplers/mcmc/mcmc.py index 1255ad28a..4a74ca0b7 100644 --- a/cobaya/samplers/mcmc/mcmc.py +++ b/cobaya/samplers/mcmc/mcmc.py @@ -47,6 +47,7 @@ class MCMC(CovmatSampler): learn_every: NumberWithUnits output_every: NumberWithUnits callback_every: NumberWithUnits + temperature: float max_tries: NumberWithUnits max_samples: int drag: bool @@ -110,7 +111,8 @@ def initialize(self): # One collection per MPI process: `name` is the MPI rank + 1 name = str(1 + mpi.rank()) self.collection = SampleCollection( - self.model, self.output, name=name, resuming=self.output.is_resuming()) + self.model, self.output, name=name, resuming=self.output.is_resuming(), + weight_temperature=self.temperature) self.current_point = OneSamplePoint(self.model) # Use standard MH steps by default self.get_new_sample = self.get_new_sample_metropolis @@ -260,6 +262,10 @@ def set_proposer_blocking(self): "Use 'oversample_power' to control the amount of " "dragging steps.") # END OF DEPRECATION BLOCK + if self.temperature: + raise LoggedError( + self.log, + "Temperature != 1 and dragging are not compatible at the moment.") self.get_new_sample = self.get_new_sample_dragging self.mpi_info("Dragging with number of interpolating steps:") max_width = len(str(self.drag_interp_steps)) @@ -526,7 +532,10 @@ def metropolis_accept(self, logp_trial, logp_current): elif logp_trial > logp_current: return True else: - return self._rng.standard_exponential() > (logp_current - logp_trial) + posterior_ratio = logp_current - logp_trial + if self.temperature: + posterior_ratio /= self.temperature + return self._rng.standard_exponential() > posterior_ratio def process_accept_or_reject(self, accept_state, trial, trial_results): """Processes the acceptance/rejection of the new point.""" @@ -597,8 +606,8 @@ def check_convergence_and_learn_proposal(self): if more_than_one_process(): # Compute and gather means and covs use_first = int(self.n() / 2) - mean = self.collection.mean(first=use_first) - cov = self.collection.cov(first=use_first) + mean = self.collection.mean(first=use_first, ignore_temperature=True) + cov = self.collection.cov(first=use_first, ignore_temperature=True) acceptance_rate = self.get_acceptance_rate(use_first) Ns, means, covs, acceptance_rates = mpi.array_gather( [self.n(), mean, cov, acceptance_rate]) @@ -610,11 +619,13 @@ def check_convergence_and_learn_proposal(self): acceptance_rate = self.get_acceptance_rate(cut) Ns = np.ones(m - 1) * cut means = np.array( - [self.collection.mean(first=i * cut, last=(i + 1) * cut - 1) for i in - range(1, m)]) + [self.collection.mean(first=i * cut, last=(i + 1) * cut - 1, + ignore_temperature=True) + for i in range(1, m)]) covs = np.array( - [self.collection.cov(first=i * cut, last=(i + 1) * cut - 1) for i in - range(1, m)]) + [self.collection.cov(first=i * cut, last=(i + 1) * cut - 1, + ignore_temperature=True) + for i in range(1, m)]) except: self.log.info("Not enough points in chain to check convergence. " "Waiting for next checkpoint.") diff --git a/cobaya/samplers/mcmc/mcmc.yaml b/cobaya/samplers/mcmc/mcmc.yaml index c3f866906..bfbe4747e 100644 --- a/cobaya/samplers/mcmc/mcmc.yaml +++ b/cobaya/samplers/mcmc/mcmc.yaml @@ -20,6 +20,8 @@ proposal_scale: 2.4 output_every: 60s # Number of distinct accepted points between proposal learn & convergence checks learn_every: 40d +# Posterior temperature: >1 for more exploratory chains +temperature: # Proposal covariance matrix learning # ----------------------------------- learn_proposal: True From 937c19d5d740a4774fbf04d483fc88d8751b582c Mon Sep 17 00:00:00 2001 From: torradocacho Date: Tue, 7 Sep 2021 21:41:15 +0200 Subject: [PATCH 02/18] more on temperatures [skip travis] --- cobaya/collection.py | 92 ++++++++++++++++++++++++++-------- cobaya/model.py | 38 ++++++++++---- cobaya/samplers/mcmc/mcmc.py | 35 ++++++------- cobaya/samplers/mcmc/mcmc.yaml | 2 +- 4 files changed, 117 insertions(+), 50 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index 5b12a8d3e..642fba5f6 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -94,6 +94,10 @@ class SampleCollection(BaseCollection): but slicing can be done on the ``SampleCollection`` itself (returns a copy, not a view). + When ``temperature`` is different from 1, weights and log-posterior (but not priors + or likelihoods' chi squared's) are those of the tempered sample, obtained assuming a + posterior raised to the power of ``1/temperature``. + Note for developers: when expanding this class or inheriting from it, always access the underlying DataFrame as `self.data` and not `self._data`, to ensure the cache has been dumped. If you really need to access the actual attribute `self._data` in a @@ -102,7 +106,7 @@ class SampleCollection(BaseCollection): def __init__(self, model, output=None, cache_size=_default_cache_size, name=None, extension=None, file_name=None, resuming=False, load=False, - weight_temperature=None, onload_skip=0, onload_thin=1): + temperature=1, onload_skip=0, onload_thin=1): super().__init__(model, name) self.cache_size = cache_size # Create/load the main data frame and the tracking indices @@ -167,7 +171,7 @@ def __init__(self, model, output=None, cache_size=_default_cache_size, name=None self._out_delete() if not resuming and not load: self.reset() - self.weight_temperature = weight_temperature + self.temperature = temperature # Prepare fast numpy cache self._icol = {col: i for i, col in enumerate(self.columns)} self._cache_reset() @@ -184,6 +188,7 @@ def add(self, values: Union[Sequence[float], np.ndarray], logpriors: Optional[Sequence[float]] = None, loglikes: Optional[Sequence[float]] = None, derived: Optional[Sequence[float]] = None, + temperature: Optional[float] = None, weight: float = 1): """ Adds a point to the collection. @@ -192,8 +197,8 @@ def add(self, values: Union[Sequence[float], np.ndarray], `logpriors`, `loglikes` are both required). """ logposterior = self._check_before_adding( - values, logpost=logpost, logpriors=logpriors, - loglikes=loglikes, derived=derived, weight=weight) + values, logpost=logpost, logpriors=logpriors, loglikes=loglikes, + derived=derived, temperature=temperature, weight=weight) self._cache_add(values, logposterior=logposterior, weight=weight) def _check_before_adding(self, values: Union[Sequence[float], np.ndarray], @@ -201,6 +206,7 @@ def _check_before_adding(self, values: Union[Sequence[float], np.ndarray], logpriors: Optional[Sequence[float]] = None, loglikes: Optional[Sequence[float]] = None, derived: Optional[Sequence[float]] = None, + temperature: Optional[float] = None, weight: float = 1 ) -> LogPosterior: """ @@ -239,11 +245,26 @@ def _check_before_adding(self, values: Union[Sequence[float], np.ndarray], self.log, "derived params not consistent with those of LogPosterior object " "passed.") + if temperature is not None: + if not np.isclose(temperature, logpost.temperature): + raise LoggedError( + self.log, + "temperature not consistent with that of LogPosterior object " + "passed.") + if not np.isclose(self.temperature, logpost.temperature): + raise LoggedError( + self.log, + "temperature of added point not consistent with that of sample.") return_logpost = logpost elif isinstance(logpost, float) or logpost is None: + if temperature is not None: + if not np.isclose(temperature, self.temperature): + raise LoggedError( + self.log, "temperature not consistent with that of sample.") try: return_logpost = LogPosterior( - logpriors=logpriors, loglikes=loglikes, derived=derived) + logpost=logpost, logpriors=logpriors, loglikes=loglikes, + derived=derived, temperature=temperature) except ValueError as valerr: # missing logpriors/loglikes if logpost is None, # or inconsistent sum if logpost given @@ -457,18 +478,46 @@ def cov(self, first=None, last=None, derived=False, pweight=False, dtype=np.float64).T, **weights_kwarg)) - def tempered_weights(self, first=None, last=None): + def reweight(self, importance_weights): + self._cache_dump() + self._data[OutPar.weight] *= importance_weights + self._data = self.data[self._data.weight > 0].reset_index(drop=True) + self._n = self._data.last_valid_index() + 1 + + def detempering_reweight_factor(self, first=None, last=None): """ - Returns the weights after applying the sampling temperature. + Returns the reweighting factor necessary to remove the effect of the sampling + temperature. """ - if self.weight_temperature: + if self.temperature != 1: + minuslogpost = self[OutPar.minuslogpost][first:last].to_numpy() maxlogpost = -self.MAP()[OutPar.minuslogpost] - reweighting = np.exp( - (-self[OutPar.minuslogpost][first:last].to_numpy() - maxlogpost) * - (1 - 1 / self.weight_temperature)) + return np.exp((-minuslogpost - maxlogpost) * (1 - 1 / self.temperature)) + else: + return 1 + + def detempered_minuslogpost(self, first=None, last=None): + """ + Returns the minus log-posterior of the original (temperature 1) pdf at the + samples, as a numpy array. + """ + if self.temperature != 1: + return (self.data[OutPar.minuslogprior][first:last].to_numpy() + + self.data[OutPar.chi2][first:last].to_numpy() / 2) else: - reweighting = 1 - return self[OutPar.weight][first:last].to_numpy() * reweighting + return self[OutPar.minuslogpost][first:last].to_numpy() + + def detemper(self): + """ + Removes the effect of sampling temperature. + + This cannot be fully undone (e.g. recovering the original integer weights). + You may want to call this method on a copy (see :func:`~SampleCollection.copy`). + """ + self._cache_dump() + self._data[OutPar.minuslogpost] = self.detempered_minuslogpost() + self.reweight(self.detempering_reweight_factor()) + self.temperature = 1 def filtered_copy(self, where) -> 'SampleCollection': return self._copy(self.data[where].reset_index(drop=True)) @@ -528,12 +577,6 @@ def sampled_to_getdist_mcsamples(self, first=None, last=None): names=names) return mcsamples - def reweight(self, importance_weights): - self._cache_dump() - self._data[OutPar.weight] *= importance_weights - self._data = self.data[self._data.weight > 0].reset_index(drop=True) - self._n = self._data.last_valid_index() + 1 - # Saving and updating def _get_driver(self, method): return getattr(self, method + derived_par_name_separator + self.driver) @@ -647,8 +690,17 @@ def add(self, values, logpost: LogPosterior): def logpost(self): return self.results.logpost + @property + def temperature(self): + return self.results.temperature + def add_to_collection(self, collection: SampleCollection): - """Adds this point at the end of a given collection.""" + """ + Adds this point at the end of a given collection. + + It is assumed that both this instance and the collection passed were + initialised with the same :class:`model.Model` (no checks are performed). + """ if self.output_thin > 1: self._added_weight += self.weight if self._added_weight >= self.output_thin: diff --git a/cobaya/model.py b/cobaya/model.py index 679832956..fcb2132cc 100644 --- a/cobaya/model.py +++ b/cobaya/model.py @@ -58,12 +58,16 @@ class LogPosterior: If ``finite=True`` (default: False), it will try to represent infinities as the largest real numbers allowed by machine precision. + + If ``temperature`` is different from 1, the log-posterior (but not priors + or likelihoods' chi squared's) is stored to the power of ``1/temperature``. """ logpost: Optional[float] = None logpriors: Optional[Sequence[float]] = None loglikes: Optional[Sequence[float]] = None derived: Optional[Sequence[float]] = None + temperature: float = 1 finite: Optional[bool] = False logprior: float = dataclasses.field(init=False, repr=False) loglike: float = dataclasses.field(init=False, repr=False) @@ -100,9 +104,11 @@ def _logpost_is_consistent(self): """ if self.finite: return np.isclose(np.nan_to_num(self.logpost), - np.nan_to_num(self.logprior + self.loglike)) + np.nan_to_num( + (self.logprior + self.loglike) / self.temperature)) else: - return np.isclose(self.logpost, self.logprior + self.loglike) + return np.isclose(self.logpost, + (self.logprior + self.loglike) / self.temperature) def make_finite(self): """ @@ -483,8 +489,8 @@ def loglike(self, def logposterior(self, params_values: Union[Dict[str, float], Sequence[float]], as_dict: bool = False, make_finite: bool = False, - return_derived: bool = True, cached: bool = True, - _no_check: bool = False + return_derived: bool = True, temperature: float = 1, + cached: bool = True, _no_check: bool = False ) -> Union[LogPosterior, dict]: """ Takes an array or dictionary of sampled parameter values. @@ -517,6 +523,9 @@ def logposterior(self, If ``make_finite=True``, it will try to represent infinities as the largest real numbers allowed by machine precision. + If ``temperature`` is different from 1, the log-posterior (but not priors + or likelihoods' chi squared's) is returned to the power of ``1/temperature``. + If ``cached=False`` (default: True), it ignores previously computed results that could be reused. """ @@ -554,8 +563,9 @@ def logposterior(self, else: loglikes = [] derived_sampler = [] - logposterior = LogPosterior(logpriors=logpriors, loglikes=loglikes, - derived=derived_sampler, finite=make_finite) + logposterior = LogPosterior( + logpriors=logpriors, loglikes=loglikes, temperature=temperature, + derived=derived_sampler, finite=make_finite) if as_dict: return logposterior.as_dict(self) else: @@ -563,7 +573,8 @@ def logposterior(self, def logpost(self, params_values: Union[Dict[str, float], Sequence[float]], - make_finite: bool = False, cached: bool = True) -> float: + make_finite: bool = False, temperature: float = 1, cached: bool = True + ) -> float: """ Takes an array or dictionary of sampled parameter values. If the argument is an array, parameters must have the same order as in the input. @@ -575,11 +586,15 @@ def logpost(self, If ``make_finite=True``, it will try to represent infinities as the largest real numbers allowed by machine precision. + If ``temperature`` is different from 1, the log-posterior is returned to the + power of ``1/temperature``. + If ``cached=False`` (default: True), it ignores previously computed results that could be reused. """ - return self.logposterior(params_values, make_finite=make_finite, - return_derived=False, cached=cached).logpost + return self.logposterior( + params_values, make_finite=make_finite, return_derived=False, + temperature=temperature, cached=cached).logpost def get_valid_point(self, max_tries: int, ignore_fixed_ref: bool = False, logposterior_as_dict: bool = False, random_state=None, @@ -596,6 +611,9 @@ def get_valid_point(self, max_tries: int, ignore_fixed_ref: bool = False, Returns (point, LogPosterior(logpost, logpriors, loglikes, derived)) + If ``temperature`` is different from 1, the log-posterior (but not priors + or likelihoods' chi squared's) is returned to the power of ``1/temperature``. + If ``logposterior_as_dict=True`` (default: False), returns for the log-posterior a dictionary containing prior names, likelihoods names and, if applicable, derived parameters names as keys, instead of a :class:`~model.LogPosterior` object. @@ -605,7 +623,7 @@ def get_valid_point(self, max_tries: int, ignore_fixed_ref: bool = False, ignore_fixed=ignore_fixed_ref, warn_if_no_ref=not loop, random_state=random_state) - results = self.logposterior(initial_point) + results = self.logposterior(initial_point, temperature=temperature) if results.logpost != -np.inf: break else: diff --git a/cobaya/samplers/mcmc/mcmc.py b/cobaya/samplers/mcmc/mcmc.py index 26e09373a..c46dd309e 100644 --- a/cobaya/samplers/mcmc/mcmc.py +++ b/cobaya/samplers/mcmc/mcmc.py @@ -110,7 +110,7 @@ def initialize(self): name = str(1 + mpi.rank()) self.collection = SampleCollection( self.model, self.output, name=name, resuming=self.output.is_resuming(), - weight_temperature=self.temperature) + temperature=self.temperature) self.current_point = OneSamplePoint(self.model) # Use standard MH steps by default self.get_new_sample = self.get_new_sample_metropolis @@ -147,14 +147,16 @@ def initialize(self): loglikes=-0.5 * (self.collection[self.collection.chi2_names] .iloc[last].to_numpy(dtype=np.float64, copy=True)), derived=(self.collection[self.collection.derived_params].iloc[last] - .to_numpy(dtype=np.float64, copy=True))) + .to_numpy(dtype=np.float64, copy=True)), + temperature=self.temperature) else: # NB: max_tries adjusted to dim instead of #cycles (blocking not computed yet) self.max_tries.set_scale(self.model.prior.d()) self.log.info("Getting initial point... (this may take a few seconds)") initial_point, results = \ - self.model.get_valid_point(max_tries=self.max_tries.value, - random_state=self._rng) + self.model.get_valid_point( + max_tries=self.max_tries.value, random_state=self._rng, + temperature=self.temperature) # If resuming but no existing chain, assume failed run and ignore blocking # if speeds measurement requested if self.output.is_resuming() and not len(self.collection) \ @@ -258,10 +260,6 @@ def set_proposer_blocking(self): "Use 'oversample_power' to control the amount" " of dragging steps.") # END OF DEPRECATION BLOCK - if self.temperature: - raise LoggedError( - self.log, - "Temperature != 1 and dragging are not compatible at the moment.") self.get_new_sample = self.get_new_sample_dragging self.mpi_info("Dragging with number of interpolating steps:") max_width = len(str(self.drag_interp_steps)) @@ -418,7 +416,7 @@ def get_new_sample_metropolis(self): """ trial = self.current_point.values.copy() self.proposer.get_proposal(trial) - trial_results = self.model.logposterior(trial) + trial_results = self.model.logposterior(trial, temperature=self.temperature) accept = self.metropolis_accept(trial_results.logpost, self.current_point.logpost) self.process_accept_or_reject(accept, trial, trial_results) return accept @@ -444,7 +442,8 @@ def get_new_sample_dragging(self): self.log.debug("Proposed slow end-point: %r", current_end_point) # Save derived parameters of delta_slow jump, in case I reject all the dragging # steps but accept the move in the slow direction only - current_end = self.model.logposterior(current_end_point) + current_end = self.model.logposterior( + current_end_point, temperature=self.temperature) if current_end.logpost == -np.inf: self.current_point.weight += 1 return False @@ -471,13 +470,13 @@ def get_new_sample_dragging(self): # point, but discard them, since they contain the starting point's fast ones, # not used later -- save the end point's ones. proposal_start_logpost = self.model.logposterior( - proposal_start_point, return_derived=derived, _no_check=True).logpost - + proposal_start_point, return_derived=derived, + temperature=self.temperature, _no_check=True).logpost if proposal_start_logpost != -np.inf: proposal_end_point = current_end_point + delta_fast proposal_end = self.model.logposterior( - proposal_end_point, return_derived=derived, _no_check=True) - + proposal_end_point, return_derived=derived, + temperature=self.temperature, _no_check=True) if proposal_end.logpost != -np.inf: # create the interpolated probability and do a Metropolis test @@ -510,8 +509,8 @@ def get_new_sample_dragging(self): start_drag_logpost_acc / n_average) if accept and not derived: # recompute with derived parameters (slow parameter ones should be cached) - current_end = self.model.logposterior(current_end_point) - + current_end = self.model.logposterior( + current_end_point, temperature=self.temperature) self.process_accept_or_reject(accept, current_end_point, current_end) self.log.debug("TOTAL step: %s", ("accepted" if accept else "rejected")) return accept @@ -529,8 +528,6 @@ def metropolis_accept(self, logp_trial, logp_current): return True else: posterior_ratio = logp_current - logp_trial - if self.temperature: - posterior_ratio /= self.temperature return self._rng.standard_exponential() > posterior_ratio def process_accept_or_reject(self, accept_state, trial, trial_results): @@ -557,7 +554,7 @@ def process_accept_or_reject(self, accept_state, trial, trial_results): self.current_point.weight += 1 # Failure criterion: chain stuck! (but be more permissive during burn_in) max_tries_now = self.max_tries.value * ( - 1 + (10 - 1) * np.sign(self.burn_in_left)) + 1 + (10 - 1) * np.sign(self.burn_in_left)) if self.current_point.weight > max_tries_now: self.collection.out_update() raise LoggedError( diff --git a/cobaya/samplers/mcmc/mcmc.yaml b/cobaya/samplers/mcmc/mcmc.yaml index 2dcccef34..b2c669a12 100644 --- a/cobaya/samplers/mcmc/mcmc.yaml +++ b/cobaya/samplers/mcmc/mcmc.yaml @@ -21,7 +21,7 @@ output_every: 60s # Number of distinct accepted points between proposal learn & convergence checks learn_every: 40d # Posterior temperature: >1 for more exploratory chains -temperature: +temperature: 1 # Proposal covariance matrix learning # ----------------------------------- learn_proposal: True From 1606031d7e3f0e81d370f6e1307d26f38b8ea1c9 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Wed, 8 Sep 2021 10:39:21 +0200 Subject: [PATCH 03/18] collection: some type hinting --- cobaya/collection.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index 642fba5f6..5d3ef644b 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -424,8 +424,8 @@ def copy(self) -> 'SampleCollection': """ return self._copy() - def mean(self, first=None, last=None, derived=False, pweight=False, - ignore_temperature=False): + def mean(self, first: Optional[int] = None, last: Optional[int] = None, + derived=False, pweight=False, ignore_temperature=False): """ Returns the (weighted) mean of the parameters in the chain, between `first` (default 0) and `last` (default last obtained), @@ -449,8 +449,8 @@ def mean(self, first=None, last=None, derived=False, pweight=False, [first:last].to_numpy(dtype=np.float64).T, weights=weights, axis=-1) - def cov(self, first=None, last=None, derived=False, pweight=False, - ignore_temperature=False): + def cov(self, first: Optional[int] = None, last: Optional[int] = None, + derived=False, pweight=False, ignore_temperature=False): """ Returns the (weighted) covariance matrix of the parameters in the chain, between `first` (default 0) and `last` (default last obtained), @@ -484,7 +484,8 @@ def reweight(self, importance_weights): self._data = self.data[self._data.weight > 0].reset_index(drop=True) self._n = self._data.last_valid_index() + 1 - def detempering_reweight_factor(self, first=None, last=None): + def detempering_reweight_factor(self, first: Optional[int] = None, + last: Optional[int] = None): """ Returns the reweighting factor necessary to remove the effect of the sampling temperature. @@ -496,7 +497,8 @@ def detempering_reweight_factor(self, first=None, last=None): else: return 1 - def detempered_minuslogpost(self, first=None, last=None): + def detempered_minuslogpost(self, first: Optional[int] = None, + last: Optional[int] = None): """ Returns the minus log-posterior of the original (temperature 1) pdf at the samples, as a numpy array. @@ -558,7 +560,8 @@ def MAP(self): """Maximum-a-posteriori (MAP) sample. Returns a copy.""" return self.data.loc[self.data[OutPar.minuslogpost].idxmin()].copy() - def sampled_to_getdist_mcsamples(self, first=None, last=None): + def sampled_to_getdist_mcsamples(self, first: Optional[int] = None, + last: Optional[int] = None): """ Basic interface with getdist -- internal use only! (For analysis and plotting use `getdist.mcsamples.MCSamplesFromCobaya From 5a12f05d5cd3770f5ba29158968c85695f8c6fa2 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Wed, 8 Sep 2021 11:29:45 +0200 Subject: [PATCH 04/18] temperature: fixes --- cobaya/collection.py | 2 +- cobaya/model.py | 16 ++++++++++------ 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index 5d3ef644b..b160de0ef 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -171,7 +171,7 @@ def __init__(self, model, output=None, cache_size=_default_cache_size, name=None self._out_delete() if not resuming and not load: self.reset() - self.temperature = temperature + self.temperature = temperature if temperature is not None else 1 # Prepare fast numpy cache self._icol = {col: i for i, col in enumerate(self.columns)} self._cache_reset() diff --git a/cobaya/model.py b/cobaya/model.py index fcb2132cc..09832f58e 100644 --- a/cobaya/model.py +++ b/cobaya/model.py @@ -83,13 +83,15 @@ def __post_init__(self): object.__setattr__( self, 'loglike', sum(self.loglikes) if self.loglikes is not None else None) + if self.temperature is None: + object.__setattr__(self, 'temperature', 1) if self.finite: self.make_finite() if self.logpost is None: if self.logpriors is None or self.loglikes is None: raise ValueError("If `logpost` not passed, both `logpriors` and " "`loglikes` must be passed.") - object.__setattr__(self, 'logpost', self.logprior + self.loglike) + object.__setattr__(self, 'logpost', self._logpost()) elif self.logpriors is not None and self.loglikes is not None: if not self._logpost_is_consistent(): raise ValueError("The given log-posterior is not equal to the " @@ -97,18 +99,19 @@ def __post_init__(self): "%g != sum(%r) + sum(%r)" % (self.logpost, self.logpriors, self.loglikes)) + def _logpost(self): + """Computes logpost from prior and likelihood product.""" + return (self.logprior + self.loglike) / self.temperature + def _logpost_is_consistent(self): """ Checks that the sum of logpriors and loglikes (if present) add up to logpost, if passed. """ if self.finite: - return np.isclose(np.nan_to_num(self.logpost), - np.nan_to_num( - (self.logprior + self.loglike) / self.temperature)) + return np.isclose(np.nan_to_num(self.logpost), np.nan_to_num(self._logpost())) else: - return np.isclose(self.logpost, - (self.logprior + self.loglike) / self.temperature) + return np.isclose(self.logpost, self._logpost()) def make_finite(self): """ @@ -598,6 +601,7 @@ def logpost(self, def get_valid_point(self, max_tries: int, ignore_fixed_ref: bool = False, logposterior_as_dict: bool = False, random_state=None, + temperature: Optional[float] = None, ) -> Union[Tuple[np.ndarray, LogPosterior], Tuple[np.ndarray, dict]]: """ From 4e803628c335e7fae90743dd2fa3c7aa103f6748 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Wed, 8 Sep 2021 17:51:07 +0200 Subject: [PATCH 05/18] collection: bugfixes and reworking [skip travis] --- cobaya/collection.py | 96 +++++++++++++++++++++++++++++--------------- 1 file changed, 63 insertions(+), 33 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index b160de0ef..de65d87ab 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -413,7 +413,7 @@ def _copy(self, data=None) -> 'SampleCollection': delattr(self, "_data") self_copy = deepcopy(self) setattr(self, "_data", current_data) - setattr(self_copy, "_data", data) + setattr(self_copy, "_data", data.copy()) setattr(self_copy, "_n", data.last_valid_index() + 1) return self_copy @@ -424,59 +424,88 @@ def copy(self) -> 'SampleCollection': """ return self._copy() - def mean(self, first: Optional[int] = None, last: Optional[int] = None, - derived=False, pweight=False, ignore_temperature=False): + def _weights_for_stats(self, first: Optional[int] = None, last: Optional[int] = None, + pweight: bool = False, ignore_temperature: bool = False + ) -> np.ndarray: """ - Returns the (weighted) mean of the parameters in the chain, - between `first` (default 0) and `last` (default last obtained), - optionally including derived parameters if `derived=True` (default `False`). + Returns weights for computation of statistical quantities such as mean and + covariance. - If `pweight=True` (default `False`) weights every point with its probability. - The estimate of the mean in this case is unstable; use carefully. + If ``ignore_temperature=True``, (default ``False``), the sample is weighted as + if the tempered logposterior were the real one (e.g. variances would be larger + for hot samples). + + If ``pweight=True`` (default ``False``) every point is weighted with its + probability, and sample weights are ignored (may lead to very inaccurate + estimates!). """ if pweight: - logps = -self[OutPar.minuslogpost][first:last].to_numpy(dtype=np.float64, - copy=True) + if ignore_temperature: + logps = -self[OutPar.minuslogpost][first:last]\ + .to_numpy(dtype=np.float64, copy=True) + else: + logps = -self.detempered_minuslogpost(first, last) logps -= max(logps) weights = np.exp(logps) else: - if ignore_temperature: - weights = self[OutPar.weight][first:last].to_numpy(dtype=np.float64) - else: - weights = self.tempered_weights(first, last) + weights = self[OutPar.weight][first:last]\ + .to_numpy(dtype=np.float64, copy=True) + if not ignore_temperature: + weights *= self.detempering_reweight_factor(first, last) + return weights + + def mean(self, first: Optional[int] = None, last: Optional[int] = None, + pweight: bool = False, ignore_temperature: bool = False, + derived: bool = False) -> np.ndarray: + """ + Returns the (weighted) mean of the parameters in the chain, + between `first` (default 0) and `last` (default last obtained), + optionally including derived parameters if `derived=True` (default `False`). + + If ``ignore_temperature=True``, (default ``False``), the sample is weighted as + if the tempered logposterior were the real one (e.g. variances would be larger + for hot samples). + + If ``pweight=True`` (default ``False``) every point is weighted with its + probability, and sample weights are ignored (may lead to very inaccurate + estimates!). + """ + weights = self._weights_for_stats(first, last, pweight=pweight, + ignore_temperature=ignore_temperature) return np.average(self[list(self.sampled_params) + (list(self.derived_params) if derived else [])] [first:last].to_numpy(dtype=np.float64).T, weights=weights, axis=-1) def cov(self, first: Optional[int] = None, last: Optional[int] = None, - derived=False, pweight=False, ignore_temperature=False): + pweight: bool = False, ignore_temperature: bool = False, + derived: bool = False) -> np.ndarray: """ Returns the (weighted) covariance matrix of the parameters in the chain, between `first` (default 0) and `last` (default last obtained), optionally including derived parameters if `derived=True` (default `False`). - If `pweight=True` (default `False`) weights every point with its probability. - The estimate of the covariance matrix in this case is unstable; use carefully. + If ``ignore_temperature=True``, (default ``False``), the sample is weighted as + if the tempered logposterior were the real one (e.g. variances would be larger + for hot samples). + + If ``pweight=True`` (default ``False``) every point is weighted with its + probability, and sample weights are ignored (may lead to very inaccurate + estimates!). """ + weights = self._weights_for_stats(first, last, pweight=pweight, + ignore_temperature=ignore_temperature) if pweight: - logps = -self[OutPar.minuslogpost][first:last].to_numpy(dtype=np.float64, - copy=True) - logps -= max(logps) - weights = np.exp(logps) - kwarg = "aweights" + weight_type_kwarg = "aweights" + elif np.allclose(np.round(weights), weights): + weight_type_kwarg = "fweights" else: - if ignore_temperature: - weights = self[OutPar.weight][first:last].to_numpy(dtype=np.float64) - else: - weights = self.tempered_weights(first, last) - kwarg = "fweights" if np.allclose(np.round(weights), weights) else "aweights" - weights_kwarg = {kwarg: weights} + weight_type_kwarg = "aweights" return np.atleast_2d(np.cov( self[list(self.sampled_params) + (list(self.derived_params) if derived else [])][first:last].to_numpy( dtype=np.float64).T, - **weights_kwarg)) + **{weight_type_kwarg: weights})) def reweight(self, importance_weights): self._cache_dump() @@ -492,8 +521,8 @@ def detempering_reweight_factor(self, first: Optional[int] = None, """ if self.temperature != 1: minuslogpost = self[OutPar.minuslogpost][first:last].to_numpy() - maxlogpost = -self.MAP()[OutPar.minuslogpost] - return np.exp((-minuslogpost - maxlogpost) * (1 - 1 / self.temperature)) + return np.exp((-minuslogpost - -np.min(minuslogpost)) * + (1 - 1 / self.temperature)) else: return 1 @@ -507,7 +536,7 @@ def detempered_minuslogpost(self, first: Optional[int] = None, return (self.data[OutPar.minuslogprior][first:last].to_numpy() + self.data[OutPar.chi2][first:last].to_numpy() / 2) else: - return self[OutPar.minuslogpost][first:last].to_numpy() + return self[OutPar.minuslogpost][first:last].to_numpy(copy=True) def detemper(self): """ @@ -517,8 +546,9 @@ def detemper(self): You may want to call this method on a copy (see :func:`~SampleCollection.copy`). """ self._cache_dump() - self._data[OutPar.minuslogpost] = self.detempered_minuslogpost() self.reweight(self.detempering_reweight_factor()) + # log-post and temperature must be changed *after* computing reweighting factors! + self._data[OutPar.minuslogpost] = self.detempered_minuslogpost() self.temperature = 1 def filtered_copy(self, where) -> 'SampleCollection': From d524abf37f4649a9ebc259ce33f38d5749d14201 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Wed, 8 Sep 2021 17:53:37 +0200 Subject: [PATCH 06/18] collection: removed old slicing made it unnecessary bc caching --- cobaya/collection.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index de65d87ab..01b10207f 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -336,7 +336,7 @@ def append(self, collection): Append another collection. Internal method: does not check for consistency! """ - self._data = pd.concat([self.data[:len(self)], collection.data], + self._data = pd.concat([self.data, collection.data], ignore_index=True) def __len__(self): @@ -352,11 +352,11 @@ def data(self): # Make the dataframe printable (but only the filled ones!) def __repr__(self): - return self.data[:len(self)].__repr__() + return self.data.__repr__() # Make the dataframe iterable over rows def __iter__(self): - return self.data[:len(self)].iterrows() + return self.data.iterrows() # Accessing the dataframe def __getitem__(self, *args): @@ -601,11 +601,11 @@ def sampled_to_getdist_mcsamples(self, first: Optional[int] = None, # No logging of warnings temporarily, so getdist won't complain unnecessarily with NoLogging(): mcsamples = MCSamples( - samples=self.data[:len(self)][names].to_numpy(dtype=np.float64)[ + samples=self.data[names].to_numpy(dtype=np.float64)[ first:last], - weights=self.data[:len(self)][OutPar.weight].to_numpy(dtype=np.float64)[ + weights=self.data[OutPar.weight].to_numpy(dtype=np.float64)[ first:last], - loglikes=self.data[:len(self)][OutPar.minuslogpost].to_numpy( + loglikes=self.data[OutPar.minuslogpost].to_numpy( dtype=np.float64)[first:last], names=names) return mcsamples From 3484e2784d1485cb7c746b033d7cd369548fed1d Mon Sep 17 00:00:00 2001 From: torradocacho Date: Wed, 8 Sep 2021 18:30:46 +0200 Subject: [PATCH 07/18] mcmc: takes temperature into account for bound computation --- cobaya/collection.py | 7 +++++-- cobaya/samplers/mcmc/mcmc.py | 6 ++++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index 01b10207f..eee135169 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -590,8 +590,9 @@ def MAP(self): """Maximum-a-posteriori (MAP) sample. Returns a copy.""" return self.data.loc[self.data[OutPar.minuslogpost].idxmin()].copy() - def sampled_to_getdist_mcsamples(self, first: Optional[int] = None, - last: Optional[int] = None): + def sampled_to_getdist_mcsamples( + self, first: Optional[int] = None, last: Optional[int] = None, + ignore_temperature: bool = False) -> MCSamples: """ Basic interface with getdist -- internal use only! (For analysis and plotting use `getdist.mcsamples.MCSamplesFromCobaya @@ -608,6 +609,8 @@ def sampled_to_getdist_mcsamples(self, first: Optional[int] = None, loglikes=self.data[OutPar.minuslogpost].to_numpy( dtype=np.float64)[first:last], names=names) + if not ignore_temperature: + mcsamples.cool(self.temperature) return mcsamples # Saving and updating diff --git a/cobaya/samplers/mcmc/mcmc.py b/cobaya/samplers/mcmc/mcmc.py index c46dd309e..16a73470a 100644 --- a/cobaya/samplers/mcmc/mcmc.py +++ b/cobaya/samplers/mcmc/mcmc.py @@ -700,7 +700,8 @@ def check_convergence_and_learn_proposal(self): # in units of the mean standard deviation of the chains if converged_means: if more_than_one_process(): - mcsamples = self.collection.sampled_to_getdist_mcsamples(first=use_first) + mcsamples = self.collection.sampled_to_getdist_mcsamples( + first=use_first, ignore_temperature=True) try: bound = np.array([[ mcsamples.confidence(i, limfrac=self.Rminus1_cl_level / 2., @@ -716,7 +717,8 @@ def check_convergence_and_learn_proposal(self): try: mcsamples_list = [ self.collection.sampled_to_getdist_mcsamples( - first=i * cut, last=(i + 1) * cut - 1) + first=i * cut, last=(i + 1) * cut - 1, + ignore_temperature=True) for i in range(1, m)] except always_stop_exceptions: raise From 0580398910d742661fec5ff018d8296b3f84fb90 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Wed, 22 Sep 2021 15:19:51 +0200 Subject: [PATCH 08/18] mcmc_tempered: bugfix in Collection, and some refactoring --- cobaya/collection.py | 20 +++++++++++--------- cobaya/samplers/mcmc/mcmc.py | 27 +++++++++++++++++++++------ 2 files changed, 32 insertions(+), 15 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index eee135169..df2da727d 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -99,9 +99,9 @@ class SampleCollection(BaseCollection): posterior raised to the power of ``1/temperature``. Note for developers: when expanding this class or inheriting from it, always access - the underlying DataFrame as `self.data` and not `self._data`, to ensure the cache has - been dumped. If you really need to access the actual attribute `self._data` in a - method, make sure to decorate it with `@ensure_cache_dumped`. + the underlying DataFrame as ``self.data`` and not ``self._data``, to ensure the cache + has been dumped. If you really need to access the actual attribute ``self._data`` in a + method, make sure to decorate it with ``@ensure_cache_dumped``. """ def __init__(self, model, output=None, cache_size=_default_cache_size, name=None, @@ -520,9 +520,9 @@ def detempering_reweight_factor(self, first: Optional[int] = None, temperature. """ if self.temperature != 1: - minuslogpost = self[OutPar.minuslogpost][first:last].to_numpy() - return np.exp((-minuslogpost - -np.min(minuslogpost)) * - (1 - 1 / self.temperature)) + logp_temp = -self[OutPar.minuslogpost][first:last].to_numpy(dtype=np.float64) + log_ratio = logp_temp * (self.temperature - 1) + return np.exp(log_ratio - np.max(log_ratio)) else: return 1 @@ -533,10 +533,12 @@ def detempered_minuslogpost(self, first: Optional[int] = None, samples, as a numpy array. """ if self.temperature != 1: - return (self.data[OutPar.minuslogprior][first:last].to_numpy() + - self.data[OutPar.chi2][first:last].to_numpy() / 2) + return ( + self.data[OutPar.minuslogprior][first:last].to_numpy(dtype=np.float64) + + self.data[OutPar.chi2][first:last].to_numpy(dtype=np.float64) / 2) else: - return self[OutPar.minuslogpost][first:last].to_numpy(copy=True) + return self[OutPar.minuslogpost][first:last].to_numpy( + dtype=np.float64, copy=True) def detemper(self): """ diff --git a/cobaya/samplers/mcmc/mcmc.py b/cobaya/samplers/mcmc/mcmc.py index 16a73470a..833f40ac5 100644 --- a/cobaya/samplers/mcmc/mcmc.py +++ b/cobaya/samplers/mcmc/mcmc.py @@ -98,6 +98,8 @@ def initialize(self): self._quants_d_units.append(number) setattr(self, q, number) self.output_every = NumberWithUnits(self.output_every, "s", dtype=int) + if self.temperature is None: + self.temperature = 1 if is_main_process(): if self.output.is_resuming() and ( max(self.mpi_size or 0, 1) != mpi.size()): @@ -169,7 +171,7 @@ def initialize(self): n = None if self.measure_speeds is True else int(self.measure_speeds) self.model.measure_and_set_speeds(n=n, discard=0, random_state=self._rng) self.set_proposer_blocking() - self.set_proposer_covmat(load=True) + self.set_proposer_initial_covmat(load=True) self.current_point.add(initial_point, results) self.log.info("Initial point: %s", self.current_point) @@ -296,10 +298,10 @@ def set_proposer_blocking(self): for number in self._quants_d_units: number.set_scale(self.cycle_length // self.current_point.output_thin) - def set_proposer_covmat(self, load=False): + def set_proposer_initial_covmat(self, load=False): if load: # Build the initial covariance matrix of the proposal, or load from checkpoint - self._covmat, where_nan = self._load_covmat( + self._intial_covmat, where_nan = self._load_covmat( prefer_load_old=self.output.is_resuming()) if np.any(where_nan) and self.learn_proposal: # We want to start learning the covmat earlier. @@ -312,11 +314,11 @@ def set_proposer_covmat(self, load=False): self.learn_proposal_Rminus1_max = self.learn_proposal_Rminus1_max_early self.log.debug( "Sampling with covmat:\n%s", - DataFrame(self._covmat, + DataFrame(self._intial_covmat, columns=self.model.parameterization.sampled_params(), index=self.model.parameterization.sampled_params()).to_string( line_width=line_width)) - self.proposer.set_covariance(self._covmat) + self.proposer.set_covariance(self._temper_covmat(self._intial_covmat)) def _get_last_nondragging_block(self, blocks, speeds): # blocks and speeds are already sorted @@ -772,6 +774,7 @@ def check_convergence_and_learn_proposal(self): return mean_of_covs = mpi.share(mean_of_covs) try: + # No need to temper this covmat: already computed from tempered post self.proposer.set_covariance(mean_of_covs) self.mpi_info(" - Updated covariance matrix of proposal pdf.") self.mpi_debug("%r", mean_of_covs) @@ -781,6 +784,18 @@ def check_convergence_and_learn_proposal(self): # Save checkpoint info self.write_checkpoint() + def _temper_covmat(self, covmat): + """ + Convert covmat to that of ``probability^(1/temperature)`` posterior. + """ + return covmat * self.temperature + + def _detemper_covmat(self, covmat): + """ + Convert covmat of ``probability^(1/temperature)`` posterior to original one. + """ + return covmat / self.temperature + def do_output(self, date_time): self.collection.out_update() msg = "Progress @ %s : " % date_time.strftime("%Y-%m-%d %H:%M:%S") @@ -794,7 +809,7 @@ def do_output(self, date_time): def write_checkpoint(self): if is_main_process() and self.output: checkpoint_filename = self.checkpoint_filename() - self.dump_covmat(self.proposer.get_covariance()) + self.dump_covmat(self._detemper_covmat(self.proposer.get_covariance())) checkpoint_info = {"sampler": {self.get_name(): dict([ ("converged", self.converged), ("Rminus1_last", self.Rminus1_last), From 0f99386b942da41848aa5ba20a2decc96120ccd9 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Wed, 22 Sep 2021 17:47:33 +0200 Subject: [PATCH 09/18] collection: bugfix (new pandas behaviour) and temp in MAP --- cobaya/collection.py | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index df2da727d..a729c22ad 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -431,7 +431,7 @@ def _weights_for_stats(self, first: Optional[int] = None, last: Optional[int] = Returns weights for computation of statistical quantities such as mean and covariance. - If ``ignore_temperature=True``, (default ``False``), the sample is weighted as + If ``ignore_temperature=True`` (default ``False``), the sample is weighted as if the tempered logposterior were the real one (e.g. variances would be larger for hot samples). @@ -462,7 +462,7 @@ def mean(self, first: Optional[int] = None, last: Optional[int] = None, between `first` (default 0) and `last` (default last obtained), optionally including derived parameters if `derived=True` (default `False`). - If ``ignore_temperature=True``, (default ``False``), the sample is weighted as + If ``ignore_temperature=True`` (default ``False``), the sample is weighted as if the tempered logposterior were the real one (e.g. variances would be larger for hot samples). @@ -485,7 +485,7 @@ def cov(self, first: Optional[int] = None, last: Optional[int] = None, between `first` (default 0) and `last` (default last obtained), optionally including derived parameters if `derived=True` (default `False`). - If ``ignore_temperature=True``, (default ``False``), the sample is weighted as + If ``ignore_temperature=True`` (default ``False``), the sample is weighted as if the tempered logposterior were the real one (e.g. variances would be larger for hot samples). @@ -586,11 +586,21 @@ def thin_samples(self, thin, inplace=False) -> 'SampleCollection': def bestfit(self): """Best fit (maximum likelihood) sample. Returns a copy.""" - return self.data.loc[self.data[OutPar.chi2].idxmin()].copy() + return self.data.loc[self.data[OutPar.chi2].astype(np.float64).idxmin()].copy() - def MAP(self): - """Maximum-a-posteriori (MAP) sample. Returns a copy.""" - return self.data.loc[self.data[OutPar.minuslogpost].idxmin()].copy() + def MAP(self, ignore_temperature: bool = False): + """ + Maximum-a-posteriori (MAP) sample. Returns a copy. + + If ``ignore_temperature=True`` (default ``False``), it returns the maximum + tempered posterior, instead of the original one. + """ + i_max_logpost = self.data[OutPar.minuslogpost].astype(np.float64).idxmin() + MAP = self.data.loc[i_max_logpost].copy() + if not ignore_temperature: + MAP[OutPar.minuslogpost] = self.detempered_minuslogpost( + first=i_max_logpost, last=i_max_logpost + 1)[0] + return MAP def sampled_to_getdist_mcsamples( self, first: Optional[int] = None, last: Optional[int] = None, From 89b31887357ea132a09e458441869c781cceb065 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Wed, 22 Sep 2021 17:49:52 +0200 Subject: [PATCH 10/18] mcmc_tempered: documentation --- docs/post.rst | 6 ++---- docs/sampler_mcmc.rst | 18 +++++++++++++++++- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/docs/post.rst b/docs/post.rst index 00a7b1143..a05e842b1 100644 --- a/docs/post.rst +++ b/docs/post.rst @@ -1,10 +1,6 @@ Importance reweighting and general ``post``-processing ====================================================== -.. note:: - - *new in 1.2* - The `post` component provides a way to post-process an existing sample in different ways: - Add/remove/recompute a prior, e.g. to impose a parameter cut. @@ -16,6 +12,8 @@ The `post` component provides a way to post-process an existing sample in differ The domain of the new pdf (after post-processing) must be well-sampled over the domain of the old one. This may not happen when the new one has a larger domain (specially after removing likelihoods or removing/relaxing priors). In that case, it may be a good idea to redo the full sampling on the new model, instead of post-processing the old one. + To prevent this scenario, if you expect to reweigh your resulting sample, you may want to produce one with a higher-than-1 temperature (see :ref:`mcmc_tempered`). + The requested operations are detailed in a ``post`` block, which may contain one ``add`` and one ``remove`` sub-block. Under ``add``, *using standard input syntax*, specify what you want to add during preprocessing, and under ``remove`` what you would like to remove (no options necessary for removed stuff: just a mention). To force an update of a prior/likelihood/derived-parameter, include it both under ``remove`` and ``add``, with the new options, if needed inside ``add``. The input sample is specified via the ``output`` option with the same value as the original sample. Cobaya will look for it and check that it is compatible with the requested operations. If multiple samples are found (e.g. from an MPI run), all of them are loaded and each processed separately. The resulting samples will have a suffix ``.post.[your_suffix]``, where ``[your_suffix]`` is specified with the ``suffix`` option of ``post`` (you can also change the original prefix [path and base for file name] by setting ``output: [new prefix]`` inside ``post``. diff --git a/docs/sampler_mcmc.rst b/docs/sampler_mcmc.rst index 28c465f3f..201ae261d 100644 --- a/docs/sampler_mcmc.rst +++ b/docs/sampler_mcmc.rst @@ -182,6 +182,23 @@ E.g. if a likelihood depends of parameters ``a``, ``b`` and ``c`` and the cost o If automatic learning of the proposal covariance is enabled, after some checkpoint the proposed steps will mix parameters from different blocks, but *always towards faster ones*. Thus, it is important to specify your blocking in **ascending order of speed**, when not prevented by the architecture of your likelihood (e.g. due to internal caching of intermediate results that require some particular order of parameter variation). + +.. _mcmc_tempered: + +Tempered MCMC +------------- + +Some times it is convenient to sample from a power-reduced (or softened), *tempered* version of the posterior. This produces a Monte Carlo sample with more points towards the tails of the distribution, which is useful when e.g. estimating a quantity by weighting with samples with their probability, or to be able to do more robust :doc:`importance reweighting `. + +By setting a value greater than 1 for the ``temperature`` option, the ``mcmc`` sampler will produce a chain sampled from :math:`p^{1/t}`, where :math:`p` is the posterior and :math:`t` is the temperature. The resulting :class:`~collection.SampleCollection` and output file will contain as weights and log-posterior those of the tempered posterior :math:`p^{1/t}` (this is done like this because it is advantageous to retain integer weights). The original prior- and likelihood-related values in the sample/output are preserved. + +To recover the corresponding weights of the original posterior, multiply the weight column by the output of the :func:`~collection.SampleCollection.detempering_reweight_factor()` method, and the original posterior values are recovered with the :func:`~collection.SampleCollection.detempered_minuslogpost()` method. To remove the temperature from a tempered :class:`~collection.SampleCollection` (i.e. have weights and log-posterior of the original posterior), call the :func:`~collection.SampleCollection.detempered_minuslogpost()` method, possibly on a copy produced with the :func:`~collection.SampleCollection.detempered_copy()` method. + +Despite storing the tempered log-posterior, methods producing statistics such as :func:`~collection.SampleCollection.mean()`, :func:`~collection.SampleCollection.cov()` and :func:`~collection.SampleCollection.MAP()` return the results corresponding to the original posterior, unless they are called with ``ignore_temperature=True``. + +[TODO] Interaction with GetDist. + + .. _mcmc_convergence: Convergence checks @@ -334,4 +351,3 @@ Proposal :members: .. autoclass:: samplers.mcmc.proposal.BlockedProposer :members: - From ad74149afadb05582bfd286b61091c3b8f1ca580 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Wed, 22 Sep 2021 18:11:59 +0200 Subject: [PATCH 11/18] mcmc_tempered: added test --- cobaya/collection.py | 2 +- tests/common_sampler.py | 3 +++ tests/test_mcmc.py | 7 +++++-- tests/test_polychord.py | 2 +- 4 files changed, 10 insertions(+), 4 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index a729c22ad..0c080e7dd 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -307,7 +307,7 @@ def _cache_add_row(self, pos: int, values: Union[Sequence[float], np.ndarray], for name, value in zip(self.chi2_names, logposterior.loglikes): self._cache[pos, self._icol[name]] = -2 * value self._cache[pos, self._icol[OutPar.chi2]] = -2 * logposterior.loglike - if logposterior.derived != []: + if len(logposterior.derived): for name, value in zip(self.derived_params, logposterior.derived): self._cache[pos, self._icol[name]] = value diff --git a/tests/common_sampler.py b/tests/common_sampler.py index 80eb90f9b..9279202cf 100644 --- a/tests/common_sampler.py +++ b/tests/common_sampler.py @@ -92,6 +92,9 @@ def body_of_sampler_test(info_sampler: SamplersDict, dimension=1, n_modes=1, tmp ignore_rows = 0 results = MCSamplesFromCobaya(updated_info, products["sample"], ignore_rows=ignore_rows, name_tag="sample") + temperature = info["sampler"][sampler_name].get("temperature", 1) + if temperature != 1: + results.cool(temperature) clusters = None if "clusters" in products: clusters = [MCSamplesFromCobaya( diff --git a/tests/test_mcmc.py b/tests/test_mcmc.py index aa71df508..324a54b6f 100644 --- a/tests/test_mcmc.py +++ b/tests/test_mcmc.py @@ -18,7 +18,8 @@ @flaky(max_runs=max_runs, min_passes=1) -def test_mcmc(tmpdir, packages_path=None): +@pytest.mark.parametrize("temperature", (1, 2)) +def test_mcmc(tmpdir, temperature, packages_path=None): dimension = 3 # Random initial proposal # cov = mpi.share( @@ -31,7 +32,9 @@ def test_mcmc(tmpdir, packages_path=None): # Bad guess for covmat, so big burn in and max_tries "max_tries": 3000, "burn_in": 100 * dimension, # Proposal, relearn quickly as poor guess - "covmat": cov, "learn_proposal_Rminus1_max": 30}} + "covmat": cov, "learn_proposal_Rminus1_max": 30, + # Sampling temperature, 1 and >1 to be tested + "temperature": temperature}} def check_gaussian(sampler_instance): proposer = KL_norm( diff --git a/tests/test_polychord.py b/tests/test_polychord.py index 23f6233ed..999c74cb5 100644 --- a/tests/test_polychord.py +++ b/tests/test_polychord.py @@ -30,7 +30,7 @@ def test_polychord_resume(packages_path, skip_not_installed, tmpdir): def callback(sampler): nonlocal dead_points - dead_points = sampler.dead[["a", "b"]].values.copy() + dead_points = sampler.dead[["a", "b"]].to_numpy(dtype=np.float64).copy() info = { "likelihood": { From 81f75bed88cdf81cafcca3d695228a5150f93373 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Tue, 23 Nov 2021 00:09:23 +0100 Subject: [PATCH 12/18] typo --- CHANGELOG.md | 5 +++++ docs/theories_and_dependencies.rst | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 83b045f02..f6003e8a5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,11 @@ - Added Boss DR16 likelihood (by @Pablo-Lemos) +#### BICEP/Keck + +- Bugfix in decorrelation function #196 (by Caterina Umilta, @umilta) +- Updated to 2021 data release (2018 data) (by Dominic Beck, @doicbek) + ## 3.1.1 – 2021-07-22 - Changes for compatibility with Pandas 1.3 (which broke convergence testing amongst other things). diff --git a/docs/theories_and_dependencies.rst b/docs/theories_and_dependencies.rst index fa4094069..8f514e637 100644 --- a/docs/theories_and_dependencies.rst +++ b/docs/theories_and_dependencies.rst @@ -31,7 +31,7 @@ The theory code also needs to tell other theory codes and likelihoods the things Use a ``get_X`` method when you need to add optional arguments to provide different outputs from the computed quantity. Quantities returned by :meth:`~.theory.Theory.get_can_provide` should be stored in the state dictionary by the calculate function -or returned by the ``get_results(X)`` for each quantity ``X`` (which by default just returns the value stored in the current state dictionary). +or returned by the ``get_result(X)`` for each quantity ``X`` (which by default just returns the value stored in the current state dictionary). The results stored by calculate for a given set of input parameters are cached, and ``self.current_state`` is set to the current state whenever ``get_X``, ``get_param`` etc are called. From 5dd8be8dbca7302b32e5bda9511c148867236849 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Tue, 23 Nov 2021 08:58:36 +0100 Subject: [PATCH 13/18] docs typo --- docs/installation.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/installation.rst b/docs/installation.rst index bc37dd4ef..a6c725e37 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -71,7 +71,7 @@ To install **cobaya** or upgrade it to the latest release, simply type in a term .. code:: bash - $ python -m pip install cobaya--upgrade + $ python -m pip install cobaya --upgrade To go on installing **cosmological requisites**, see :doc:`installation_cosmo`. From 6dc60cdd1ada134ff5b99f3f9ad0a0824ab77a4e Mon Sep 17 00:00:00 2001 From: torradocacho Date: Tue, 23 Nov 2021 13:21:44 +0100 Subject: [PATCH 14/18] typos and clarifications --- cobaya/theories/classy/classy.py | 5 ++--- cobaya/theories/cosmo/boltzmannbase.py | 13 +++++++------ 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/cobaya/theories/classy/classy.py b/cobaya/theories/classy/classy.py index 890054cd0..1a43b2757 100644 --- a/cobaya/theories/classy/classy.py +++ b/cobaya/theories/classy/classy.py @@ -251,7 +251,6 @@ def must_provide(self, **requirements): # (default: 0.1). But let's leave it like this in case this changes # in the future. self.add_z_for_matter_power(v.pop("z")) - if v["nonlinear"] and "non linear" not in self.extra_args: self.extra_args["non linear"] = non_linear_default_code pair = k[2:] @@ -484,8 +483,8 @@ def get_can_provide_params(self): return names def get_can_support_params(self): - # non-exhaustive list of supported input parameters that will be assigne do classy - # if they are varied + # non-exhaustive list of supported input parameters that will be assigned to + # classy if they are varied return ['H0'] def get_version(self): diff --git a/cobaya/theories/cosmo/boltzmannbase.py b/cobaya/theories/cosmo/boltzmannbase.py index cf83ae190..21fc270f8 100644 --- a/cobaya/theories/cosmo/boltzmannbase.py +++ b/cobaya/theories/cosmo/boltzmannbase.py @@ -314,9 +314,9 @@ def get_Pk_interpolator(self, var_pair=("delta_tot", "delta_tot"), nonlinear=Tru Get a :math:`P(z,k)` bicubic interpolation object (:class:`PowerSpectrumInterpolator`). - In the interpolator returned, both the input :math:`k` and resulting - :math:`P(z,k)` values are in units of :math:`1/\mathrm{Mpc}` (not :math:`h^{-1}` - units). + In the interpolator returned, the input :math:`k` and resulting + :math:`P(z,k)` are in units of :math:`1/\mathrm{Mpc}` and + :math:`\mathrm{Mpc}^3` respectively (not in :math:`h^{-1}` units). :param var_pair: variable pair for power spectrum :param nonlinear: non-linear spectrum (default True) @@ -353,9 +353,10 @@ def get_Pk_grid(self, var_pair=("delta_tot", "delta_tot"), nonlinear=True): Returned arrays may be bigger or more densely sampled than requested, but will include required values. - In the grid returned, both :math:`k` and :math:`P(z,k)` values are in units of - :math:`1/\mathrm{Mpc}` (not :math:`h^{-1}` units), and :math:`z` and :math:`k` - are in **ascending** order. + In the grid returned, :math:`k` and :math:`P(z,k)` are in units of + :math:`1/\mathrm{Mpc}` and :math:`\mathrm{Mpc}^3` respectively + (not in :math:`h^{-1}` units), and :math:`z` and :math:`k` are in + **ascending** order. :param nonlinear: whether the linear or non-linear spectrum :param var_pair: which power spectrum From 4597f507363ce4d0baff2082c6d346942a82c596 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Wed, 24 Nov 2021 16:55:29 +0100 Subject: [PATCH 15/18] model and collection: more typing --- cobaya/collection.py | 4 ++-- cobaya/model.py | 19 +++++++++++++------ 2 files changed, 15 insertions(+), 8 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index 0c080e7dd..320134133 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -13,7 +13,7 @@ import pandas as pd from copy import deepcopy from typing import Union, Sequence, Optional -from getdist import MCSamples, chains +from getdist import MCSamples, chains # type: ignore # Local from cobaya.conventions import OutPar, minuslogprior_names, chi2_names, \ @@ -562,7 +562,7 @@ def thin_samples(self, thin, inplace=False) -> 'SampleCollection': if thin != int(thin) or thin < 1: raise LoggedError(self.log, "Thin factor must be an positive integer, got %s", thin) - from getdist.chains import WeightedSamples, WeightedSampleError + from getdist.chains import WeightedSamples, WeightedSampleError # type: ignore thin = int(thin) try: if hasattr(WeightedSamples, "thin_indices_and_weights"): diff --git a/cobaya/model.py b/cobaya/model.py index 09832f58e..c995cc019 100644 --- a/cobaya/model.py +++ b/cobaya/model.py @@ -63,12 +63,19 @@ class LogPosterior: or likelihoods' chi squared's) is stored to the power of ``1/temperature``. """ - logpost: Optional[float] = None - logpriors: Optional[Sequence[float]] = None - loglikes: Optional[Sequence[float]] = None - derived: Optional[Sequence[float]] = None + # A note on typing: + # Though None is allowed for some arguments, after initialisation everything should + # be not None. So we can either (a) use Optional, and then get A LOT of typing errors + # or (b) not use it (use dataclasses.field(default=None) instead) and get fewer errors + # (only wherever LogPosterior is initialised). + # Let's opt for (b) and suppress errors there. + + logpost: float = dataclasses.field(default=None) # type: ignore + logpriors: Sequence[float] = dataclasses.field(default=None) # type: ignore + loglikes: Sequence[float] = dataclasses.field(default=None) # type: ignore + derived: Sequence[float] = dataclasses.field(default=None) # type: ignore temperature: float = 1 - finite: Optional[bool] = False + finite: bool = dataclasses.field(default=False) logprior: float = dataclasses.field(init=False, repr=False) loglike: float = dataclasses.field(init=False, repr=False) @@ -128,7 +135,7 @@ def make_finite(self): object.__setattr__(self, 'loglikes', np.nan_to_num(self.loglikes)) object.__setattr__(self, 'loglike', np.nan_to_num(self.loglike)) - def as_dict(self, model: "Model") -> Dict[str, float]: + def as_dict(self, model: "Model") -> Dict[str, Union[float, Dict[str, float]]]: """ Given a :class:`~model.Model`, returns a more informative version of itself, containing the names of priors, likelihoods and derived parameters. From b169ec681411f8e16d6f8d93b7195982a25b02c2 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Thu, 25 Nov 2021 17:32:16 +0100 Subject: [PATCH 16/18] collection: covmat computation robustness --- cobaya/collection.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index 320134133..3c24bd16c 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -498,13 +498,14 @@ def cov(self, first: Optional[int] = None, last: Optional[int] = None, if pweight: weight_type_kwarg = "aweights" elif np.allclose(np.round(weights), weights): + weights = np.round(weights).astype(int) weight_type_kwarg = "fweights" else: weight_type_kwarg = "aweights" return np.atleast_2d(np.cov( self[list(self.sampled_params) + (list(self.derived_params) if derived else [])][first:last].to_numpy( - dtype=np.float64).T, + dtype=np.float64).T, **{weight_type_kwarg: weights})) def reweight(self, importance_weights): From 48fc2e6eb6db38062c2f54885fab9683df669f18 Mon Sep 17 00:00:00 2001 From: torradocacho Date: Sun, 5 Dec 2021 07:27:51 +0100 Subject: [PATCH 17/18] prior: added method set_reference to update reference pdf, MPI-aware --- CHANGELOG.md | 5 +- cobaya/model.py | 2 +- cobaya/prior.py | 112 +++++++++++++++++++++++++++++------------- docs/sampler_mcmc.rst | 2 +- tests/test_ref.py | 86 ++++++++++++++++++++++++++++++++ 5 files changed, 171 insertions(+), 36 deletions(-) create mode 100644 tests/test_ref.py diff --git a/CHANGELOG.md b/CHANGELOG.md index f6003e8a5..69d14d6c8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,7 +1,10 @@ ## 3.X.Y – YYYY-MM-DD -- Documented uses of `Model` class in general contexts (previously only cosmo) +### General + +- Documented uses of `Model` class in general contexts (previously only cosmo). - `Model` methods to compute log-probabilities and derived parameters now have an `as_dict` keyword (default `False`), for more informative return value. +- `Prior` now has method `set_reference`, to update the reference pdf's if needed (MPI-aware, documented). ### Cosmological likelihoods and theory codes diff --git a/cobaya/model.py b/cobaya/model.py index c995cc019..d773f3a7e 100644 --- a/cobaya/model.py +++ b/cobaya/model.py @@ -638,7 +638,7 @@ def get_valid_point(self, max_tries: int, ignore_fixed_ref: bool = False, if results.logpost != -np.inf: break else: - if self.prior.reference_is_pointlike(): + if self.prior.reference_is_pointlike: raise LoggedError(self.log, "The reference point provided has null " "likelihood. Set 'ref' to a different point " "or a pdf.") diff --git a/cobaya/prior.py b/cobaya/prior.py index ccfe7278f..4801d8a90 100644 --- a/cobaya/prior.py +++ b/cobaya/prior.py @@ -29,16 +29,21 @@ (e.g. to guess the pdf tails correctly if the derived quantity needs to be positive or negative -- defaulted to ``-.inf``, ``.inf`` resp.). -The (optional) **reference** pdf (``ref``) for **sampled** parameters -defines the region of the prior which is of most -interest (e.g. where most of the prior mass is expected); -samplers needing an initial point -for a chain will attempt to draw it from the ``ref`` if it has been defined (otherwise -from the prior). A good reference pdf will avoid a long *burn-in* stage during the -sampling. If you assign a single value to ``ref``, samplers will always start from -that value; however this makes convergence tests less reliable as each chain will -start from the same point (so all chains could be stuck near the same point). +The (optional) **reference** pdf (``ref``) for **sampled** parameters defines the region +of the prior which is of most interest (e.g. where most of the prior mass is expected); +samplers needing an initial point for a chain will attempt to draw it from the ``ref`` +if it has been defined (otherwise from the prior). A good reference pdf will avoid a long +*burn-in* stage during the sampling. If you assign a single value to ``ref``, samplers +will always start from that value; however this makes convergence tests less reliable as +each chain will start from the same point (so all chains could be stuck near the same +point). +.. note:: + + When running in parallel with MPI you can give ``ref`` different values for each of + the different parallel processes (either at initialisation or by calling + :func:`Prior.set_reference`). They will be taken into account e.g. for starting + parallel MCMC chains at different points/distributions of your choice. The syntax for priors and ref's has the following fields: @@ -376,8 +381,6 @@ def __init__(self, parameterization: Parameterization, # in principle, separable: one per parameter self.params = [] self.pdf = [] - self.ref_pdf = [] - self._ref_is_pointlike = True self._bounds = np.zeros((len(sampled_params_info), 2)) for i, p in enumerate(sampled_params_info): self.params += [p] @@ -386,26 +389,6 @@ def __init__(self, parameterization: Parameterization, fast_logpdf = fast_logpdfs.get(self.pdf[-1].dist.name) if fast_logpdf: self.pdf[-1].logpdf = MethodType(fast_logpdf, self.pdf[-1]) - # Get the reference (1d) pdf - ref = sampled_params_info[p].get("ref") - # Cases: number, pdf (something, but not a number), nothing - if isinstance(ref, Sequence) and len(ref) == 2 and all( - isinstance(n, numbers.Number) for n in ref): - ref = {"dist": "norm", "loc": ref[0], "scale": ref[1]} - if isinstance(ref, numbers.Real): - self.ref_pdf += [float(ref)] - elif isinstance(ref, Mapping): - self.ref_pdf += [get_scipy_1d_pdf({p: ref})] - self._ref_is_pointlike = False - elif ref is None: - self.ref_pdf += [np.nan] - self._ref_is_pointlike = False - else: - raise LoggedError(self.log, - "'ref' for starting position should be None or a number" - ", a list of two numbers for normal mean and deviation," - "or a dict with parameters for a scipy distribution.") - self._bounds[i] = [-np.inf, np.inf] try: self._bounds[i] = self.pdf[-1].interval(1) @@ -423,7 +406,8 @@ def __init__(self, parameterization: Parameterization, self._lower_limits = self._bounds[:, 0].copy() self._uniform_logp = -np.sum(np.log(self._upper_limits[self._uniform_indices] - self._lower_limits[self._uniform_indices])) - + # Set the reference pdf's + self.set_reference({p: v.get("ref") for p, v in sampled_params_info.items()}) # Process the external prior(s): self.external = {} self.external_dependence = set() @@ -590,7 +574,69 @@ def covmat(self, ignore_external=False): "It is not possible to get the covariance matrix from an external prior.") return np.diag([pdf.var() for pdf in self.pdf]).T + def set_reference(self, ref_info): + """ + Sets or updates the reference pdf with the given parameter input info. + + ``ref_info`` should be a dict ``{parameter_name: [ref definition]}``, not + ``{parameter_name: {"ref": [ref definition]}}``. + + When called after prior initialisation, not mentioning a parameter leaves + its reference pdf unchanged, whereas explicitly setting ``ref: None`` sets + the prior as the reference pdf. + + You can set different reference pdf's/values for different MPI processes, + e.g. for fixing different starting points for parallel MCMC chains. + """ + if not hasattr(self, "ref_pdf"): + # Initialised with nan's in case ref==None: no ref -> uses prior + self.ref_pdf = [np.nan] * self.d() + unknown = set(ref_info).difference(self.params) + if unknown: + raise LoggedError(self.log, + f"Cannot set reference pdf for parameter(s) {unknown}: " + "not sampled parameters.") + for i, p in enumerate(self.params): + # The next if ensures correct behavious in "update call", + # where not mentioning a parameter and making its ref None are different + # (not changing vs setting to prior) + if p not in ref_info: + continue # init: use prior; update: don't change + ref = ref_info[p] + # [number, number] interpreted as Gaussian + if isinstance(ref, Sequence) and len(ref) == 2 and all( + isinstance(n, numbers.Number) for n in ref): + ref = {"dist": "norm", "loc": ref[0], "scale": ref[1]} + if isinstance(ref, numbers.Real): + self.ref_pdf[i] = float(ref) + elif isinstance(ref, Mapping): + self.ref_pdf[i] = get_scipy_1d_pdf({p: ref}) + elif ref is None: + # We only get here if explicit `param: None` mention! + self.ref_pdf[i] = np.nan + else: + raise LoggedError(self.log, + "'ref' for starting position should be None or a number" + ", a list of two numbers for normal mean and deviation," + "or a dict with parameters for a scipy distribution.") + # Re-set the pointlike-ref property + if hasattr(self, "_ref_is_pointlike"): + delattr(self, "_ref_is_pointlike") + self.reference_is_pointlike + +# TODO: check for errors due to the change into property +# TODO: check for errors in test for pointlike-ness + @property def reference_is_pointlike(self) -> bool: + """ + Whether there is a fixed reference point for all parameters, such that calls to + :func:`Prior.reference` would always return the same. + """ + if not hasattr(self, "_ref_is_pointlike"): + self._ref_is_pointlike = all( + # np.nan is a numbers.Number instance, but not a fixed ref (uses prior) + (isinstance(ref, numbers.Number) and ref is not np.nan) + for ref in self.ref_pdf) return self._ref_is_pointlike def reference(self, max_tries=np.inf, warn_if_tries="10d", ignore_fixed=False, @@ -643,7 +689,7 @@ def reference(self, max_tries=np.inf, warn_if_tries="10d", ignore_fixed=False, "If stuck here, maybe it is not possible to sample from the " "reference pdf a point with non-null prior. Check that the reference " "pdf and the prior are consistent.") - if self.reference_is_pointlike(): + if self.reference_is_pointlike: raise LoggedError(self.log, "The reference point provided has null prior. " "Set 'ref' to a different point or a pdf.") raise LoggedError( diff --git a/docs/sampler_mcmc.rst b/docs/sampler_mcmc.rst index 201ae261d..db95fa34f 100644 --- a/docs/sampler_mcmc.rst +++ b/docs/sampler_mcmc.rst @@ -208,7 +208,7 @@ Convergence of an MCMC run is assessed in terms a generalized version of the :ma In particular, given a small number of chains from the same run, the :math:`R-1` statistic measures (from the last half of each chain) the variance between the means of the different chains in units of the covariance of the chains (in other words, that all chains are centered around the same point, not deviating from it a significant fraction of the standard deviation of the posterior). When that number becomes smaller than ``Rminus1_stop`` twice in a row, a second :math:`R-1` check is also performed on the bounds of the ``Rminus1_cl_level`` % confidence level interval, which, if smaller than ``Rminus1_cl_stop``, stops the run. -The default settings have *Rminus1_stop = 0.01*, *Rminus1_cl_level = 0.95* and *Rminus1_cl_level = 0.2*; the stop values can be decreased for better convergence. +The default settings are ``Rminus1_stop = 0.01``, ``Rminus1_cl_level = 0.95`` and ``Rminus1_cl_level = 0.2``; the stop values can be decreased for better convergence. .. note:: diff --git a/tests/test_ref.py b/tests/test_ref.py new file mode 100644 index 000000000..68722d3e3 --- /dev/null +++ b/tests/test_ref.py @@ -0,0 +1,86 @@ +""" +Tests the setting and updateing of the reference pdf, including re-checking if point-like. +""" +import pytest +import numpy as np + +from cobaya.model import get_model +from cobaya.sampler import get_sampler +from cobaya import mpi + + +def test_ref(): + val = 1 + mean, std = 0.5, 0.1 + info = { + "params": { + "a": {"prior": [0, 1]}, + "b": {"prior": [0, 1], "ref": None}, + "c": {"prior": [0, 1], "ref": val}, + "d": {"prior": [0, 1], "ref": [mean, std]}, + "e": {"prior": [0, 1], "ref": {"dist": "norm", "loc": mean, "scale": std}}}, + "likelihood": {"one": None}} + model = get_model(info) + for i in [3, 4]: + assert model.prior.ref_pdf[i].dist.name == "norm" + assert model.prior.ref_pdf[i].mean() == mean + assert model.prior.ref_pdf[i].std() == std + # Not point-like: 2 nan's (use prior) and 2 norms + assert not model.prior.reference_is_pointlike + # Let's update it -- remove norms + upd1_ref_info = { + "d": val + 2, + "e": val + 3} + model.prior.set_reference(upd1_ref_info) + # Updated: + assert model.prior.ref_pdf[3] == upd1_ref_info["d"] + assert model.prior.ref_pdf[4] == upd1_ref_info["e"] + # Unchanged: + assert model.prior.ref_pdf[0] is np.nan + assert model.prior.ref_pdf[1] is np.nan + assert model.prior.ref_pdf[2] == val + # Still not pointlike: there are nan's --> use prior + assert not model.prior.reference_is_pointlike + # Let's update it -- remove nan's + upd2_ref_info = { + "a": val - 2, + "b": val - 1} + model.prior.set_reference(upd2_ref_info) + # Updated: + assert model.prior.ref_pdf[0] == upd2_ref_info["a"] + assert model.prior.ref_pdf[1] == upd2_ref_info["b"] + # Unchanged: + assert model.prior.ref_pdf[2] == val + assert model.prior.ref_pdf[3] == upd1_ref_info["d"] + assert model.prior.ref_pdf[4] == upd1_ref_info["e"] + # Should be point-like now! + assert model.prior.reference_is_pointlike + # Let's update it -- back to one norm + upd3_ref_info = { + "a": [mean, std]} + model.prior.set_reference(upd3_ref_info) + # Updated: + assert model.prior.ref_pdf[0].dist.name == "norm" + assert model.prior.ref_pdf[0].mean() == mean + assert model.prior.ref_pdf[0].std() == std + # Unchanged: + assert model.prior.ref_pdf[1] == upd2_ref_info["b"] + assert model.prior.ref_pdf[2] == val + assert model.prior.ref_pdf[3] == upd1_ref_info["d"] + assert model.prior.ref_pdf[4] == upd1_ref_info["e"] + # Should be non-point-like again + assert not model.prior.reference_is_pointlike + + +# Tests MPI-awareness of the referece, when using it to get the initial point of mcmc +@pytest.mark.mpi +@mpi.sync_errors +def test_ref_mcmc_initial_point(): + val = 0.5 + info = { + "params": {"a": {"prior": [0, 1 + mpi.size()], "ref": val + mpi.rank()}}, + "likelihood": {"one": None}} + model = get_model(info) + mcmc_sampler = get_sampler({"mcmc": None}, model) + initial_point = mcmc_sampler.current_point.values[0] + assert initial_point == val + mpi.rank() From 8e8a5f42c723dc86c3ff9371f39dc30a220f7a56 Mon Sep 17 00:00:00 2001 From: Antony Lewis Date: Thu, 6 Oct 2022 18:16:01 +0100 Subject: [PATCH 18/18] trivial --- cobaya/collection.py | 11 ++++++----- cobaya/theories/camb/camb.py | 2 +- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/cobaya/collection.py b/cobaya/collection.py index ce989def0..5fcef081f 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -448,14 +448,14 @@ def _weights_for_stats(self, first: Optional[int] = None, last: Optional[int] = """ if pweight: if ignore_temperature: - logps = -self[OutPar.minuslogpost][first:last]\ + logps = -self[OutPar.minuslogpost][first:last] \ .to_numpy(dtype=np.float64, copy=True) else: logps = -self.detempered_minuslogpost(first, last) logps -= max(logps) weights = np.exp(logps) else: - weights = self[OutPar.weight][first:last]\ + weights = self[OutPar.weight][first:last] \ .to_numpy(dtype=np.float64, copy=True) if not ignore_temperature: weights *= self.detempering_reweight_factor(first, last) @@ -512,7 +512,7 @@ def cov(self, first: Optional[int] = None, last: Optional[int] = None, return np.atleast_2d(np.cov( self[list(self.sampled_params) + (list(self.derived_params) if derived else [])][first:last].to_numpy( - dtype=np.float64).T, + dtype=np.float64).T, **{weight_type_kwarg: weights})) def reweight(self, importance_weights): @@ -542,8 +542,9 @@ def detempered_minuslogpost(self, first: Optional[int] = None, """ if self.temperature != 1: return ( - self.data[OutPar.minuslogprior][first:last].to_numpy(dtype=np.float64) + - self.data[OutPar.chi2][first:last].to_numpy(dtype=np.float64) / 2) + self.data[OutPar.minuslogprior][first:last].to_numpy( + dtype=np.float64) + + self.data[OutPar.chi2][first:last].to_numpy(dtype=np.float64) / 2) else: return self[OutPar.minuslogpost][first:last].to_numpy( dtype=np.float64, copy=True) diff --git a/cobaya/theories/camb/camb.py b/cobaya/theories/camb/camb.py index 3e6d675b8..d1890712e 100644 --- a/cobaya/theories/camb/camb.py +++ b/cobaya/theories/camb/camb.py @@ -810,7 +810,7 @@ def set(self, params_values_dict, state): else: raise LoggedError( self.log, - "Some of the attributes to be set manually were not " + "Some attributes to be set manually were not " "recognized: %s=%s", attr, value) # Sources if getattr(self, "sources", None):