Summary

Inconsistency between class' internal $X_e$ and the one form 21cmFAST

In [1]:
from classy import Class
/home/kolen/.conda/envs/all3-intel/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 216, got 192
  return f(*args, **kwds)
In [2]:
import os
import io
In [3]:
import pandas as pd
import numpy as np
In [4]:
import matplotlib.pyplot as plt
In [5]:
%matplotlib inline
In [6]:
import plotly.offline as py
py.init_notebook_mode()
In [7]:
from dautil.plot import iplot_column_slider
In [8]:
def read_txt(path):
    with open(path, 'r') as f:
        # remove beginning `#`
        text = f.read()[1:]
    with io.StringIO(text) as f:
        # input has 6 sig. fig. which fits in float32
        return pd.read_csv(f, delim_whitespace=True, index_col=0, dtype=np.float32)
In [9]:
def parse_class_comment(comment):
    '''parse class' comment line and return a header list

    class' comment line example:
      1:z  2:conf. time 3:...
    '''
    # remove beginning '#'
    comment = comment[1:]
    # remove whitespace at boundary
    comemnt = comment.strip()
    # split on :
    comment_list = comment.split(':')
    # throw away first, which is ' 1'
    comment_list = comment_list[1:]

    n = len(comment_list)
    # last doesn't has a number
    for i in range(n - 1):
        # remove last number after white space
        comment_list[i] = ' '.join(comment_list[i].split()[:-1])
    return comment_list
In [10]:
def read_class_txt(path, camb=False):
    # read once
    with open(path, 'r') as f:
        text = f.read()

    # get last comment line
    comment = None
    for line in text.split('\n'):
        if line.startswith('#'):
            comment = line

    names = parse_class_comment(comment)

    with io.StringIO(text) as f:
        df = pd.read_csv(f, delim_whitespace=True, index_col=0, comment='#', header=None, names=names)
        return df if camb else df * 1.e12
In [11]:
def plot_compare(df1, df2, keys, log=False, relative=False):
    '''plot to compare columns from `df1` and `df2`,
    using columns from `df1`.

    `keys`: tuple of str of names of `df1` and `df2`.
    `relative`: if True, plot the relative error w.r.t. `df1` instead.
    '''
    for col in df1.columns:
#         pd.merge(df1[col], df2[col], left_index=True, right_index=True).plot()
        df_temp = pd.concat((df1[col], df2[col]), axis=1, join='inner', keys=(' '.join((key, col)) for key in keys))
        if relative:
            temp = df_temp.values
            if log:
                plt.loglog(df_temp.index, np.abs((temp[:, 0] - temp[:, 1]) / temp[:, 0]), label=col)
            else:
                plt.plot(df_temp.index, np.abs((temp[:, 0] - temp[:, 1]) / temp[:, 0]), label=col)
            plt.legend()
            plt.show()
        else:
            if log:
                df_temp.plot(logx=True, logy=True)
            else:
                df_temp.plot()
In [12]:
!tail /home/kolen/git/fork/class/base_2018_plikHM_TTTEEE_lowl_lowE_lensing_custom_x_e.ini
binned_reio_num = 7
binned_reio_z = 6.665454200000001,8.344146,10.390461700000001,12.8849089,15.9256246,19.6322398,24.150583
binned_reio_xe = 0.00055323428,0.00030092037,0.00022954238000000003,0.00018775178999999999,0.00017818018,0.00018369758,0.00019151189999999998
binned_reio_step_sharpness = 0.1

#reio_inter_num = 79
#reio_inter_z = 0,5.95,5.96404,6.103321,6.245388,6.390296,6.538101,6.688863,6.842639999999999,6.999492999999999,7.159483,7.322673,7.489126,7.658908,7.832087,8.008728999999999,8.188903999999999,8.372682000000001,8.560135000000001,8.751338000000001,8.946363999999999,9.145291,9.348197000000001,9.555161,9.766264,9.981589,10.201221,10.425246,10.653751,10.886827,11.124563,11.367055,11.614396,11.866684,12.124018,12.386498,12.654228,12.927313,13.205859,13.489976,13.779776,14.075371,14.376878,14.684415,14.998103,15.318067,15.644429,15.977318,16.316864,16.6632,17.016464,17.376793,17.744329,18.119217,18.5016,18.891632,19.289465,19.695253,20.109158,20.531342,20.961967,21.401207,21.849232,22.306217,22.772341,23.247787,23.732742,24.227398,24.731945,25.246584,25.771515,26.306946,26.853085,27.410147,27.97835,28.557917,29.149075,29.752056,30.
#reio_inter_xe = -2,-1,0.001117013,0.001023735,0.0009362602,0.000854464,0.0007782149,0.0007073621,0.00064172,0.0005811004,0.0005262021,0.0004780447,0.0004365587000000001,0.0004013149,0.0003715573000000001,0.0003463328,0.000325201,0.0003074503,0.000292537,0.0002800515,0.0002696148,0.0002608876,0.0002536635,0.0002476625,0.0002427183,0.0002386412,0.0002349889,0.0002309913,0.0002263293,0.0002210947,0.0002155082,0.0002098274,0.0002042996,0.0001991403,0.0001944928,0.0001904641,0.0001870542,0.0001842557,0.0001820293,0.0001803089,0.0001790365,0.000178146,0.0001775812,0.0001772763,0.0001771834,0.0001772654,0.0001774832,0.000177809,0.0001782208,0.0001786995,0.0001792314,0.0001798065,0.0001804161,0.0001810538,0.0001817147,0.0001823958,0.0001830941,0.0001838079,0.0001845358,0.0001852772,0.0001860313,0.0001867979,0.0001875768,0.000188368,0.0001891715,0.0001899874,0.0001908159,0.0001916572,0.0001925115,0.0001933791,0.0001942602,0.0001951552,0.0001960643,0.0001969879,0.0001979262,0.0001988797,0.0001998485,0.0002008332,0

