Merge branch 'master' into bump_version

.gitignore

@@ -13,3 +13,5 @@ plikHM_TTTEEE_lowl_lowE_lensing_NS/
 *~
 .pytest_cache/*
 .coverage
+
+.DS_Store

README.rst

@@ -2,7 +2,7 @@
 anesthetic: nested sampling post-processing
 ===========================================
 :Authors: Will Handley and Lukas Hergt
-:Version: 2.6.1
+:Version: 2.7.1
 :Homepage: https://github.com/handley-lab/anesthetic
 :Documentation: http://anesthetic.readthedocs.io/
 

anesthetic/__init__.py

@@ -9,7 +9,6 @@
 import pandas.plotting._misc
 from anesthetic._format import _DataFrameFormatter
 from anesthetic._version import __version__  # noqa: F401
-
 # TODO: remove this when conda pandas version catches up
 from packaging.version import parse
 assert parse(pandas.__version__) >= parse('2.0.0')

anesthetic/_version.py

@@ -1 +1 @@
-__version__ = '2.6.1'
+__version__ = '2.7.1'

anesthetic/examples/perfect_ns.py

@@ -5,7 +5,8 @@
 from anesthetic.samples import merge_nested_samples
 
 
-def gaussian(nlive, ndims, sigma=0.1, R=1, logLmin=-1e-2):
+def gaussian(nlive, ndims, sigma=0.1, R=1, logLmin=-1e-2, logLmax=0,
+             *args, **kwargs):
     """Perfect nested sampling run for a spherical Gaussian & prior.
 
     Up to normalisation this is identical to the example in John Skilling's
@@ -29,14 +30,20 @@ def gaussian(nlive, ndims, sigma=0.1, R=1, logLmin=-1e-2):
     logLmin : float
         loglikelihood at which to terminate
 
+    logLmax : float
+        maximum loglikelihood
+
+    The remaining arguments are passed to the
+    :class:`anesthetic.samples.NestedSamples` constructor.
+
     Returns
     -------
     samples : :class:`anesthetic.samples.NestedSamples`
         Nested sampling run
     """
 
     def logLike(x):
-        return -(x**2).sum(axis=-1)/2/sigma**2
+        return logLmax - (x**2).sum(axis=-1)/2/sigma**2
 
     def random_sphere(n):
         return random_ellipsoid(np.zeros(ndims), np.eye(ndims), n)
@@ -48,7 +55,8 @@ def random_sphere(n):
     while logL.min() < logLmin:
         points = r * random_sphere(nlive)
         logL = logLike(points)
-        samples.append(NestedSamples(points, logL=logL, logL_birth=logL_birth))
+        samples.append(NestedSamples(points, logL=logL, logL_birth=logL_birth,
+                                     *args, **kwargs))
         logL_birth = logL.copy()
         r = (points**2).sum(axis=-1, keepdims=True)**0.5
 
@@ -57,14 +65,14 @@ def random_sphere(n):
     return samples.loc[samples.logL_birth < logLend].recompute()
 
 
-def correlated_gaussian(nlive, mean, cov, bounds=None):
+def correlated_gaussian(nlive, mean, cov, bounds=None, logLmax=0,
+                        *args, **kwargs):
     """Perfect nested sampling run for a correlated gaussian likelihood.
 
-    This produces a perfect nested sampling run with a uniform prior over the
-    unit hypercube, with a likelihood gaussian in the parameters normalised so
-    that the evidence is unity. The algorithm proceeds by simultaneously
-    rejection sampling from the prior and sampling exactly and uniformly from
-    the known ellipsoidal contours.
+    This produces a perfect nested sampling run with a uniform prior over
+    the unit hypercube. The algorithm proceeds by simultaneously rejection
+    sampling from the prior and sampling exactly and uniformly from the
+    known ellipsoidal contours.
 
     This can produce perfect runs in up to around d~15 dimensions. Beyond
     this rejection sampling from a truncated gaussian in the early stage
@@ -85,6 +93,12 @@ def correlated_gaussian(nlive, mean, cov, bounds=None):
     bounds : 2d array-like, shape (ndims, 2)
         bounds of a gaussian, default ``[[0, 1]]*ndims``
 
+    logLmax : float
+        maximum loglikelihood
+
+    The remaining arguments are passed to the
+    :class:`anesthetic.samples.NestedSamples` constructor.
+
     Returns
     -------
     samples : :class:`anesthetic.samples.NestedSamples`
@@ -95,7 +109,7 @@ def correlated_gaussian(nlive, mean, cov, bounds=None):
     invcov = np.linalg.inv(cov)
 
     def logLike(x):
-        return -0.5 * ((x-mean) @ invcov * (x-mean)).sum(axis=-1)
+        return logLmax - 0.5 * ((x - mean) @ invcov * (x - mean)).sum(axis=-1)
 
     ndims = len(mean)
 
@@ -104,10 +118,9 @@ def logLike(x):
 
     bounds = np.array(bounds, dtype=float)
 
-    logLmax = logLike(mean)
-
     points = np.random.uniform(*bounds.T, (2*nlive, ndims))
-    samples = NestedSamples(points, logL=logLike(points), logL_birth=-np.inf)
+    samples = NestedSamples(points, logL=logLike(points), logL_birth=-np.inf,
+                            *args, **kwargs)
 
     while (1/samples.nlive.iloc[:-nlive]).sum() < samples.D_KL()*2:
         logLs = samples.logL.iloc[-nlive]
@@ -120,19 +133,21 @@ def logLike(x):
         points = np.random.uniform(*bounds.T, (nlive, ndims))
         logL = logLike(points)
         i = logL > logLs
-        samps_1 = NestedSamples(points[i], logL=logL[i], logL_birth=logLs)
+        samps_1 = NestedSamples(points[i], logL=logL[i], logL_birth=logLs,
+                                *args, **kwargs)
 
         # Ellipsoidal round
         points = random_ellipsoid(mean, cov*2*(logLmax - logLs), nlive)
         logL = logLike(points)
         i = ((points > bounds.T[0]) & (points < bounds.T[1])).all(axis=1)
-        samps_2 = NestedSamples(points[i], logL=logL[i], logL_birth=logLs)
+        samps_2 = NestedSamples(points[i], logL=logL[i], logL_birth=logLs,
+                                *args, **kwargs)
         samples = merge_nested_samples([samples, samps_1, samps_2])
 
     return samples
 
 
-def wedding_cake(nlive, ndims, sigma=0.01, alpha=0.5):
+def wedding_cake(nlive, ndims, sigma=0.01, alpha=0.5, *args, **kwargs):
     """Perfect nested sampling run for a wedding cake likelihood.
 
     This is a likelihood with nested hypercuboidal plateau regions of constant
@@ -157,6 +172,9 @@ def wedding_cake(nlive, ndims, sigma=0.01, alpha=0.5):
 
     alpha : float
         volume compression between plateau regions
+
+    The remaining arguments are passed to the
+    :class:`anesthetic.samples.NestedSamples` constructor.
     """
 
     def i(x):
@@ -192,7 +210,8 @@ def logL(x):
         live_points[j] = x_
         live_likes[j] = logL(x_)
 
-        samps = NestedSamples(points, logL=death_likes, logL_birth=birth_likes)
+        samps = NestedSamples(points, logL=death_likes, logL_birth=birth_likes,
+                              *args, **kwargs)
         weights = samps.get_weights()
         if weights[-nlive:].sum() < 0.001 * weights.sum():
             break
@@ -201,7 +220,8 @@ def logL(x):
     birth_likes = np.concatenate([birth_likes, live_birth_likes])
     points = np.concatenate([points, live_points])
 
-    return NestedSamples(points, logL=death_likes, logL_birth=birth_likes)
+    return NestedSamples(points, logL=death_likes, logL_birth=birth_likes,
+                         *args, **kwargs)
 
 
 def planck_gaussian(nlive=500):
@@ -244,14 +264,7 @@ def planck_gaussian(nlive=500):
         [8.00e-01, 1.20e+00]])
 
     logL_mean = -1400.35
-
-    samples = correlated_gaussian(nlive, mean, cov, bounds)
-
-    data = samples.iloc[:, :len(columns)].to_numpy()
-    logL = samples['logL'].to_numpy()
-    logL_birth = samples['logL_birth'].to_numpy()
-    logL_birth += logL_mean - samples.logL.mean()
-    logL += logL_mean - samples.logL.mean()
-    samples = NestedSamples(data=data, columns=columns, labels=labels,
-                            logL=logL, logL_birth=logL_birth)
-    return samples
+    d = len(mean)
+    logLmax = logL_mean + d/2
+    return correlated_gaussian(nlive, mean, cov, bounds, logLmax,
+                               columns=columns, labels=labels)

pyproject.toml

@@ -25,6 +25,7 @@ authors = [
     { name="Aleksandr Petrosyan", email="[email protected]" },
     { name="Liangliang Su", email="[email protected]"},
     { name="David Yallup", email="[email protected]" },
+    { name="Dily Ong", email="[email protected]" }
 ]
 description = "nested sampling post-processing"
 readme = "README.rst"

tests/test_examples.py

@@ -144,8 +144,25 @@ def test_planck():
     assert (planck[params] < bounds.T[1]).all().all()
     assert_allclose(planck[params].mean(), mean, atol=1e-3)
     assert_allclose(planck[params].cov(), cov, atol=1e-3)
-    assert_allclose(planck.logL.mean(), logL_mean, atol=1e-3)
+    assert_allclose(planck.logL.mean(), logL_mean,
+                    atol=3*planck.logL_P(12).std())
     assert_allclose(planck.logZ(), logZ, atol=3*planck.logZ(12).std())
     assert_allclose(planck.D_KL(), D_KL, atol=3*planck.D_KL(12).std())
     assert_allclose(planck.d_G(), len(params), atol=3*planck.d_G(12).std())
     assert_allclose(planck.logL_P(), logL_mean, atol=3*planck.logL_P(12).std())
+
+
+def test_logLmax():
+    nlive = 1000
+    mean = [0.1, 0.3, 0.5]
+    # cov = covariance matrix
+    cov = np.array([[.01, 0.009, 0], [0.009, .01, 0], [0, 0, 0.1]])*0.01
+    bounds = [[0, 1], [0, 1], [0, 1]]
+    logLmax = 10
+    # d = number of parameters
+    d = len(mean)
+    samples = correlated_gaussian(nlive, mean, cov, bounds=bounds,
+                                  logLmax=logLmax)
+    abs_err = samples.logL.std()/np.sqrt(samples.neff())
+    atol = abs_err*3
+    assert_allclose(samples.logL.mean(), logLmax-d/2, atol=atol)