From f9f8f7436e55a585121d6da0416aa17c1fe377fd Mon Sep 17 00:00:00 2001 From: Larry Bradley Date: Fri, 8 Dec 2023 16:17:03 -0500 Subject: [PATCH] Add custom function to read MultiNest output data --- bagpipes/fitting/fit.py | 39 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/bagpipes/fitting/fit.py b/bagpipes/fitting/fit.py index 6ff6014..8f5fe5e 100644 --- a/bagpipes/fitting/fit.py +++ b/bagpipes/fitting/fit.py @@ -2,6 +2,7 @@ import numpy as np import os +import re import time import warnings import h5py @@ -44,6 +45,41 @@ from .posterior import posterior +def _read_multinest_data(filename): + """ + Read MultiNest data. + + By default, Fortran drops the "E" symbol for 3-digit exponent output + (e.g., '0.148232-104'). This impacts the output files currently + being written by MultiNest. For such caes, this reader inserts + the "E" symbol into the number string so that the number can be + converted to a float. + + Parameters + ---------- + filename : str + The filename to read. + """ + # count the columns in the first row of data; + # without a converter, genfromtxt will read 3-digit exponents as np.nan + ncolumns = np.genfromtxt(filename, max_rows=1).shape[0] + + # insert "E" before the "+" or "-" exponent symbol if it is missing, + # so the string can be converted to a float + # '1.148232-104' -> 1.148232e-104 + # '1.148232+104' -> 1.148232e+104 + # '-1.148232-104' -> -1.148232e-104 + # '+1.148232+104' -> 1.148232e+104 + # '0.148232-104' -> 1.482320e-105 + # '0.148232E-10' -> 1.482320e-011 + # '1.148232' -> 1.48232e+000 + convert = lambda s: float(re.sub(r'(\d)([\+\-])(\d)', r'\1E\2\3', + s.decode())) + converters = dict(zip(range(ncolumns), [convert] * ncolumns)) + + return np.genfromtxt(filename, converters=converters) + + class fit(object): """ Top-level class for fitting models to observational data. Interfaces with MultiNest to sample from the posterior distribution @@ -169,7 +205,8 @@ def fit(self, verbose=False, n_live=400, use_MPI=True, sampler="multinest"): # Load MultiNest outputs and save basic quantities to file. if sampler == "multinest": - samples2d = np.loadtxt(self.fname + "post_equal_weights.dat") + multinest_fname = self.fname + 'post_equal_weights.dat' + samples2d = _read_multinest_data(multinest_fname) lnz_line = open(self.fname + "stats.dat").readline().split() self.results["samples2d"] = samples2d[:, :-1] self.results["lnlike"] = samples2d[:, -1]