write thermodynamics = yes
In [13]:
!ls ~/git/fork/class/output/
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_cl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_cl_lensed.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_cl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_cl_lensed.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_parameters.ini
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_pk_cb.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_pk_cb_nl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_pk.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_pk_nl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_thermodynamics.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_unused_parameters
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_external_cl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_external_cl_lensed.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_external_parameters.ini
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_external_pk_cb.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_external_pk_cb_nl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_external_pk.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_external_pk_nl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_external_unused_parameters
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_parameters.ini
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_pk_cb.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_pk_cb_nl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_pk.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_pk_nl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_unused_parameters
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_cl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_cl_lensed.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_parameters.ini
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_pk_cb.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_pk_cb_nl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_pk.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_pk_nl.dat
base_2018_plikHM_TTTEEE_lowl_lowE_lensing_unused_parameters
explanatory00_cl.dat
explanatory00_cl_lensed.dat
explanatory00_parameters.ini
explanatory00_unused_parameters
reference
In [14]:
df = read_class_txt(os.path.expanduser('~/git/fork/class/output/base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_cl_lensed.dat'), camb=True)
In [15]:
df2 = read_class_txt(os.path.expanduser('~/git/fork/class/output/base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_cl_lensed.dat'), camb=True)
In [16]:
plot_compare(df, df2, ('default', 'custom X_e'), relative=False)

Thermo

In [17]:
path = '~/git/fork/class/output/base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_thermodynamics.dat'
In [18]:
!head -n 15 $path
# Table of selected thermodynamics quantities
# The following notation is used in column titles:
#    x_e = electron ionization fraction
# -kappa = optical depth
# kappa' = Thomson scattering rate, prime denotes conformal time derivatives
#      g = kappa' e^-kappa = visibility function 
#     Tb = baryon temperature 
#  c_b^2 = baryon sound speed squared 
#  tau_d = baryon drag optical depth 
#    1:z                      2:conf. time [Mpc]       3:x_e                    4:kappa' [Mpc^-1]        5:exp(-kappa)            6:g [Mpc^-1]             7:Tb [K]                 8:c_b^2                  9:tau_d              
       0.000000000000e+00       1.415008256241e+04       1.081885046089e+00       4.209937145613e-07       1.000000000000e+00       4.209937145613e-07       2.824481002508e-01       6.984592211472e-14       0.000000000000e+00 
       1.007782391480e-03       1.414559596211e+04       1.081885046089e+00       4.218426822381e-07       9.999981092663e-01       4.218418846459e-07       2.830025812987e-01       7.062110751657e-14       2.786728606699e-09 
       1.426905918761e-02       1.408675623346e+04       1.081885046089e+00       4.330937998900e-07       9.999729564413e-01       4.330820874924e-07       2.903846560033e-01       7.243831081895e-14       4.012965499360e-08 
       2.775825945696e-02       1.402728836678e+04       1.081885046089e+00       4.446902045629e-07       9.999468563821e-01       4.446665721165e-07       2.979818011260e-01       7.430704742560e-14       7.939428281323e-08 
       4.148112112022e-02       1.396719025852e+04       1.081885046089e+00       4.566446939905e-07       9.999197727299e-01       4.566080586333e-07       3.058012782515e-01       7.622897960252e-14       1.206827379361e-07 
In [19]:
df_thermo = read_class_txt(os.path.expanduser(path), camb=True)
In [33]:
df_thermo_z_restricted = df_thermo.loc[(7 < df_thermo.index) & (df_thermo.index < 30)]
In [34]:
py.iplot(iplot_column_slider(df_thermo_z_restricted))
In [ ]: