LeR Short Example
This notebook is created by Phurailatpam Hemantakumar
ler is a comprehensive framework for simulating gravitational wave (GW) events and calculating their detection rates, including gravitational lensing effects. - The ``LeR`` class is the primary interface for these simulations. - The ``GWRATES`` class focuses on standard (unlensed) Compact Binary Coalescence (CBC) events. Refer to the GWRATES complete example for more details.
This notebook demonstrates how to simulate both lensed and unlensed CBC populations and compare their detection rates using the LeR class.
Table of Contents
LeR initialization with default arguments
Basic GW Event Simulation and Rate Calculation
2.1 Simulate Unlensed GW Population
2.2 Calculate Unlensed Detection Rates
2.3 Inspect Generated Unlensed Parameters
Lensed GW Population Simulation and Rates
3.1 Simulate Lensed GW Population
3.2 Calculate Lensed Detection Rates
3.3 Inspect Generated Lensed Parameters
Rate Comparison and Visualization
4.1 Compare Lensed vs Unlensed Rates
4.2 Access Saved Data Files
4.3 Load and Examine Saved Parameters
4.4 Visualize Parameter Distributions
Available Internal Functions and Their Usage
5.1 Explore Functions and Parameters
5.2 Short Example on Using Internal Functions
Summary
1. LeR initialization with default arguments
The LeR class is the main interface for simulating unlensed and lensed GW events and calculating detection rates. By default, it uses:
Event type: BBH (Binary Black Hole).
Lens galaxy model: EPL+Shear (Elliptical Power Law galaxy with external shears).
Detectors: H1, L1, V1 (LIGO Hanford, LIGO Livingston, Virgo) with O4 design sensitivity.
All outputs will be saved in the ./ler_data directory by default.
[6]:
# Import LeR
from ler import LeR
# Initialize LeR with default settings
ler = LeR()
Initializing LeR class...
Initializing LensGalaxyParameterDistribution class...
Initializing OpticalDepth class
comoving_distance interpolator will be loaded from ./interpolator_json/comoving_distance/comoving_distance_0.json
angular_diameter_distance interpolator will be loaded from ./interpolator_json/angular_diameter_distance/angular_diameter_distance_0.json
angular_diameter_distance interpolator will be loaded from ./interpolator_json/angular_diameter_distance/angular_diameter_distance_0.json
differential_comoving_volume interpolator will be loaded from ./interpolator_json/differential_comoving_volume/differential_comoving_volume_0.json
using ler available velocity dispersion function : velocity_dispersion_ewoud
velocity_dispersion_ewoud interpolator will be loaded from ./interpolator_json/velocity_dispersion/velocity_dispersion_ewoud_0.json
using ler available axis_ratio function : axis_ratio_rayleigh
axis_ratio_rayleigh interpolator will be loaded from ./interpolator_json/axis_ratio/axis_ratio_rayleigh_1.json
using ler available axis_rotation_angle function : axis_rotation_angle_uniform
using ler available density_profile_slope function : density_profile_slope_normal
using ler available external_shear function : external_shear_normal
Cross section interpolation data points loaded from ./interpolator_json/cross_section_function/cross_section_function_0.json
using ler available lens_redshift function : lens_redshift_strongly_lensed_numerical
Numerically solving the lens redshift distribution...
lens_redshift_strongly_lensed_numerical_epl_shear_galaxy interpolator will be loaded from ./interpolator_json/lens_redshift/lens_redshift_strongly_lensed_numerical_epl_shear_galaxy_1.json
lens_redshift_intrinsic interpolator will be loaded from ./interpolator_json/lens_redshift_intrinsic/lens_redshift_intrinsic_1.json
using ler available density_profile_slope_sl function : density_profile_slope_normal
using ler available external_shear_sl function : external_shear_normal
using ler available optical depth function : optical_depth_numerical
optical_depth_numerical interpolator will be loaded from ./interpolator_json/optical_depth/optical_depth_numerical_1.json
Initializing CBCSourceRedshiftDistribution class...
luminosity_distance interpolator will be loaded from ./interpolator_json/luminosity_distance/luminosity_distance_0.json
differential_comoving_volume interpolator will be loaded from ./interpolator_json/differential_comoving_volume/differential_comoving_volume_0.json
using ler available merger rate density model: merger_rate_density_madau_dickinson_belczynski_ng
merger_rate_density_madau_dickinson_belczynski_ng interpolator will be loaded from ./interpolator_json/merger_rate_density/merger_rate_density_madau_dickinson_belczynski_ng_0.json
merger_rate_density_detector_frame interpolator will be loaded from ./interpolator_json/merger_rate_density/merger_rate_density_detector_frame_1.json
source_redshift interpolator will be loaded from ./interpolator_json/source_redshift/source_redshift_0.json
Initializing CBCSourceParameterDistribution class...
using ler available zs function : source_redshift
using ler available source_frame_masses function : binary_masses_BBH_powerlaw_gaussian
using ler available geocent_time function : sampler_uniform
using ler available ra function : sampler_uniform
using ler available dec function : sampler_cosine
using ler available phase function : sampler_uniform
using ler available psi function : sampler_uniform
using ler available theta_jn function : sampler_sine
using ler available a_1 function : sampler_uniform
using ler available a_2 function : sampler_uniform
Faster, njitted and importance sampling based lens parameter sampler will be used.
Initializing GWSNR class...
psds not given. Choosing bilby's default psds
Interpolator will be loaded for L1 detector from ./interpolator_json/L1/partialSNR_dict_0.json
Interpolator will be loaded for H1 detector from ./interpolator_json/H1/partialSNR_dict_0.json
Interpolator will be loaded for V1 detector from ./interpolator_json/V1/partialSNR_dict_0.json
Chosen GWSNR initialization parameters:
npool: 4
snr type: interpolation_aligned_spins
waveform approximant: IMRPhenomD
sampling frequency: 2048.0
minimum frequency (fmin): 20.0
reference frequency (f_ref): 20.0
mtot=mass1+mass2
min(mtot): 9.96
max(mtot) (with the given fmin=20.0): 500.0
detectors: ['L1', 'H1', 'V1']
psds: [[array([ 10.21659, 10.23975, 10.26296, ..., 4972.81 ,
4984.081 , 4995.378 ], shape=(2736,)), array([4.43925574e-41, 4.22777986e-41, 4.02102594e-41, ...,
6.51153524e-46, 6.43165104e-46, 6.55252996e-46],
shape=(2736,)), <scipy.interpolate._interpolate.interp1d object at 0x304993ab0>], [array([ 10.21659, 10.23975, 10.26296, ..., 4972.81 ,
4984.081 , 4995.378 ], shape=(2736,)), array([4.43925574e-41, 4.22777986e-41, 4.02102594e-41, ...,
6.51153524e-46, 6.43165104e-46, 6.55252996e-46],
shape=(2736,)), <scipy.interpolate._interpolate.interp1d object at 0x3040aa5c0>], [array([ 10. , 10.02306 , 10.046173, ...,
9954.0389 , 9976.993 , 10000. ], shape=(3000,)), array([1.22674387e-42, 1.20400299e-42, 1.18169466e-42, ...,
1.51304203e-43, 1.52010157e-43, 1.52719372e-43],
shape=(3000,)), <scipy.interpolate._interpolate.interp1d object at 0x304160cc0>]]
#-------------------------------------
# LeR initialization input arguments:
#-------------------------------------
# LeR set up input arguments:
npool = 4,
z_min = 0.0,
z_max = 10.0,
event_type = 'BBH',
lens_type = 'epl_shear_galaxy',
cosmology = LambdaCDM(H0=70.0 km / (Mpc s), Om0=0.3, Ode0=0.7, Tcmb0=0.0 K, Neff=3.04, m_nu=None, Ob0=0.0),
pdet_finder = <bound method GWSNR.pdet of <gwsnr.core.gwsnr.GWSNR object at 0x313d6be80>>,
json_file_names = dict(
ler_params = 'ler_params.json',
unlensed_param = 'unlensed_param.json',
unlensed_param_detectable = 'unlensed_param_detectable.json',
lensed_param = 'lensed_param.json',
lensed_param_detectable = 'lensed_param_detectable.json',
),
interpolator_directory = './interpolator_json',
ler_directory = './ler_data',
create_new_interpolator = dict(
merger_rate_density = {'create_new': False, 'resolution': 500},
redshift_distribution = {'create_new': False, 'resolution': 500},
luminosity_distance = {'create_new': False, 'resolution': 500},
differential_comoving_volume = {'create_new': False, 'resolution': 500},
source_frame_masses = {'create_new': False, 'resolution': 500},
geocent_time = {'create_new': False, 'resolution': 500},
ra = {'create_new': False, 'resolution': 500},
dec = {'create_new': False, 'resolution': 500},
phase = {'create_new': False, 'resolution': 500},
psi = {'create_new': False, 'resolution': 500},
theta_jn = {'create_new': False, 'resolution': 500},
a_1 = {'create_new': False, 'resolution': 500},
a_2 = {'create_new': False, 'resolution': 500},
tilt_1 = {'create_new': False, 'resolution': 500},
tilt_2 = {'create_new': False, 'resolution': 500},
phi_12 = {'create_new': False, 'resolution': 500},
phi_jl = {'create_new': False, 'resolution': 500},
velocity_dispersion = {'create_new': False, 'resolution': 500, 'zl_resolution': 48},
axis_ratio = {'create_new': False, 'resolution': 500, 'sigma_resolution': 48},
lens_redshift = {'create_new': False, 'resolution': 48, 'zl_resolution': 48},
lens_redshift_intrinsic = {'create_new': False, 'resolution': 500},
optical_depth = {'create_new': False, 'resolution': 48},
comoving_distance = {'create_new': False, 'resolution': 500},
angular_diameter_distance = {'create_new': False, 'resolution': 500},
angular_diameter_distance_z1z2 = {'create_new': False, 'resolution': 500},
density_profile_slope = {'create_new': False, 'resolution': 100},
lens_parameters_kde_sl = {'create_new': False, 'resolution': 5000},
cross_section = {'create_new': False, 'resolution': [25, 25, 45, 15, 15]},
),
# LeR also takes other CBCSourceParameterDistribution class input arguments as kwargs, as follows:
source_priors = dict(
merger_rate_density = 'merger_rate_density_madau_dickinson_belczynski_ng',
zs = 'source_redshift',
source_frame_masses = 'binary_masses_BBH_powerlaw_gaussian',
geocent_time = 'sampler_uniform',
ra = 'sampler_uniform',
dec = 'sampler_cosine',
phase = 'sampler_uniform',
psi = 'sampler_uniform',
theta_jn = 'sampler_sine',
a_1 = 'sampler_uniform',
a_2 = 'sampler_uniform',
),
source_priors_params = dict(
merger_rate_density = {'R0': 1.9e-08, 'alpha_F': 2.57, 'beta_F': 5.83, 'c_F': 3.36},
zs = None,
source_frame_masses = {'mminbh': 4.98, 'mmaxbh': 112.5, 'alpha': 3.78, 'mu_g': 32.27, 'sigma_g': 3.88, 'lambda_peak': 0.03, 'delta_m': 4.8, 'beta': 0.81},
geocent_time = {'xmin': 1238166018, 'xmax': 1269702018},
ra = {'xmin': 0.0, 'xmax': 6.283185307179586},
dec = None,
phase = {'xmin': 0.0, 'xmax': 6.283185307179586},
psi = {'xmin': 0.0, 'xmax': 3.141592653589793},
theta_jn = None,
a_1 = {'xmin': -0.8, 'xmax': 0.8},
a_2 = {'xmin': -0.8, 'xmax': 0.8},
),
spin_zero = False,
spin_precession = False,
# LeR also takes other LensGalaxyParameterDistribution class input arguments as kwargs, as follows:
lens_functions = dict(
param_sampler_type = 'sample_all_routine_epl_shear_sl',
cross_section_based_sampler = 'importance_sampling_with_cross_section',
optical_depth = 'optical_depth_numerical',
cross_section = 'cross_section_epl_shear_interpolation',
),
lens_functions_params = dict(
param_sampler_type = None,
cross_section_based_sampler = {'n_prop': 200},
optical_depth = None,
cross_section = None,
),
lens_param_samplers = dict(
source_redshift_sl = 'strongly_lensed_source_redshifts',
lens_redshift = 'lens_redshift_strongly_lensed_numerical',
velocity_dispersion = 'velocity_dispersion_ewoud',
axis_ratio = 'axis_ratio_rayleigh',
axis_rotation_angle = 'axis_rotation_angle_uniform',
external_shear = 'external_shear_normal',
density_profile_slope = 'density_profile_slope_normal',
external_shear_sl = 'external_shear_normal',
density_profile_slope_sl = 'density_profile_slope_normal',
),
lens_param_samplers_params = dict(
source_redshift_sl = None,
lens_redshift = {'integration_size': 25000, 'use_multiprocessing': False},
velocity_dispersion = {'sigma_min': 100.0, 'sigma_max': 400.0, 'alpha': 0.94, 'beta': 1.85, 'phistar': np.float64(0.02099), 'sigmastar': 113.78, 'name': 'velocity_dispersion_ewoud'},
axis_ratio = {'q_min': 0.2, 'q_max': 1.0, 'name': 'axis_ratio_rayleigh'},
axis_rotation_angle = {'phi_min': 0.0, 'phi_max': 6.283185307179586, 'name': 'axis_rotation_angle_uniform'},
external_shear = {'mean': 0.0, 'std': 0.05, 'name': 'external_shear_normal'},
density_profile_slope = {'mean': 1.99, 'std': 0.149, 'name': 'density_profile_slope_normal'},
external_shear_sl = {'mean': 0.0, 'std': 0.05},
density_profile_slope_sl = {'mean': 2.078, 'std': 0.16},
),
# LeR also takes other ImageProperties class input arguments as kwargs, as follows:
n_min_images = 2,
n_max_images = 4,
time_window = 630720000,
lens_model_list = ['EPL_NUMBA', 'SHEAR'],
# LeR also takes other gwsnr.GWSNR input arguments as kwargs, as follows:
npool = 4,
snr_method = 'interpolation_aligned_spins',
snr_type = 'optimal_snr',
gwsnr_verbose = True,
multiprocessing_verbose = True,
pdet_kwargs = dict(
snr_th = 10.0,
snr_th_net = 10.0,
pdet_type = 'boolean',
distribution_type = 'noncentral_chi2',
include_optimal_snr = False,
include_observed_snr = False,
),
mtot_min = 9.96,
mtot_max = 500.0,
ratio_min = 0.1,
ratio_max = 1.0,
spin_max = 0.99,
mtot_resolution = 200,
ratio_resolution = 20,
spin_resolution = 10,
batch_size_interpolation = 1000000,
interpolator_dir = './interpolator_json',
create_new_interpolator = False,
sampling_frequency = 2048.0,
waveform_approximant = 'IMRPhenomD',
frequency_domain_source_model = 'lal_binary_black_hole',
minimum_frequency = 20.0,
reference_frequency = None,
duration_max = None,
duration_min = None,
fixed_duration = None,
mtot_cut = False,
psds = None,
ifos = None,
noise_realization = None,
ann_path_dict = None,
snr_recalculation = False,
snr_recalculation_range = [6, 14],
snr_recalculation_waveform_approximant = 'IMRPhenomXPHM',
psds_list = [[array([ 10.21659, 10.23975, 10.26296, ..., 4972.81 ,
4984.081 , 4995.378 ], shape=(2736,)), array([4.43925574e-41, 4.22777986e-41, 4.02102594e-41, ...,
6.51153524e-46, 6.43165104e-46, 6.55252996e-46],
shape=(2736,)), <scipy.interpolate._interpolate.interp1d object at 0x304993ab0>], [array([ 10.21659, 10.23975, 10.26296, ..., 4972.81 ,
4984.081 , 4995.378 ], shape=(2736,)), array([4.43925574e-41, 4.22777986e-41, 4.02102594e-41, ...,
6.51153524e-46, 6.43165104e-46, 6.55252996e-46],
shape=(2736,)), <scipy.interpolate._interpolate.interp1d object at 0x3040aa5c0>], [array([ 10. , 10.02306 , 10.046173, ...,
9954.0389 , 9976.993 , 10000. ], shape=(3000,)), array([1.22674387e-42, 1.20400299e-42, 1.18169466e-42, ...,
1.51304203e-43, 1.52010157e-43, 1.52719372e-43],
shape=(3000,)), <scipy.interpolate._interpolate.interp1d object at 0x304160cc0>]],
# To print all initialization input arguments, use:
ler._print_all_init_args()
2. Basic GW Event Simulation and Rate Calculation
This section demonstrates how to simulate unlensed binary black hole mergers and calculate their detection rates.
2.1 Simulate Unlensed GW Population
Generate a population of unlensed Compact Binary Coalescence (CBC) events. This step: - Samples intrinsic (masses and spins) and extrinsic (redshift, sky location, inclination angle, etc.) GW parameters from initialized priors. - Calculates the probability of detection (Pdet) for each event based on the detector network sensitivity. - Stores the output in ./ler_data/unlensed_param.json.
Note: For realistic results, use size >= 1,000,000 with batch_size = 100,000
[12]:
# Simulate 100,000 unlensed GW events in batches of 50,000
# use 'print(ler.unlensed_cbc_statistics.__doc__)' to see all the input args
unlensed_param = ler.unlensed_cbc_statistics(size=100000, resume=False)
print(f"\nTotal unlensed events simulated: {len(unlensed_param['zs'])}")
print(f"Sampled source redshift values (first 5): {unlensed_param['zs'][:5]}")
unlensed params will be stored in ./ler_data/unlensed_param.json
removing ./ler_data/unlensed_param.json if it exists
Batch no. 1
sampling gw source params...
calculating pdet...
Batch no. 2
sampling gw source params...
calculating pdet...
saving all unlensed parameters in ./ler_data/unlensed_param.json
Total unlensed events simulated: 100000
Sampled source redshift values (first 5): [2.38659538 1.87375022 3.05753123 2.03724768 3.06881079]
2.2 Calculate Unlensed Detection Rates
Select detectable events and calculate the detection rate. This function: - Filters events using a Pdet threshold. By default, Pdet is based on observed detector network SNR > 10 (gwsnr’s default), where the SNR follows a non-central chi-squared distribution centered at the optimal SNR (Essick et al. 2023). - Returns the rate in detectable events per year. - Saves detectable events to ./ler_data/unlensed_param_detectable.json.
[15]:
# Calculate the detection rate and extract detectable unlensed events
# use 'print(ler.unlensed_rate.__doc__)' to see all the input args
rate_unlensed, unlensed_param_detectable = ler.unlensed_rate()
print(f"\n=== Unlensed Detection Rate Summary ===")
print(f"Detectable event rate: {rate_unlensed:.2e} events per year")
print(f"Total event rate: {ler.normalization_pdf_z:.2e} events per year")
print(f"Percentage fraction of the detectable events: {rate_unlensed/ler.normalization_pdf_z*100:.2e}%")
Getting unlensed_param from json file ./ler_data/unlensed_param.json...
total unlensed rate (yr^-1): 294.89984863034954
number of simulated unlensed detectable events: 322
number of simulated all unlensed events: 100000
storing detectable params in ./ler_data/unlensed_param_detectable.json
=== Unlensed Detection Rate Summary ===
Detectable event rate: 2.95e+02 events per year
Total event rate: 9.16e+04 events per year
Percentage fraction of the detectable events: 3.22e-01%
2.3 Inspect Generated Unlensed Parameters
View the available parameters in the generated event population.
[8]:
# List all parameters available in the detectable event population
print("Detectable unlensed event parameters:")
print(list(unlensed_param_detectable.keys()))
print("\nExample values for mass_1_source (first 5 events):")
print(unlensed_param_detectable['mass_1_source'][:5])
Detectable unlensed event parameters:
['zs', 'geocent_time', 'ra', 'dec', 'phase', 'psi', 'theta_jn', 'a_1', 'a_2', 'luminosity_distance', 'mass_1_source', 'mass_2_source', 'mass_1', 'mass_2', 'pdet_L1', 'pdet_H1', 'pdet_V1', 'pdet_net']
Example values for mass_1_source (first 5 events):
[13.76440519 36.07826189 29.91132175 36.52243878 77.44086226]
3. Lensed GW Population Simulation and Rates
This part demonstrates the simulation of lensed GW events and rate calculation. Lensing includes additional parameters such as lens galaxy properties (lens redshift, velocity dispersion, etc.) and image characteristics (magnification, time delays, etc.).
3.1 Simulate Lensed GW Population
Generate a population of lensed CBC events including lens galaxy properties and image parameters: - Source parameters (redshift, masses, spins) - Lens parameters (redshift, Einstein radius, ellipticity, shear components) - Image parameters (magnifications, time delays)
This step stores output in ./ler_data/lensed_param.json
Note: The simulation includes: - Lens galaxy population sampling - Selection based on lensing cross-section - Lens equation solving for multiple image generation - Pdet calculation for each image
[16]:
# Simulate 100,000 unlensed GW events in batches of 50,000
# use 'ler.lensed_cbc_statistics.__doc__' to see all the input args
lensed_param = ler.lensed_cbc_statistics(size=100000, resume=False)
print(f"\nTotal lensed events simulated: {len(lensed_param['zs'])}")
print(f"Sampled source redshift values (first 5): {lensed_param['zs'][:5]}")
lensed params will be stored in ./ler_data/lensed_param.json
removing ./ler_data/lensed_param.json if it exists
Batch no. 1
sampling lensed params...
sampling lens parameters with sample_all_routine_epl_shear_sl...
solving lens equations...
100%|███████████████████████████████████████████████████████| 50000/50000 [00:11<00:00, 4524.78it/s]
calculating pdet...
Batch no. 2
sampling lensed params...
sampling lens parameters with sample_all_routine_epl_shear_sl...
solving lens equations...
100%|███████████████████████████████████████████████████████| 50000/50000 [00:10<00:00, 4563.30it/s]
calculating pdet...
saving all lensed parameters in ./ler_data/lensed_param.json
Total lensed events simulated: 100000
Sampled source redshift values (first 5): [2.96150209 3.62643698 2.22925302 6.69295117 3.6017759 ]
3.2 Calculate Lensed Detection Rates
Calculate the lensed detection rate, with each event requiring at least two detectable images for a valid detection.
The detection rate is stored in
./ler_data/lensed_param_detectable.json.
[20]:
# Calculate the detection rate for lensed events
# use 'print(ler.lensed_rate.__doc__)' to see all the input args
rate_lensed, lensed_param_detectable = ler.lensed_rate()
print(f"\n=== Lensed Detection Rate Summary ===")
print(f"Detectable event rate: {rate_lensed:.2e} events per year")
print(f"Total event rate: {ler.normalization_pdf_z_lensed:.2e} events per year")
print(f"Percentage fraction of the detectable events: {rate_lensed/ler.normalization_pdf_z_lensed*100:.2e}%")
Getting lensed_param from json file ./ler_data/lensed_param.json...
total lensed rate (yr^-1): 0.10082199896856103
number of simulated lensed detectable events: 89
number of simulated all lensed events: 100000
storing detectable params in ./ler_data/lensed_param_detectable.json
=== Lensed Detection Rate Summary ===
Detectable event rate: 1.01e-01 events per year
Total event rate: 1.13e+02 events per year
Percentage fraction of the detectable events: 8.90e-02%
3.3 Inspect Generated Lensed Parameters
[9]:
# List all parameters available in the detectable lensed event population
print("Detectable lensed event parameters:")
print(list(lensed_param_detectable.keys()))
print("\nLens-specific parameters include:")
lens_params = ['zl', 'sigma', 'theta_E', 'q', 'phi', 'gamma', 'gamma1', 'gamma2']
for param in lens_params:
if param in lensed_param_detectable:
print(f" {param}: {lensed_param_detectable[param][:3]}")
Detectable lensed event parameters:
['zl', 'zs', 'sigma', 'theta_E', 'q', 'phi', 'gamma', 'gamma1', 'gamma2', 'geocent_time', 'ra', 'dec', 'phase', 'psi', 'theta_jn', 'a_1', 'a_2', 'luminosity_distance', 'mass_1_source', 'mass_2_source', 'mass_1', 'mass_2', 'x0_image_positions', 'x1_image_positions', 'magnifications', 'time_delays', 'image_type', 'n_images', 'x_source', 'y_source', 'effective_luminosity_distance', 'effective_geocent_time', 'effective_phase', 'pdet_net', 'L1', 'H1', 'V1']
Lens-specific parameters include:
zl: [0.54337231 1.87777572 0.64609008]
sigma: [181.65920777 128.87706255 196.8531673 ]
theta_E: [2.03764065e-06 6.01901250e-07 2.62310898e-06]
q: [0.66697535 0.4918578 0.77095448]
phi: [3.25526037 6.01826202 3.75664554]
gamma: [1.88747336 2.18793729 2.0281945 ]
gamma1: [ 0.05469759 -0.06132396 -0.02247388]
gamma2: [-0.1018271 0.0365925 -0.02447439]
4. Rate Comparison and Visualization
Compare the lensed and unlensed detection rates, access saved data files, and visualize parameter distributions.
4.1 Compare Lensed vs Unlensed Rates
[21]:
# Compare lensed and unlensed rates
rate_ratio = ler.rate_ratio()
print(f"\n=== Rate Comparison ===")
print(f"Unlensed detection rate: {rate_unlensed:.4e} events/yr")
print(f"Lensed detection rate: {rate_lensed:.4e} events/yr")
print(f"Rate ratio (lensed/unlensed): {rate_ratio:.4f}")
unlensed_rate: 294.89984863034954
lensed_rate: 0.10082199896856103
ratio: 2924.9553832225356
=== Rate Comparison ===
Unlensed detection rate: 2.9490e+02 events/yr
Lensed detection rate: 1.0082e-01 events/yr
Rate ratio (lensed/unlensed): 2924.9554
Bonus: The following command simultaneously selects detectable events, calculates rates, and compares unlensed and lensed events.
unlensed_rate, lensed_rate, rate_ratio, unlensed_param_detectable, lensed_param_detectable = ler.rate_comparison_with_rate_calculation()
4.2 Access Saved Data Files
All simulation results are saved in JSON files for future reference and analysis. View the saved file locations and load parameters from disk.
[8]:
# View saved file locations and names
print(f"Output directory: {ler.ler_directory}")
print(f"\nSaved JSON files:")
for key, filename in ler.json_file_names.items():
print(f" {key}: {filename}")
Output directory: ./ler_data
Saved JSON files:
ler_params: ler_params.json
unlensed_param: unlensed_param.json
unlensed_param_detectable: unlensed_param_detectable.json
lensed_param: lensed_param.json
lensed_param_detectable: lensed_param_detectable.json
4.3 Load and Examine Saved Parameters
Reload parameters from JSON files for further analysis. This example also hightlights some of the useful functions from ler.utils module.
[9]:
# Load parameters from saved JSON files
from ler.utils import get_param_from_json, load_json
# Load detectable parameters from files
unlensed_param_from_file = get_param_from_json(
ler.ler_directory + '/' + ler.json_file_names['unlensed_param_detectable']
)
lensed_param_from_file = get_param_from_json(
ler.ler_directory + '/' + ler.json_file_names['lensed_param_detectable']
)
print(f"Unlensed parameters loaded: {list(unlensed_param_from_file.keys())}")
print(f"Lensed parameters loaded: {list(lensed_param_from_file.keys())}")
# Load initialization parameters and results
ler_params = load_json(ler.ler_directory + '/ler_params.json')
print(f"\n=== Rates from saved file ===")
print(f"Detectable unlensed rate: {ler_params['detectable_unlensed_rate_per_year']:.4e} events/yr")
print(f"Detectable lensed rate: {ler_params['detectable_lensed_rate_per_year']:.4e} events/yr")
print(f"Rate ratio: {ler_params['rate_ratio']:.4f}")
Unlensed parameters loaded: ['zs', 'geocent_time', 'ra', 'dec', 'phase', 'psi', 'theta_jn', 'a_1', 'a_2', 'luminosity_distance', 'mass_1_source', 'mass_2_source', 'mass_1', 'mass_2', 'pdet_L1', 'pdet_H1', 'pdet_V1', 'pdet_net']
Lensed parameters loaded: ['zl', 'zs', 'sigma', 'theta_E', 'q', 'phi', 'gamma', 'gamma1', 'gamma2', 'geocent_time', 'ra', 'dec', 'phase', 'psi', 'theta_jn', 'a_1', 'a_2', 'luminosity_distance', 'mass_1_source', 'mass_2_source', 'mass_1', 'mass_2', 'x0_image_positions', 'x1_image_positions', 'magnifications', 'time_delays', 'image_type', 'n_images', 'x_source', 'y_source', 'effective_luminosity_distance', 'effective_geocent_time', 'effective_phase', 'pdet_net', 'L1', 'H1', 'V1']
=== Rates from saved file ===
Detectable unlensed rate: 2.9582e+02 events/yr
Detectable lensed rate: 8.4985e-02 events/yr
Rate ratio: 3480.8049
4.4 Visualize Parameter Distributions
Create KDE (Kernel Density Estimation) plots to compare the distributions of source redshift for lensed and unlensed populations.
[22]:
# Visualize parameter distributions
import matplotlib.pyplot as plt
from ler.utils import plots as lerplt
# Plot source redshift distributions
plt.figure(figsize=(6, 4))
# Unlensed populations
lerplt.param_plot(
param_name='zs',
param_dict='ler_data/unlensed_param_detectable.json', # can also provide the dict directly
plot_label='unlensed-detectable',
histogram=False,
kde=True,
kde_bandwidth=0.5,
)
lerplt.param_plot(
param_name='zs',
param_dict='ler_data/unlensed_param.json',
plot_label='unlensed-all',
histogram=False,
kde=True,
)
# Lensed populations
lerplt.param_plot(
param_name='zs',
param_dict='ler_data/lensed_param_detectable.json',
plot_label='lensed-detectable',
histogram=False,
kde=True,
kde_bandwidth=0.5,
)
lerplt.param_plot(
param_name='zs',
param_dict='ler_data/lensed_param.json',
plot_label='lensed-all',
histogram=False,
kde=True,
)
plt.xlim(0.001, 8)
plt.grid(alpha=0.4)
plt.xlabel('Source redshift (zs)', fontsize=12)
plt.ylabel('Probability Density', fontsize=12)
plt.title('Source Redshift Distribution: Lensed vs Unlensed', fontsize=14, fontweight='bold')
plt.legend(fontsize=10, loc='upper right')
plt.tight_layout()
plt.show()
getting gw_params from json file ler_data/unlensed_param_detectable.json...
getting gw_params from json file ler_data/unlensed_param.json...
getting gw_params from json file ler_data/lensed_param_detectable.json...
getting gw_params from json file ler_data/lensed_param.json...
5. Available Internal Functions and Their Usage
5.1 Explore Functions and Parameters
Inspect the internal functions and parameters used for LeR initialization.
[11]:
# List all available GW prior are in ler.available_gw_prior
print("# ----------------------------------------------------")
print("# GW prior sampler functions and it's input arguments:")
print("# ----------------------------------------------------\n")
for key, value in ler.available_gw_prior.items():
print(f"{key} = dict(")
if isinstance(value, dict):
for k, v in value.items():
print(f" {k} = {v},")
else:
print(f" {value},")
print(")")
# list all available lens prior are in ler.available_lens_samplers
print("\n# ----------------------------------------------------")
print("# Lens prior sampler functions and it's input arguments:")
print("# ----------------------------------------------------\n")
for key, value in ler.available_lens_samplers.items():
print(f"{key} = dict(")
if isinstance(value, dict):
for k, v in value.items():
print(f" {k} = {v},")
else:
print(f" {value},")
print(")")
# list all available lens functions are in ler.available_lens_functions
print("\n# ----------------------------------------------------")
print("# Other lens related functions and it's input arguments:")
print("# ----------------------------------------------------\n")
for key, value in ler.available_lens_functions.items():
print(f"{key} = dict(")
if isinstance(value, dict):
for k, v in value.items():
print(f" {k} = {v},")
else:
print(f" {value},")
print(")")
# ----------------------------------------------------
# GW prior sampler functions and it's input arguments:
# ----------------------------------------------------
merger_rate_density = dict(
merger_rate_density_bbh_oguri2018 = {'R0': 1.9e-08, 'b2': 1.6, 'b3': 2.1, 'b4': 30},
merger_rate_density_madau_dickinson2014 = {'R0': 1.9e-08, 'a': 0.015, 'b': 2.7, 'c': 2.9, 'd': 5.6},
merger_rate_density_madau_dickinson_belczynski_ng = {'R0': 1.9e-08, 'alpha_F': 2.57, 'beta_F': 5.83, 'c_F': 3.36},
sfr_with_time_delay = {'R0': 1.9e-08, 'a': 0.01, 'b': 2.6, 'c': 3.2, 'd': 6.2, 'td_min': 0.01, 'td_max': 10.0},
merger_rate_density_bbh_popIII_ken2022 = {'R0': 1.92e-08, 'aIII': 0.66, 'bIII': 0.3, 'zIII': 11.6},
merger_rate_density_bbh_primordial_ken2022 = {'R0': 4.4e-11, 't0': 13.786885302009708},
)
zs = dict(
source_redshift = None,
)
source_frame_masses = dict(
binary_masses_BBH_powerlaw_gaussian = {'mminbh': 4.98, 'mmaxbh': 112.5, 'alpha': 3.78, 'mu_g': 32.27, 'sigma_g': 3.88, 'lambda_peak': 0.03, 'delta_m': 4.8, 'beta': 0.81},
binary_masses_BBH_popIII_lognormal = {'m_min': 5.0, 'm_max': 150.0, 'Mc': 30.0, 'sigma': 0.3},
binary_masses_BBH_primordial_lognormal = {'m_min': 1.0, 'm_max': 100.0, 'Mc': 20.0, 'sigma': 0.3},
binary_masses_NSBH_broken_powerlaw = {'mminbh': 26, 'mmaxbh': 125, 'alpha_1': 6.75, 'alpha_2': 6.75, 'b': 0.5, 'delta_m': 5, 'mminns': 1.0, 'mmaxns': 3.0, 'alphans': 0.0},
binary_masses_uniform = {'m_min': 1.0, 'm_max': 3.0},
binary_masses_BNS_bimodal = {'w': 0.643, 'muL': 1.352, 'sigmaL': 0.08, 'muR': 1.88, 'sigmaR': 0.3, 'mmin': 1.0, 'mmax': 2.3},
)
a_1 = dict(
constant_values_n_size = {'value': 0.0},
sampler_uniform = {'xmin': -0.8, 'xmax': 0.8},
)
a_2 = dict(
constant_values_n_size = {'value': 0.0},
sampler_uniform = {'xmin': -0.8, 'xmax': 0.8},
)
tilt_1 = dict(
constant_values_n_size = {'value': 0.0},
sampler_sine = None,
)
tilt_2 = dict(
constant_values_n_size = {'value': 0.0},
sampler_sine = None,
)
phi_12 = dict(
constant_values_n_size = {'value': 0.0},
sampler_uniform = {'xmin': 0.0, 'xmax': 6.283185307179586},
)
phi_jl = dict(
constant_values_n_size = {'value': 0.0},
sampler_uniform = {'xmin': 0.0, 'xmax': 6.283185307179586},
)
geocent_time = dict(
sampler_uniform = {'xmin': 1238166018, 'xmax': 1269723618.0},
constant_values_n_size = {'value': 1238166018},
)
ra = dict(
sampler_uniform = {'xmin': 0.0, 'xmax': 6.283185307179586},
constant_values_n_size = {'value': 0.0},
)
dec = dict(
sampler_cosine = None,
constant_values_n_size = {'value': 0.0},
sampler_uniform = {'xmin': -1.5707963267948966, 'xmax': 1.5707963267948966},
)
phase = dict(
sampler_uniform = {'xmin': 0.0, 'xmax': 6.283185307179586},
constant_values_n_size = {'value': 0.0},
)
psi = dict(
sampler_uniform = {'xmin': 0.0, 'xmax': 3.141592653589793},
constant_values_n_size = {'value': 0.0},
)
theta_jn = dict(
sampler_sine = None,
constant_values_n_size = {'value': 0.0},
sampler_uniform = {'xmin': 0.0, 'xmax': 3.141592653589793},
)
# ----------------------------------------------------
# Lens prior sampler functions and it's input arguments:
# ----------------------------------------------------
source_redshift_sl = dict(
strongly_lensed_source_redshifts = None,
)
lens_redshift = dict(
lens_redshift_strongly_lensed_sis_haris = None,
lens_redshift_strongly_lensed_numerical = {'integration_size': 25000, 'use_multiprocessing': False},
lens_redshift_strongly_lensed_hemanta = None,
)
velocity_dispersion = dict(
velocity_dispersion_gengamma = {'sigma_min': 100.0, 'sigma_max': 400.0, 'alpha': 0.94, 'beta': 1.85, 'phistar': np.float64(0.02099), 'sigmastar': 113.78},
velocity_dispersion_choi = {'sigma_min': 100.0, 'sigma_max': 400.0, 'alpha': 2.32, 'beta': 2.67, 'phistar': np.float64(0.0027439999999999995), 'sigmastar': 161.0},
velocity_dispersion_bernardi = {'sigma_min': 100.0, 'sigma_max': 400.0, 'alpha': 0.94, 'beta': 1.85, 'phistar': np.float64(0.02099), 'sigmastar': 113.78},
velocity_dispersion_ewoud = {'sigma_min': 100.0, 'sigma_max': 400.0, 'alpha': 0.94, 'beta': 1.85, 'phistar': np.float64(0.02099), 'sigmastar': 113.78},
)
axis_ratio = dict(
axis_ratio_rayleigh = {'q_min': 0.2, 'q_max': 1.0},
axis_ratio_padilla_strauss = {'q_min': 0.2, 'q_max': 1.0},
axis_ratio_uniform = {'q_min': 0.2, 'q_max': 1.0},
)
axis_rotation_angle = dict(
axis_rotation_angle_uniform = {'phi_min': 0.0, 'phi_max': 6.283185307179586},
)
external_shear = dict(
external_shear_normal = {'mean': 0.0, 'std': 0.05},
)
external_shear_sl = dict(
external_shear_normal = {'mean': 0.0, 'std': 0.05},
external_shear_sl_numerical_hemanta = {'external_shear_normal': {'mean': 0.0, 'std': 0.05}},
)
density_profile_slope = dict(
density_profile_slope_normal = {'mean': 1.99, 'std': 0.149},
)
density_profile_slope_sl = dict(
density_profile_slope_normal = {'mean': 2.091, 'std': 0.133},
density_profile_slope_sl_numerical_hemanta = {'density_profile_slope_normal': {'mean': 1.99, 'std': 0.149}},
)
source_parameters = dict(
sample_gw_parameters = None,
)
# ----------------------------------------------------
# Other lens related functions and it's input arguments:
# ----------------------------------------------------
cross_section_based_sampler = dict(
rejection_sampling_with_cross_section = {'safety_factor': 1.2},
importance_sampling_with_cross_section = {'n_prop': 200},
)
optical_depth = dict(
optical_depth_sis_analytic = None,
optical_depth_epl_shear_hemanta = None,
optical_depth_numerical = None,
)
param_sampler_type = dict(
sample_all_routine_epl_shear_sl = None,
)
cross_section = dict(
cross_section_sie_feixu = None,
cross_section_sis = None,
cross_section_epl_shear_numerical = None,
cross_section_epl_shear_interpolation = None,
)
5.2 Short Example on Using Internal Functions
Following is a code snippet to initialize LeR with internal functions. Refer to the dedicated LeR with custom functions example for more details.
ler = LeR(
source_priors=dict(
merger_rate_density='sfr_madau_fragos2017',
),
source_priors_params=dict(
merger_rate_density={'R0': 2.39e-08, 'b2': 1.6, 'b3': 2.1, 'b4': 30}
),
lens_param_samplers=dict(
velocity_dispersion='velocity_dispersion_bernardi',
)
lens_param_samplers_params={'sigma_min': 100.0, 'sigma_max': 400.0, 'alpha': 0.94, 'beta': 1.85, 'phistar': np.float64(0.02099), 'sigmastar': 113.78}
)
6. Summary
This notebook demonstrates gravitational wave detection rate simulations using the LeR framework for both unlensed and strongly lensed compact binary coalescence (CBC) events.
Workflow: Initialize LeR with default BBH sources, EPL+Shear lens model, and O4-sensitivity detectors (H1, L1, V1) → simulate 100,000 unlensed CBC events and compute detection rates → sample lens galaxy properties, solve lens equations, and compute lensed rates (requiring ≥2 detectable images) → compare rates and visualize redshift distributions.
Sampled Parameters: - GW source: redshift, masses, spins, sky location, orientation, luminosity distance, detection probability - Lens: redshift, velocity dispersion, Einstein radius, axis ratio, density slope, shear - Images: magnifications, time delays, image types, effective luminosity distances
Output Files (saved to ./ler_data/): unlensed_param.json, unlensed_param_detectable.json, lensed_param.json, lensed_param_detectable.json, ler_params.json
Customization: Access built-in samplers via ler.available_gw_prior, ler.available_lens_samplers, and ler.available_lens_functions. See the LeR with custom functions example for advanced configurations.