LeR with Custom Functions and Parameters

This notebook is created by Phurailatpam Hemantakumar

Documentation

This notebook demonstrates how to customize the LeR class by defining custom source priors, lens parameter samplers, and detection criteria. We’ll simulate both unlensed and lensed binary neutron star (BNS) events with fully custom configurations, then compare the results with default models.

Table of Contents

  1. LeR initialization

    • 1.1 Code snippet for initialization with default input arguments

    • 1.2 Initialize LeR with default settings

  2. Examine Available Prior Functions

    • 2.1 Accessing functions as Class Attributes (Unlensed)

    • 2.2 Testing source related parameter prior functions

    • 2.3 Accessing functions as Class Attributes (Lensed)

    • 2.4 Accessing functions from the ler.gw_source_population module (Unlensed)

    • 2.5 Accessing functions from the ler.lens_galaxy_population module (Lensed)

  3. Using Custom Functions in LeR Initialization

    • 3.1 Custom Event Type with non-spinning configuration

    • 3.2 Custom Merger Rate Density

    • 3.3 Custom Source Frame Masses

    • 3.4 Custom Lens Model

      • 3.4.1 Custom Velocity Dispersion

      • 3.4.2 Custom Axis Ratio

    • 3.5 Custom Detection Criteria

    • 3.6 LeR initialization with custom functions and parameters

    • 3.7 Simulate Unlensed Population

    • 3.8 Simulate Lensed Population

    • 3.9 Calculate Rates and Compare Results

  4. Compare Custom vs Default Models

    • 4.1 Mass Distribution Comparison

    • 4.2 Axis-Ratio Distribution Comparison

  5. Summary

1. LeR initialization

1.1. Code snippet for initialization with default input arguments

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

