brinckmann / montepython_public

Public repository for the Monte Python Code
MIT License
92 stars 78 forks source link

BK18 likelihood added #370

Open ThomasTram opened 3 months ago

ThomasTram commented 3 months ago

I implemented the BK18 likelihood based on the previous BK14 likelihood. I initially wanted to do it by inheriting from the BK14 likelihood class, but in the end a lot of small things had changed so I decided to copy-paste it and then edit it. There are actually two likelihoods, one using only B-modes (which is the default) called BK18lf_dust and one that includes E-modes which is called BK18lf_dust_incEE. The two likelihoods are identical and differ only in the datafiles and the maps used inside the .data file.

To validate the likelihood I recreated Figure 4 of BKXIII in 2110.00483. Using this script:

import getdist.plots as gdplt
from getdist import MCSamples, loadMCSamples
from pathlib import Path

# Load chain file
samples = loadMCSamples('chains/BK18_fig4_20/2024-04-15_20000_', settings={"ignore_rows": 0.3})
samples.paramNames.setLabelsFromParamNames('chains/BK18_fig4_20/2024-04-15_20000_.paramnames')
# Plot settings
plot_settings = {
    "smooth_scale_1D": 0.15,
    "smooth_scale_2D": 2,
    "num_bins": 100,
    "num_bins_2D": 40,
    "contours": [0.68, 0.95],
    "triangle_plot": True,
    "plot_2D_param": 0,
    "plot_2D_num": 1,
    "params": ["r", "BBdust", "BBsync"],
    "triangle_params": ["r", "BBdust", "BBsync"],
    "param_limits": {"r": [0, 0.18], "BBdust": [0, 11], "BBsync": [0, 8], "BBdustsynccorr": [-1, 1], "BBbetadust": [1, 2], 
                     "BBbetasync": [-4.2, -1.8], "BBalphadust": [-1, 0], "BBalphasync": [-1, 0]},
}

# Make triangle plot of figure 4 of BKXIII
gdplt.get_subplot_plotter()
gdplot = gdplt.getSubplotPlotter()
gdplot.settings.num_plot_contours = 2
gdplot.settings.alpha_filled_add = 0.6
gdplot.triangle_plot(samples, **plot_settings, filled=True)

# Make 1d posterior plots of figure 4 of BKXIII
plot_settings['params'] = ["BBbetadust", "BBbetasync", "BBdustsynccorr", "BBalphadust", "BBalphasync"]
plot_settings['xlims'] = [plot_settings['param_limits'][p] for p in plot_settings['params']]
gdplot.plots_1d(samples, **plot_settings)

I get these figures: image

image

which at least by eye agree.

brinckmann commented 2 months ago

Hi Thomas,

Thanks a lot for doing this, would you like to go ahead and add it to the private devel branch? I can also do it, but will be quite busy the next three weeks and might not get around to it soon. Then once we have everything there (e.g. DESI BAO and new WMAP Python 3 wrapper from @schoeneberg , NPIPE, ground based CMB, etc) we can put out a new release.

If there are any other fixes or improvements you've done that haven't been included yet (e.g. PR #372 ) you can also add those. Thanks a lot!

Best, Thejs

ThomasTram commented 2 months ago

@brinckmann I actually forgot that I had access to the private repo. I can add the changes yes.