Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Compatibility update with GNFW #44

Merged
merged 9 commits into from
Sep 13, 2024
Merged

Compatibility update with GNFW #44

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
944ef77
Merge pull request #5 from sibirrer/main
williyamshoe Apr 24, 2024
68e2239
Merge pull request #6 from williyamshoe/wsheu_dev
williyamshoe Apr 24, 2024
0e163bb
Merge pull request #7 from sibirrer/main
williyamshoe Sep 12, 2024
4bace44
added functionality to input alpha_Rs for GNFW profile
williyamshoe Sep 12, 2024
2bcd21a
added test cases
williyamshoe Sep 12, 2024
0ca8de4
added docstrings
williyamshoe Sep 12, 2024
6d4ef83
formatting changes
williyamshoe Sep 12, 2024
6e5e0f3
changed alpha_Rs_or_kappa_s to halo_normalization
williyamshoe Sep 12, 2024
20c2f76
changed _alpha_Rs_or_kappa_s_bool to _is_normalization_alpha_Rs
williyamshoe Sep 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 49 additions & 19 deletions hierarc/LensPosterior/kin_constraints_composite.py
View file
Open in desktop
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
Loading