Comparing 21cmFAST with 21cmEMU directly

In this tutorial, we will use 21cmFAST to simulate a lightcone and compare it with the output from 21cmEMU. We will compare the EoR history and global signal.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import py21cmfast as p21c

from py21cmemu import Emulator
from py21cmemu.properties import (
    COSMO_PARAMS,
    FLAG_OPTIONS,
    USER_PARAMS,
    emulator_properties,
)

limits = emulator_properties.limits

/home/dani/anaconda3/envs/21cmEMU/lib/python3.9/site-packages/py21cmfast/_cfg.py:57: UserWarning: Your configuration file is out of date. Updating...
  warnings.warn(
/home/dani/anaconda3/envs/21cmEMU/lib/python3.9/site-packages/py21cmfast/_cfg.py:41: UserWarning: Your configuration file is out of date. Updating...
  warnings.warn("Your configuration file is out of date. Updating...")

p21c.__version__

'3.3.1.dev6+g95523bc'

# Let's take a test param that is already in 21cmFAST-friendly units
test_param = [
    -0.98454527,
    0.84028646,
    -1.01608287,
    0.03414988,
    9.02499104,
    0.45168016,
    40.0,
    500.0,
    1.0,
]
keys = [
    "F_STAR10",
    "ALPHA_STAR",
    "F_ESC10",
    "ALPHA_ESC",
    "M_TURN",
    "t_STAR",
    "L_X",
    "NU_X_THRESH",
    "X_RAY_SPEC_INDEX",
]

input_dict = {k: v for k, v in zip(keys, test_param, strict=False)}

input_dict

{'F_STAR10': -0.98454527,
 'ALPHA_STAR': 0.84028646,
 'F_ESC10': -1.01608287,
 'ALPHA_ESC': 0.03414988,
 'M_TURN': 9.02499104,
 't_STAR': 0.45168016,
 'L_X': 40.0,
 'NU_X_THRESH': 500.0,
 'X_RAY_SPEC_INDEX': 1.0}

ap = p21c.AstroParams(input_dict)

lightcone = p21c.run_lightcone(
    redshift=6.0,
    max_redshift=15.0,
    astro_params=ap,
    user_params=USER_PARAMS,
    lightcone_quantities=("brightness_temp", "density"),
    global_quantities=("brightness_temp", "density", "xH_box"),
    flag_options=FLAG_OPTIONS,
    cosmo_params=COSMO_PARAMS,
    direc="_cache",
)

np.savez(
    "lightcone_example",
    redshifts=lightcone.node_redshifts,
    xHI=lightcone.global_xHI,
    Tb=lightcone.global_brightness_temp,
    pars=input_dict,
)

lightcone = np.load("lightcone_example.npz")

emu = Emulator()

theta, out, errs = emu.predict(input_dict)

zs = out.redshifts
xHI = out.xHI
Tb = out.Tb

plt.figure(figsize=(10, 6))
fs = 25
plt.plot(zs, xHI, ls="-.", color="k", label="21cmEMU")
# plt.plot(lightcone.node_redshifts, lightcone.global_xHI, color = 'k', label = '21cmFAST')
plt.plot(lightcone["redshifts"], lightcone["xHI"], color="k", label="21cmFAST")

plt.legend(fontsize=fs)
plt.xticks(fontsize=fs)
plt.yticks(fontsize=fs)
plt.ylabel(r"$\overline{x}$HI", fontsize=fs)
plt.xlabel(r"Redshift z", fontsize=fs)
plt.show()

# Directly compare Tb from test set and Tb predicted by emulator
N = 7
cs = ["k", "r", "lime", "b", "orange", "cyan", "magenta", "coral"]
plt.figure(figsize=(10, 6))

plt.plot(zs, Tb, ls="-.", color="k", label="21cmEMU")
# plt.plot(lightcone.node_redshifts, lightcone.global_brightness_temp, color = 'k', label = '21cmFAST')
plt.plot(lightcone["redshifts"], lightcone["Tb"], color="k", label="21cmFAST")
plt.legend(fontsize=fs)
plt.xticks(fontsize=fs)
plt.yticks(fontsize=fs)
plt.ylabel(r"$\overline{T}_b$ [mK]", fontsize=fs)
plt.xlabel(r"Redshift z", fontsize=fs)
plt.show()