handley-lab · williamjameshandley · Aug 19, 2020 · Aug 18, 2020 · Aug 18, 2020 · Aug 18, 2020
diff --git a/README.rst b/README.rst
@@ -3,7 +3,7 @@ anesthetic: nested sampling visualisation
 =========================================
 :anesthetic: nested sampling visualisation
 :Author: Will Handley and Lukas Hergt
-:Version: 2.0.0-beta.1
+:Version: 2.0.0-beta.2
 :Homepage: https://github.com/williamjameshandley/anesthetic
 :Documentation: http://anesthetic.readthedocs.io/
 

diff --git a/anesthetic/gui/plot.py b/anesthetic/gui/plot.py
@@ -245,8 +245,7 @@ def points(self, label):
 
     def update(self, _):
         """Update all the plots upon slider changes."""
-        with np.errstate(divide='ignore'):
-            logX = np.log(self.samples.nlive / (self.samples.nlive+1)).cumsum()
+        logX = np.log(self.samples.nlive / (self.samples.nlive+1)).cumsum()
         kT = self.temperature()
         LX = self.samples.logL/kT + logX
         LX = np.exp(LX-LX.max())

diff --git a/anesthetic/samples.py b/anesthetic/samples.py
@@ -590,14 +590,13 @@ def dlogX(self, nsamples=None):
             distribution. (Default: None)
 
         """
-        with np.errstate(divide='ignore'):
-            if np.ndim(nsamples) > 0:
-                return nsamples
-            elif nsamples is None:
-                t = np.log(self.nlive/(self.nlive+1)).to_frame()
-            else:
-                r = np.log(np.random.rand(len(self), nsamples))
-                t = pandas.DataFrame(r, self.index).divide(self.nlive, axis=0)
+        if np.ndim(nsamples) > 0:
+            return nsamples
+        elif nsamples is None:
+            t = np.log(self.nlive/(self.nlive+1)).to_frame()
+        else:
+            r = np.log(np.random.rand(len(self), nsamples))
+            t = pandas.DataFrame(r, self.index).divide(self.nlive, axis=0)
 
         logX = t.cumsum()
         logXp = logX.shift(1, fill_value=0)

diff --git a/anesthetic/utils.py b/anesthetic/utils.py
@@ -182,13 +182,13 @@ def compute_nlive(death, birth):
     nlive: np.array
         number of live points at each contour
     """
-    birth_index = death.searchsorted(birth)
+    birth_index = death.searchsorted(birth, side='right')
     births = pandas.Series(+1, index=birth_index).sort_index()
-    index = np.arange(death.size)
+    index = np.arange(death.size)+1
     deaths = pandas.Series(-1, index=index)
     nlive = pandas.concat([births, deaths]).sort_index()
     nlive = nlive.groupby(nlive.index).sum().cumsum()
-    return nlive.values
+    return nlive.values[:-1]
 
 
 def compute_insertion_indexes(death, birth):

diff --git a/tests/test_gui.py b/tests/test_gui.py
@@ -6,10 +6,10 @@ def test_gui():
     plotter = NestedSamples(root='./tests/example_data/pc').gui()
 
     # Type buttons
-    plotter.type.buttons.set_active(0)
-    assert(plotter.type() == 'live')
     plotter.type.buttons.set_active(1)
     assert(plotter.type() == 'posterior')
+    plotter.type.buttons.set_active(0)
+    assert(plotter.type() == 'live')
 
     # Parameter choice buttons
     plotter.param_choice.buttons.set_active(1)
@@ -25,11 +25,17 @@ def test_gui():
     plotter.evolution.slider.set_val(100)
     assert(plotter.evolution() == 100)
     plotter.type.buttons.set_active(1)
+
     plotter.temperature.slider.set_val(0)
     assert(plotter.temperature() == 1)
-
     plotter.temperature.slider.set_val(1)
     assert(plotter.temperature() == 10)
     plotter.temperature.slider.set_val(2)
     assert(plotter.temperature() == 100)
     plotter.type.buttons.set_active(0)
+
+    # Reload button
+    plotter.reload.button.observers[0](None)
+
+    # Reset button
+    plotter.reset.button.observers[0](None)
diff --git a/tests/test_samples.py b/tests/test_samples.py
@@ -12,6 +12,7 @@
                            assert_array_less)
 from matplotlib.colors import to_hex
 from scipy.stats import ks_2samp, kstest
+from wedding_cake import WeddingCake
 try:
     import montepython  # noqa: F401
 except ImportError:
@@ -218,6 +219,13 @@ def test_plot_2d_types():
     fig, axes = ns.plot_2d(params, types={'lower': 'kde', 'diagonal': 'kde',
                                           'upper': 'scatter'})
     assert((~axes.isnull()).sum().sum() == 12)
+
+    with pytest.raises(NotImplementedError):
+        fig, axes = ns.plot_2d(params, types={'lower': 'not a plot type'})
+
+    with pytest.raises(NotImplementedError):
+        fig, axes = ns.plot_2d(params, types={'diagonal': 'not a plot type'})
+
     plt.close("all")
 
 
@@ -539,9 +547,9 @@ def test_live_points():
     np.random.seed(4)
     pc = NestedSamples(root="./tests/example_data/pc")
 
