Emulated termination criteria

README.rst

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

anesthetic/_version.py

@@ -1 +1 @@
-__version__ = '2.8.3'
+__version__ = '2.9.0'

anesthetic/read/csv.py

@@ -5,7 +5,7 @@
 
 
 def read_csv(filename, *args, **kwargs):
-    """Read a CSV file into a :class:`Samples` object."""
+    """Read a CSV file into a :class:`anesthetic.samples.Samples` object."""
     filename = Path(filename)
     kwargs['label'] = kwargs.get('label', filename.stem)
     wldf = wl_read_csv(filename.with_suffix('.csv'))

anesthetic/samples.py

@@ -18,6 +18,7 @@
 from anesthetic.plot import (make_1d_axes, make_2d_axes,
                              AxesSeries, AxesDataFrame)
 from anesthetic.utils import adjust_docstrings
+from anesthetic import termination
 
 
 class Samples(WeightedLabelledDataFrame):
@@ -1193,6 +1194,26 @@ def truncate(self, logL=None):
         index = np.concatenate([dead_points.index, live_points.index])
         return self.loc[index].recompute()
 
+    def terminated(self, criterion='logZ', *args, **kwargs):
+        """Check if a simulated run has terminated.
+
+        Parameters
+        ----------
+        criterion : str, optional
+            The termination criterion to choose.
+            This should be one of
+            {'logZ', 'D_KL', 'logX', 'ndead', 'logL'}.
+            Default is logZ.
+            See the documentation for the specific criterion in
+            ``anesthetic.termination`` for more details.
+
+        The remaining arguments are passed to the termination criterion.
+        """
+        try:
+            return termination.__dict__[criterion](self, *args, **kwargs)
+        except KeyError:
+            raise KeyError("Unknown termination criterion: %s" % criterion)
+
     def posterior_points(self, beta=1):
         """Get equally weighted posterior points at temperature beta."""
         return self.set_beta(beta).compress('equal')

anesthetic/termination.py

@@ -0,0 +1,41 @@
+"""Nested sampling termination criteria."""
+from scipy.special import logsumexp
+import numpy as np
+
+
+def logZ(samples, eps=1e-3, nsamples=None, beta=None):
+    """Terminate when evidence in the live points is fraction eps of total."""
+    i_live = samples.live_points().index
+    logw = samples.logw(nsamples, beta)
+    logZ_live = logsumexp(logw[i_live])
+    logZ = logsumexp(logw)
+    return logZ_live - logZ < np.log(eps)
+
+
+def D_KL(samples, eps=1e-3, nsamples=None, beta=None):
+    """Terminate when D_KL in the live points is fraction eps of total."""
+    i_live = samples.live_points().index
+    logw = samples.logw(nsamples, beta)
+    logZ = logsumexp(logw)
+    betalogL = samples._betalogL(beta)
+    S = (logw*0).add(betalogL, axis=0) - logZ
+    w = np.exp(logw-logZ)
+    D_KL_live = (S*w).loc[i_live].sum()
+    D_KL = (S*w).sum()
+    return D_KL_live / D_KL < eps
+
+
+def logX(samples, max_logX, nsamples=None):
+    """Terminate when the log-volume of the live points reaches a threshold."""
+    i_live = samples.live_points().index
+    return samples.logX(nsamples)[i_live[0]] < max_logX
+
+
+def ndead(samples, max_ndead):
+    """Terminate if the number of dead points exceeds a maximum."""
+    return len(samples.dead_points()) > max_ndead
+
+
+def logL(samples, max_logL):
+    """Terminate if the lowest live likelihood exceeds a threshold."""
+    return samples.contour() > max_logL

docs/source/anesthetic.read.rst

@@ -25,6 +25,15 @@ anesthetic.read.cobaya module
    :show-inheritance:
 
 
+anesthetic.read.csv module
+--------------------------
+
+.. automodule:: anesthetic.read.csv
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
 anesthetic.read.getdist module
 ------------------------------
 

docs/source/anesthetic.rst

@@ -83,6 +83,15 @@ anesthetic.scripts module
    :show-inheritance:
 
 
+anesthetic.termination module
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. automodule:: anesthetic.termination
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
 anesthetic.testing module
 ~~~~~~~~~~~~~~~~~~~~~~~~~
 

tests/test_samples.py

@@ -1180,6 +1180,24 @@ def test_truncate(cut):
         assert_array_equal(pc, truncated_run)
 
 
+@pytest.mark.parametrize("criterion,args,terminated",
+                         [("logZ", (1e-3,), False),
+                          ("logZ", (1e-1,), True),
+                          ("D_KL", (1e-1,), True),
+                          ("D_KL", (1e-3,), False),
+                          ("logX", (-10.0), True),
+                          ("logX", (-15.0), False),
+                          ("ndead", (1000), True),
+                          ("ndead", (2000), False),
+                          ("logL", (5.0,), True),
+                          ("logL", (6.0,), False),
+                          ])
+def test_terminated(criterion, args, terminated):
+    np.random.seed(4)
+    pc = read_chains("./tests/example_data/pc")
+    assert pc.terminated(criterion, *args) is terminated
+
+
 def test_hist_range_1d():
     """Test to provide a solution to #89"""
     np.random.seed(3)