Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/main'
Browse files Browse the repository at this point in the history
  • Loading branch information
@sibirrer
sibirrer committed Sep 16, 2024
2 parents a73eb84 + 5bc9cac commit 5a1cb4c
Show file tree
Hide file tree
Showing 2 changed files with 276 additions and 19 deletions.
68 changes: 49 additions & 19 deletions hierarc/LensPosterior/kin_constraints_composite.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
from lenstronomy.Util import constants as const
from lenstronomy.Analysis.light_profile import LightProfileAnalysis
from lenstronomy.LightModel.light_model import LightModel
from lenstronomy.LensModel.Profiles.gnfw import GNFW


class KinConstraintsComposite(KinConstraints):
Expand All @@ -16,7 +17,7 @@ def __init__(
z_source,
gamma_in_array,
log_m2l_array,
kappa_s_array,
alpha_Rs_array,
r_s_angle_array,
theta_E,
theta_E_error,
Expand All @@ -42,6 +43,7 @@ def __init__(
num_kin_sampling=1000,
multi_observations=False,
rho0_array=None,
kappa_s_array=None,
r_s_array=None,
is_m2l_population_level=True,
):
Expand All @@ -53,7 +55,7 @@ def __init__(
:param log_m2l_array: array of log10(mass-to-light ratios) of the stellar component,
needs to be in the unit/scaling such that m2l / sigma_crit * amp in the
kwargs_lens_light provides the convergence amplitude of the stars
:param kappa_s_array: array of generalized NFW profile's convergence normalization at the scale radius
:param alpha_Rs_array: array of the deflection (angular units) at projected Rs
:param r_s_angle_array: array of halo scale radii in arcsecond
:param theta_E: Einstein radius (in arc seconds)
:param theta_E_error: 1-sigma error on Einstein radius
Expand Down Expand Up @@ -88,6 +90,7 @@ def __init__(
:param multi_observations: bool, if True, interprets kwargs_aperture and
kwargs_seeing as lists of multiple observations
:param rho0_array: array of halo mass normalizations in M_sun / Mpc^3
:param kappa_s_array: array of generalized NFW profile's convergence normalization at the scale radius
:param r_s_array: array of halo scale radii in Mpc
"""

Expand Down Expand Up @@ -152,16 +155,23 @@ def __init__(
log_m2l_scaling=log_m2l_scaling,
)

if self._check_arrays(kappa_s_array, r_s_angle_array):
self._kappa_s_array = kappa_s_array
if self._check_arrays(alpha_Rs_array, r_s_angle_array):
self._halo_normalization_array = alpha_Rs_array
self._is_normalization_alpha_Rs = True
self._r_scale_angle_array = r_s_angle_array
elif self._check_arrays(kappa_s_array, r_s_angle_array):
self._halo_normalization_array = kappa_s_array
self._is_normalization_alpha_Rs = False
self._r_scale_angle_array = r_s_angle_array
elif self._check_arrays(rho0_array, r_s_array):
self._kappa_s_array, self._r_scale_angle_array = self.get_kappa_s_r_s_angle(
kappa_s_array, self._r_scale_angle_array = self.get_kappa_s_r_s_angle(
rho0_array, r_s_array
)
self._halo_normalization_array = kappa_s_array
self._is_normalization_alpha_Rs = False
else:
raise ValueError(
"Both kappa_s_array and r_s_angle_array, or rho0_array and r_s_array must be arrays of the same length!"
"Both alpha_Rs_array and r_s_angle_array, kappa_s_array and r_s_angle_array, or rho0_array and r_s_array must be arrays of the same length!"
)

self.gamma_in_array = gamma_in_array
Expand All @@ -172,10 +182,10 @@ def __init__(
self._gamma_in_prior_std = gamma_in_prior_std

if not is_m2l_population_level and not self._check_arrays(
self._kappa_s_array, log_m2l_array
self._halo_normalization_array, log_m2l_array
):
raise ValueError(
"log_m2l_array must have the same length as rho0_array or kappa_s_array!"
"log_m2l_array must have the same length as alpha_Rs_array, kappa_s_array, or rho0_array!"
)

@staticmethod
Expand Down Expand Up @@ -210,22 +220,24 @@ def draw_lens(self, no_error=False):
if no_error is True:
if self._is_m2l_population_level:
return (
np.mean(self._kappa_s_array),
np.mean(self._halo_normalization_array),
np.mean(self._r_scale_angle_array),
self._r_eff,
1,
)
else:
return (
np.mean(self._kappa_s_array),
np.mean(self._halo_normalization_array),
np.mean(self._r_scale_angle_array),
np.mean(self.log_m2l_array),
self._r_eff,
1,
)

random_index = np.random.randint(low=0, high=len(self._kappa_s_array))
kappa_s_draw = self._kappa_s_array[random_index]
random_index = np.random.randint(
low=0, high=len(self._halo_normalization_array)
)
halo_normalization_draw = self._halo_normalization_array[random_index]
r_scale_angle_draw = self._r_scale_angle_array[random_index]

# we make sure no negative r_eff are being sampled
Expand All @@ -235,11 +247,11 @@ def draw_lens(self, no_error=False):
r_eff_draw = delta_r_eff * self._r_eff

if self._is_m2l_population_level:
return kappa_s_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff
return halo_normalization_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff
else:
log_m2l_draw = self.log_m2l_array[random_index]
return (
kappa_s_draw,
halo_normalization_draw,
r_scale_angle_draw,
log_m2l_draw,
r_eff_draw,
Expand Down Expand Up @@ -290,8 +302,8 @@ def j_kin_draw_composite(
the mean values instead
:return: dimensionless kinematic component J() Birrer et al. 2016, 2019
"""
kappa_s_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff = self.draw_lens(
no_error=no_error
halo_normalization_draw, r_scale_angle_draw, r_eff_draw, delta_r_eff = (
self.draw_lens(no_error=no_error)
)

kwargs_lens_stars = copy.deepcopy(self._kwargs_lens_light[0])
Expand All @@ -304,11 +316,20 @@ def j_kin_draw_composite(
for kwargs in kwargs_light:
kwargs["sigma"] *= delta_r_eff

# Input is alpha_Rs
if self._is_normalization_alpha_Rs:
alpha_Rs_draw = halo_normalization_draw
# Input is kappa_s
else:
alpha_Rs_draw = GNFW().kappa_s_to_alpha_Rs(
halo_normalization_draw, r_scale_angle_draw, gamma_in
)

kwargs_lens = [
{
"Rs": r_scale_angle_draw,
"gamma_in": gamma_in,
"kappa_s": kappa_s_draw,
"alpha_Rs": alpha_Rs_draw,
"center_x": 0,
"center_y": 0,
},
Expand Down Expand Up @@ -336,7 +357,7 @@ def j_kin_draw_composite_m2l(self, kwargs_anisotropy, gamma_in, no_error=False):
:return: dimensionless kinematic component J() Birrer et al. 2016, 2019
"""
(
kappa_s_draw,
halo_normalization_draw,
r_scale_angle_draw,
log_m2l_draw,
r_eff_draw,
Expand All @@ -353,11 +374,20 @@ def j_kin_draw_composite_m2l(self, kwargs_anisotropy, gamma_in, no_error=False):
for kwargs in kwargs_light:
kwargs["sigma"] *= delta_r_eff

# Input is alpha_Rs
if self._is_normalization_alpha_Rs:
alpha_Rs_draw = halo_normalization_draw
# Input is kappa_s
else:
alpha_Rs_draw = GNFW().kappa_s_to_alpha_Rs(
halo_normalization_draw, r_scale_angle_draw, gamma_in
)

kwargs_lens = [
{
"Rs": r_scale_angle_draw,
"gamma_in": gamma_in,
"kappa_s": kappa_s_draw,
"alpha_Rs": alpha_Rs_draw,
"center_x": 0,
"center_y": 0,
},
Expand Down
Loading

0 comments on commit 5a1cb4c

Please sign in to comment.