Optical depth SIE model validation

This notebook shows example how to use the OpticalDepth class to create an optical depth distribution for a ‘SIE’ model and make comparison with various models.
You can change model or model parameters and compare the results.

[1]:

import pycbc
import numpy as np
import matplotlib.pyplot as plt
from astropy.cosmology import LambdaCDM
cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7)
from ler.lens_galaxy_population import OpticalDepth

[6]:

od = OpticalDepth()

z_to_Dc interpolator will be loaded from ./interpolator_pickle/z_to_Dc/z_to_Dc_1.pickle
Dc_to_z interpolator will be loaded from ./interpolator_pickle/Dc_to_z/Dc_to_z_1.pickle
angular_diameter_distance interpolator will be loaded from ./interpolator_pickle/angular_diameter_distance/angular_diameter_distance_1.pickle
differential_comoving_volume interpolator will be loaded from ./interpolator_pickle/differential_comoving_volume/differential_comoving_volume_1.pickle
velocity_dispersion_bernardi interpolator will be loaded from ./interpolator_pickle/velocity_dispersion_bernardi/velocity_dispersion_bernardi_0.pickle
optical_depth_SIE_hemanta interpolator will be loaded from ./interpolator_pickle/optical_depth_SIE_hemanta/optical_depth_SIE_hemanta_0.pickle

Let’s look at the available functions and priors

[7]:

od.available_optical_depth_list_and_its_params

[7]:

{'optical_depth_SIS_haris': None,
 'optical_depth_SIS_hemanta': None,
 'optical_depth_SIE_hemanta': None}

Choice 1

optical_depth_SIS_haris (analytical)
velocity_dispersion_gengamma
axis_ratio_rayleigh
cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7)

[10]:

od1 = OpticalDepth(
    npool=4,
    z_min=0.001,
    z_max=10.0,
    optical_depth_function="optical_depth_SIS_haris",
    sampler_priors=dict(
        velocity_dispersion="velocity_dispersion_gengamma",
        axis_ratio="axis_ratio_rayleigh",
    ),
    sampler_priors_params=dict(
        velocity_dispersion=dict(vd_min=0, vd_max=600),
        axis_ratio=dict(q_min=0.2, q_max=1),
    ),
    cosmology=cosmo,
    directory="./interpolator_pickle",
    #create_new_interpolator=True,
    )

z_to_Dc interpolator will be loaded from ./interpolator_pickle/z_to_Dc/z_to_Dc_1.pickle
Dc_to_z interpolator will be loaded from ./interpolator_pickle/Dc_to_z/Dc_to_z_1.pickle
angular_diameter_distance interpolator will be loaded from ./interpolator_pickle/angular_diameter_distance/angular_diameter_distance_1.pickle
differential_comoving_volume interpolator will be loaded from ./interpolator_pickle/differential_comoving_volume/differential_comoving_volume_1.pickle
velocity_dispersion_gengamma interpolator will be loaded from ./interpolator_pickle/velocity_dispersion_gengamma/velocity_dispersion_gengamma_1.pickle

Choice 2

optical_depth_SIS_hemanta (numerical)
velocity_dispersion_bernardi
axis_ratio_rayleigh
cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7)

[12]:

od2 = OpticalDepth(
    npool=4,
    z_min=0.001,
    z_max=10.0,
    optical_depth_function="optical_depth_SIS_hemanta",
    sampler_priors=dict(
        velocity_dispersion="velocity_dispersion_bernardi",
        axis_ratio="axis_ratio_rayleigh",
    ),
    sampler_priors_params=dict(
        velocity_dispersion=dict(vd_min=0., vd_max=600),
        axis_ratio=dict(q_min=0.2, q_max=1),
    ),
    cosmology=cosmo,
    directory="./interpolator_pickle",
    #create_new_interpolator=True,
    )

z_to_Dc interpolator will be loaded from ./interpolator_pickle/z_to_Dc/z_to_Dc_1.pickle
Dc_to_z interpolator will be loaded from ./interpolator_pickle/Dc_to_z/Dc_to_z_1.pickle
angular_diameter_distance interpolator will be loaded from ./interpolator_pickle/angular_diameter_distance/angular_diameter_distance_1.pickle
differential_comoving_volume interpolator will be loaded from ./interpolator_pickle/differential_comoving_volume/differential_comoving_volume_1.pickle
velocity_dispersion_bernardi interpolator will be loaded from ./interpolator_pickle/velocity_dispersion_bernardi/velocity_dispersion_bernardi_0.pickle
optical_depth_SIS_hemanta interpolator will be generated at ./interpolator_pickle/optical_depth_SIS_hemanta/optical_depth_SIS_hemanta_0.pickle

Choice 3

optical_depth_SIE_hemanta (numerical)
velocity_dispersion_bernardi
axis_ratio_rayleigh
cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7)

[15]:

od3 = OpticalDepth(
    npool=4,
    z_min=0.001,
    z_max=10.0,
    optical_depth_function="optical_depth_SIE_hemanta",
    sampler_priors=dict(
        velocity_dispersion="velocity_dispersion_bernardi",
        axis_ratio="axis_ratio_rayleigh",
    ),
    sampler_priors_params=dict(
        velocity_dispersion=dict(vd_min=0, vd_max=600),
        axis_ratio=dict(q_min=0.2, q_max=1),
    ),
    cosmology=cosmo,
    directory="./interpolator_pickle",
    #create_new_interpolator=True,
    )

z_to_Dc interpolator will be loaded from ./interpolator_pickle/z_to_Dc/z_to_Dc_1.pickle
Dc_to_z interpolator will be loaded from ./interpolator_pickle/Dc_to_z/Dc_to_z_1.pickle
angular_diameter_distance interpolator will be loaded from ./interpolator_pickle/angular_diameter_distance/angular_diameter_distance_1.pickle
differential_comoving_volume interpolator will be loaded from ./interpolator_pickle/differential_comoving_volume/differential_comoving_volume_1.pickle
velocity_dispersion_bernardi interpolator will be loaded from ./interpolator_pickle/velocity_dispersion_bernardi/velocity_dispersion_bernardi_0.pickle
optical_depth_SIE_hemanta interpolator will be loaded from ./interpolator_pickle/optical_depth_SIE_hemanta/optical_depth_SIE_hemanta_0.pickle

Choice 4

optical_depth_SIE_hemanta (numerical)
velocity_dispersion_ewoud
axis_ratio_rayleigh
cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7)

[16]:

od4 = OpticalDepth(
    npool=2,
    z_min=0.001,
    z_max=10.0,
    optical_depth_function="optical_depth_SIE_hemanta",
    sampler_priors=dict(
        velocity_dispersion="velocity_dispersion_ewoud",
        axis_ratio="axis_ratio_rayleigh",
    ),
    sampler_priors_params=dict(
        velocity_dispersion=dict(vd_min=0, vd_max=600),
        axis_ratio=dict(q_min=0.2, q_max=1),
    ),
    cosmology=cosmo,
    directory="./interpolator_pickle",
    #create_new_interpolator=True,
    )

z_to_Dc interpolator will be loaded from ./interpolator_pickle/z_to_Dc/z_to_Dc_1.pickle
Dc_to_z interpolator will be loaded from ./interpolator_pickle/Dc_to_z/Dc_to_z_1.pickle
angular_diameter_distance interpolator will be loaded from ./interpolator_pickle/angular_diameter_distance/angular_diameter_distance_1.pickle
differential_comoving_volume interpolator will be loaded from ./interpolator_pickle/differential_comoving_volume/differential_comoving_volume_1.pickle
velocity_dispersion_ewoud interpolator will be generated at ./interpolator_pickle/velocity_dispersion_ewoud/velocity_dispersion_ewoud_0.pickle
optical_depth_SIE_hemanta interpolator will be generated at ./interpolator_pickle/optical_depth_SIE_hemanta/optical_depth_SIE_hemanta_1.pickle

Choice 4

optical_depth_SIE_hemanta (numerical)
velocity_dispersion_ewoud
axis_ratio_rayleigh
cosmo = Planck18

[17]:

from astropy.cosmology import Planck18

od5 = OpticalDepth(
    npool=2,
    z_min=0.001,
    z_max=10.0,
    optical_depth_function="optical_depth_SIE_hemanta",
    sampler_priors=dict(
        velocity_dispersion="velocity_dispersion_ewoud",
        axis_ratio="axis_ratio_rayleigh",
    ),
    sampler_priors_params=dict(
        velocity_dispersion=dict(vd_min=0, vd_max=600),
        axis_ratio=dict(q_min=0.2, q_max=1),
    ),
    cosmology=Planck18,
    directory="./interpolator_pickle",
    #create_new_interpolator=True,
    )

z_to_Dc interpolator will be generated at ./interpolator_pickle/z_to_Dc/z_to_Dc_2.pickle
Dc_to_z interpolator will be generated at ./interpolator_pickle/Dc_to_z/Dc_to_z_2.pickle
angular_diameter_distance interpolator will be generated at ./interpolator_pickle/angular_diameter_distance/angular_diameter_distance_2.pickle
differential_comoving_volume interpolator will be generated at ./interpolator_pickle/differential_comoving_volume/differential_comoving_volume_2.pickle
velocity_dispersion_ewoud interpolator will be generated at ./interpolator_pickle/velocity_dispersion_ewoud/velocity_dispersion_ewoud_1.pickle
optical_depth_SIE_hemanta interpolator will be generated at ./interpolator_pickle/optical_depth_SIE_hemanta/optical_depth_SIE_hemanta_2.pickle

Velocity dispersion

[18]:

