Merge branch 'master' into post_2_deprecate

README.rst

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

anesthetic/_version.py

@@ -1 +1 @@
-__version__ = '2.2.0'
+__version__ = '2.4.0'

anesthetic/plot.py

@@ -25,7 +25,7 @@
 from anesthetic.utils import (sample_compression_1d, quantile,
                               triangular_sample_compression_2d,
                               iso_probability_contours,
-                              match_contour_to_contourf)
+                              match_contour_to_contourf, histogram_bin_edges)
 from anesthetic.boundary import cut_and_normalise_gaussian
 
 
@@ -948,6 +948,12 @@ def hist_plot_1d(ax, data, *args, **kwargs):
     xmax = quantile(data, q[-1], weights)
     range = kwargs.pop('range', (xmin, xmax))
 
+    if isinstance(bins, str) and bins in ['fd', 'scott', 'sqrt']:
+        bins = histogram_bin_edges(data,
+                                   weights=weights,
+                                   bins=bins,
+                                   beta=kwargs.pop('beta', 'equal'),
+                                   range=range)
     if isinstance(bins, str) and bins in ['knuth', 'freedman', 'blocks']:
         try:
             from astropy.visualization import hist

anesthetic/plotting/_matplotlib/hist.py

@@ -1,3 +1,4 @@
+from __future__ import annotations
 from matplotlib import pyplot as plt
 import numpy as np
 from pandas.plotting._matplotlib.hist import (HistPlot as _HistPlot,
@@ -21,7 +22,7 @@
     hist_plot_1d,
     quantile_plot_interval,
 )
-from anesthetic.utils import quantile
+from anesthetic.utils import quantile, histogram_bin_edges
 
 
 class HistPlot(_WeightedMPLPlot, _HistPlot):
@@ -139,19 +140,32 @@ class Hist1dPlot(HistPlot):
     def _kind(self) -> Literal["hist_1d"]:
         return "hist_1d"
 
-    def _calculate_bins(self, data):
-        if "range" not in self.kwds:
+    def __init__(
+            self,
+            data,
+            bins: str | int | np.ndarray | list[np.ndarray] = 'fd',
+            bottom: int | np.ndarray = 0,
+            **kwargs,
+    ) -> None:
+        super().__init__(data, bins=bins, bottom=bottom, **kwargs)
+
+    def _args_adjust(self) -> None:
+        if 'range' not in self.kwds:
             q = self.kwds.get('q', 5)
             q = quantile_plot_interval(q=q)
