diff --git a/autogalaxy/config/priors/mass/total/dual_pseudo_isothermal_mass.yaml b/autogalaxy/config/priors/mass/total/dual_pseudo_isothermal_mass.yaml index 7b5298aca..9c568f16a 100644 --- a/autogalaxy/config/priors/mass/total/dual_pseudo_isothermal_mass.yaml +++ b/autogalaxy/config/priors/mass/total/dual_pseudo_isothermal_mass.yaml @@ -1,126 +1,288 @@ -dPIEMass: - centre_0: - type: Gaussian - mean: 0.0 - sigma: 0.1 - width_modifier: - type: Absolute - value: 0.05 - limits: - lower: -inf - upper: inf - centre_1: - type: Gaussian - mean: 0.0 - sigma: 0.1 - width_modifier: - type: Absolute - value: 0.05 - limits: - lower: -inf - upper: inf - ell_comps_0: - type: TruncatedGaussian - mean: 0.0 - sigma: 0.3 - lower_limit: -1.0 - upper_limit: 1.0 - width_modifier: - type: Absolute - value: 0.2 - limits: - lower: -1.0 - upper: 1.0 - ell_comps_1: - type: TruncatedGaussian - mean: 0.0 - sigma: 0.3 - lower_limit: -1.0 - upper_limit: 1.0 - width_modifier: - type: Absolute - value: 0.2 - limits: - lower: -1.0 - upper: 1.0 - ra: - type: Uniform - lower_limit: 0.0 - upper_limit: 10.0 - width_modifier: - type: Relative - value: 0.25 - limits: - lower: 0.0 - upper: inf - rs: - type: Uniform - lower_limit: 0.0 - upper_limit: 10.0 - width_modifier: - type: Relative - value: 0.25 - limits: - lower: 0.0 - upper: inf - b0: - type: Uniform - lower_limit: 0.0 - upper_limit: 10.0 - width_modifier: - type: Relative - value: 0.25 - limits: - lower: 0.0 - upper: inf -dPIEMassSph: - centre_0: - type: Gaussian - mean: 0.0 - sigma: 0.1 - width_modifier: - type: Absolute - value: 0.05 - limits: - lower: -inf - upper: inf - centre_1: - type: Gaussian - mean: 0.0 - sigma: 0.1 - width_modifier: - type: Absolute - value: 0.05 - limits: - lower: -inf - upper: inf - ra: - type: Uniform - lower_limit: 0.0 - upper_limit: 10.0 - width_modifier: - type: Relative - value: 0.25 - limits: - lower: 0.0 - upper: inf - rs: - type: Uniform - lower_limit: 0.0 - upper_limit: 10.0 - width_modifier: - type: Relative - value: 0.25 - limits: - lower: 0.0 - upper: inf - b0: - type: Uniform - lower_limit: 0.0 - upper_limit: 10.0 - width_modifier: - type: Relative - value: 0.25 - limits: - lower: 0.0 - upper: inf \ No newline at end of file +dPIEMass: + centre_0: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + centre_1: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + ell_comps_0: + type: TruncatedGaussian + mean: 0.0 + sigma: 0.3 + lower_limit: -1.0 + upper_limit: 1.0 + width_modifier: + type: Absolute + value: 0.2 + limits: + lower: -1.0 + upper: 1.0 + ell_comps_1: + type: TruncatedGaussian + mean: 0.0 + sigma: 0.3 + lower_limit: -1.0 + upper_limit: 1.0 + width_modifier: + type: Absolute + value: 0.2 + limits: + lower: -1.0 + upper: 1.0 + ra: + type: Uniform + lower_limit: 0.0 + upper_limit: 10.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf + rs: + type: Uniform + lower_limit: 0.0 + upper_limit: 10.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf + b0: + type: Uniform + lower_limit: 0.0 + upper_limit: 10.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf +dPIEMassSph: + centre_0: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + centre_1: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + ra: + type: Uniform + lower_limit: 0.0 + upper_limit: 10.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf + rs: + type: Uniform + lower_limit: 0.0 + upper_limit: 10.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf + b0: + type: Uniform + lower_limit: 0.0 + upper_limit: 10.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf +dPIEMassLenstool: + centre_0: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + centre_1: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + ellipticity: + type: Uniform + lower_limit: 0.0 + upper_limit: 0.9 + width_modifier: + type: Absolute + value: 0.2 + limits: + lower: 0.0 + upper: 1.0 + angle_pos: + type: Uniform + lower_limit: 0.0 + upper_limit: 180.0 + width_modifier: + type: Absolute + value: 30.0 + limits: + lower: -inf + upper: inf + sigma: + type: Uniform + lower_limit: 0.0 + upper_limit: 1000.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf + r_core: + type: Uniform + lower_limit: 0.0 + upper_limit: 10.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf + r_cut: + type: Uniform + lower_limit: 0.0 + upper_limit: 100.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf + redshift_object: + type: Uniform + lower_limit: 0.0 + upper_limit: 1.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: inf + redshift_source: + type: Uniform + lower_limit: 0.0 + upper_limit: 1.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: inf +dPIEMassLenstoolSph: + centre_0: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + centre_1: + type: Gaussian + mean: 0.0 + sigma: 0.1 + width_modifier: + type: Absolute + value: 0.05 + limits: + lower: -inf + upper: inf + sigma: + type: Uniform + lower_limit: 0.0 + upper_limit: 1000.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf + r_core: + type: Uniform + lower_limit: 0.0 + upper_limit: 10.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf + r_cut: + type: Uniform + lower_limit: 0.0 + upper_limit: 100.0 + width_modifier: + type: Relative + value: 0.25 + limits: + lower: 0.0 + upper: inf + redshift_object: + type: Uniform + lower_limit: 0.0 + upper_limit: 1.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: inf + redshift_source: + type: Uniform + lower_limit: 0.0 + upper_limit: 1.0 + width_modifier: + type: Relative + value: 0.5 + limits: + lower: 0.0 + upper: inf diff --git a/autogalaxy/profiles/mass/__init__.py b/autogalaxy/profiles/mass/__init__.py index 89f029b36..560a4283e 100644 --- a/autogalaxy/profiles/mass/__init__.py +++ b/autogalaxy/profiles/mass/__init__.py @@ -4,6 +4,8 @@ from .total import ( dPIEMass, dPIEMassSph, + dPIEMassLenstool, + dPIEMassLenstoolSph, PIEMass, dPIEPotential, dPIEPotentialSph, diff --git a/autogalaxy/profiles/mass/total/__init__.py b/autogalaxy/profiles/mass/total/__init__.py index 88ad4d76f..897ea34cb 100644 --- a/autogalaxy/profiles/mass/total/__init__.py +++ b/autogalaxy/profiles/mass/total/__init__.py @@ -1,5 +1,11 @@ from .dual_pseudo_isothermal_potential import dPIEPotential, dPIEPotentialSph -from .dual_pseudo_isothermal_mass import PIEMass, dPIEMass, dPIEMassSph +from .dual_pseudo_isothermal_mass import ( + PIEMass, + dPIEMass, + dPIEMassSph, + dPIEMassLenstool, + dPIEMassLenstoolSph, +) from .isothermal import Isothermal, IsothermalSph from .isothermal_core import IsothermalCore, IsothermalCoreSph from .power_law import PowerLaw, PowerLawSph diff --git a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py index 72db43108..2fd78dca0 100644 --- a/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py +++ b/autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py @@ -2,8 +2,44 @@ import numpy as np import autoarray as aa +from autogalaxy import convert +from autogalaxy.cosmology.model import LensingCosmology from autogalaxy.profiles.mass.abstract.abstract import MassProfile + +def _b0_from_lenstool_sigma( + sigma: float, + redshift_object: float, + redshift_source: float, + cosmology: LensingCosmology, +) -> float: + """ + Convert Lenstool's fiducial velocity dispersion ``v_disp`` (sigma_LT, in km/s) to the + dPIE lens strength ``b0`` in arcseconds. + + Lenstool stores ``b0 = 6 * pia_c2 * sigma_LT^2`` with ``pia_c2 = 648000 / c^2`` + (``constant.h``; c in km/s) and applies the distance ratio D_LS / D_S separately when + computing deflections (``e_grad.c``). PyAutoGalaxy's ``b0`` is fully normalized, so the + ratio is folded in here: + + b0 [arcsec] = 6 * 648000 * (sigma_LT / c)^2 * (D_LS / D_S) + + This is identical to ``4 * pi * (sigma_0 / c)^2 * D_LS / D_S`` (in radians) for the + central velocity dispersion ``sigma_0 = sqrt(3/2) * sigma_LT`` — the fiducial-vs-central + distinction of Eliasdottir et al. (2007) App. A / Bergamini et al. (2019). + """ + c_km_s = 299792.458 + + d_s = cosmology.angular_diameter_distance_to_earth_in_kpc_from( + redshift=redshift_source + ) + d_ls = cosmology.angular_diameter_distance_between_redshifts_in_kpc_from( + redshift_0=redshift_object, redshift_1=redshift_source + ) + + return 6.0 * 648000.0 * (sigma / c_km_s) ** 2 * (d_ls / d_s) + + # Within this profile family, PIEMass, dPIEMass, and dPIEMassSph are directly ported from Lenstool's C code, and have been thoroughly annotated and adapted for PyAutoLens. # The dPIEPotential and dPIEPotentialSph profiles are modified from the original `dPIEPotential` and `dPIEPotentialSph`, which were implemented to PyAutoLens by Jackson O'Donnell. @@ -215,6 +251,45 @@ def _mdci05(x, y, eps, rcore, b0, xp=np): return a, b, c, d +def _pi05(x, y, eps, rcore, xp=np): + """ + Returns the lensing potential psi / b0 of the single-core PIEMass (Kassiola & + Kovner 1993 I0.5c model) at positions (x, y), ported from Lenstool's ``pi05`` + (``e_pcpx.c``). Note b0 is outside ``_pi05``, mirroring ``_ci05``. + + The dPIE potential is the two-component difference (Lenstool ``e_pot.c`` case 81): + psi = b0 * rs / (rs - ra) * (_pi05(rcore=ra) - _pi05(rcore=rs)). + + Parameters + ---------- + eps + The ellipticity of the corresponding profiles. + rcore + The core radius of the corresponding profiles. + """ + sqe = xp.sqrt(eps) + ci = 0.5 * (1.0 - eps**2) / sqe + cxro = (1.0 + eps) ** 2 + cyro = (1.0 - eps) ** 2 + rem2 = x**2 / cxro + y**2 / cyro + e1 = 2.0 * sqe / (1.0 - eps) + e2 = 2.0 * sqe / (1.0 + eps) + z = xp.sqrt(x**2 + y**2) + + eta = -0.5 * xp.arcsinh(e1 * y / z) + 0.5j * xp.arcsin(e2 * x / z) + zeta = 0.5 * xp.log((xp.sqrt(rem2) + xp.sqrt(rcore**2 + rem2)) / rcore) + 0.0j + + b1 = xp.cosh(eta + zeta) + b2 = xp.cosh(eta - zeta) + a1 = xp.log(xp.cosh(eta) ** 2 / (b1 * b2)) + a2 = xp.log(b1 / b2) + c1 = xp.sinh(2.0 * eta) * a1 + c2 = xp.sinh(2.0 * zeta) * a2 + ckk = c1 + c2 + + return ci * rcore / xp.sqrt(rem2) * (ckk.imag * x - ckk.real * y) + + class PIEMass(MassProfile): def __init__( self, @@ -259,8 +334,11 @@ def __init__( def _ellip(self, xp=np): ellip = xp.sqrt(self.ell_comps[0] ** 2 + self.ell_comps[1] ** 2) + # The ci05 deflection integral is degenerate (NaN) at exactly zero ellipticity; + # Lenstool clamps to 1e-5 at setup for the same reason (set_lens.c). + MIN_ELLIP = 0.00001 MAX_ELLIP = 0.99999 - return xp.min(xp.array([ellip, MAX_ELLIP])) + return xp.clip(ellip, MIN_ELLIP, MAX_ELLIP) @aa.decorators.to_vector_yx @aa.decorators.transform(rotate_back=True) @@ -455,10 +533,102 @@ def __init__( self.rs = rs self.b0 = b0 + @classmethod + def from_lenstool( + cls, + centre: Tuple[float, float] = (0.0, 0.0), + ellipticity: float = 0.0, + angle_pos: float = 0.0, + sigma: float = 200.0, + r_core: float = 0.1, + r_cut: float = 20.0, + redshift_object: float = 0.5, + redshift_source: float = 1.0, + cosmology: LensingCosmology = None, + ) -> "dPIEMass": + """ + Construct a ``dPIEMass`` from Lenstool's native dPIE / PIEMD parameterization, as + read directly out of a Lenstool ``.par`` file (``potentiel`` profil 81) or the + parameter tables of Lenstool-based papers. + + Three Lenstool conventions are converted (each verified against the Lenstool C + source): + + - ``sigma`` is Lenstool's **fiducial** velocity dispersion ``v_disp`` (sigma_LT), + *not* the central velocity dispersion sigma_0 of the dPIE profile. They differ + by sigma_0 = sqrt(3/2) * sigma_LT (Eliasdottir et al. 2007, App. A; Bergamini + et al. 2019, Eq. 5) — quoting a measured central/aperture dispersion here + overestimates the mass by 50%. The lens strength is + b0 = 6 * 648000 * (sigma_LT / c)^2 * (D_LS / D_S) arcsec, where Lenstool's + stored ``b0 = 6 * pia_c2 * sigma^2`` (``set_potfile.c``) carries no distance + ratio — Lenstool applies D_LS / D_S separately at deflection time + (``e_grad.c``), whereas PyAutoGalaxy's ``b0`` is fully normalized. + - ``ellipticity`` is Lenstool's ``ellipticite`` for mass-type profiles, + emass = (a^2 - b^2) / (a^2 + b^2). Lenstool converts it internally + (``set_lens.c``) to epot = (1 - q) / (1 + q) before evaluating deflections; + that epot is exactly the magnitude of PyAutoGalaxy's ``ell_comps``, so the + conversion here is emass -> q = sqrt((1 - e) / (1 + e)) -> ``ell_comps``. + - ``r_core`` / ``r_cut`` (Lenstool ``core_radius`` / ``cut_radius``, arcsec) map + one-to-one onto ``ra`` / ``rs``. For ``.par`` files using the kpc variants + (``core_radius_kpc`` / ``cut_radius_kpc``), pre-convert with + ``r_core = r_core_kpc / cosmology.kpc_per_arcsec_from(redshift=redshift_object)``. + + Parameters + ---------- + centre + The (y,x) arc-second coordinates of the profile centre. + ellipticity + Lenstool mass ellipticity, (a^2 - b^2) / (a^2 + b^2). + angle_pos + Position angle in degrees, counter-clockwise from the positive x-axis + (Lenstool ``angle_pos`` in its tangent plane; axis flips from WCS + conventions must be handled when ingesting real-data catalogues). + sigma + Lenstool fiducial velocity dispersion ``v_disp`` (sigma_LT) in km/s. + r_core + Lenstool ``core_radius`` in arcseconds (becomes ``ra``). + r_cut + Lenstool ``cut_radius`` in arcseconds (becomes ``rs``). + redshift_object + The redshift of the lens, used for the D_LS / D_S normalization of ``b0``. + redshift_source + The redshift of the source used to normalize ``b0``. For multi-plane cluster + models this is the reference source plane the Lenstool model was normalized + to. + cosmology + The cosmology used to compute the distance ratio (default ``Planck15``; pass + the cosmology of the Lenstool run for exact comparisons). + """ + if cosmology is None: + from autogalaxy.cosmology.model import Planck15 + + cosmology = Planck15() + + axis_ratio = np.sqrt((1.0 - ellipticity) / (1.0 + ellipticity)) + ell_comps = convert.ell_comps_from(axis_ratio=axis_ratio, angle=angle_pos) + + b0 = _b0_from_lenstool_sigma( + sigma=sigma, + redshift_object=redshift_object, + redshift_source=redshift_source, + cosmology=cosmology, + ) + + return cls( + centre=centre, + ell_comps=ell_comps, + ra=r_core, + rs=r_cut, + b0=b0, + ) + def _ellip(self, xp=np): ellip = xp.sqrt(self.ell_comps[0] ** 2 + self.ell_comps[1] ** 2) + # The ci05 deflection integral is degenerate (NaN) at exactly zero ellipticity; + # Lenstool clamps to 1e-5 at setup for the same reason (set_lens.c). + MIN_ELLIP = 0.00001 MAX_ELLIP = 0.99999 - return xp.min(xp.array([ellip, MAX_ELLIP])) + return xp.clip(ellip, MIN_ELLIP, MAX_ELLIP) @aa.decorators.to_vector_yx @aa.decorators.transform(rotate_back=True) @@ -583,18 +753,44 @@ def analytical_magnification_2d_from( @aa.over_sample @aa.decorators.to_array @aa.decorators.transform + @aa.decorators.to_array + @aa.decorators.transform def potential_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): - from autogalaxy.profiles.mass.abstract.mge import MGEDecomposer - - radii_min = max(self.ra, 0.001) / 10.0 - radii_max = self.rs * 20.0 - sigmas = xp.exp(xp.linspace(xp.log(radii_min), xp.log(radii_max), 30)) - mge_decomp = MGEDecomposer(mass_profile=self) - return mge_decomp.potential_2d_via_mge_from( - grid=grid, xp=xp, sigma_log_list=sigmas, - ellipticity_convention="major", three_D=False, + """ + Returns the two dimensional projected lensing potential on a grid of (y,x) arc-second + coordinates. + + The analytic Kassiola & Kovner (1993) I0.5 potential of the dPIE is the same + two-component difference as the deflections, ported from Lenstool's C code + (``e_pot.c`` case 81, ``pi05`` in ``e_pcpx.c``): + + psi = b0 * rs / (rs - ra) * (pi05(ra) - pi05(rs)) + + Parameters + ---------- + grid + The grid of (y,x) arc-second coordinates the potential is computed on. + """ + ellip = self._ellip(xp) + factor = self.b0 * self.rs / (self.rs - self.ra) + + pot_core = _pi05( + x=grid.array[:, 1], + y=grid.array[:, 0], + eps=ellip, + rcore=self.ra, + xp=xp, + ) + pot_cut = _pi05( + x=grid.array[:, 1], + y=grid.array[:, 0], + eps=ellip, + rcore=self.rs, + xp=xp, ) + return factor * (pot_core - pot_cut) + class dPIEMassSph(dPIEMass): r"""Spherical dual pseudo-isothermal mass distribution (dPIE, mass parameterisation). @@ -681,6 +877,59 @@ def __init__( self.rs = rs self.b0 = b0 + @classmethod + def from_lenstool( + cls, + centre: Tuple[float, float] = (0.0, 0.0), + sigma: float = 200.0, + r_core: float = 0.1, + r_cut: float = 20.0, + redshift_object: float = 0.5, + redshift_source: float = 1.0, + cosmology: LensingCosmology = None, + ) -> "dPIEMassSph": + """ + Construct a ``dPIEMassSph`` from Lenstool's native dPIE / PIEMD parameterization + (circular case). See ``dPIEMass.from_lenstool`` for the full conversion + conventions; the ellipticity and angle inputs are absent here. + + Parameters + ---------- + centre + The (y,x) arc-second coordinates of the profile centre. + sigma + Lenstool fiducial velocity dispersion ``v_disp`` (sigma_LT) in km/s — not the + central velocity dispersion sigma_0 = sqrt(3/2) * sigma_LT. + r_core + Lenstool ``core_radius`` in arcseconds (becomes ``ra``). + r_cut + Lenstool ``cut_radius`` in arcseconds (becomes ``rs``). + redshift_object + The redshift of the lens, used for the D_LS / D_S normalization of ``b0``. + redshift_source + The redshift of the source used to normalize ``b0``. + cosmology + The cosmology used to compute the distance ratio (default ``Planck15``). + """ + if cosmology is None: + from autogalaxy.cosmology.model import Planck15 + + cosmology = Planck15() + + b0 = _b0_from_lenstool_sigma( + sigma=sigma, + redshift_object=redshift_object, + redshift_source=redshift_source, + cosmology=cosmology, + ) + + return cls( + centre=centre, + ra=r_core, + rs=r_cut, + b0=b0, + ) + @aa.decorators.to_vector_yx @aa.decorators.transform(rotate_back=True) def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs): @@ -805,3 +1054,138 @@ def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", xp=np, **kwargs hessian_yy = t05 * (p * Y + z * X / R2) return hessian_yy, hessian_xy, hessian_yx, hessian_xx + + +class dPIEMassLenstool(dPIEMass): + """ + The dPIE mass profile in Lenstool's native parameterization, supporting model-fitting + with priors placed directly on the Lenstool parameters. + + This is a thin wrapper around :class:`dPIEMass` whose free parameters are + (``ellipticity``, ``angle_pos``, ``sigma``, ``r_core``, ``r_cut``) exactly as they appear + in a Lenstool ``.par`` file (``potentiel`` profil 81), rather than the internal + (``ell_comps``, ``ra``, ``rs``, ``b0``). Use it to fit a model whose posteriors read + like a Lenstool results table; see ``dPIEMass.from_lenstool`` for the full + conversion conventions (verified against the Lenstool C source). + + The distance ratio D_LS / D_S entering ``b0`` uses the ``Planck15`` cosmology + (matching the ``NFWMCRLudlow`` convention); for a different cosmology construct via + ``dPIEMass.from_lenstool(..., cosmology=...)`` instead. + + Parameters + ---------- + centre : (float, float) + (y, x) arc-second coordinates of the profile centre. + ellipticity : float + Lenstool mass ellipticity ``ellipticite`` = (a^2 - b^2) / (a^2 + b^2). + angle_pos : float + Position angle in degrees, counter-clockwise from the positive x-axis. + sigma : float + Lenstool fiducial velocity dispersion ``v_disp`` (sigma_LT) in km/s — not the + central velocity dispersion sigma_0 = sqrt(3/2) * sigma_LT. + r_core : float + Lenstool ``core_radius`` in arcseconds (the dPIE ``ra``). + r_cut : float + Lenstool ``cut_radius`` in arcseconds (the dPIE ``rs``). + redshift_object : float + The redshift of the lens, used for the D_LS / D_S normalization of ``b0``. + redshift_source : float + The redshift of the source used to normalize ``b0``. + """ + + def __init__( + self, + centre: Tuple[float, float] = (0.0, 0.0), + ellipticity: float = 0.0, + angle_pos: float = 0.0, + sigma: float = 200.0, + r_core: float = 0.1, + r_cut: float = 20.0, + redshift_object: float = 0.5, + redshift_source: float = 1.0, + ): + from autogalaxy.cosmology.model import Planck15 + + cosmology = Planck15() + + axis_ratio = np.sqrt((1.0 - ellipticity) / (1.0 + ellipticity)) + ell_comps = convert.ell_comps_from(axis_ratio=axis_ratio, angle=angle_pos) + + b0 = _b0_from_lenstool_sigma( + sigma=sigma, + redshift_object=redshift_object, + redshift_source=redshift_source, + cosmology=cosmology, + ) + + super().__init__( + centre=centre, + ell_comps=ell_comps, + ra=r_core, + rs=r_cut, + b0=b0, + ) + + self.ellipticity = ellipticity + self.angle_pos = angle_pos + self.sigma = sigma + self.r_core = r_core + self.r_cut = r_cut + self.redshift_object = redshift_object + self.redshift_source = redshift_source + + +class dPIEMassLenstoolSph(dPIEMassSph): + """ + The spherical dPIE mass profile in Lenstool's native parameterization, supporting + model-fitting with priors placed directly on the Lenstool parameters + (``sigma``, ``r_core``, ``r_cut``). See :class:`dPIEMassLenstool`. + + Parameters + ---------- + centre : (float, float) + (y, x) arc-second coordinates of the profile centre. + sigma : float + Lenstool fiducial velocity dispersion ``v_disp`` (sigma_LT) in km/s. + r_core : float + Lenstool ``core_radius`` in arcseconds (the dPIE ``ra``). + r_cut : float + Lenstool ``cut_radius`` in arcseconds (the dPIE ``rs``). + redshift_object : float + The redshift of the lens, used for the D_LS / D_S normalization of ``b0``. + redshift_source : float + The redshift of the source used to normalize ``b0``. + """ + + def __init__( + self, + centre: Tuple[float, float] = (0.0, 0.0), + sigma: float = 200.0, + r_core: float = 0.1, + r_cut: float = 20.0, + redshift_object: float = 0.5, + redshift_source: float = 1.0, + ): + from autogalaxy.cosmology.model import Planck15 + + cosmology = Planck15() + + b0 = _b0_from_lenstool_sigma( + sigma=sigma, + redshift_object=redshift_object, + redshift_source=redshift_source, + cosmology=cosmology, + ) + + super().__init__( + centre=centre, + ra=r_core, + rs=r_cut, + b0=b0, + ) + + self.sigma = sigma + self.r_core = r_core + self.r_cut = r_cut + self.redshift_object = redshift_object + self.redshift_source = redshift_source diff --git a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_mass.py b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_mass.py index af12334f9..6f26d5073 100644 --- a/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_mass.py +++ b/test_autogalaxy/profiles/mass/total/test_dual_pseudo_isothermal_mass.py @@ -1,74 +1,308 @@ -import pytest - -import autogalaxy as ag - -grid = ag.Grid2DIrregular([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0], [2.0, 4.0]]) - - -def test__deflections_yx_2d_from__sph_config_1(): - mp = ag.mp.dPIEMassSph(centre=(-0.7, 0.5), b0=5.2, ra=2.0, rs=3.0) - - deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) - - assert deflections[0, 0] == pytest.approx(1.033080741, 1e-4) - assert deflections[0, 1] == pytest.approx(-0.39286169026, 1e-4) - - -def test__deflections_yx_2d_from__sph_config_2(): - mp = ag.mp.dPIEMassSph(centre=(-0.1, 0.1), b0=20.0, ra=2.0, rs=3.0) - - deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) - - assert deflections[0, 0] == pytest.approx(1.4212977207, 1e-4) - assert deflections[0, 1] == pytest.approx(0.308977765378, 1e-4) - - -def test__deflections_yx_2d_from__elliptical(): - # First deviation from potential case due to ellipticity - - mp = ag.mp.dPIEMass( - centre=(0, 0), ell_comps=(0.0, 0.333333), b0=4.0, ra=2.0, rs=3.0 - ) - - deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) - - assert deflections[0, 0] == pytest.approx(0.21461366, 1e-3) - assert deflections[0, 1] == pytest.approx(0.10753914, 1e-3) - - -def test__deflections_yx_2d_from__elliptical_vs_spherical(): - elliptical = ag.mp.dPIEMass( - centre=(1.1, 1.1), ell_comps=(0.000001, 0.0000001), b0=12.0, ra=2.0, rs=3.0 - ) - spherical = ag.mp.dPIEMassSph(centre=(1.1, 1.1), b0=12.0, ra=2.0, rs=3.0) - - assert elliptical.deflections_yx_2d_from(grid=grid).array == pytest.approx( - spherical.deflections_yx_2d_from(grid=grid).array, 1e-1 - ) - - -def test__convergence_func__matches_private_helper(): - """Regression: dPIEMass must override the abstract `convergence_func` - so MGEDecomposer.decompose_convergence_via_mge (which walks the - convergence radially during MGE potential decomposition) doesn't - fall through to the abstract NotImplementedError. The shim delegates - to the existing `_convergence` radial helper that `convergence_2d_from` - already uses.""" - - import numpy as np - - mp = ag.mp.dPIEMass(centre=(0.0, 0.0), b0=5.2, ra=2.0, rs=3.0) - - # Scalar radius: equals the _convergence formula directly. - assert mp.convergence_func(1.5) == pytest.approx(mp._convergence(1.5), 1e-12) - - # 1-D array of radii: shape preserved, element-wise equality. - radii = np.array([0.1, 0.5, 1.0, 2.5, 5.0]) - expected = mp._convergence(radii) - actual = mp.convergence_func(radii) - assert actual.shape == radii.shape - assert actual == pytest.approx(expected, 1e-12) - - # dPIEMassSph inherits the override from dPIEMass. - sph = ag.mp.dPIEMassSph(centre=(0.0, 0.0), b0=5.2, ra=2.0, rs=3.0) - assert sph.convergence_func(1.5) == pytest.approx(sph._convergence(1.5), 1e-12) +import pytest + +import autogalaxy as ag + +grid = ag.Grid2DIrregular([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0], [2.0, 4.0]]) + + +def test__deflections_yx_2d_from__sph_config_1(): + mp = ag.mp.dPIEMassSph(centre=(-0.7, 0.5), b0=5.2, ra=2.0, rs=3.0) + + deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) + + assert deflections[0, 0] == pytest.approx(1.033080741, 1e-4) + assert deflections[0, 1] == pytest.approx(-0.39286169026, 1e-4) + + +def test__deflections_yx_2d_from__sph_config_2(): + mp = ag.mp.dPIEMassSph(centre=(-0.1, 0.1), b0=20.0, ra=2.0, rs=3.0) + + deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1875, 0.1625]])) + + assert deflections[0, 0] == pytest.approx(1.4212977207, 1e-4) + assert deflections[0, 1] == pytest.approx(0.308977765378, 1e-4) + + +def test__deflections_yx_2d_from__elliptical(): + # First deviation from potential case due to ellipticity + + mp = ag.mp.dPIEMass( + centre=(0, 0), ell_comps=(0.0, 0.333333), b0=4.0, ra=2.0, rs=3.0 + ) + + deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.1625, 0.1625]])) + + assert deflections[0, 0] == pytest.approx(0.21461366, 1e-3) + assert deflections[0, 1] == pytest.approx(0.10753914, 1e-3) + + +def test__deflections_yx_2d_from__elliptical_vs_spherical(): + elliptical = ag.mp.dPIEMass( + centre=(1.1, 1.1), ell_comps=(0.000001, 0.0000001), b0=12.0, ra=2.0, rs=3.0 + ) + spherical = ag.mp.dPIEMassSph(centre=(1.1, 1.1), b0=12.0, ra=2.0, rs=3.0) + + assert elliptical.deflections_yx_2d_from(grid=grid).array == pytest.approx( + spherical.deflections_yx_2d_from(grid=grid).array, 1e-1 + ) + + +def test__convergence_func__matches_private_helper(): + """Regression: dPIEMass must override the abstract `convergence_func` + so MGEDecomposer.decompose_convergence_via_mge (which walks the + convergence radially during MGE potential decomposition) doesn't + fall through to the abstract NotImplementedError. The shim delegates + to the existing `_convergence` radial helper that `convergence_2d_from` + already uses.""" + + import numpy as np + + mp = ag.mp.dPIEMass(centre=(0.0, 0.0), b0=5.2, ra=2.0, rs=3.0) + + # Scalar radius: equals the _convergence formula directly. + assert mp.convergence_func(1.5) == pytest.approx(mp._convergence(1.5), 1e-12) + + # 1-D array of radii: shape preserved, element-wise equality. + radii = np.array([0.1, 0.5, 1.0, 2.5, 5.0]) + expected = mp._convergence(radii) + actual = mp.convergence_func(radii) + assert actual.shape == radii.shape + assert actual == pytest.approx(expected, 1e-12) + + # dPIEMassSph inherits the override from dPIEMass. + sph = ag.mp.dPIEMassSph(centre=(0.0, 0.0), b0=5.2, ra=2.0, rs=3.0) + assert sph.convergence_func(1.5) == pytest.approx(sph._convergence(1.5), 1e-12) + + +def test__from_lenstool__b0_conversion__matches_isothermal_relation(): + # b0 = 6 * 648000 * (sigma_LT / c)^2 * (D_LS / D_S); cross-checked via the independent + # cosmology.velocity_dispersion_from (isothermal theta_E -> sigma_0), which must + # recover the central dispersion sigma_0 = sqrt(3/2) * sigma_LT. + + cosmology = ag.cosmo.Planck15() + + mp = ag.mp.dPIEMassSph.from_lenstool( + sigma=1000.0, + r_core=0.1, + r_cut=20.0, + redshift_object=0.5, + redshift_source=2.0, + cosmology=cosmology, + ) + + sigma_0 = cosmology.velocity_dispersion_from( + redshift_0=0.5, redshift_1=2.0, einstein_radius=mp.b0 + ) + + assert sigma_0 == pytest.approx(1000.0 * (1.5**0.5), rel=1e-6) + + +def test__from_lenstool__b0_explicit_value(): + # Pure prefactor check with the distance ratio divided out: + # b0 / (D_LS / D_S) = 6 * 648000 * (sigma / c)^2. + + cosmology = ag.cosmo.Planck15() + + mp = ag.mp.dPIEMassSph.from_lenstool( + sigma=250.0, + redshift_object=0.3, + redshift_source=1.5, + cosmology=cosmology, + ) + + d_s = cosmology.angular_diameter_distance_to_earth_in_kpc_from(redshift=1.5) + d_ls = cosmology.angular_diameter_distance_between_redshifts_in_kpc_from( + redshift_0=0.3, redshift_1=1.5 + ) + + assert mp.b0 / (d_ls / d_s) == pytest.approx( + 6.0 * 648000.0 * (250.0 / 299792.458) ** 2, rel=1e-8 + ) + + +def test__from_lenstool__ellipticity_conversion__matches_lenstool_internal(): + # Lenstool converts emass = (a^2-b^2)/(a^2+b^2) to epot = (1 - sqrt(1 - emass^2)) / emass + # (set_lens.c) and passes epot to ci05f; the port's _ellip() = |ell_comps| must equal it. + + emass = 0.4 + + mp = ag.mp.dPIEMass.from_lenstool(ellipticity=emass, angle_pos=0.0) + + epot_lenstool = (1.0 - (1.0 - emass**2) ** 0.5) / emass + + assert mp._ellip() == pytest.approx(epot_lenstool, rel=1e-10) + + axis_ratio = ((1.0 - emass) / (1.0 + emass)) ** 0.5 + + assert mp._ellip() == pytest.approx( + (1.0 - axis_ratio) / (1.0 + axis_ratio), rel=1e-10 + ) + + +def test__from_lenstool__angle_and_radii_passthrough(): + mp = ag.mp.dPIEMass.from_lenstool( + centre=(0.1, -0.2), + ellipticity=0.3, + angle_pos=45.0, + sigma=200.0, + r_core=0.15, + r_cut=12.0, + ) + + assert mp.centre == (0.1, -0.2) + assert mp.ra == 0.15 + assert mp.rs == 12.0 + + axis_ratio = ((1.0 - 0.3) / (1.0 + 0.3)) ** 0.5 + ell_comps = ag.convert.ell_comps_from(axis_ratio=axis_ratio, angle=45.0) + + assert mp.ell_comps[0] == pytest.approx(ell_comps[0], rel=1e-10) + assert mp.ell_comps[1] == pytest.approx(ell_comps[1], rel=1e-10) + + +def test__from_lenstool__sph_matches_elliptical_zero_ellipticity(): + kwargs = dict( + sigma=300.0, r_core=0.2, r_cut=15.0, redshift_object=0.4, redshift_source=1.8 + ) + + sph = ag.mp.dPIEMassSph.from_lenstool(**kwargs) + ell = ag.mp.dPIEMass.from_lenstool(ellipticity=0.0, angle_pos=0.0, **kwargs) + + grid_check = ag.Grid2DIrregular([[0.5, 0.8], [1.5, -0.3]]) + + sph_deflections = sph.deflections_yx_2d_from(grid=grid_check) + ell_deflections = ell.deflections_yx_2d_from(grid=grid_check) + + assert sph_deflections[0, 0] == pytest.approx(ell_deflections[0, 0], rel=1e-4) + assert sph_deflections[1, 1] == pytest.approx(ell_deflections[1, 1], rel=1e-4) + + +def test__from_lenstool__isothermal_limit_deflection_equals_b0(): + # ra -> 0, rs -> inf, e = 0: the dPIE tends to a SIS whose deflection magnitude is b0 + # everywhere, so b0 is the Einstein radius in this limit. + + mp = ag.mp.dPIEMassSph.from_lenstool( + sigma=800.0, r_core=1e-5, r_cut=1e6, redshift_object=0.5, redshift_source=2.0 + ) + + deflections = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.0, 3.0]])) + + assert deflections[0, 1] == pytest.approx(mp.b0, rel=1e-3) + + +def test__potential_2d_from__gradient_matches_deflections(): + # The analytic potential (Lenstool pi05 port) must be exactly consistent with the + # ci05f deflections: finite-difference grad(psi) = alpha, including rotation. + + mp = ag.mp.dPIEMass( + centre=(0.1, -0.3), ell_comps=(0.15, 0.1), ra=1.5, rs=30.0, b0=10.0 + ) + + eps_fd = 1e-6 + y, x = 0.8, 1.3 + + grid_fd = ag.Grid2DIrregular( + [ + [y, x + eps_fd], + [y, x - eps_fd], + [y + eps_fd, x], + [y - eps_fd, x], + ] + ) + + psi = mp.potential_2d_from(grid=grid_fd) + + alpha_fd_x = (psi[0] - psi[1]) / (2 * eps_fd) + alpha_fd_y = (psi[2] - psi[3]) / (2 * eps_fd) + + alpha = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[y, x]])) + + assert alpha_fd_y == pytest.approx(alpha[0, 0], rel=1e-5) + assert alpha_fd_x == pytest.approx(alpha[0, 1], rel=1e-5) + + +def test__potential_2d_from__spherical_limit_matches_quadrature(): + # For a near-circular profile psi(R) = 2 * int_0^R kappa(r) r ln(R/r) dr + ... ; + # simpler and equivalent: d(psi)/dR must equal the spherical deflection alpha(R). + + import numpy as np + + mp = ag.mp.dPIEMassSph(centre=(0.0, 0.0), ra=2.0, rs=20.0, b0=5.0) + + eps_fd = 1e-6 + R = 4.0 + + psi = mp.potential_2d_from( + grid=ag.Grid2DIrregular([[0.0, R + eps_fd], [0.0, R - eps_fd]]) + ) + alpha_fd = (psi[0] - psi[1]) / (2 * eps_fd) + + alpha = mp.deflections_yx_2d_from(grid=ag.Grid2DIrregular([[0.0, R]])) + + assert alpha_fd == pytest.approx(alpha[0, 1], rel=1e-5) + + +def test__lenstool_wrapper__matches_from_lenstool_constructor(): + # The model-fittable wrapper class must produce the identical profile to the + # from_lenstool classmethod with the (fixed Planck15) cosmology. + + kwargs = dict( + centre=(0.1, -0.2), + ellipticity=0.3, + angle_pos=45.0, + sigma=350.0, + r_core=0.15, + r_cut=12.0, + redshift_object=0.4, + redshift_source=1.8, + ) + + wrapper = ag.mp.dPIEMassLenstool(**kwargs) + constructed = ag.mp.dPIEMass.from_lenstool(**kwargs) + + assert wrapper.b0 == pytest.approx(constructed.b0, rel=1e-10) + assert wrapper.ra == constructed.ra + assert wrapper.rs == constructed.rs + assert wrapper.ell_comps[0] == pytest.approx(constructed.ell_comps[0], rel=1e-10) + assert wrapper.ell_comps[1] == pytest.approx(constructed.ell_comps[1], rel=1e-10) + + grid_check = ag.Grid2DIrregular([[0.5, 0.8]]) + + assert wrapper.deflections_yx_2d_from(grid=grid_check)[0, 0] == pytest.approx( + constructed.deflections_yx_2d_from(grid=grid_check)[0, 0], rel=1e-10 + ) + + sph_kwargs = { + k: v for k, v in kwargs.items() if k not in ("ellipticity", "angle_pos") + } + + wrapper_sph = ag.mp.dPIEMassLenstoolSph(**sph_kwargs) + constructed_sph = ag.mp.dPIEMassSph.from_lenstool(**sph_kwargs) + + assert wrapper_sph.b0 == pytest.approx(constructed_sph.b0, rel=1e-10) + + +def test__lenstool_wrapper__supports_model_composition(): + # Fitting in Lenstool parameters: af.Model must resolve priors for every __init__ + # arg from the config, and fixing the redshifts must leave the Lenstool free + # parameters (centre_0, centre_1, ellipticity, angle, sigma, r_core, r_cut). + + import autofit as af + + model = af.Model(ag.mp.dPIEMassLenstool) + + assert model.prior_count == 9 + + model.redshift_object = 0.5 + model.redshift_source = 2.0 + + assert model.prior_count == 7 + + instance = model.instance_from_unit_vector([0.5] * model.prior_count) + + assert isinstance(instance, ag.mp.dPIEMassLenstool) + assert instance.b0 > 0.0 + + model_sph = af.Model(ag.mp.dPIEMassLenstoolSph) + + assert model_sph.prior_count == 7