Merge branch 'polychord' into 'master'

Polychord

See merge request !264

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Unreleased
 
 ### Added
+- Added the PolyChord sampler, which can be accessed by using `sampler='pypolychord'` in `run_sampler`
 - `emcee` now writes all progress to disk and can resume from a previous run.
 
 ### Changed

Pipfile

@@ -6,11 +6,12 @@ name = "pypi"
 [packages]
 future = "*"
 corner = "*"
-numpy = ">=1.9"
+numpy = "==1.15.2"
+ligotimegps = "<=1.2.3"
 matplotlib = "<3"
 scipy = ">=0.16"
-pandas = "*"
-deepdish = "*"
+pandas = "==0.23.0"
+deepdish = "==0.3.6"
 mock = "*"
 astropy = "<3"
 gwpy = "*"
@@ -23,7 +24,6 @@ emcee = "*"
 nestle = "*"
 ptemcee = "*"
 pymc3 = "*"
-pymultinest = "*"
 
 [requires]
 

bilby/core/result.py

@@ -157,7 +157,15 @@ class Result(object):
             else:
                 filename = result_file_name(outdir, label)
         if os.path.isfile(filename):
-            return cls(**deepdish.io.load(filename))
+            dictionary = deepdish.io.load(filename)
+            # Some versions of deepdish/pytables return the dictionanary as
+            # a dictionary with a kay 'data'
+            if len(dictionary) == 1 and 'data' in dictionary:
+                dictionary = dictionary['data']
+            try:
+                return cls(**dictionary)
+            except TypeError as e:
+                raise IOError("Unable to load dictionary, error={}".format(e))
         else:
             raise IOError("No result '{}' found".format(filename))
 

bilby/core/sampler/__init__.py

@@ -11,20 +11,21 @@ from .cpnest import Cpnest
 from .dynesty import Dynesty
 from .emcee import Emcee
 from .nestle import Nestle
+from .polychord import PyPolyChord
 from .ptemcee import Ptemcee
 from .ptmcmc import PTMCMCSampler
 from .pymc3 import Pymc3
 from .pymultinest import Pymultinest
 
-implemented_samplers = {
+IMPLEMENTED_SAMPLERS = {
     'cpnest': Cpnest, 'dynesty': Dynesty, 'emcee': Emcee, 'nestle': Nestle,
     'ptemcee': Ptemcee,'ptmcmcsampler' : PTMCMCSampler,
-    'pymc3': Pymc3, 'pymultinest': Pymultinest }
+    'pymc3': Pymc3, 'pymultinest': Pymultinest, 'pypolychord': PyPolyChord }
 
 if command_line_args.sampler_help:
     sampler = command_line_args.sampler_help
-    if sampler in implemented_samplers:
-        sampler_class = implemented_samplers[sampler]
+    if sampler in IMPLEMENTED_SAMPLERS:
+        sampler_class = IMPLEMENTED_SAMPLERS[sampler]
         print('Help for sampler "{}":'.format(sampler))
         print(sampler_class.__doc__)
     else:
@@ -33,7 +34,7 @@ if command_line_args.sampler_help:
                   'the name of the sampler')
         else:
             print('Requested sampler {} not implemented'.format(sampler))
-        print('Available samplers = {}'.format(implemented_samplers))
+        print('Available samplers = {}'.format(IMPLEMENTED_SAMPLERS))
 
     sys.exit()
 
@@ -106,7 +107,7 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
     if command_line_args.clean:
         kwargs['resume'] = False
 
-    from . import implemented_samplers
+    from . import IMPLEMENTED_SAMPLERS
 
     if priors is None:
         priors = dict()
@@ -128,15 +129,15 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
     if isinstance(sampler, Sampler):
         pass
     elif isinstance(sampler, str):
-        if sampler.lower() in implemented_samplers:
-            sampler_class = implemented_samplers[sampler.lower()]
+        if sampler.lower() in IMPLEMENTED_SAMPLERS:
+            sampler_class = IMPLEMENTED_SAMPLERS[sampler.lower()]
             sampler = sampler_class(
                 likelihood, priors=priors, outdir=outdir, label=label,
                 injection_parameters=injection_parameters, meta_data=meta_data,
                 use_ratio=use_ratio, plot=plot, result_class=result_class,
                 **kwargs)
         else:
-            print(implemented_samplers)
+            print(IMPLEMENTED_SAMPLERS)
             raise ValueError(
                 "Sampler {} not yet implemented".format(sampler))
     elif inspect.isclass(sampler):
@@ -148,7 +149,7 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
     else:
         raise ValueError(
             "Provided sampler should be a Sampler object or name of a known "
-            "sampler: {}.".format(', '.join(implemented_samplers.keys())))
+            "sampler: {}.".format(', '.join(IMPLEMENTED_SAMPLERS.keys())))
 
     if sampler.cached_result:
         logger.warning("Using cached result")

bilby/core/sampler/polychord.py

