Add custom function to read MultiNest output data

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]