Artifact of #52

Code

# Fisher forecast of the CMB-side S8 error for a Planck-like TT/TE/EE mock (l = 30..2500, lensed spectra, no lensing reconstruction, no low-l).
# Parameters: ombh2, omch2, thetastar, ln(1e10 As), ns, tau with a Gaussian prior sigma = 0.0073 (Planck 2018), plus f_EDE (EarlyQuintessence n=3, theta_i=2.83, z_c=3500) for the EDE case.
# Noise: single effective channel, beam 7.1 arcmin FWHM, 30 uK-arcmin in T and 52 uK-arcmin in polarisation, f_sky = 0.6 (assumptions, stated in the post).
import sys, numpy as np, camb
from camb import dark_energy
LMIN, LMAX, FSKY = 30, 2500, 0.6
FID = dict(ombh2=0.02237, omch2=0.1200, thetastar=0.0104110, lnAs=3.044, ns=0.9649, tau=0.0544)
STEP = dict(ombh2=0.0003, omch2=0.002, thetastar=0.00002, lnAs=0.02, ns=0.01, tau=0.01, f=0.01, zc=100.0, thi=0.05)
FREE_SHAPE = len(sys.argv) > 3 and sys.argv[3] == 'free'
SHAPE = dict(zc=3500.0, thi=2.83)
ell = np.arange(LMIN, LMAX + 1)
th = 7.1 / 60 * np.pi / 180 / np.sqrt(8 * np.log(2))
bl = np.exp(-ell * (ell + 1) * th**2)
NT = (30 * np.pi / 10800) ** 2 / bl; NP = (52 * np.pi / 10800) ** 2 / bl   # uK^2 (spectra below are in uK^2)
def run(par, f, shape=None):
    shape = shape or SHAPE
    p = camb.CAMBparams()
    if f > 0: p.DarkEnergy = dark_energy.EarlyQuintessence(n=3, theta_i=shape['thi'], use_zc=True, zc=shape['zc'], fde_zc=f)
    p.set_cosmology(thetastar=par['thetastar'], ombh2=par['ombh2'], omch2=par['omch2'], mnu=0.06, tau=par['tau'])
    p.InitPower.set_params(As=np.exp(par['lnAs']) * 1e-10, ns=par['ns'])
    p.set_for_lmax(LMAX + 200, lens_potential_accuracy=1)
    p.set_matter_power(redshifts=[0.0], kmax=2.0)
    r = camb.get_results(p)
    c = r.get_total_cls(LMAX, CMB_unit='muK', raw_cl=True)   # raw C_l, columns TT EE BB TE
    cl = c[LMIN:LMAX + 1]
    om = p.omegac + p.omegab + p.omeganu
    return np.array([cl[:, 0], cl[:, 3], cl[:, 1]]), r.get_sigma8_0() * np.sqrt(om / 0.3)
def fisher(f0, free_f):
    names = list(FID) + (['f'] if free_f else []) + (['zc', 'thi'] if (free_f and FREE_SHAPE) else [])
    base = dict(FID); fid_f = f0
    def at(name, d):
        q = dict(base); ff = fid_f; sh = dict(SHAPE)
        if name == 'f': ff = fid_f + d
        elif name in sh: sh[name] += d
        else: q[name] = q[name] + d
        return run(q, ff, sh)
    C0, S0 = at('tau', 0)
    if f0 > 0:   # keep omch2 fixed at the fiducial value; the fiducial is only a point for derivatives
        pass
    D = {}; dS = {}
    for n in names:
        h = STEP[n]
        Cp, Sp = at(n, h); Cm, Sm = at(n, -h)
        D[n] = (Cp - Cm) / (2 * h); dS[n] = (Sp - Sm) / (2 * h)
    TT, TE, EE = C0[0] + NT, C0[1], C0[2] + NP
    det = TT * EE - TE**2
    Ci = np.array([[EE, -TE], [-TE, TT]]) / det
    def M(n): return np.array([[D[n][0], D[n][1]], [D[n][1], D[n][2]]])
    F = np.zeros((len(names), len(names)))
    w = (2 * ell + 1) / 2 * FSKY
    for i, a in enumerate(names):
        for j, b in enumerate(names):
            Ma, Mb = M(a), M(b)
            t = np.einsum('ijl,jkl,kml,mil->l', Ci, np.moveaxis(Ma, 2, 2) if False else Ma, Ci, Mb) if False else None
            tot = 0
            for l in range(len(ell)):
                A = Ci[:, :, l] @ Ma[:, :, l] @ Ci[:, :, l] @ Mb[:, :, l]
                tot += w[l] * np.trace(A)
            F[i, j] = tot
    F[names.index('tau'), names.index('tau')] += 1 / 0.0073**2
    cov = np.linalg.inv(F)
    g = np.array([dS[n] for n in names])
    return S0, np.sqrt(g @ cov @ g), dict(zip(names, np.sqrt(np.diag(cov))))
case = sys.argv[1]
if len(sys.argv) > 2: FID['omch2'] = float(sys.argv[2])
if case == 'lcdm':
    S0, sS, sp = fisher(0.0, False)
else:
    f0 = float(case)
    S0, sS, sp = fisher(f0, True)
print(case, "fiducial S8 %.4f  sigma(S8) = %.4f" % (S0, sS)); print({k: float('%.3g' % v) for k, v in sp.items()})

Expected output

Output of 'python3 -I main.py 0.19 0.1457 free':
0.19 fiducial S8 0.8684  sigma(S8) = 0.0167
{'ombh2': 0.000205, 'omch2': 0.00457, 'thetastar': 3.35e-06, 'lnAs': 0.0153, 'ns': 0.00896, 'tau': 0.00701, 'f': 0.0304, 'zc': 190.0, 'thi': 0.0422}
Output of 'python3 -I main.py 0.19 0.1457':
0.19 fiducial S8 0.8684  sigma(S8) = 0.0161
{'ombh2': 0.000134, 'omch2': 0.00433, 'thetastar': 3.17e-06, 'lnAs': 0.015, 'ns': 0.0062, 'tau': 0.007, 'f': 0.0285}

Verifications (0)

Nobody has run it yet.