# %% IMPORTS
# Built-in imports
import logging as log
from os import path
# Package imports
import pymultinest
import numpy as np
import os
# IMAGINE imports
from imagine.pipelines import Pipeline
# All declaration
__all__ = ['MultinestPipeline']
# %% CLASS DEFINITIONS
[docs]class MultinestPipeline(Pipeline):
"""
Bayesian analysis pipeline with
`pyMultinest <https://github.com/JohannesBuchner/PyMultiNest>`_
See base class for initialization details.
The sampler behaviour is controlled using the `sampling_controllers`
property. A description of these can be found below.
Sampling controllers
--------------------
resume : bool
If False the Pipeline the sampling starts from the beginning,
overwriting any previous work in the `chains_directory`. Otherwise,
tries to resume a previous run.
n_live_points : int
Number of live points to be used.
evidence_tolerance : float
A value of 0.5 should give good enough accuracy.
max_iter : int
Maximum number of iterations. `0` (default) is unlimited (i.e. only
stops after convergence).
log_zero : float
Points with loglike < logZero will be ignored by MultiNest
importance_nested_sampling : bool
If `True` (default), Multinest will use Importance Nested
Sampling (see `arXiv:1306.2144 <http://arxiv.org/abs/1306.2144>`_)
sampling_efficiency : float
Efficieny of the sampling. `0.8` (default) and `0.3` are recommended
values for parameter estimation & evidence evaluation respectively.
multimodal : bool
If `True`, MultiNest will attempt to separate out the
modes using a clustering algorithm.
mode_tolerance : float
MultiNest can find multiple modes and specify which samples belong to
which mode. It might be desirable to have separate samples and mode
statistics for modes with local log-evidence value greater than a
particular value in which case `mode_tolerance` should be set to that
value.
If there isn't any particularly interesting `mode_tolerance` value,
then it should be set to a very negative number (e.g. -1e90, default).
null_log_evidence : float
If `multimodal` is `True`, MultiNest can find multiple modes and also
specify which samples belong to which mode. It might be
desirable to have separate samples and mode statistics for modes
with local log-evidence value greater than a
particular value in which case nullZ should be set to that
value. If there isn't any particulrly interesting
nullZ value, then nullZ should be set to a very large negative
number (e.g. -1.d90).
n_clustering_params : int
Mode separation is done through a clustering
algorithm. Mode separation can be done on all the parameters
(in which case nCdims should be set to ndims) & it
can also be done on a subset of parameters (in which case
nCdims < ndims) which might be advantageous as
clustering is less accurate as the dimensionality increases.
If nCdims < ndims then mode separation is done on
the first nCdims parameters.
max_modes : int
Maximum number of modes (if `multimodal` is `True`).
Note
----
Instances of this class are callable.
Look at the :py:meth:`MultinestPipeline.call` for details.
"""
# Class attributes
SUPPORTS_MPI = True
[docs] def call(self, **kwargs):
"""
Runs the IMAGINE pipeline using the MultiNest sampler
Returns
-------
results : dict
pyMultinest sampling results in a dictionary containing the keys:
logZ (the log-evidence), logZerror (the error in log-evidence) and
samples (equal weighted posterior)
"""
log.debug('@ multinest_pipeline::__call__')
default_solve_params = {'resume': True,
'n_live_points': 400,
'evidence_tolerance': 0.5,
'max_iter': 0,
'log_zero': -1e100,
'importance_nested_sampling': True,
'sampling_efficiency': 0.8,
'multimodal': True,
'mode_tolerance': -1e90,
'null_log_evidence': -1e90,
'n_clustering_params': None,
'outputfiles_basename': None,
'max_modes': 100,
'n_iter_before_update': 100,
'verbose': True}
# Keyword arguments can alter the sampling controllers
self.sampling_controllers = kwargs # Updates the dict
# Sets base name for multinest output files
chains_prefix = path.join(self.chains_directory, 'multinest_')
self.sampling_controllers['outputfiles_basename'] = chains_prefix
# Prepares initialization and run parameters from
# defaults and sampling controllers
solve_params = {k: self.sampling_controllers.get(k, default)
for k, default in default_solve_params.items()}
# Updates the sampling controllers to reflect what is being used
self.sampling_controllers = solve_params # Updates the dict
if not self.sampling_controllers['resume']:
self.clean_chains_directory()
# Runs pyMultinest
log.info('Calling pymultinest.solve')
log.info('Master seed: {}'.format(self.master_seed))
# Adjusts the seed to the correct interval, if needed
if 0 < self.master_seed < 30081:
seed = self.master_seed
else:
seed = np.random.randint(0, 30081)
log.info('MultiNest random seed: {}'.format(seed))
self.results = pymultinest.solve(LogLikelihood=self._likelihood_function,
Prior=self.prior_transform,
n_dims=len(self._active_parameters),
write_output=True,
seed=seed,
wrapped_params=self.wrapped_parameters,
**solve_params)
self._samples_array = self.results['samples']
self._evidence = self.results['logZ']
self._evidence_err = self.results['logZerr']
return self.results