-    for i, logL in pc.logL.iteritems():
+    for i, logL in pc.logL.iloc[:-1].iteritems():
         live_points = pc.live_points(logL)
-        assert len(live_points) == int(pc.nlive[i])
+        assert len(live_points) == int(pc.nlive[i[0]+1])
 
         live_points_from_int = pc.live_points(i[0])
         assert_array_equal(live_points_from_int, live_points)
@@ -573,7 +581,7 @@ def test_limit_assignment():
     # limits for logL, weights, nlive
     assert ns.limits['logL'][0] == -777.0115456428716
     assert ns.limits['logL'][1] == 5.748335384373301
-    assert ns.limits['nlive'][0] == 0
+    assert ns.limits['nlive'][0] == 1
     assert ns.limits['nlive'][1] == 125
 
 
@@ -638,3 +646,15 @@ def test_posterior_points():
     ns = NestedSamples(root='./tests/example_data/pc')
     assert_array_equal(ns.posterior_points(), ns.posterior_points())
     assert_array_equal(ns.posterior_points(0.5), ns.posterior_points(0.5))
+
+
+def test_wedding_cake():
+    np.random.seed(3)
+    wc = WeddingCake(4, 0.5, 0.01)
+    nlive = 500
+    samples = wc.sample(nlive)
+    assert samples.nlive.iloc[0] == nlive
+    assert samples.nlive.iloc[-1] == 1
+    assert (samples.nlive <= nlive).all()
+    out = samples.logZ(100)
+    assert abs(out.mean()-wc.logZ()) < out.std()*3
diff --git a/tests/test_utils.py b/tests/test_utils.py
@@ -37,8 +37,8 @@ def test_compute_nlive():
     # Check the first half are constant
     assert_array_equal(nlives[:len(nlives)//2], nlive)
 
-    # Check no points at the end
-    assert(nlives[-1] == 0)
+    # Check one point at the end
+    assert(nlives[-1] == 1)
 
     # Check never more than nlive
     assert(nlives.max() <= nlive)

diff --git a/tests/wedding_cake.py b/tests/wedding_cake.py
@@ -0,0 +1,86 @@
+import numpy as np
+from scipy.special import logsumexp
+from anesthetic import NestedSamples
+
+
+class WeddingCake():
+    """Class for generating samples from a wedding cake 'likelihood'.
+
+    This is a likelihood with nested hypercuboidal plateau regions of constant
+    likelihood centered on 0.5, with geometrically decreasing volume by a
+    factor of alpha. The value of the likelihood in these plateau regions has a
+    gaussian profile with width sigma.
+
+    logL = - alpha^(2 floor(D*log_alpha(2|x-0.5|_infinity))/D) / (8 sigma^2)
+
+    Parameters
+    ----------
+    D: int
+        dimensionality (number of parameters) of the likelihood
+
+    alpha: float
+        volume compression between plateau regions
+
+    sigma: float
+        width of gaussian profile
+    """
+    def __init__(self, D=4, alpha=0.5, sigma=0.01):
+        self.D = D
+        self.alpha = alpha
+        self.sigma = sigma
+
+    def logZ(self):
+        """Numerically compute the true evidence."""
+        mean = np.sqrt(self.D/2.) * (np.log(4*self.D*self.sigma**2)-1
+                                     ) / np.log(self.alpha)
+        std = -np.sqrt(self.D/2.)/np.log(self.alpha)
+        i = np.arange(mean + std*10)
+
+        return logsumexp(-self.alpha**(2*i/self.D)/8/self.sigma**2
+                         + i*np.log(self.alpha) + np.log(1-self.alpha))
+
+    def i(self, x):
+        """Plateau number of a parameter point."""
+        r = np.max(abs(x-0.5), axis=-1)
+        return np.floor(self.D*np.log(2*r)/np.log(self.alpha))
+
+    def logL(self, x):
+        """Gaussian log-likelihood."""
+        ri = self.alpha**(self.i(x)/self.D)/2
+        return - ri**2/2/self.sigma**2
+
+    def sample(self, nlive=500):
+        """Generate samples from a perfect nested sampling run."""
+        points = np.zeros((0, self.D))
+        death_likes = np.zeros(0)
+        birth_likes = np.zeros(0)
+
+        live_points = np.random.rand(nlive, self.D)
+        live_likes = self.logL(live_points)
+        live_birth_likes = np.ones(nlive)*-np.inf
+
+        while True:
+            logL_ = live_likes.min()
+            j = live_likes == logL_
+
+            death_likes = np.concatenate([death_likes, live_likes[j]])
+            birth_likes = np.concatenate([birth_likes, live_birth_likes[j]])
+            points = np.concatenate([points, live_points[j]])
+
+            i_ = self.i(live_points[j][0])+1
+            r_ = self.alpha**(i_/self.D)/2
+            x_ = np.random.uniform(0.5-r_, 0.5+r_, size=(j.sum(), self.D))
+            live_birth_likes[j] = logL_
+            live_points[j] = x_
+            live_likes[j] = self.logL(x_)
+
+            samps = NestedSamples(points, logL=death_likes,
+                                  logL_birth=birth_likes)
+            if samps.live_points().weights.sum()/samps.weights.sum() < 0.001:
+                break
+
+        death_likes = np.concatenate([death_likes, live_likes])
+        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)