-            weights = self.kwds.get("weights", None)
-            xmin = quantile(data, q[0], weights)
-            xmax = quantile(data, q[-1], weights)
-            self.kwds["range"] = (xmin, xmax)
-            result = super()._calculate_bins(data)
-            self.kwds.pop("range")
-        else:
-            result = super()._calculate_bins(data)
-        return result
+            weights = self.kwds.get('weights', None)
+            xmin = quantile(self.data, q[0], weights)
+            xmax = quantile(self.data, q[-1], weights)
+            self.kwds['range'] = (xmin, xmax)
+        if isinstance(self.bins, str) and self.bins in ['fd', 'scott', 'sqrt']:
+            self.bins = histogram_bin_edges(
+                self.data,
+                weights=self.kwds.get('weights', None),
+                bins=self.bins,
+                beta=self.kwds.pop('beta', 'equal'),
+                range=self.kwds.get('range', None)
+            )
+        super()._args_adjust()
 
     @classmethod
     def _plot(

anesthetic/read/chain.py

@@ -4,6 +4,7 @@
 from anesthetic.read.cobaya import read_cobaya
 from anesthetic.read.multinest import read_multinest
 from anesthetic.read.ultranest import read_ultranest
+from anesthetic.read.nestedfit import read_nestedfit
 
 
 def read_chains(root, *args, **kwargs):
@@ -14,6 +15,7 @@ def read_chains(root, *args, **kwargs):
         * `PolyChord <https://github.com/PolyChord/PolyChordLite>`_,
         * `MultiNest <https://github.com/farhanferoz/MultiNest>`_,
         * `UltraNest <https://github.com/JohannesBuchner/UltraNest>`_,
+        * `Nested_fit <https://github.com/martinit18/Nested_Fit>`_,
         * `CosmoMC <https://github.com/cmbant/CosmoMC>`_,
         * `Cobaya <https://github.com/CobayaSampler/cobaya>`_,
         * or anything `GetDist <https://github.com/cmbant/getdist>`_
@@ -49,7 +51,7 @@ def read_chains(root, *args, **kwargs):
     errors = []
     readers = [
         read_polychord, read_multinest, read_cobaya,
-        read_ultranest, read_getdist
+        read_ultranest, read_nestedfit, read_getdist
     ]
     for read in readers:
         try:

anesthetic/read/nestedfit.py

@@ -0,0 +1,35 @@
+"""Read NestedSamples from Nested_Fit chains."""
+import os
+import numpy as np
+from anesthetic.read.getdist import read_paramnames
+from anesthetic.samples import NestedSamples
+
+
+def read_nestedfit(root, *args, **kwargs):
+    """Read Nested_Fit chain files.
+
+    Parameters
+    ----------
+    root : str
+        root specify the directory only, no specific roots,
+        The files read files are ``nf_output_points.txt``
+        and ``nf_output_diag.txt``.
+
+    """
+    dead_file = os.path.join(root, 'nf_output_points.txt')
+    birth_file = os.path.join(root, 'nf_output_diag.dat')
+    data_dead = np.loadtxt(dead_file)
+    data_birth = np.loadtxt(birth_file)
+    weight, logL, data = np.split(data_dead, [1, 2], axis=1)
+    logL_birth = data_birth[:, 0]
+    root_getdist = os.path.join(root, 'nf_output_points')
+    columns, labels = read_paramnames(root_getdist)
+    # No specific labeling is implemented in nested_fit
+    labels = columns
+    columns = kwargs.pop('columns', columns)
+    labels = kwargs.pop('labels', labels)
+    kwargs['label'] = kwargs.get('label', os.path.basename(root))
+
+    return NestedSamples(data=data, columns=columns,
+                         logL=logL, logL_birth=logL_birth,
+                         labels=labels, *args, **kwargs)

anesthetic/utils.py

@@ -219,6 +219,72 @@ def histogram(a, **kwargs):
     return xpath, ypath
 
 
+def histogram_bin_edges(samples, weights, bins='fd', range=None, beta='equal'):
+    """Compute a good number of bins dynamically from weighted samples.
+
+    Parameters
+    ----------
+    samples : array_like
+        Input data.
+
+    weights : array-like
+        Array of sample weights.
+
+    bins : str, default='fd'
+        String defining the rule used to automatically compute a good number
+        of bins for the weighted samples:
+
+            * 'fd'    : Freedman--Diaconis rule (modified for weighted data)
+            * 'scott' : Scott's rule (modified for weighted data)
+            * 'sqrt'  : Square root estimator (modified for weighted data)
+
+    range : (float, float), optional
+        The lower and upper range of the bins. If not provided, range is simply
+        ``(a.min(), a.max())``. Values outside the range are ignored. The first
+        element of the range must be less than or equal to the second.
+
+    beta : float, default='equal'
+        The value of beta>0 used to calculate the number of effective
+        samples via :func:`neff`.
+
+    Returns
+    -------
+    bin_edges : array of dtype float
+        The edges to pass to :func:`numpy.histogram`.
+
+    """
+    if weights is None:
+        weights = np.ones_like(samples)
+    if range is None:
+        minimum = np.min(samples)
+        maximum = np.max(samples)
+        data_range = maximum - minimum
+    else:
+        minimum = range[0]
+        maximum = range[1]
+        data_range = maximum - minimum
+        data_mask = (samples >= minimum) & (samples <= maximum)
+        samples = samples[data_mask]
+        weights = weights[data_mask]
+    w = weights / np.sum(weights)
+    N_eff = neff(w=w, beta=beta)
+    if bins == 'fd':  # Freedman--Diaconis rule
+        q1, q3 = quantile(samples, [1/4, 3/4], w=w)
+        bin_width = 2 * (q3 - q1) * N_eff**(-1/3)
+    elif bins == 'scott':  # Scott's rule
+        weighted_mean = np.average(samples, weights=w)
+        weighted_var = np.average((samples - weighted_mean)**2, weights=w)
+        weighted_std = np.sqrt(weighted_var)
+        bin_width = (24 * np.pi**0.5 / N_eff)**(1/3) * weighted_std
+    elif bins == 'sqrt':  # Square root estimator
+        samples_eff, _ = sample_compression_1d(samples, w=w, ncompress=N_eff)
+        data_range_eff = np.max(samples_eff) - np.min(samples_eff)
+        bin_width = data_range_eff / np.sqrt(N_eff)
+    nbins = int(np.ceil(data_range / bin_width))
+    bin_edges = np.linspace(minimum, maximum, nbins+1)
+    return bin_edges
+
+
 def compute_nlive(death, birth):
     """Compute number of live points from birth and death contours.
 

docs/source/reading_writing.rst

@@ -5,14 +5,15 @@ Reading and writing
 
 .. _reading chains:
 
-Reading chain files from PolyChord, MultiNest, UltraNest, CosmoMC, or Cobaya
-============================================================================
+Reading chain files from PolyChord, MultiNest, UltraNest, NestedFit, CosmoMC, or Cobaya
+=======================================================================================
 
 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
+* NestedFit: https://github.com/martinit18/nested_fit 
 * `CosmoMC <https://cosmologist.info/cosmomc/readme.html>`_: https://github.com/cmbant/CosmoMC
 * `Cobaya <https://cobaya.readthedocs.io>`_: https://github.com/CobayaSampler/cobaya
 
@@ -38,6 +39,14 @@ out the examples listed here.
       from anesthetic import read_chains
       samples = read_chains("anesthetic/tests/example_data/un")
 
+* NestedFit 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/nf")
+
 * Cobaya samples, which will be an instance of the
   :class:`anesthetic.samples.MCMCSamples` class:
 

pyproject.toml

@@ -21,6 +21,7 @@ authors = [
     { name="Stefan Heimersheim", email="[email protected]" },
     { name="Pablo Lemos", email="[email protected]" },
     { name="Toby Lovick", email="[email protected]"},
+    { name="Deborah Odunuyi", email="[email protected]"},
     { name="Aleksandr Petrosyan", email="[email protected]" },
     { name="Liangliang Su", email="[email protected]"},
     { name="David Yallup", email="[email protected]" },
@@ -58,7 +59,7 @@ classifiers = [
 "JOSS paper" = "https://joss.theoj.org/papers/10.21105/joss.01414"
 
 [project.optional-dependencies]
-docs = ["sphinx", "sphinx_rtd_theme", "sphinx-copybutton", "numpydoc"]
+docs = ["sphinx>=4.2.0", "sphinx_rtd_theme>=1.2.2", "sphinx-copybutton", "numpydoc"]
 test = ["pytest", "pytest-cov", "flake8", "pydocstyle", "packaging", "pre-commit"]
 ultranest = ["h5py"]
 astropy = ["astropy"]