# Uncomment the below code if you need to change the default arguments.
ler = LeR(
    # LeR setup arguments
    npool=4, # number of processors to use
    z_min=0.0, # minimum redshift
    z_max=10.0, # maximum redshift
    event_type='BBH', # event type
    lens_type = 'epl_shear_galaxy', # lens type
    cosmology=cosmo, # cosmology
    pdet_finder=None,  # if None, the pdet_finder will be calculated using the gwsnr package.
    json_file_names=dict(
        ler_params="ler_params.json", # to store initialization parameters and important results
        unlensed_param="unlensed_param.json", # to store all unlensed events
        unlensed_param_detectable="unlensed_param_detectable.json", # to store only detectable unlensed events
        lensed_param="lensed_param.json", # to store all lensed events
        lensed_param_detectable="lensed_param_detectable.json"), # to store only detectable lensed events
    interpolator_directory='./interpolator_json', # directory to store the interpolator pickle files. 'ler' uses interpolation to get values of various functions to speed up the calculations (relying on numba njit).
    ler_directory='./ler_data', # directory to store all the outputs
    verbose=True, # if True, will print all information at initialization

    # CBCSourceParameterDistribution class arguments
    gw_functions = dict(
        merger_rate_density = 'merger_rate_density_madau_dickinson_belczynski_ng', # function to calculate the merger rate density at a given redshift
        param_sampler_type = 'gw_parameters_rvs', # function to sample the source parameters. This will be used to sample the source parameters for all events (lensed and unlensed)
    ),
    gw_functions_params = dict(
        merger_rate_density = dict(param_name='merger_rate_density', function_type='merger_rate_density_madau_dickinson_belczynski_ng', R0=1.9e-08, alpha_F=2.57, beta_F=5.83, c_F=3.36),
        param_sampler_type = None,
    ),
    gw_priors = dict(
        zs = 'merger_rate_density_based_source_redshift',
        mass_1_source = 'broken_powerlaw_plus_2peaks',
        mass_ratio = 'powerlaw_with_smoothing',
        mass_2_source = None,
        geocent_time = 'uniform',
        ra = 'uniform',
        dec = 'sampler_cosine',
        phase = 'uniform',
        psi = 'uniform',
        theta_jn = 'sampler_sine',
        a_1 = 'uniform',
        a_2 = 'uniform',
        tilt_1 = None,
        tilt_2 = None,
        phi_12 = None,
        phi_jl = None,
    ),
    gw_priors_params = dict(
        zs = None,
        mass_1_source = dict(param_name='mass_1_source', sampler_type='broken_powerlaw_plus_2peaks', lam_0=0.361, lam_1=0.586, mpp_1=9.764, sigpp_1=0.649, mpp_2=32.763, sigpp_2=3.918, mlow_1=5.059, delta_m_1=4.321, break_mass=35.622, alpha_1=1.728, alpha_2=4.512, mmax=300.0, normalization_size=500),
        mass_ratio = dict(param_name='mass_ratio', sampler_type='powerlaw_with_smoothing', q_min=0.01, q_max=1.0, mlow_2=3.551, mmax=300.0, beta=1.171, delta_m=4.91, mmin=3.551),
        mass_2_source = None,
        geocent_time = dict(param_name='geocent_time', sampler_type='uniform', x_min=1238166018, x_max=1269702018),
        ra = dict(param_name='ra', sampler_type='uniform', x_min=0.0, x_max=6.283185307179586),
        dec = dict(param_name='dec', sampler_type='sampler_cosine', dec_min=-1.5707963267948966, dec_max=1.5707963267948966),
        phase = dict(param_name='phase', sampler_type='uniform', x_min=0.0, x_max=6.283185307179586),
        psi = dict(param_name='psi', sampler_type='uniform', x_min=0.0, x_max=3.141592653589793),
        theta_jn = dict(param_name='theta_jn', sampler_type='sampler_sine', theta_jn_min=0.0, theta_jn_max=3.141592653589793),
        a_1 = dict(param_name='a_1', sampler_type='uniform', x_min=-1.0, x_max=1.0),
        a_2 = dict(param_name='a_2', sampler_type='uniform', x_min=-1.0, x_max=1.0),
        tilt_1 = None,
        tilt_2 = None,
        phi_12 = None,
        phi_jl = None,
    ),
    spin_zero = False,
    spin_precession = False,

    # LensGalaxyParameterDistribution class arguments
    lens_functions = dict(
        param_sampler_type = 'epl_shear_sl_parameters_rvs',
        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 = dict(n_prop_factor=10),
        optical_depth = None,
        cross_section = dict(num_th=500, maginf=-100.0),
    ),
    lens_param_samplers = dict(
        zs_sl = 'strongly_lensed_source_redshift',
        lens_redshift_sl = 'lens_redshift_strongly_lensed_numerical',
        lens_redshift = 'lens_redshift_intrinsic_numerical',
        velocity_dispersion = 'velocity_dispersion_ewoud',
        axis_ratio = 'rayleigh',
        axis_rotation_angle = 'uniform',
        external_shear1 = 'normal',
        external_shear2 = 'normal',
        density_profile_slope = 'normal',
    ),
    lens_param_samplers_params = dict(
        zs_sl = None,
        lens_redshift_sl = dict(integration_size=25000, use_multiprocessing=False),
        lens_redshift = None,
        velocity_dispersion = dict(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(q_min=0.2, q_max=1.0),
        axis_rotation_angle = dict(x_min=0.0, x_max=6.283185307179586),
        external_shear1 = dict(mu=0.0, sigma=0.05, x_min=-np.inf, x_max=np.inf),
        external_shear2 = dict(mu=0.0, sigma=0.05, x_min=-np.inf, x_max=np.inf),
        density_profile_slope = dict(mu=1.99, sigma=0.149, x_min=-np.inf, x_max=np.inf),
    ),

    # ImageProperties class arguments
    n_min_images = 2,
    n_max_images = 4,
    time_window = 365. * 24. * 3600. * 1.,  # 1 year
    include_effective_parameters = False,
    lens_model_list = ['EPL_NUMBA', 'SHEAR'],
    image_properties_function = 'image_properties_epl_shear_njit',
    multiprocessing_verbose = True,
    include_redundant_parameters = False,

    # gwsnr package arguments
    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 = 1.0,
    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',
    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,  # will consider the default bilby psds of 'L1', 'H1', 'V1'
    ifos = None,  # will consider the default bilby interferometer objects of 'L1', 'H1', 'V1'
    snr_recalculation = False,
    snr_recalculation_range = [6, 14],
    snr_recalculation_waveform_approximant = 'IMRPhenomXPHM',

    # common arguments, to generate interpolator
    create_new_interpolator = dict(
        velocity_dispersion = dict(create_new=False, resolution=200, zl_resolution=48),
        axis_ratio = dict(create_new=False, resolution=200, sigma_resolution=48),
        lens_redshift_sl = dict(create_new=False, resolution=16, zl_resolution=48),
        lens_redshift = dict(create_new=False, resolution=16, zl_resolution=48),
        optical_depth = dict(create_new=False, resolution=16),
        comoving_distance = dict(create_new=False, resolution=500),
        angular_diameter_distance = dict(create_new=False, resolution=500),
        angular_diameter_distance_z1z2 = dict(create_new=False, resolution=500),
        differential_comoving_volume = dict(create_new=False, resolution=500),
        zs_sl = dict(create_new=False, resolution=200),
        cross_section = dict(create_new=False, resolution=[25, 25, 45, 15, 15], spacing_config=dict(e1=dict(mode='two_sided_mixed_grid', power_law_part='lower', spacing_trend='increasing', power=2.5, value_transition_fraction=0.7, num_transition_fraction=0.9, auto_match_slope=True), e2=dict(mode='two_sided_mixed_grid', power_law_part='lower', spacing_trend='increasing', power=2.5, value_transition_fraction=0.7, num_transition_fraction=0.9, auto_match_slope=True), gamma=dict(mode='two_sided_mixed_grid', power_law_part='lower', spacing_trend='increasing', power=2.5, value_transition_fraction=0.7, num_transition_fraction=0.9, auto_match_slope=True), gamma1=dict(mode='two_sided_mixed_grid', power_law_part='lower', spacing_trend='increasing', power=2.5, value_transition_fraction=0.7, num_transition_fraction=0.9, auto_match_slope=True), gamma2=dict(mode='two_sided_mixed_grid', power_law_part='lower', spacing_trend='increasing', power=2.5, value_transition_fraction=0.7, num_transition_fraction=0.9, auto_match_slope=True))),
        merger_rate_density = dict(create_new=False, resolution=200),
        redshift_distribution = dict(create_new=False, resolution=200),
        luminosity_distance = dict(create_new=False, resolution=500),
        mass_1_source = dict(create_new=False, resolution=200),
        mass_ratio = dict(create_new=False, resolution=100, m1_resolution=200),
        mass_2_source = dict(create_new=False, resolution=200),
        geocent_time = dict(create_new=False, resolution=200),
        ra = dict(create_new=False, resolution=200),
        dec = dict(create_new=False, resolution=200),
        phase = dict(create_new=False, resolution=200),
        psi = dict(create_new=False, resolution=200),
        theta_jn = dict(create_new=False, resolution=200),
        a_1 = dict(create_new=False, resolution=200),
        a_2 = dict(create_new=False, resolution=200),
        tilt_1 = dict(create_new=False, resolution=200),
        tilt_2 = dict(create_new=False, resolution=200, tilt_1_resolution=100),
        phi_12 = dict(create_new=False, resolution=200),
        phi_jl = dict(create_new=False, resolution=200),
    ),
)

1.2. Initialize LeR with default settings

Set verbose=False to suppress lengthy output

[1]:

from ler import LeR
import numpy as np

# use this gwsnr's input args if you want SNR values in the output, besides the boolean detection probability values. This will increase the output dictionary size and the runtime of the code. For details, refer to https://gwsnr.hemantaph.com/examples/pdet_generation.html.:
# pdet_kwargs=dict(snr_th=10.0, snr_th_net=10.0, pdet_type="boolean", distribution_type="noncentral_chi2", include_optimal_snr=True, include_observed_snr=True)
ler = LeR(verbose=False)

2. Examine Available Prior Functions

There are two ways of accessing the built-in GW parameter prior functions and their default parameters.

2.1. Accessing functions as Class Attributes (Unlensed)

[2]:

# Display all available GW prior functions and sampler functions and their parameters
print("Built-in GW parameter functions and sampler functions and parameters:\n")
for func_name, func_params in ler.available_gw_functions.items():
    print(f"{func_name}:")

    if isinstance(func_params, dict):
        for param_name, param_value in func_params.items():
            print(f"  {param_name}: {param_value}")
    else:
        print(f"  {func_params}")
    print()

for func_name, func_params in ler.available_gw_prior.items():
    print(f"{func_name}:")

    if isinstance(func_params, dict):
        for param_name, param_value in func_params.items():
            print(f"  {param_name}: {param_value}")
    else:
        print(f"  {func_params}")
    print()

Built-in GW parameter functions and sampler functions and parameters:

merger_rate_density:
  merger_rate_density_madau_dickinson_belczynski_ng: {'param_name': 'merger_rate_density', 'function_type': 'merger_rate_density_madau_dickinson_belczynski_ng', 'R0': 1.9e-08, 'alpha_F': 2.57, 'beta_F': 5.83, 'c_F': 3.36}
  merger_rate_density_bbh_oguri2018: {'param_name': 'merger_rate_density', 'function_type': 'merger_rate_density_bbh_oguri2018', 'R0': 1.9e-08, 'b2': 1.6, 'b3': 2.1, 'b4': 30}
  merger_rate_density_madau_dickinson2014: {'param_name': 'merger_rate_density', 'function_type': 'merger_rate_density_madau_dickinson2014', 'R0': 8.9e-08, 'a': 0.015, 'b': 2.7, 'c': 2.9, 'd': 5.6}
  sfr_with_time_delay: {'param_name': 'merger_rate_density', 'function_type': '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: {'param_name': 'merger_rate_density', 'function_type': '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: {'param_name': 'merger_rate_density', 'function_type': 'merger_rate_density_bbh_primordial_ken2022', 'R0': 4.4e-11, 't0': 13.786885302009708}

param_sampler_type:
  gw_parameters_rvs: None
  gw_parameters_rvs_njit: None

zs:
  merger_rate_density_based_source_redshift: None

mass_1_source:
  broken_powerlaw_plus_2peaks: {'param_name': 'mass_1_source', 'sampler_type': 'broken_powerlaw_plus_2peaks', 'lam_0': 0.361, 'lam_1': 0.586, 'mpp_1': 9.764, 'sigpp_1': 0.649, 'mpp_2': 32.763, 'sigpp_2': 3.918, 'mlow_1': 5.059, 'delta_m_1': 4.321, 'break_mass': 35.622, 'alpha_1': 1.728, 'alpha_2': 4.512, 'mmax': 300.0, 'normalization_size': 500}
  powerlaw_plus_peak: {'param_name': 'mass_1_source', 'sampler_type': 'powerlaw_plus_peak', '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, 'normalization_size': 500}
  truncated_normal: {'param_name': 'mass_1_source', 'sampler_type': 'truncated_normal', 'mu': 1.4, 'sigma': 0.68, 'x_min': 1.0, 'x_max': 2.5}
  uniform: {'param_name': 'mass_1_source', 'sampler_type': 'uniform', 'x_min': 1.0, 'x_max': 2.5}
  powerlaw: {'param_name': 'mass_1_source', 'sampler_type': 'powerlaw', 'x_min': 1.0, 'x_max': 2.5, 'alpha': 7.7}
  bimodal: {'param_name': 'mass_1_source', 'sampler_type': 'bimodal', 'w': 0.643, 'muL': 1.352, 'sigmaL': 0.08, 'muR': 1.88, 'sigmaR': 0.3, 'mmin': 1.0, 'mmax': 2.3}
  broken_powerlaw: {'mminbh': 26, 'mmaxbh': 125, 'alpha_1': 6.75, 'alpha_2': 6.75, 'b': 0.5, 'delta_m': 5}

mass_ratio:
  powerlaw_with_smoothing: {'param_name': 'mass_ratio', 'sampler_type': 'powerlaw_with_smoothing', 'q_min': 0.01, 'q_max': 1.0, 'mlow_2': 3.551, 'mmax': 300.0, 'beta': 1.171, 'delta_m': 4.91}

mass_2_source:
  truncated_normal: {'param_name': 'mass_2_source', 'sampler_type': 'truncated_normal', 'mu': 1.4, 'sigma': 0.68, 'x_min': 1.0, 'x_max': 2.5}
  uniform: {'param_name': 'mass_2_source', 'sampler_type': 'uniform', 'x_min': 1.0, 'x_max': 2.5}
  powerlaw: {'param_name': 'mass_2_source', 'sampler_type': 'powerlaw', 'x_min': 1.0, 'x_max': 2.5, 'alpha': 7.7}
  bimodal: {'param_name': 'mass_2_source', 'sampler_type': 'bimodal', 'w': 0.643, 'muL': 1.352, 'sigmaL': 0.08, 'muR': 1.88, 'sigmaR': 0.3, 'mmin': 1.0, 'mmax': 2.3}

a_1:
  constant_values_n_size: {'param_name': 'a_1', 'sampler_type': 'constant_values_n_size', 'value': 0.0}
  uniform: {'param_name': 'a_1', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 0.8}
  truncated_normal: {'param_name': 'a_1', 'sampler_type': 'truncated_normal', 'x_min': 0.0, 'x_max': 1.0, 'mu': 0.085, 'sigma': 0.33}

a_2:
  constant_values_n_size: {'param_name': 'a_2', 'sampler_type': 'constant_values_n_size', 'value': 0.0}
  uniform: {'param_name': 'a_2', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 0.8}
  truncated_normal: {'param_name': 'a_2', 'sampler_type': 'truncated_normal', 'x_min': 0.0, 'x_max': 1.0, 'mu': 0.085, 'sigma': 0.33}

tilt_1:
  constant_values_n_size: {'param_name': 'tilt_1', 'sampler_type': 'constant_values_n_size', 'value': 0.0}
  sampler_sine: {'param_name': 'tilt_1', 'sampler_type': 'sampler_sine', 'tilt_1_min': 0.0, 'tilt_1_max': 3.141592653589793}
  gaussian_plus_isotropic: {'param_name': 'tilt_1', 'sampler_type': 'gaussian_plus_isotropic', 'tilt_1_min': 0.0, 'tilt_1_max': 3.141592653589793, 'mu_t': 0.426, 'sigma_t': 1.222, 'zeta': 0.652}

tilt_2:
  constant_values_n_size: {'param_name': 'tilt_2', 'sampler_type': 'constant_values_n_size', 'value': 0.0}
  sampler_sine: {'param_name': 'tilt_2', 'sampler_type': 'sampler_sine', 'tilt_2_min': 0.0, 'tilt_2_max': 3.141592653589793}
  gaussian_plus_isotropic_joint: {'param_name': 'tilt_2', 'sampler_type': 'gaussian_plus_isotropic_joint', 'tilt_2_min': 0.0, 'tilt_2_max': 3.141592653589793, 'mu_t': 0.426, 'sigma_t': 1.222, 'zeta': 0.652}

phi_12:
  constant_values_n_size: {'param_name': 'phi_12', 'sampler_type': 'constant_values_n_size', 'value': 0.0}
  uniform: {'param_name': 'phi_12', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 6.283185307179586}

phi_jl:
  constant_values_n_size: {'param_name': 'phi_jl', 'sampler_type': 'constant_values_n_size', 'value': 0.0}
  uniform: {'param_name': 'phi_jl', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 6.283185307179586}

geocent_time:
  uniform: {'param_name': 'geocent_time', 'sampler_type': 'uniform', 'x_min': 1238166018, 'x_max': 1269723618.0}
  constant_values_n_size: {'param_name': 'geocent_time', 'sampler_type': 'constant_values_n_size', 'value': 1238166018}

ra:
  uniform: {'param_name': 'ra', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 6.283185307179586}
  constant_values_n_size: {'param_name': 'ra', 'sampler_type': 'constant_values_n_size', 'value': 0.0}

dec:
  sampler_cosine: {'param_name': 'dec', 'sampler_type': 'sampler_cosine', 'x_min': -1.5707963267948966, 'x_max': 1.5707963267948966}
  constant_values_n_size: {'param_name': 'dec', 'sampler_type': 'constant_values_n_size', 'value': 0.0}
  uniform: {'param_name': 'dec', 'sampler_type': 'uniform', 'x_min': -1.5707963267948966, 'x_max': 1.5707963267948966}

phase:
  uniform: {'param_name': 'phase', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 6.283185307179586}
  constant_values_n_size: {'param_name': 'phase', 'sampler_type': 'constant_values_n_size', 'value': 0.0}

psi:
  uniform: {'param_name': 'psi', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 3.141592653589793}
  constant_values_n_size: {'param_name': 'psi', 'sampler_type': 'constant_values_n_size', 'value': 0.0}

theta_jn:
  sampler_sine: {'param_name': 'theta_jn', 'sampler_type': 'sampler_sine', 'x_min': 0.0, 'x_max': 3.141592653589793}
  constant_values_n_size: {'param_name': 'theta_jn', 'sampler_type': 'constant_values_n_size', 'value': 0.0}
  uniform: {'param_name': 'theta_jn', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 3.141592653589793}

2.3. Accessing functions as Class Attributes (Lensed)

[4]:

# Display all available lens parameter sampler functions and their parameters
print("Built-in lens parameter sampler functions and parameters:\n")
for func_name, func_params in ler.available_lens_priors.items():
    print(f"{func_name}:")

    if isinstance(func_params, dict):
        for param_name, param_value in func_params.items():
            print(f"  {param_name}: {param_value}")
    else:
        print(f"  {func_params}")
    print()

Built-in lens parameter sampler functions and parameters:

zs_sl:
  strongly_lensed_source_redshift: {'tau_approximation': True}

lens_redshift_sl:
  lens_redshift_strongly_lensed_sis_analytical: None
  lens_redshift_strongly_lensed_numerical: {'param_name': 'lens_redshift_sl', 'sampler_type': 'lens_redshift_strongly_lensed_numerical', 'lens_type': 'epl_shear_galaxy', 'integration_size': 50000, 'use_multiprocessing': False, 'cross_section_epl_shear_interpolation': False}

lens_redshift:
  lens_redshift_intrinsic_numerical: None

velocity_dispersion:
  gengamma: {'param_name': 'velocity_dispersion', 'sampler_type': '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: {'param_name': 'velocity_dispersion', 'sampler_type': '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: {'param_name': 'velocity_dispersion', 'sampler_type': '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: {'param_name': 'velocity_dispersion', 'sampler_type': '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:
  rayleigh: {'param_name': 'axis_ratio', 'sampler_type': 'rayleigh', 'q_min': 0.2, 'q_max': 1.0}
  axis_ratio_padilla_strauss: {'param_name': 'axis_ratio', 'sampler_type': 'axis_ratio_padilla_strauss', 'q_min': 0.2, 'q_max': 1.0}
  uniform: {'param_name': 'axis_ratio', 'sampler_type': 'uniform', 'x_min': 0.2, 'x_max': 1.0}
  constant_values_n_size: {'param_name': 'axis_ratio', 'sampler_type': 'constant_values_n_size', 'value': 1.0}

axis_rotation_angle:
  uniform: {'param_name': 'axis_rotation_angle', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 6.283185307179586}
  constant_values_n_size: {'param_name': 'axis_rotation_angle', 'sampler_type': 'constant_values_n_size', 'value': 0.0}

density_profile_slope:
  normal: {'param_name': 'density_profile_slope', 'sampler_type': 'normal', 'mu': 1.99, 'sigma': 0.149}
  constant_values_n_size: {'param_name': 'density_profile_slope', 'sampler_type': 'constant_values_n_size', 'value': 2.0}

external_shear1:
  normal: {'param_name': 'external_shear1', 'sampler_type': 'normal', 'mu': 0.0, 'sigma': 0.05}
  constant_values_n_size: {'param_name': 'external_shear1', 'sampler_type': 'constant_values_n_size', 'value': 0.0}

external_shear2:
  normal: {'param_name': 'external_shear2', 'sampler_type': 'normal', 'mu': 0.0, 'sigma': 0.05}
  constant_values_n_size: {'param_name': 'external_shear2', 'sampler_type': 'constant_values_n_size', 'value': 0.0}


[5]:

# use the following code to inspect one of the velocity dispersion sampler functions
# print(ler.gengamma.__doc__)

# Test one of the velocity dispersion sampler functions
print("\nTesting gengamma sampler function")
size = 5
print(f"Velocity dispersion: {ler.gengamma(size)} km/s")


Testing gengamma sampler function
gengamma interpolator will be generated at ./interpolator_json/velocity_dispersion/gengamma_1.json
Velocity dispersion: [130.32498575 152.65861243 153.17794808 100.9058325
 182.76520196] km/s

2.4. Accessing functions from the ler.gw_source_population module (Unlensed)

[6]:

import ler.gw_source_population as gsp

for prior in gsp.available_prior_list():
    print(prior)

merger_rate_density_bbh_oguri2018_function
merger_rate_density_bbh_popIII_ken2022_function
merger_rate_density_madau_dickinson2014_function
merger_rate_density_madau_dickinson_belczynski_ng_function
merger_rate_density_bbh_primordial_ken2022_function
sfr_madau_fragos2017_with_bbh_td
sfr_madau_dickinson2014_with_bbh_td
sfr_madau_fragos2017_with_bns_td
sfr_madau_dickinson2014_with_bns_td
sfr_madau_fragos2017
sfr_madau_dickinson2014
ng2022_lognormal_joint_pdf
binary_masses_BBH_popIII_lognormal_rvs
binary_masses_BBH_primordial_lognormal_rvs
bimodal_pdf
binary_masses_BNS_bimodal_rvs
broken_powerlaw_pdf
gaussian_plus_isotropic_pdf
gaussian_plus_isotropic_joint_pdf
powerlaw_pdf
powerlaw_rvs
truncated_normal_pdf
truncated_normal_rvs
powerlaw_with_smoothing
powerlaw_plus_peak_pdf
powerlaw_plus_peak_function
powerlaw_plus_peak_rvs
broken_powerlaw_plus_2peaks_pdf
broken_powerlaw_plus_2peaks_function
broken_powerlaw_plus_2peaks_rvs
mass_ratio_powerlaw_with_smoothing_pdf
mass_ratio_powerlaw_with_smoothing_rvs

[7]:

# use the following code to inspect one of the merger rate density function
# print(gsp.merger_rate_density_bbh_oguri2018.__doc__)

# Test one of the merger rate density function
print("\nTesting merger_rate_density_bbh_oguri2018 function")
zs =  np.array([0.1, 0.2, 0.3])
print(f"Redshifts: {zs}")
print(f"Merger Rate Denisty: {gsp.merger_rate_density_bbh_oguri2018_function(zs)} Mpc^-3 yr^-1")


Testing merger_rate_density_bbh_oguri2018 function
Redshifts: [0.1 0.2 0.3]
Merger Rate Denisty: [2.21298914e-08 2.57321630e-08 2.98600744e-08] Mpc^-3 yr^-1

2.5. Accessing functions from the ler.lens_galaxy_population module (Lensed)

[8]:

import ler.lens_galaxy_population as lgp

for sampler in lgp.available_sampler_list():
    print(sampler)

lens_redshift_strongly_lensed_sis_analytical_pdf
lens_redshift_strongly_lensed_sis_analytical_rvs
velocity_dispersion_ewoud_denisty_function
velocity_dispersion_bernardi_denisty_function
gengamma_function
gengamma_pdf
gengamma_rvs
rayleigh_rvs
rayleigh_pdf
axis_ratio_padilla_strauss_rvs
axis_ratio_padilla_strauss_pdf
bounded_normal_sample
rejection_sampler
importance_sampler
importance_sampler_mp

[9]:

# use the following code to inspect one of the velocity dispersion sampler functions
# print(lgp.gengamma_rvs.__doc__)

# Test one of the velocity dispersion sampler functions
print("\nTesting gengamma sampler function")
size = 5
print(f"Velocity dispersion: {lgp.gengamma_rvs(size)} km/s")


Testing gengamma sampler function
Velocity dispersion: [151.09757508 182.54117392 201.69941703 180.6957201
 110.82457502] km/s

3. Using Custom Functions in LeR Initialization

The ler package allows full customization of sampling functions and detection criteria. This section demonstrates a Binary Neutron Star (BNS) configuration with custom settings:

Component

Custom Configuration

Default (BBH)

Event Type

BNS (non-spinning)

BBH (spinning, aligned)

Merger Rate

GWTC-3 based

GWTC-4 based

Source Masses

Uniform 1.0-2.5 \(M_{\odot}\)

truncated normal

Lens Model

SIE (Singular Isothermal)

EPL+Shear

Velocity Dispersion

\(\sigma_* = 161\) km/s

\(\sigma_* = 113.78\) km/s

Axis Ratio

Padilla & Strauss (2008)

Rayleigh distribution

Detectors

3G (ET, CE), \(\rho_{\rm opt}\) > 12

O4 (H1, L1, V1), \(\rho_{\rm obs}\) > 10

Notes:

  • GW priors/samplers must be a function with size as the only input argument, or a ler.utils.FunctionConditioning class object (use cases shown in later sections). Use numba.njit decorator for prior/sampler functions when possible. ler.utils.FunctionConditioning creates interpolators for the custom functions to speed up the calculations relying on numba njit.

  • Lens parameter samplers must be ler.utils.FunctionConditioning class object. Generally it takes it the function (1D or 2D) and create njit interpolators pdf and inverse rvs (sampler).

  • Merger rate density must be a function of redshift, i.e., \(F(z_s)\).

  • Velocity dispersion function (galaxy number density) must use a function of velocity dispersion, i.e., \(F(\sigma)\) or \(F(\sigma, z_l)\), to create the ler.utils.FunctionConditioning class object.

3.1. Custom Event Type with non-spinning configuration

Using event_type='BNS' in the LeR class initialization will default to the following GW parameter priors corresponding to BNS. Other allowed event types are ‘BBH’ and ‘NSBH’.

gw_functions = dict(
    merger_rate_density = 'merger_rate_density_madau_dickinson2014',
),
gw_functions_params = dict(
    merger_rate_density = dict(
        param_name='merger_rate_density',
        function_type='merger_rate_density_madau_dickinson2014',
        'R0': 89 * 1e-9, 'a': 0.015, 'b': 2.7, 'c': 2.9, 'd': 5.6
    ),
),
gw_priors = dict(
    mass_1_source = 'truncated_normal',
    mass_2_source = 'truncated_normal',
    mass_ratio = None,
    a_1 = 'uniform',
    a_2 = 'uniform',
),
gw_priors_params= dict(
    mass_1_source = dict(
        param_name='mass_1_source', sampler_type='truncated_normal',
        mu=1.4, sigma=0.68, x_min=1.0, x_max=2.5
    ),
    mass_2_source = dict(
        param_name='mass_2_source', sampler_type='truncated_normal',
        mu=1.4, sigma=0.68, x_min=1.0, x_max=2.5
    ),
    mass_ratio = None,
    a_1 = dict(
        param_name='a_1', sampler_type='uniform',
        xmin=-0.4, xmax=0.4
    ),
    a_2 = dict(
        param_name='a_2', sampler_type='uniform',
        xmin=-0.4, xmax=0.4
    ),
),

We will change some of these priors with our custom ones in the next sections.

For non-spinning configuration (for faster calculation in our example), we can set:

spin_zero=True,
spin_precession=False,

3.2. Custom Merger Rate Density

Using the default BNS merger rate density prior model with the local merger rate density change from the default value of \(R_0 = 89 \times 10^{-9} \, \text{Mpc}^{-3}\text{yr}^{-1}\) (GWTC-4) to \(R_0 = 105.5 \times 10^{-9} \, \text{Mpc}^{-3}\text{yr}^{-1}\) (GWTC-3).

[10]:

merger_rate_density_function = 'merger_rate_density_madau_dickinson2014'
merger_rate_density_input_args = dict(
    param_name='merger_rate_density',
    function_type='merger_rate_density_madau_dickinson2014',
    R0=105.5e-9,
    a=0.015,
    b=2.7,
    c=2.9,
    d=5.6,
)

print("Merger rate density function:", merger_rate_density_function)
print("Parameters:", merger_rate_density_input_args)

Merger rate density function: merger_rate_density_madau_dickinson2014
Parameters: {'param_name': 'merger_rate_density', 'function_type': 'merger_rate_density_madau_dickinson2014', 'R0': 1.055e-07, 'a': 0.015, 'b': 2.7, 'c': 2.9, 'd': 5.6}

3.3. Custom Source Frame Masses (\(m_1^{src}\), \(m_2^{src}\))

Using a uniform distribution to sample the binary masses mass_1_source and mass_2_source. Swapping of values if mass_1_source < mass_2_source will be done internally.

Note: Alternatively, in LeR initialization, you can provide mass_ratio instead of mass_2_source. In that case, mass_2_source will be internally calculated as mass_1_source * mass_ratio.

[11]:

import numpy as np
from numba import njit

mmin=1.0
mmax=2.5

# define your custom function
# it should have 'size' as the only argument
# the same function will be used to sample both mass_1_source and mass_2_source. The code will swap the values if mass_1_source < mass_2_source.
@njit
def source_frame_masses_uniform(size):
    """
    Function to sample mass1 or mass2 from a uniform distribution between mmin and mmax.

    Parameters
    ----------
    size : `int`
        Number of samples to draw

    Returns
    -------
    mass_source : `numpy.ndarray`
        Array of mass samples in the source frame (Msun)
    """

    mass_source = np.random.uniform(mmin, mmax, size)

    return mass_source

# test
mass_source = source_frame_masses_uniform(size=5)
print(f"mass: {mass_source} M_sun")

mass: [2.47061195 1.22504778 2.47096762 1.35593826 1.96734825] M_sun

3.4. Custom Lens Model

Using lens_model='sie_galaxy' in the LeR initialization will use the following default settings. Other allowed lens models are ‘epl_shear_galaxy’ and ‘sis_galaxy’.

lens_priors = dict(
    zs_sl="strongly_lensed_source_redshift",
    lens_redshift_sl="lens_redshift_strongly_lensed_numerical",
    lens_redshift="lens_redshift_intrinsic_numerical",
    velocity_dispersion="velocity_dispersion_ewoud",
    axis_ratio="rayleigh",
    axis_rotation_angle="uniform",
    external_shear1="constant_values_n_size",
    external_shear2="constant_values_n_size",
    density_profile_slope="constant_values_n_size",
    external_shear1_sl="constant_values_n_size",
    external_shear2_sl="constant_values_n_size",
    density_profile_slope_sl="constant_values_n_size",
)
lens_priors_params = dict(
    zs_sl=None,
    lens_redshift_sl=dict(
        param_name = "lens_redshift_sl",
        sampler_type = "lens_redshift_strongly_lensed_numerical",
        lens_type = 'sie_galaxy',
        integration_size=25000, use_multiprocessing=False
    ),
    lens_redshift=None,
    velocity_dispersion=dict(
        param_name='velocity_dispersion', sampler_type='velocity_dispersion_ewoud',
        sigma_min=100.0,
        sigma_max=400.0,
        alpha=0.94,
        beta=1.85,
        phistar=2.099e-2 * (self.cosmo.h/0.7)**3,
        sigmastar=113.78,
    ),
    axis_ratio=dict(
        param_name='axis_ratio', sampler_type='rayleigh',
        q_min=0.2, q_max=1.0
    ),
    axis_rotation_angle=dict(
        param_name='axis_rotation_angle', sampler_type='uniform',
        x_min=0.0, x_max=2 * np.pi
    ),
    external_shear1=dict(param_name='external_shear1', sampler_type= 'constant_values_n_size', value=0.0),
    external_shear2=dict(param_name='external_shear2', sampler_type= 'constant_values_n_size', value=0.0),
    density_profile_slope=dict(param_name='density_profile_slope', sampler_type='constant_values_n_size', value=2.0),
    external_shear1_sl=dict(param_name='external_shear1_sl', sampler_type= 'constant_values_n_size', value=0.0),
    external_shear2_sl=dict(param_name='external_shear2_sl', sampler_type= 'constant_values_n_size', value=0.0),
    density_profile_slope_sl=dict(param_name='density_profile_slope', sampler_type='constant_values_n_size', value=2.0),
)
lens_functions = dict(
    param_sampler_type="epl_shear_sl_parameters_rvs",
    cross_section_based_sampler="importance_sampling_with_cross_section",
    optical_depth="optical_depth_numerical",
    cross_section="cross_section_sie_feixu",
)
lens_functions_params = dict(
    param_sampler_type=None,
    cross_section_based_sampler=dict(n_prop_factor=10),
    optical_depth=dict(param_name="optical_depth", function_type="optical_depth_numerical"),
    cross_section=None,
)

We will change the velocity dispersion and axis ratio samplers to our custom functions.

3.4.1. Custom Velocity Dispersion

Unlike other lens parameter, custom velocity dispersion needs to be a density function (galaxy number density as a function of velocity dispersion).

[12]:

import numpy as np
from ler.utils import FunctionConditioning, generate_mixed_grid
# import number density function wrt velocity dispersion and redshift
from ler.lens_galaxy_population import velocity_dispersion_ewoud_denisty_function

velocity_dispersion_args = dict(
    sigma_min=100,  # default ler sigma_min=100 km/s
    sigma_max=400,  # default ler sigma_max=400 km/s
    alpha=0.94,
    beta=1.85,
    phistar=2.099e-2,
    sigmastar=161.0,
)

# identifier_dict dict allows for easy tracking of the generated interpolator in future usage
identifier_dict = {"name": "velocity_dispersion_custom"}
identifier_dict["sigma_min"] = velocity_dispersion_args["sigma_min"]
identifier_dict["sigma_max"] = velocity_dispersion_args["sigma_max"]
identifier_dict["resolution"] = 200
identifier_dict["zl_resolution"] = 48

# setting up inputs for the interpolator
# Note: sigma_array and zl_array will form a 2D grid where the function is evaluated. Gird points are use for cubic spline interpolation.
sigma_array = np.linspace(
    identifier_dict["sigma_min"],
    identifier_dict["sigma_max"],
    identifier_dict["resolution"],
)
z_min = 0.001
z_max = 10.0
z_resolution = identifier_dict["zl_resolution"]
zl_array = generate_mixed_grid(z_min, z_max, z_resolution)

# define the function
number_density_function = lambda sigma, zl: velocity_dispersion_ewoud_denisty_function(  # noqa: E731
    sigma,
    zl,
    alpha=velocity_dispersion_args["alpha"],
    beta=velocity_dispersion_args["beta"],
    phistar=velocity_dispersion_args["phistar"],
    sigmastar=velocity_dispersion_args["sigmastar"],
)

sigma_object = FunctionConditioning(
    function=number_density_function,
    x_array=sigma_array,
    non_negative_function=True,
    conditioned_y_array=zl_array,
    identifier_dict=identifier_dict,
    directory="./interpolator_json",
    sub_directory="velocity_dispersion",
    name=identifier_dict["name"],
    create_new=False,
    create_function_inverse=False,
    create_function=True,
    create_pdf=True,
    create_rvs=True,
    callback="rvs",
    cdf_size=500,
)

# Test
sigma = np.array([100., 200., 300.])
zl = np.array([0.1, 0.2, 0.3])
print(f"Velocity dispersions: {sigma} km/s at lens redshifts: {zl}")
print(f"Velocity dispersion density function values: {sigma_object.function(sigma, zl)} Mpc^-3")
print(f"Random velocity dispersion samples: {sigma_object.rvs(len(zl), zl)} Mpc^-3")

velocity_dispersion_custom interpolator will be generated at ./interpolator_json/velocity_dispersion/velocity_dispersion_custom_2.json
Velocity dispersions: [100. 200. 300.] km/s at lens redshifts: [0.1 0.2 0.3]
Velocity dispersion density function values: [9.75545717e-05 3.09054572e-05 5.78713823e-06] Mpc^-3
Random velocity dispersion samples: [187.89865267 239.40749587 130.15268119] Mpc^-3

3.4.2. Custom Axis Ratio

FunctionConditioning class object creates pdf and rvs, which is required by ler for sampling.

[13]:

from scipy.interpolate import CubicSpline  # noqa: E402
from ler.utils import FunctionConditioning  # noqa: E402

axis_ratio_args = dict(
    q_min=0.2,
    q_max=1.0,
)

identifier_dict = {"name": "axis_ratio_padilla_strauss_custom"}
identifier_dict["q_min"] = axis_ratio_args["q_min"]
identifier_dict["q_max"] = axis_ratio_args["q_max"]
identifier_dict["resolution"] = 500

# Using Padilla and Strauss 2008 distribution for axis ratio
q_array = np.array(
    [0.04903276402927845,0.09210526315789469,0.13596491228070173,0.20789473684210524,0.2899703729522482,0.3230132450331126,0.35350877192982455,0.37946148483792264,0.4219298245614036,0.4689525967235971,0.5075026141512723,0.5226472638550018,0.5640350877192983,0.6096491228070177,0.6500000000000001,0.6864848379226213,0.7377192982456142,0.7787295224817011,0.8007581038689441,0.822786685256187,0.8668438480306729,0.8973684210526317,0.9254385964912283,
    ]
)
pdf = np.array(
    [0.04185262687135349,0.06114520695141845,0.096997499638376,0.1932510900336828,0.39547914337673706,0.49569751276216234,0.6154609137685201,0.7182049959882812,0.920153741243567,1.1573982157399754,1.3353263628106684,1.413149656448315,1.5790713532948977,1.7280185150744938,1.8132994441344819,1.8365803753840484,1.8178662203211204,1.748929843583365,1.688182592496342,1.6274353414093188,1.4948487090314488,1.402785526832393,1.321844068356993,
    ]
)

# Interpolate the pdf
spline = CubicSpline(q_array, pdf, extrapolate=True)
q_array = np.linspace(identifier_dict["q_min"], identifier_dict["q_max"], identifier_dict["resolution"])
pdf = spline(q_array)

q_object = FunctionConditioning(
    function=pdf,  # it also allows precomputed values, besides function
    x_array=q_array,
    non_negative_function=True,
    conditioned_y_array=None,
    identifier_dict=identifier_dict,
    directory="./interpolator_json",
    sub_directory="axis_ratio",
    name="axis_ratio_padilla_strauss",
    create_new=False,
    create_function_inverse=False,
    create_function=True,
    create_pdf=True,
    create_rvs=True,
    callback="rvs",
    cdf_size=500,
)

# test
# sampling
q_samples = q_object.rvs(5)
print(f"Axis ratio samples: {q_samples}")
# pdf
q_pdf = q_object.pdf(q_samples)
print(f"Axis ratio pdf values: {q_pdf}")

axis_ratio_padilla_strauss interpolator will be generated at ./interpolator_json/axis_ratio/axis_ratio_padilla_strauss_1.json
Axis ratio samples: [0.61392959 0.38463519 0.90852885 0.67962473 0.82826483]
Axis ratio pdf values: [1.76358827 0.75046076 1.38887165 1.8601851  1.63387105]

3.5. Custom Detection Criteria

Custom pdet_finder should be of this format:

def custom_pdet_finder(gw_param_dict):
    """
    Parameters
    ----------
    gw_param_dict: dict
        Dictionary containing the GW parameters of the event.

    Returns
    -------
    pdet_dict: dict
        Dictionary with key 'pdet_net'.

  • For our example, we will use 3G detectors (Einstein Telescope and Cosmic Explorer) with SNR threshold of 12.

Note: Even though ler uses pdet_net to determine detectability, we can use standard SNR (or other detector sensitivity arguments) instead, and set a corresponding threshold later.

[21]:

# Define a function that sets detection criteria
from gwsnr import GWSNR

# 3G detectors: Einstein Telescope (ET) and Cosmic Explorer (CE)
mmin, mmax = 1.0, 2.5
zmin, zmax = 0.0, 10.0

# for spin precessing waveform
gwsnr_3g = GWSNR(
    npool=6,
    ifos=['CE'],  # 3G detector network. Can put ifos=['ET','CE'] for multiple detectors
    snr_method='interpolation_aligned_spins',  # BNS have no spins
    mtot_min=2*mmin*(1+zmin),
    mtot_max=2*mmax*(1+zmax),
    sampling_frequency=1024.0,
    waveform_approximant='IMRPhenomXPHM',
    minimum_frequency=20.0,
    gwsnr_verbose=False,
    snr_recalculation = True,
    snr_recalculation_range = [6, 14],
    snr_recalculation_waveform_approximant = 'IMRPhenomXPHM',
)

detection_threshold=10

def detection_criteria(gw_param_dict):
    """
    Custom detection criteria for 3G detectors with optimal SNR as output.
    """

    dict_ = {}
    dict_['pdet_net'] = gwsnr_3g.optimal_snr(gw_param_dict=gw_param_dict)['optimal_snr_net']
    return dict_

# test
gw_param_dict = dict(
    mass_1 = np.array([20.0, 20.0]),
    mass_2 = np.array([10.0, 10.0]),
    luminosity_distance = np.array([1000.0, 40000.0]),
    a_1 = np.array([0.5, 0.5]),
    a_2 = np.array([0.5, 0.5]),
    tilt_1 = np.array([0.5, 0.5]),
    tilt_2 = np.array([0.5, 0.5]),
    phi_12 = np.array([0.5, 0.5]),
    phi_jl = np.array([0.5, 0.5]),
)
detection_dict = detection_criteria(gw_param_dict)
print(f"SNR results: {detection_dict}")
pdet_net = detection_dict['pdet_net']>= detection_threshold
print(f"Detection results: {pdet_net}")


Initializing GWSNR class...

Intel processor has trouble allocating memory when the data is huge. So, by default for IMRPhenomXPHM, duration_max = 64.0. Otherwise, set to some max value like duration_max = 600.0 (10 mins)
Interpolator will be loaded for CE detector from ./interpolator_json/CE/partialSNR_dict_1.json


Recalculating SNR for 1 out of 2 samples in the SNR range of 6 to 14

100%|█████████████████████████████████████████████████████████████████| 1/1 [00:01<00:00,  1.32s/it]

SNR results: {'pdet_net': array([412.10519461,  10.26345904])}
Detection results: [ True  True]

3.6. LeR initialization with custom functions and parameters

Create a LeR instance with custom source priors, lens samplers, and detection criteria for BNS events.

[22]:

from ler import LeR

ler = LeR(
    # Core setup
    npool=6,
    event_type='BNS',
    lens_type='sie_galaxy',

    # Source priors
    gw_functions=dict(
        merger_rate_density=merger_rate_density_function,
    ),
    gw_functions_params=dict(
        merger_rate_density=merger_rate_density_input_args,
    ),
    gw_priors=dict(
        mass_1_source=source_frame_masses_uniform,
        mass_2_source=source_frame_masses_uniform,
    ),
    gw_priors_params=dict(
        mass_1_source=dict(
            param_name='mass_1_source',
            sampler_type='source_frame_masses_uniform',
            mmin=1.0,
            mmax=2.5,
        ),
        mass_2_source=dict(
            param_name='mass_2_source',
            sampler_type='source_frame_masses_uniform',
            mmin=1.0,
            mmax=2.5,
        ),
    ),

    # Lens samplers
    lens_priors=dict(
        velocity_dispersion=sigma_object,
        axis_ratio=q_object,
    ),
    lens_priors_params=dict(
        velocity_dispersion=velocity_dispersion_args,
        axis_ratio=axis_ratio_args,
    ),

    # Custom detection
    pdet_finder=detection_criteria,
    ler_directory='./ler_data_custom',
)


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 user provided custom velocity_dispersion class/object of type ler.utils.FunctionConditioning
using user provided custom axis_ratio class/object of type ler.utils.FunctionConditioning
using ler available axis_rotation_angle function : uniform
using ler available density_profile_slope function : constant_values_n_size
using ler available external_shear1 function : constant_values_n_size
using ler available external_shear2 function : constant_values_n_size
using ler available cross_section function : cross_section_sie_feixu
using ler available lens_redshift_sl function : lens_redshift_strongly_lensed_numerical
Numerically solving the lens redshift distribution...
Using multiprocessing
Computing lens redshift distribution with multiprocessing...

100%|███████████████████████████████████████████████████████████████| 16/16 [00:15<00:00,  1.05it/s]

lens_redshift_strongly_lensed_numerical interpolator will be generated at ./interpolator_json/lens_redshift_sl/lens_redshift_strongly_lensed_numerical_1.json
using ler available lens_redshift function : lens_redshift_intrinsic_numerical
lens_redshift_intrinsic_numerical interpolator will be generated at ./interpolator_json/lens_redshift/lens_redshift_intrinsic_numerical_3.json
using ler available optical depth function : optical_depth_numerical
optical_depth_numerical interpolator will be generated at ./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_dickinson2014
merger_rate_density_madau_dickinson2014 interpolator will be generated at ./interpolator_json/merger_rate_density/merger_rate_density_madau_dickinson2014_3.json
merger_rate_density_madau_dickinson2014_detector_frame interpolator will be generated at ./interpolator_json/merger_rate_density/merger_rate_density_madau_dickinson2014_detector_frame_4.json
merger_rate_density_based_source_redshift interpolator will be generated at ./interpolator_json/merger_rate_density_based_source_redshift/merger_rate_density_based_source_redshift_1.json

Initializing CBCSourceParameterDistribution class...

using ler available zs function : merger_rate_density_based_source_redshift
using user defined custom mass_1_source function
using user defined custom mass_2_source function
No mass_ratio prior provided. mass_ratio = mass_2_source / mass_1_source
using ler available geocent_time function : uniform
using ler available ra function : uniform
using ler available dec function : sampler_cosine
using ler available phase function : uniform
using ler available psi function : uniform
using ler available theta_jn function : sampler_sine
using ler available a_1 function : uniform
using ler available a_2 function : uniform
No tilt_1 prior provided. tilt_1 prior will be set to None
No tilt_2 prior provided. tilt_2 prior will be set to None
No phi_12 prior provided. phi_12 prior will be set to None
No phi_jl prior provided. phi_jl prior will be set to None
Initializing GW parameter samplers...

initializing epl shear njit solver
using ler available zs_sl function : strongly_lensed_source_redshift
strongly_lensed_source_redshift interpolator will be generated at ./interpolator_json/zs_sl/strongly_lensed_source_redshift_1.json
Slower, non-njit and importance sampling based lens parameter sampler will be used.
#-------------------------------------
# LeR initialization input arguments:
#-------------------------------------

    # LeR set up input arguments:
    npool = 6,
    z_min = 0.0,
    z_max = 10.0,
    event_type = 'BNS',
    lens_type = 'sie_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 = <function detection_criteria at 0x129311120>,
    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_custom',
    create_new_interpolator = dict(
        velocity_dispersion = {'create_new': False, 'resolution': 200, 'zl_resolution': 48, 'cdf_size': 400},
        axis_ratio = {'create_new': False, 'resolution': 200, 'sigma_resolution': 48, 'cdf_size': 400},
        lens_redshift_sl = {'create_new': False, 'resolution': 16, 'zl_resolution': 48, 'cdf_size': 400},
        lens_redshift = {'create_new': False, 'resolution': 16, 'zl_resolution': 48, 'cdf_size': 400},
        optical_depth = {'create_new': False, 'resolution': 16, 'cdf_size': 500},
        comoving_distance = {'create_new': False, 'resolution': 500},
        angular_diameter_distance = {'create_new': False, 'resolution': 500},
        angular_diameter_distance_z1z2 = {'create_new': False, 'resolution': 500},
        differential_comoving_volume = {'create_new': False, 'resolution': 500},
        zs_sl = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        cross_section = {'create_new': False, 'resolution': [25, 25, 5, 5, 5], 'spacing_config': {'e1': {'mode': 'two_sided_mixed_grid', 'power_law_part': 'lower', 'spacing_trend': 'increasing', 'power': 2.5, 'value_transition_fraction': 0.7, 'num_transition_fraction': 0.9, 'auto_match_slope': True}, 'e2': {'mode': 'two_sided_mixed_grid', 'power_law_part': 'lower', 'spacing_trend': 'increasing', 'power': 2.5, 'value_transition_fraction': 0.7, 'num_transition_fraction': 0.9, 'auto_match_slope': True}, 'gamma': {'mode': 'two_sided_mixed_grid', 'power_law_part': 'lower', 'spacing_trend': 'increasing', 'power': 2.5, 'value_transition_fraction': 0.7, 'num_transition_fraction': 0.9, 'auto_match_slope': True}, 'gamma1': {'mode': 'two_sided_mixed_grid', 'power_law_part': 'lower', 'spacing_trend': 'increasing', 'power': 2.5, 'value_transition_fraction': 0.7, 'num_transition_fraction': 0.9, 'auto_match_slope': True}, 'gamma2': {'mode': 'two_sided_mixed_grid', 'power_law_part': 'lower', 'spacing_trend': 'increasing', 'power': 2.5, 'value_transition_fraction': 0.7, 'num_transition_fraction': 0.9, 'auto_match_slope': True}}},
        merger_rate_density = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        redshift_distribution = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        luminosity_distance = {'create_new': False, 'resolution': 500},
        mass_1_source = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        mass_ratio = {'create_new': False, 'resolution': 100, 'm1_resolution': 200, 'cdf_size': 200},
        mass_2_source = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        geocent_time = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        ra = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        dec = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        phase = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        psi = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        theta_jn = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        a_1 = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        a_2 = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        tilt_1 = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        tilt_2 = {'create_new': False, 'resolution': 200, 'tilt_1_resolution': 100, 'cdf_size': 200},
        phi_12 = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
        phi_jl = {'create_new': False, 'resolution': 200, 'cdf_size': 500},
    ),

    # LeR also takes other CBCSourceParameterDistribution class input arguments as kwargs, as follows:
    gw_functions = dict(
        merger_rate_density = 'merger_rate_density_madau_dickinson2014',
        param_sampler_type = 'gw_parameters_rvs',
    ),
    gw_functions_params = dict(
        merger_rate_density = {'param_name': 'merger_rate_density', 'function_type': 'merger_rate_density_madau_dickinson2014', 'R0': 1.055e-07, 'a': 0.015, 'b': 2.7, 'c': 2.9, 'd': 5.6},
        param_sampler_type = None,
    ),
    gw_priors = dict(
        zs = 'merger_rate_density_based_source_redshift',
        mass_1_source = CPUDispatcher(<function source_frame_masses_uniform at 0x1290bd800>),
        mass_ratio = None,
        mass_2_source = CPUDispatcher(<function source_frame_masses_uniform at 0x1290bd800>),
        geocent_time = 'uniform',
        ra = 'uniform',
        dec = 'sampler_cosine',
        phase = 'uniform',
        psi = 'uniform',
        theta_jn = 'sampler_sine',
        a_1 = 'uniform',
        a_2 = 'uniform',
        tilt_1 = None,
        tilt_2 = None,
        phi_12 = None,
        phi_jl = None,
    ),
    gw_priors_params = dict(
        zs = None,
        mass_1_source = {'param_name': 'mass_1_source', 'sampler_type': 'source_frame_masses_uniform', 'mmin': 1.0, 'mmax': 2.5},
        mass_ratio = None,
        mass_2_source = {'param_name': 'mass_2_source', 'sampler_type': 'source_frame_masses_uniform', 'mmin': 1.0, 'mmax': 2.5},
        geocent_time = {'param_name': 'geocent_time', 'sampler_type': 'uniform', 'x_min': 1238166018, 'x_max': 1269702018},
        ra = {'param_name': 'ra', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 6.283185307179586},
        dec = {'param_name': 'dec', 'sampler_type': 'sampler_cosine', 'x_min': -1.5707963267948966, 'x_max': 1.5707963267948966},
        phase = {'param_name': 'phase', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 6.283185307179586},
        psi = {'param_name': 'psi', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 3.141592653589793},
        theta_jn = {'param_name': 'theta_jn', 'sampler_type': 'sampler_sine', 'x_min': 0.0, 'x_max': 3.141592653589793},
        a_1 = {'param_name': 'a_1', 'sampler_type': 'uniform', 'x_min': -0.4, 'x_max': 0.4},
        a_2 = {'param_name': 'a_2', 'sampler_type': 'uniform', 'x_min': -0.4, 'x_max': 0.4},
        tilt_1 = None,
        tilt_2 = None,
        phi_12 = None,
        phi_jl = None,
    ),
    spin_zero = False,
    spin_precession = False,

    # LeR also takes other LensGalaxyParameterDistribution class input arguments as kwargs, as follows:
    lens_functions = dict(
        param_sampler_type = 'epl_shear_sl_parameters_rvs',
        cross_section_based_sampler = 'importance_sampler_partial',
        optical_depth = 'optical_depth_numerical',
        cross_section = 'cross_section_sie_feixu',
    ),
    lens_functions_params = dict(
        param_sampler_type = None,
        cross_section_based_sampler = {'n_prop': 50, 'threshold_factor': 0.0, 'sigma_min': 100.0, 'sigma_max': 400.0, 'q_min': 0.2, 'q_max': 1.0, 'phi_min': 0.0, 'phi_max': 6.283185307179586, 'gamma_min': 2.0, 'gamma_max': 2.0, 'shear_min': 0.0, 'shear_max': 0.0},
        optical_depth = {'param_name': 'optical_depth', 'function_type': 'optical_depth_numerical'},
        cross_section = None,
    ),
    lens_priors = dict(
        zs_sl = 'strongly_lensed_source_redshift',
        lens_redshift_sl = 'lens_redshift_strongly_lensed_numerical',
        lens_redshift = 'lens_redshift_intrinsic_numerical',
        velocity_dispersion = <ler.utils.function_interpolation.FunctionConditioning object at 0x1291c2310>,
        axis_ratio = <ler.utils.function_interpolation.FunctionConditioning object at 0x12d1502d0>,
        axis_rotation_angle = 'uniform',
        external_shear1 = 'constant_values_n_size',
        external_shear2 = 'constant_values_n_size',
        density_profile_slope = 'constant_values_n_size',
    ),
    lens_priors_params = dict(
        zs_sl = {'tau_approximation': True},
        lens_redshift_sl = {'param_name': 'lens_redshift_sl', 'sampler_type': 'lens_redshift_strongly_lensed_numerical', 'lens_type': 'sie_galaxy', 'integration_size': 50000, 'use_multiprocessing': False, 'cross_section_epl_shear_interpolation': False},
        lens_redshift = None,
        velocity_dispersion = {'param_name': 'velocity_dispersion', 'sampler_type': 'velocity_dispersion_ewoud', 'sigma_min': 100, 'sigma_max': 400, 'alpha': 0.94, 'beta': 1.85, 'phistar': 0.02099, 'sigmastar': 161.0},
        axis_ratio = {'param_name': 'axis_ratio', 'sampler_type': 'rayleigh', 'q_min': 0.2, 'q_max': 1.0},
        axis_rotation_angle = {'param_name': 'axis_rotation_angle', 'sampler_type': 'uniform', 'x_min': 0.0, 'x_max': 6.283185307179586},
        external_shear1 = {'param_name': 'external_shear1', 'sampler_type': 'constant_values_n_size', 'value': 0.0},
        external_shear2 = {'param_name': 'external_shear2', 'sampler_type': 'constant_values_n_size', 'value': 0.0},
        density_profile_slope = {'param_name': 'density_profile_slope', 'sampler_type': 'constant_values_n_size', 'value': 2.0},
    ),

    # LeR also takes other ImageProperties class input arguments as kwargs, as follows:
    n_min_images = 2,
    n_max_images = 4,
    time_window = 31536000.0,
    include_effective_parameters = False,
    lens_model_list = ['EPL_NUMBA', 'SHEAR'],
    image_properties_function = image_properties_epl_shear_njit,
    multiprocessing_verbose = True,
    include_redundant_parameters = False,

3.7. Simulate Unlensed Population

Generate unlensed BNS events using the custom mass distribution and merger rate density.

[23]:

unlensed_params = ler.unlensed_cbc_statistics(size=20000, batch_size=10000, resume=False)

print(f"\nTotal unlensed events simulated: {len(unlensed_params['zs'])}")
print(f"Sample source redshift values (first 5): {unlensed_params['zs'][:5]}")
print(f"Sample masses (first 3):")
print(f"  mass_1: {unlensed_params['mass_1_source'][:3]}")
print(f"  mass_2: {unlensed_params['mass_2_source'][:3]}")

unlensed params will be stored in ./ler_data_custom/unlensed_param.json
removing ./ler_data_custom/unlensed_param.json if it exists
Batch no. 1
sampling gw source params...
calculating pdet...
Recalculating SNR for 4361 out of 10000 samples in the SNR range of 6 to 14

100%|██████████████████████████████████████████████████████████| 4361/4361 [00:12<00:00, 349.15it/s]

Batch no. 2
sampling gw source params...
calculating pdet...
Recalculating SNR for 4308 out of 10000 samples in the SNR range of 6 to 14

100%|██████████████████████████████████████████████████████████| 4308/4308 [00:11<00:00, 364.34it/s]

saving all unlensed parameters in ./ler_data_custom/unlensed_param.json

Total unlensed events simulated: 20000
Sample source redshift values (first 5): [2.38335663 1.66229721 2.0978767  1.30300035 1.71605284]
Sample masses (first 3):
  mass_1: [1.13298159 2.28624153 1.71188966]
  mass_2: [1.5889242  2.15626486 1.23868748]

3.8. Simulate Lensed Population

Generate lensed BNS events with custom lens parameters.

[24]:

lensed_params = ler.lensed_cbc_statistics(size=20000, batch_size=10000, resume=False)

print(f"\nTotal lensed events simulated: {len(lensed_params['zs'])}")
print(f"Sample source redshift values (first 5): {lensed_params['zs'][:5]}")
print(f"Lens parameters (first 3):")
lens_params = ['zl', 'sigma', 'q']
for param in lens_params:
    if param in lensed_params:
        print(f"  {param}: {lensed_params[param][:3]}")

lensed params will be stored in ./ler_data_custom/lensed_param.json
removing ./ler_data_custom/lensed_param.json if it exists
Batch no. 1
sampling lensed params...
sampling lens parameters with epl_shear_sl_parameters_rvs...
computing image properties using ler's epl+shear (analytical, njit) solver...
Invalid sample found. Resampling 661 lensed events...
sampling lens parameters with epl_shear_sl_parameters_rvs...
computing image properties using ler's epl+shear (analytical, njit) solver...
Invalid sample found. Resampling 31 lensed events...
sampling lens parameters with epl_shear_sl_parameters_rvs...
computing image properties using ler's epl+shear (analytical, njit) solver...
calculating pdet...
Recalculating SNR for 1 out of 1 samples in the SNR range of 6 to 14

100%|█████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00,  1.07it/s]

Recalculating SNR for 4598 out of 10000 samples in the SNR range of 6 to 14


100%|██████████████████████████████████████████████████████████| 4598/4598 [00:09<00:00, 472.00it/s]

Recalculating SNR for 1836 out of 7765 samples in the SNR range of 6 to 14


100%|██████████████████████████████████████████████████████████| 1836/1836 [00:04<00:00, 372.57it/s]

Recalculating SNR for 201 out of 515 samples in the SNR range of 6 to 14

100%|████████████████████████████████████████████████████████████| 201/201 [00:01<00:00, 144.79it/s]

Recalculating SNR for 156 out of 492 samples in the SNR range of 6 to 14


100%|████████████████████████████████████████████████████████████| 156/156 [00:01<00:00, 114.71it/s]

Batch no. 2
sampling lensed params...
sampling lens parameters with epl_shear_sl_parameters_rvs...
computing image properties using ler's epl+shear (analytical, njit) solver...
Invalid sample found. Resampling 668 lensed events...
sampling lens parameters with epl_shear_sl_parameters_rvs...
computing image properties using ler's epl+shear (analytical, njit) solver...
Invalid sample found. Resampling 37 lensed events...
sampling lens parameters with epl_shear_sl_parameters_rvs...
computing image properties using ler's epl+shear (analytical, njit) solver...
Invalid sample found. Resampling 1 lensed events...
sampling lens parameters with epl_shear_sl_parameters_rvs...
computing image properties using ler's epl+shear (analytical, njit) solver...
Invalid sample found. Resampling 1 lensed events...
sampling lens parameters with epl_shear_sl_parameters_rvs...
computing image properties using ler's epl+shear (analytical, njit) solver...
calculating pdet...
Recalculating SNR for 4581 out of 10000 samples in the SNR range of 6 to 14

100%|██████████████████████████████████████████████████████████| 4581/4581 [00:09<00:00, 475.12it/s]

Recalculating SNR for 1975 out of 7728 samples in the SNR range of 6 to 14


100%|██████████████████████████████████████████████████████████| 1975/1975 [00:05<00:00, 377.59it/s]

Recalculating SNR for 222 out of 515 samples in the SNR range of 6 to 14


100%|████████████████████████████████████████████████████████████| 222/222 [00:01<00:00, 153.29it/s]

Recalculating SNR for 170 out of 490 samples in the SNR range of 6 to 14


100%|████████████████████████████████████████████████████████████| 170/170 [00:01<00:00, 123.38it/s]

saving all lensed parameters in ./ler_data_custom/lensed_param.json

Total lensed events simulated: 20000
Sample source redshift values (first 5): [4.393585   4.43653119 7.74806238 5.82014538 4.66031067]
Lens parameters (first 3):
  zl: [0.54230372 0.33748345 1.94065031]
  sigma: [362.33283443 354.76713343 266.12957698]
  q: [0.6649898  0.7713353  0.44194038]

[25]:

# include other useful parameters in the output dictionary.
# It is omitted by default to save runtime and memory.
# For theta_E, n_images, mass_1, mass_2, luminosity_distance:
lensed_param = ler.recover_redundant_parameters(lensed_params)
# For effective_luminosity_distance, effective_geocent_time, effective_phase, effective_ra, effective_dec:
lensed_param = ler.produce_effective_params(lensed_param)

# or you can simply set 'include_redundant_parameters=True, include_effective_parameters=True' in the input args of LeR initialization to include these parameters in the output dictionary by default.

3.9. Calculate Rates and Compare Results

Calculate unlensed and lensed detection rates using the custom detection criteria.

[26]:

# Calculate detection rates
rate_unlensed, unlensed_param_detectable = ler.unlensed_rate(
    unlensed_param=unlensed_params,
    pdet_threshold=detection_threshold
)
rate_lensed, lensed_param_detectable = ler.lensed_rate(
    lensed_param=lensed_params,
    pdet_threshold=[detection_threshold, detection_threshold],
    num_img=[1, 1]
)
rate_ratio = ler.rate_ratio()

print(f"\n=== Detection Rates (BNS with Custom Configuration) ===")
print(f"Unlensed rate: {rate_unlensed:.4e} events/year")
print(f"Lensed rate:   {rate_lensed:.4e} events/year")
print(f"Rate ratio (lensed/unlensed): {rate_ratio:.4f}")

Using provided unlensed_param dict...
total unlensed rate (yr^-1): 118984.60266071376
number of simulated unlensed detectable events: 6054
number of simulated all unlensed events: 20000
storing detectable params in ./ler_data_custom/unlensed_param_detectable.json
Using provided lensed_param dict...
total lensed rate (yr^-1): 146.82448281660734
number of simulated lensed detectable events: 1952
number of simulated all lensed events: 20000
storing detectable params in ./ler_data_custom/lensed_param_detectable.json
unlensed_rate: 118984.60266071376
lensed_rate: 146.82448281660734
ratio: 810.3866628928141

=== Detection Rates (BNS with Custom Configuration) ===
Unlensed rate: 1.1898e+05 events/year
Lensed rate:   1.4682e+02 events/year
Rate ratio (lensed/unlensed): 810.3867

4. Compare Custom vs Internal Models

4.1. Mass Distribution Comparison

Compare the custom uniform mass distribution with an internal BNS mass model (normal distribution).

[27]:

import ler.utils as lerplt
import matplotlib.pyplot as plt

# Custom mass distributions for comparison (uniform)
m1_custom = source_frame_masses_uniform(size=5000)
m2_custom = source_frame_masses_uniform(size=5000)
# swapping to ensure mass_1 >= mass_2
idx = m1_custom < m2_custom
m1_custom[idx], m2_custom[idx] = m2_custom[idx], m1_custom[idx]

# internal BNS normal distribution in LeR
kwargs = dict(
    param_name = "mass_1_source",
    sampler_type = "truncated_normal",
    mu=1.4, sigma=0.68, x_min=1.0, x_max=2.5
)
m1_default = ler.truncated_normal(size=5000, **kwargs)
m2_default = ler.truncated_normal(size=5000, **kwargs)
# swapping to ensure mass_1 >= mass_2
idx = m1_default < m2_default
m1_default[idx], m2_default[idx] = m2_default[idx], m1_default[idx]

custom_dict = dict(mass_1=m1_custom)
default_dict = dict(mass_1=m1_default)

# Plot comparison
plt.figure(figsize=(6, 4))
lerplt.param_plot(
    param_name="mass_1",
    param_dict=custom_dict, # or the json file name
    plot_label='custom model',
);
lerplt.param_plot(
    param_name="mass_1",
    param_dict=default_dict,
    plot_label='default model',
);
plt.xlabel(r'Primary Mass $m_1$ ($M_{\odot}$)', fontsize=11)
plt.ylabel(r'Probability Density', fontsize=11)
plt.title('BNS Mass Distribution (Source Frame)', fontsize=13, fontweight='bold')
plt.xlim(1.0, 2.3)
plt.grid(alpha=0.3)
plt.legend(fontsize=10)
plt.tight_layout()
plt.show()
../_images/examples_LeR_custom_functions_41_0.png

4.2. Axis-Ratio Distribution Comparison

Compare the ler default axis-ratio distribution (Rayleigh distribution, Collet et al. 2018) with custom axis-ratio distribution (Padilla and Strauss 2008)

[28]:

# Compare axis-ratio distributions
size = 5000
zl = np.ones(size)
sigma = ler.velocity_dispersion(size, zl)
rayleigh = ler.rayleigh(size, sigma)
padilla_strauss = q_object(size)

axis_ratio_dict = dict(rayleigh=rayleigh, padilla_strauss=padilla_strauss)

# Plot comparison
plt.figure(figsize=(6, 4))
lerplt.param_plot(
    param_name="padilla_strauss",
    param_dict=axis_ratio_dict,
    plot_label='padilla_strauss',
)
lerplt.param_plot(
    param_name="rayleigh",
    param_dict=axis_ratio_dict,
    plot_label='rayleigh',
)
plt.xlabel(r'Axis Ratio $q$', fontsize=11)
plt.ylabel(r'Probability Density', fontsize=11)
plt.title('Lens Axis Ratio Distribution: Custom vs Default', fontsize=13, fontweight='bold')
plt.xlim(0.2, 1.)
plt.grid(alpha=0.3)
plt.legend(fontsize=10)
plt.tight_layout()
plt.show()

rayleigh interpolator will be generated at ./interpolator_json/axis_ratio/rayleigh_2.json
../_images/examples_LeR_custom_functions_43_1.png

5. Summary

This notebook demonstrated how to customize the LeR class for gravitational wave population studies:

  • Custom source priors: Replaced the default BNS bimodal mass distribution with a uniform distribution

  • Custom merger rate density: Modified the local merger rate density parameter

  • Custom lens parameters: Implemented Padilla & Strauss (2008) axis-ratio distribution and modified velocity dispersion

  • Custom detection criteria: Used 3G detectors (ET, CE) with SNR threshold of 12

These customizations enable flexible modeling of different astrophysical scenarios while leveraging LeR’s efficient sampling and rate calculation infrastructure.