+from __future__ import absolute_import
+
+import numpy as np
+
+from .base_sampler import NestedSampler
+
+
+class PyPolyChord(NestedSampler):
+
+    """
+    Bilby wrapper of PyPolyChord
+    https://arxiv.org/abs/1506.00171
+
+    PolyChordLite is available at:
+    https://github.com/PolyChord/PolyChordLite
+
+    Follow the installation instructions at their github page.
+
+    Keyword arguments will be passed into `pypolychord.run_polychord` into the `settings`
+    argument. See the PolyChord documentation for what all of those mean.
+
+    To see what the keyword arguments are for, see the docstring of PyPolyChordSettings
+    """
+
+    default_kwargs = dict(use_polychord_defaults=False, nlive=None, num_repeats=None,
+                          nprior=-1, do_clustering=True, feedback=1, precision_criterion=0.001,
+                          logzero=-1e30, max_ndead=-1, boost_posterior=0.0, posteriors=True,
+                          equals=True, cluster_posteriors=True, write_resume=True,
+                          write_paramnames=False, read_resume=True, write_stats=True,
+                          write_live=True, write_dead=True, write_prior=True,
+                          compression_factor=np.exp(-1), base_dir='outdir',
+                          file_root='polychord', seed=-1, grade_dims=None, grade_frac=None, nlives={})
+
+    def run_sampler(self):
+        import pypolychord
+        from pypolychord.settings import PolyChordSettings
+        if self.kwargs['use_polychord_defaults']:
+            settings = PolyChordSettings(nDims=self.ndim, nDerived=self.ndim, base_dir=self.outdir,
+                                         file_root=self.label)
+        else:
+            self._setup_dynamic_defaults()
+            pc_kwargs = self.kwargs.copy()
+            pc_kwargs['base_dir'] = self.outdir
+            pc_kwargs['file_root'] = self.label
+            pc_kwargs.pop('use_polychord_defaults')
+            settings = PolyChordSettings(nDims=self.ndim, nDerived=self.ndim, **pc_kwargs)
+        self._verify_kwargs_against_default_kwargs()
+
+        pypolychord.run_polychord(loglikelihood=self.log_likelihood, nDims=self.ndim,
+                                  nDerived=self.ndim, settings=settings, prior=self.prior_transform)
+
+        self.result.log_evidence, self.result.log_evidence_err = self._read_out_stats_file()
+        self.result.samples = self._read_sample_file()
+        return self.result
+
+    def _setup_dynamic_defaults(self):
+        """ Sets up some interdependent default argument if none are given by the user """
+        if not self.kwargs['grade_dims']:
+            self.kwargs['grade_dims'] = [self.ndim]
+        if not self.kwargs['grade_frac']:
+            self.kwargs['grade_frac'] = [1.0] * len(self.kwargs['grade_dims'])
+        if not self.kwargs['nlive']:
+            self.kwargs['nlive'] = self.ndim * 25
+        if not self.kwargs['num_repeats']:
+            self.kwargs['num_repeats'] = self.ndim * 5
+
+    def _translate_kwargs(self, kwargs):
+        if 'nlive' not in kwargs:
+            for equiv in self.npoints_equiv_kwargs:
+                if equiv in kwargs:
+                    kwargs['nlive'] = kwargs.pop(equiv)
+
+    def log_likelihood(self, theta):
+        """ Overrides the log_likelihood so that PolyChord understands it """
+        return super(PyPolyChord, self).log_likelihood(theta), theta
+
+    def _read_out_stats_file(self):
+        statsfile = self.outdir + '/' + self.label + '.stats'
+        with open(statsfile) as f:
+            for line in f:
+                if line.startswith('log(Z)'):
+                    line = line.replace('log(Z)', '')
+                    line = line.replace('=', '')
+                    line = line.replace(' ', '')
+                    print(line)
+                    z = line.split('+/-')
+                    log_z = float(z[0])
+                    log_z_err = float(z[1])
+                    return log_z, log_z_err
+
+    def _read_sample_file(self):
+        sample_file = self.outdir + '/' + self.label + '_equal_weights.txt'
+        samples = np.loadtxt(sample_file)
+        return samples[:, -self.ndim:]  # extract last ndim columns

docs/samplers.txt

@@ -48,7 +48,7 @@ you would use
    is given in the API information below.
 
 Below, we give the detailed API for the samplers. Remember, this API is not
-recomended for direct use by the user, rather it should be accessed via the
+recommended for direct use by the user, rather it should be accessed via the
 :code:`run_sampler`.
 
 ---------------
@@ -60,6 +60,7 @@ Dynesty
 
 .. autoclass:: bilby.core.sampler.dynesty.Dynesty
 
+
 Nestle
 ======
 
@@ -71,11 +72,19 @@ CPNest
 
 .. autoclass:: bilby.core.sampler.cpnest.Cpnest
 
+
 PyMultinest
 ===========
 
 .. autoclass:: bilby.core.sampler.pymultinest.Pymultinest
 
+
+PyPolyChord
+===========
+
+.. autoclass:: bilby.core.sampler.pypolychord.PyPolyChord
+
+
 -------------
 MCMC samplers
 -------------
