camb_vs_cobaya.py

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import cobaya
import camb
print("     Numpy :", np.__version__)
print("Matplotlib :", mpl.__version__)
print("      CAMB :", camb.__version__)
print("    Cobaya :", cobaya.__version__)

cosmo_params = {
    "cosmomc_theta": 0.0104085,
    "As": 2.1e-9,
    "ombh2": 0.02237,
    "omch2": 0.1200,
    "ns": 0.9649,
    "Alens": 1.0,
    "tau": 0.0544
}

info = {
    "params": cosmo_params,
    "likelihood": {"one": None},
    # "theory": {"camb": None},
    "theory": {"camb": {"extra_args": {"lens_potential_accuracy": 1}}},
}

lmin, lmax = 2, 1000

# CAMB via cobaya
from cobaya.model import get_model
model = get_model(info)
Cl = {"tt": lmax, "ee": lmax, "te": lmax, "bb":lmax}
if cobaya.__version__ < "2.1.0":
    model.likelihood.theory.needs(Cl=Cl)
    model.logposterior({})
    Dls = model.likelihood.theory.get_Cl(ell_factor=True)
else:
    model.theory["camb"].needs(Cl=Cl)
    model.logposterior({})
    Dls = model.theory["camb"].get_Cl(ell_factor=True)

l = Dls["ell"]
dls_cobaya = {cl: Dls[cl] for cl in Cl.keys()}

# Standalone CAMB
camb_params = cosmo_params.copy()
camb_params.update({"lmax": lmax, "lens_potential_accuracy": 1})
pars = camb.set_params(**camb_params)
results = camb.get_results(pars)
powers = results.get_cmb_power_spectra(pars, CMB_unit="muK")
dls_camb = {k: powers["total"][:, v]
            for k, v in {"tt": 0, "ee": 1, "te": 3, "bb": 2}.items()}

# Plot Cl and delta Cl
fig, axes = plt.subplots(2, 2, sharex=True, figsize=(8, 8))
axes = axes.flatten()
for i, cl in enumerate(Cl.keys()):
    ax = axes[2] if cl == "te" else axes[0]
    ax.plot(l, dls_cobaya[cl], "-C{}".format(i), label=cl.upper())
    ax.plot(l, dls_camb[cl], "--C{}".format(i))
    ax = axes[3] if cl == "te" else axes[1]
    ax.plot(l, dls_camb[cl]-dls_cobaya[cl], "-C{}".format(i), label=cl.upper())

axes[0].set_yscale("log")
axes[0].set_ylabel(r"$D_\ell$")
axes[2].set_ylabel(r"$D_\ell$")
axes[1].set_ylabel(r"$\Delta D_\ell$")
axes[3].set_ylabel(r"$\Delta D_\ell$")
axes[2].set_xlabel(r"$\ell$")
axes[3].set_xlabel(r"$\ell$")
for ax in axes:
    ax.legend()

plt.tight_layout()
plt.subplots_adjust(top=0.9)
fig.suptitle("CAMB {} - COBAYA {}".format(camb.__version__, cobaya.__version__))
plt.show()

Results from cobaya 2.0.3 with the temperature fixed to 2.7255

Results from cobaya devel branch (at revision CobayaSampler/cobaya@db9b6e6)