add ultranest import

README.rst

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

anesthetic/_version.py

@@ -1 +1 @@
-__version__ = '2.0.0'
+__version__ = '2.1.0'

anesthetic/plot.py

@@ -355,7 +355,8 @@ def _set_labels(axes, labels, **kwargs):
         all_params = list(axes.columns) + list(axes.index)
         if labels is None:
             labels = {}
-        labels = {p: labels[p] if p in labels else p for p in all_params}
+        labels = {p: labels[p] if isinstance(labels, dict) and p in labels
+                  else p for p in all_params}
 
         for y, axes_row in axes.iterrows():
             if axes_row.size:

anesthetic/read/chain.py

@@ -3,6 +3,7 @@
 from anesthetic.read.getdist import read_getdist
 from anesthetic.read.cobaya import read_cobaya
 from anesthetic.read.multinest import read_multinest
+from anesthetic.read.ultranest import read_ultranest
 
 
 def read_chains(root, *args, **kwargs):
@@ -45,7 +46,11 @@ def read_chains(root, *args, **kwargs):
             "anesthetic.html#anesthetic.samples.MCMCSamples.remove_burn_in"
         )
     errors = []
-    for read in [read_polychord, read_multinest, read_cobaya, read_getdist]:
+    readers = [
+        read_polychord, read_multinest, read_cobaya,
+        read_ultranest, read_getdist
+    ]
+    for read in readers:
         try:
             samples = read(root, *args, **kwargs)
             samples.root = root

anesthetic/read/ultranest.py

@@ -0,0 +1,33 @@
+"""Read NestedSamples from ultranest results."""
+import os
+import json
+from anesthetic.samples import NestedSamples
+try:
+    import h5py
+except ImportError:
+    pass
+
+
+def read_ultranest(root, *args, **kwargs):
+    """Read <root>ev.dat and <root>phys_live.points in MultiNest format."""
+    with open(os.path.join(root, 'info', 'results.json')) as infofile:
+        labels = json.load(infofile)['paramnames']
+    num_params = len(labels)
+
+    filepath = os.path.join(root, 'results', 'points.hdf5')
+    with h5py.File(filepath, 'r') as fileobj:
+        points = fileobj['points']
+        _, ncols = points.shape
+        x_dim = ncols - 3 - num_params
+        logL_birth = points[:, 0]
+        logL = points[:, 1]
+        # u_samples = points[:,3:x_dim]
+        samples = points[:, 3 + x_dim:3 + x_dim + num_params]
+
+    kwargs['label'] = kwargs.get('label', os.path.basename(root))
+    columns = kwargs.pop('columns', labels)
+    labels = kwargs.pop('labels', labels)
+    data = samples
+
+    return NestedSamples(data=data, logL=logL, logL_birth=logL_birth,
+                         columns=columns, labels=labels, *args, **kwargs)

docs/source/reading_writing.rst

@@ -5,16 +5,17 @@ Reading and writing
 
 .. _reading chains:
 
-Reading chain files from PolyChord, MultiNest, CosmoMC, or Cobaya
-=================================================================
+Reading chain files from PolyChord, MultiNest, UltraNest, CosmoMC, or Cobaya
+============================================================================
 
-`MultiNest <https://github.com/farhanferoz/MultiNest>`_
+If you have finished nested sampling or MCMC runs from one of:
+
+* `PolyChord <https://polychord.io>`_: https://github.com/PolyChord/PolyChordLite
+* MultiNest: https://github.com/farhanferoz/MultiNest
+* `UltraNest <https://johannesbuchner.github.io/UltraNest/index.html>`_: https://github.com/JohannesBuchner/UltraNest
+* `CosmoMC <https://cosmologist.info/cosmomc/readme.html>`_: https://github.com/cmbant/CosmoMC
+* `Cobaya <https://cobaya.readthedocs.io>`_: https://github.com/CobayaSampler/cobaya
 
-If you have finished nested sampling or MCMC runs from
-`PolyChord <https://github.com/PolyChord/PolyChordLite>`_,
-`MultiNest <https://github.com/farhanferoz/MultiNest>`_,
-`CosmoMC <https://github.com/cmbant/CosmoMC>`_, or
-`Cobaya <https://github.com/CobayaSampler/cobaya>`_,
 then you should be able to read in the chain files directly, by passing the
 ``root`` to the :func:`anesthetic.read.chain.read_chains` function.
 
@@ -29,6 +30,14 @@ out the examples listed here.
       from anesthetic import read_chains
       samples = read_chains("anesthetic/tests/example_data/pc")
 
+* UltraNest samples, which will be an instance of the
+  :class:`anesthetic.samples.NestedSamples` class:
+
+  ::
+
+      from anesthetic import read_chains
+      samples = read_chains("anesthetic/tests/example_data/un")
+
 * Cobaya samples, which will be an instance of the
   :class:`anesthetic.samples.MCMCSamples` class:
 

pyproject.toml

@@ -14,6 +14,7 @@ authors = [
     { name="Lukas Hergt", email="[email protected]" },
     { name="Adam Ormondroyd", email="[email protected]" },
     { name="Harry Bevins", email="[email protected]" },
+    { name="Johannes Buchner", email="[email protected]" },
     { name="Ethan Carragher", email="[email protected]" },
     { name="Andrew Fowlie", email="[email protected]" },
     { name="Thomas Gessey-Jones", email="[email protected]"},
@@ -59,11 +60,12 @@ classifiers = [
 [project.optional-dependencies]
 docs = ["sphinx", "sphinx_rtd_theme", "numpydoc"]
 test = ["pytest", "pytest-cov", "flake8", "pydocstyle", "packaging", "pre-commit"]
+ultranest = ["h5py"]
 astropy = ["astropy"]
 fastkde = ["fastkde"]
 getdist = ["getdist"]
 hdf5 = ["tables"]
-all = ["astropy", "fastkde", "getdist", "tables"]
+all = ["h5py", "astropy", "fastkde", "getdist", "tables"]
 
 [project.scripts]
 anesthetic = "anesthetic.scripts:gui"

tests/example_data/un/info/results.json

@@ -0,0 +1,77 @@
+{
+    "niter": 16650,
+    "logz": -12.414914594805953,
+    "logzerr": 0.7649516021328628,
+    "logz_bs": -12.378097297348035,
+    "logz_single": -12.414914594805953,
+    "logzerr_tail": 0.405716490463611,
+    "logzerr_bs": 0.6484944741256662,
+    "ess": 988.8766189106832,
+    "H": 38.639675527033575,
+    "Herr": 0.22373181701346842,
+    "posterior": {
+        "mean": [
+            0.5025783329483371,
+            0.7370934893843076,
+            0.9206375311748793,
+            0.9987450408249887
+        ],
+        "stdev": [
+            0.09929557150158441,
+            0.0021313377831688133,
+            4.6033000067476255e-05,
+            1.0182458288564012e-06
+        ],
+        "median": [
+            0.506654697190454,
+            0.7370548529653416,
+            0.920636543102102,
+            0.9987450545520016
+        ],
+        "errlo": [
+            0.4016698354319299,
+            0.7350305173203235,
+            0.9205913816528906,
+            0.9987439938546689
+        ],
+        "errup": [
+            0.601282668275072,
+            0.7393071000969336,
+            0.9206855255078583,
+            0.9987460211457329
+        ],
+        "information_gain_bits": [
+            3.7847253159407837,
+            4.14647919764959,
+            4.14647919764959,
+            4.14647919764959
+        ]
+    },
+    "maximum_likelihood": {
+        "logl": 28.526361154649237,
+        "point": [
+            0.4912994020358319,
+            0.7373000808358334,
+            0.9206344043582959,
+            0.9987452206617764
+        ],
+        "point_untransformed": [
+            0.5245649701017916,
+            0.5368650040417917,
+            0.5460317202179148,
+            0.5499372610330888
+        ]
+    },
+    "ncall": 39885,
+    "paramnames": [
+        "a",
+        "b",
+        "c",
+        "d"
+    ],
+    "logzerr_single": 0.31080410038734035,
+    "insertion_order_MWW_test": {
+        "independent_iterations": Infinity,
+        "converged": true
+    }
+}

tests/example_data/un_gen_data.py

@@ -0,0 +1,31 @@
+import numpy as np
+from ultranest import ReactiveNestedSampler
+
+ndim = 4
+sigma = np.logspace(-1, -6, ndim)
+width = 1 - 5 * sigma
+width[width < 1e-20] = 1e-20
+centers = (np.sin(np.arange(ndim) / 2.) * width + 1.) / 2.
+
+
+def loglike(theta):
+    """compute log-likelihood."""
+    like = - 0.5 * (((theta - centers) / sigma)**2).sum(axis=1) \
+           - 0.5 * np.log(2 * np.pi * sigma**2).sum()
+    return like
+
+
+def transform(x):
+    """transform according to prior."""
+    return x * 20 - 10
+
+
+paramnames = ['a', 'b', 'c', 'd']
+
+sampler = ReactiveNestedSampler(
+    paramnames, loglike, transform=transform,
+    log_dir='un', resume=True, vectorized=True)
+
+sampler.run(frac_remain=0.5, min_num_live_points=400)
+sampler.print_results()
+sampler.plot()

tests/test_gui.py

@@ -1,16 +1,25 @@
+import sys
 import anesthetic.examples._matplotlib_agg  # noqa: F401
 from anesthetic import read_chains
 import pytest
 import pandas._testing as tm
+try:
+    import h5py  # noqa: F401
+except ImportError:
+    pass
 
 
 @pytest.fixture(autouse=True)
 def close_figures_on_teardown():
     tm.close()
 
 
-def test_gui():
-    samples = read_chains('./tests/example_data/pc')
+@pytest.mark.parametrize('root', ["./tests/example_data/pc",
+                                  "./tests/example_data/un"])
+def test_gui(root):
+    if 'un' in root and 'h5py' not in sys.modules:
+        pytest.skip("`h5py` not in sys.modules, but needed for ultranest.")
+    samples = read_chains(root)
     plotter = samples.gui()
 
     # Type buttons
@@ -30,8 +39,11 @@ def test_gui():
     assert len(plotter.triangle.ax) == 4
 
     # Sliders
+    old = plotter.evolution()
     plotter.evolution.slider.set_val(5)
-    assert plotter.evolution() == 627
+    assert plotter.evolution() != old
+    plotter.evolution.slider.set_val(0)
+    assert plotter.evolution() == old
     plotter.type.buttons.set_active(1)
 
     plotter.beta.slider.set_val(0)
@@ -50,15 +62,25 @@ def test_gui():
     plotter.reset.button.on_clicked(plotter.reset_range(None))
 
 
-def test_gui_params():
-    plotter = read_chains('./tests/example_data/pc').gui()
-    assert len(plotter.param_choice.buttons.labels) == 8
+@pytest.mark.parametrize('root', ["./tests/example_data/pc",
+                                  "./tests/example_data/un"])
+def test_gui_params(root):
+    if 'un' in root and 'h5py' not in sys.modules:
+        pytest.skip("`h5py` not in sys.modules, but needed for ultranest.")
+    samples = read_chains(root)
+    params = samples.columns.get_level_values(0).to_list()
+    plotter = samples.gui()
+    assert len(plotter.param_choice.buttons.labels) == len(params)
 
-    plotter = read_chains('./tests/example_data/pc').gui(params=['x0', 'x1'])
+    plotter = samples.gui(params=params[:2])
     assert len(plotter.param_choice.buttons.labels) == 2
 
 
-def test_slider_reset_range():
-    plotter = read_chains('./tests/example_data/pc').gui()
+@pytest.mark.parametrize('root', ["./tests/example_data/pc",
+                                  "./tests/example_data/un"])
+def test_slider_reset_range(root):
+    if 'un' in root and 'h5py' not in sys.modules:
+        pytest.skip("`h5py` not in sys.modules, but needed for ultranest.")
+    plotter = read_chains(root).gui()
     plotter.evolution.reset_range(-3, 3)
     assert plotter.evolution.ax.get_xlim() == (-3.0, 3.0)

tests/test_reader.py

@@ -11,8 +11,13 @@
 from anesthetic.read.getdist import read_getdist
 from anesthetic.read.cobaya import read_cobaya
 from anesthetic.read.multinest import read_multinest
+from anesthetic.read.ultranest import read_ultranest
 import pandas._testing as tm
 from anesthetic.read.hdf import HDFStore, read_hdf
+try:
+    import h5py  # noqa: F401
+except ImportError:
+    pass
 try:
     import getdist
 except ImportError:
@@ -184,6 +189,28 @@ def test_read_multinest():
     ns.plot_1d(['x0', 'x1', 'x2', 'x3'])
 
 
+@pytest.mark.xfail('h5py' not in sys.modules,
+                   raises=NameError,
+                   reason="ultranest reading requires `h5py` package")
+def test_read_ultranest():
+    np.random.seed(3)
+    ns = read_ultranest('./tests/example_data/un')
+    params = ['a', 'b', 'c', 'd', 'logL', 'logL_birth', 'nlive']
+    assert_array_equal(ns.drop_labels().columns, params)
+    labels = ['a',
+              'b',
+              'c',
+              'd',
+              r'$\ln\mathcal{L}$',
+              r'$\ln\mathcal{L}_\mathrm{birth}$',
+              r'$n_\mathrm{live}$']
+    assert_array_equal(ns.get_labels(), labels)
+
+    assert isinstance(ns, NestedSamples)
+    ns.plot_2d(['a', 'b', 'c', 'd'])
+    ns.plot_1d(['a', 'b', 'c', 'd'])
+
+
 def test_read_polychord():
     np.random.seed(3)
     ns = read_polychord('./tests/example_data/pc')