diff --git a/slsim/Deflectors/DeflectorPopulation/generate_cluster_deflectors.py b/slsim/Deflectors/DeflectorPopulation/generate_cluster_deflectors.py new file mode 100644 index 000000000..4567031ea --- /dev/null +++ b/slsim/Deflectors/DeflectorPopulation/generate_cluster_deflectors.py @@ -0,0 +1,455 @@ +"""The GeneratedDeflector class implements the analytic subhalo model from Han +et al. (2016) Some other choices are made following Abe et al. (2025) + +Model summary: + +1. First, the accretion subhalo masses are sampled from the SHMF given in Han et al. (2018), table 1, fit to the Millenium 2 simulation +The SHMF given in the original paper is inaccurate for higher subhalo masses + +2. Host halo concentration is set by the Diemer19 mass-concentration relation with log scatter 0.33 +The accreted subhalos are placed randomly in the host halo (within R200) tracing the host halo density + +3. The evolved subhalo masses are calculated according to eq. 7 in Han 2016, with some subhalos being completely stripped +Subhalo concentration is chosen the same way using the accretion mass +Subhalos use a truncated NFW profile with truncation radius from eq. 6 in Gilman et al. (2016) + +4. Galaxy masses are calculated from subhalo accretion mass from the Behroozi SMHM relation, with log scatter 0.2 + +5. Galaxies are randomly selected as red or blue according to their distance to the center. Hennig et al. (2017) describes the density of galaxies +in the cluster as NFW profiles with red galaxies being more concentrated in the center, as well as the average fraction of red galaxies +which is a function of host mass and redshift. + +6. The closest matching red / blue galaxy in terms of stellar mass and redshift is then selected from the skypy catalog + +TODO better truncation radius, different subhalo profile, eccentricity for halos + +References: +Han et al. (2016): https://academic.oup.com/mnras/article/457/2/1208/965286 +Han et al. (2018): https://ui.adsabs.harvard.edu/abs/2018MNRAS.474..604H/abstract +Abe et al. (2025): https://arxiv.org/abs/2411.07509 +Gilman et al. (2019): https://arxiv.org/abs/1908.06983 +Hennig et al. (2017): https://ui.adsabs.harvard.edu/abs/2017MNRAS.467.4015H/abstract +""" + +from colossus.halo import profile_nfw +from colossus.halo import concentration + +import numpy as np + +# from scipy.stats import truncnorm + +from slsim.Deflectors.MassLightConnection.galaxy_population import ( + gals_init, + stellarmass_halomass, +) +from slsim.Deflectors.deflector_util import set_colossus_cosmo + +from slsim.Deflectors.deflector_group import DeflectorGroup + +from slsim.Sources.SourcePopulation.galaxies import galaxy_projected_eccentricity + +from slsim.Halos.halo_population import dNhalodzdlnM_lens + +import random + + +class GeneratedDeflector: + def __init__(self, M200h, z, cosmo): + """ + :param M: host halo mass (200c) in solar masses / h + :param z: deflector redshift + :param cosmo: astropy cosmology instance + """ + + set_colossus_cosmo(cosmo) + + c = concentration.concentration(M200h, "200c", z, model="diemer19") + c = np.random.lognormal( + mean=np.log(c), sigma=0.33 + ) # concentration distribution 0.16 dex + self.p_nfw = profile_nfw.NFWProfile(M=M200h, c=c, z=z, mdef="200c") + + self.M200h = M200h # solar masses / h + self.R200h = self.p_nfw.RDelta(z, "200c") # kpc / h + + self.z = z + self.cosmo = cosmo + + # fraction of red galaxies, function of redshift and mass + f_red_200 = ( + 0.68 + * (self.M200h / self.cosmo.h / (6 * 10**14)) ** -0.10 + * ((1 + self.z) / (1 + 0.46)) ** -0.65 + ) + self.f_red_200 = np.random.normal(loc=f_red_200, scale=0.14 * f_red_200) + + def accreted_subhalo_density_pdf(self, m_acc_h): + """Returns accreted subhalo dN / dlnm. + + :param m_acc: subhalo accretion mass [solar masses / h] + """ + + # Han et al. 2018, table 1 + a1 = 0.11 + al1 = 0.95 + a2 = 0.20 + al2 = 0.30 + b = 7.6 + beta = 2.1 + + mu = m_acc_h / self.M200h + + return (a1 * mu**-al1 + a2 * mu**-al2) * np.exp(-b * mu**beta) + + def unevolved_spatial_distribution_pdf(self, r): + """ + r: kpc/h + + In the model, when selecting by accretion mass the final position traces host halo density + This function returns normalized dN / d3r as a function of radius, so the total probability over the host halo is 1 + """ + + return self.p_nfw.density(r) / self.M200h + + def evolved_subhalo_mass_fraction_of_accretion(self, r): + """Returns m_evolved / m_acc. + + :param r: distance to halo center [kpc / h] + """ + + M = self.M200h / 10**10 + + mustar = 0.5 * M**-0.03 # stripping function amplitude + beta = 1.7 * M**-0.04 # stripping function slope + sigma = 1.1 # log scatter of m / m_acc + fs = 0.55 # fraction of survived subhaloes + + mubar = mustar * (r / self.R200h) ** beta # average evolved mass fraction + mu = np.random.lognormal(mean=np.log(mubar), sigma=sigma) # with scatter + + return np.where(np.random.rand(*r.shape) < fs, mu, 0) # 0 mass when stripped + + def generate_subhalos(self, m_acc_min, m_acc_max): + """Randomly samples subhalos from m_acc_min to m_acc_max (solar masses + / h) + + Returns list of tuples (accretion mass [solar masses / h], + evolved mass [solar masses / h], distance to center [kpc / h]) + """ + + m_acc = np.logspace(np.log10(m_acc_min), np.log10(m_acc_max), num=1000, base=10) + r = self.R200h * np.linspace(0.001, 1, 1000) + + d3r = 4 / 3 * np.pi * (r[1:] ** 3 - r[:-1] ** 3) + + # probability that the accreted subhalo ended up at each radius (traces halo density profile) + spatial_distribution_multiplier = ( + self.unevolved_spatial_distribution_pdf(r[:-1]) * d3r + ) + + dlnm = np.log(m_acc[1]) - np.log(m_acc[0]) + + # mass_acc, mass_evolved, radius + subhalos = [] + + for m in m_acc: + # expected number of accreted subhalos of this mass + accreted_subhalo_expected = self.accreted_subhalo_density_pdf(m) * dlnm + + # function of radius, expected number of subhalos with this accretion mass at this radius + probability_present = ( + spatial_distribution_multiplier * accreted_subhalo_expected + ) + + # poisson unnecessary + present_r_indices = np.where( + np.random.rand(*probability_present.shape) < probability_present + )[0] + + mass_acc = [m] * len(present_r_indices) + mass_evolved = m * self.evolved_subhalo_mass_fraction_of_accretion( + r[present_r_indices] + ) + radius = r[present_r_indices] + + subhalos += list(zip(mass_acc, mass_evolved, radius)) + + return subhalos + + def get_deflector_data( + self, + red_galaxies, + blue_galaxies, + min_subhalo_accretion_mass=None, + crop_subhalo_dist=1000, + ): + """ + :param galaxy_list: list of galaxies to assign as deflectors + :param min_subhalo_accretion_mass: min subhalo accretion mass [solar masses / h] + :param crop_subhalo_dist: only return subhalos closer than this distance to the BCG [arcsecs] + + returns kwargs_mass_list, kwargs_light_list, center_x_deflector_list, center_y_deflector_list to be passed into DeflectorGroup + """ + + if min_subhalo_accretion_mass is None: + min_subhalo_accretion_mass = 10**-3 * self.M200h + + paramc, params = gals_init() + angular_diameter_dist = self.cosmo.angular_diameter_distance(self.z).to_value( + "kpc" + ) # distance corresponding to 1 radian + + mean_position_angle = np.random.rand() * 2 * np.pi + + light_dicts = [] + mass_dicts = [] + center_x_list = [] + center_y_list = [] + + log_skypy_red_galaxy_masses = np.log10(red_galaxies["stellar_mass"]) + log_skypy_blue_galaxy_masses = np.log10(blue_galaxies["stellar_mass"]) + + halos = [(self.M200h, self.M200h, 0)] + self.generate_subhalos( + min_subhalo_accretion_mass, self.M200h + ) + for subhalo_m_acc, subhalo_m_evolved, dist_to_center in halos: + # Mo/h, Mo/h, kpc/h + + # place randomly in 3d and cast to 2d + pos = np.random.normal(size=(3)) + pos /= np.linalg.norm(pos) + pos_2d = ( + pos[:2] * dist_to_center / self.cosmo.h / angular_diameter_dist * 206265 + ) # 2d coordinate in arcsecs + + if not ( + (-crop_subhalo_dist < pos_2d[0] < crop_subhalo_dist) + and (-crop_subhalo_dist < pos_2d[1] < crop_subhalo_dist) + ): + continue + + # compute stellar mass from halo mass, use paramc if host halo, else params for SMHM relation + galaxy_mass = stellarmass_halomass( + subhalo_m_acc, self.z, paramc if dist_to_center == 0 else params + ) # solar masses / h + galaxy_mass = np.random.lognormal( + mean=np.log(galaxy_mass), sigma=0.2 + ) # scatter + + # TODO remove weighing, make skypy generate more massive galaxies + if np.random.rand() < self.fraction_red_galaxies(dist_to_center): # red + # Find sample galaxy with closest redshift and stellar mass #skypy returns in physical mass + # weigh redshift higher + closest_real_galaxy_index = np.argmin( + np.hypot( + 5 * (red_galaxies["z"] - self.z), + log_skypy_red_galaxy_masses + - np.log10(galaxy_mass / self.cosmo.h), + ) + ) + + light_dict = dict(red_galaxies[closest_real_galaxy_index]) + else: # blue + # Find sample galaxy with closest redshift and stellar mass #skypy returns in physical mass + # weigh redshift higher + closest_real_galaxy_index = np.argmin( + np.hypot( + 5 * (blue_galaxies["z"] - self.z), + log_skypy_blue_galaxy_masses + - np.log10(galaxy_mass / self.cosmo.h), + ) + ) + + light_dict = dict(blue_galaxies[closest_real_galaxy_index]) + + del light_dict["z"] + light_dict["extended_source_type"] = "hernquist" + + # eccentricity + light_dict["e1"], light_dict["e2"] = galaxy_projected_eccentricity( + light_dict["ellipticity"], + np.random.normal(loc=mean_position_angle, scale=35.4 * np.pi / 180), + ) + + if subhalo_m_evolved > 0: + # use accretion mass for concentration + c = concentration.concentration( + subhalo_m_acc, "200c", self.z, model="diemer19" + ) + c = np.random.lognormal( + mean=np.log(c), sigma=0.33 + ) # concentration distribution 0.16 dex + + # slsim wants physical masses, not / h + mass_dict = { + "mass_type": "NFW_HERNQUIST", + "halo_mass": subhalo_m_evolved / self.cosmo.h, + "concentration": c, + "e1": 0, + "e2": 0, + } + + if dist_to_center != 0: # subhalo + mass_dict["truncation_radius"] = ( + 1.4 + * (subhalo_m_evolved / self.cosmo.h / 10**7) ** (1 / 3) + * (dist_to_center / self.cosmo.h / 50) ** (2 / 3) + ) + else: # subhalo is completely stripped + mass_dict = {"mass_type": "HERNQUIST"} + + light_dicts.append(light_dict) + mass_dicts.append(mass_dict) + center_x_list.append(pos_2d[0]) + center_y_list.append(pos_2d[1]) + + return { + "kwargs_mass_list": mass_dicts, + "kwargs_light_list": light_dicts, + "center_x_deflector_list": center_x_list, + "center_y_deflector_list": center_y_list, + } + + def get_deflector( + self, + red_galaxies, + blue_galaxies, + min_subhalo_accretion_mass=None, + crop_subhalo_dist=1000, + ): + """ + :param galaxy_list: list of galaxies to assign as deflectors + :param min_subhalo_accretion_mass: min subhalo accretion mass [solar masses / h] + :param crop_subhalo_dist: only return subhalos closer than this distance to the BCG [arcsecs] + + returns DeflectorGroup + """ + + return DeflectorGroup( + self.z, + **self.get_deflector_data( + red_galaxies, + blue_galaxies, + min_subhalo_accretion_mass, + crop_subhalo_dist, + ) + ) + + def fraction_red_galaxies(self, r): + """Returns probability that a galaxy is red at a distance r (kpc / h) + + Red and blue galaxy densities in the cluster are modeled as NFW + profiles with different concentrations + + f_red_200 is calculated earlier as a function of mass and + redshift + + Data taken from Hennig et al. (2017) + """ + + # From f_red_200, we need to calculate the NFW densities describing the galaxy distributions from their concentrations + + c_red = 5.37 + c_blue = 1.38 + + def nfw(x, c): + y = c * x + y = max(y, 0.0001) + return 1 / (y * (1 + y) ** 2) + + def enclosed(y): + y = max(y, 0.0001) + return np.log(1 + y) - y / (1 + y) + + I_red = enclosed(c_red) / c_red**3 + I_blue = enclosed(c_blue) / c_blue**3 + + A_red = self.f_red_200 / I_red + A_blue = (1 - self.f_red_200) / I_blue + + rho_red = A_red * nfw(r / self.R200h, c_red) + rho_blue = A_blue * nfw(r / self.R200h, c_blue) + + return rho_red / (rho_red + rho_blue) + + +class GeneratedDeflectorPopulation: + def __init__( + self, + Mmin, + Mmax, + zmin, + zmax, + red_galaxies, + blue_galaxies, + sky_area, + cosmo, + crop_subhalo_dist=100, + min_subhalo_accretion_mass=None, + ): + """ + :param Mmin: Minimum host halo mass (solar masses / h) + :param Mmax: Maximum host halo mass (solar masses / h) + :param zmin: min redshift + :param zmax: max redshift + :param red_galaxies: list of red galaxies from skypy + :param blue_galaxies: list of blue galaxies from skypy + :param sky_area: sky area over which to draw deflectors, set very high to get a realistic population + :type sky_area: astropy.Quantity + :param cosmo: astropy cosmology instance + :param crop_subhalo_dist: maximum distance in arcsecs from the center that satellites will be generated + :param min_subhalo_accretion_mass: minimum mass of accreted subhalos to sample (solar masses / h). Default 10^-3 M_hh, but can be set explicitly if generating smaller groups + """ + + self.Mmin = Mmin + self.Mmax = Mmax + self.zmin = zmin + self.zmax = zmax + self.red_galaxies = red_galaxies + self.blue_galaxies = blue_galaxies + self.crop_subhalo_dist = crop_subhalo_dist + self.min_subhalo_accretion_mass = min_subhalo_accretion_mass + self.cosmo = cosmo + self.sky_area = sky_area + + self.deflectors = [] + + cosmo_col = set_colossus_cosmo(cosmo) + + MM_h = np.logspace(np.log10(Mmin), np.log10(Mmax), 1000) + dlnm = np.log(MM_h[1]) - np.log(MM_h[0]) + + dz = 0.001 + for z in np.arange(zmin, zmax, dz): + counts = np.random.poisson( + dNhalodzdlnM_lens(MM_h, z, cosmo_col, mdef="200c", model="tinker08") + * sky_area.to_value("deg2") + * dlnm + * dz + ) + + for m, count in zip(MM_h, counts): + if count > 0: + self.deflectors += [(m, z)] * count + + def draw_deflector(self): + """Draw a random deflector. + + Returns DeflectorGroup object + """ + + M, z = random.choice(self.deflectors) + + # print(f"Drew M {np.log10(M):.4f}/h z {z:.4f}") + + generated_deflector = GeneratedDeflector(M, z, self.cosmo) + + return generated_deflector.get_deflector( + self.red_galaxies, + self.blue_galaxies, + crop_subhalo_dist=self.crop_subhalo_dist, + min_subhalo_accretion_mass=self.min_subhalo_accretion_mass, + ) + + def deflector_number(self): + return len(self.deflectors) diff --git a/slsim/Deflectors/MassTypes/nfw.py b/slsim/Deflectors/MassTypes/nfw.py index 64efaa92b..78f075386 100644 --- a/slsim/Deflectors/MassTypes/nfw.py +++ b/slsim/Deflectors/MassTypes/nfw.py @@ -6,7 +6,16 @@ class NFW(MassBase): """Class of a NFW lens model.""" - def __init__(self, light, halo_mass, concentration, e1=0, e2=0, vel_disp=None): + def __init__( + self, + light, + halo_mass, + concentration, + e1=0, + e2=0, + vel_disp=None, + truncation_radius=None, + ): """ :param light: light model (used for position of deflector and stellar mass density profile) @@ -18,11 +27,13 @@ def __init__(self, light, halo_mass, concentration, e1=0, e2=0, vel_disp=None): :param e2: halo eccentricity component 2 :param vel_disp: velocity dispersion [km/s], optional as pre-computed value. ATTENTION: consistency is not checked with mass profile. + :param truncation_radius: if not None, uses the truncated nfw profile with this truncation radius """ super().__init__(light=light, vel_disp=vel_disp) self._halo_mass = halo_mass self._concentration = concentration self._e1_mass, self._e2_mass = e1, e2 + self._truncation_radius = truncation_radius def velocity_dispersion(self, cosmo=None): """Velocity dispersion of deflector. Simplified assumptions on @@ -56,7 +67,9 @@ def mass_model_lenstronomy(self, lens_cosmo, spherical=False): else: _spherical = False - if _spherical is True: + if self._truncation_radius is not None: + lens_mass_model_list = ["TNFW"] + elif _spherical is True: lens_mass_model_list = ["NFW"] else: lens_mass_model_list = ["NFW_ELLIPSE_CSE"] @@ -74,6 +87,11 @@ def mass_model_lenstronomy(self, lens_cosmo, spherical=False): "center_y": center_y, } ] + if self._truncation_radius is not None: + kwargs_lens_mass[0]["r_trunc"] = lens_cosmo.phys2arcsec_lens( + self._truncation_radius / 1000 + ) # function converts mpc to arcsec + if _spherical is False: e1_mass_lenstronomy, e2_mass_lenstronomy = ellipticity_slsim_to_lenstronomy( e1_slsim=self._e1_mass, e2_slsim=self._e2_mass diff --git a/slsim/Deflectors/MassTypes/nfw_hernquist.py b/slsim/Deflectors/MassTypes/nfw_hernquist.py index d0b7840a8..b166a6256 100644 --- a/slsim/Deflectors/MassTypes/nfw_hernquist.py +++ b/slsim/Deflectors/MassTypes/nfw_hernquist.py @@ -10,7 +10,16 @@ class NFWHernquist(MassBase): """Class of a NFW+Hernquist lens model with a Hernquist light mode.""" - def __init__(self, light, halo_mass, concentration, e1=0, e2=0, vel_disp=None): + def __init__( + self, + light, + halo_mass, + concentration, + e1=0, + e2=0, + vel_disp=None, + truncation_radius=None, + ): """ :param light: light model (used for position of deflector and stellar mass density profile) @@ -22,10 +31,16 @@ def __init__(self, light, halo_mass, concentration, e1=0, e2=0, vel_disp=None): :param e2: halo eccentricity component 2 :param vel_disp: velocity dispersion [km/s], optional as pre-computed value. ATTENTION: consistency is not checked with mass profile. + :param truncation_radius: if not None, uses the truncated NFW profile with this truncation radius """ super().__init__(light=light, vel_disp=vel_disp) self._nfw = NFW( - light=light, halo_mass=halo_mass, concentration=concentration, e1=e1, e2=e2 + light=light, + halo_mass=halo_mass, + concentration=concentration, + e1=e1, + e2=e2, + truncation_radius=truncation_radius, ) self._hernquist = Hernquist(light=light) self.num_mass_models = 2 diff --git a/slsim/Deflectors/deflector_util.py b/slsim/Deflectors/deflector_util.py index e37d6b67b..4f9a8ce57 100644 --- a/slsim/Deflectors/deflector_util.py +++ b/slsim/Deflectors/deflector_util.py @@ -6,6 +6,9 @@ from slsim.Util import param_util import numpy as np from colossus.cosmology import cosmology as colossus_cosmo +from lenstronomy.LensModel.lens_model import LensModel +from lenstronomy.LensModel.lens_model_extensions import LensModelExtensions +from lenstronomy.Cosmo.lens_cosmo import LensCosmo def deflector_from_table(table, mass_type, extended_source_type, cosmo=None): @@ -191,4 +194,46 @@ def set_colossus_cosmo(cosmo): sigma8=0.8102, ns=0.9660499, ) - colossus_cosmo.setCosmology(cosmo_name="halo_cosmo", **params) + return colossus_cosmo.setCosmology(cosmo_name="halo_cosmo", **params) + + +def critical_curves_caustics_list( + deflector, z_source, cosmo, kwargs_critical_curve_caustics=None +): + """Returns list of critical curves and caustics for a source at `z_source` + + :param deflector: `DeflectorGroup` or `Deflector` object + :param z_source: redshift at which to compute curves + :param cosmo: astropy.cosmology instance + :param kwargs_critical_curve_caustics: arguments passed into the `critical_curve_caustics` function from LensModelExtensions. Keys: + - compute_window: window size in arcsec where the critical curve is computed + - grid_scale: numerical grid spacing of the computation of the critical curves + - center_x: float, center of the window to compute critical curves and caustics + - center_y: float, center of the window to compute critical curves and caustics + - kwargs_lens: lens model kwargs + """ + + lens_cosmo = LensCosmo( + z_lens=deflector.redshift, + z_source=z_source, + cosmo=cosmo, + ) + + lens_mass_model_list, model_params = deflector.mass_model_lenstronomy(lens_cosmo) + + lens_model = LensModel( + lens_model_list=lens_mass_model_list, + cosmo=cosmo, + z_lens=deflector.redshift, + z_source=z_source, + multi_plane=False, + ) + + lens_model_ext = LensModelExtensions(lens_model) + ra_crit_list, dec_crit_list, ra_caustic_list, dec_caustic_list = ( + lens_model_ext.critical_curve_caustics( + model_params, **kwargs_critical_curve_caustics + ) + ) + + return ra_crit_list, dec_crit_list, ra_caustic_list, dec_caustic_list diff --git a/slsim/Halos/halo_population.py b/slsim/Halos/halo_population.py index f5ead0613..e51efd6c8 100644 --- a/slsim/Halos/halo_population.py +++ b/slsim/Halos/halo_population.py @@ -79,11 +79,11 @@ def calc_vol(z, cosmo_col): return (dis * dis / 3282.806350011744) * drdz * (1.0 + z) * (1.0 + z) * (1.0 + z) -def dNhalodzdlnM_lens(M, z, cosmo_col): +def dNhalodzdlnM_lens(M, z, cosmo_col, **massfunction_kwargs): """Compute the differential number density of halos with respect to redshift and log halo mass, per a unit of solid angle [deg^2] - :param M: The masses of the dark matter halos. + :param M: The masses of the dark matter halos (solar masses / h). :type M: ndarray, float :param z: The redshift at which to compute the differential number density. @@ -94,11 +94,14 @@ def dNhalodzdlnM_lens(M, z, cosmo_col): unit redshift per natural log mass interval per unit area in units of #/deg^2/dlnM[M_sol/h]. """ + + massfunction_params = {"mdef": "fof", "model": "sheth99"} + massfunction_params.update(**massfunction_kwargs) + dvoldzdO = calc_vol(z, cosmo_col) hhh = (cosmo_col.H0 / 100.0) ** 3 mfunc_so = ( - mass_function.massFunction(M, z, mdef="fof", model="sheth99", q_out="dndlnM") - * hhh + mass_function.massFunction(M, z, q_out="dndlnM", **massfunction_params) * hhh ) return dvoldzdO * mfunc_so diff --git a/slsim/Lenses/LensPopulation/lens_pop.py b/slsim/Lenses/LensPopulation/lens_pop.py index bd2a03d9c..5602d07bd 100644 --- a/slsim/Lenses/LensPopulation/lens_pop.py +++ b/slsim/Lenses/LensPopulation/lens_pop.py @@ -8,7 +8,8 @@ from slsim.LOS.los_pop import LOSPop from slsim.Deflectors.DeflectorPopulation.deflectors_base import DeflectorsBase from slsim.Lenses.LensPopulation.lensed_population_base import LensedPopulationBase - +from matplotlib.path import Path +from slsim.Deflectors.deflector_util import critical_curves_caustics_list from tqdm import tqdm @@ -107,6 +108,132 @@ def select_lens_at_random(self, test_area=None, verbose=False, **kwargs_lens_cut return gg_lens n += 1 + def select_lens_at_random_multi_source( + self, + source_area, + verbose=False, + min_num_sources=1, + return_only_multiply_imaged_sources=False, + **kwargs_lens_cut + ): + """Draw a random lens with an entire source field spanning + `source_area`, with at least min_num_sources satisfying + kwargs_lens_cut. + + :param sky_area: Sky area to draw sources from + :type sky_area: (astropy.units.Quantity) + :param kwargs_lens_cut: Dictionary of cuts that one wants to apply to the lens. + Example: kwargs_lens_cut = { + "min_image_separation": 0.5, + "max_image_separation": 10, + "mag_arc_limit": {"i": 24}, + "second_brightest_image_cut": {"i": 24}}. + All these cuts are optional. + :type kwargs_lens_cut: dict + :param min_num_sources: Minimum number of lensed sources that must satisfy kwargs_lens_cut + :type min_num_sources: int + :param return_only_multiply_imaged_sources: Return a lens class only with sources that are multiply imaged. Default False. + :type return_only_multiply_imaged_sources: bool + :param verbose: print statements added + :type verbose: bool + :return: Lens() instance with parameters of the deflector and lens and source field. + :rtype: Lens + """ + + # utility function - returns True if point is inside one of the paths in the ra and dec lists, or within 3 arcseconds outside. Used to cut sources for image checking. + def _in_caustic_or_close_outside(ra_caustic_list, dec_caustic_list, s): + point = s.extended_source_position + + for ra, dec in zip(ra_caustic_list, dec_caustic_list): + curve_path = Path(np.column_stack((ra, dec))) + + area = (np.max(ra) - np.min(ra)) * (np.max(dec) - np.min(dec)) + if area < 1: + continue + + if curve_path.contains_point(point): + return True + + if np.min((ra - point[0]) ** 2 + (dec - point[1]) ** 2) < 3**2: + return True + + n = 0 + while True: + n += 1 + + # draw random deflector + _deflector = self._lens_galaxies.draw_deflector() + + ### compute caustics at high redshift to filter source galaxies for validity checking + _, _, ra_caustic_list, dec_caustic_list = critical_curves_caustics_list( + _deflector, + 10, + self.cosmo, + { + "compute_window": np.sqrt(source_area.to_value("arcsec2") / np.pi) + * 2, + "grid_scale": 0.5, + }, + ) + + if len(ra_caustic_list) == 0: + continue + + # draw all sources, and filter onces near caustics for validity checking + _source = self._sources.draw_galaxies(source_area) + _source_cut = [ + s + for s in _source + if s.redshift > _deflector.redshift + and _in_caustic_or_close_outside(ra_caustic_list, dec_caustic_list, s) + ] + + if len(_source_cut) < min_num_sources: + continue + + # lens only with sources near caustics to speed up validity checking + test_lens = Lens( + deflector_class=_deflector, + source_class=_source_cut, + cosmo=self.cosmo, + use_jax=self._use_jax, + multi_plane="Source", + create_field_galaxies=True, + ) + + test_res = test_lens.validity_test(**kwargs_lens_cut) + if not isinstance(test_res, dict): + test_res = {0: test_res} + + if len([x for x in test_res.values() if x]) >= min_num_sources: + if verbose is True: + print("selected lens after %s tries." % n) + + if not return_only_multiply_imaged_sources: + # final lens with all sources + return Lens( + deflector_class=_deflector, + source_class=_source, + cosmo=self.cosmo, + use_jax=self._use_jax, + multi_plane="Source", + create_field_galaxies=True, + ) + else: + # only sources that are multiply imaged + return Lens( + deflector_class=_deflector, + source_class=[ + _source_cut[i] + for i in range(len(_source_cut)) + if test_res[i] + ], + cosmo=self.cosmo, + use_jax=self._use_jax, + multi_plane="Source", + create_field_galaxies=True, + ) + def _draw_source(self, mag_arc_limit=None, magnification_limit=2, **kwargs): """Draw from source population considering some additional constraints to be fulfilled. diff --git a/slsim/Lenses/lens.py b/slsim/Lenses/lens.py index cca7eabda..9b84dd6f4 100644 --- a/slsim/Lenses/lens.py +++ b/slsim/Lenses/lens.py @@ -42,6 +42,7 @@ def __init__( shear=True, convergence=True, field_galaxies=None, + create_field_galaxies=False, ): """ @@ -74,8 +75,21 @@ def __init__( Instances should be generated via :meth:`slsim.Sources.SourcePopulation.Galaxies.draw_field_galaxies`. If None, no field galaxies are included. :type field_galaxies: list[`slsim.Sources.source.Source`] or None - + :param create_field_galaxies: If True, make all source galaxies in front of `deflector_class.redshift` into field galaxies + :type create_field_galaxies: bool """ + + if create_field_galaxies and isinstance(source_class, list): + if field_galaxies is None: + field_galaxies = [] + + field_galaxies += [ + s for s in source_class if s.redshift < deflector_class.redshift + ] + source_class = [ + s for s in source_class if s.redshift > deflector_class.redshift + ] + LensedSystemBase.__init__( self, source_class=source_class, @@ -196,7 +210,12 @@ def _image_position_from_source(self, x_source, y_source, source_index): list], [DEC list] """ lens_model_class, kwargs_lens = self.deflector_mass_model_lenstronomy( - source_index=source_index + source_index=source_index, + multi_plane=self.multi_plane + in [ + "Deflector", + "Both", + ], # unnecessary for source, makes calculations much faster ) lens_eq_solver = LensEquationSolver(lens_model_class) point_source_pos_x, point_source_pos_y = x_source, y_source @@ -226,6 +245,7 @@ def validity_test( self, min_image_separation=0, max_image_separation=10, + min_num_of_images=2, mag_arc_limit=None, second_brightest_image_cut=None, snr_limit=None, @@ -235,6 +255,7 @@ def validity_test( :param min_image_separation: minimum image separation :param max_image_separation: maximum image separation + :param min_num_of_images: minimum number of images :param mag_arc_limit: dictionary with key of bands and values of magnitude limits of integrated lensed arc :type mag_arc_limit: dict with key of bands and values of @@ -254,6 +275,7 @@ def validity_test( validity_results[index] = self._validity_test( min_image_separation=min_image_separation, max_image_separation=max_image_separation, + min_num_of_images=min_num_of_images, mag_arc_limit=mag_arc_limit, second_brightest_image_cut=second_brightest_image_cut, snr_limit=snr_limit, @@ -268,6 +290,7 @@ def _validity_test( self, min_image_separation=0, max_image_separation=10, + min_num_of_images=2, mag_arc_limit=None, second_brightest_image_cut=None, snr_limit=None, @@ -278,6 +301,7 @@ def _validity_test( :param min_image_separation: minimum image separation :param max_image_separation: maximum image separation + :param min_num_of_images: minimum number of images :param mag_arc_limit: dictionary with key of bands and values of magnitude limits of integrated lensed arc :type mag_arc_limit: dict with key of bands and values of @@ -317,9 +341,9 @@ def _validity_test( if np.sum((center_lens - center_source) ** 2) > einstein_radius**2 * 2: return False - # Criteria 4: The lensing configuration must produce at least two SL images. + # Criteria 4: The lensing configuration must produce at least min_num_of_images SL images. image_positions = self.point_source_image_positions()[source_index] - if len(image_positions[0]) < 2: + if len(image_positions[0]) < min_num_of_images: return False # Criteria 5: The maximum separation between any two image positions must be @@ -1440,7 +1464,7 @@ def lenstronomy_kwargs( return kwargs_model, kwargs_params - def deflector_mass_model_lenstronomy(self, source_index=None): + def deflector_mass_model_lenstronomy(self, source_index=None, multi_plane=None): """Returns lens model instance and parameters in lenstronomy conventions. @@ -1508,7 +1532,7 @@ def deflector_mass_model_lenstronomy(self, source_index=None): z_source=z_source, z_source_convention=self.max_redshift_source_class.redshift, use_jax=use_jax, - multi_plane=bool(self.multi_plane), + multi_plane=bool(self.multi_plane) if multi_plane is None else multi_plane, ) return lens_model, self._kwargs_lens diff --git a/slsim/Pipelines/skypy_pipeline.py b/slsim/Pipelines/skypy_pipeline.py index 187e5c4b4..37c0f931c 100644 --- a/slsim/Pipelines/skypy_pipeline.py +++ b/slsim/Pipelines/skypy_pipeline.py @@ -2,6 +2,7 @@ from skypy.pipeline import Pipeline import tempfile import slsim.Util.param_util as util +from astropy.table import vstack import astropy.units as u @@ -119,3 +120,12 @@ def red_galaxies(self): :rtype: list of dict """ return self._pipeline["red"] + + @property + def all_galaxies(self): + """Skypy pipeline for red and blue galaxies. + + :return: list of red and blue galaxies + :rtype: list of dict + """ + return vstack([self.red_galaxies, self.blue_galaxies]) diff --git a/tests/test_Deflectors/test_deflector.py b/tests/test_Deflectors/test_deflector.py index f190f7753..6fdbd04b9 100644 --- a/tests/test_Deflectors/test_deflector.py +++ b/tests/test_Deflectors/test_deflector.py @@ -2,6 +2,7 @@ import numpy.testing as npt from slsim.Deflectors.deflector import Deflector from lenstronomy.Cosmo.lens_cosmo import LensCosmo +from slsim.Deflectors.deflector_util import critical_curves_caustics_list class TestDeflector(object): @@ -208,3 +209,21 @@ def test_theta_e_when_source_infinity(self): npt.assert_almost_equal(theta_E_infinity, 1.8024, decimal=2) npt.assert_almost_equal(theta_E_infinity, theta_E_infinity_new, decimal=5) + + def test_critical_curves_caustics(self): + ra_crit_list, dec_crit_list, ra_caustic_list, dec_caustic_list = ( + critical_curves_caustics_list( + self.deflector_nfw_her, + 10, + None, + { + "compute_window": 40, + "grid_scale": 1.0, + }, + ) + ) + + assert len(ra_crit_list) > 0 + assert len(dec_crit_list) > 0 + assert len(ra_caustic_list) > 0 + assert len(dec_caustic_list) > 0 diff --git a/tests/test_Lenses/test_LensPopulation/test_lens_pop.py b/tests/test_Lenses/test_LensPopulation/test_lens_pop.py index 5d5fb7585..cf2a6cb19 100644 --- a/tests/test_Lenses/test_LensPopulation/test_lens_pop.py +++ b/tests/test_Lenses/test_LensPopulation/test_lens_pop.py @@ -13,6 +13,10 @@ CompoundLensHalosGalaxies, ) from slsim.Deflectors.DeflectorPopulation.cluster_deflectors import ClusterDeflectors +from slsim.Deflectors.DeflectorPopulation.generate_cluster_deflectors import ( + GeneratedDeflectorPopulation, +) + from slsim.Sources.SourcePopulation.galaxies import Galaxies from slsim.Sources.SourceCatalogues.SupernovaeCatalog.supernovae_sample import ( SupernovaeCatalog, @@ -252,6 +256,57 @@ def test_cluster_lens_pop_instance(): assert pes_lens_class.deflector_velocity_dispersion() > 100 +def test_cluster_lens_pop_instance_multi_source(): + sky_area = Quantity(value=0.05, unit="deg2") + + # one of the other tests adds a column to the galaxies table, which causes a bug when creating a Hernquist source + _galaxy_simulation_pipeline = pipelines.SkyPyPipeline( + skypy_config=None, + sky_area=sky_area, + filters=None, + ) + + cosmo = FlatLambdaCDM(H0=70, Om0=0.3) + + kwargs_source_cut = {"band": "g", "band_max": 26, "z_min": 0.21, "z_max": 5.0} + + source_galaxies = Galaxies( + galaxy_list=_galaxy_simulation_pipeline.all_galaxies, + kwargs_cut=kwargs_source_cut, + cosmo=cosmo, + sky_area=sky_area, + catalog_type="skypy", + ) + + deflectors = GeneratedDeflectorPopulation( + 10**14.5, + 10**15, + 0.3, + 0.4, + _galaxy_simulation_pipeline.red_galaxies, + _galaxy_simulation_pipeline.blue_galaxies, + Quantity(100, "deg2"), + cosmo, + ) + + lenspop = LensPop( + deflector_population=deflectors, + source_population=source_galaxies, + cosmo=cosmo, + sky_area=sky_area, + ) + + kwargs_lens_cut_plot = { + "min_image_separation": 1.0, + "max_image_separation": 100.0, + } + + pes_lens_class = lenspop.select_lens_at_random_multi_source( + source_area=Quantity(40**2, "arcsec2"), **kwargs_lens_cut_plot + ) + assert isinstance(pes_lens_class, Lens) + + def test_galaxies_lens_pop_instance(): cosmo = FlatLambdaCDM(H0=70, Om0=0.3) sky_area = Quantity(value=0.001, unit="deg2")