[Draft] Support for converting from dynesty and emcee samples

[Stefan-Heimersheim]

Description

[This is a draft, mostly for feedback of how to best implement this. I didn't check for typos/functionality, but I have been using these methods myself.]

Import dynesty and emcee sample objects (in python, not from files) easily. Short version:

dynesty_to_anesthetic = lambda s: anesthetic.samples.NestedSamples(data=s['samples'], weights=np.exp(s['logwt']), logL=s['logl'], columns=columns, tex=tex, limits=prior_dict)
emcee_to_anesthetic = lambda s: anesthetic.samples.MCMCSamples(data=s.flatchain, columns=columns, tex=tex)

Todo

I will probably have time to finalize this in ~2 weeks

• Edit: I just noticed the nested samples currently don't work correctly, e.g. logZ() fails
• Check if any of the unused arguments are relevant (e.g. dynesty_results['logz'], or not passing logL_birth)
• Figure out how/if to do checks like isinstance(dynesty_sampler, dynesty.dynesty.NestedSampler): without adding dependency on dynesty? (maybe something like Iterable)
• There should be a way to obtain the limits from the samplers (DynamicNestedSampler even requires as input a function that returns the limits), but for now it's manual

Checklist:

• [ ] I have performed a self-review of my own code
• [ ] My code is PEP8 compliant (flake8 anesthetic tests)
• [ ] My code contains compliant docstrings (pydocstyle --convention=numpy anesthetic)
• [ ] New and existing unit tests pass locally with my changes (python -m pytest)
• [ ] I have added tests that prove my fix is effective or that my feature works

[williamjameshandley]

Edit: I just noticed the nested samples currently don't work correctly, e.g. logZ() fails Check if any of the unused arguments are relevant (e.g. dynesty_results['logz'], or not passing logL_birth)

The first thing to get NestedSamples working is to make sure that you supply a logL_birth argument. At the moment, this can either be a list of birth contours (which are provided by MultiNest and PolyChord in the dead-birth files, but I don't know if dynesty records these), or a single integer to represent the number of live points. If you're using dynamic nested sampling from dynesty, then we might need to extend this to allow for a list of integers. If you provide these, then logZ() etc should work as normal. You also then won't need to provide a weights argument (since these are computed automatically).

Figure out how/if to do checks like isinstance(dynesty_sampler, dynesty.dynesty.NestedSampler): without adding dependency on dynesty? (maybe something like Iterable)

Pythonically when duck typing you should probably aim for a try-except style framework, rather than specific isinstance checking

There should be a way to obtain the limits from the samplers (DynamicNestedSampler even requires as input a function that returns the limits), but for now it's manual

If you don't supply these then they are computed automatically from the samples #80, and we're thinking about removing this altogether (although I now can't find that conversation).

[Stefan-Heimersheim]

I will probably have time to finalize this in ~2 weeks

Well this was wildly optimistic (also turned out to be more complicated that initially thought) but I aim to work on this as soon as I get a bit of free time

[codecov[bot]]

Codecov Report

Patch coverage: 15.78% and project coverage change: -0.57 :warning:

Comparison is base (d3aafc2) 100.00% compared to head (c03fc1f) 99.43%.

Additional details and impacted files

@@             Coverage Diff             @@
##            master     #158      +/-   ##
===========================================
- Coverage   100.00%   99.43%   -0.57%     
===========================================
  Files           33       33              
  Lines         2790     2809      +19     
===========================================
+ Hits          2790     2793       +3     
- Misses           0       16      +16     

Impacted Files	Coverage Δ
anesthetic/convert.py	38.46% <15.78%> (-61.54%)	:arrow_down:

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

[williamjameshandley]

Hi @Stefan-Heimersheim -- if you can commit an example of typical dynesty/emcee output (in any form), I can write the nlive interface and plumb it into the suite in the same way as is done in #313

[Stefan-Heimersheim]

Here's a typical emcee output. I never ended up using dynesty, and since switched to using polychord only (no longer emcee), so I don't need this functionality, but I'll leave the file & details here in case you want to add it anyway.

File flatchain.csv generated with

import numpy as np
import emcee
import anesthetic
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

def log_prob(p):
    return multivariate_normal.logpdf(p, mean=[0, 0], cov=[[1, 0.5], [0.5, 1]])

def sample_from_prior(nsamples):
    return np.random.uniform(-10, 10, size=(nsamples, 2))

nwalkers = 20
ndim = 2
burnin = 100
nsteps = 100

sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob)
p0 = sample_from_prior(nwalkers)

# Burn in
state = sampler.run_mcmc(p0, burnin)
sampler.reset()

# Sample
sampler.run_mcmc(state, nsteps)
chains = sampler.get_chain()
flatchain = sampler.flatchain

# Save flatchain as csv
np.savetxt("flatchain.csv", flatchain, delimiter=",")

and loaded with

flatchain = np.loadtxt("flatchain.csv", delimiter=",")

columns = ["alpha", "beta"]
samples = anesthetic.samples.MCMCSamples(data=flatchain, columns=columns)

samples.plot_2d(columns)
plt.show()

The example chains are a bit shorter than I would usually choose, but I wanted to minimize file size for the demo.