Intialize
[1]:
import classy_sz
[3]:
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from classy_sz import Class
import os
import time
from scipy.integrate import simpson
import numpy as np
font = {'size' : 16, 'family':'STIXGeneral'}
axislabelfontsize='large'
matplotlib.rc('font', **font)
plt.rcParams.update({
"text.usetex": True,
"font.family": "sans-serif",
"font.sans-serif": ["Helvetica"]})
def l_to_dl(lp):
if type(lp) == list:
lp = np.asarray(lp)
return lp*(lp+1.)/2./np.pi
[54]:
import os
path_to_class_sz = os.environ['PATH_TO_CLASS_SZ_DATA']
path_to_class_sz_files = path_to_class_sz+"/class_sz/class-sz/class_sz_auxiliary_files/includes/"
set cosmological parameters
Let us first set the cosmological parameters
[5]:
cosmo_params = {
'omega_b': 0.022,
'omega_cdm': 0.122,
'H0': 67.5,
'tau_reio': 0.0561,
'ln10^{10}A_s': 3.047,
'n_s': 0.965,
}
deal with dndz
put the relevant file in
path_to_class_sz_files
[6]:
zs = np.loadtxt(path_to_class_sz_files+'z_zsplit2.txt')
nzs = np.loadtxt(path_to_class_sz_files+'nz_zsplit2bin0.txt')
check normalization
[10]:
if np.trapezoid(nzs,zs) != 1:
print(np.trapezoid(nzs,zs))
# normalize if not normalized
nzs_norm = nzs/np.trapezoid(nzs,zs)
else:
nzs_norm
1355882.8763808226
save in relevant format
[11]:
np.savetxt(path_to_class_sz_files+'/dndz_quaia.txt',np.c_[zs,nzs_norm])
quaia bias redshift dependence
Inferred b(z) for quaia from Alonso et al .
[12]:
@np.vectorize
def bquaia(z,bl):
return bl*(0.278*((1 + z)**2 - 6.565) + 2.393)
compute power spectra
[13]:
class_sz = Class()
class_sz.set(cosmo_params)
common_params = {
# here gal_gal_1h and gal_gal_2h correspond to the HOD calculation
# while gal_gal_hf is the linear bias calculation
# dndlnm ensures we can access the halo-function after computation
'output': 'gal_gal_1h,gal_gal_2h,gal_gal_hf,dndlnm',
'galaxy_sample' : 'custom', # custom made galaxy sample
'full_path_to_dndz_gal':path_to_class_sz_files+'/dndz_quaia.txt', # redshift distribution of custom made galaxy sample
# pick redshift range
'z_min' : 0.005,
'z_max' : 3.0,
'has_b_custom1': 1, # do you want a cutom function for b(z)? yes: 1, no: 0
'array_b_custom1_n_z': 200, # precision parameters for z sampling of b(z)
'non_linear' : 'hmcode', # non-linear model for matter P(k), e.g., 'hmcode' or 'halofit'
'effective_galaxy_bias': 2.,
# for the custom bias function we want:
# b(z) = bquaia(z,bl) as defined above and bl is effective_galaxy_bias
# We code this in path_to_classy_szfast/classy_szfast/custom_bias/custom_bias.py
# pick mass range (relevant for HOD prediction)
'M_min' : 1e11,
'M_max' : 1e15,
# Other parameters relevant for HOD predictions
'mass_function' : 'T08M200c', # halo mass function, e.g., Tinker et al 2008 defined at M200c
'concentration_parameter' : 'B13', # halo concentration parameter
# main HOD parameters controlling Nc, Ns
'sigma_log10M_HOD': 0.68,
'alpha_s_HOD': 1.30,
'M1_prime_HOD': 10**12.7, # msun/h
'M_min_HOD': 10**11.8, # msun/h
'M0_HOD' :0,
# these are not important (maybe for fine-tuning HOD)
'x_out_truncated_nfw_profile_satellite_galaxies': 1.,
'f_cen_HOD' : 1.,
'hm_consistency' : 0,
}
# if you need other cosmological parameters
# do:
cosmo_params.update(
{
'omega_cdm':0.11,
'cosmo_model': 1,
'skip_background_and_thermo': 0,
'N_ncdm': 1,
'N_ur': 2.0308,
'm_ncdm': 0.06,
'deg_ncdm': 1,
'skip_input': 0,
}
)
class_sz.set(common_params)
[13]:
True
[14]:
%%time
class_sz.compute_class_szfast()
/Users/boris/venvdir/class_sz_312_brew/lib/python3.12/site-packages/mcfit/mcfit.py:130: UserWarning: use backend='jax' if desired
warnings.warn("use backend='jax' if desired")
CPU times: user 3.34 s, sys: 113 ms, total: 3.45 s
Wall time: 812 ms
[15]:
## check why returns 0?
# class_sz.get_current_derived_parameters(['m_ncdm_tot'])
[16]:
class_sz.Neff()
[16]:
3.0441685911057843
collect the output
[17]:
# factor for converting to cl
# linear bias calculations
l, cl_gg_linear_bias = np.asarray(class_sz.cl_gg()['ell']),np.asarray(class_sz.cl_gg()['hf'])
cl_gg_linear_bias /= l_to_dl(l)
# hod calculation
l, cl_gg_hod_1h, cl_gg_hod_2h = np.asarray(class_sz.cl_gg()['ell']),np.asarray(class_sz.cl_gg()['1h']),np.asarray(class_sz.cl_gg()['2h'])
cl_gg_hod_tot = cl_gg_hod_1h + cl_gg_hod_2h
cl_gg_hod_tot /= l_to_dl(l)
cl_gg_hod_1h /= l_to_dl(l)
cl_gg_hod_2h /= l_to_dl(l)
plot
[18]:
label_size = 17
title_size = 18
legend_size = 13
handle_length = 1.5
fig, (ax1) = plt.subplots(1,1,figsize=(5,5))
ax = ax1
ax.tick_params(axis = 'x',which='both',length=5,direction='in', pad=10)
ax.tick_params(axis = 'y',which='both',length=5,direction='in', pad=5)
ax.xaxis.set_ticks_position('both')
ax.yaxis.set_ticks_position('both')
plt.setp(ax.get_yticklabels(), rotation='horizontal', fontsize=label_size)
plt.setp(ax.get_xticklabels(), fontsize=label_size)
ax.grid( visible=True, which="both", alpha=0.2, linestyle='--')
ax.set_ylabel(r"$C_\ell$ ",size=title_size)
ax.set_xlabel(r"$\ell$",size=title_size)
ax.set_xscale('log')
ax.set_yscale('log')
plt.loglog(l, cl_gg_hod_tot,ls='-', label='Total HOD spectrum',c='blue')
plt.loglog(l, cl_gg_hod_1h,ls='--', label='1-halo term',c='blue')
plt.loglog(l, cl_gg_hod_2h,ls='-.', label='2-halo term',c='blue')
plt.loglog(l, cl_gg_linear_bias, label='Total linear-bias spectrum',c='r')
# plt.ylim(1e-2,1e4)
ax.legend(fontsize=11,ncol=1,frameon=False)
ax.legend(fontsize=13,ncol=1,frameon=False)
fig.tight_layout()
plot HOD quantities
In the code the HOD is coded in these functions:
def get_N_satellites(self,double z,double M_halo,double Nc_mean,double M_min,double alpha_s,double M1_prime):
Ns = HOD_mean_number_of_satellite_galaxies(z,M_halo,Nc_mean,M_min,alpha_s,M1_prime,&self.tsz,&self.ba)
return Ns
def get_N_centrals(self,double z,double M_halo,double M_min,double sigma_log10M,double fc):
Nc = HOD_mean_number_of_central_galaxies(z,M_halo,M_min,sigma_log10M,fc,&self.tsz,&self.ba)
return Nc
[19]:
zx = 1.
ngbar = class_sz.get_ng_bar_at_z(zx)
print(f'ngbar = {ngbar}')
ngbar = 0.0119082187262307
[20]:
marr = np.geomspace(common_params['M_min'],common_params['M_max'],500)
M_min = common_params['M_min_HOD']
sigma_log10M = common_params['sigma_log10M_HOD']
fc = common_params['f_cen_HOD']
get_nc = np.vectorize(class_sz.get_N_centrals)
nc = get_nc(zx,marr,M_min,sigma_log10M,fc).squeeze()
get_ns = np.vectorize(class_sz.get_N_satellites)
alpha_s = common_params['alpha_s_HOD']
m1p = common_params['M1_prime_HOD']
ns = get_ns(zx,marr,nc,M_min,alpha_s,m1p)
[21]:
plt.plot(marr,nc)
plt.plot(marr,ns)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('mass')
plt.ylabel('Nc,Ns')
[21]:
Text(0, 0.5, 'Nc,Ns')
Halo bias (Tinker 2010)
[25]:
# "Tinker bias", i.e., halo bias
# here we pick z = zx defined ealier
bh = class_sz.get_first_order_bias_at_z_and_m(zx,marr)[0]
[26]:
bh.shape
[26]:
(500,)
[27]:
plt.plot(marr,bh)
plt.yscale('log')
plt.xlabel('mass msun/h')
plt.ylabel('b(m,z)')
[27]:
Text(0, 0.5, 'b(m,z)')
get halo mass function
[28]:
dndlnm = class_sz.get_dndlnm(zx, marr) # halo mass function
dndm = dndlnm/marr
compute formula if needed
[38]:
# assume Nc = 1
# need to fin Mmin such that beff(z) = RHS
@np.vectorize
def beff(mmin):
integrand = (dndm*bh)[0]
integrand[marr<mmin] = 0.
return np.trapezoid(integrand,x=marr)/ngbar
[39]:
mmin = np.geomspace(1e11,1e15,200)
bs = beff(mmin)
Define bl for that sample
[40]:
bl = 1.3
Compare
[42]:
fig, ax = plt.subplots(figsize=(10,5))
plt.plot(mmin,bs)
ax.axhline(bquaia(zx,bl),label=f'measured quia bias at z = {zx} and bl = {bl}',c='r')
plt.xscale('log')
plt.xlabel(r'$M_\mathrm{min} [M_\mathrm{sun}/h]$')
plt.ylabel(f'Quaia bias estmate at z = {zx}')
plt.legend()
[42]:
<matplotlib.legend.Legend at 0x388b2ab10>
Quaia Six bins
[48]:
for i in range(6):
zs = np.loadtxt(path_to_class_sz_files+'zs_zsplit6.txt')
nzs = np.loadtxt(path_to_class_sz_files+f'nzs_zsplit6bin{i}.txt')
if np.trapezoid(nzs,zs) != 1:
print(np.trapezoid(nzs,zs))
# normalize if not normalized
nzs_norm = nzs/np.trapezoid(nzs,zs)
else:
nzs_norm
np.savetxt(path_to_class_sz_files+f'/dndz_quaia_{i}.txt',np.c_[zs,nzs_norm])
71481.04066872306
473301.79462171864
668156.2756827383
580528.2260208138
280667.9358349544
140526.58722148684
[49]:
import classy_szfast
import classy_sz
[51]:
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from classy_sz import Class
import os
import time
from scipy.integrate import simpson
import numpy as np
font = {'size' : 16, 'family':'STIXGeneral'}
axislabelfontsize='large'
matplotlib.rc('font', **font)
plt.rcParams.update({
"text.usetex": True,
"font.family": "sans-serif",
"font.sans-serif": ["Helvetica"]})
def l_to_dl(lp):
if type(lp) == list:
lp = np.asarray(lp)
return lp*(lp+1.)/2./np.pi
[55]:
# path_to_class_sz = "/your/path/class_sz/class-sz/"
# path_to_class_sz_files = path_to_class_sz+"class_sz_auxiliary_files/excludes/galaxies/"
The history saving thread hit an unexpected error (OperationalError('unable to open database file')).History will not be written to the database.
[56]:
path_to_class_sz_files+f'dndz_quaia_'
[56]:
'/Users/boris/class_sz_data_directory/class_sz/class-sz/class_sz_auxiliary_files/includes/dndz_quaia_'
[57]:
# for z dependence
# set 'has_tracer_bias_zdependence = 1'
# this is hard coded in perturbations.c
cosmo_params = {
'omega_b': 0.022,
'omega_cdm': 0.122,
'H0': 67.5,
'tau_reio': 0.0561,
'ln10^{10}A_s': 3.047,
'n_s': 0.965,
}
class_sz = Class()
class_sz.set(cosmo_params)
class_sz.set({
'output' : 'galn_galn_hf,galn_lens_hf',
'galaxy_samples_list_num' : 6, # the number of galaxy samples
'galaxy_samples_list' : '0,1,2,3,4,5', # the id string of each sample, can be any integer
'full_path_and_prefix_to_dndz_ngal': path_to_class_sz_files+f'dndz_quaia_',
})
class_sz.set({# class_sz parameters:
'z_min' : 0.,
'z_max' : 4.635,
'effective_galaxy_bias_ngal_0' : 1.,
'effective_galaxy_bias_ngal_1' : 1.,
'effective_galaxy_bias_ngal_2' : 1.,
'effective_galaxy_bias_ngal_3' : 1.,
'effective_galaxy_bias_ngal_4' : 1.,
'effective_galaxy_bias_ngal_5' : 1.,
'dlogell' : 0.1,
'ell_max' : 5000.0,
'ell_min' : 2.0,
'redshift_epsrel' : 1e-4,
'ndim_redshifts': 30,
'ngal_ngal_auto_only' : 1,
# 'has_tracer_bias_zdependence': 1,
'non_linear' : 'hmcode',
'cosmo_model' : 1,
})
[57]:
True
[58]:
# class_sz.pars
[59]:
%%time
class_sz.compute_class_szfast()
# class_sz.compute()
cl_galn_galn = class_sz.cl_galn_galn()
cl_galn_lens = class_sz.cl_galn_lens()
CPU times: user 373 ms, sys: 61 ms, total: 434 ms
Wall time: 318 ms
/Users/boris/venvdir/class_sz_312_brew/lib/python3.12/site-packages/mcfit/mcfit.py:130: UserWarning: use backend='jax' if desired
warnings.warn("use backend='jax' if desired")
[67]:
fig, axes = plt.subplots(figsize=(10, 7), sharex=True,
sharey=True,
ncols=3,
nrows=2,
)
plt.subplots_adjust(left=0.1,
bottom=0.1,
right=0.9,
top=0.9,
wspace=0.4,
hspace=0.4)
ik = 0
for i in range(2):
for j in range(3):
# axes[i, j].plot(x, np.sin((i+j) *x))
axes[i, j].set_title(ik,size=15)
axes[i, j].set_xlabel(r"$\ell$",fontsize=18)
axes[i, j].set_ylabel(r"$C_l^{gg}$",fontsize=18)
axes[i, j].grid()
#
cl = cl_galn_galn[f'{ik}x{ik}']
fac = np.asarray(cl['ell'])*(np.asarray(cl['ell'])+1.)/2./np.pi
axes[i, j].loglog(cl['ell'],np.asarray(cl['hf'])/fac,'r-',label=r'hf')
axes[i, j].legend(loc=1,fontsize=12)
ik+=1
[71]:
plot_dim = len(cl_galn_lens.keys())
fig, axes = plt.subplots(figsize=(10, 7),
sharex=True,
sharey=True,
ncols=3,
nrows=2,
)
plt.subplots_adjust(
hspace=0.4)
ik = 0
for i in range(2):
for j in range(3):
kk = list(cl_galn_lens.keys())[ik]
ik+=1
# axes[i, j].plot(x, np.sin((i+j) *x))
axes[i,j].set_title(kk,size=9)
axes[i,j].set_xlabel(r"$\ell$",fontsize=18)
axes[i,j].set_ylabel(r"$C_l^{gk} \,\,\, [\mathrm{units}]$",fontsize=18)
axes[i,j].grid()
cl = cl_galn_lens[kk]
fac = np.asarray(cl['ell'])*(np.asarray(cl['ell'])+1.)/2./np.pi
axes[i,j].loglog(cl['ell'],np.asarray(cl['hf'])/fac,'r-',label=r'hf')
axes[i,j].legend(loc=1,fontsize=12)