# plot hist
sigma1 = od1.sample_velocity_dispersion(size=5000,)
sigma2 = od2.sample_velocity_dispersion(size=5000)
sigma3 = od3.sample_velocity_dispersion(size=5000)
sigma4 = od4.sample_velocity_dispersion(size=5000, zl=1)
sigma5 = od4.sample_velocity_dispersion(size=5000, zl=6)
sigma6 = od5.sample_velocity_dispersion(size=5000, zl=6)
plt.hist(sigma1, bins=50, density=True, alpha=0.5, label=r"$\sigma_v$ haris", histtype="step")
plt.hist(sigma2, bins=50, density=True, alpha=0.5, label="$\sigma_v$ bernardi", histtype="step")
plt.hist(sigma4, bins=50, density=True, alpha=0.5, label="$\sigma_v$ ewoud $z_l=1$", histtype="step")
plt.hist(sigma5, bins=50, density=True, alpha=0.5, label="$\sigma_v$ ewoud $z_l=6$", histtype="step")
plt.hist(sigma6, bins=50, density=True, alpha=0.5, label="$\sigma_v$ $z_l=6$ Planck18", histtype="step")
plt.xlim(0, 600)
plt.xlabel(r"$\sigma_v$")
plt.ylabel(r"$p(\sigma_v)$")
plt.legend()
plt.grid(alpha=0.5)
plt.show()

Optical depth

[20]:

z = np.linspace(0.001, 10, 1000)
# plot od.tau
plt.figure(figsize=(8,4))
plt.plot(z, od1.strong_lensing_optical_depth(z), label=r"SIS $\sigma_v$ haris LambdaCDM")
plt.plot(z, od2.strong_lensing_optical_depth(z), label=r"SIS $\sigma_v$ bernardi LambdaCDM")
plt.plot(z, od3.strong_lensing_optical_depth(z), label=r"SIE $\sigma_v$ bernardi LambdaCDM")
plt.plot(z, od4.strong_lensing_optical_depth(z), label=r"SIE $\sigma(z_l)$ ewoud LambdaCDM")
plt.plot(z, od5.strong_lensing_optical_depth(z), label=r"SIE $\sigma(z_l)$ ewoud planck18")
# plt.plot(z, od4.optical_depth(z), label=r"SIE $\sigma(z_l)$ ewoud planck18")
plt.xlabel(r"$z$")
plt.ylabel(r"$\tau$")
#plt.xlim(0, 2)
#plt.ylim(0, 0.002)
plt.legend()
plt.grid(alpha=0.5)
plt.title("optical depth")
plt.show()

../../_images/examples_optical_depth_validation_SIE_18_0.png

[21]:

z = np.linspace(0.001, 2, 1000)
# plot od.tau
plt.figure(figsize=(8,4))
plt.plot(z, od1.strong_lensing_optical_depth(z), label=r"SIS $\sigma_v$ gengamma LambdaCDM")
plt.plot(z, od2.strong_lensing_optical_depth(z), label=r"SIS $\sigma_v$ bernardi LambdaCDM")
plt.plot(z, od3.strong_lensing_optical_depth(z), label=r"SIE $\sigma_v$ bernardi LambdaCDM")
plt.plot(z, od4.strong_lensing_optical_depth(z), label=r"SIE $\sigma(z_l)$ ewoud LambdaCDM")
plt.plot(z, od5.strong_lensing_optical_depth(z), label=r"SIE $\sigma(z_l)$ ewoud planck18")
# plt.plot(z, od4.optical_depth(z), label=r"SIE $\sigma(z_l)$ ewoud planck18")
plt.xlabel(r"$z$")
plt.ylabel(r"$\tau$")
#plt.xlim(0, 2)
#plt.ylim(0, 0.002)
plt.legend()
plt.grid(alpha=0.5)
plt.title("optical depth")
plt.show()

../../_images/examples_optical_depth_validation_SIE_19_0.png

Distribution of velocity dispersion wrt redshift

reproduction of Oguri et al. (2018b) Fig. 1 with Ewoud’s (code) velocity dispersion model Wempe et al. (2022)

[1]:

from ler.lens_galaxy_population import pdf_phi_z_div_0, phi_loc_bernardi, phi
import numpy as np
import matplotlib.pyplot as plt

[2]:

# combining sigma and q in a function
nsamples_z = 100
nsamples_sigma = 100
z = np.ones(nsamples_sigma)
sigma = np.linspace(50, 400, nsamples_sigma)

[3]:

num1 = pdf_phi_z_div_0(sigma, z*1.0)
num2 = pdf_phi_z_div_0(sigma, z*2.0)
num3 = pdf_phi_z_div_0(sigma, z*3.0)
num4 = pdf_phi_z_div_0(sigma, z*4.0)
num5 = pdf_phi_z_div_0(sigma, z*5.0)
num6 = pdf_phi_z_div_0(sigma, z*6.0)

[4]:

plt.plot(sigma, num1, label='z=1')
plt.plot(sigma, num2, label='z=2')
plt.plot(sigma, num3, label='z=3')
plt.plot(sigma, num4, label='z=4')
plt.plot(sigma, num5, label='z=5')
plt.plot(sigma, num6, label='z=6')
plt.xlabel('sigma')
plt.ylabel(r'$\frac{\phi(\sigma,z)}{\phi(\sigma,z=0)}$')
plt.yscale('log')
plt.ylim(1e-1,0.5e1)
plt.grid(alpha=0.5)
plt.legend()
plt.show()

../../_images/examples_optical_depth_validation_SIE_24_0.png

[ ]: