Nested Sampling posterior sampler

examples/nested_sampling_hds_interface.py

@@ -0,0 +1,23 @@
+"""
+Example of a posterior sampling experiment. Implemented procedure is explained
+at https://en.wikipedia.org/wiki/Rejection_sampling
+"""
+import high_dimensional_sampling as hds
+import numpy as np
+from high_dimensional_sampling.posterior.nestedsampling import PolyChord
+
+procedure = PolyChord()
+experiment = hds.PosteriorSamplingExperiment(procedure, './hds')
+feeder = hds.functions.FunctionFeeder()
+feeder.load_function_group(
+    'posterior', {
+        "Block": {
+            "block_size": 8
+        },
+        "MultivariateNormal": {
+            "covariance": [[4, 0], [0, 4]]
+        }
+    })
+
+for function in feeder:
+    experiment.run(function, finish_line=200)

high_dimensional_sampling/functions.py

@@ -1063,22 +1063,26 @@ class Block(TestFunction):
             dimension as [-block_size, block_size]. Default: 1.
         block_value: Value that the function takes *inside* the ranges spanned
             by block_size. Default: 1.
+        global_size: Defines the domain for which the block function should be
+        searched [-global_size, global_size]. Default: 10.
         global_value: Value that the function takes outside of the ranges
             spanned by block_size. Default: 0.
 
     """
 
     def __init__(self,
                  dimensionality=3,
-                 block_size=1,
-                 block_value=1,
-                 global_value=0,
+                 block_size=1.,
+                 block_value=1.,
+                 global_size=10.,
+                 global_value=0.,
                  **kwargs):
         self.dimensionality = dimensionality
         self.block_size = block_size
         self.block_value = block_value
         self.global_value = global_value
-        self.ranges = self.construct_ranges(dimensionality, -np.inf, np.inf)
+        self.ranges = self.construct_ranges(dimensionality,
+                                            -global_size, global_size)
         super(Block, self).__init__(**kwargs)
 
     def _evaluate(self, x):

high_dimensional_sampling/posterior/nestedsampling.py

@@ -0,0 +1,65 @@
+import os
+import numpy as np
+import high_dimensional_sampling as hds
+import pypolychord as pc
+from anesthetic import NestedSamples
+from pypolychord.settings import PolyChordSettings
+
+
+class PolyChord(hds.Procedure):
+    def __init__(self):
+        self.store_parameters = []
+        self._is_finished = False
+        self._reset = False
+
+    def __call__(self, function):
+        ranges = function.get_ranges()
+
+        def prior(cube):
+            return cube * ranges[:, 1] + (1-cube) * ranges[:, 0]
+
+        def likelihood(theta):
+            return function([theta])[0, 0], []
+
+        def dumper(live, dead, logweights, logZ, logZerr):
+            pass
+
+        nDims = function.get_dimensionality()
+        nDerived = 0
+        settings = PolyChordSettings(nDims, nDerived)
+        settings.base_dir = '.polychord_chains'
+        settings.file_root = function.name
+        settings.posteriors = False
+        settings.equals = False
+        settings.write_live = False
+        settings.write_prior = False
+        settings.write_paramnames = False
+        settings.write_stats = True
+        settings.write_resume = True
+        settings.read_resume = True
+
+        if self._reset:
+            self._reset = False
+            settings.read_resume = False
+
+        if not self.is_finished():
+            pc.run_polychord(likelihood, nDims, nDerived,
+                             settings, prior, dumper)
+
+        self._is_finished = True
+        root = os.path.join(settings.base_dir, settings.file_root)
+        samples = NestedSamples(root=root)
+        x = samples.iloc[:, :nDims].to_numpy()
+        y = np.array([samples.logL.to_numpy()]).T
+
+        return (x, y)
+
+    def reset(self):
+        self._is_finished = False
+        self._reset = True
+
+    def is_finished(self):
+        return self._is_finished
+
+    def check_testfunction(self, function):
+        return True

manual/01_creating_a_procedure.md

@@ -29,7 +29,7 @@ This class should implement at least the following four methods:
 This is Initialisation method in which anything can be done. At the very least
 it should define the `self.store_parameters` property as a list, containing
 the names of the parameters of the Procedure. These parameters will then 
-autometically be stored when an experiment is run. If no procedure parameters
+automatically be stored when an experiment is run. If no procedure parameters
 exist, `self.store_parameters` should be an empty list.
 
 Your implementation of this initialisation method is allowed to include extra
@@ -51,7 +51,7 @@ should be a numpy array of shape `(nDatapoints, nOutputVariables)`. Note that
 `nOutputVariables` will often be 1.
 
 Although it is not required, it is advisable to implement the `__call__` method
-in such a way that it returns a single data point at a time. Whether of not
+in such a way that it returns a single data point at a time. Whether or not
 this is possible depends of course on the implemented procedure, but following
 this guideline makes sure that the experiment is stopped as soon as possible
 and no further (possibly expensive) iterations are run. This makes comparison
@@ -64,7 +64,7 @@ stopped, or `False` if it should continue. If `True`, this overrides the
 `finish_line` condition set in [the Experiment](04_running_an_experiment.md).
 
 ### `reset(self)`
-This method should reset all internal variables conserning previous runs. For
+This method should reset all internal variables concerning previous runs. For
 instance: if your implementation keeps track of the number of points sampled
 in a property called `n_sampled`, the `reset()` method should set `n_sampled`
 to 0. It is called every time a new testfunction is provided to the Procedure.
@@ -78,4 +78,4 @@ In the [examples](../examples) folder two example procedures are provided: one
 for [rejection sampling](../examples/rejection_sampling.py), a procedure for
 posterior sampling, the other for 
 [random optimisation](../examples/random_optimisation.py), a procedure for 
-optimum finding.
+optimum finding.

manual/02_using_testfunctions.md

@@ -70,7 +70,7 @@ defined.
 - **_evaluate(self, x)**: method that evaluates the test function at coordinate
 `x`. Be aware that `x` is a numpy array of shape `(nDatapoints, nVariables)`
 and can contain multiple rows (i.e. data points). This method should return
-the function values of your test fucntion in the form of a numpy array of shape
+the function values of your test function in the form of a numpy array of shape
 `(nDatapoints, nOutputVariables)`.
 - **_derivative(self, x)**: same as `_evaluate`, but now for the derivative of
 the testfunction. If no derivative is implemented, this method should raise
@@ -80,4 +80,4 @@ a `functions.NoDerivativeError`.
 Have a look at
 [Using and Looping over TestFunctions](03_selecting_and_looping_over_testfunctions.md)
 for more information on how to automatically select and loop over test
-functions.
+functions.