@@ -123,6 +132,27 @@ be found in the at the top-level of `the repository
 <https://git.ligo.org/lscsoft/bilby>`_ (Note: if you installed from pip, you
 can simply download that file and use the command above).
 
+
+Installing PyPolyChord
+======================
+
+If you want to use the `PyPolyChord` sampler, you first need the
+PolyChord library to be installed to work properly. An image of PolyChord can be found on github.
+Clone the following repository onto your system. Navigate to the folder you want to install PolyChord in and run:
+
+.. code-block:: console
+
+   $ git clone https://github.com/PolyChord/PolyChordLite.git
+
+Then navigate into the PolyChord directory and install PolyChord/PyPolyChord with
+
+.. code-block:: console
+
+   $ make pypolychord MPI=
+   $ python setup.py install --user
+
+Add a number after `MPI=` to compile with `MPI`. Leave it like it is if you don't wish to compile with MPI.
+
 Installing pymultinest
 ======================
 

test/sampler_test.py

@@ -5,6 +5,8 @@ import unittest
 from mock import MagicMock
 import numpy as np
 import os
+import sys
+import shutil
 import copy
 
 
@@ -247,6 +249,53 @@ class TestNestle(unittest.TestCase):
             self.assertDictEqual(expected, self.sampler.kwargs)
 
 
+class TestPolyChord(unittest.TestCase):
+
+    def setUp(self):
+        self.likelihood = MagicMock()
+        self.priors = dict(a=bilby.prior.Uniform(0, 1))
+        self.sampler = bilby.core.sampler.PyPolyChord(self.likelihood, self.priors,
+                                                      outdir='outdir', label='polychord',
+                                                      use_ratio=False, plot=False,
+                                                      skip_import_verification=True)
+
+    def tearDown(self):
+        del self.likelihood
+        del self.priors
+        del self.sampler
+
+    def test_default_kwargs(self):
+        expected = dict(use_polychord_defaults=False, nlive=self.sampler.ndim*25, num_repeats=self.sampler.ndim*5,
+                        nprior=-1, do_clustering=True, feedback=1, precision_criterion=0.001,
+                        logzero=-1e30, max_ndead=-1, boost_posterior=0.0, posteriors=True,
+                        equals=True, cluster_posteriors=True, write_resume=True,
+                        write_paramnames=False, read_resume=True, write_stats=True,
+                        write_live=True, write_dead=True, write_prior=True,
+                        compression_factor=np.exp(-1), base_dir='outdir',
+                        file_root='polychord', seed=-1, grade_dims=list([self.sampler.ndim]),
+                        grade_frac=list([1.0]*len([self.sampler.ndim])), nlives={})
+        self.sampler._setup_dynamic_defaults()
+        self.assertDictEqual(expected, self.sampler.kwargs)
+
+    def test_translate_kwargs(self):
+        expected = dict(use_polychord_defaults=False, nlive=123, num_repeats=self.sampler.ndim*5,
+                        nprior=-1, do_clustering=True, feedback=1, precision_criterion=0.001,
+                        logzero=-1e30, max_ndead=-1, boost_posterior=0.0, posteriors=True,
+                        equals=True, cluster_posteriors=True, write_resume=True,
+                        write_paramnames=False, read_resume=True, write_stats=True,
+                        write_live=True, write_dead=True, write_prior=True,
+                        compression_factor=np.exp(-1), base_dir='outdir',
+                        file_root='polychord', seed=-1, grade_dims=list([self.sampler.ndim]),
+                        grade_frac=list([1.0]*len([self.sampler.ndim])), nlives={})
+        self.sampler._setup_dynamic_defaults()
+        for equiv in bilby.core.sampler.base_sampler.NestedSampler.npoints_equiv_kwargs:
+            new_kwargs = self.sampler.kwargs.copy()
+            del new_kwargs['nlive']
+            new_kwargs[equiv] = 123
+            self.sampler.kwargs = new_kwargs
+            self.assertDictEqual(expected, self.sampler.kwargs)
+
+
 class TestPTEmcee(unittest.TestCase):
 
     def setUp(self):
@@ -397,11 +446,13 @@ class TestRunningSamplers(unittest.TestCase):
 
         self.priors = dict(
             m=bilby.core.prior.Uniform(0, 5), c=bilby.core.prior.Uniform(-2, 2))
+        bilby.core.utils.check_directory_exists_and_if_not_mkdir('outdir')
 
     def tearDown(self):
         del self.likelihood
         del self.priors
         bilby.core.utils.command_line_args.test = False
+        shutil.rmtree('outdir')
 
     def test_run_cpnest(self):
         _ = bilby.run_sampler(
@@ -423,6 +474,11 @@ class TestRunningSamplers(unittest.TestCase):
             likelihood=self.likelihood, priors=self.priors, sampler='nestle',
             nlive=100, save=False)
 
+    def test_run_pypolychord(self):
+        _ = bilby.run_sampler(
+            likelihood=self.likelihood, priors=self.priors,
+            sampler='pypolychord', nlive=100, save=False)
+
     def test_run_ptemcee(self):
         _ = bilby.run_sampler(
             likelihood=self.likelihood, priors=self.priors, sampler='ptemcee',