<jupyter_start><jupyter_text>cosmoprimo examples  Cosmology<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

from cosmoprimo import *
# Start with some vanilla cosmology
cosmo = Cosmology()
# Planck2018 TT, TE, EE, lowE, lensing and BAO data
cosmo_planck = fiducial.Planck2018FullFlatLCDM()
# Custom cosmology
cosmo_custom = Cosmology(omega_cdm=0.2,sigma8=0.7)
print('You can provide the following parameters (including conflicts):',
      cosmo.get_default_parameters(include_conflicts=True))
# Input and easily derived parameters can be accessed through indexing
print('h',cosmo['h'],cosmo_custom['h'])
print('Omega_cdm',cosmo['Omega_cdm'],cosmo_custom['Omega_cdm'])<jupyter_output>h 0.7 0.7
Omega_cdm 0.25 0.4081632653061225<jupyter_text>Engines & sections<jupyter_code># Cosmological calculation is performed by engines, by section 
# (Background, Thermodynamics, Primordial, Perturbations, Transfer, Harmonic, Fourier)
# engines can be set once for all sections as:
cosmo = Cosmology(engine='class')
# or:
cosmo.set_engine('camb')
# or you can pass an engine instance directly
cosmo.set_engine(engine=cosmo.engine) # CLASS engine was attached to cosmo with the previous call
# then you don't need to specify it when calling specific section, e.g.
ba = cosmo.get_background()
# or:
ba = Background(cosmo)
# you can temporarily switch to another engine with:
ba = cosmo.get_background(engine='camb',set_engine=False)
# or:
ba = Background(cosmo,engine='camb',set_engine=False)
assert cosmo.engine is not ba.engine # class engine<jupyter_output><empty_output><jupyter_text>Shortcuts<jupyter_code># If you are lazy, and you want to access directly a quantity which will be in only one section, e.g.:
cosmo = Cosmology(engine='class')
z = 0.2
d = cosmo.comoving_radial_distance(z)
# is equivalent to cosmo.get_background().comoving_radial_distance(z)
assert np.all(d == cosmo.get_background().comoving_radial_distance(z))
# however, you will get into troubles if some attributes/methods 
# of e.g. section "Fourier" are the same as section "Fourier" (e.g. "table")
assert 'table' in dir(cosmo.get_transfer())
assert 'table' in dir(cosmo.get_fourier())
# note that available attributes/methods are accessible through dir():
print('Attributes/methods available with class engine are {}'.format(dir(cosmo)))
# and their doc is given by .__doc__, e.g.
print(cosmo.comoving_radial_distance.__doc__)
# available attributes/methods can be different with a different engine, e.g.
assert set(dir(cosmo)) != set(dir(Cosmology(engine='eisenstein_hu')))<jupyter_output>Attributes/methods available with class engine are ['A_s', 'H0', 'N_eff', 'N_ncdm', 'N_ur', 'Omega0_Lambda', 'Omega0_b', 'Omega0_cdm', 'Omega0_dcdm', 'Omega0_fld', 'Omega0_g', 'Omega0_k', 'Omega0_m', 'Omega0_ncdm', 'Omega0_ncdm_tot', 'Omega0_r', 'Omega0_ur', 'Omega_Lambda', 'Omega_b', 'Omega_cdm', 'Omega_de', 'Omega_fld', 'Omega_g', 'Omega_k', 'Omega_m', 'Omega_ncdm', 'Omega_r', 'Omega_ur', 'T0_cmb', 'T0_ncdm', 'T_cmb', 'T_ncdm', '_RH0_', '__bool__', '__class__', '__copy__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattr__', '__getattribute__', '__getitem__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_check_pk_of', '_default_calculation_parameters', '_default_cosmological_parameters', '_engine', '_get_A_s_fid', '_get_p_ncdm'[...]<jupyter_text>Save/load<jupyter_code># You can save/load any cosmology instance this way:
import os
import tempfile

with tempfile.TemporaryDirectory() as tmp_dir: # in a temporary file to avoid messing up our directory
    fn = os.path.join(tmp_dir, 'cosmo.npy')
    cosmo.save(fn)
    cosmo2 = Cosmology.load(fn)

for param in cosmo.params:
    if isinstance(cosmo[param], np.ndarray):
        assert np.all(cosmo2[param] == cosmo[param])
    else:
        assert cosmo2[param] == cosmo[param]
# Note however that this will not save any calculated quantity<jupyter_output><empty_output><jupyter_text>Background<jupyter_code># you can choose an engine ('class','camb','eisenstein_hu' for approximations) for the calculation
# do ba = Background(cosmo,engine='class') or:
ba = cosmo.get_background(engine='class') # e.g. precision parameters can be passed as other arguments
ba_camb = Background(cosmo,engine='camb') # alternative way of getting background
z = np.linspace(0.,10.,1000)
plt.plot(z,ba.comoving_radial_distance(z),label='class')
plt.plot(z,ba_camb.comoving_radial_distance(z),label='camb')
plt.plot(z,cosmo_planck.get_background(engine='class').comoving_radial_distance(z),label='Planck2018')
# redoing cosmo_planck.get_background(engine='class') will not recompute the background, it's saved one for all!
# (cosmology instance is immutable)
plt.xlabel('$z$')
plt.ylabel('comoving radial distance [$\mathrm{Mpc}/h$]')
plt.legend()
plt.show()
# many quantities are available: densities, distances;
# look at the doc (for class, look at http://github.com/adematti/pyclass)
plt.plot(z,ba.Omega_cdm(z),label='class')
plt.plot(z,ba.Omega_cdm(z),label='Planck2018')
plt.xlabel('$z$')
plt.ylabel('density parameter of cold dark matter')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>You can also access thermodynamics quantities (``cosmo.get_thermodynamics()``), primordial quantities, etc. Let's move to...  Perturbations<jupyter_code>hr = Harmonic(cosmo,engine='class')
cls = hr.unlensed_cl()
ells_factor = (cls['ell'] + 1) * cls['ell'] / (2 * np.pi)
plt.plot(cls['ell'],ells_factor*cls['tt'],label='class')
hr = cosmo.get_harmonic(engine='camb') # second initialisation style
cls = hr.unlensed_cl() # same number of cells are output than in CLASS
plt.plot(cls['ell'],ells_factor*cls['tt'],label='camb')
plt.xlabel('$\ell$')
plt.ylabel('$\ell(\ell+1)C_{\ell}/(2\pi)$')
plt.legend()
plt.show()
# if you want lensed_cl(), you must set lensing=True as part of the Cosmology arguments
fo = Fourier(cosmo,engine='class')
pk = fo.pk_interpolator()
k = np.logspace(-3,2,1000)
z = 0
plt.loglog(k,pk(k,z=z),label='class')
fo_camb = cosmo.get_fourier(engine='camb') # second initialisation style
pk = fo_camb.pk_interpolator() # same number of cells are output than in CLASS
plt.loglog(k,pk(k,z=z),label='camb')
# cross density-velocity CDM + baryon power spectrum
plt.loglog(k,fo.pk_interpolator(of=('delta_cb','theta_cb'))(k,z=z),label=r'$P_{\delta_{cb}\theta_{cb}}$')
plt.xlabel('$k$ [$h/\mathrm{Mpc}$]')
plt.ylabel('$P(k)$ [$(\mathrm{Mpc}/h)^{3}$]')
plt.legend()
plt.show()
# fast approximations
pk = fo.pk_interpolator()
plt.loglog(k,pk(k,z=z),label='class')
plt.loglog(k,Fourier(cosmo,'eisenstein_hu').pk_interpolator()(k,z=z),label='EH1998')
plt.loglog(k,Fourier(cosmo,'eisenstein_hu_nowiggle').pk_interpolator()(k,z=z),label='EH1998 no wiggle')
plt.loglog(k,Fourier(cosmo,'bbks').pk_interpolator()(k,z=z),label='BBKS')
plt.xlabel('$k$ [$h/\mathrm{Mpc}$]')
plt.ylabel('$P(k)$ [$(\mathrm{Mpc}/h)^{3}$]')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Note about $f\sigma_{8}$: should be better computed as the r.m.s. of velocity perturbations<jupyter_code>fo = Fourier(cosmo,engine='class')
print('fsigma8 is',fo.sigma8_z(z=0,of='theta_cb'))
# as a function of z
k = np.logspace(-3,2,1000)
z = np.linspace(0.,4.,10)
colors = plt.cm.jet(np.linspace(0,1,len(z)))
for z_,color in zip(z,colors):
    plt.loglog(k,pk(k,z=z_),color=color)
plt.xlabel('$k$ [$h/\mathrm{Mpc}$]')
plt.ylabel('$P(k)$ [$(\mathrm{Mpc}/h)^{3}$]')
plt.show()<jupyter_output><empty_output><jupyter_text>Interpolators: easy handling of $P(k)$/$\xi(s)$<jupyter_code># 2D (k,z) interpolator (vectorized, of course)
pk = fo.pk_interpolator()
k, z = 0.1, 1.
print('2D power spectrum at k = {}, z = {} is {}'.format(k,z,pk(k,z)))
# cast to 1D z interpolator (tiny interpolation error)
pk1d = pk.to_1d(z=z)
print('1D power spectrum at k = {}, z = {} is {}'.format(k,z,pk1d(k)))
# interpolators can be built from callables, e.g.
pk1dc = PowerSpectrumInterpolator1D.from_callable(pk_callable=pk1d)
print('1D power spectrum from callable at k = {}, z = {} is {}'.format(k,z,pk1dc(k)))
# get sigma_rz() (vectorised)
r,z = 2.,1.
print('quadpack integration sigmarz = {} v.s. FFTlog sigmarz = {}'.format(pk.sigma_rz(r,z,nk=None),pk.sigma_rz(r,z)))
# estimate fsigma_rz()
print('quadpack integration fsigmarz = {} v.s. FFTlog fsigmarz = {}'.format(pk.growth_rate_rz(r,z,nk=None),pk.growth_rate_rz(r,z)))
# FFTlog transform (vectorised) to xi...
xi = pk.to_xi()
s = np.logspace(-2.,2.5,1000)
z = np.linspace(0.,4.,10)
colors = plt.cm.jet(np.linspace(0,1,len(z)))
for z_,color in zip(z,colors):
    plt.plot(s,s**2*xi(s,z=z_),color=color)
plt.xlabel('$s$ [$\mathrm{Mpc}/h$]')
plt.ylabel(r'$s^{2}\xi(s)$ [$(\mathrm{Mpc}/h)^{2}$]')
plt.show()<jupyter_output><empty_output><jupyter_text>BAO filters<jupyter_code># several engines are available for the power spectrum, but wallish2018 is the most accurate one
# note one can also provide a 1D power spectrum interpolator
pknow = PowerSpectrumBAOFilter(pk,engine='wallish2018').smooth_pk_interpolator()
k = np.logspace(-3,2,1000)
z = np.linspace(0.,4.,10)
colors = plt.cm.jet(np.linspace(0,1,len(z)))
for z_,color in zip(z,colors):
    plt.loglog(k,pk(k,z=z_)/pknow(k,z=z_),color=color)
plt.xlabel('$k$ [$h/\mathrm{Mpc}$]')
plt.ylabel('$P(k)/P_{\mathrm{now}}(k)$')
plt.show()
# for the correlation function, the Kirkby2013 approach
xinow = CorrelationFunctionBAOFilter(xi,engine='kirkby2013').smooth_xi_interpolator()
s = np.logspace(-2.,2.5,1000)
z = 1
plt.plot(s,s**2*xi(s,z=z),label='class')
plt.plot(s,s**2*xinow(s,z=z),label='kirkby2013')
plt.xlabel('$s$ [$\mathrm{Mpc}/h$]')
plt.ylabel(r'$s^{2}\xi(s)$ [$(\mathrm{Mpc}/h)^{2}$]')
plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1750100928542592, 'fraction_numerical': 0.0190754945498587, 'mean_word_length': 7.135467980295567, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '649778', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3964, 'n_tokens_neox': 3758, 'n_words': 782, 'path': 'nb/examples.ipynb', 'repo_name': 'adematti/cosmoprimo', 'size': '9908'}