Inconsistency between class' internal $X_e$ and the one form 21cmFAST
from classy import Class
import os
import io
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import plotly.offline as py
py.init_notebook_mode()
from dautil.plot import iplot_column_slider
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)
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
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
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()
!tail /home/kolen/git/fork/class/base_2018_plikHM_TTTEEE_lowl_lowE_lensing_custom_x_e.ini
!ls ~/git/fork/class/output/
df = read_class_txt(os.path.expanduser('~/git/fork/class/output/base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_cl_lensed.dat'), camb=True)
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)
plot_compare(df, df2, ('default', 'custom X_e'), relative=False)
path = '~/git/fork/class/output/base_2018_plikHM_TTTEEE_lowl_lowE_lensing_camb_custom_x_e_thermodynamics.dat'
!head -n 15 $path
df_thermo = read_class_txt(os.path.expanduser(path), camb=True)
df_thermo_z_restricted = df_thermo.loc[(7 < df_thermo.index) & (df_thermo.index < 30)]
py.iplot(iplot_column_slider(df_thermo_z_restricted))